Warning, /utils/matlab/cs_grid/bk_line/gener_bk_line.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
323aa12fe2 Jean*0001 % main script to generate broken lines that folows the cubic grid
0002 % and stay close as possible to a given latitude
0003 %-----------------------
0004 %- load definition of the grid (needs to be done at the 1rst call):
0005
0006 if size(who('krd'),1) > 0,
0007 fprintf('krd is defined and = %i \n',krd);
0008 else
0009 fprintf('krd undefined ; set to 1 \n'); krd=1 ;
0c7ba451ba Jean*0010 end
0011 if krd==1,
9194d1ae0e Jean*0012 more off ;
0c7ba451ba Jean*0013 %- set ncdf=1 to load MNC (NetCDF) grid-files ;
0014 % or ncdf=0 to load MDS (binary) grid-files :
0015 ncdf=0;
323aa12fe2 Jean*0016 %rac='/home/jmc/grid_cs32/';
0017 rac='grid_files/';
0c7ba451ba Jean*0018 G=load_grid(rac,10+ncdf);
0019 xcs=G.xC; ycs=G.yC; xcg=G.xG; ycg=G.yG; arc=G.rAc;
323aa12fe2 Jean*0020 %- 1rst try:
0c7ba451ba Jean*0021 ydim=1; jl=1; yl=35; kwr=0;
323aa12fe2 Jean*0022 else
0023 %- do the real calculation and write to file:
0024 ydim=2; kwr=1;
0025 end
0026
0027 %------------
0028 nc=size(xcs,2) ; ncp=nc+1 ;
0029
0030 %- output arrays:
0031 ylat=-87:3:87;
0032 if ydim > 1, ydim=length(ylat); end
0033 savNpts=zeros(ydim,1);
28aaff1409 Jean*0034 savFlg=zeros(6*ncp,ydim);
0035 savIuv=zeros(6*ncp,ydim); savJuv=zeros(6*ncp,ydim);
0036 savXsg=zeros(6*ncp,ydim); savYsg=zeros(6*ncp,ydim);
323aa12fe2 Jean*0037
0038 kplot=zeros(ydim,1);
28aaff1409 Jean*0039 kplot(1:ydim)=1;
0040 %kplot(1)=1;
323aa12fe2 Jean*0041 nf1=1; nf2=6; %- allow to treat only few faces (debug) : n=nf1:nf2,
0042
0043 %------------------------------
0044 %- face 3 & 6 : long cover ]-180;180]
0c7ba451ba Jean*0045 %- face 4 : long in 2 part ]-180;-135] + ]135;180];
323aa12fe2 Jean*0046 % => put together (in xx2) with xMid=180
0c7ba451ba Jean*0047 xMid=zeros(6,1);
0048 xMid(4)=180;
323aa12fe2 Jean*0049
0050 %-- define large value for X,Y position :
0051 XYout=1000;
0052
9194d1ae0e Jean*0053 nPxy=nc*6*nc; nPp2=nPxy+2;
0054 yg2=zeros(nPp2,1); yg2([1:nPxy],1)=reshape(ycg,[nPxy 1]);
0055 xg2=zeros(nPp2,1); xg2([1:nPxy],1)=reshape(xcg,[nPxy 1]);
0056 %-- add missing corner:
0057 xg2(nPxy+1)=xg2(1); yg2(nPxy+1)=yg2(1+2*nc);
0058 xg2(nPxy+2)=xg2(1+3*nc); yg2(nPxy+2)=yg2(1);
323aa12fe2 Jean*0059
9194d1ae0e Jean*0060 %------------------------------
0061 yy2=split_Z_cub(yg2);
323aa12fe2 Jean*0062 %-- define dylat = width of Lat. band that contains the broken line:
0c7ba451ba Jean*0063 ddy=zeros(ncp,ncp,6);
323aa12fe2 Jean*0064 ddy(1:nc,:,:)=yy2(1:nc,:,:)-yy2(2:ncp,:,:); ddy=abs(ddy);
0065 dyImx=max(max(max(ddy))) ;
0c7ba451ba Jean*0066 ddy=zeros(ncp,ncp,6);
323aa12fe2 Jean*0067 ddy(:,1:nc,:)=yy2(:,1:nc,:)-yy2(:,2:ncp,:); ddy=abs(ddy);
0068 dyJmx=max(max(max(ddy)));
0069 dylat=sqrt(2)*max(dyImx,dyJmx);
0070 fprintf('dyImx,dyJmx = %7.3f %7.3f ; ==> select dylat= %9.6f \n', ...
0071 dyImx,dyJmx,dylat);
9194d1ae0e Jean*0072 %-- set long precision ; deal with long=-180 / +180 separation
0073 epsil=dylat*1.e-10;
0074 dd0=abs(xg2+180); ddd=dd0; ddd(find(ddd<epsil))=200; dd1=min(ddd(:));
0075 ddd=abs(xg2-180); ddd(find(ddd<epsil))=200; dd2=min(ddd(:));
0076 if ( dd1 < dylat/10 | dd2 < dylat/10 ),
0077 fprintf(' date change line is poorly defined : dd1= %g , dd2= %g\n',dd1,dd2);
0078 fprintf(' needs to be +/-180 (+/- esilon=%10.3e) ==> Stop here\n',epsil);
0079 return
0080 end
0081 %-- change long=-180 (+/- epsil) to +180
0082 [II]=find(dd0<epsil); xg2(II)=xg2(II)+360; clear II;
0083 %---------
0084 xx1=split_Z_cub(xg2);
0085 %- dx,dy not used but would be necessary to compute the length of broken-lines
0086 % [dy6t,dx6t]=split_UV_cub(dyg,dxg);
0087 xx2=xx1;
0088 for n=1:6,
0089 if xMid(n) ~= 0,
0090 dec=xMid(n)-180;
0091 xx2(:,:,n)=rem( xx2(:,:,n)+360-dec, 360) +dec ;
0092 end
0093 fprintf('face %i : x1 min,max= %8.3f %8.3f', ...
0094 n, min(min(xx1(:,:,n))), max(max(xx1(:,:,n))) );
0095 fprintf(' ; xMid= %8.3f ; x2 min,max= %8.3f %8.3f \n', ...
0096 xMid(n), min(min(xx2(:,:,n))), max(max(xx2(:,:,n))) );
0097 end
323aa12fe2 Jean*0098 %------------------------------
0099
0100 %-----------------------------------------
0101 %-- illustrate with a 2.D plot :
0102 if max(kplot) > 0,
0103 face_color=zeros(7,3); lineThick=3*ones(7);
28aaff1409 Jean*0104 face_color(1,:)=[1 0 0]; face_color(4,:)=face_color(1,:);
0105 face_color(2,:)=[0 1 0]; face_color(5,:)=face_color(2,:);
0106 face_color(3,:)=[0 0 1]; face_color(6,:)=face_color(3,:);
0107 face_color(7,:)=[1 0 1]; lineThick(7)=2;
323aa12fe2 Jean*0108
0109 figure(1); clf;
0110 %subplot(212);
9194d1ae0e Jean*0111 ccB=[-25 40]; shift=-1; cbV=2; AxBx=[-180 180 -90 90]; kEnv=0;
0112 if ydim == 1, AxBx(3:4)=[-15 15]+yl ; end
0113 grph_CS(zeros(6*nc,nc),xcs,ycs,xg2,yg2,ccB(1),ccB(2),shift,cbV,AxBx,kEnv);
0c7ba451ba Jean*0114 AA=axis;
323aa12fe2 Jean*0115 if ydim == 1,
0116 h(1)=line(AA(1:2),[yl yl]);
0117 h(2)=line(AA(1:2),[yl-dylat/2 yl-dylat/2]);
0118 h(3)=line(AA(1:2),[yl+dylat/2 yl+dylat/2]);
0119 set(h,'Color',[0 0 0],'LineWidth',1,'LineStyle','--');
0120 end
0121 hold on
0122 end
0123 %------------------------------
0124
0c7ba451ba Jean*0125 %-- Start loop on jl :
323aa12fe2 Jean*0126 for jl=1:ydim,
0127 if ydim > 1, yl=ylat(jl) ; end
0128 fprintf(' Define broken line closest to yl= %8.3f \n',yl);
0129 %======================================================================
0130
0131 %--------------------------------------
0132 % 1rst step : fort each "yl" and for each face;
0133 % find the broken-line closest to yl, starting from the West side = min(x)
0134
28aaff1409 Jean*0135 %nMx6t=zeros(6,1);
0136 %savI=zeros(nc*6,6); savJ=zeros(nc*6,6); savF=zeros(nc*6,6);
0137 %isav=zeros(nc*6,6); jsav=zeros(nc*6,6); xsav=zeros(nc*6,6);
0138 %find_bk_line
0139 [savI,savJ,savF,isav,jsav,xsav,nMx6t] = ...
0140 find_bk_line( nf1,nf2,nc,ydim,yl,dylat,XYout,xMid,xx1,xx2,yy2,xcs,ycs );
323aa12fe2 Jean*0141
0142 %--------------------------------------
0143 %- define "segments" = continuous part of the broken line :
0144
28aaff1409 Jean*0145 %ncut=zeros(6,1); icut=zeros(nc,6,6); xcut=zeros(nc,4,6); ycut=zeros(nc,4,6);
0146 %clean_bk_line
0147 [ncut,icut,xcut,ycut,misfit,xyfit] = ....
0148 clean_bk_line( nf1,nf2,nc,ydim,yl,dylat,xMid,xx1,xx2,yy2, ...
0149 savI,savJ,savF,isav,jsav,xsav,nMx6t );
323aa12fe2 Jean*0150
0151 %-- plot the results:
0c7ba451ba Jean*0152 if kplot(jl) == 2 | ( max(kplot) > 0 & misfit > 0 ),
323aa12fe2 Jean*0153 for n=nf1:nf2, if ncut(n) > 0,
0154 yloc=zeros(6*nc,1);
0155 for i=1:nMx6t(n), yloc(i)=yy2(isav(i,n),jsav(i,n),n); end;
0c7ba451ba Jean*0156 clear P ; np=0;
0157 for in=1:ncut(n), if icut(in,6,n) >= 0,
323aa12fe2 Jean*0158 % is=3; ie=4;
0c7ba451ba Jean*0159 % if icut(in,is,n)*icut(in,ie,n)==0 | icut(in,is,n) >= icut(in,ie,n),
323aa12fe2 Jean*0160 % is=1; ie=2; end
0161 is=3; if icut(in,is,n) <= 0, is=1 ; end
0162 ie=4; if icut(in,ie,n) <= 0, ie=2 ; end
0163 if icut(in,is,n) > icut(in,ie,n), is=1 ; ie=2 ; end
0c7ba451ba Jean*0164 if icut(in,is,n) > 0 & icut(in,ie,n) > 0 & icut(in,is,n) < icut(in,ie,n),
323aa12fe2 Jean*0165 np=np+1;
0166 P(np)=plot(xsav(icut(in,is,n):icut(in,ie,n),n), ...
0167 yloc(icut(in,is,n):icut(in,ie,n)),'*');
9194d1ae0e Jean*0168 %fprintf( ...
0169 %'n,in= %i %i ; is,ict,x,y= %i %i %6.2f %6.2f ; ie,ict,x,y= %i %i %6.2f %6.2f\n', ...
0170 % n,in,is,icut(in,is,n),xsav(icut(in,is,n),n),yloc(icut(in,is,n)), ...
0171 % ie,icut(in,ie,n),xsav(icut(in,ie,n),n),yloc(icut(in,ie,n)) );
323aa12fe2 Jean*0172 end ; end; end
0173 if np > 0,
0174 set(P,'Color',face_color(n,:),'LineWidth',lineThick(n),'LineStyle','-'); end
0c7ba451ba Jean*0175 end; end
323aa12fe2 Jean*0176 end
0177
0c7ba451ba Jean*0178 if misfit > 0,
323aa12fe2 Jean*0179 fprintf('misfit= %i , xyfit= %i ; ==> must do something ! \n', ...
0c7ba451ba Jean*0180 misfit,xyfit);
323aa12fe2 Jean*0181 return
0182 end
0183
0184 %--------------------------------------
28aaff1409 Jean*0185 %- output : put together the pieces of bk-lines from the 6 faces :
0186
0c7ba451ba Jean*0187 %save_bk_line
6f2602a7ee Jean*0188 [svNpts,svFlg,svIuv,svJuv,svXsg,svYsg,svXx1,svYy1]= ...
28aaff1409 Jean*0189 save_bk_line( nf1,nf2,nc,ydim,yl,dylat,XYout,xMid,xx1,yy2,yy2, ...
0190 savI,savJ,savF,isav,jsav,xsav,ncut,icut,xcut,ycut );
0191 %- easier to debug this way:
0192 savNpts(jl)=svNpts;
0193 savFlg(:,jl)=svFlg;
0194 savIuv(:,jl)=svIuv;
0195 savJuv(:,jl)=svJuv;
6f2602a7ee Jean*0196 savXsg(:,jl)=svXx1;
0197 savYsg(:,jl)=svYy1;
28aaff1409 Jean*0198
0199 %-- check that this broken-line is different from any previous one :
0200 check_bk_line( nc,ydim,jl,ylat,savNpts,savFlg,savIuv,savJuv,savXsg,savYsg );
0201
0202 if kplot(jl) == 1 | ( max(kplot) > 0 & misfit > 0 ),
323aa12fe2 Jean*0203 titr=sprintf('nMx6t= %i %i %i %i %i %i ; yLat= %8.3f',nMx6t,yl);
0204 title(titr);
0205 if lineThick(7) > 0,
28aaff1409 Jean*0206 clear Pu ict
0207 np=1; ict(1)=0;
0208 sNp = savNpts(jl)+1;
0209 xloc=savXsg(1:sNp,jl);
0210 yloc=savYsg(1:sNp,jl);
0211 ddMx=0; % ddMx=dylat/cos(yl*pi/180);
0212 if ddMx > 0,
0213 dd=abs(xloc(2:sNp)-xloc(1:sNp-1));
0214 fprintf(' np, Max(dd) : %i %8.3f\n',np,max(dd(:)) );
0215 while max(dd(:)) > ddMx,
0216 [I]=find(dd==max(dd(:)));
0217 ict(np+1)=I(1);
0218 xloc(I+1:sNp+1)=xloc(I:sNp);
0219 yloc(I+1:sNp+1)=yloc(I:sNp);
0220 sNp=sNp+1; np=np+1;
0221 dec=xloc(I+2)-180;
0222 xloc(I+1)=dec+rem(xloc(I+1)-dec+360,360);
0223 for j=1:np-1, if ict(j) > I, ict(j)=ict(j)+1; end; end
0224 dd=abs(xloc(2:sNp)-xloc(1:sNp-1));
0225 for j=2:np, dd(ict(j))=0; end
0226 fprintf(' np, Max(dd) : %i %8.3f\n',np,max(dd(:)) );
0227 if np > 10, dd=0; end
0228 end
0229 end
0230 ict(np+1)=sNp;
0231 for j=1:np,
0232 Pu(j)=plot(xloc(1+ict(j):ict(j+1)),yloc(1+ict(j):ict(j+1)));
0233 end
0234 %Pu=plot(savXsg(1:nPts,jl),savYsg(1:nPts,jl));
0235 set(Pu,'Color',face_color(7,:),'LineWidth',lineThick(7),'LineStyle','-');
323aa12fe2 Jean*0236 end ; end
0237
0c7ba451ba Jean*0238 if misfit > 0,
323aa12fe2 Jean*0239 fprintf('misfit= %i , Pb in connecting segment jl,yLat= %i %8.3f \n', ...
0240 misfit,jl,yl);
0241 return
0242 end
0243
0244 %======================================================================
0245 end %- end jl loop
0246
0247 if max(kplot) > 0, hold off ; end
0248
0249 %------------------------------
0250 %-- save to a matlab file :
0251 if ydim > 1 & kwr == 1,
0252
0253 %- write smaller arrays : reduce size from nc*6 to Max(savNpts) :
0c7ba451ba Jean*0254
323aa12fe2 Jean*0255 mxNpt=1+max(savNpts);
0256 bkl_Flg=zeros(mxNpt,ydim);
0257 bkl_Iuv=zeros(mxNpt,ydim); bkl_Juv=zeros(mxNpt,ydim);
0258 bkl_Xsg=zeros(mxNpt,ydim); bkl_Ysg=zeros(mxNpt,ydim);
0259
0260 bkl_Ylat=ylat ;
0c7ba451ba Jean*0261 bkl_Npts=savNpts ;
323aa12fe2 Jean*0262 bkl_Flg=savFlg(1:mxNpt,:);
0263 bkl_Iuv=savIuv(1:mxNpt,:); bkl_Juv=savJuv(1:mxNpt,:);
0264 bkl_Xsg=savXsg(1:mxNpt,:); bkl_Ysg=savYsg(1:mxNpt,:);
0265
0266 dyMean=round((ylat(ydim)-ylat(1))/ydim) ;
0267 namf=['isoLat_cube',int2str(nc),'_',int2str(ydim)];
0268
0269 save(namf,'bkl_Ylat','bkl_Npts','bkl_Flg','bkl_Iuv','bkl_Juv','bkl_Xsg','bkl_Ysg');
0270
0271 fprintf([' ==> Save %i lines on matlab file: ',namf,'\n'],ydim);
0272
0273 end
0274 %------------------------------
0275
0276 return
0277 %------------------------------
