/[MITgcm]/MITgcm_contrib/gael/matlab_class/gcmfaces_calc/gcmfaces_lines_transp.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/matlab_class/gcmfaces_calc/gcmfaces_lines_transp.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.1 - (hide annotations) (download)
Fri Jun 24 17:02:05 2011 UTC (14 years, 1 month ago) by gforget
Branch: MAIN 
- add headers to gecmfaces_calc routines
- move line_greatC_TUV_mask.m, line_zonal_TUV_MASKS.m, line_zonal_TUV_mask.m to gcmfaces_legacy
- add replacements gcmfaces_calc/gcmfaces_lines_transp.m and gcmfaces_calc/gcmfaces_lines_zonal.m
- now the zonal and transport 'lines' masks are added to mygrid as LATS_MASKS and LINES_MASKS
- accordingly, the old LATS_MASKS is removed from arguments list in calc_overturn.m etc.
  (calc_transports and basic_diags_display_transport still have LINES_MASKS args. for now)
  and basic_diags_compute_v3_or_v4.m etc. now call gcmfaces_lines_zonal and gcmfaces_lines_transp.
- line_greatC_TUV_MASKS_v3.m and line_greatC_TUV_MASKS_v4.m now are in sample_analysis
  and they simply output the 'lines' definitions (rather than the masks off gcmfaces_lines_transp)
1 gforget 1.1 function [line_out]=gcmfaces_lines_transp(varargin);
2    
3     global mygrid;
4    
5     if nargin>0;
6     lonPairs=varargin{1};
7     latPairs=varargin{2};
8     names=varargin{3};
9     else;
10     lonPairs=[-68 -63];
11     latPairs=[-54 -66];
12     names={'Drake Passage'};
13     end;
14    
15     for iy=1:length(names);
16    
17     lonPair=lonPairs(iy,:);
18     latPair=latPairs(iy,:);
19     name=names{iy};
20    
21     %get carthesian coordinates:
22     %... of grid
23     lon=mygrid.XC; lat=mygrid.YC;
24     x=cos(lat*pi/180).*cos(lon*pi/180);
25     y=cos(lat*pi/180).*sin(lon*pi/180);
26     z=sin(lat*pi/180);
27     %... and of end points
28     x0=cos(latPair*pi/180).*cos(lonPair*pi/180);
29     y0=cos(latPair*pi/180).*sin(lonPair*pi/180);
30     z0=sin(latPair*pi/180);
31    
32     %get the rotation matrix:
33     %1) rotate around x axis to put first point at z=0
34     theta=atan2(-z0(1),y0(1));
35     R1=[[1;0;0] [0;cos(theta);sin(theta)] [0;-sin(theta);cos(theta)]];
36     tmp0=[x0;y0;z0]; tmp1=R1*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:);
37     x0=x1; y0=y1; z0=z1;
38     %2) rotate around z axis to put first point at y=0
39     theta=atan2(x0(1),y0(1));
40     R2=[[cos(theta);sin(theta);0] [-sin(theta);cos(theta);0] [0;0;1]];
41     tmp0=[x0;y0;z0]; tmp1=R2*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:);
42     x0=x1; y0=y1; z0=z1;
43     %3) rotate around y axis to put second point at z=0
44     theta=atan2(-z0(2),-x0(2));
45     R3=[[cos(theta);0;-sin(theta)] [0;1;0] [sin(theta);0;cos(theta)]];
46     tmp0=[x0;y0;z0]; tmp1=R3*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:);
47     x0=x1; y0=y1; z0=z1;
48    
49     %apply rotation to grid:
50     tmpx=convert2array(x); tmpy=convert2array(y); tmpz=convert2array(z);
51     tmp1=find(~isnan(tmpx));
52     tmpx2=tmpx(tmp1); tmpy2=tmpy(tmp1); tmpz2=tmpz(tmp1);
53     tmp2=[tmpx2';tmpy2';tmpz2'];
54     tmp3=R3*R2*R1*tmp2;
55     tmpx2=tmp3(1,:); tmpy2=tmp3(2,:); tmpz2=tmp3(3,:);
56     tmpx(tmp1)=tmpx2; tmpy(tmp1)=tmpy2; tmpz(tmp1)=tmpz2;
57     x=convert2array(tmpx); y=convert2array(tmpy); z=convert2array(tmpz);
58    
59     %compute the great circle mask:
60     zz=exch_T_N(z);
61     mm=zz<=0; mm(isnan(mm))=0;
62     for iF=1:y.nFaces;
63     eval(['tmp1=mm.f' num2str(iF) ';']);
64     tmp2=circshift(tmp1,[-1 0])+circshift(tmp1,[1 0])+...
65     circshift(tmp1,[0 -1])+circshift(tmp1,[0 1]);
66     tmp2(tmp1>0)=0; tmp2(tmp2~=0)=1; tmp2(tmp2==0)=NaN;
67     eval(['mm.f' num2str(iF) '=tmp2(2:end-1,2:end-1);']);
68     end;
69    
70     %UV part:
71     mmt=mm;
72     nn=exch_T_N(mmt);
73     mm=zz<=0; mm(isnan(mm))=0;
74     mmu=mmt; mmv=mmt;
75     for iFace=1:y.nFaces;
76     iF=num2str(iFace);
77     eval(['tmp_u=mm.f' iF '(1:end-1,:)==1&nn.f' iF '(2:end,:)==1;']);
78     eval(['tmp_u2=mm.f' iF '(2:end,:)==1&nn.f' iF '(1:end-1,:)==1;']);
79     tmp_u=1*(tmp_u(1:end-1,2:end-1)-tmp_u2(1:end-1,2:end-1)); eval(['mmu.f' iF '=tmp_u;']);
80     eval(['tmp_v=mm.f' iF '(:,1:end-1)==1&nn.f' iF '(:,2:end)==1;']);
81     eval(['tmp_v2=mm.f' iF '(:,2:end)==1&nn.f' iF '(:,1:end-1)==1;']);
82     tmp_v=1*(tmp_v(2:end-1,1:end-1)-tmp_v2(2:end-1,1:end-1)); eval(['mmv.f' iF '=tmp_v;']);
83     end;
84     mmu(mmu==0)=NaN; mmv(mmv==0)=NaN;
85    
86     for kk=1:3;
87     %select field to treat:
88     switch kk;
89     case 1; mm=mmt;
90     case 2; mm=mmu;
91     case 3; mm=mmv;
92     end;
93     %split in two contours:
94     theta=[];
95     theta(1)=atan2(y0(1),x0(1));
96     theta(2)=atan2(y0(2),x0(2));
97    
98     tmpx=convert2array(x); tmpy=convert2array(y); tmpz=convert2array(z);
99     tmptheta=atan2(tmpy,tmpx);
100     tmpm=convert2array(mm); tmpmIn=tmpm; tmpmOut=tmpmIn;
101     if theta(2)<0;
102     tmp00=find(tmptheta<=theta(2)); tmptheta(tmp00)=tmptheta(tmp00)+2*pi;
103     theta(2)=theta(2)+2*pi;
104     end;
105     tmpmIn(find(tmptheta>theta(2)|tmptheta<theta(1)))=NaN;
106     tmpmOut(find(tmptheta<=theta(2)&tmptheta>=theta(1)))=NaN;
107     mmIn=convert2array(tmpmIn); mmOut=convert2array(tmpmOut);
108     %ensure that we select the shorther segment as mmIn:
109     if theta(2)-theta(1)>pi; tmp1=mmIn; mmIn=mmOut; mmOut=tmp1; end;
110     %store result:
111     switch kk;
112     case 1; mmtIn=mmIn; mmtOut=mmOut;
113     case 2; mmuIn=mmIn; mmuOut=mmOut;
114     case 3; mmvIn=mmIn; mmvOut=mmOut;
115     end;
116     end;
117    
118     %store so-defined line:
119     line_cur=struct('mmt',mmt,'mmtIn',mmtIn,'mmtOut',mmtOut,...
120     'mmu',mmu,'mmuIn',mmuIn,'mmuOut',mmuOut,...
121     'mmv',mmv,'mmvIn',mmvIn,'mmvOut',mmvOut,...
122     'lonPair',lonPair,'latPair',latPair,'name',name);
123    
124     if 0;
125     iFace=4;
126     figure;
127     tmp1=mygrid.hFacW{iFace}(:,:,1); tmp2=line_cur.mmuIn{iFace}; tmp1(~isnan(tmp2))=-1;
128     subplot(2,2,1); imagesc(tmp1);
129     tmp1=mygrid.hFacS{iFace}(:,:,1); tmp2=line_cur.mmvIn{iFace}; tmp1(~isnan(tmp2))=-1;
130     subplot(2,2,2); imagesc(tmp1);
131    
132     subplot(2,1,2);
133     [X,Y,MSK]=convert2pcol(mygrid.XC,mygrid.YC,mygrid.Depth); MSK(MSK==0)=NaN;
134     pcolor(X,Y,MSK); axis([-180 180 -90 90]); shading flat; caxis([0 6000]);
135    
136     hold on; plot(lonPair(1),latPair(1),'r.','MarkerSize',32);
137     plot(lonPair(2),latPair(2),'k.','MarkerSize',32);
138    
139     [X,Y,mmtIn]=convert2pcol(mygrid.XC,mygrid.YC,line_cur.mmtIn);
140     kk=find(mmtIn==1); plot(X(kk),Y(kk),'m.');
141     % [X,Y,mmtOut]=convert2pcol(mygrid.XC,mygrid.YC,line_cur.mmtOut);
142     % kk=find(mmtOut==1); plot(X(kk),Y(kk),'cx');
143     end;
144    
145     %add to lines:
146     if iy==1;
147     LINES_MASKS=line_cur;
148     else;
149     LINES_MASKS(iy)=line_cur;
150     end;
151    
152     end;
153    
154     mygrid.LINES_MASKS=LINES_MASKS;
155    
156    
157    
  ViewVC Help
Powered by ViewVC 1.1.22