proj.c

/*
*				proj.c
*
* Spherical map projections.
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
*	This file part of:	AstrOmatic WCS library
*
*	Copyright:		(C) 2000-2010 IAP/CNRS/UPMC
*				(C) 1995-1999 Mark Calabretta
*
*	Authors:		Emmanuel Bertin (this version)
*				Mark Calabretta (original version)
*
*	Licenses:		GNU General Public License
*
*	AstrOmatic software is free software: you can redistribute it and/or
*	modify it under the terms of the GNU General Public License as
*	published by the Free Software Foundation, either version 3 of the
*	License, or (at your option) any later version.
*	AstrOmatic software is distributed in the hope that it will be useful,
*	but WITHOUT ANY WARRANTY; without even the implied warranty of
*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*	GNU General Public License for more details.
*	You should have received a copy of the GNU General Public License
*	along with AstrOmatic software.
*	If not, see <http://www.gnu.org/licenses/>.
*
*	Last modified:		10/10/2010
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/*============================================================================
*
*   WCSLIB - an implementation of the FITS WCS proposal.
*   Copyright (C) 1995-1999, Mark Calabretta
*
*   This library is free software; you can redistribute it and/or modify it
*   under the terms of the GNU Library General Public License as published
*   by the Free Software Foundation; either version 2 of the License, or (at
*   your option) any later version.
*
*   This library is distributed in the hope that it will be useful, but
*   WITHOUT ANY WARRANTY; without even the implied warranty of
*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library
*   General Public License for more details.
*
*   You should have received a copy of the GNU Library General Public License
*   along with this library; if not, write to the Free Software Foundation,
*   Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
*   Correspondence concerning WCSLIB may be directed to:
*      Internet email: mcalabre@atnf.csiro.au
*      Postal address: Dr. Mark Calabretta,
*                      Australia Telescope National Facility,
*                      P.O. Box 76,
*                      Epping, NSW, 2121,
*                      AUSTRALIA
*
*=============================================================================
*
*   C implementation of the spherical map projections recognized by the FITS
*   "World Coordinate System" (WCS) convention.
*
*   Summary of routines
*   -------------------
*   Each projection is implemented via separate functions for the forward,
*   *fwd(), and reverse, *rev(), transformation.
*
*   Initialization routines, *set(), compute intermediate values from the
*   projection parameters but need not be called explicitly - see the
*   explanation of prj.flag below.
*
*      azpset azpfwd azprev   AZP: zenithal/azimuthal perspective
*      tanset tanfwd tanrev   TAN: gnomonic
*      sinset sinfwd sinrev   SIN: orthographic/synthesis
*      stgset stgfwd stgrev   STG: stereographic
*      arcset arcfwd arcrev   ARC: zenithal/azimuthal equidistant
*      zpnset zpnfwd zpnrev   ZPN: zenithal/azimuthal polynomial
*      zeaset zeafwd zearev   ZEA: zenithal/azimuthal equal area
*      airset airfwd airrev   AIR: Airy
*      cypset cypfwd cyprev   CYP: cylindrical perspective
*      carset carfwd carrev   CAR: Cartesian
*      merset merfwd merrev   MER: Mercator
*      ceaset ceafwd cearev   CEA: cylindrical equal area
*      copset copfwd coprev   COP: conic perspective
*      codset codfwd codrev   COD: conic equidistant
*      coeset coefwd coerev   COE: conic equal area
*      cooset coofwd coorev   COO: conic orthomorphic
*      bonset bonfwd bonrev   BON: Bonne
*      pcoset pcofwd pcorev   PCO: polyconic
*      glsset glsfwd glsrev   GLS: Sanson-Flamsteed (global sinusoidal)
*      parset parfwd parrev   PAR: parabolic
*      aitset aitfwd aitrev   AIT: Hammer-Aitoff
*      molset molfwd molrev   MOL: Mollweide
*      cscset cscfwd cscrev   CSC: COBE quadrilateralized spherical cube
*      qscset qscfwd qscrev   QSC: quadrilateralized spherical cube
*      tscset tscfwd tscrev   TSC: tangential spherical cube
*      tnxset tnxfwd tnxrev   TNX: IRAF's gnomonic
*
*
*   Initialization routine; *set()
*   ------------------------------
*   Initializes members of a prjprm data structure which hold intermediate
*   values.  Note that this routine need not be called directly; it will be
*   invoked by prjfwd() and prjrev() if the "flag" structure member is
*   anything other than a predefined magic value.
*
*   Given and/or returned:
*      prj      prjprm*  Projection parameters (see below).
*
*   Function return value:
*               int      Error status
*                           0: Success.
*                           1: Invalid projection parameters.
*
*   Forward transformation; *fwd()
*   -----------------------------
*   Compute (x,y) coordinates in the plane of projection from native spherical
*   coordinates (phi,theta).
*
*   Given:
*      phi,     const double
*      theta             Longitude and latitude of the projected point in
*                        native spherical coordinates, in degrees.
*
*   Given and returned:
*      prj      prjprm*  Projection parameters (see below).
*
*   Returned:
*      x,y      double*  Projected coordinates.
*
*   Function return value:
*               int      Error status
*                           0: Success.
*                           1: Invalid projection parameters.
*                           2: Invalid value of (phi,theta).
*
*   Reverse transformation; *rev()
*   -----------------------------
*   Compute native spherical coordinates (phi,theta) from (x,y) coordinates in
*   the plane of projection.
*
*   Given:
*      x,y      const double
*                        Projected coordinates.
*
*   Given and returned:
*      prj      prjprm*  Projection parameters (see below).
*
*   Returned:
*      phi,     double*  Longitude and latitude of the projected point in
*      theta             native spherical coordinates, in degrees.
*
*   Function return value:
*               int      Error status
*                           0: Success.
*                           1: Invalid projection parameters.
*                           2: Invalid value of (x,y).
*
*   Projection parameters
*   ---------------------
*   The prjprm struct consists of the following:
*
*      int flag
*         This flag must be set to zero whenever any of p[10] or r0 are set
*         or changed.  This signals the initialization routine to recompute
*         intermediaries.  flag may also be set to -1 to disable strict bounds
*         checking for the AZP, TAN, SIN, ZPN, and COP projections.
*      double r0
*         r0; The radius of the generating sphere for the projection, a linear
*         scaling parameter.  If this is zero, it will be reset to the default
*         value of 180/pi (the value for FITS WCS).
*      double p[10]
*         The first 10 elements contain projection parameters which correspond
*         to the PROJPn keywords in FITS, so p[0] is PROJP0, and p[9] is
*         PROJP9.  Many projections use p[1] (PROJP1) and some also use p[2]
*         (PROJP2).  ZPN is the only projection which uses any of the others.
*
*   The remaining members of the prjprm struct are maintained by the
*   initialization routines and should not be modified.  This is done for the
*   sake of efficiency and to allow an arbitrary number of contexts to be
*   maintained simultaneously.
*
*      int n
*      double w[10]
*         Intermediate values derived from the projection parameters.
*
*   Usage of the p[] array as it applies to each projection is described in
*   the prologue to each trio of projection routines.
*
*   Argument checking
*   -----------------
*   Forward routines:
*
*      The values of phi and theta (the native longitude and latitude)
*      normally lie in the range [-180,180] for phi, and [-90,90] for theta.
*      However, all forward projections will accept any value of phi and will
*      not normalize it.
*
*      The forward projection routines do not explicitly check that theta lies
*      within the range [-90,90].  They do check for any value of theta which
*      produces an invalid argument to the projection equations (e.g. leading
*      to division by zero).  The forward routines for AZP, TAN, SIN, ZPN, and
*      COP also return error 2 if (phi,theta) corresponds to the overlapped
*      (far) side of the projection but also return the corresponding value of
*      (x,y).  This strict bounds checking may be relaxed by setting prj->flag
*      to -1 (rather than 0) when these projections are initialized.
*
*   Reverse routines:
*
*      Error checking on the projected coordinates (x,y) is limited to that
*      required to ascertain whether a solution exists.  Where a solution does
*      exist no check is made that the value of phi and theta obtained lie
*      within the ranges [-180,180] for phi, and [-90,90] for theta.
*
*   Accuracy
*   --------
*   Closure to a precision of at least 1E-10 degree of longitude and latitude
*   has been verified for typical projection parameters on the 1 degree grid
*   of native longitude and latitude (to within 5 degrees of any latitude
*   where the projection may diverge).
*
*   Author: Mark Calabretta, Australia Telescope National Facility
*   IRAF's TNX added by E.Bertin 2000/08/23
*   Behaviour of Cartesian-like projections modified by E.Bertin 2005/08/29
*   $Id: proj.c,v 1.1.1.1 2008/01/10 18:00:00 bertin Exp $
*===========================================================================*/

#ifdef HAVE_CONFIG_H
#include	"config.h"
#endif

#ifdef HAVE_MATHIMF_H
#include <mathimf.h>
#else
#include <math.h>
#endif
#include <stdio.h>
#include <stdlib.h>
#include "poly.h"
#include "proj.h"
#include "tnx.h"
#include "wcsmath.h"
#include "wcstrig.h"

/* Map error number to error message for each function. */
const char *prjset_errmsg[] = {
   0,
   "Invalid projection parameters"};

const char *prjfwd_errmsg[] = {
   0,
   "Invalid projection parameters",
   "Invalid value of (phi,theta)"};

const char *prjrev_errmsg[] = {
   0,
   "Invalid projection parameters",
   "Invalid value of (x,y)"};

#define wcs_copysign(X, Y) ((Y) < 0.0 ? -fabs(X) : fabs(X))

/*============================================================================
*   AZP: zenithal/azimuthal perspective projection.
*
*   Given:
*      prj->p[1]   AZP distance parameter, mu in units of r0.
*
*   Given and/or returned:
*      prj->r0     r0; reset to 180/pi if 0.
*      prj->w[0]   r0*(mu+1)
*      prj->w[1]   1/prj->w[0]
*      prj->w[2]   Boundary parameter, -mu    for |mu| <= 1,
*                                      -1/mu  for |mu| >= 1.
*===========================================================================*/

int azpset(prj)

struct prjprm *prj;

{
   if (prj->r0 == 0.0) prj->r0 = R2D;

   prj->w[0] = prj->r0*(prj->p[1] + 1.0);
   if (prj->w[0] == 0.0) {
      return 1;
   }

   prj->w[1] = 1.0/prj->w[0];
   if (fabs(prj->p[1]) <= 1.0) {
      prj->w[2] = -prj->p[1];
   } else {
      prj->w[2] = -1.0/prj->p[1];
   }

   if (prj->flag == -1) {
      prj->flag = -PRJSET;
   } else {
      prj->flag = PRJSET;
   }

   return 0;
}

/*--------------------------------------------------------------------------*/

int azpfwd(phi, theta, prj, x, y)

const double phi, theta;
struct prjprm *prj;
double *x, *y;

{
   double r, s, sthe;

   if (abs(prj->flag) != PRJSET) {
      if (azpset(prj)) return 1;
   }

   sthe = wcs_sind(theta);

   s = prj->p[1] + sthe;
   if (s == 0.0) {
      return 2;
   }

   r =  prj->w[0]*wcs_cosd(theta)/s;
   *x =  r*wcs_sind(phi);
   *y = -r*wcs_cosd(phi);

   if (prj->flag == PRJSET && sthe < prj->w[2]) {
      return 2;
   }

   return 0;
}

/*--------------------------------------------------------------------------*/

int azprev(x, y, prj, phi, theta)

const double x, y;
struct prjprm *prj;
double *phi, *theta;

{
   double r, rho, s;
   const double tol = 1.0e-13;

   if (abs(prj->flag) != PRJSET) {
      if (azpset(prj)) return 1;
   }

   r = sqrt(x*x + y*y);
   if (r == 0.0) {
      *phi = 0.0;
   } else {
      *phi = wcs_atan2d(x, -y);
   }

   rho = r*prj->w[1];
   s = rho*prj->p[1]/sqrt(rho*rho+1.0);
   if (fabs(s) > 1.0) {
      if (fabs(s) > 1.0+tol) {
         return 2;
      }
      *theta = wcs_atan2d(1.0,rho) - wcs_copysign(90.0,s);
   } else {
      *theta = wcs_atan2d(1.0,rho) - wcs_asind(s);
   }

   return 0;
}

/*============================================================================
*   TAN: gnomonic projection.
*
*   Given and/or returned:
*      prj->r0     r0; reset to 180/pi if 0.
*===========================================================================*/

int tanset(prj)

struct prjprm *prj;

{
   int	k;

   if (prj->r0 == 0.0) prj->r0 = R2D;

   if (prj->flag == -1) {
      prj->flag = -PRJSET;
   } else {
      prj->flag = PRJSET;
   } 

   for (k = 99; k >= 0 && prj->p[k] == 0.0 && prj->p[k+100] == 0.0; k--);
   if (k < 0)
     k = 0;

   prj->n = k;

   return 0;
}

/*--------------------------------------------------------------------------*/

int tanfwd(phi, theta, prj, x, y)

const double phi, theta;
struct prjprm *prj;
double *x, *y;

{
   double r, s, xp[2];

   if (abs(prj->flag) != PRJSET) {
      if(tanset(prj)) return 1;
   }

   s = wcs_sind(theta);
   if (s == 0.0) return 2;

   r =  prj->r0*wcs_cosd(theta)/s;
   xp[0] =  r*wcs_sind(phi);
   xp[1] = -r*wcs_cosd(phi);
   *x = prj->inv_x? poly_func(prj->inv_x, xp) : xp[0];
   *y = prj->inv_y? poly_func(prj->inv_y, xp) : xp[1];

   if (prj->flag == PRJSET && s < 0.0) {
      return 2;
   }

   return 0;
}

/*--------------------------------------------------------------------------*/

int tanrev(x, y, prj, phi, theta)

const double x, y;
struct prjprm *prj;
double *phi, *theta;

{
   double	xp,yp, rp;

   if (abs(prj->flag) != PRJSET) {
      if (tanset(prj)) return 1;
   }

   if (prj->n)
     raw_to_pv(prj, x,y, &xp, &yp);
   else
     {
     xp = x;
     yp = y;
     }
   rp = sqrt(xp*xp+yp*yp);
   if (rp == 0.0) {
      *phi = 0.0;
   } else {
      *phi = wcs_atan2d(xp, -yp);
   }
   *theta = wcs_atan2d(prj->r0, rp);

   return 0;
}

/*============================================================================
*   SIN: orthographic/synthesis projection.
*
*   Given:
*      prj->p[1:2] SIN obliqueness parameters, alpha and beta.
*
*   Given and/or returned:
*      prj->r0     r0; reset to 180/pi if 0.
*      prj->w[0]   1/r0
*      prj->w[1]   alpha**2 + beta**2
*      prj->w[2]   2*(alpha**2 + beta**2)
*      prj->w[3]   2*(alpha**2 + beta**2 + 1)
*      prj->w[4]   alpha**2 + beta**2 - 1
*===========================================================================*/

int sinset(prj)

struct prjprm *prj;

{
   if (prj->r0 == 0.0) prj->r0 = R2D;

   prj->w[0] = 1.0/prj->r0;
   prj->w[1] = prj->p[1]*prj->p[1] + prj->p[2]*prj->p[2];
   prj->w[2] = 2.0*prj->w[1];
   prj->w[3] = prj->w[2] + 2.0;
   prj->w[4] = prj->w[1] - 1.0;

   if (prj->flag == -1) {
      prj->flag = -PRJSET;
   } else {
      prj->flag = PRJSET;
   }

   return 0;
}

/*--------------------------------------------------------------------------*/

int sinfwd(phi, theta, prj, x, y)

const double phi, theta;
struct prjprm *prj;
double *x, *y;

{
   double cphi, cthe, sphi, t, z;

   if (abs(prj->flag) != PRJSET) {
      if (sinset(prj)) return 1;
   }

   t = (90.0 - fabs(theta))*D2R;
   if (t < 1.0e-5) {
      if (theta > 0.0) {
         z = -t*t/2.0;
      } else {
         z = -2.0 + t*t/2.0;
      }
      cthe = t;
   } else {
      z =  wcs_sind(theta) - 1.0;
      cthe = wcs_cosd(theta);
   }

   cphi = wcs_cosd(phi);
   sphi = wcs_sind(phi);
   *x =  prj->r0*(cthe*sphi + prj->p[1]*z);
   *y = -prj->r0*(cthe*cphi + prj->p[2]*z);
   /* Validate this solution. */
   if (prj->flag == PRJSET) {
      if (prj->w[1] == 0.0) {
         /* Orthographic projection. */
         if (theta < 0.0) {
            return 2;
         }
      } else {
         /* "Synthesis" projection. */
         t = wcs_atand(prj->p[1]*sphi + prj->p[2]*cphi);
         if (theta < t) {
            return 2;
         }
      }
   }

   return 0;
}

/*--------------------------------------------------------------------------*/

int sinrev (x, y, prj, phi, theta)

const double x, y;
struct prjprm *prj;
double *phi, *theta;

{
   const double tol = 1.0e-13;
   double a, b, c, d, r2, sth, sth1, sth2, sxy, x0, xp, y0, yp, z;

   if (abs(prj->flag) != PRJSET) {
      if (sinset(prj)) return 1;
   }

   /* Compute intermediaries. */
   x0 = x*prj->w[0];
   y0 = y*prj->w[0];
   r2 = x0*x0 + y0*y0;

   if (prj->w[1] == 0.0) {
      /* Orthographic projection. */
      if (r2 != 0.0) {
         *phi = wcs_atan2d(x0, -y0);
      } else {
         *phi = 0.0;
      }

      if (r2 < 0.5) {
         *theta = wcs_acosd(sqrt(r2));
      } else if (r2 <= 1.0) {
         *theta = wcs_asind(sqrt(1.0 - r2));
      } else {
         return 2;
      }

   } else {
      /* "Synthesis" projection. */
      if (r2 < 1.0e-10) {
         /* Use small angle formula. */
         z = -r2/2.0;
         *theta = 90.0 - R2D*sqrt(r2/(1.0 - x0*prj->p[1] + y0*prj->p[2]));

      } else {
         sxy = 2.0*(prj->p[1]*x0 - prj->p[2]*y0);

         a = prj->w[3];
         b = -(sxy + prj->w[2]);
         c = r2 + sxy + prj->w[4];
         d = b*b - 2.0*a*c;

         /* Check for a solution. */
         if (d < 0.0) {
            return 2;
         }
         d = sqrt(d);

         /* Choose solution closest to pole. */
         sth1 = (-b + d)/a;
         sth2 = (-b - d)/a;
         sth = (sth1>sth2) ? sth1 : sth2;
         if (sth > 1.0) {
            if (sth-1.0 < tol) {
               sth = 1.0;
            } else {
               sth = (sth1<sth2) ? sth1 : sth2;
            }
         }
         if (sth > 1.0 || sth < -1.0) {
            return 2;
         }

         *theta = wcs_asind(sth);
         z = sth - 1.0;
      }

      xp = -y0 - prj->p[2]*z;
      yp =  x0 - prj->p[1]*z;
      if (xp == 0.0 && yp == 0.0) {
         *phi = 0.0;
      } else {
         *phi = wcs_atan2d(yp,xp);
      }
   }

   return 0;
}

/*============================================================================
*   STG: stereographic projection.
*
*   Given and/or returned:
*      prj->r0     r0; reset to 180/pi if 0.
*      prj->w[0]   2*r0
*      prj->w[1]   1/(2*r0)
*===========================================================================*/

int stgset(prj)

struct prjprm *prj;

{
   if (prj->r0 == 0.0) {
      prj->r0 = R2D;
      prj->w[0] = 360.0/PI;
      prj->w[1] = PI/360.0;
   } else {
      prj->w[0] = 2.0*prj->r0;
      prj->w[1] = 1.0/prj->w[0];
   }

   prj->flag = PRJSET;
   return 0;
}

/*--------------------------------------------------------------------------*/

int stgfwd(phi, theta, prj, x, y)

const double phi, theta;
struct prjprm *prj;
double *x, *y;

{
   double r, s;

   if (prj->flag != PRJSET) {
      if (stgset(prj)) return 1;
   }

   s = 1.0 + wcs_sind(theta);
   if (s == 0.0) {
      return 2;
   }

   r =  prj->w[0]*wcs_cosd(theta)/s;
   *x =  r*wcs_sind(phi);
   *y = -r*wcs_cosd(phi);

   return 0;
}

/*--------------------------------------------------------------------------*/

int stgrev(x, y, prj, phi, theta)

const double x, y;
struct prjprm *prj;
double *phi, *theta;

{
   double r;

   if (prj->flag != PRJSET) {
      if (stgset(prj)) return 1;
   }

   r = sqrt(x*x + y*y);
   if (r == 0.0) {
      *phi = 0.0;
   } else {
      *phi = wcs_atan2d(x, -y);
   }
   *theta = 90.0 - 2.0*wcs_atand(r*prj->w[1]);

   return 0;
}

/*============================================================================
*   ARC: zenithal/azimuthal equidistant projection.
*
*   Given and/or returned:
*      prj->r0     r0; reset to 180/pi if 0.
*      prj->w[0]   r0*(pi/180)
*      prj->w[1]   (180/pi)/r0
*===========================================================================*/

int arcset(prj)

struct prjprm *prj;

{
   if (prj->r0 == 0.0) {
      prj->r0 = R2D;
      prj->w[0] = 1.0;
      prj->w[1] = 1.0;
   } else {
      prj->w[0] = prj->r0*D2R;
      prj->w[1] = 1.0/prj->w[0];
   }

   prj->flag = PRJSET;

   return 0;
}

/*--------------------------------------------------------------------------*/

int arcfwd(phi, theta, prj, x, y)

const double phi, theta;
struct prjprm *prj;
double *x, *y;

{
   double r;

   if (prj->flag != PRJSET) {
      if (arcset(prj)) return 1;
   }

   r =  prj->w[0]*(90.0 - theta);
   *x =  r*wcs_sind(phi);
   *y = -r*wcs_cosd(phi);

   return 0;
}

/*--------------------------------------------------------------------------*/

int arcrev(x, y, prj, phi, theta)

const double x, y;
struct prjprm *prj;
double *phi, *theta;

{
   double r;

   if (prj->flag != PRJSET) {
      if (arcset(prj)) return 1;
   }

   r = sqrt(x*x + y*y);
   if (r == 0.0) {
      *phi = 0.0;
   } else {
      *phi = wcs_atan2d(x, -y);
   }
   *theta = 90.0 - r*prj->w[1];

   return 0;
}

/*============================================================================
*   ZPN: zenithal/azimuthal polynomial projection.
*
*   Given:
*      prj->p[0:99] Polynomial coefficients.
*
*   Given and/or returned:
*      prj->r0     r0; reset to 180/pi if 0.
*      prj->n      Degree of the polynomial, N.
*      prj->w[0]   Co-latitude of the first point of inflection (N > 2).
*      prj->w[1]   Radius of the first point of inflection (N > 2).
*===========================================================================*/

int zpnset(prj)

struct prjprm *prj;

{
   int   i, j, k;
   double d, d1, d2, r, zd, zd1, zd2;
   const double tol = 1.0e-13;

   if (prj->r0 == 0.0) prj->r0 = R2D;

   /* Find the highest non-zero coefficient. */
   for (k = 99; k >= 0 && prj->p[k] == 0.0; k--);
   if (k < 0) return 1;

   prj->n = k;

   if (k >= 3) {
      /* Find the point of inflection closest to the pole. */
      zd1 = 0.0;
      d1  = prj->p[1];
      if (d1 <= 0.0) {
         return 1;
      }

      /* Find the point where the derivative first goes negative. */
      for (i = 0; i < 180; i++) {
         zd2 = i*D2R;
         d2  = 0.0;
         for (j = k; j > 0; j--) {
            d2 = d2*zd2 + j*prj->p[j];
         }

         if (d2 <= 0.0) break;
         zd1 = zd2;
         d1  = d2;
      }

      if (i == 180) {
         /* No negative derivative -> no point of inflection. */
         zd = PI;
      } else {
         /* Find where the derivative is zero. */
         for (i = 1; i <= 10; i++) {
            zd = zd1 - d1*(zd2-zd1)/(d2-d1);

            d = 0.0;
            for (j = k; j > 0; j--) {
               d = d*zd + j*prj->p[j];
            }

            if (fabs(d) < tol) break;

            if (d < 0.0) {
               zd2 = zd;
               d2  = d;
            } else {
               zd1 = zd;
               d1  = d;
            }
         }
      }

      r = 0.0;
      for (j = k; j >= 0; j--) {
         r = r*zd + prj->p[j];
      }
      prj->w[0] = zd;
      prj->w[1] = r;
   }

   if (prj->flag == -1) {
      prj->flag = -PRJSET;
   } else {
      prj->flag = PRJSET;
   }

   return 0;
}

/*--------------------------------------------------------------------------*/

int zpnfwd(phi, theta, prj, x, y)

const double phi, theta;
struct prjprm *prj;
double *x, *y;

{
   int   j;
   double r, s;

   if (abs(prj->flag) != PRJSET) {
      if (zpnset(prj)) return 1;
   }

   s = (90.0 - theta)*D2R;
   r = 0.0;
   for (j = prj->n; j >= 0; j--) {
      r = r*s + prj->p[j];
   }
   r = prj->r0*r;

   *x =  r*wcs_sind(phi);
   *y = -r*wcs_cosd(phi);

   if (prj->flag == PRJSET && s > prj->w[0]) {
      return 2;
   }

   return 0;
}

/*--------------------------------------------------------------------------*/

int zpnrev(x, y, prj, phi, theta)

const double x, y;
struct prjprm *prj;
double *phi, *theta;

{
   int   i, j, k;
   double a, b, c, d, lambda, r, r1, r2, rt, zd, zd1, zd2;
   const double tol = 1.0e-13;

   if (abs(prj->flag) != PRJSET) {
      if (zpnset(prj)) return 1;
   }

   k = prj->n;

   r = sqrt(x*x + y*y)/prj->r0;

   if (k < 1) {
      /* Constant - no solution. */
      return 1;
   } else if (k == 1) {
      /* Linear. */
      zd = (r - prj->p[0])/prj->p[1];
   } else if (k == 2) {
      /* Quadratic. */
      a = prj->p[2];
      b = prj->p[1];
      c = prj->p[0] - r;

      d = b*b - 4.0*a*c;
      if (d < 0.0) {
         return 2;
      }
      d = sqrt(d);

      /* Choose solution closest to pole. */
      zd1 = (-b + d)/(2.0*a);
      zd2 = (-b - d)/(2.0*a);
      zd  = (zd1<zd2) ? zd1 : zd2;
      if (zd < -tol) zd = (zd1>zd2) ? zd1 : zd2;
      if (zd < 0.0) {
         if (zd < -tol) {
            return 2;
         }
         zd = 0.0;
      } else if (zd > PI) {
         if (zd > PI+tol) {
            return 2;
         }
         zd = PI;
      }
   } else {
      /* Higher order - solve iteratively. */
      zd1 = 0.0;
      r1  = prj->p[0];
      zd2 = prj->w[0];
      r2  = prj->w[1];

      if (r < r1) {
         if (r < r1-tol) {
            return 2;
         }
         zd = zd1;
      } else if (r > r2) {
         if (r > r2+tol) {
            return 2;
         }