high_res_cube/matlab-grid-generator/map_xyz2lonlat.m

function [lon,lat]=map_xyz2lonlat(x,y,z)
% Convert 3-D coordinates [x,y,z] to [lon,lat]
%
% Assumes "lat" is positive with "z", equatorial plane
% falls at z=0  and "lon" is measured anti-clockwise (eastward)
% from x-axis (y=0) about z-axis.

% Latitude
Req=sqrt( x.^2 + y.^2 );
Z=z;
ii=find(Req == 0);
Z(ii)=Z(ii)*Inf;
Req(ii)=1;
lat=atan( Z ./ Req );

% Longitude
X=x;
Y=y;
ii=find(x == 0);
Y(ii)=Inf;
X(ii)=1;
lon=atan( Y ./ X );

ii=find(x<0 & y>=0);
lon(ii)=pi+lon(ii);
ii=find(x<=0 & y<0);
lon(ii)=-pi+lon(ii);
1	cnh	1.1	function [lon,lat]=map_xyz2lonlat(x,y,z)
2			% Convert 3-D coordinates [x,y,z] to [lon,lat]
3			%
4			% Assumes "lat" is positive with "z", equatorial plane
5			% falls at z=0 and "lon" is measured anti-clockwise (eastward)
6			% from x-axis (y=0) about z-axis.
7
8			% Latitude
9			Req=sqrt( x.^2 + y.^2 );
10			Z=z;
11			ii=find(Req == 0);
12			Z(ii)=Z(ii)*Inf;
13			Req(ii)=1;
14			lat=atan( Z ./ Req );
15
16			% Longitude
17			X=x;
18			Y=y;
19			ii=find(x == 0);
20			Y(ii)=Inf;
21			X(ii)=1;
22			lon=atan( Y ./ X );
23
24			ii=find(x<0 & y>=0);
25			lon(ii)=pi+lon(ii);
26			ii=find(x<=0 & y<0);
27			lon(ii)=-pi+lon(ii);