# Documentation of stereo2

Global Index (all files) (short | long) | Local Index (files in subdir) (short | long)

## Function Synopsis

`[out1,out2,savepts] = stereo(mstruct,in1,in2,object,direction,savepts)`

## Help text

```STEREO  Stereographic Azimuthal Projection

This is a perspective projection on a plane tangent at the center
point from the point antipodal to the center point.  The center
point is a pole in the common polar aspect, but it can be any
point.  This projection has two significant properties.  It is
conformal, being free from angular distortion.  Additionally, all
great and small circles are either straight lines or circular
arcs on this projection.  Scale is true only at the center point,
and is constant along any circle having the center point as its
center.  This projection is not equal area.

The polar aspect of this projection appears to have been developed
by the Egyptians and Greeks by the second century B.C.
```

## Listing of function stereo2

```function [out1,out2,savepts] = stereo(mstruct,in1,in2,object,direction,savepts)

%  Copyright (c) 1996-98 by Systems Planning and Analysis, Inc. and The MathWorks, Inc.
% \$Revision: 1.5 \$
%  Written by:  E. Byrns, E. Brown

if nargin == 1                  %  Set the default structure entries
mstruct.mapparallels = [];
mstruct.nparallels   = 0;
mstruct.fixedorient  = [];
mstruct.trimlat = angledim([-inf 90],'degrees',mstruct.angleunits);
mstruct.trimlon = angledim([-180 180],'degrees',mstruct.angleunits);
out1 = mstruct;     return
elseif nargin ~= 5 & nargin ~= 6
error('Incorrect number of arguments')
end

%  Extract the necessary projection data and convert to radians

units   = mstruct.angleunits;
a       = mstruct.geoid(1);
e       = mstruct.geoid(2);

%  Compute projection parameters

a = mstruct.geoid(1);    e = mstruct.geoid(2);

den1  = (1 - (e*sin(origin(1)))^2);
m1    = cos(origin(1)) / sqrt(den1);
if abs(pi/2 - abs(origin(1))) <= epsilon
fact1 = 2*a/sqrt(den1);
else
fact1 = 2*a*m1 / cos(latcnf);
end

%  Adjust the origin latitude to the auxiliary sphere

origin(1)  = latcnf;

switch direction
case 'forward'

%  Back off of the +/- 180 degree points.  This allows
%  the differentiation of longitudes at the dateline in
%  the inverse calculation.  For example, if not performed, -180 degree
%  points may become 180 degrees in the inverse calculation.

indx = find( abs(pi - abs(long)) <= epsilon);

if ~isempty(indx)
long(indx) = long(indx) - sign(long(indx))*epsilon;
end

%  Calculate the azimuths and distances on the rotated auxiliary sphere

[lat,long] = rotatem(lat,long,origin,direction);   %  Rotate to new origin

zeromat = zeros(size(lat));

%  Trim data exceeding a specified range from the origin

[rng,az,trimmed] = trimdata(rng,[-inf max(trimlat)],...
az,[-inf inf],object);
if size(trimmed,2) == 3
indx = trimmed(:,1);
trimmed(:,[2:3]) = [lat(indx) long(indx)];
elseif size(trimmed,2) == 4
indx = trimmed(:,1) + (trimmed(:,2) - 1) * size(lat,1);
trimmed(:,[3:4]) = [lat(indx) long(indx)];
end

%  Construct the structure of altered points

savepts.trimmed = trimmed;
savepts.clipped = [];

%  Perform the projection calculations

A = fact1 ./ ( 1 + cos(rng) );

x = A .* sin(rng) .* sin(az);
y = A .* sin(rng) .* cos(az);

%  Set the output variables

out1 = x;   out2 = y;

case 'inverse'

x = in1;   y = in2;

rho = (x.^2 + y.^2) / fact1^2;

indx1 = find(rho == 0);
indx2 = find(rho ~= 0);

lat = x;    long = y;

if ~isempty(indx1);   lat(indx1)  = 0;    long(indx1) = 0;   end

if ~isempty(indx2)
az = atan2(x(indx2),y(indx2));
rng = acos( (1-rho(indx2)) ./ (1+rho(indx2)));

zeromat = zeros(size(lat(indx2)));
end

%  Undo trims and clips

[lat,long] = undotrim(lat,long,savepts.trimmed,object);
[lat,long] = undoclip(lat,long,savepts.clipped,object);

%  Rotate to Greenwich and transform to desired units

[lat,long] = rotatem(lat,long,origin,direction);

out1 = angledim(lat, 'radians', units);   %  Transform units

otherwise
error('Unrecognized direction string')
end

%  Some operations on NaNs produce NaN + NaNi.  However operations
%  outside the map may product complex results and we don't want
%  to destroy this indicator.

if isieee == 1
indx = find(isnan(out1) | isnan(out2));
out1(indx) = NaN;   out2(indx) = NaN;
end
```