ell2cart.m

Set of tools to perform transformations between projection, ellipsoidal and cartesian coordinates in

function CART=ell2cart(ELL,ellips)

% ELL2CART performs transformation from ellipsoidal coordinates to cartesian coordinates on that
%          ellipsoid
% 
% CART=ell2cart(ELL,ellips)
% 
% Also necessary:   Ellipsoids.mat   (see beneath)
% 
% Inputs:   ELL coordinates on ellipsoid as nx3-matrix (longitude, latitude, height) [degree, m]
%               3xn-matrices are allowed. Be careful with 3x3-matrices!
%               nx2-matrices are allowed, all heights are set to 0 in that case.
%               Southern hemisphere is signalled by negative latitude.
%        ellips the underlying ellipsoid as string in lower case letters, default is 'besseldhdn'
%               The ellipsoids are stored in the mat-File "Ellipsoids.mat" also as
%               cell-arrays named by the ellipsoid, e.g. 'bessel1841' or 'wgs84'.
%               The ellipsoid cell array fields are
%                    a  semimajor axis
%                    b  semiminor axis
%                    f  flattening
%                   
% Outputs: CART 3xn-matrix with right-handed cartesian coordinates
% 
% ell2cart is meant as second step for geodetic transformations, when projection coordinates have to
% be transformed either to global coordinates or other projections.
%         
% Ellipsoids.mat are used because several of my files and GUIs need that parameters. 
% I considered it a better way to add those definition files rather to hardcode it
% several times. Feel free to enter your own ellipsoids.
%
% See also: ell2proj, proj2ell, cart2ell, d3trafo, Ellipsoids, Projections

% Author:
% Peter Wasmeier, Technical University of Munich
% [e-mail]
% Jan 18, 2006

% Do some input checking
if ~any(ismember(size(ELL),[2 3]))
    error('Coordinate list ELL must be a nx3- or nx2-matrix!')
elseif (ismember(size(ELL,1),[2 3]))&&(~ismember(size(ELL,2),[2 3]))
    ELL=ELL';
end
if size(ELL,2)==2,ELL(:,3)=zeros(size(ELL,1),1);end
if nargin<2,ellips='besseldhdn';end

% Load ellipsoids
load Ellipsoids;
if ~exist(ellips,'var'), error(['Ellipsoid ',ellips,' is not defined in Ellipsoids.mat - check your definitions!.'])
end
eval(['ell=',ellips,';']);

% Do calculations
CART=zeros(size(ELL));
rho=180/pi;
B=ELL(:,2)/rho;
L=ELL(:,1)/rho;

% 1. numerical eccentricity
e2=(ell.a^2-ell.b^2)/ell.a^2;

% norm radius
N=ell.a./sqrt(1-e2*sin(B).^2);

% cartesian coordinates
CART(:,1)=(N+ELL(:,3)).*cos(B).*cos(L);
CART(:,2)=(N+ELL(:,3)).*cos(B).*sin(L);
CART(:,3)=(N.*(1-e2)+ELL(:,3)).*sin(B);