| 1 |
gforget |
1.1 |
function [line_out]=gcmfaces_lines_transp(varargin); |
| 2 |
gforget |
1.4 |
%object: compute geat circle transport lines, defined |
| 3 |
|
|
% by pair of coordinates, for use in calc_transports |
| 4 |
|
|
%inputs: lonPairs (e.g. [-68 -63] by default) is the pair of longitudes |
| 5 |
|
|
% latPairs (e.g. [-54 -66] by default) is the pair of latgitudes |
| 6 |
|
|
% names (e.g. {'Drake Passage'} by default) is the transport line name |
| 7 |
|
|
%output: (optional) line_out is a strcuture with the following fields |
| 8 |
|
|
% lonPair, latPair, name, mskCedge, mskWedge, mskSedge |
| 9 |
|
|
% If no output is specified then line_out is copied to |
| 10 |
|
|
% mygrid (global variable) as 'mygrid.LINES_MASKS' |
| 11 |
gforget |
1.1 |
|
| 12 |
|
|
global mygrid; |
| 13 |
|
|
|
| 14 |
|
|
if nargin>0; |
| 15 |
|
|
lonPairs=varargin{1}; |
| 16 |
|
|
latPairs=varargin{2}; |
| 17 |
|
|
names=varargin{3}; |
| 18 |
|
|
else; |
| 19 |
|
|
lonPairs=[-68 -63]; |
| 20 |
|
|
latPairs=[-54 -66]; |
| 21 |
|
|
names={'Drake Passage'}; |
| 22 |
|
|
end; |
| 23 |
|
|
|
| 24 |
|
|
for iy=1:length(names); |
| 25 |
|
|
|
| 26 |
|
|
lonPair=lonPairs(iy,:); |
| 27 |
|
|
latPair=latPairs(iy,:); |
| 28 |
|
|
name=names{iy}; |
| 29 |
|
|
|
| 30 |
|
|
%get carthesian coordinates: |
| 31 |
|
|
%... of grid |
| 32 |
|
|
lon=mygrid.XC; lat=mygrid.YC; |
| 33 |
|
|
x=cos(lat*pi/180).*cos(lon*pi/180); |
| 34 |
|
|
y=cos(lat*pi/180).*sin(lon*pi/180); |
| 35 |
|
|
z=sin(lat*pi/180); |
| 36 |
|
|
%... and of end points |
| 37 |
|
|
x0=cos(latPair*pi/180).*cos(lonPair*pi/180); |
| 38 |
|
|
y0=cos(latPair*pi/180).*sin(lonPair*pi/180); |
| 39 |
|
|
z0=sin(latPair*pi/180); |
| 40 |
|
|
|
| 41 |
|
|
%get the rotation matrix: |
| 42 |
|
|
%1) rotate around x axis to put first point at z=0 |
| 43 |
|
|
theta=atan2(-z0(1),y0(1)); |
| 44 |
|
|
R1=[[1;0;0] [0;cos(theta);sin(theta)] [0;-sin(theta);cos(theta)]]; |
| 45 |
|
|
tmp0=[x0;y0;z0]; tmp1=R1*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:); |
| 46 |
|
|
x0=x1; y0=y1; z0=z1; |
| 47 |
|
|
%2) rotate around z axis to put first point at y=0 |
| 48 |
|
|
theta=atan2(x0(1),y0(1)); |
| 49 |
|
|
R2=[[cos(theta);sin(theta);0] [-sin(theta);cos(theta);0] [0;0;1]]; |
| 50 |
|
|
tmp0=[x0;y0;z0]; tmp1=R2*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:); |
| 51 |
|
|
x0=x1; y0=y1; z0=z1; |
| 52 |
|
|
%3) rotate around y axis to put second point at z=0 |
| 53 |
|
|
theta=atan2(-z0(2),-x0(2)); |
| 54 |
|
|
R3=[[cos(theta);0;-sin(theta)] [0;1;0] [sin(theta);0;cos(theta)]]; |
| 55 |
|
|
tmp0=[x0;y0;z0]; tmp1=R3*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:); |
| 56 |
|
|
x0=x1; y0=y1; z0=z1; |
| 57 |
|
|
|
| 58 |
|
|
%apply rotation to grid: |
| 59 |
|
|
tmpx=convert2array(x); tmpy=convert2array(y); tmpz=convert2array(z); |
| 60 |
|
|
tmp1=find(~isnan(tmpx)); |
| 61 |
|
|
tmpx2=tmpx(tmp1); tmpy2=tmpy(tmp1); tmpz2=tmpz(tmp1); |
| 62 |
|
|
tmp2=[tmpx2';tmpy2';tmpz2']; |
| 63 |
|
|
tmp3=R3*R2*R1*tmp2; |
| 64 |
|
|
tmpx2=tmp3(1,:); tmpy2=tmp3(2,:); tmpz2=tmp3(3,:); |
| 65 |
|
|
tmpx(tmp1)=tmpx2; tmpy(tmp1)=tmpy2; tmpz(tmp1)=tmpz2; |
| 66 |
|
|
x=convert2array(tmpx); y=convert2array(tmpy); z=convert2array(tmpz); |
| 67 |
gforget |
1.2 |
|
| 68 |
gforget |
1.1 |
%compute the great circle mask: |
| 69 |
gforget |
1.2 |
mskCint=1*(z>0); |
| 70 |
|
|
[mskCedge,mskWedge,mskSedge]=gcmfaces_edge_mask(mskCint); |
| 71 |
|
|
|
| 72 |
|
|
%select the shorther segment: |
| 73 |
gforget |
1.1 |
for kk=1:3; |
| 74 |
|
|
%select field to treat: |
| 75 |
|
|
switch kk; |
| 76 |
gforget |
1.2 |
case 1; mm=mskCedge; |
| 77 |
|
|
case 2; mm=mskWedge; |
| 78 |
|
|
case 3; mm=mskSedge; |
| 79 |
gforget |
1.1 |
end; |
| 80 |
gforget |
1.2 |
%split in two segments: |
| 81 |
gforget |
1.1 |
theta=[]; |
| 82 |
|
|
theta(1)=atan2(y0(1),x0(1)); |
| 83 |
|
|
theta(2)=atan2(y0(2),x0(2)); |
| 84 |
|
|
|
| 85 |
|
|
tmpx=convert2array(x); tmpy=convert2array(y); tmpz=convert2array(z); |
| 86 |
|
|
tmptheta=atan2(tmpy,tmpx); |
| 87 |
gforget |
1.2 |
tmpm=convert2array(mm); |
| 88 |
gforget |
1.1 |
if theta(2)<0; |
| 89 |
|
|
tmp00=find(tmptheta<=theta(2)); tmptheta(tmp00)=tmptheta(tmp00)+2*pi; |
| 90 |
|
|
theta(2)=theta(2)+2*pi; |
| 91 |
|
|
end; |
| 92 |
gforget |
1.2 |
%select the shorther segment: |
| 93 |
|
|
if theta(2)-theta(1)<=pi; |
| 94 |
|
|
tmpm(find(tmptheta>theta(2)|tmptheta<theta(1)))=NaN; |
| 95 |
|
|
else; |
| 96 |
|
|
tmpm(find(tmptheta<=theta(2)&tmptheta>=theta(1)))=NaN; |
| 97 |
|
|
end; |
| 98 |
|
|
mm=convert2array(tmpm); |
| 99 |
gforget |
1.1 |
%store result: |
| 100 |
|
|
switch kk; |
| 101 |
gforget |
1.2 |
case 1; mskCedge=mm; |
| 102 |
|
|
case 2; mskWedge=mm; |
| 103 |
|
|
case 3; mskSedge=mm; |
| 104 |
gforget |
1.1 |
end; |
| 105 |
|
|
end; |
| 106 |
|
|
|
| 107 |
|
|
%store so-defined line: |
| 108 |
gforget |
1.2 |
line_cur=struct('lonPair',lonPair,'latPair',latPair,'name',name,... |
| 109 |
|
|
'mskCedge',mskCedge,'mskWedge',mskWedge,'mskSedge',mskSedge); |
| 110 |
gforget |
1.1 |
|
| 111 |
|
|
%add to lines: |
| 112 |
|
|
if iy==1; |
| 113 |
|
|
LINES_MASKS=line_cur; |
| 114 |
|
|
else; |
| 115 |
|
|
LINES_MASKS(iy)=line_cur; |
| 116 |
|
|
end; |
| 117 |
|
|
|
| 118 |
|
|
end; |
| 119 |
|
|
|
| 120 |
gforget |
1.3 |
if nargout==0; %add to mygrid |
| 121 |
|
|
mygrid.LINES_MASKS=LINES_MASKS; |
| 122 |
|
|
else; %output to line_out |
| 123 |
|
|
line_out=line_cur; |
| 124 |
|
|
end; |
| 125 |
gforget |
1.1 |
|
| 126 |
|
|
|
Parent Directory
| 
Revision Log
| 
Revision Graph