matlab_class/gcmfaces_calc/gcmfaces_lines_transp.m

function [line_out]=gcmfaces_lines_transp(varargin);
%object:    compute geat circle transport lines, defined 
%           by pair of coordinates, for use in calc_transports
%inputs:    lonPairs (e.g. [-68 -63] by default) is the pair of longitudes
%           latPairs (e.g. [-54 -66] by default) is the pair of latgitudes
%           names (e.g. {'Drake Passage'} by default) is the transport line name
%output:    (optional) line_out is a strcuture with the following fields 
%           lonPair, latPair, name, mskCedge, mskWedge, mskSedge
%           If no output is specified then line_out is copied to 
%           mygrid (global variable) as 'mygrid.LINES_MASKS'

global mygrid;

if nargin>0;
    lonPairs=varargin{1};
    latPairs=varargin{2};
    names=varargin{3};
else;
    lonPairs=[-68 -63];
    latPairs=[-54 -66];
    names={'Drake Passage'};
end;

for iy=1:length(names);
    
    lonPair=lonPairs(iy,:);
    latPair=latPairs(iy,:);
    name=names{iy};
    
    %get carthesian coordinates:
    %... of grid
    lon=mygrid.XC; lat=mygrid.YC;
    x=cos(lat*pi/180).*cos(lon*pi/180);
    y=cos(lat*pi/180).*sin(lon*pi/180);
    z=sin(lat*pi/180);
    %... and of end points
    x0=cos(latPair*pi/180).*cos(lonPair*pi/180);
    y0=cos(latPair*pi/180).*sin(lonPair*pi/180);
    z0=sin(latPair*pi/180);
    
    %get the rotation matrix:
    %1) rotate around x axis to put first point at z=0
    theta=atan2(-z0(1),y0(1));
    R1=[[1;0;0] [0;cos(theta);sin(theta)] [0;-sin(theta);cos(theta)]];
    tmp0=[x0;y0;z0]; tmp1=R1*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:);
    x0=x1; y0=y1; z0=z1;
    %2) rotate around z axis to put first point at y=0
    theta=atan2(x0(1),y0(1));
    R2=[[cos(theta);sin(theta);0] [-sin(theta);cos(theta);0] [0;0;1]];
    tmp0=[x0;y0;z0]; tmp1=R2*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:);
    x0=x1; y0=y1; z0=z1;
    %3) rotate around y axis to put second point at z=0
    theta=atan2(-z0(2),-x0(2));
    R3=[[cos(theta);0;-sin(theta)] [0;1;0] [sin(theta);0;cos(theta)]];
    tmp0=[x0;y0;z0]; tmp1=R3*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:);
    x0=x1; y0=y1; z0=z1;
    
    %apply rotation to grid:
    tmpx=convert2array(x); tmpy=convert2array(y); tmpz=convert2array(z);
    tmp1=find(~isnan(tmpx));
    tmpx2=tmpx(tmp1); tmpy2=tmpy(tmp1); tmpz2=tmpz(tmp1);
    tmp2=[tmpx2';tmpy2';tmpz2'];
    tmp3=R3*R2*R1*tmp2;
    tmpx2=tmp3(1,:); tmpy2=tmp3(2,:); tmpz2=tmp3(3,:);
    tmpx(tmp1)=tmpx2; tmpy(tmp1)=tmpy2; tmpz(tmp1)=tmpz2;
    x=convert2array(tmpx); y=convert2array(tmpy); z=convert2array(tmpz);

    %compute the great circle mask:
    mskCint=1*(z>0);
    [mskCedge,mskWedge,mskSedge]=gcmfaces_edge_mask(mskCint);
        
    %select the shorther segment:
    for kk=1:3;
        %select field to treat:
        switch kk;
            case 1; mm=mskCedge;
            case 2; mm=mskWedge;
            case 3; mm=mskSedge;
        end;
        %split in two segments:
        theta=[];
        theta(1)=atan2(y0(1),x0(1));
        theta(2)=atan2(y0(2),x0(2));
        
        tmpx=convert2array(x); tmpy=convert2array(y); tmpz=convert2array(z);
        tmptheta=atan2(tmpy,tmpx);
        tmpm=convert2array(mm);
        if theta(2)<0;
            tmp00=find(tmptheta<=theta(2)); tmptheta(tmp00)=tmptheta(tmp00)+2*pi;
            theta(2)=theta(2)+2*pi;
        end;
        %select the shorther segment:
        if theta(2)-theta(1)<=pi;
            tmpm(find(tmptheta>theta(2)|tmptheta<theta(1)))=NaN;
        else;
            tmpm(find(tmptheta<=theta(2)&tmptheta>=theta(1)))=NaN;
        end;
        mm=convert2array(tmpm);
        %store result:
        switch kk;
            case 1; mskCedge=mm;
            case 2; mskWedge=mm;
            case 3; mskSedge=mm;
        end;
    end;
    
    %store so-defined line:
    line_cur=struct('lonPair',lonPair,'latPair',latPair,'name',name,...
        'mskCedge',mskCedge,'mskWedge',mskWedge,'mskSedge',mskSedge);
        
    %add to lines:
    if iy==1;
        LINES_MASKS=line_cur;
    else;
        LINES_MASKS(iy)=line_cur;
    end;
    
end;

if nargout==0; %add to mygrid
  mygrid.LINES_MASKS=LINES_MASKS;
else;          %output to line_out 
  line_out=line_cur;
end;


1	function [line_out]=gcmfaces_lines_transp(varargin);
2	%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
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(latpi/180).cos(lon*pi/180);
34	y=cos(latpi/180).sin(lon*pi/180);
35	z=sin(lat*pi/180);
36	%... and of end points
37	x0=cos(latPairpi/180).cos(lonPair*pi/180);
38	y0=cos(latPairpi/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=R3R2R1*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
68	%compute the great circle mask:
69	mskCint=1*(z>0);
70	[mskCedge,mskWedge,mskSedge]=gcmfaces_edge_mask(mskCint);
71
72	%select the shorther segment:
73	for kk=1:3;
74	%select field to treat:
75	switch kk;
76	case 1; mm=mskCedge;
77	case 2; mm=mskWedge;
78	case 3; mm=mskSedge;
79	end;
80	%split in two segments:
81	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	tmpm=convert2array(mm);
88	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	%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	%store result:
100	switch kk;
101	case 1; mskCedge=mm;
102	case 2; mskWedge=mm;
103	case 3; mskSedge=mm;
104	end;
105	end;
106
107	%store so-defined line:
108	line_cur=struct('lonPair',lonPair,'latPair',latPair,'name',name,...
109	'mskCedge',mskCedge,'mskWedge',mskWedge,'mskSedge',mskSedge);
110
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	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
126