Warning, /utils/matlab/cs_grid/bk_line/find_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
28aaff1409 Jean*0001 function [savI,savJ,savF,isav,jsav,xsav,nMx6t]=find_bk_line( ...
0002 nf1,nf2,nc,ydim,yl,dylat,XYout,xMid,xx1,xx2,yy2,xcs,ycs);
0003 % 1rst step : fort each "yl" and for each face;
0004 % find the broken-line closest to yl, starting from the West side = min(x)
0005 %--------------------------------------
323aa12fe2 Jean*0006
28aaff1409 Jean*0007 if yl == 1, fprintf( ...
0008 '--- find_bk_line: nf1,nf2,nc,ydim,yl= %i %i %i %i %8.3f\n', ...
0009 nf1,nf2,nc,ydim,yl); end
0010 ncp=nc+1;
0011 nMx6t=zeros(6,1);
0012 savI=zeros(nc*6,6); savJ=zeros(nc*6,6); savF=zeros(nc*6,6);
0013 isav=zeros(nc*6,6); jsav=zeros(nc*6,6); xsav=zeros(nc*6,6);
323aa12fe2 Jean*0014
0015 for n=nf1:nf2,
0016
0017 %---
0018 % to find the 1rst point :
0019 % start from the midle point ; If no point => give up.
0020 % decrease X until encounter boundary :
0021 % set 1rst point = 1rst point over boundary
0022
0023 nPts=0; yl1=yl-dylat/2 ; yl2=yl+dylat/2 ;
0024 [Il,Jl]=find( yy2(:,:,n) >= yl1 & yy2(:,:,n) < yl2) ; Ntloc=size(Il,1);
0025 if Ntloc > 0,
0026
0027 xloc=zeros(Ntloc,1);
0028 for i=1:Ntloc, xloc(i)=xx2(Il(i),Jl(i),n) ; end
0029 Xmidle=mean(xloc); xloc=abs(xloc-Xmidle); Xmn=min(xloc);
0030 ix=find(xloc <= Xmn+dylat); Nxloc=size(ix,1) ;
0031 yloc=zeros(Nxloc,1);
0032 for j=1:Nxloc, yloc(j)=yy2(Il(ix(j)),Jl(ix(j)),n) ; end
0033 yloc=abs(yloc-yl); Ymn=min(yloc) ; jx=find(yloc==Ymn);
0034 i0=Il(ix(jx(1))); j0=Jl(ix(jx(1)));
0035 x0=xx2(i0,j0,n) ; y0=yy2(i0,j0,n) ;
0036
0037 if ydim == 1, fprintf('n=%i, i,j,X,Y(-1)= %2i %2i %8.3f %8.3f ;', ...
0038 n,i0,j0,x0,y0); end
0039 %-------------------------------------------------
0040 while Ymn < XYout/2,
0041
0042 %- save the 1rst point :
0043 if nPts == 0,
0044 nPts = 1 ;
0045 else
0046 %- save the next point :
0047 i1 = i0 ; j1 = j0 ;
0048 if i4c(1) == 1, i1 = i0 -1 ; end
0049 if i4c(1) == 2, i1 = i0 +1 ; end
0050 if i4c(1) == 3, j1 = j0 -1 ; end
0051 if i4c(1) == 4, j1 = j0 +1 ; end
0052 nPts=nPts+1; x1=xx2(i1,j1,n);
0053 if x1 > x0, x1=x1-360 ; Ymn=XYout; end
0054 i0 = i1 ; j0 = j1 ; x0=x1;
0055 end
0056
0057 if Ymn < XYout/2,
0058 %- look for the next point : try the 4 direction i-1,i+1,j-1,j+1 :
0059 y4c=XYout*ones(4,1); y0=yy2(i0,j0,n) ;
0060
0061 if i0 > 1, xx = xx2(i0-1,j0,n);
0062 ddx= xx-x0+3*180 ; ddx=rem(ddx,360)-180 ;
0063 if ddx < 0, y4c(1)=yy2(i0-1,j0,n) ; end ; end
0064 if i0 < ncp, xx = xx2(i0+1,j0,n);
0065 ddx= xx-x0+3*180 ; ddx=rem(ddx,360)-180 ;
0066 if ddx < 0, y4c(2)=yy2(i0+1,j0,n) ; end ; end
0067 if j0 > 1, xx = xx2(i0,j0-1,n);
0068 ddx= xx-x0+3*180 ; ddx=rem(ddx,360)-180 ;
0069 if ddx < 0, y4c(3)=yy2(i0,j0-1,n) ; end ; end
0070 if j0 < ncp, xx = xx2(i0,j0+1,n);
0071 ddx= xx-x0+3*180 ; ddx=rem(ddx,360)-180 ;
0072 if ddx < 0, y4c(4)=yy2(i0,j0+1,n) ; end ; end
0073 %-- avoid the North & South poles (since Long is not well define there):
0074 y4c(find( abs(y4c) == 90)) = XYout;
0075 y4c=abs(y4c - yl); Ymn=min(y4c) ; i4c=find(y4c==Ymn);
0076 end
0077 %if nPts > 0,
0078 %fprintf('nPts,y4c : %i %8.3f %8.3f %8.3f %8.3f ; x0,y0= %8.3f %8.3f\n', ...
0079 % nPts,y4c,x0,y0); end
0080 %if nPts > 53, return ; end
0081
0082 end; % (while)
0083
0084 % x0=xx2(i0,j0,n) ;
0085 %- x0 is already set to xx2(i0) or to xx2(i0)-360
0086 y0=yy2(i0,j0,n) ; Ymn=abs(y0-yl);
0087 if ydim == 1,
0088 fprintf(' nPts=%3i, i0,j0,x0,y0= %2i %2i %8.3f %8.3f \n', ...
0089 nPts,i0,j0,x0,y0); end
0090
0091 %-------------------------------------------------
0092
0093 while Ymn < XYout/2,
0094
0095 %- save the 1rstpoint :
0096 if nMx6t(n) == 0,
0097 nMx6t(n) = 1 ; isav(1,n)=i0 ; jsav(1,n) = j0 ;
28aaff1409 Jean*0098 xsav(1,n) = x0 ; if xMid(n) ~= 0, xsav(1,n)=xx1(i0,j0,n) ; end
323aa12fe2 Jean*0099
0100 else
0101 %- save the next point :
0102 i1 = i0 ; j1 = j0 ; ii=i0 ; jj=j0 ; is=nMx6t(n) ;
0103 if i4c(1) == 1, i1 = i0 -1 ; savF(is,n)=2 ; ii=i1 ; end
0104 if i4c(1) == 2, i1 = i0 +1 ; savF(is,n)=2 ; end
0105 if i4c(1) == 3, j1 = j0 -1 ; savF(is,n)=1 ; jj=j1 ; end
0106 if i4c(1) == 4, j1 = j0 +1 ; savF(is,n)=1 ; end
0107 x1=xx2(i1,j1,n); if x1 < x0, Ymn=XYout; savF(is,n)=0; end
0108 if Ymn < XYout/2,
0109 savI(is,n)=ii ; savJ(is,n)=jj ;
0110 %- cut !
0111 if ii == ncp | jj == ncp, savF(is,n)=0 ; end
0112 %- suppress point outside yl1,yl2 :
0113 if ( mod(i0,nc) == 1 | mod(j0,nc) == 1 ) & ...
0114 abs(yy2(i0,j0,n)-yl) > dylat/2, savF(is,n)=0 ; end
0115 if ( mod(i1,nc) == 1 | mod(j1,nc) == 1 ) & ...
0116 abs(yy2(i1,j1,n)-yl) > dylat/2, savF(is,n)=0 ; end
0117
0118 %- if jump in xx1 => cut and add 1 point :
28aaff1409 Jean*0119 if xMid(n) ~= 0 & xx1(i1,j1,n) < xx1(i0,j0,n),
323aa12fe2 Jean*0120 savF(is+1,n)=savF(is,n); savF(is,n)=0; is=1+is;
0121 isav(is,n)=i0 ; jsav(is,n)=j0 ; xsav(is,n)=xx1(i0,j0,n)-360;
0122 savI(is,n)=ii ; savJ(is,n)=jj ;
0123 end
0124
0125 xPmid=(xx1(i1,j1,n) + xsav(is,n))/2;
0126 yPmid=(yy2(i1,j1,n)+yy2(i0,j0,n))/2;
0127 ic=min(nc,ii); jc=min(nc,jj);
0128 xCenter=xcs(ic+(n-1)*nc,jc) ; yCenter=ycs(ic+(n-1)*nc,jc) ;
0129 %- compute vector product p1,p2 x Pmid,Center :
0130 dirUV = (yCenter - yPmid)*(xx1(i1,j1,n) - xsav(is,n)) ...
0131 - (xCenter - xPmid)*(yy2(i1,j1,n)-yy2(i0,j0,n)) ;
0132 if dirUV < 0, savF(is,n)=-savF(is,n) ; end
0133 % if yCenter <= yMid, savF(is,n)=-savF(is,n) ; end
0134 % if ydim == 1 & dirUV*(yCenter - yPmid) <= 0 ,
28aaff1409 Jean*0135 if dirUV*(yCenter - yPmid) < 0 ,
0136 fprintf(['- Sign - : n,is,i0,i1,j0,j1= %i %i %i %i %i %i ', ...
0137 ' ; dirUV= %e \n'], n,is,i0,i1,j0,j1,dirUV);
323aa12fe2 Jean*0138 fprintf(' %8.3f %8.3f %8.3f %8.3f \n', ...
0139 xsav(is,n),xx1(i1,j1,n),xPmid,xCenter);
0140 fprintf(' %8.3f %8.3f %8.3f %8.3f \n', ...
0141 yy2(i0,j0,n),yy2(i1,j1,n),yPmid,yCenter);
0142 end
0143
0144 isav(1+is,n)=i1 ; jsav(1+is,n)=j1 ; xsav(1+is,n)=xx1(i1,j1,n);
0145 nMx6t(n)=1+is ; x1=xx2(i1,j1,n);
0146 i0 = i1 ; j0 = j1 ; x0=x1 ;
0147 end
0148
0149 end
0150
0151 if Ymn < XYout/2,
0152 %- look for the next point : try the 4 direction i-1,i+1,j-1,j+1 :
0153 y4c=XYout*ones(4,1); y0=yy2(i0,j0,n) ;
0154 %if xMid(n) == 1, x0=xx2(i0,j0,n) ; end
0155
0156 if i0 > 1, xx = xx2(i0-1,j0,n);
0157 ddx= xx-x0+3*180 ; ddx=rem(ddx,360)-180 ;
0158 if ddx > 0, y4c(1)=yy2(i0-1,j0,n) ; end ; end
0159 if i0 < ncp, xx = xx2(i0+1,j0,n);
0160 ddx= xx-x0+3*180 ; ddx=rem(ddx,360)-180 ;
0161 if ddx > 0, y4c(2)=yy2(i0+1,j0,n) ; end ; end
0162 if j0 > 1, xx = xx2(i0,j0-1,n);
0163 ddx= xx-x0+3*180 ; ddx=rem(ddx,360)-180 ;
0164 if ddx > 0, y4c(3)=yy2(i0,j0-1,n) ; end ; end
0165 if j0 < ncp, xx = xx2(i0,j0+1,n);
0166 ddx= xx-x0+3*180 ; ddx=rem(ddx,360)-180 ;
0167 if ddx > 0, y4c(4)=yy2(i0,j0+1,n) ; end ; end
0168 %-- avoid the North & South poles (since Long is not well define there):
0169 y4c(find( abs(y4c) == 90)) = XYout;
0170 y4c=abs(y4c - yl); Ymn=min(y4c) ; i4c=find(y4c==Ymn);
0171 end
0172
0173 end; % end(while)
0174
0175 end ; end
0176 return
0177
