function varargout = polyconstd(varargin)
%POLYCONSTD  Polyconic Projection -- Standard
%
%  POLYCONSTD  Implements the Polyconic projection directly on a
%  reference ellipsoid, consistent with the industry-standard definition
%  of this projection.  See POLYCON for an alternative implementation
%  based on rotating the rectifying sphere. 
%
%  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 2006-2021 The MathWorks, Inc.

mproj.default = @polyconstdDefault;
mproj.forward = @polyconstdFwd;
mproj.inverse = @polyconstdInv;
mproj.auxiliaryLatitudeType = 'rectifying';
mproj.classCode = 'Pstd';

varargout = applyProjection(mproj, varargin{:});

%--------------------------------------------------------------------------

function mstruct = polyconstdDefault(mstruct)

mstruct.mapparallels = [];
mstruct.nparallels   = 0;
mstruct.fixedorient  = [];
[mstruct.trimlat, mstruct.trimlon] ...
          = fromDegrees(mstruct.angleunits, [-90 90], [-75 75]);

%--------------------------------------------------------------------------

function [x, y] = polyconstdFwd(mstruct, mu, lambda)

[a, e, M0, tolerance, radius] = deriveParameters(mstruct);

% 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(lambda)) <= tolerance);
if ~isempty(indx)
    lambda(indx) = lambda(indx) - sign(lambda(indx))*tolerance;
end

% Convert back to geodetic latitude -- both types are needed

phi = convertlat([a e], mu, 'rectifying', 'geodetic', 'nocheck');

% 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(mu)) <= tolerance);
if ~isempty(indx)
    mu(indx)  = (pi/2 - tolerance) * sign(mu(indx));
    phi(indx) = (pi/2 - tolerance) * sign(phi(indx));
end

% Pick up NaN place holders

x = lambda;
y = mu;

% Eliminate singularities in transformations at 0 latitude.

q = (phi == 0);
x(q) = a * lambda(q);
y(q) = -M0;

% Points at non-zero latitude

N = a ./ sqrt(1 - (e*sin(phi(~q))).^2);
E = lambda(~q) .* sin(phi(~q));
x(~q) = N .* cot(phi(~q)) .* sin(E);
y(~q) = radius*mu(~q) - M0 + N .* cot(phi(~q)) .* (1-cos(E));

%--------------------------------------------------------------------------

function [mu, lambda] = polyconstdInv(mstruct, x, y)

[a, e, M0, tolerance] = deriveParameters(mstruct);

%  Eliminate singularities in transformations at 0 latitude.

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

A = (M0 + y) / a;
B = (x / a).^2 + A.^2;
convergence = 1E-10;
maxsteps = 100;
steps = 1;
phiNew = A;
converged = 0;

while ~converged && (steps <= maxsteps)
    steps = steps + 1;
    phiOld = phiNew;

    C = sqrt(1 - (e*sin(phiOld)).^2) .* tan(phiOld);

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

    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*phiOld) + ...
        4 * (15*e^4/256 + 45*e^6/1024) * cos(4*phiOld) - ...
        6 * (35*e^6/3072) * cos(6*phiOld);

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

    if max(abs(deltaPhi(:))) <= convergence
        converged = 1;
    else
        phiNew = phiOld - deltaPhi;
    end
end

lambda = asin(x.*C/a) ./ sin(phiNew);

phiNew(indx) = 0;       %  Correct for y = 0 points

% Convert to rectifying latitude for use in rotatem
mu = convertlat([a e], phiNew, 'geodetic', 'rectifying', 'nocheck');

%--------------------------------------------------------------------------

function [a, e, M0, tolerance, radius] = deriveParameters(mstruct)

[a, e, radius] = ellipsoidpropsRectifying(mstruct);
tolerance = 5*toleranceInRadians;

phi0 = toRadians(mstruct.angleunits, mstruct.origin(1));
M0 = meridianarc(0, phi0, [a e]);
