Warning, /utils/matlab/cs_grid/bk_line/grt_circ_bkl.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
28aaff1409 Jean*0001 % Script to generate the broken line that folows the cubic grid and
0002 % stay as close as possible to the great circle arc between 2 given points A,B.
0003 % Transport computed from this broken line (across great circle arc AB) is
0004 % positive when the flow goes from right to left of arc A -> B.
0005 % (e.g.: if lat_A = lat_B and long_A < long_B, northward flow -> transport > 0 )
0006 %-----------------------
0007 %- kwr=1 : => write broken-line to file: filnam(.mat)
0008 %- kplot : to check different steps of making bk-line:
0009 % < 0 : no plot ; 0 : just plot cos(longitude) after rotation
0010 % = 1 : add bk-line after "find_bk_line" ;
0011 % = 2 : same after "clean_bk_line" ;
0012 % = 3 : same after "save_bk_line" ;
0013 % = 4 : same after "shift_bk_line" (= final step);
0014
0015 kwr=0; kplot=4;
0016
0c7ba451ba Jean*0017 %- set ncdf=1 to load MNC (NetCDF) grid-files ;
0018 % or ncdf=0 to load MDS (binary) grid-files :
0019 ncdf=0;
0020 rac='grid_files/';
0021 G=load_grid(rac,10+ncdf);
0022 xc0=G.xC; yc0=G.yC; xg0=G.xG; yg0=G.yG; arc=G.rAc;
0023 nc=G.dims(2); ; ncp=nc+1 ;
28aaff1409 Jean*0024 nXY=6*nc*nc;
0025 nXYz=nXY+2;
0026
0027 xa0=-10;
0028 ya0=-45;
0029 xb0=80;
0030 yb0=45;
0031 filnam='bkl_AB';
0032
0033 %xa0=10;
0034 %ya0=2;
0035 %xb0=30;
0036 %yb0=2;
0037
0038 %----------------------------------
0039 rad=pi/180;
0040
0041 xg0=reshape(xg0,[nXY 1]);
0042 yg0=reshape(yg0,[nXY 1]);
0043 %- add th 2 missing corner:
0044 xg0(nXY+1)=xg0(1); yg0(nXY+1)=yg0(1+2*nc);
0045 xg0(nXY+2)=xg0(1+3*nc); yg0(nXY+2)=yg0(1);
0046
0047 %----------------------------------
0048 dx=(xg0-xa0);
0049 dy=(yg0-ya0);
0050 dr=dx.*dx+dy.*dy;
0051 drMn=min(dr(:)); [I J]=find(dr==drMn);
0052 if length(I) > 1,
0053 fprintf(' find more than 1 point next to A\n')
0054 else
0c7ba451ba Jean*0055 if length(I) ==0,
28aaff1409 Jean*0056 fprintf(' did not find any point next to A\n')
0057 return
0058 end
0059 end
0060 xa=xg0(I(1),J(1));
0061 ya=yg0(I(1),J(1));
0062
0063 dx=(xg0-xb0);
0064 dy=(yg0-yb0);
0065 dr=dx.*dx+dy.*dy;
0066 drMn=min(dr(:)); [I J]=find(dr==drMn);
0067 if length(I) > 1,
0068 fprintf(' find more than 1 point next to B\n')
0069 else
0c7ba451ba Jean*0070 if length(I) ==0,
28aaff1409 Jean*0071 fprintf(' did not find any point next to B\n')
0072 return
0073 end
0074 end
0075 xb=xg0(I(1),J(1));
0076 yb=yg0(I(1),J(1));
0077
0078 fprintf('A point: %8.3f %8.3f\n',xa,ya);
0079 fprintf('B point: %8.3f %8.3f\n',xb,yb);
0080 %----
0081
0082 [xcg,ycg,xcs,ycs,alp,bet]=rotate_xy(xg0,yg0,xc0,yc0,xa,ya,xb,yb);
0083
0084 %-----------------
0085
0086 xx1=split_Z_cub(xcg);
0087 yy2=split_Z_cub(ycg);
0088 %-- define dylat = width of Lat. band that contains the broken line:
0089 ddy=zeros(ncp,ncp,6);
0090 ddy(1:nc,:,:)=yy2(1:nc,:,:)-yy2(2:ncp,:,:); ddy=abs(ddy);
0091 dyImx=max(ddy(:)) ;
0092 ddy=zeros(ncp,ncp,6);
0093 ddy(:,1:nc,:)=yy2(:,1:nc,:)-yy2(:,2:ncp,:); ddy=abs(ddy);
0094 dyJmx=max(ddy(:));
0095 dylat=sqrt(2)*max(dyImx,dyJmx);
0096 fprintf('dyImx,dyJmx = %7.3f %7.3f ; ==> select dylat= %9.6f \n', ...
0097 dyImx,dyJmx,dylat);
0098
0099 %-----------------
0100 % to avoid -180,+180 jump in long, center longitude arround xMid:
0101 xMid=zeros(6,1);
0102 for n=1:6,
0103 var=reshape(yy2(:,:,n),[ncp*ncp 1]);
0104 [I]=find(abs(var) <= dylat);
0105 fprintf(' face: %i, %i points in dylat band',n,length(I));
0106 if length(I) > 0,
0107 var=reshape(xx1(:,:,n),[ncp*ncp 1]);
0108 xx=var(I)*rad; xM=mean(cos(xx)); yM=mean(sin(xx));
0c7ba451ba Jean*0109 xMid(n)=atan2(yM,xM)/rad;
28aaff1409 Jean*0110 fprintf(', xMid= %8.3f\n',xMid(n));
0111 end
0112 end
0113 xMid=round(xMid);
0114 xx2=xx1;
0115 for n=1:6,
0116 if xMid(n) ~= 0,
0117 dec=xMid(n)-180;
0118 xx2(:,:,n)=rem( xx2(:,:,n)+360-dec, 360) +dec ;
0119 end
0120 fprintf('face %i : xMid= %8.2f ; Xmin,Xmax = %8.3f %8.3f \n', ...
0121 n, xMid(n), min(min(xx2(:,:,n))), max(max(xx2(:,:,n))) );
0122 end
0123
0124 yl=0; nf1=1; nf2=6; ydim=1; jl=1;
0125 XYout=1000;
0126 %nMx6t=zeros(6,1);
0127 %savI=zeros(nc*6,6); savJ=zeros(nc*6,6); savF=zeros(nc*6,6);
0128 %isav=zeros(nc*6,6); jsav=zeros(nc*6,6); xsav=zeros(nc*6,6);
0c7ba451ba Jean*0129 %find_bk_line
28aaff1409 Jean*0130 [savI,savJ,savF,isav,jsav,xsav,nMx6t] = ...
0131 find_bk_line( nf1,nf2,nc,ydim,yl,dylat,XYout,xMid,xx1,xx2,yy2,xcs,ycs );
0132
0133 %- define "segments" = continuous part of the broken line :
0134 %ncut=zeros(6,1); icut=zeros(nc,6,6); xcut=zeros(nc,4,6); ycut=zeros(nc,4,6);
0135 %clean_bk_line
0136 [ncut,icut,xcut,ycut,misfit,xyfit] = ...
0137 clean_bk_line( nf1,nf2,nc,ydim,yl,dylat,xMid,xx1,xx2,yy2, ...
0138 savI,savJ,savF,isav,jsav,xsav,nMx6t );
0139
0140 if misfit > 0,
0141 fprintf('misfit= %i , xyfit= %i ; ==> must do something ! \n', ...
0142 misfit,xyfit);
0143 return
0144 end
0145
0146 %-----------------
0147
0148 xx0=split_Z_cub(xg0);
0149 yy0=split_Z_cub(yg0);
0150
0151 %- output : put together the pieces of bk-lines from the 6 faces :
0152 %save_bk_line
0153 savNpts=zeros(ydim,1);
0154 savFlg=zeros(6*ncp,ydim);
0155 savIuv=zeros(6*ncp,ydim); savJuv=zeros(6*ncp,ydim);
0156 savXsg=zeros(6*ncp,ydim); savYsg=zeros(6*ncp,ydim);
0157
6f2602a7ee Jean*0158 [svNpts,svFlg,svIuv,svJuv,svXsg,svYsg,svXx1,svYy1]= ...
28aaff1409 Jean*0159 save_bk_line( nf1,nf2,nc,ydim,yl,dylat,XYout,xMid,xx0,yy0,yy2, ...
0160 savI,savJ,savF,isav,jsav,xsav,ncut,icut,xcut,ycut );
0161 %- easier to debug this way:
0162 savNpts(jl)=svNpts;
0163 savFlg(:,jl)=svFlg;
0164 savIuv(:,jl)=svIuv;
0165 savJuv(:,jl)=svJuv;
0166 savXsg(:,jl)=svXsg;
0167 savYsg(:,jl)=svYsg;
0168
0169 [ab_Npts,ab_Flg,ab_IJuv,ab_Xsg,ab_Ysg] = ...
0170 shift_bk_line( nc,ydim,jl,xa,ya,xb,yb,savNpts,savFlg,savIuv,savJuv,savXsg,savYsg );
0171
0172 if kwr == 1,
0173 %- save to matlab file "filnam":
0174 filnam='bkl_AB';
0175 namf=[filnam,'_cs',int2str(nc)];
0176 fprintf(['save bk-line to file: ',namf,'\n']);
0177 save(namf,'ab_Npts','ab_Flg','ab_IJuv','ab_Xsg','ab_Ysg');
0178 end
0179 %-----------------
0180 if kplot < 0 , return ; end
0181
0182 ccB=[0 0]; shift=-1; cbV=1; AxBx=[-180 180 -90 90];
0183 %AxBx=[-50 30 -50 -10];
0184 %var=ycs;
0185 var=cos(xcs*rad); ccB=[-1 1]*3;
0186 figure(1);clf;
0187 grph_CS(var,xc0,yc0,xg0,yg0,ccB(1),ccB(2),shift,cbV,AxBx);
0188 hold on;
0189 [L]=line([xa xb],[ya yb]); set(L,'Color',[1 1 1]);
0190 plot(alp/rad,bet/rad,'ro-');
0191 %nj=21; plot(xx0(isav(nj,1),jsav(nj,1),1)+1,yy0(isav(nj,1),jsav(nj,1),1)+1,'cx');
0192 if kplot == 1,
0193 for n=nf1:nf2, if nMx6t(n) > 0,
0194 xloc=zeros(nMx6t(n),1); yloc=xloc;
0195 for i=1:nMx6t(n),
0196 xloc(i)=xx0(isav(i,n),jsav(i,n),n);
0197 yloc(i)=yy0(isav(i,n),jsav(i,n),n);
0198 end
0199 plot(xloc,yloc,'g*-');
0200 end ; end
0201 end
0202 if kplot == 2,
0203 clrs='bscdmvr^y<g>kp';
0204 clear P ; np=0;
0205 for n=nf1:nf2, if ncut(n) > 0,
0206 xloc=zeros(6*nc,1); yloc=zeros(6*nc,1);
0207 for i=1:nMx6t(n),
0208 xloc(i)=xx0(isav(i,n),jsav(i,n),n);
0209 yloc(i)=yy0(isav(i,n),jsav(i,n),n);
0210 end
0211 for in=1:ncut(n), if icut(in,6,n) >= 0,
0212 is=3; if icut(in,is,n) <= 0, is=1 ; end
0213 ie=4; if icut(in,ie,n) <= 0, ie=2 ; end
0214 if icut(in,is,n) > icut(in,ie,n), is=1 ; ie=2 ; end
0215 if icut(in,is,n) > 0 & icut(in,ie,n) > 0 & icut(in,is,n) < icut(in,ie,n),
0216 np=np+1;
0217 % P(np)=plot(xsav(icut(in,is,n):icut(in,ie,n),n), ...
0218 P(np)=plot(xloc(icut(in,is,n):icut(in,ie,n)), ...
0219 yloc(icut(in,is,n):icut(in,ie,n)),[clrs(1+2*(np-1):2*np),'-']);
0220 % yloc(icut(in,is,n):icut(in,ie,n)),'g*-');
0221 end
0222 end; end
0223 end; end
0224 set(P,'LineWidth',4);
0225 end
0226 if kplot == 3,
0227 clear P ; np=1; ict(1)=0;
0228 sNp=savNpts(jl)+1;
0229 xloc=savXsg(1:sNp,jl);
0230 yloc=savYsg(1:sNp,jl);
0231 %xloc=ab_Xsg; yloc=ab_Ysg;
0232 dd=abs(xloc(2:sNp)-xloc(1:sNp-1));
0233 fprintf(' np, Max(dd) : %i %8.3f\n',np,max(dd(:)) );
0234 while max(dd(:)) > 2*dylat,
0235 [I]=find(dd==max(dd(:)));
0236 ict(np+1)=I;
0237 xloc(I+1:sNp+1)=xloc(I:sNp);
0238 yloc(I+1:sNp+1)=yloc(I:sNp);
0c7ba451ba Jean*0239 sNp=sNp+1; np=np+1;
28aaff1409 Jean*0240 dec=xloc(I+2)-180;
0241 xloc(I+1)=dec+rem(xloc(I+1)-dec+360,360);
0242 for j=1:np-1, if ict(j) > I, ict(j)=ict(j)+1; end; end
0243 dd=abs(xloc(2:sNp)-xloc(1:sNp-1));
0244 for j=2:np, dd(ict(j))=0; end
0245 fprintf(' np, Max(dd) : %i %8.3f\n',np,max(dd(:)) );
0246 end
0247 ict(np+1)=sNp;
0248 for j=1:np,
0249 P(j)=plot(xloc(1+ict(j):ict(j+1)),yloc(1+ict(j):ict(j+1)),'r*-');
0250 end
0251 set(P,'LineWidth',4);
0252 end
0253 if kplot == 4,
0254 clear P ; sNp=ab_Npts+1;
0255 xloc=ab_Xsg(1:sNp);
0256 yloc=ab_Ysg(1:sNp);
0257 P=plot(xloc,yloc,'b*-');
0258 set(P,'LineWidth',4);
0259 end
0260 hold off
0261 title('COS( new-long )');
