Warning, /utils/matlab/cs_grid/mk_psiLine_CS.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 4ec37fd8 on 2023-10-03 22:00:32 UTC
323aa12fe2 Jean*0001 kgr=3; nRgr=3; kwr=1;
0002 dbug=1;
0003
0c7ba451ba Jean*0004 %- set ncdf=1 to load MNC (NetCDF) grid-files ;
0005 % or ncdf=0 to load MDS (binary) grid-files :
0006 ncdf=0;
4ec37fd829 Jean*0007 gDir='grid_files/';
0008 G=load_grid(gDir,10+ncdf);
0009 xcs=G.xC; ycs=G.yC; xcg=G.xG; ycg=G.yG;
0010
0011 %------------
0012 n1h=size(xcs,1); n2h=size(xcs,2);
0013 if n1h == 6*n2h, nc=n2h;
0014 elseif n1h*6 == n2h, nc=n1h;
0015 else
0016 error([' grid var size: ',int2str(n1h),' x ',int2str(n2h),' does not fit regular cube !']);
0017 end
0018 nPg=nc*nc*6; ncp=nc+1; nPp2=nPg+2; n6c=nc*6;
0019 %------------
0020
0021 %--------------------------------
0022 %- when stored in long-vector, use "compact" convention (i.e., 1 face after the other)
0023 if n2h == nc,
0024 xcg=permute(reshape(xcg,[nc 6 nc]),[1 3 2]);
0025 ycg=permute(reshape(ycg,[nc 6 nc]),[1 3 2]);
0026 end
0027 xcg=reshape(xcg,[nPg 1]); ycg=reshape(ycg,[nPg 1]);
0028 %- add the 2 missing corners:
0029 xcg(nPg+1)=xcg(1); ycg(nPg+1)=ycg(1+2*nc*nc);
0030 xcg(nPg+2)=xcg(1+3*nc*nc); ycg(nPg+2)=ycg(1);
0031 fprintf(' add 1rst corner: %i %8.3f %8.3f\n',nPg+1,xcg(nPg+1,1),ycg(nPg+1,1));
0032 fprintf(' add 2nd corner: %i %8.3f %8.3f\n',nPg+2,xcg(nPg+2,1),ycg(nPg+2,1));
0033 %--------------------------------
0034
0035 %- make a "small list" of Z point (= with no double point)
0036 xx1=xcg; yy1=ycg;
323aa12fe2 Jean*0037
0c7ba451ba Jean*0038 xx2=split_Z_cub(xcg);
323aa12fe2 Jean*0039 yy2=split_Z_cub(ycg);
0040 xx3=reshape(xx2,ncp*ncp*6,1); yy3=reshape(yy2,ncp*ncp*6,1);
0041
4ec37fd829 Jean*0042 %- store number of "round" (nRd) and for each round:
0043 % list_N = number of Strean-Function (Psi) to assign
0044 % list_C = list of Current Psi point to assign
0045 % list_P = list of Previously assigned Psi point to add to
0046 % listUV = list of transport (U or V) to use from P to C
0047 % list_T = list to select sign and component (U or V) to use
0048 % tagC = just to keep track of which Psi point has been assigned
0049 %--
0050 list_C=zeros(n6c,n6c); list_P=zeros(n6c,n6c);
0051 listUV=zeros(n6c,n6c); list_T=zeros(n6c,n6c);
0052 list_N=zeros(1,n6c); tagC=zeros(nPp2,1);
323aa12fe2 Jean*0053
0054 %-- start from the N.Pole:
0c7ba451ba Jean*0055 [IJ]=find(yy1==max(yy1));
4ec37fd829 Jean*0056 nRd=1; nbp=length(IJ); list_N(1,nRd)=nbp;
0057 list_C(1:nbp,nRd)=IJ; tagC(IJ(1:nbp),1)=nRd;
323aa12fe2 Jean*0058
0059 doRd=1;
0060 fprintf('nRd= %i ; nbp= %i \n',nRd,nbp);
0061 while doRd > 0,
0062 %- next round :
0063 n2p=0;
4ec37fd829 Jean*0064 list_c=zeros(n6c,1); list_p=zeros(n6c,1);
0065 listuv=zeros(n6c,1); list_t=zeros(n6c,1);
323aa12fe2 Jean*0066
0067 %- for all point of the previous list, find the neighbour point not yet tagged:
0068 for n=1:nbp,
0069 %- find the 4 neighbours:
0c7ba451ba Jean*0070 ij=list_C(n,nRd);
4ec37fd829 Jean*0071 i=1+rem(ij-1,nc); j=fix((ij-1)/nc); k=1+fix(j/nc); j=1+rem(j,nc);
0072 if ij > nPg | i*j==1,
323aa12fe2 Jean*0073 %-- Corner point:
0074 nbv=3;
0075 [I]=find( xx3==xx1(ij) & yy3==yy1(ij) ) ;
0076 if dbug>0, fprintf('Corner: i,j,k= %2i %2i %i ',i,j,k); end
0077 if length(I) ~= 3, fprintf('Pb C: stop\n'); return; end
0078 for l=1:length(I),
0079 i2=I(l);
4ec37fd829 Jean*0080 j2=fix((i2-1)/ncp); k2=1+fix(j2/ncp); j2=1+rem(j2,ncp); i2=1+rem(i2-1,ncp);
323aa12fe2 Jean*0081 if dbug>0, fprintf('; l=%i: %2i %2i %i ',l,i2,j2,k2); end
4ec37fd829 Jean*0082 if ij==nPg+2,
0083 %- 3 times (nc,1) :
0084 lc(l)=nc+(j2-1)*nc+(k2-1)*nc*nc; luv(l)=lc(l); typ(l)=-2;
0085 elseif ij==nPg+1,
0086 %- 3 times (1,nc) :
0087 lc(l)=i2+(nc-1)*nc+(k2-1)*nc*nc; luv(l)=lc(l); typ(l)=1;
323aa12fe2 Jean*0088 elseif i2==ncp & j2==1,
4ec37fd829 Jean*0089 lc(2)=nc+(j2-1)*nc+(k2-1)*nc*nc; luv(2)=lc(2); typ(2)=-2;
323aa12fe2 Jean*0090 elseif i2==1 & j2==ncp,
4ec37fd829 Jean*0091 lc(2)=i2+(nc-1)*nc+(k2-1)*nc*nc; luv(2)=lc(2); typ(2)=1;
323aa12fe2 Jean*0092 end
0093 end
0094 if dbug>0, fprintf('\n'); end
4ec37fd829 Jean*0095 if ij <= nPg,
323aa12fe2 Jean*0096 lc(1)=ij+1; luv(1)=ij; typ(1)=2;
4ec37fd829 Jean*0097 lc(3)=ij+nc; luv(3)=ij; typ(3)=-1;
323aa12fe2 Jean*0098 end
0099 % fprintf('Pb C: stop\n'); return;
0100 else
0101 %-- Not a corner point:
0102 nbv=4;
0103 %--- set 1rst neighbour i+1 :
0104 lc(1)=ij+1; luv(1)=ij; typ(1)=2;
0105 if i==nc,
0106 %- edge point:
0107 [I]=find( xx1==xx2(i+1,j,k) & yy1==yy2(i+1,j,k) ) ;
0108 if length(I) ~= 1, fprintf('Pb A-i: stop\n'); return; end
0109 lc(1)=I(1);
0110 end
0111 %--- set 2nd neighbour i-1 :
0112 if i==1,
0113 %- edge point:
0114 [I]=find( xx3==xx1(ij) & yy3==yy1(ij) ) ;
0115 if dbug>1, fprintf('Edge: i,j,k= %2i %2i %i ',i,j,k); end
0116 if length(I) ~= 2, fprintf('Pb E-i: stop\n'); return; end
0117 for l=1:length(I),
0118 i2=I(l);
4ec37fd829 Jean*0119 j2=fix((i2-1)/ncp); k2=1+fix(j2/ncp); j2=1+rem(j2,ncp); i2=1+rem(i2-1,ncp);
323aa12fe2 Jean*0120 if dbug>1, fprintf('; l= %i : i,j,k_2= %2i %2i %i ',l,i2,j2,k2); end
0c7ba451ba Jean*0121 if k2 ~= k & i2==ncp,
4ec37fd829 Jean*0122 lc(2)=nc+(j2-1)*nc+(k2-1)*nc*nc; luv(2)=lc(2); typ(2)=-2;
323aa12fe2 Jean*0123 elseif k2 ~= k & j2==ncp,
4ec37fd829 Jean*0124 lc(2)=i2+(nc-1)*nc+(k2-1)*nc*nc; luv(2)=lc(2); typ(2)=1;
323aa12fe2 Jean*0125 end
0126 end
0127 if dbug>1, fprintf('\n'); end
0128 else
0129 lc(2)=ij-1; luv(2)=ij-1; typ(2)=-2;
0130 end
0131 %--- set 3rd neighbour j+1 :
4ec37fd829 Jean*0132 lc(3)=ij+nc; luv(3)=ij; typ(3)=-1;
323aa12fe2 Jean*0133 if j==nc,
0134 %- edge point:
0135 [I]=find( xx1==xx2(i,j+1,k) & yy1==yy2(i,j+1,k) ) ;
0136 if length(I) ~= 1, fprintf('Pb A-j: stop\n'); return; end
0137 lc(3)=I(1);
0138 end
0139 %--- set 4th neighbour j-1 :
0140 if j==1,
0141 %- edge point:
0142 [I]=find( xx3==xx1(ij) & yy3==yy1(ij) ) ;
0143 if dbug>1, fprintf('Edge: i,j,k= %2i %2i %i ',i,j,k); end
0144 if length(I) ~= 2, fprintf('Pb E-j: stop\n'); return; end
0145 for l=1:length(I),
0146 i2=I(l);
4ec37fd829 Jean*0147 j2=fix((i2-1)/ncp); k2=1+fix(j2/ncp); j2=1+rem(j2,ncp); i2=1+rem(i2-1,ncp);
323aa12fe2 Jean*0148 if dbug>1, fprintf('; l= %i : i,j,k_2= %2i %2i %i ',l,i2,j2,k2); end
0c7ba451ba Jean*0149 if k2 ~= k & j2==ncp,
4ec37fd829 Jean*0150 lc(4)=i2+(nc-1)*nc+(k2-1)*nc*nc; luv(4)=lc(4); typ(4)=1;
323aa12fe2 Jean*0151 elseif k2 ~= k & i2==ncp,
4ec37fd829 Jean*0152 lc(4)=nc+(j2-1)*nc+(k2-1)*nc*nc; luv(4)=lc(4); typ(4)=-2;
323aa12fe2 Jean*0153 end
0154 end
0155 if dbug>1, fprintf('\n'); end
0156 else
4ec37fd829 Jean*0157 lc(4)=ij-nc; luv(4)=ij-nc; typ(4)=1;
323aa12fe2 Jean*0158 end
0159 end
0160 %- keep the untagged points in a temp list:
0c7ba451ba Jean*0161 for l=1:nbv, if tagC(lc(l))==0,
0162 n2p=n2p+1;
0163 list_p(n2p)=ij; list_c(n2p)=lc(l);
323aa12fe2 Jean*0164 listuv(n2p)=luv(l); list_t(n2p)=typ(l);
0165 end; end;
0166 end
0167 if kgr == 1 & nRd==nRgr,
0168 for n=1:n2p,
0169 xloc=[xx1(list_p(n)) xx1(list_c(n))];
4ec37fd829 Jean*0170 xloc(2)=xloc(1)+rem(xloc(2)-xloc(1)+540,360)-180;
323aa12fe2 Jean*0171 yloc=[yy1(list_p(n)) yy1(list_c(n))];
4ec37fd829 Jean*0172 plot(xloc,yloc,'bx-');
323aa12fe2 Jean*0173 if n==1, hold on ; end
0174 end
0175 hold off ; grid
0176 elseif kgr == 2 & nRd==nRgr,
0177 xloc=xx1(list_c(1:n2p)); yloc=yy1(list_c(1:n2p));
0178 plot(xloc,yloc,'ro'); hold on ;
0179 xloc=xx1(list_p(1:n2p)); yloc=yy1(list_p(1:n2p));
0180 plot(xloc,yloc,'bo');
0c7ba451ba Jean*0181 hold off ; grid ;
323aa12fe2 Jean*0182 end
0183 if n2p == 0, doRd=0; else nRd=nRd+1 ; end
0184 %- deal with non-single points:
0185 nbp=0;
0186 for n=1:n2p,
0187 ij=list_c(n); dbl=0;
0188 if n > 1, [I]=find(list_C(1:nbp,nRd)==ij); dbl=length(I); end
0189 if dbl==0,
0190 nbp=nbp+1;
0191 list_C(nbp,nRd)=list_c(n);
0192 list_P(nbp,nRd)=list_p(n);
0193 listUV(nbp,nRd)=listuv(n);
0194 list_T(nbp,nRd)=list_t(n);
0195 tagC(ij,1)=nRd;
0196 elseif dbl==1,
0197 %- choose between 2 origin points: select the smaller Delta-X
0c7ba451ba Jean*0198 xx=xx1(ij);
323aa12fe2 Jean*0199 dx1=xx1(list_P(I,nRd)); dx1=rem(dx1-xx+540,360)-180;
0200 dx2=xx1(list_p(n)); dx2=rem(dx2-xx+540,360)-180;
0201 if abs(dx2) < abs(dx1),
0202 %- replace I with n :
0203 list_C(I,nRd)=list_c(n);
0204 list_P(I,nRd)=list_p(n);
0205 listUV(I,nRd)=listuv(n);
0206 list_T(I,nRd)=list_t(n);
0207 else
0208 %- remove point n :
0209 list_c(n)=0;
0210 end
0211 else
0212 fprintf('Pb D: stop\n'); return;
0213 end
0214 end
4ec37fd829 Jean*0215 if nbp > 0,
323aa12fe2 Jean*0216 list_N(1,nRd)=nbp;
0217 fprintf('nRd= %i ; nbp= %i \n',nRd,nbp);
0218 end
0219 %if nRd > 64, doRd=0; end
0220 end
0221
0222 if kgr == 3,
0223 for n=2:nRd,
4ec37fd829 Jean*0224 if n == nRgr, linC='r-o'; else linC='b-'; end
323aa12fe2 Jean*0225 for i=1:list_N(1,n),
0226 xloc=[xx1(list_P(i,n)) xx1(list_C(i,n))];
0227 yloc=[yy1(list_P(i,n)) yy1(list_C(i,n))];
0228 plot(xloc,yloc,linC);
0229 if i==1 & n==2, hold on ; end
0230 end
0231 end
0c7ba451ba Jean*0232 hold off ; AA=axis;
0233 AA(1)=max(-180,AA(1)); AA(2)=min(183,AA(2));
0234 AA(3)=max(-90,AA(3)); AA(4)=min(90,AA(4));
323aa12fe2 Jean*0235 axis(AA); grid
0236 elseif kgr == 4,
0237 for n=1:nRd,
0238 if n == nRd, linC='ro'; else linC='ko'; end
0239 nbp=list_N(1,n);
0240 xloc=xx1(list_C(1:nbp,n)); yloc=yy1(list_C(1:nbp,n));
0241 plot(xloc,yloc,linC);
0242 if n==1, hold on ; end
0243 end
0244 xloc=xx1(list_P(1:nbp,nRd)); yloc=yy1(list_P(1:nbp,nRd));
0245 plot(xloc,yloc,'bx');
0c7ba451ba Jean*0246 hold off ; grid ;
323aa12fe2 Jean*0247 end
0248
0249 if kwr==1,
0250 %-- copy to the right sixe arrays:
0251 psiNx=max(list_N);
0c7ba451ba Jean*0252 psiNy=nRd;
4ec37fd829 Jean*0253 %- store number of pass (nRd) and for each pass:
0254 % psi_N = number of Strean-Function (Psi) to assign
0255 % psi_C = list of Current Psi point to assign
0256 % psi_P = list of Previously assigned Psi point to add to
0257 % psiUV = list of transport (U or V) to use from P to C
0258 % psi_T = list to select sign and component (U or V) to use
0259 %--
323aa12fe2 Jean*0260 psi_N=zeros(1,psiNy);
0261 psi_C=zeros(psiNx,psiNy); psi_P=zeros(psiNx,psiNy);
0262 psiUV=zeros(psiNx,psiNy); psi_T=zeros(psiNx,psiNy);
0263 psi_N=list_N(1,1:psiNy);
0264 psi_C=list_C(1:psiNx,1:psiNy);
0265 psi_P=list_P(1:psiNx,1:psiNy);
0266 psiUV=listUV(1:psiNx,1:psiNy);
0267 psi_T=list_T(1:psiNx,1:psiNy);
0268 %- write to a file :
0269 namf=['psiLine_N2S_cs',int2str(nc)];
0270 fprintf(['write on file : ',namf,'.mat ']);
0271 save(namf,'psi_N','psi_C','psi_P','psiUV','psi_T');
0272 fprintf(' <== O.K. \n');
0273 end
