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	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