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

%POLYCON  Polyconic Projection
%
%  For this projection, each parallel has a curvature identical to
%  its curvature on a cone tangent at that latitude.  Since each
%  parallel would have its own cone, this is a "polyconic" projection.
%  Scale is true along the central meridian and along each parallel.
%  This projection is free of distortion only along the central meridian;
%  distortion can be severe at extreme longitudes.  This projection is
%  neither conformal nor equal area.
%
%  This projection was apparently originated about 1820 by Ferdinand
%  Rudolph Hassler.  It is also known as the American Polyconic and
%  the Ordinary Polyconic projections.

%  Copyright 1996-1998 by Systems Planning and Analysis, Inc. and The MathWorks, Inc.
% $Revision: 1.7 $
%  Written by:  E. Byrns, E. Brown


if nargin == 1                  %  Set the default structure entries
	  mstruct.mapparallels = [];
      mstruct.nparallels   = 0;
	  mstruct.fixedorient  = [];
	  mstruct.trimlat = angledim([-90  90],'degrees',mstruct.angleunits);
	  mstruct.trimlon = angledim([-75  75],'degrees',mstruct.angleunits);

%  Longitude limits are set well within the [-90 90] algorithm limit
%  because the algorithm sometimes does not converge when approaching
%  this boundary.

	  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;
aspect  = mstruct.aspect;
geoid   = mstruct.geoid;
radius  = rsphere('rectifying',mstruct.geoid);
origin  = angledim(mstruct.origin,units,'radians');
trimlat = angledim(mstruct.flatlimit,units,'radians');
trimlon = angledim(mstruct.flonlimit,units,'radians');
epsilon = 5*epsm('radians');
scalefactor = mstruct.scalefactor;
falseeasting = mstruct.falseeasting;
falsenorthing = mstruct.falsenorthing;

%  Adjust the origin latitude to the auxiliary sphere

origin(1) = geod2rec(origin(1),mstruct.geoid,'radians');
trimlat   = geod2rec(trimlat  ,mstruct.geoid,'radians');

switch direction
case 'forward'

     lat  = angledim(in1,units,'radians');
     lat  = geod2rec(lat,mstruct.geoid,'radians');
     long = angledim(in2,units,'radians');

%  Back off of the +/- 180 degree points.  The inverse
%  algorithm has trouble distinguishing points at -180 degrees
%  when they are near the pole.

     indx = find( abs(pi - abs(long)) <= epsilon);
     if ~isempty(indx)
	      long(indx) = long(indx) - sign(long(indx))*epsilon;
     end

%  Continue with the transformation

     [lat,long] = rotatem(lat,long,origin,direction);   %  Rotate to new origin
     [lat,long,clipped] = clipdata(lat,long,object);    %  Clip at the date line
     latgeod = rec2geod(lat,mstruct.geoid,'radians');
     [lat,long,trimmed] = trimdata(lat,trimlat,long,trimlon,object);

%  Construct the structure of altered points

     savepts.trimmed = trimmed;
     savepts.clipped = clipped;

%  Back off of the +/- 90 degree points.  This allows
%  the differentiation of longitudes at the poles of the transformed
%  coordinate system.

     indx = find(abs(pi/2 - abs(lat)) <= epsilon);
     if ~isempty(indx)
	    lat(indx) = (pi/2 - epsilon) * sign(lat(indx));
	    latgeod(indx) = (pi/2 - epsilon) * sign(latgeod(indx));
     end

%  Pick up NaN place holders

      x = long;   y = lat;

%  Eliminate singularities in transformations at 0 latitude.

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

%  Points at zero latitude

     if ~isempty(indx1)
          x(indx1) = geoid(1) * long(indx1);
		  y(indx1) = 0;
     end

%  Points at non-zero latitude

	 if ~isempty(indx2)
          N = geoid(1) ./ sqrt(1 - (geoid(2)*sin(latgeod(indx2))).^2);
	      E = long(indx2) .* sin(latgeod(indx2));
          x(indx2) = N .* cot(latgeod(indx2)) .* sin(E);
          y(indx2) = radius*lat(indx2) + ...
		             N .* cot(latgeod(indx2)) .* (1-cos(E));
     end

%  Apply scale factor, false easting, northing

	x = x*scalefactor+falseeasting;
	y = y*scalefactor+falsenorthing;

%  Set the output variables

     switch  aspect
	    case 'normal',         out1 = x;      out2 = y;
	    case 'transverse',	   out1 = y;      out2 = -x;
        otherwise,             error('Unrecognized aspect string')
     end


case 'inverse'

     switch  aspect
	    case 'normal',         x = in1;    y = in2;
	    case 'transverse',	   x = -in2;   y = in1;
        otherwise,             error('Unrecognized aspect string')
     end

%  Apply scale factor, false easting, northing

	x = (x-falseeasting)/scalefactor;
	y = (y-falsenorthing)/scalefactor;

% Inverse projection
%  Eliminate singularities in transformations at 0 latitude.

     indx = find(y == 0);
	 if ~isempty(indx);   y(indx) = epsilon;   end


     A = y / geoid(1);                  B = (x / geoid(1)).^2 + A.^2;
     convergence = 1E-10;             maxsteps = 100;
     steps = 1;                       latnew = A;
     converged = 0;                   e = geoid(2);

     while ~converged & steps <= maxsteps
         steps = steps + 1;
		 latold = latnew;

         C = sqrt(1 - (e*sin(latold)).^2) .* tan(latold);
         Ma = (1 - e^2/4 - 3*e^4/64 - 5*e^6/256) * latold - ...
		      (3*e^2/8 + 3*e^4/32 + 45*e^6/1024) * sin(2*latold) + ...
			  (15*e^4/256 + 45*e^6/1024) * sin(4*latold) - ...
			  (35*e^6/3072) * sin(6*latold);
		 M = Ma * geoid(1);

         Mp = (1 - e^2/4 - 3*e^4/64 - 5*e^6/256) - ...
		      2 * (3*e^2/8 + 3*e^4/32 + 45*e^6/1024) * cos(2*latold) + ...
			  4 * (15*e^4/256 + 45*e^6/1024) * cos(4*latold) - ...
			  6 * (35*e^6/3072) * cos(6*latold);


         num = A.*(C.*Ma + 1) - Ma - 0.5*(Ma.^2 + B).*C;
		 fact1 = (e^2)*sin(2*latold) .* (Ma.^2 + B - 2*A.*Ma) ./ (4*C);
		 fact2 = (A-Ma) .* (C.*Mp - 2./sin(2*latold)) - Mp;
		 dellat = num ./ (fact1 + fact2);

         if max(abs(dellat(:))) <= convergence;     converged = 1;
		      else;                      latnew = latold - dellat;
	     end
     end

     long = asin(x.*C/geoid(1)) ./ sin(latnew);

     latnew(indx) = 0;       %  Correct for y = 0 points
     lat  = geod2rec(latnew,mstruct.geoid,'radians');

%  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);
     lat        = rec2geod(lat,mstruct.geoid,'radians');

%  Transform to desired units

     out1 = angledim(lat, 'radians', units);
     out2 = angledim(long,'radians', 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
