matlab_class/gcmfaces_maps/gcmfaces_stereoproj.m

function [xx,yy]=gcmfaces_stereoproj(XC0,YC0,XC,YC);
%[xx,yy]=gcmfaces_stereoproj(XC0,YC0);
%object:    compute stereographic projection putting XC0,YC0 at 0,0
%inputs:    XC0,YC0 are the origin longitude,latitude
%(optional) XC,YC are the lon,lat points to project (XC,YC by default)
%outputs:   xx,yy are the projected points

%for additional information see :
% http://www4.ncsu.edu/~franzen/public_html/CH795Z/math/lab_frame/lab_frame.html
% http://physics.unm.edu/Courses/Finley/p503/handouts/SphereProjectionFINALwFig.pdf

gcmfaces_global;

if isempty(who('XC'))|isempty(who('YC'));
    XC=mygrid.XC; YC=mygrid.YC;
end;

%compute spherical coordinates:
phi=pi/180*XC; theta=pi/180*(90-YC);
phi0=pi/180*XC0; theta0=pi/180*(90-YC0);

%compute cartesian coordinates:
X=sin(theta).*cos(phi);
Y=sin(theta).*sin(phi);
Z=cos(theta);

x=X; y=Y; z=Z;

%bring chosen point to the north pole:
xx=x; yy=y; zz=z;
x=cos(phi0)*xx+sin(phi0)*yy;
y=-sin(phi0)*xx+cos(phi0)*yy;
z=zz;
%
xx=x; yy=y; zz=z;
x=cos(theta0)*xx-sin(theta0)*zz;
y=yy;
z=sin(theta0)*xx+cos(theta0)*zz;

%stereographic projection from the south pole:
xx=x./(1+z);
yy=y./(1+z);

% nrm=sqrt(xx.^2+yy.^2); 
% msk=1+0*nrm; msk(nrm>tan(pi/4/2))=NaN;%mask points outside of pi/4 cone

1	gforget	1.1	function [xx,yy]=gcmfaces_stereoproj(XC0,YC0,XC,YC);
2			%[xx,yy]=gcmfaces_stereoproj(XC0,YC0);
3			%object: compute stereographic projection putting XC0,YC0 at 0,0
4			%inputs: XC0,YC0 are the origin longitude,latitude
5			%(optional) XC,YC are the lon,lat points to project (XC,YC by default)
6			%outputs: xx,yy are the projected points
7
8			%for additional information see :
9			% http://www4.ncsu.edu/~franzen/public_html/CH795Z/math/lab_frame/lab_frame.html
10			% http://physics.unm.edu/Courses/Finley/p503/handouts/SphereProjectionFINALwFig.pdf
11
12			gcmfaces_global;
13
14			if isempty(who('XC'))\|isempty(who('YC'));
15			XC=mygrid.XC; YC=mygrid.YC;
16			end;
17
18			%compute spherical coordinates:
19			phi=pi/180XC; theta=pi/180(90-YC);
20			phi0=pi/180XC0; theta0=pi/180(90-YC0);
21
22			%compute cartesian coordinates:
23			X=sin(theta).*cos(phi);
24			Y=sin(theta).*sin(phi);
25			Z=cos(theta);
26
27			x=X; y=Y; z=Z;
28
29			%bring chosen point to the north pole:
30			xx=x; yy=y; zz=z;
31			x=cos(phi0)xx+sin(phi0)yy;
32			y=-sin(phi0)xx+cos(phi0)yy;
33			z=zz;
34			%
35			xx=x; yy=y; zz=z;
36			x=cos(theta0)xx-sin(theta0)zz;
37			y=yy;
38			z=sin(theta0)xx+cos(theta0)zz;
39
40			%stereographic projection from the south pole:
41			xx=x./(1+z);
42			yy=y./(1+z);
43
44			% nrm=sqrt(xx.^2+yy.^2);
45			% msk=1+0*nrm; msk(nrm>tan(pi/4/2))=NaN;%mask points outside of pi/4 cone
46