head	2.17;
access;
symbols;
locks; strict;
comment	@c @;


2.17
date	2002.04.03.01.24.44;	author mcalabre;	state Exp;
branches;
next	2.16;

2.16
date	2001.11.16.03.52.15;	author mcalabre;	state Exp;
branches;
next	2.15;

2.15
date	2000.12.06.23.53.57;	author mcalabre;	state Exp;
branches;
next	2.14;

2.14
date	2000.12.06.00.50.47;	author mcalabre;	state Exp;
branches;
next	2.13;

2.13
date	2000.05.10.04.51.30;	author mcalabre;	state Exp;
branches;
next	2.12;

2.12
date	2000.02.01.04.17.25;	author mcalabre;	state Exp;
branches;
next	2.11;

2.11
date	99.12.13.03.05.34;	author mcalabre;	state Exp;
branches;
next	2.10;

2.10
date	99.11.19.05.54.32;	author mcalabre;	state Exp;
branches;
next	2.9;

2.9
date	99.11.19.01.17.18;	author mcalabre;	state Exp;
branches;
next	2.8;

2.8
date	99.04.09.05.09.27;	author mcalabre;	state Exp;
branches;
next	2.7;

2.7
date	99.04.06.07.58.16;	author mcalabre;	state Exp;
branches;
next	2.6;

2.6
date	98.02.25.07.02.23;	author mcalabre;	state Exp;
branches;
next	2.5;

2.5
date	97.03.19.05.34.13;	author mcalabre;	state Exp;
branches;
next	2.3;

2.3
date	96.05.07.21.06.38;	author mcalabre;	state Exp;
branches;
next	2.2;

2.2
date	95.12.22.02.40.14;	author mcalabre;	state Exp;
branches;
next	2.1;

2.1
date	95.11.09.03.14.18;	author mcalabre;	state Exp;
branches;
next	2.0;

2.0
date	95.09.11.08.35.22;	author mcalabre;	state Exp;
branches;
next	1.2;

1.2
date	95.09.11.08.34.50;	author mcalabre;	state Exp;
branches;
next	1.1;

1.1
date	95.01.31.03.20.39;	author mcalabre;	state Exp;
branches;
next	;


desc
@Implementation of the WCS spherical projection equations.
@


2.17
log
@Updated date in copyright notice.
@
text
@*=======================================================================
*
*   WCSLIB - an implementation of the FITS WCS proposal.
*   Copyright (C) 1995-2002, 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
*
*=======================================================================
*
*   FORTRAN routines for the spherical map projections recognized by the
*   FITS "World Coordinate System" (WCS) convention.
*
*   Summary of routines
*   -------------------
*   Each projection is implemented via separate subroutines 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(11) below.
*
*      PRJSET PRJFWD PRJREV   Driver routines (see below).
*
*      AZPSET AZPFWD AZPREV   AZP: zenithal/azimuthal perspective
*      SZPSET SZPFWD SZPREV   SZP: slant zenithal perspective
*      TANSET TANFWD TANREV   TAN: gnomonic
*      STGSET STGFWD STGREV   STG: stereographic
*      SINSET SINFWD SINREV   SIN: orthographic/synthesis
*      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
*      CEASET CEAFWD CEAREV   CEA: cylindrical equal area
*      CARSET CARFWD CARREV   CAR: Cartesian
*      MERSET MERFWD MERREV   MER: Mercator
*      SFLSET SFLFWD SFLREV   SFL: Sanson-Flamsteed
*      PARSET PARFWD PARREV   PAR: parabolic
*      MOLSET MOLFWD MOLREV   MOL: Molweide
*      AITSET AITFWD AITREV   AIT: Hammer-Aitoff
*      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
*      TSCSET TSCFWD TSCREV   TSC: tangential spherical cube
*      CSCSET CSCFWD CSCREV   CSC: COBE quadrilateralized spherical cube
*      QSCSET QSCFWD QSCREV   QSC: quadrilateralized spherical cube
*
*
*   Driver routines; PRJSET, PRJFWD & PRJREV
*   ----------------------------------------
*   A set of driver routines are available for use as a generic
*   interface to the specific projection routines.  The interfaces to
*   PRJFWD and PRJREV are the same as those of the forward and reverse
*   transformation routines for the specific projections.  They use the
*   value of PRJ(11) to identify the projection.
*
*   The interface to PRJSET differs from that of the initialization
*   routines for the specific projections and unlike them it must be
*   invoked explicitly to use PRJFWD and PRJREV.
*
*   Given and/or returned:
*      PCODE    C*3      WCS projection code.
*
*   Given and/or returned:
*      PRJ      D(0:20)  Projection parameters (see below).
*
*   Returned:
*      PHI0     D        Native longitude and latitude of the reference
*      THETA0   D        point, in degrees.
*
*      IERR     I        Error status
*                           0: Success.
*                           1: Invalid projection parameters.
*
*
*   Initialization routine; *SET
*   ----------------------------
*   Initializes elements of a PRJ array which hold intermediate values.
*   Note that this routine need not be called directly; it will be
*   invoked by the forward or reverse projection routine if PRJ(11) is
*   zero.
*
*   Given and/or returned:
*      PRJ      D(0:20)  Projection parameters (see below).
*
*   Returned:
*      IERR     I        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,     D        Longitude and latitude of the projected point
*      THETA             in native spherical coordinates, in degrees.
*
*   Given and returned:
*      PRJ      D(0:20)  Projection parameters (see below).
*
*   Returned:
*      X,Y      D        Projected coordinates.
*      IERR     I        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      D        Projected coordinates.
*
*   Given and returned:
*      PRJ      D(0:20)  Projection parameters (see below).
*
*   Returned:
*      PHI,     D        Longitude and latitude of the projected point
*      THETA             in native spherical coordinates, in degrees.
*      IERR     I        Error status
*                           0: Success.
*                           1: Invalid projection parameters.
*                           2: Invalid value of (x,y).
*
*   Projection parameters
*   ---------------------
*   The PRJ array consists of the following:
*
*      PRJ(0:9)
*         The first 10 elements contain projection parameters which
*         correspond to the PROJPn keywords in FITS, so PRJ(0) is
*         PROJP0, and PRJ(9) is PROJP9.  Many projections use PRJ(1)
*         PROJP1 and some also use PRJ(2) (PROJP2 ).  ZPN is the only
*         projection which uses any of the others.
*      PRJ(10)
*         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).
*      PRJ(11)
*         This must be set to zero whenever any of PRJ(0:10) are set or
*         changed.  This signals the initialization routines to
*         recompute intermediaries.  PRJ(11) may also be set to -1 to
*         disable strict bounds checking for the AZP, SZP, TAN, SIN,
*         ZPN, and COP projections.  The initialization routine for each
*         projection resets PRJ(11) to a code unique to the projection.
*
*   The remaining elements of the PRJ array 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.
*
*      PRJ(12:20)
*         Intermediate values derived from the projection parameters.
*
*   Usage of the PRJ 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, SZP, 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(11) 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 1D-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
*   $Id: proj.f,v 2.16 2001/11/16 03:52:15 mcalabre Exp mcalabre $
*=======================================================================
      SUBROUTINE PRJSET (PCODE, PRJ, PHI0, THETA0, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR, IREF(26), NDX
      DOUBLE PRECISION PHI0, PRJ(0:20), THETA0
      CHARACTER PCODES*103, PCODE*3

      DATA PCODES(01:40)  /'AZP SZP TAN STG SIN ARC ZPN ZEA AIR CYP '/
      DATA PCODES(41:80)  /'CEA CAR MER SFL PAR MOL AIT COP COE COD '/
      DATA PCODES(81:103) /'COO BON PCO TSC CSC QSC'/
      DATA IREF /2,2,2,2,2,2,2,2,2,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0/
*-----------------------------------------------------------------------
*     Determine the projection type.
      NDX = (INDEX(PCODES,PCODE) - 1)/4 + 1
      IF (NDX.EQ.0) THEN
*        Unrecognized projection code.
         IERR = 1
         RETURN
      END IF

*     Native longitude and latitude of the reference point.
      PHI0 = 0D0
      IF (IREF(NDX).EQ.0) THEN
         THETA0 = 0D0
      ELSE IF (IREF(NDX).EQ.1) THEN
         THETA0 = PRJ(1)
      ELSE IF (IREF(NDX).EQ.2) THEN
         THETA0 = 90D0
      ELSE
         IERR = 1
         RETURN
      END IF

*     The computed GOTO is used instead of an IF-block for efficiency.
      GO TO (101, 102, 103, 104, 105, 106, 107, 108, 109,
     *       201, 202, 203, 204,
     *       301, 302, 303,
     *       401,
     *       501, 502, 503, 504,
     *       601, 602,
     *       701, 702, 703) NDX

*     AZP: zenithal/azimuthal perspective.
 101  CALL AZPSET (PRJ, IERR)
      RETURN

*     SZP: slant zenithal perspective.
 102  CALL SZPSET (PRJ, IERR)
      RETURN

*     TAN: gnomonic (tan).
 103  CALL TANSET (PRJ, IERR)
      RETURN

*     SIN: orthographic (sine).
 104  CALL STGSET (PRJ, IERR)
      RETURN

*     STG: stereographic.
 105  CALL SINSET (PRJ, IERR)
      RETURN

*     ARC: zenithal/azimuthal equidistant.
 106  CALL ARCSET (PRJ, IERR)
      RETURN

*     ZPN: zenithal/azimuthal polynomial.
 107  CALL ZPNSET (PRJ, IERR)
      RETURN

*     ZEA: zenithal/azimuthal equal area.
 108  CALL ZEASET (PRJ, IERR)
      RETURN

*     AIR: Airy's zenithal projection.
 109  CALL AIRSET (PRJ, IERR)
      RETURN

*     CYP: cylindrical perspective.
 201  CALL CYPSET (PRJ, IERR)
      RETURN

*     CAR: Cartesian.
 202  CALL CEASET (PRJ, IERR)
      RETURN

*     MER: Mercator's.
 203  CALL CARSET (PRJ, IERR)
      RETURN

*     CEA: cylindrical equal area.
 204  CALL MERSET (PRJ, IERR)
      RETURN

*     SFL: Sanson-Flamsteed.
 301  CALL SFLSET (PRJ, IERR)
      RETURN

*     PAR: parabolic.
 302  CALL PARSET (PRJ, IERR)
      RETURN

*     MOL: Mollweide's projection.
 303  CALL MOLSET (PRJ, IERR)
      RETURN

*     AIT: Hammer-Aitoff.
 401  CALL AITSET (PRJ, IERR)
      RETURN

*     COP: conic perspective.
 501  CALL COPSET (PRJ, IERR)
      RETURN

*     COE: conic equal area.
 502  CALL COESET (PRJ, IERR)
      RETURN

*     COD: conic equidistant.
 503  CALL CODSET (PRJ, IERR)
      RETURN

*     COO: conic orthomorphic.
 504  CALL COOSET (PRJ, IERR)
      RETURN

*     BON: Bonne's projection.
 601  CALL BONSET (PRJ, IERR)
      RETURN

*     PCO: polyconic.
 602  CALL PCOSET (PRJ, IERR)
      RETURN

*     TSC: tangential spherical cube.
 701  CALL TSCSET (PRJ, IERR)
      RETURN

*     CSC: COBE quadrilateralized spherical cube.
 702  CALL CSCSET (PRJ, IERR)
      RETURN

*     QSC: quadrilateralized spherical cube.
 703  CALL QSCSET (PRJ, IERR)
      RETURN

      END

*-----------------------------------------------------------------------
      SUBROUTINE PRJFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   I, IERR, J, K, NDX(7,9)
      DOUBLE PRECISION PHI, PRJ(0:20), THETA, X, Y

      DATA (NDX(1,J),J=1,9) / 1,  2,  3,  4,  5,  6,  7,  8,  9/
      DATA (NDX(2,J),J=1,9) /10, 11, 12, 13,  0,  0,  0,  0,  0/
      DATA (NDX(3,J),J=1,9) /14, 15, 16,  0,  0,  0,  0,  0,  0/
      DATA (NDX(4,J),J=1,9) /17,  0,  0,  0,  0,  0,  0,  0,  0/
      DATA (NDX(5,J),J=1,9) /18, 19, 20, 21,  0,  0,  0,  0,  0/
      DATA (NDX(6,J),J=1,9) /22, 23,  0,  0,  0,  0,  0,  0,  0/
      DATA (NDX(7,J),J=1,9) /24, 25, 26,  0,  0,  0,  0,  0,  0/
*-----------------------------------------------------------------------
      K = NINT(ABS(PRJ(11)))
      I = K/100
      J = MOD(K, 100)

      IF (I.LT.1 .OR. I.GT.7 .OR. J.LT.1 .OR. J.GT.9) THEN
         IERR = 1
         RETURN
      END IF


*     The computed GOTO is used instead of an IF-block for efficiency.
      GO TO (101, 102, 103, 104, 105, 106, 107, 108, 109,
     *       201, 202, 203, 204,
     *       301, 302, 303,
     *       401,
     *       501, 502, 503, 504,
     *       601, 602,
     *       701, 702, 703) NDX(I,J)

      IERR = 1
      RETURN

*     AZP: zenithal/azimuthal perspective.
 101  CALL AZPFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     SZP: slant zenithal perspective.
 102  CALL SZPFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     TAN: gnomonic.
 103  CALL TANFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     STG: stereographic.
 104  CALL STGFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     SIN: orthographic.
 105  CALL SINFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     ARC: zenithal/azimuthal equidistant.
 106  CALL ARCFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     ZPN: zenithal/azimuthal polynomial.
 107  CALL ZPNFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     ZEA: zenithal/azimuthal equal area.
 108  CALL ZEAFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     AIR: Airy's zenithal projection.
 109  CALL AIRFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     CYP: cylindrical perspective.
 201  CALL CYPFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     CEA: cylindrical equal area.
 202  CALL CEAFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     CAR: Cartesian.
 203  CALL CARFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     MER: Mercator's.
 204  CALL MERFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     SFL: Sanson-Flamsteed.
 301  CALL SFLFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     PAR: parabolic.
 302  CALL PARFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     MOL: Mollweide's projection.
 303  CALL MOLFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     AIT: Hammer-Aitoff.
 401  CALL AITFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     COP: conic perspective.
 501  CALL COPFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     COE: conic equal area.
 502  CALL COEFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     COD: conic equidistant.
 503  CALL CODFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     COO: conic orthomorphic.
 504  CALL COOFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     BON: Bonne's projection.
 601  CALL BONFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     PCO: polyconic.
 602  CALL PCOFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     TSC: tangential spherical cube.
 701  CALL TSCFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     CSC: COBE quadrilateralized spherical cube.
 702  CALL CSCFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

*     QSC: quadrilateralized spherical cube.
 703  CALL QSCFWD (PHI, THETA, PRJ, X, Y, IERR)
      RETURN

      END

*-----------------------------------------------------------------------
      SUBROUTINE PRJREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   I, IERR, J, K, NDX(7,9)
      DOUBLE PRECISION PHI, PRJ(0:20), THETA, X, Y

      DATA (NDX(1,J),J=1,9) / 1,  2,  3,  4,  5,  6,  7,  8,  9/
      DATA (NDX(2,J),J=1,9) /10, 11, 12, 13,  0,  0,  0,  0,  0/
      DATA (NDX(3,J),J=1,9) /14, 15, 16,  0,  0,  0,  0,  0,  0/
      DATA (NDX(4,J),J=1,9) /17,  0,  0,  0,  0,  0,  0,  0,  0/
      DATA (NDX(5,J),J=1,9) /18, 19, 20, 21,  0,  0,  0,  0,  0/
      DATA (NDX(6,J),J=1,9) /22, 23,  0,  0,  0,  0,  0,  0,  0/
      DATA (NDX(7,J),J=1,9) /24, 25, 26,  0,  0,  0,  0,  0,  0/
*-----------------------------------------------------------------------
      K = NINT(ABS(PRJ(11)))
      I = K/100
      J = MOD(K, 100)

      IF (I.LT.1 .OR. I.GT.7 .OR. J.LT.1 .OR. J.GT.9) THEN
         IERR = 1
         RETURN
      END IF


*     The computed GOTO is used instead of an IF-block for efficiency.
      GO TO (101, 102, 103, 104, 105, 106, 107, 108, 109,
     *       201, 202, 203, 204,
     *       301, 302, 303,
     *       401,
     *       501, 502, 503, 504,
     *       601, 602,
     *       701, 702, 703) NDX(I,J)

      IERR = 1
      RETURN

*     AZP: zenithal/azimuthal perspective.
 101  CALL AZPREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     SZP: slant zenithal perspective.
 102  CALL SZPREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     TAN: gnomonic.
 103  CALL TANREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     STG: stereographic.
 104  CALL STGREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     SIN: orthographic.
 105  CALL SINREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     ARC: zenithal/azimuthal equidistant.
 106  CALL ARCREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     ZPN: zenithal/azimuthal polynomial.
 107  CALL ZPNREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     ZEA: zenithal/azimuthal equal area.
 108  CALL ZEAREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     AIR: Airy's zenithal projection.
 109  CALL AIRREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     CYP: cylindrical perspective.
 201  CALL CYPREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     CEA: cylindrical equal area.
 202  CALL CEAREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     CAR: Cartesian.
 203  CALL CARREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     MER: Mercator's.
 204  CALL MERREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     SFL: Sanson-Flamsteed.
 301  CALL SFLREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     PAR: parabolic.
 302  CALL PARREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     MOL: Mollweide's projection.
 303  CALL MOLREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     AIT: Hammer-Aitoff.
 401  CALL AITREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     COP: conic perspective.
 501  CALL COPREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     COE: conic equal area.
 502  CALL COEREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     COD: conic equidistant.
 503  CALL CODREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     COO: conic orthomorphic.
 504  CALL COOREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     BON: Bonne's projection.
 601  CALL BONREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     PCO: polyconic.
 602  CALL PCOREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     TSC: tangential spherical cube.
 701  CALL TSCREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     CSC: COBE quadrilateralized spherical cube.
 702  CALL CSCREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

*     QSC: quadrilateralized spherical cube.
 703  CALL QSCREV (X, Y, PRJ, PHI, THETA, IERR)
      RETURN

      END

*=======================================================================
*   AZP: zenithal/azimuthal perspective projection.
*
*   Given:
*      PRJ(1)    Distance parameter, mu in units of r0.
*      PRJ(2)    Tilt angle, gamma in degrees.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*      PRJ(11)   State flag, 101 or -101 if given < 0.
*
*   Returned:
*      PRJ(12)   r0*(mu+1)
*      PRJ(13)   tan(gamma)
*      PRJ(14)   sec(gamma)
*      PRJ(15)   cos(gamma)
*      PRJ(16)   sin(gamma)
*      PRJ(17)   asin(-1/mu) for |mu| >= 1, -90 otherwise
*      PRJ(18)   mu*cos(gamma)
*      PRJ(19)   1 if |mu*cos(gamma)| < 1, 0 otherwise
*=======================================================================
      SUBROUTINE AZPSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      DOUBLE PRECISION ASIND, COSD, SIND
*-----------------------------------------------------------------------
      PRJ(11) = SIGN(101.0, PRJ(11))

      IF (PRJ(10).EQ.0D0) PRJ(10) = R2D

      PRJ(12) = PRJ(10)*(PRJ(1) + 1D0)
      IF (PRJ(12).EQ.0D0) THEN
         IERR = 1
         RETURN
      END IF

      PRJ(15) = COSD(PRJ(2))
      IF (PRJ(15).EQ.0D0) THEN
         IERR = 1
         RETURN
      END IF

      PRJ(14) = 1D0/PRJ(15)
      PRJ(16) = SIND(PRJ(2))
      PRJ(13) = PRJ(16)/PRJ(15)

      IF (ABS(PRJ(1)).GT.1D0) THEN
         PRJ(17) = ASIND(-1D0/PRJ(1))
      ELSE
         PRJ(17) = -90D0
      END IF

      PRJ(18) = PRJ(1)*PRJ(15)

      IF (ABS(PRJ(18)).LT.1D0) THEN
         PRJ(19) = 1D0
      ELSE
         PRJ(19) = 0D0
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE AZPFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION A, B, CPHI, CTHE, PHI, PRJ(0:20), R, S, T, THETA,
     *          X, Y

      DOUBLE PRECISION ASIND, ATAND, COSD, SIND
*-----------------------------------------------------------------------
      IF (ABS(PRJ(11)).NE.101D0) THEN
         CALL AZPSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      CPHI = COSD(PHI)
      CTHE = COSD(THETA)

      S = PRJ(13)*CPHI
      T = (PRJ(1) + SIND(THETA)) + CTHE*S
      IF (T.EQ.0D0) THEN
         IERR = 2
         RETURN
      END IF

      R =  PRJ(12)*CTHE/T
      X =  R*SIND(PHI)
      Y = -R*CPHI*PRJ(14)

*     Bounds checking.
      IF (PRJ(11).GT.0D0) THEN
*        Overlap.
         IF (THETA.LT.PRJ(17)) THEN
            IERR = 2
            RETURN
         END IF

*        Divergence.
         IF (PRJ(19).GT.0D0) THEN
            T = PRJ(1)/SQRT(1d0 + S*S)

            IF (ABS(T).LE.1D0) THEN
               S = ATAND(-S)
               T = ASIND(T)
               A = S - T
               B = S + T + 180D0

               IF (A.GT.90D0) A = A - 360D0
               IF (B.GT.90D0) B = B - 360D0

               IF (THETA.LT.MAX(A,B)) THEN
                  IERR = 2
                  RETURN
               END IF
            END IF
         END IF
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE AZPREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION A, B, PHI, PRJ(0:20), R, S, T, THETA, TOL, X, Y,
     *          YCOSG

      PARAMETER (TOL = 1D-13)

      DOUBLE PRECISION ASIND, ATAN2D
*-----------------------------------------------------------------------
      IF (ABS(PRJ(11)).NE.101D0) THEN
         CALL AZPSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      YCOSG = Y*PRJ(15)

      R = SQRT(X*X + YCOSG*YCOSG)
      IF (R.EQ.0D0) THEN
         PHI   =  0D0
         THETA = 90D0
      ELSE
         PHI = ATAN2D(X, -YCOSG)

         S = R/(PRJ(12) + Y*PRJ(16))
         T = S*PRJ(1)/SQRT(S*S + 1D0)

         S = ATAN2D(1D0, S)

         IF (ABS(T).GT.1D0) THEN
            T = SIGN(90D0, T)
            IF (ABS(T).GT.1D0+TOL) THEN
               IERR = 2
               RETURN
            END IF
         ELSE
            T = ASIND(T)
         END IF

         A = S - T
         B = S + T + 180D0

         IF (A.GT.90D0) A = A - 360D0
         IF (B.GT.90D0) B = A - 360D0

         THETA = MAX(A, B)
      END IF

      IERR = 0
      RETURN
      END

*=======================================================================
*   SZP: slant zenithal perspective projection.
*
*   Given:
*      PRJ(1)    Distance of the point of projection from the centre of
*                the generating sphere, mu in units of r0.
*      PRJ(2)    Native longitude, phi_c, and ...
*      PRJ(3)    Native latitude, theta_c, on the planewards side of
*                the intersection of the line through the point of
*                projection and the centre of the generating sphere,
*                phi_c in degrees.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*      PRJ(11)   State flag, 102 or -102 if given < 0.
*
*   Returned:
*      PRJ(12)   1/r0
*      PRJ(13)   xp = -mu*cos(theta_c)*sin(phi_c)
*      PRJ(14)   yp =  mu*cos(theta_c)*cos(phi_c)
*      PRJ(15)   zp =  mu*sin(theta_c) + 1
*      PRJ(16)   r0*xp
*      PRJ(17)   r0*yp
*      PRJ(18)   r0*zp
*      PRJ(19)   (zp - 1)^2
*      PRJ(20)   asin(1-zp) if |1 - zp| < 1, -90 otherwise
*=======================================================================
      SUBROUTINE SZPSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      DOUBLE PRECISION ASIND, COSD, SIND
*-----------------------------------------------------------------------
      PRJ(11) = SIGN(102.0, PRJ(11))

      IF (PRJ(10).EQ.0D0) PRJ(10) = R2D

      PRJ(12) = 1D0/PRJ(10)

      PRJ(15) = PRJ(1)*SIND(PRJ(3)) + 1D0
      IF (PRJ(15).EQ.0D0) THEN
         IERR = 1
         RETURN
      END IF

      PRJ(13) =  -PRJ(1)*COSD(PRJ(3))*SIND(PRJ(2))
      PRJ(14) =   PRJ(1)*COSD(PRJ(3))*COSD(PRJ(2))
      PRJ(16) =  PRJ(10)*PRJ(13)
      PRJ(17) =  PRJ(10)*PRJ(14)
      PRJ(18) =  PRJ(10)*PRJ(15)
      PRJ(19) = (PRJ(15) - 1D0)**2

      IF (ABS(PRJ(15)-1D0).LT.1D0) THEN
         PRJ(20) = ASIND(1D0 - PRJ(15))
      ELSE
         PRJ(20) = -90D0
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE SZPFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION A, B, CPHI, CTHE, PHI, PRJ(0:20), S, SPHI, T,
     *          THETA, X, Y

      DOUBLE PRECISION ASIND, ATAN2D, COSD, SIND
*-----------------------------------------------------------------------
      IF (ABS(PRJ(11)).NE.102D0) THEN
         CALL SZPSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      CPHI = COSD(PHI)
      SPHI = SIND(PHI)
      CTHE = COSD(THETA)
      S = 1D0 - SIND(THETA)

      T = PRJ(15) - S
      IF (T.EQ.0D0) THEN
         IERR = 2
         RETURN
      END IF

      X =  (PRJ(18)*CTHE*SPHI - PRJ(16)*S)/T
      Y = -(PRJ(18)*CTHE*CPHI + PRJ(17)*S)/T

*     Bounds checking.
      IF (PRJ(11).GT.0D0) THEN
*        Divergence.
         IF (THETA.LT.PRJ(20)) THEN
            IERR = 2
            RETURN
         END IF

*        Overlap.
         IF (ABS(PRJ(1)).GT.1D0) THEN
            S = PRJ(13)*SPHI - PRJ(14)*CPHI
            T = 1D0/SQRT(PRJ(19) + S*S)

            IF (ABS(T).LE.1D0) THEN
               S = ATAN2D(S, PRJ(15)-1D0)
               T = ASIND(T)
               A = S - T
               B = S + T + 180D0

               IF (A.GT.90D0) A = A - 360D0
               IF (B.GT.90D0) B = B - 360D0

               IF (THETA.LT.MAX(A,B)) THEN
                  IERR = 2
                  RETURN
               END IF
            END IF
         END IF
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE SZPREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION A, B, C, D, PHI, PRJ(0:20), R2, STH1, STH2, STHE,
     *          SXY, T, THETA, TOL, X, X1, XP, Y, Y1, YP, Z

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      PARAMETER (TOL = 1D-13)

      DOUBLE PRECISION ASIND, ATAN2D
*-----------------------------------------------------------------------
      IF (ABS(PRJ(11)).NE.102D0) THEN
         CALL SZPSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      XP = X*PRJ(12)
      YP = Y*PRJ(12)
      R2 = XP*XP + YP*YP

      X1 = (XP - PRJ(13))/PRJ(15)
      Y1 = (YP - PRJ(14))/PRJ(15)
      SXY = XP*X1 + YP*Y1

      IF (R2.LT.1D-10) THEN
*        Use small angle formula.
         Z = R2/2D0
         THETA = 90D0 - R2D*SQRT(R2/(1D0 + SXY))

      ELSE
         T = X1*X1 + Y1*Y1
         A = T + 1D0
         B = SXY - T
         C = R2 - SXY - SXY + T - 1D0
         D = B*B - A*C

*        Check for a solution.
         IF (D.LT.0D0) THEN
            IERR = 2
            RETURN
         END IF
         D = SQRT(D)

*        Choose solution closest to pole.
         STH1 = (-B + D)/A
         STH2 = (-B - D)/A
         STHE = MAX(STH1, STH2)
         IF (STHE.GT.1D0) THEN
            IF (STHE-1D0.LT.TOL) THEN
               STHE = 1D0
            ELSE
               STHE = MIN(STH1, STH2)
            END IF
         END IF

         IF (STHE.LT.-1D0) THEN
            IF (STHE+1D0.GT.-TOL) THEN
               STHE = -1D0
            END IF
         END IF

         IF (STHE.GT.1D0 .OR. STHE.LT.-1D0) THEN
            IERR = 2
            RETURN
         END IF

         THETA = ASIND(STHE)

         Z = 1D0 - STHE
      END IF

      PHI = ATAN2D(XP - X1*Z, -(YP - Y1*Z))

      IERR = 0
      RETURN
      END

*=======================================================================
*   TAN: gnomonic projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*      PRJ(11)   State flag, 103 or -103 if given < 0.
*=======================================================================
      SUBROUTINE TANSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      PRJ(11) = SIGN(103.0, PRJ(11))

      IF (PRJ(10).EQ.0D0) PRJ(10) = R2D

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE TANFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), R, S, THETA, X, Y

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (ABS(PRJ(11)).NE.103D0) THEN
         CALL TANSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      S = SIND(THETA)
      IF (S.EQ.0D0) THEN
         IERR = 2
         RETURN
      END IF

      R =  PRJ(10)*COSD(THETA)/S
      X =  R*SIND(PHI)
      Y = -R*COSD(PHI)

      IF (PRJ(11).GT.0D0 .AND. S.LT.0D0) THEN
         IERR = 2
         RETURN
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE TANREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), R, THETA, X, Y

      DOUBLE PRECISION ATAN2D
*-----------------------------------------------------------------------
      IF (ABS(PRJ(11)).NE.103D0) THEN
         CALL TANSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      R = SQRT(X*X + Y*Y)
      IF (R.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = ATAN2D(X, -Y)
      END IF
      THETA = ATAN2D(PRJ(10), R)

      IERR = 0
      RETURN
      END

*=======================================================================
*   STG: stereographic projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 104.
*      PRJ(12)   2*r0
*      PRJ(13)   1/(2*r0)
*=======================================================================
      SUBROUTINE STGSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      PRJ(11) = 104D0

      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D
         PRJ(12) = 360D0/PI
         PRJ(13) = PI/360D0
      ELSE
         PRJ(12) = 2D0*PRJ(10)
         PRJ(13) = 1D0/PRJ(12)
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE STGFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), R, S, THETA, X, Y

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.104D0) THEN
         CALL STGSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      S = 1D0 + SIND(THETA)
      IF (S.EQ.0D0) THEN
         IERR = 2
         RETURN
      END IF

      R =  PRJ(12)*COSD(THETA)/S
      X =  R*SIND(PHI)
      Y = -R*COSD(PHI)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE STGREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), R, THETA, X, Y

      DOUBLE PRECISION ATAN2D, ATAND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.104D0) THEN
         CALL STGSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      R = SQRT(X*X + Y*Y)
      IF (R.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = ATAN2D(X, -Y)
      END IF
      THETA = 90D0 - 2D0*ATAND(R*PRJ(13))

      IERR = 0
      RETURN
      END

*=======================================================================
*   SIN: orthographic/synthesis projection.
*
*   Given:
*      PRJ(1:2)  SIN obliqueness parameters, xi and eta.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*      PRJ(11)   State flag, 105 or -105 if given < 0.
*
*   Returned:
*      PRJ(12)   1/r0
*      PRJ(13)   xi**2 + eta**2
*      PRJ(14)   xi**2 + eta**2 + 1
*      PRJ(15)   xi**2 + eta**2 - 1
*=======================================================================
      SUBROUTINE SINSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      PRJ(11) = SIGN(105.0, PRJ(11))

      IF (PRJ(10).EQ.0D0) PRJ(10) = R2D

      PRJ(12) = 1D0/PRJ(10)
      PRJ(13) = PRJ(1)*PRJ(1) + PRJ(2)*PRJ(2)
      PRJ(14) = PRJ(13) + 1D0
      PRJ(15) = PRJ(13) - 1D0

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE SINFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION CPHI, CTHE, PHI, PRJ(0:20), SPHI, T, THETA, X, Y,
     *          Z

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      DOUBLE PRECISION ATAND, COSD, SIND
*-----------------------------------------------------------------------
      IF (ABS(PRJ(11)).NE.105D0) THEN
         CALL SINSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      T = (90D0 - ABS(THETA))*D2R
      IF (T.LT.1D-4) THEN
         IF (THETA.GT.0D0) THEN
            Z = T*T/2D0
         ELSE
            Z = 2D0 - T*T/2D0
         END IF
         CTHE = T
      ELSE
         Z = 1D0 - SIND(THETA)
         CTHE = COSD(THETA)
      END IF

      CPHI = COSD(PHI)
      SPHI = SIND(PHI)
      X =  PRJ(10)*(CTHE*SPHI + PRJ(1)*Z)
      Y = -PRJ(10)*(CTHE*CPHI - PRJ(2)*Z)

*     Validate this solution.
      IF (PRJ(11).GT.0D0) THEN
         IF (PRJ(13).EQ.0D0) THEN
*           Orthographic projection.
            IF (THETA.LT.0D0) THEN
               IERR = 2
               RETURN
            END IF
         ELSE
*           "Synthesis" projection.
            T = -ATAND(PRJ(1)*SPHI - PRJ(2)*CPHI)
            IF (THETA.LT.T) THEN
               IERR = 2
               RETURN
            END IF
         END IF
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE SINREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION A, B, C, D, PHI, PRJ(0:20), R2, STH1, STH2, STHE,
     *          SXY, THETA, TOL, X, X0, X1, XP, Y, Y0, Y1, YP, Z

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      PARAMETER (TOL = 1D-13)

      DOUBLE PRECISION ACOSD, ASIND, ATAN2D
*-----------------------------------------------------------------------
      IF (ABS(PRJ(11)).NE.105D0) THEN
         CALL SINSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

*     Compute intermediaries.
      X0 = X*PRJ(12)
      Y0 = Y*PRJ(12)
      R2 = X0*X0 + Y0*Y0

      IF (PRJ(13).EQ.0D0) THEN
*        Orthographic projection.
         IF (R2.NE.0D0) THEN
            PHI = ATAN2D(X0, -Y0)
         ELSE
            PHI = 0D0
         END IF

         IF (R2.LT.0.5D0) THEN
            THETA = ACOSD(SQRT(R2))
         ELSE IF (R2.LE.1D0) THEN
            THETA = ASIND(SQRT(1D0 - R2))
         ELSE
            IERR = 2
            RETURN
         END IF

      ELSE
*        "Synthesis" projection.
         X1 = PRJ(1)
         Y1 = PRJ(2)
         SXY = X0*X1 + Y0*Y1

         IF (R2.LT.1D-10) THEN
*           Use small angle formula.
            Z = R2/2D0
            THETA = 90D0 - R2D*SQRT(R2/(1D0 + SXY))

         ELSE
            A = PRJ(14)
            B = SXY - PRJ(13)
            C = R2 - SXY - SXY + PRJ(15)
            D = B*B - A*C

*           Check for a solution.
            IF (D.LT.0D0) THEN
               IERR = 2
               RETURN
            END IF
            D = SQRT(D)

*           Choose solution closest to pole.
            STH1 = (-B + D)/A
            STH2 = (-B - D)/A
            STHE = MAX(STH1,STH2)
            IF (STHE.GT.1D0) THEN
               IF (STHE-1D0.LT.TOL) THEN
                  STHE = 1D0
               ELSE
                  STHE = MIN(STH1,STH2)
               END IF
            END IF

            IF (STHE.LT.-1D0) THEN
               IF (STHE+1D0.GT.-TOL) THEN
                  STHE = -1D0
               END IF
            END IF

            IF (STHE.GT.1D0 .OR. STHE.LT.-1D0) THEN
               IERR = 2
               RETURN
            END IF

            THETA = ASIND(STHE)
            Z = 1D0 - STHE
         END IF

         XP = -Y0 + PRJ(2)*Z
         YP =  X0 - PRJ(1)*Z
         IF (XP.EQ.0D0 .AND. YP.EQ.0D0) THEN
            PHI = 0D0
         ELSE
            PHI = ATAN2D(YP,XP)
         END IF
      END IF

      IERR = 0
      RETURN
      END

*=======================================================================
*   ARC: zenithal/azimuthal equidistant projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 106.
*      PRJ(12)   r0*(pi/180)
*      PRJ(13)   (180/pi)/r0
*=======================================================================
      SUBROUTINE ARCSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      PRJ(11) = 106D0

      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D
         PRJ(12) = 1D0
         PRJ(13) = 1D0
      ELSE
         PRJ(12) = PRJ(10)*D2R
         PRJ(13) = 1D0/PRJ(12)
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE ARCFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), R, THETA, X, Y

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.106D0) THEN
         CALL ARCSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      R =  PRJ(12)*(90D0 - THETA)
      X =  R*SIND(PHI)
      Y = -R*COSD(PHI)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE ARCREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), R, THETA, X, Y

      DOUBLE PRECISION ATAN2D
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.106D0) THEN
         CALL ARCSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      R = SQRT(X*X + Y*Y)
      IF (R.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = ATAN2D(X, -Y)
      END IF
      THETA = 90D0 - R*PRJ(13)

      IERR = 0
      RETURN
      END

*=======================================================================
*   ZPN: zenithal/azimuthal polynomial projection.
*
*   Given:
*      PRJ(0:9)  Polynomial coefficients.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*      PRJ(11)   State flag, 107 or -107 if given < 0.
*
*   Returned:
*      PRJ(12)   Degree of the polynomial, N.
*      PRJ(13)   Co-latitude of the first point of inflection (N > 2).
*      PRJ(14)   Radius of the first point of inflection (N > 2).
*=======================================================================
      SUBROUTINE ZPNSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR, I, J, K
      DOUBLE PRECISION D, D1, D2, PRJ(0:20), R, TOL, ZD, ZD1, ZD2

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      PARAMETER (TOL = 1D-13)
*-----------------------------------------------------------------------
      PRJ(11) = SIGN(107.0, PRJ(11))

      IF (PRJ(10).EQ.0D0) PRJ(10) = R2D

*     Find the highest non-zero coefficient.
      DO 10 K = 9, 0, -1
         IF (PRJ(K).NE.0D0) GO TO 20
 10   CONTINUE
      IERR = 1
      RETURN

 20   PRJ(12) = DBLE(K)

      IF (K.GE.3) THEN
*        Find the point of inflection closest to the pole.
         ZD1 = 0D0
         D1  = PRJ(1)
         IF (D1.LE.0D0) THEN
            IERR = 1
            RETURN
         END IF

*        Find the point where the derivative first goes negative.
         DO 40 I = 1, 180
            ZD2 = PI*DBLE(I)/180D0
            D2  = 0D0
            DO 30 J = K, 1, -1
               D2 = D2*ZD2 + J*PRJ(J)
 30         CONTINUE

            IF (D2.LE.0D0) GO TO 50
            ZD1 = ZD2
            D1  = D2
 40      CONTINUE

*        No negative derivative -> no point of inflection.
         ZD = PI
         GO TO 80

*        Find where the derivative is zero.
 50      DO 70 I = 1, 10
            ZD = ZD1 - D1*(ZD2-ZD1)/(D2-D1)

            D = 0D0
            DO 60 J = K, 1, -1
               D = D*ZD + J*PRJ(J)
 60         CONTINUE

            IF (ABS(D).LT.TOL) GO TO 80

            IF (D.LT.0D0) THEN
               ZD2 = ZD
               D2  = D
            ELSE
               ZD1 = ZD
               D1  = D
            END IF
 70      CONTINUE

 80      R = 0D0
         DO 90 J = K, 0, -1
            R = R*ZD + PRJ(J)
 90      CONTINUE
         PRJ(13) = ZD
         PRJ(14) = R
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE ZPNFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR, J
      DOUBLE PRECISION PHI, PRJ(0:20), R, S, THETA, X, Y

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (ABS(PRJ(11)).NE.107D0) THEN
         CALL ZPNSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      S = (90D0 - THETA)*D2R

      R = 0D0
      DO 10 J = 9, 0, -1
         R = R*S + PRJ(J)
 10   CONTINUE
      R = PRJ(10)*R

      X =  R*SIND(PHI)
      Y = -R*COSD(PHI)

      IF (PRJ(11).GT.0D0 .AND. S.GT.PRJ(13)) THEN
         IERR = 2
         RETURN
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE ZPNREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR, I, J, K
      DOUBLE PRECISION A, B, C, D, LAMBDA, PHI, PRJ(0:20), R, R1, R2,
     *          RT, THETA, TOL, X, Y, ZD, ZD1, ZD2

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      PARAMETER (TOL = 1D-13)

      DOUBLE PRECISION ATAN2D
*-----------------------------------------------------------------------
      IF (ABS(PRJ(11)).NE.107D0) THEN
         CALL ZPNSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      K = NINT(PRJ(12))

      R = SQRT(X*X + Y*Y)/PRJ(10)

      IF (K.LT.1) THEN
*        Constant - no solution.
         IERR = 1
         RETURN
      ELSE IF (K.EQ.1) THEN
*        Linear.
         ZD = (R - PRJ(0))/PRJ(1)
      ELSE IF (K.EQ.2) THEN
*        Quadratic.
         A = PRJ(2)
         B = PRJ(1)
         C = PRJ(0) - R

         D = B*B - 4.0*A*C
         IF (D.LT.0D0) THEN
            IERR = 2
            RETURN
         END IF
         D = SQRT(D)

*        Choose solution closest to pole.
         ZD1 = (-B + D)/(2D0*A)
         ZD2 = (-B - D)/(2D0*A)
         ZD  = MIN(ZD1,ZD2)
         IF (ZD.LT.-TOL) ZD = MAX(ZD1,ZD2)
         IF (ZD.LT.0D0) THEN
            IF (ZD.LT.-TOL) THEN
               IERR = 2
               RETURN
            END IF
            ZD = 0D0
         ELSE IF (ZD.GT.PI) THEN
            IF (ZD.GT.PI+TOL) THEN
               IERR = 2
               RETURN
            END IF
            ZD = PI
         END IF
      ELSE
*        Higher order - solve iteratively.
         ZD1 = 0D0
         R1  = PRJ(0)
         ZD2 = PRJ(13)
         R2  = PRJ(14)

         IF (R.LT.R1) THEN
            IF (R.LT.R1-TOL) THEN
               IERR = 2
               RETURN
            END IF
            ZD = ZD1
            GO TO 120
         ELSE IF (R.GT.R2) THEN
            IF (R.GT.R2+TOL) THEN
               IERR = 2
               RETURN
            END IF
            ZD = ZD2
            GO TO 120
         ELSE
*           Disect the interval.
            DO 110 J = 1, 100
               LAMBDA = (R2 - R)/(R2 - R1)
               IF (LAMBDA.LT.0.1D0) THEN
                  LAMBDA = 0.1D0
               ELSE IF (LAMBDA.GT.0.9D0) THEN
                  LAMBDA = 0.9D0
               END IF

               ZD = ZD2 - LAMBDA*(ZD2 - ZD1)

               RT = 0D0
               DO 100 I = K, 0, -1
                   RT = (RT * ZD) + PRJ(I)
 100           CONTINUE

               IF (RT.LT.R) THEN
                   IF (R-RT.LT.TOL) GO TO 120
                   R1 = RT
                   ZD1 = ZD
               ELSE
                   IF (RT-R.LT.TOL) GO TO 120
                   R2 = RT
                   ZD2 = ZD
               END IF

               IF (ABS(ZD2-ZD1).LT.TOL) GO TO 120
 110        CONTINUE
         END IF
      END IF

 120  IF (R.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = ATAN2D(X, -Y)
      END IF
      THETA = 90D0 - ZD*R2D

      IERR = 0
      RETURN
      END

*=======================================================================
*   ZEA: zenithal/azimuthal equal area projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 108.
*      PRJ(12)   2*r0
*      PRJ(13)   1/(2*r0)
*=======================================================================
      SUBROUTINE ZEASET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      PRJ(11) = 108D0

      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D
         PRJ(12) = 360D0/PI
         PRJ(13) = PI/360D0
      ELSE
         PRJ(12) = 2D0*PRJ(10)
         PRJ(13) = 1D0/PRJ(12)
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE ZEAFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), R, THETA, X, Y

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.108D0) THEN
         CALL ZEASET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      R =  PRJ(12)*SIND((90D0 - THETA)/2D0)
      X =  R*SIND(PHI)
      Y = -R*COSD(PHI)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE ZEAREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), R, S, THETA, TOL, X, Y

      PARAMETER (TOL = 1D-12)

      DOUBLE PRECISION ASIND, ATAN2D
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.108D0) THEN
         CALL ZEASET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      R = SQRT(X*X + Y*Y)
      IF (R.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = ATAN2D(X, -Y)
      END IF

      S = R*PRJ(13)
      IF (ABS(S).GT.1D0) THEN
         IF (ABS(R-PRJ(12)).LT.TOL) THEN
            THETA = -90D0
         ELSE
            IERR = 2
            RETURN
         END IF
      ELSE
         THETA = 90D0 - 2D0*ASIND(R*PRJ(13))
      END IF

      IERR = 0
      RETURN
      END

*=======================================================================
*   AIR: Airy's projection.
*
*   Given:
*      PRJ(1)    Latitude theta_b within which the error is minimized,
*                in degrees.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 109.
*      PRJ(12)   2*r0
*      PRJ(13)   ln(cos(xi_b))/tan(xi_b)**2, where xi_b = (90-theta_b)/2
*      PRJ(14)   1/2 - PRJ(13)
*      PRJ(15)   2*r0*PRJ(14)
*      PRJ(16)   tol, cutoff for using small angle approximation, in
*                radians.
*      PRJ(17)   PRJ(14)*tol
*      PRJ(18)   (180/pi)/PRJ(14)
*=======================================================================
      SUBROUTINE AIRSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION CXI, PRJ(0:20), TOL

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      PARAMETER (TOL = 1D-4)

      DOUBLE PRECISION COSD
*-----------------------------------------------------------------------
      PRJ(11) = 109D0

      IF (PRJ(10).EQ.0D0) PRJ(10) = R2D

      PRJ(12) = 2D0*PRJ(10)
      IF (PRJ(1).EQ.90D0) THEN
         PRJ(13) = -0.5D0
         PRJ(14) = 1D0
      ELSE IF (PRJ(1).GT.-90D0) THEN
         CXI = COSD((90D0 - PRJ(1))/2D0)
         PRJ(13) = LOG(CXI)*(CXI*CXI)/(1D0-CXI*CXI)
         PRJ(14) = 0.5D0 - PRJ(13)
      ELSE
         IERR = 1
         RETURN
      END IF

      PRJ(15) = PRJ(12)*PRJ(14)
      PRJ(16) = TOL
      PRJ(17) = PRJ(14)*TOL
      PRJ(18) = R2D/PRJ(14)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE AIRFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION CXI, PHI, PRJ(0:20), R, THETA, TXI, X, XI, Y

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.109D0) THEN
         CALL AIRSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      IF (THETA.EQ.90D0) THEN
         R = 0D0
      ELSE IF (THETA.GT.-90D0) THEN
         XI = D2R*(90D0 - THETA)/2D0
         IF (XI.LT.PRJ(16)) THEN
            R = XI*PRJ(15)
         ELSE
            CXI = COS(XI)
            TXI = SQRT(1D0-CXI*CXI)/CXI
            R = -PRJ(12)*(LOG(CXI)/TXI + PRJ(13)*TXI)
         END IF
      ELSE
         IERR = 2
         RETURN
      END IF

      X =  R*SIND(PHI)
      Y = -R*COSD(PHI)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE AIRREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR, J
      DOUBLE PRECISION CXI, LAMBDA, PHI, PRJ(0:20), R, R1, R2, RT,
     *          THETA, TOL, TXI, X, X1, X2, XI, Y

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      PARAMETER (TOL = 1D-12)

      DOUBLE PRECISION ACOSD, ATAN2D
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.109D0) THEN
         CALL AIRSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      R = SQRT(X*X + Y*Y)/PRJ(12)

      IF (R.EQ.0D0) THEN
         XI = 0D0
      ELSE IF (R.LT.PRJ(17)) THEN
         XI = R*PRJ(18)
      ELSE
*        Find a solution interval.
         X1 = 1D0
         R1 = 0D0
         DO 10 J = 1, 30
            X2 = X1/2D0
            TXI = SQRT(1D0-X2*X2)/X2
            R2 = -(LOG(X2)/TXI + PRJ(13)*TXI)

            IF (R2.GE.R) GO TO 20
            X1 = X2
            R1 = R2
 10      CONTINUE
         IERR = 2
         RETURN

 20      DO 30 J = 1, 100
*           Weighted division of the interval.
            LAMBDA = (R2-R)/(R2-R1)
            IF (LAMBDA.LT.0.1D0) THEN
               LAMBDA = 0.1D0
            ELSE IF (LAMBDA.GT.0.9D0) THEN
               LAMBDA = 0.9D0
            END IF
            CXI = X2 - LAMBDA*(X2-X1)

            TXI = SQRT(1D0-CXI*CXI)/CXI
            RT = -(LOG(CXI)/TXI + PRJ(13)*TXI)

            IF (RT.LT.R) THEN
               IF (R-RT.LT.TOL) GO TO 40
               R1 = RT
               X1 = CXI
            ELSE
               IF (RT-R.LT.TOL) GO TO 40
               R2 = RT
               X2 = CXI
            END IF
 30      CONTINUE
         IERR = 2
         RETURN

 40      XI = ACOSD(CXI)
      END IF

      IF (R.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = ATAN2D(X, -Y)
      END IF
      THETA = 90D0 - 2D0*XI

      IERR = 0
      RETURN
      END

*=======================================================================
*   CYP: cylindrical perspective projection.
*
*   Given:
*      PRJ(1)    Distance of point of projection from the centre of the
*                generating sphere, mu, in units of r0.
*      PRJ(2)    Radius of the cylinder of projection, lambda, in units
*                of r0.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 201.
*      PRJ(12)   r0*lambda*(pi/180)
*      PRJ(13)   (180/pi)/(r0*lambda)
*      PRJ(14)   r0*(mu + lambda)
*      PRJ(15)   1/(r0*(mu + lambda))
*=======================================================================
      SUBROUTINE CYPSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      PRJ(11) = 201D0

      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D

         PRJ(12) = PRJ(2)
         IF (PRJ(12).EQ.0D0) THEN
            IERR = 1
            RETURN
         END IF

         PRJ(13) = 1D0/PRJ(12)

         PRJ(14) = R2D*(PRJ(1) + PRJ(2))
         IF (PRJ(14).EQ.0D0) THEN
            IERR = 1
            RETURN
         END IF

         PRJ(15) = 1D0/PRJ(14)
      ELSE
         PRJ(12) = PRJ(10)*PRJ(2)*D2R
         IF (PRJ(12).EQ.0D0) THEN
            IERR = 1
            RETURN
         END IF

         PRJ(13) = 1D0/PRJ(12)

         PRJ(14) = PRJ(10)*(PRJ(1) + PRJ(2))
         IF (PRJ(14).EQ.0D0) THEN
            IERR = 1
            RETURN
         END IF

         PRJ(15) = 1D0/PRJ(14)
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE CYPFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), S, THETA, X, Y

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.201D0) THEN
         CALL CYPSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      S = PRJ(1) + COSD(THETA)
      IF (S.EQ.0D0) THEN
            IERR = 2
            RETURN
         END IF

      X = PRJ(12)*PHI
      Y = PRJ(14)*SIND(THETA)/S

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE CYPREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION ETA, PHI, PRJ(0:20), THETA, X, Y

      DOUBLE PRECISION ASIND, ATAN2D
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.201D0) THEN
         CALL CYPSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      PHI   = X*PRJ(13)
      ETA   = Y*PRJ(15)
      THETA = ATAN2D(ETA,1D0) + ASIND(ETA*PRJ(1)/SQRT(ETA*ETA+1D0))

      IERR = 0
      RETURN
      END

*=======================================================================
*   CEA: cylindrical equal area projection.
*
*   Given:
*      PRJ(1)    Square of the cosine of the latitude at which the
*                projection is conformal, lambda.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 202.
*      PRJ(12)   r0*(pi/180)
*      PRJ(13)   (180/pi)/r0
*      PRJ(14)   r0/lambda
*      PRJ(15)   lambda/r0
*=======================================================================
      SUBROUTINE CEASET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      PRJ(11) = 202D0

      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D
         PRJ(12) = 1D0
         PRJ(13) = 1D0
         IF (PRJ(1).LE.0D0 .OR. PRJ(1).GT.1D0) THEN
            IERR = 1
            RETURN
         END IF
         PRJ(14) = PRJ(10)/PRJ(1)
         PRJ(15) = PRJ(1)/PRJ(10)
      ELSE
         PRJ(12) = PRJ(10)*D2R
         PRJ(13) = R2D/PRJ(10)
         IF (PRJ(1).LE.0D0 .OR. PRJ(1).GT.1D0) THEN
            IERR = 1
            RETURN
         END IF
         PRJ(14) = PRJ(10)/PRJ(1)
         PRJ(15) = PRJ(1)/PRJ(10)
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE CEAFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), THETA, X, Y

      DOUBLE PRECISION SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.202D0) THEN
         CALL CEASET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      X = PRJ(12)*PHI
      Y = PRJ(14)*SIND(THETA)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE CEAREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), S, THETA, TOL, X, Y

      PARAMETER (TOL = 1D-13)

      DOUBLE PRECISION ASIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.202D0) THEN
         CALL CEASET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      S = Y*PRJ(15)
      IF (ABS(S).GT.1D0) THEN
         IF (ABS(S).GT.1D0+TOL) THEN
            IERR = 2
            RETURN
         END IF

         S = SIGN(1D0,S)
      END IF

      PHI   = X*PRJ(13)
      THETA = ASIND(S)

      IERR = 0
      RETURN
      END

*=======================================================================
*   CAR: Cartesian projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 203.
*      PRJ(12)   r0*(pi/180)
*      PRJ(13)   (180/pi)/r0
*=======================================================================
      SUBROUTINE CARSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      PRJ(11) = 203D0

      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D
         PRJ(12) = 1D0
         PRJ(13) = 1D0
      ELSE
         PRJ(12) = PRJ(10)*D2R
         PRJ(13) = 1D0/PRJ(12)
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE CARFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), THETA, X, Y
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.203D0) THEN
         CALL CARSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      X = PRJ(12)*PHI
      Y = PRJ(12)*THETA

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE CARREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), THETA, X, Y
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.203D0) THEN
         CALL CARSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      PHI   = PRJ(13)*X
      THETA = PRJ(13)*Y

      IERR = 0
      RETURN
      END

*=======================================================================
*   MER: Mercator's projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 204.
*      PRJ(12)   r0*(pi/180)
*      PRJ(13)   (180/pi)/r0
*=======================================================================
      SUBROUTINE MERSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      PRJ(11) = 204D0

      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D
         PRJ(12) = 1D0
         PRJ(13) = 1D0
      ELSE
         PRJ(12) = PRJ(10)*D2R
         PRJ(13) = 1D0/PRJ(12)
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE MERFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), THETA, X, Y

      DOUBLE PRECISION TAND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.204D0) THEN
         CALL MERSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      IF (THETA.LE.-90D0 .OR. THETA.GE.90D0) THEN
         IERR = 2
         RETURN
      END IF

      X = PRJ(12)*PHI
      Y = PRJ(10)*LOG(TAND((90D0+THETA)/2D0))

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE MERREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), THETA, X, Y

      DOUBLE PRECISION ATAND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.204D0) THEN
         CALL MERSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      PHI   = X*PRJ(13)
      THETA = 2D0*ATAND(EXP(Y/PRJ(10))) - 90D0

      IERR = 0
      RETURN
      END

*=======================================================================
*   SFL: Sanson-Flamsteed ("global sinusoid") projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 301.
*      PRJ(12)   r0*(pi/180)
*      PRJ(13)   (180/pi)/r0
*=======================================================================
      SUBROUTINE SFLSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      PRJ(11) = 301D0

      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D
         PRJ(12) = 1D0
         PRJ(13) = 1D0
      ELSE
         PRJ(12) = PRJ(10)*D2R
         PRJ(13) = 1D0/PRJ(12)
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE SFLFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), THETA, X, Y

      DOUBLE PRECISION COSD
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.301D0) THEN
         CALL SFLSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      X = PRJ(12)*PHI*COSD(THETA)
      Y = PRJ(12)*THETA

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE SFLREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), THETA, W, X, Y
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.301D0) THEN
         CALL SFLSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      W = COS(Y/PRJ(10))
      IF (W.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = X*PRJ(13)/COS(Y/PRJ(10))
      END IF
      THETA = Y*PRJ(13)

      IERR = 0
      RETURN
      END

*=======================================================================
*   PAR: parabolic projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 302.
*      PRJ(12)   r0*(pi/180)
*      PRJ(13)   (180/pi)/r0
*      PRJ(14)   pi*r0
*      PRJ(15)   1/(pi*r0)
*=======================================================================
      SUBROUTINE PARSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      PRJ(11) = 302D0

      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D
         PRJ(12) = 1D0
         PRJ(13) = 1D0
         PRJ(14) = 180D0
         PRJ(15) = 1D0/PRJ(14)
      ELSE
         PRJ(12) = PRJ(10)*D2R
         PRJ(13) = 1D0/PRJ(12)
         PRJ(14) = PI*PRJ(10)
         PRJ(15) = 1D0/PRJ(14)
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE PARFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), S, THETA, X, Y

      DOUBLE PRECISION SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.302D0) THEN
         CALL PARSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      S = SIND(THETA/3D0)
      X = PRJ(12)*PHI*(1D0 - 4D0*S*S)
      Y = PRJ(14)*S

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE PARREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), S, T, THETA, X, Y

      DOUBLE PRECISION ASIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.302D0) THEN
         CALL PARSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      S = Y*PRJ(15)
      IF (S.GT.1D0 .OR. S.LT.-1D0) THEN
         IERR = 2
         RETURN
      END IF

      T = 1D0 - 4D0*S*S
      IF (T.EQ.0D0) THEN
         IF (X.EQ.0D0) THEN
            PHI = 0D0
         ELSE
            IERR = 2
            RETURN
         END IF
      ELSE
         PHI = PRJ(13)*X/T
      END IF

      THETA = 3D0*ASIND(S)

      IERR = 0
      RETURN
      END

*=======================================================================
*   MOL: Mollweide's projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 303.
*      PRJ(12)   sqrt(2)*r0
*      PRJ(13)   sqrt(2)*r0/90
*      PRJ(14)   1/(sqrt(2)*r0)
*      PRJ(15)   90/r0
*      PRJ(16)   2/pi
*=======================================================================
      SUBROUTINE MOLSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      DOUBLE PRECISION SQRT2
      PARAMETER (SQRT2 = 1.4142135623730950488D0)
*-----------------------------------------------------------------------
      PRJ(11) = 303D0

      IF (PRJ(10).EQ.0D0) PRJ(10) = R2D

      PRJ(12) = SQRT2*PRJ(10)
      PRJ(13) = PRJ(12)/90D0
      PRJ(14) = 1D0/PRJ(12)
      PRJ(15) = 90D0/PRJ(10)
      PRJ(16) = 2D0/PI

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE MOLFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR, J
      DOUBLE PRECISION GAMMA, PHI, PRJ(0:20), RESID, THETA, TOL, U,
     *          V, V0, V1, X, Y

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      PARAMETER (TOL = 1D-13)

      DOUBLE PRECISION SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.303D0) THEN
         CALL MOLSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      IF (ABS(THETA).EQ.90D0) THEN
         X = 0D0
         Y = SIGN(PRJ(12),THETA)
      ELSE IF (THETA.EQ.0D0) THEN
         X = PRJ(13)*PHI
         Y = 0D0
      ELSE
         U  = PI*SIND(THETA)
         V0 = -PI
         V1 = +PI
         V  = U
         DO 10 J = 1, 100
            RESID = (V - U) + SIN(V)
            IF (RESID.LT.0D0) THEN
               IF (RESID.GT.-TOL) GO TO 20
               V0 = V
            ELSE
               IF (RESID.LT.TOL) GO TO 20
               V1 = V
            END IF
            V = (V0 + V1)/2D0
 10      CONTINUE

 20      GAMMA = V/2D0
         X = PRJ(13)*PHI*COS(GAMMA)
         Y = PRJ(12)*SIN(GAMMA)
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE MOLREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), S, THETA, TOL, X, Y, Y0, Z

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      PARAMETER (TOL = 1D-12)

      DOUBLE PRECISION ASIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.303D0) THEN
         CALL MOLSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      Y0 = Y/PRJ(10)
      S  = 2D0 - Y0*Y0
      IF (S.LE.TOL) THEN
         IF (S.LT.-TOL) THEN
            IERR = 2
            RETURN
         END IF
         S = 0D0

         IF (ABS(X).GT.TOL) THEN
            IERR = 2
            RETURN
         END IF
         PHI = 0D0
      ELSE
         S = SQRT(S)
         PHI = PRJ(15)*X/S
      END IF

      Z = Y*PRJ(14)
      IF (ABS(Z).GT.1D0) THEN
         IF (ABS(Z).GT.1D0+TOL) THEN
            IERR = 2
            RETURN
         END IF
         Z = SIGN(1D0,Z) + Y0*S/PI
      ELSE
         Z = ASIN(Z)*PRJ(16) + Y0*S/PI
      END IF

      IF (ABS(Z).GT.1D0) THEN
         IF (ABS(Z).GT.1D0+TOL) THEN
            IERR = 2
            RETURN
         END IF
         Z = SIGN(1D0,Z)
      END IF

      THETA = ASIND(Z)

      IERR = 0
      RETURN
      END

*=======================================================================
*   AIT: Hammer-Aitoff projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 401.
*      PRJ(12)   2*r0**2
*      PRJ(13)   1/(2*r0)**2
*      PRJ(14)   1/(4*r0)**2
*      PRJ(15)   1/(2*r0)
*=======================================================================
      SUBROUTINE AITSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      PRJ(11) = 401D0

      IF (PRJ(10).EQ.0D0) PRJ(10) = R2D

      PRJ(12) = 2D0*PRJ(10)*PRJ(10)
      PRJ(13) = 1D0/(2D0*PRJ(12))
      PRJ(14) = PRJ(13)/4D0
      PRJ(15) = 1D0/(2D0*PRJ(10))

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE AITFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION CTHE, PHI, PRJ(0:20), THETA, W, X, Y

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.401D0) THEN
         CALL AITSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      CTHE = COSD(THETA)
      W = SQRT(PRJ(12)/(1D0 + CTHE*COSD(PHI/2D0)))
      X = 2D0*W*CTHE*SIND(PHI/2D0)
      Y = W*SIND(THETA)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE AITREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), S, THETA, TOL, U, X, XP, Y, YP, Z

      PARAMETER (TOL = 1D-13)

      DOUBLE PRECISION ASIND, ATAN2D
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.401D0) THEN
         CALL AITSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      U = 1D0 - X*X*PRJ(14) - Y*Y*PRJ(13)
      IF (U.LT.0D0) THEN
         IF (U.LT.-TOL) THEN
            IERR = 2
            RETURN
         END IF

         U = 0D0
      END IF

      Z = SQRT(U)
      S = Z*Y/PRJ(10)
      IF (ABS(S).GT.1D0) THEN
         IF (ABS(S).GT.1D0+TOL) THEN
            IERR = 2
            RETURN
         END IF

         S = SIGN(1D0,S)
      END IF

      XP = 2D0*Z*Z - 1D0
      YP = Z*X*PRJ(15)
      IF (XP.EQ.0D0 .AND. YP.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = 2D0*ATAN2D(YP, XP)
      END IF
      THETA = ASIND(S)

      IERR = 0
      RETURN
      END

*=======================================================================
*   COP: conic perspective projection.
*
*   Given:
*      PRJ(1)    sigma = (theta2+theta1)/2
*      PRJ(2)    delta = (theta2-theta1)/2, where theta1 and theta2 are
*                the latitudes of the standard parallels, in degrees.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*      PRJ(11)   State flag, 501 or -501 if given < 0.
*
*   Returned:
*      PRJ(12)   C  = sin(sigma)
*      PRJ(13)   1/C
*      PRJ(14)   Y0 = r0*cos(delta)*cot(sigma)
*      PRJ(15)   r0*cos(delta)
*      PRJ(16)   1/(r0*cos(delta))
*      PRJ(17)   cot(sigma)
*=======================================================================
      SUBROUTINE COPSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      DOUBLE PRECISION COSD, SIND, TAND
*-----------------------------------------------------------------------
      PRJ(11) = SIGN(501.0, PRJ(11))

      IF (PRJ(10).EQ.0D0) PRJ(10) = R2D

      PRJ(12) = SIND(PRJ(1))
      IF (PRJ(12).EQ.0D0) THEN
         IERR = 1
         RETURN
      END IF

      PRJ(13) = 1D0/PRJ(12)

      PRJ(15) = PRJ(10)*COSD(PRJ(2))
      IF (PRJ(15).EQ.0D0) THEN
         IERR = 1
         RETURN
      END IF

      PRJ(16) = 1D0/PRJ(15)
      PRJ(17) = 1D0/TAND(PRJ(1))

      PRJ(14) = PRJ(15)*PRJ(17)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE COPFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION A, PHI, PRJ(0:20), R, S, T, THETA, X, Y

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (ABS(PRJ(11)).NE.501D0) THEN
         CALL COPSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      T = THETA - PRJ(1)
      S = COSD(T)
      IF (S.EQ.0D0) THEN
         IERR = 2
         RETURN
      END IF

      A = PRJ(12)*PHI
      R = PRJ(14) - PRJ(15)*SIND(T)/S

      X =           R*SIND(A)
      Y = PRJ(14) - R*COSD(A)

      IF (PRJ(11).GT.0D0 .AND. R*PRJ(12).LT.0D0) THEN
         IERR = 2
         RETURN
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE COPREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION A, DY, PHI, PRJ(0:20), R, THETA, X, Y

      DOUBLE PRECISION ATAN2D, ATAND
*-----------------------------------------------------------------------
      IF (ABS(PRJ(11)).NE.501D0) THEN
         CALL COPSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      DY = PRJ(14) - Y
      R  = SQRT(X*X + DY*DY)
      IF (PRJ(1).LT.0D0) R = -R

      IF (R.EQ.0D0) THEN
         A = 0D0
      ELSE
         A = ATAN2D(X/R, DY/R)
      END IF

      PHI   = A*PRJ(13)
      THETA = PRJ(1) + ATAND(PRJ(17) - R*PRJ(16))

      IERR = 0
      RETURN
      END

*=======================================================================
*   COE: conic equal area projection.
*
*   Given:
*      PRJ(1)    sigma = (theta2+theta1)/2
*      PRJ(2)    delta = (theta2-theta1)/2, where theta1 and theta2 are
*                the latitudes of the standard parallels, in degrees.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 502.
*      PRJ(12)   C = (sin(theta1) + sin(theta2))/2
*      PRJ(13)   1/C
*      PRJ(14)   Y0 = chi*sqrt(psi - 2C*sind(sigma))
*      PRJ(15)   chi = r0/C
*      PRJ(16)   psi = 1 + sin(theta1)*sin(theta2)
*      PRJ(17)   2C
*      PRJ(18)   (1 + sin(theta1)*sin(theta2))*(r0/C)**2
*      PRJ(19)   C/(2*r0**2)
*      PRJ(20)   Y0 = chi*sqrt(psi + 2C)
*=======================================================================
      SUBROUTINE COESET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20), THETA1, THETA2

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      DOUBLE PRECISION SIND
*-----------------------------------------------------------------------
      PRJ(11) = 502D0

      IF (PRJ(10).EQ.0D0) PRJ(10) = R2D

      THETA1 = PRJ(1) - PRJ(2)
      THETA2 = PRJ(1) + PRJ(2)

      PRJ(12) = (SIND(THETA1) + SIND(THETA2))/2D0
      IF (PRJ(12).EQ.0D0) THEN
         IERR = 1
         RETURN
      END IF

      PRJ(13) = 1D0/PRJ(12)

      PRJ(15) = PRJ(10)/PRJ(12)
      PRJ(16) = 1D0 + SIND(THETA1)*SIND(THETA2)
      PRJ(17) = 2D0*PRJ(12)
      PRJ(18) = PRJ(15)*PRJ(15)*PRJ(16)
      PRJ(19) = 1D0/(2D0*PRJ(10)*PRJ(15))
      PRJ(20) = PRJ(15)*SQRT(PRJ(16) + PRJ(17))

      PRJ(14) = PRJ(15)*SQRT(PRJ(16) - PRJ(17)*SIND(PRJ(1)))

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE COEFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION A, PHI, PRJ(0:20), R, THETA, X, Y

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.502D0) THEN
         CALL COESET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      A = PHI*PRJ(12)
      IF (THETA.EQ.-90D0) THEN
         R = PRJ(20)
      ELSE
         R = PRJ(15)*SQRT(PRJ(16) - PRJ(17)*SIND(THETA))
      END IF

      X =           R*SIND(A)
      Y = PRJ(14) - R*COSD(A)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE COEREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION A, DY, PHI, PRJ(0:20), R, THETA, TOL, W, X, Y

      PARAMETER (TOL = 1D-12)

      DOUBLE PRECISION ASIND, ATAN2D
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.502D0) THEN
         CALL COESET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      DY = PRJ(14) - Y
      R  = SQRT(X*X + DY*DY)
      IF (PRJ(1).LT.0D0) R = -R

      IF (R.EQ.0D0) THEN
         A = 0D0
      ELSE
         A = ATAN2D(X/R, DY/R)
      END IF

      PHI = A*PRJ(13)
      IF (ABS(R - PRJ(20)).LT.TOL) THEN
         THETA = -90D0
      ELSE
         W = (PRJ(18) - R*R)*PRJ(19)
         IF (ABS(W).GT.1D0) THEN
            IF (ABS(W-1D0).LT.TOL) THEN
               THETA = 90D0
            ELSE IF (ABS(W+1D0).LT.TOL) THEN
               THETA = -90D0
            ELSE
               IERR = 2
               RETURN
            END IF
         ELSE
            THETA = ASIND(W)
         END IF
      END IF

      IERR = 0
      RETURN
      END

*=======================================================================
*   COD: conic equidistant projection.
*
*   Given:
*      PRJ(1)    sigma = (theta2+theta1)/2
*      PRJ(2)    delta = (theta2-theta1)/2, where theta1 and theta2 are
*                the latitudes of the standard parallels, in degrees.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 503.
*      PRJ(12)   C = r0*sin(sigma)*sin(delta)/delta
*      PRJ(13)   1/C
*      PRJ(14)   Y0 = delta*cot(delta)*cot(sigma)
*      PRJ(15)   Y0 + sigma
*=======================================================================
      SUBROUTINE CODSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      PRJ(11) = 503D0

      IF (PRJ(10).EQ.0D0) PRJ(10) = R2D

      IF (PRJ(2).EQ.0D0) THEN
         PRJ(12) = PRJ(10)*SIND(PRJ(1))*D2R
      ELSE
         PRJ(12) = PRJ(10)*SIND(PRJ(1))*SIND(PRJ(2))/PRJ(2)
      END IF

      IF (PRJ(12).EQ.0D0) THEN
         IERR = 1
         RETURN
      END IF

      PRJ(13) = 1D0/PRJ(12)
      PRJ(14) = PRJ(10)*COSD(PRJ(2))*COSD(PRJ(1))/PRJ(12)
      PRJ(15) = PRJ(14) + PRJ(1)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE CODFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION A, PHI, PRJ(0:20), R, THETA, X, Y

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.503D0) THEN
         CALL CODSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      A = PRJ(12)*PHI
      R = PRJ(15) - THETA

      X =           R*SIND(A)
      Y = PRJ(14) - R*COSD(A)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE CODREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION A, DY, PHI, PRJ(0:20), R, THETA, X, Y

      DOUBLE PRECISION ATAN2D
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.503D0) THEN
         CALL CODSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      DY = PRJ(14) - Y
      R  = SQRT(X*X + DY*DY)
      IF (PRJ(1).LT.0D0) R = -R

      IF (R.EQ.0D0) THEN
         A = 0D0
      ELSE
         A = ATAN2D(X/R, DY/R)
      END IF

      PHI   = A*PRJ(13)
      THETA = PRJ(15) - R

      IERR = 0
      RETURN
      END

*=======================================================================
*   COO: conic orthomorphic projection.
*
*   Given:
*      PRJ(1)    sigma = (theta2+theta1)/2
*      PRJ(2)    delta = (theta2-theta1)/2, where theta1 and theta2 are
*                the latitudes of the standard parallels, in degrees.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 504.
*      PRJ(12)   C = ln(cos(theta2)/cos(theta1))/ln(tan(tau2)/tan(tau1))
*                where tau1 = (90 - theta1)/2
*                      tau2 = (90 - theta2)/2
*      PRJ(13)   1/C
*      PRJ(14)   Y0 = psi*tan((90-sigma)/2)**C
*      PRJ(15)   psi = (r0*cos(theta1)/C)/tan(tau1)**C
*      PRJ(16)   1/psi
*=======================================================================
      SUBROUTINE COOSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION COS1, COS2, PRJ(0:20), TAN1, TAN2, THETA1, THETA2

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      DOUBLE PRECISION COSD, SIND, TAND
*-----------------------------------------------------------------------
      PRJ(11) = 504D0

      IF (PRJ(10).EQ.0D0) PRJ(10) = R2D

      THETA1 = PRJ(1) - PRJ(2)
      THETA2 = PRJ(1) + PRJ(2)

      TAN1 = TAND((90D0 - THETA1)/2D0)
      COS1 = COSD(THETA1)

      IF (THETA1.EQ.THETA2) THEN
         PRJ(12) = SIND(THETA1)
      ELSE
         TAN2 = TAND((90D0 - THETA2)/2D0)
         COS2 = COSD(THETA2)
         PRJ(12) = LOG(COS2/COS1)/LOG(TAN2/TAN1)
      END IF
      IF (PRJ(12).EQ.0D0) THEN
         IERR = 1
         RETURN
      END IF

      PRJ(13) = 1D0/PRJ(12)

      PRJ(15) = PRJ(10)*(COS1/PRJ(12))/(TAN1)**PRJ(12)
      IF (PRJ(15).EQ.0D0) THEN
         IERR = 1
         RETURN
      END IF
      PRJ(16) = 1D0/PRJ(15)
      PRJ(14) = PRJ(15)*TAND((90D0-PRJ(1))/2D0)**PRJ(12)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE COOFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION A, PHI, PRJ(0:20), R, THETA, X, Y

      DOUBLE PRECISION COSD, SIND, TAND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.504D0) THEN
         CALL COOSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      A = PRJ(12)*PHI
      IF (THETA.EQ.-90D0) THEN
         IF (PRJ(12).LT.0D0) THEN
            R = 0D0
         ELSE
            IERR = 2
            RETURN
         END IF
      ELSE
         R = PRJ(15)*TAND((90D0 - THETA)/2D0)**PRJ(12)
      END IF

      X =           R*SIND(A)
      Y = PRJ(14) - R*COSD(A)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE COOREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION A, DY, PHI, PRJ(0:20), R, THETA, X, Y

      DOUBLE PRECISION ATAN2D, ATAND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.504D0) THEN
         CALL COOSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      DY = PRJ(14) - Y
      R  = SQRT(X*X + DY*DY)
      IF (PRJ(1).LT.0D0) R = -R

      IF (R.EQ.0D0) THEN
         A = 0D0
      ELSE
         A = ATAN2D(X/R, DY/R)
      END IF

      PHI = A*PRJ(13)
      IF (R.EQ.0D0) THEN
         IF (PRJ(12).LT.0D0) THEN
            THETA = -90D0
         ELSE
            RETURN
         END IF
      ELSE
         THETA = 90D0 - 2D0*ATAND((R*PRJ(16))**PRJ(13))
      END IF

      IERR = 0
      RETURN
      END

*=======================================================================
*   BON: Bonne's projection.
*
*   Given:
*      PRJ(1)    Bonne conformal latitude, theta1, in degrees.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 601.
*      PRJ(12)   r0*pi/180
*      PRJ(13)   Y0 = r0*(cot(theta1) + theta1*pi/180)
*=======================================================================
      SUBROUTINE BONSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      PRJ(11) = 601D0

      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D
         PRJ(12) = 1D0
         PRJ(13) = PRJ(10)*COSD(PRJ(1))/SIND(PRJ(1)) + PRJ(1)
      ELSE
         PRJ(12) = PRJ(10)*D2R
         PRJ(13) = PRJ(10)*(COSD(PRJ(1))/SIND(PRJ(1)) + PRJ(1)*D2R)
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE BONFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION A, PHI, PRJ(0:20), R, THETA, X, Y

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(1).EQ.0D0) THEN
*        Sanson-Flamsteed.
         CALL SFLFWD (PHI, THETA, PRJ, X, Y, IERR)
         RETURN
      END IF

      IF (PRJ(11).NE.601D0) THEN
         CALL BONSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      R = PRJ(13) - THETA*PRJ(12)
      A = PRJ(10)*PHI*COSD(THETA)/R

      X =           R*SIND(A)
      Y = PRJ(13) - R*COSD(A)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE BONREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION A, CTHE, DY, PHI, PRJ(0:20), R, THETA, X, Y

      DOUBLE PRECISION ATAN2D, COSD
*-----------------------------------------------------------------------
      IF (PRJ(1).EQ.0D0) THEN
*        Sanson-Flamsteed.
         CALL SFLREV (X, Y, PRJ, PHI, THETA, IERR)
         RETURN
      END IF

      IF (PRJ(11).NE.601D0) THEN
         CALL BONSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      DY = PRJ(13) - Y
      R = SQRT(X*X + DY*DY)
      IF (PRJ(1).LT.0D0) R = -R

      IF (R.EQ.0D0) THEN
         A = 0D0
      ELSE
         A = ATAN2D(X/R, DY/R)
      END IF

      THETA = (PRJ(13) - R)/PRJ(12)
      CTHE = COSD(THETA)
      IF (CTHE.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = A*(R/PRJ(10))/CTHE
      END IF

      IERR = 0
      RETURN
      END

*=======================================================================
*   PCO: polyconic projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 602.
*      PRJ(12)   r0*(pi/180)
*      PRJ(13)   1/r0
*      PRJ(14)   2*r0
*=======================================================================
      SUBROUTINE PCOSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      PRJ(11) = 602D0

      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D
         PRJ(12) = 1D0
         PRJ(13) = 1D0
         PRJ(14) = 360D0/PI
      ELSE
         PRJ(12) = PRJ(10)*D2R
         PRJ(13) = 1D0/PRJ(12)
         PRJ(14) = 2D0*PRJ(10)
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE PCOFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION A, COTTHE, CTHE, PHI, PRJ(0:20), STHE, THETA, X,
     *          Y

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.602D0) THEN
         CALL PCOSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      CTHE = COSD(THETA)
      STHE = SIND(THETA)
      A = PHI*STHE

      IF (STHE.EQ.0D0) THEN
         X = PRJ(12)*PHI
         Y = 0D0
      ELSE
         COTTHE = CTHE/STHE
         X = PRJ(10)*COTTHE*SIND(A)
         Y = PRJ(10)*(COTTHE*(1D0 - COSD(A)) + THETA*D2R)
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE PCOREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR, J
      DOUBLE PRECISION F, FNEG, FPOS, LAMBDA, PHI, PRJ(0:20), TANTHE,
     *          THENEG, THEPOS, THETA, TOL, W, X, XP, XX, Y, YMTHE, YP

      PARAMETER (TOL = 1D-12)

      DOUBLE PRECISION ATAN2D, SIND, TAND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.602D0) THEN
         CALL PCOSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      W = ABS(Y*PRJ(13))
      IF (W.LT.TOL) THEN
         PHI = X*PRJ(13)
         THETA = 0D0
      ELSE IF (ABS(W-90D0).LT.TOL) THEN
         PHI = 0D0
         THETA = SIGN(90D0,Y)
      ELSE
*        Iterative solution using weighted division of the interval.
         IF (Y.GT.0D0) THEN
            THEPOS =  90D0
         ELSE
            THEPOS = -90D0
         END IF
         THENEG = 0D0

         XX = X*X
         YMTHE = Y - PRJ(12)*THEPOS
         FPOS = XX + YMTHE*YMTHE
         FNEG = -999D0

         DO 10 J = 1, 64
            IF (FNEG.LT.-100D0) THEN
*              Equal division of the interval.
               THETA = (THEPOS+THENEG)/2D0
            ELSE
*              Weighted division of the interval.
               LAMBDA = FPOS/(FPOS-FNEG)
               IF (LAMBDA.LT.0.1D0) THEN
                  LAMBDA = 0.1D0
               ELSE IF (LAMBDA.GT.0.9D0) THEN
                  LAMBDA = 0.9D0
               END IF
               THETA = THEPOS - LAMBDA*(THEPOS-THENEG)
            END IF

*           Compute the residue.
            YMTHE = Y - PRJ(12)*THETA
            TANTHE = TAND(THETA)
            F = XX + YMTHE*(YMTHE - PRJ(14)/TANTHE)

*           Check for convergence.
            IF (ABS(F).LT.TOL) GO TO 20
            IF (ABS(THEPOS-THENEG).LT.TOL) GO TO 20

*           Redefine the interval.
            IF (F.GT.0D0) THEN
               THEPOS = THETA
               FPOS = F
            ELSE
               THENEG = THETA
               FNEG = F
            END IF
 10      CONTINUE

 20      XP = PRJ(10) - YMTHE*TANTHE
         YP = X*TANTHE
         IF (XP.EQ.0D0 .AND. YP.EQ.0D0) THEN
            PHI = 0D0
         ELSE
            PHI = ATAN2D(YP, XP)/SIND(THETA)
         END IF
      END IF

      IERR = 0
      RETURN
      END

*=======================================================================
*   TSC: tangential spherical cube projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 701.
*      PRJ(12)   r0*(pi/4)
*      PRJ(13)   (4/pi)/r0
*=======================================================================
      SUBROUTINE TSCSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      PRJ(11) = 701D0

      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D
         PRJ(12) = 45D0
         PRJ(13) = 1D0/45D0
      ELSE
         PRJ(12) = PRJ(10)*PI/4D0
         PRJ(13) = 1D0/PRJ(12)
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE TSCFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   FACE, IERR
      DOUBLE PRECISION CTHE, L, M, N, PHI, PRJ(0:20), RHO, THETA, TOL,
     *          X, X0, XF, Y, Y0, YF

      PARAMETER (TOL = 1D-12)

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.701D0) THEN
         CALL TSCSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      CTHE = COSD(THETA)
      L = CTHE*COSD(PHI)
      M = CTHE*SIND(PHI)
      N = SIND(THETA)

      FACE = 0
      RHO  = N
      IF (L.GT.RHO) THEN
         FACE = 1
         RHO  = L
      END IF
      IF (M.GT.RHO) THEN
         FACE = 2
         RHO  = M
      END IF
      IF (-L.GT.RHO) THEN
         FACE = 3
         RHO  = -L
      END IF
      IF (-M.GT.RHO) THEN
         FACE = 4
         RHO  = -M
      END IF
      IF (-N.GT.RHO) THEN
         FACE = 5
         RHO  = -N
      END IF

      IF (FACE.EQ.0) THEN
         XF =  M/RHO
         YF = -L/RHO
         X0 =  0D0
         Y0 =  2D0
      ELSE IF (FACE.EQ.1) THEN
         XF =  M/RHO
         YF =  N/RHO
         X0 =  0D0
         Y0 =  0D0
      ELSE IF (FACE.EQ.2) THEN
         XF = -L/RHO
         YF =  N/RHO
         X0 =  2D0
         Y0 =  0D0
      ELSE IF (FACE.EQ.3) THEN
         XF = -M/RHO
         YF =  N/RHO
         X0 =  4D0
         Y0 =  0D0
      ELSE IF (FACE.EQ.4) THEN
         XF =  L/RHO
         YF =  N/RHO
         X0 =  6D0
         Y0 =  0D0
      ELSE IF (FACE.EQ.5) THEN
         XF =  M/RHO
         YF =  L/RHO
         X0 =  0D0
         Y0 = -2D0
      END IF

      IF (ABS(XF).GT.1D0) THEN
         IF (ABS(XF).GT.1D0+TOL) THEN
            IERR = 2
            RETURN
         END IF
         XF = SIGN(1D0,XF)
      END IF
      IF (ABS(YF).GT.1D0) THEN
         IF (ABS(YF).GT.1D0+TOL) THEN
            IERR = 2
            RETURN
         END IF
         YF = SIGN(1D0,YF)
      END IF

      X = PRJ(12)*(XF + X0)
      Y = PRJ(12)*(YF + Y0)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE TSCREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION L, M, N, PHI, PRJ(0:20), THETA, X, XF, Y, YF

      DOUBLE PRECISION ASIND, ATAN2D
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.701D0) THEN
         CALL TSCSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      XF = X*PRJ(13)
      YF = Y*PRJ(13)

*     Check bounds.
      IF (ABS(XF).LE.1D0) THEN
         IF (ABS(YF).GT.3D0) THEN
            IERR = 2
            RETURN
         END IF
      ELSE
         IF (ABS(XF).GT.7D0 .OR. ABS(YF).GT.1D0) THEN
            IERR = 2
            RETURN
         END IF
      END IF

*     Map negative faces to the other side.
      IF (XF.LT.-1D0) XF = XF + 8D0

*     Determine the face.
      IF (XF.GT.5D0) THEN
*        FACE = 4
         XF = XF - 6D0
         M  = -1D0/SQRT(1D0 + XF*XF + YF*YF)
         L  = -M*XF
         N  = -M*YF
      ELSE IF (XF.GT.3D0) THEN
*        FACE = 3
         XF = XF - 4D0
         L  = -1D0/SQRT(1D0 + XF*XF + YF*YF)
         M  =  L*XF
         N  = -L*YF
      ELSE IF (XF.GT.1D0) THEN
*        FACE = 2
         XF = XF - 2D0
         M  =  1D0/SQRT(1D0 + XF*XF + YF*YF)
         L  = -M*XF
         N  =  M*YF
      ELSE IF (YF.GT.1D0) THEN
*        FACE = 0
         YF = YF - 2D0
         N  = 1D0/SQRT(1D0 + XF*XF + YF*YF)
         L  = -N*YF
         M  =  N*XF
      ELSE IF (YF.LT.-1D0) THEN
*        FACE = 5
         YF = YF + 2D0
         N  = -1D0/SQRT(1D0 + XF*XF + YF*YF)
         L  = -N*YF
         M  = -N*XF
      ELSE
*        FACE = 1
         L  =  1D0/SQRT(1D0 + XF*XF + YF*YF)
         M  =  L*XF
         N  =  L*YF
      END IF

      IF (L.EQ.0D0 .AND. M.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = ATAN2D(M, L)
      END IF
      THETA = ASIND(N)

      IERR = 0
      RETURN
      END

*=======================================================================
*   CSC: COBE quadrilateralized spherical cube projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 702.
*      PRJ(12)   r0*(pi/4)
*      PRJ(13)   (4/pi)/r0
*=======================================================================
      SUBROUTINE CSCSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      PRJ(11) = 702D0

      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D
         PRJ(12) = 45D0
         PRJ(13) = 1D0/45D0
      ELSE
         PRJ(12) = PRJ(10)*PI/4D0
         PRJ(13) = 1D0/PRJ(12)
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE CSCFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   FACE, IERR
      DOUBLE PRECISION CTHE, ETA, L, M, N, PHI, PRJ(0:20), RHO, THETA,
     *          X, XI, Y

      REAL      A, A2, A2B2, A4, AB, B, B2, B4, CA2, CB2, TOL, X0, XF,
     *          Y0, YF
      REAL      C00, C10, C01, C11, C20, C02, D0, D1, MM, GAMMA, GSTAR,
     *          OMEGA1
      PARAMETER (GSTAR  =  1.37484847732)
      PARAMETER (MM     =  0.004869491981)
      PARAMETER (GAMMA  = -0.13161671474)
      PARAMETER (OMEGA1 = -0.159596235474)
      PARAMETER (D0  =  0.0759196200467)
      PARAMETER (D1  = -0.0217762490699)
      PARAMETER (C00 =  0.141189631152)
      PARAMETER (C10 =  0.0809701286525)
      PARAMETER (C01 = -0.281528535557)
      PARAMETER (C11 =  0.15384112876)
      PARAMETER (C20 = -0.178251207466)
      PARAMETER (C02 =  0.106959469314)

      PARAMETER (TOL = 1E-7)

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.702D0) THEN
         CALL CSCSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      CTHE = COSD(THETA)
      L = CTHE*COSD(PHI)
      M = CTHE*SIND(PHI)
      N = SIND(THETA)

      FACE = 0
      RHO  = N
      IF (L.GT.RHO) THEN
         FACE = 1
         RHO  = L
      END IF
      IF (M.GT.RHO) THEN
         FACE = 2
         RHO  = M
      END IF
      IF (-L.GT.RHO) THEN
         FACE = 3
         RHO  = -L
      END IF
      IF (-M.GT.RHO) THEN
         FACE = 4
         RHO  = -M
      END IF
      IF (-N.GT.RHO) THEN
         FACE = 5
         RHO  = -N
      END IF

      IF (FACE.EQ.0) THEN
         XI  =  M
         ETA = -L
         X0  =  0.0
         Y0  =  2.0
      ELSE IF (FACE.EQ.1) THEN
         XI  =  M
         ETA =  N
         X0  =  0.0
         Y0  =  0.0
      ELSE IF (FACE.EQ.2) THEN
         XI  = -L
         ETA =  N
         X0  =  2.0
         Y0  =  0.0
      ELSE IF (FACE.EQ.3) THEN
         XI  = -M
         ETA =  N
         X0  =  4.0
         Y0  =  0.0
      ELSE IF (FACE.EQ.4) THEN
         XI  =  L
         ETA =  N
         X0  =  6.0
         Y0  =  0.0
      ELSE IF (FACE.EQ.5) THEN
         XI  =  M
         ETA =  L
         X0  =  0.0
         Y0  = -2.0
      END IF

      A =  XI/RHO
      B = ETA/RHO

      A2 = A*A
      B2 = B*B
      CA2 = 1.0 - A2
      CB2 = 1.0 - B2

*     Avoid floating underflows.
      A4 = 0.0
      B4 = 0.0
      AB = ABS(A*B)
      A2B2 = 0.0
      IF (A2.GT.1E-16) A4 = A2*A2
      IF (B2.GT.1E-16) B4 = B2*B2
      IF (AB.GT.1E-16) A2B2 = AB*AB

      XF = A*(A2 + CA2*(GSTAR + B2*(GAMMA*CA2 + MM*A2 +
     *       CB2*(C00 + C10*A2 + C01*B2 + C11*A2B2 + C20*A4 + C02*B4)) +
     *       A2*(OMEGA1 - CA2*(D0 + D1*A2))))
      YF = B*(B2 + CB2*(GSTAR + A2*(GAMMA*CB2 + MM*B2 +
     *       CA2*(C00 + C10*B2 + C01*A2 + C11*A2B2 + C20*B4 + C02*A4)) +
     *       B2*(OMEGA1 - CB2*(D0 + D1*B2))))

      IF (ABS(XF).GT.1.0) THEN
         IF (ABS(XF).GT.1.0+TOL) THEN
            IERR = 2
            RETURN
         END IF
         XF = SIGN(1.0,XF)
      END IF
      IF (ABS(YF).GT.1.0) THEN
         IF (ABS(YF).GT.1.0+TOL) THEN
            IERR = 2
            RETURN
         END IF
         YF = SIGN(1.0,YF)
      END IF

      X = PRJ(12)*(X0 + XF)
      Y = PRJ(12)*(Y0 + YF)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE CSCREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   FACE, IERR
      DOUBLE PRECISION L, M, N, PHI, PRJ(0:20), THETA, X, Y

      REAL      A, B, XF, XX, YF, YY
      REAL      P00, P01, P02, P03, P04, P05, P06, P10, P11, P12,
     *          P13, P14, P15, P20, P21, P22, P23, P24, P30, P31, P32,
     *          P33, P40, P41, P42, P50, P51, P60
      PARAMETER (P00 = -0.27292696)
      PARAMETER (P10 = -0.07629969)
      PARAMETER (P20 = -0.22797056)
      PARAMETER (P30 =  0.54852384)
      PARAMETER (P40 = -0.62930065)
      PARAMETER (P50 =  0.25795794)
      PARAMETER (P60 =  0.02584375)
      PARAMETER (P01 = -0.02819452)
      PARAMETER (P11 = -0.01471565)
      PARAMETER (P21 =  0.48051509)
      PARAMETER (P31 = -1.74114454)
      PARAMETER (P41 =  1.71547508)
      PARAMETER (P51 = -0.53022337)
      PARAMETER (P02 =  0.27058160)
      PARAMETER (P12 = -0.56800938)
      PARAMETER (P22 =  0.30803317)
      PARAMETER (P32 =  0.98938102)
      PARAMETER (P42 = -0.83180469)
      PARAMETER (P03 = -0.60441560)
      PARAMETER (P13 =  1.50880086)
      PARAMETER (P23 = -0.93678576)
      PARAMETER (P33 =  0.08693841)
      PARAMETER (P04 =  0.93412077)
      PARAMETER (P14 = -1.41601920)
      PARAMETER (P24 =  0.33887446)
      PARAMETER (P05 = -0.63915306)
      PARAMETER (P15 =  0.52032238)
      PARAMETER (P06 =  0.14381585)

      DOUBLE PRECISION ASIND, ATAN2D
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.702D0) THEN
         CALL CSCSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      XF = X*PRJ(13)
      YF = Y*PRJ(13)

*     Check bounds.
      IF (ABS(XF).LE.1D0) THEN
         IF (ABS(YF).GT.3D0) THEN
            IERR = 2
            RETURN
         END IF
      ELSE
         IF (ABS(XF).GT.7D0 .OR. ABS(YF).GT.1D0) THEN
            IERR = 2
            RETURN
         END IF
      END IF

*     Map negative faces to the other side.
      IF (XF.LT.-1D0) XF = XF + 8D0

*     Determine the face.
      IF (XF.GT.5.0) THEN
         FACE = 4
         XF = XF - 6.0
      ELSE IF (XF.GT.3.0) THEN
         FACE = 3
         XF = XF - 4.0
      ELSE IF (XF.GT.1.0) THEN
         FACE = 2
         XF = XF - 2.0
      ELSE IF (YF.GT.1.0) THEN
         FACE = 0
         YF = YF - 2.0
      ELSE IF (YF.LT.-1.0) THEN
         FACE = 5
         YF = YF + 2.0
      ELSE
         FACE = 1
      END IF

      XX  =  XF*XF
      YY  =  YF*YF

      A =   (P00+XX*(P10+XX*(P20+XX*(P30+XX*(P40+XX*(P50+XX*(P60)))))) +
     *   YY*(P01+XX*(P11+XX*(P21+XX*(P31+XX*(P41+XX*(P51))))) +
     *   YY*(P02+XX*(P12+XX*(P22+XX*(P32+XX*(P42)))) +
     *   YY*(P03+XX*(P13+XX*(P23+XX*(P33))) +
     *   YY*(P04+XX*(P14+XX*(P24)) +
     *   YY*(P05+XX*(P15) +
     *   YY*(P06)))))))
      A = XF + XF*(1.0 - XX)*A

      B =   (P00+YY*(P10+YY*(P20+YY*(P30+YY*(P40+YY*(P50+YY*(P60)))))) +
     *   XX*(P01+YY*(P11+YY*(P21+YY*(P31+YY*(P41+YY*(P51))))) +
     *   XX*(P02+YY*(P12+YY*(P22+YY*(P32+YY*(P42)))) +
     *   XX*(P03+YY*(P13+YY*(P23+YY*(P33))) +
     *   XX*(P04+YY*(P14+YY*(P24)) +
     *   XX*(P05+YY*(P15) +
     *   XX*(P06)))))))
      B = YF + YF*(1.0 - YY)*B

      IF (FACE.EQ.0) THEN
         N =  1D0/SQRT(A*A + B*B + 1D0)
         L = -B*N
         M =  A*N
      ELSE IF (FACE.EQ.1) THEN
         L =  1D0/SQRT(A*A + B*B + 1D0)
         M =  A*L
         N =  B*L
      ELSE IF (FACE.EQ.2) THEN
         M =  1D0/SQRT(A*A + B*B + 1D0)
         L = -A*M
         N =  B*M
      ELSE IF (FACE.EQ.3) THEN
         L = -1D0/SQRT(A*A + B*B + 1D0)
         M =  A*L
         N = -B*L
      ELSE IF (FACE.EQ.4) THEN
         M = -1D0/SQRT(A*A + B*B + 1D0)
         L = -A*M
         N = -B*M
      ELSE IF (FACE.EQ.5) THEN
         N = -1D0/SQRT(A*A + B*B + 1D0)
         L = -B*N
         M = -A*N
      END IF

      IF (L.EQ.0D0 .AND. M.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = ATAN2D(M, L)
      END IF
      THETA = ASIND(N)

      IERR = 0
      RETURN
      END

*=======================================================================
*   QSC: quadrilaterilized spherical cube projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*
*   Returned:
*      PRJ(11)   State flag, 703.
*      PRJ(12)   r0*(pi/4)
*      PRJ(13)   (4/pi)/r0
*=======================================================================
      SUBROUTINE QSCSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      PRJ(11) = 703D0

      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D
         PRJ(12) = 45D0
         PRJ(13) = 1D0/45D0
      ELSE
         PRJ(12) = PRJ(10)*PI/4D0
         PRJ(13) = 1D0/PRJ(12)
      END IF

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE QSCFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   FACE, IERR
      DOUBLE PRECISION CTHE, ETA, L, M, N, OMEGA, P, PHI, PRJ(0:20),
     *          RHO, RHU, T, TAU, THETA, TOL, X, X0, XF, XI, Y, Y0, YF

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      PARAMETER (TOL = 1D-12)

      DOUBLE PRECISION ASIND, ATAND, COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.703D0) THEN
         CALL QSCSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      IF (ABS(THETA).EQ.90D0) THEN
         X = 0D0
         Y = SIGN(2D0*PRJ(12),THETA)
         IERR = 0
         RETURN
      END IF

      CTHE = COSD(THETA)
      L = CTHE*COSD(PHI)
      M = CTHE*SIND(PHI)
      N = SIND(THETA)

      FACE = 0
      RHO  = N
      IF (L.GT.RHO) THEN
         FACE = 1
         RHO  = L
      END IF
      IF (M.GT.RHO) THEN
         FACE = 2
         RHO  = M
      END IF
      IF (-L.GT.RHO) THEN
         FACE = 3
         RHO  = -L
      END IF
      IF (-M.GT.RHO) THEN
         FACE = 4
         RHO  = -M
      END IF
      IF (-N.GT.RHO) THEN
         FACE = 5
         RHO  = -N
      END IF

      RHU = 1D0 - RHO

      IF (FACE.EQ.0) THEN
         XI  =  M
         ETA = -L
         IF (RHU.LT.1D-8) THEN
*           Small angle formula.
            T = (90D0 - THETA)*D2R
            RHU = T*T/2D0
         END IF
         X0  =  0D0
         Y0  =  2D0
      ELSE IF (FACE.EQ.1) THEN
         XI  =  M
         ETA =  N
         IF (RHU.LT.1D-8) THEN
*           Small angle formula.
            T = THETA*D2R
            P = MOD(PHI,360D0)
            IF (P.LT.-180D0) P = P + 360D0
            IF (P.GT.+180D0) P = P - 360D0
            P = P*D2R
            RHU = (P*P + T*T)/2D0
         END IF
         X0  =  0D0
         Y0  =  0D0
      ELSE IF (FACE.EQ.2) THEN
         XI  = -L
         ETA =  N
         IF (RHU.LT.1D-8) THEN
*           Small angle formula.
            T = THETA*D2R
            P = MOD(PHI,360D0)
            IF (P.LT.-180D0) P = P + 360D0
            P = (90D0 - P)*D2R
            RHU = (P*P + T*T)/2D0
         END IF
         X0  =  2D0
         Y0  =  0D0
      ELSE IF (FACE.EQ.3) THEN
         XI  = -M
         ETA =  N
         IF (RHU.LT.1D-8) THEN
*           Small angle formula.
            T = THETA*D2R
            P = MOD(PHI,360D0)
            IF (P.LT.0D0) P = P + 360D0
            P = (180D0 - P)*D2R
            RHU = (P*P + T*T)/2D0
         END IF
         X0  =  4D0
         Y0  =  0D0
      ELSE IF (FACE.EQ.4) THEN
         XI  =  L
         ETA =  N
         IF (RHU.LT.1D-8) THEN
*           Small angle formula.
            T = THETA*D2R
            P = MOD(PHI,360D0)
            IF (P.GT.+180D0) P = P - 360D0
            P = (90D0 + P)*D2R
            RHU = (P*P + T*T)/2D0
         END IF
         X0  =  6
         Y0  =  0D0
      ELSE IF (FACE.EQ.5) THEN
         XI  =  M
         ETA =  L
         IF (RHU.LT.1D-8) THEN
*           Small angle formula.
            T = (90D0 + THETA)*D2R
            RHU = T*T/2D0
         END IF
         X0  =  0D0
         Y0  = -2
      END IF

      IF (XI.EQ.0D0 .AND. ETA.EQ.0D0) THEN
         XF  = 0D0
         YF  = 0D0
      ELSE IF (-XI.GE.ABS(ETA)) THEN
         OMEGA = ETA/XI
         TAU = 1D0 + OMEGA*OMEGA
         XF  = -SQRT(RHU/(1D0-1D0/SQRT(1D0+TAU)))
         YF  = (XF/15D0)*(ATAND(OMEGA) - ASIND(OMEGA/SQRT(TAU+TAU)))
      ELSE IF (XI.GE.ABS(ETA)) THEN
         OMEGA = ETA/XI
         TAU = 1D0 + OMEGA*OMEGA
         XF  =  SQRT(RHU/(1D0-1D0/SQRT(1D0+TAU)))
         YF  = (XF/15D0)*(ATAND(OMEGA) - ASIND(OMEGA/SQRT(TAU+TAU)))
      ELSE IF (-ETA.GT.ABS(XI)) THEN
         OMEGA = XI/ETA
         TAU = 1D0 + OMEGA*OMEGA
         YF  = -SQRT(RHU/(1D0-1D0/SQRT(1D0+TAU)))
         XF  = (YF/15D0)*(ATAND(OMEGA) - ASIND(OMEGA/SQRT(TAU+TAU)))
      ELSE IF (ETA.GT.ABS(XI)) THEN
         OMEGA = XI/ETA
         TAU = 1D0 + OMEGA*OMEGA
         YF  =  SQRT(RHU/(1D0-1D0/SQRT(1D0+TAU)))
         XF  = (YF/15D0)*(ATAND(OMEGA) - ASIND(OMEGA/SQRT(TAU+TAU)))
      END IF

      IF (ABS(XF).GT.1D0) THEN
         IF (ABS(XF).GT.1D0+TOL) THEN
            IERR = 2
            RETURN
         END IF
         XF = SIGN(1D0,XF)
      END IF
      IF (ABS(YF).GT.1D0) THEN
         IF (ABS(YF).GT.1D0+TOL) THEN
            IERR = 2
            RETURN
         END IF
         YF = SIGN(1D0,YF)
      END IF

      X = PRJ(12)*(XF + X0)
      Y = PRJ(12)*(YF + Y0)


      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE QSCREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      LOGICAL   DIRECT
      INTEGER   FACE, IERR
      DOUBLE PRECISION L, M, N, OMEGA, PHI, PRJ(0:20), RHO, RHU, TAU,
     *          THETA, TOL, X, XF, Y, YF, W

      PARAMETER (TOL = 1D-12)

      DOUBLE PRECISION C2
      PARAMETER (C2 = 1D0/1.4142135623730950488D0)

      DOUBLE PRECISION ASIND, ATAN2D, COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.703D0) THEN
         CALL QSCSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      XF = X*PRJ(13)
      YF = Y*PRJ(13)

*     Check bounds.
      IF (ABS(XF).LE.1D0) THEN
         IF (ABS(YF).GT.3D0) THEN
            IERR = 2
            RETURN
         END IF
      ELSE
         IF (ABS(XF).GT.7D0 .OR. ABS(YF).GT.1D0) THEN
            IERR = 2
            RETURN
         END IF
      END IF

*     Map negative faces to the other side.
      IF (XF.LT.-1D0) XF = XF + 8D0

*     Determine the face.
      IF (XF.GT.5D0) THEN
         FACE = 4
         XF = XF - 6D0
      ELSE IF (XF.GT.3D0) THEN
         FACE = 3
         XF = XF - 4D0
      ELSE IF (XF.GT.1D0) THEN
         FACE = 2
         XF = XF - 2D0
      ELSE IF (YF.GT.1D0) THEN
         FACE = 0
         YF = YF - 2D0
      ELSE IF (YF.LT.-1D0) THEN
         FACE = 5
         YF = YF + 2D0
      ELSE
         FACE = 1
      END IF

      DIRECT = ABS(XF).GT.ABS(YF)
      IF (DIRECT) THEN
         IF (XF.EQ.0D0) THEN
            OMEGA = 0D0
            TAU = 1D0
            RHO = 1D0
            RHU = 0D0
         ELSE
            W = 15D0*YF/XF
            OMEGA = SIND(W)/(COSD(W) - C2)
            TAU = 1D0 + OMEGA*OMEGA
            RHU = XF*XF*(1D0 - 1D0/SQRT(1D0 + TAU))
            RHO = 1D0 - RHU
         END IF
      ELSE
         IF (YF.EQ.0D0) THEN
            OMEGA = 0D0
            TAU = 1D0
            RHO = 1D0
            RHU = 0D0
         ELSE
            W = 15D0*XF/YF
            OMEGA = SIND(W)/(COSD(W) - C2)
            TAU = 1D0 + OMEGA*OMEGA
            RHU = YF*YF*(1D0 - 1D0/SQRT(1D0 + TAU))
            RHO = 1D0 - RHU
         END IF
      END IF

      IF (RHO.LT.-1D0) THEN
         IF (RHO.LT.-1D0-TOL) THEN
            IERR = 2
            RETURN
         END IF

         RHO = -1D0
         RHU =  2D0
         W   =  0D0
      ELSE
         W = SQRT(RHU*(2D0-RHU)/TAU)
      END IF

      IF (FACE.EQ.0) THEN
         N = RHO
         IF (DIRECT) THEN
            M = W
            IF (XF.LT.0D0) M = -M
            L = -M*OMEGA
         ELSE
            L = W
            IF (YF.GT.0D0) L = -L
            M = -L*OMEGA
         END IF
      ELSE IF (FACE.EQ.1) THEN
         L = RHO
         IF (DIRECT) THEN
            M = W
            IF (XF.LT.0D0) M = -M
            N = M*OMEGA
         ELSE
            N = W
            IF (YF.LT.0D0) N = -N
            M = N*OMEGA
         END IF
      ELSE IF (FACE.EQ.2) THEN
         M =  RHO
         IF (DIRECT) THEN
            L = W
            IF (XF.GT.0D0) L = -L
            N = -L*OMEGA
         ELSE
            N = W
            IF (YF.LT.0D0) N = -N
            L = -N*OMEGA
         END IF
      ELSE IF (FACE.EQ.3) THEN
         L = -RHO
         IF (DIRECT) THEN
            M = W
            IF (XF.GT.0D0) M = -M
            N = -M*OMEGA
         ELSE
            N = W
            IF (YF.LT.0D0) N = -N
            M = -N*OMEGA
         END IF
      ELSE IF (FACE.EQ.4) THEN
         M = -RHO
         IF (DIRECT) THEN
            L = W
            IF (XF.LT.0D0) L = -L
            N = L*OMEGA
         ELSE
            N = W
            IF (YF.LT.0D0) N = -N
            L = N*OMEGA
         END IF
      ELSE IF (FACE.EQ.5) THEN
         N = -RHO
         IF (DIRECT) THEN
            M = W
            IF (XF.LT.0D0) M = -M
            L = M*OMEGA
         ELSE
            L = W
            IF (YF.LT.0D0) L = -L
            M = L*OMEGA
         END IF
      END IF

      IF (L.EQ.0D0 .AND. M.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = ATAN2D(M, L)
      END IF
      THETA = ASIND(N)

      IERR = 0
      RETURN
      END
@


2.16
log
@Generalized AZP to support the tilt parameter, gamma; added support for the
SZP projection with SZPSET, SZPFWD and SZPREV; changed the call sequence of
PRJSET to return PHI0 along with THETA0; revised usage of the SIN parameters
to accord with changes to the paper and likewise renamed some variables.
@
text
@d4 1
a4 1
*   Copyright (C) 1995-2001, Mark Calabretta
d221 1
a221 1
*   $Id: proj.f,v 2.15 2000/12/06 23:53:57 mcalabre Exp mcalabre $
@


2.15
log
@Provide driver routines PRJSET, PRJFWD & PRJREV which key on the projection-
specific value of PRJ(11); renamed glsset(), glsfwd() & glsrev() to sflset(),
sflfwd() & sflrev(); reordered projections to match paper.
@
text
@d4 1
a4 1
*   Copyright (C) 1995-2000, Mark Calabretta
d45 1
d91 2
a92 2
*      THETA0   D        Native latitude of the reference point, in
*                        degrees.
d170 2
a171 2
*         disable strict bounds checking for the AZP, TAN, SIN, ZPN, and
*         COP projections.  The initialization routine for each
d198 1
a198 1
*      forward routines for AZP, TAN, SIN, ZPN, and COP also return
d221 1
a221 1
*   $Id: proj.f,v 2.14 2000/12/06 00:50:47 mcalabre Exp mcalabre $
d223 1
a223 1
      SUBROUTINE PRJSET (PCODE, PRJ, THETA0, IERR)
d225 8
a232 8
      INTEGER   IERR, IREF(25), NDX
      DOUBLE PRECISION PRJ(0:20), THETA0
      CHARACTER PCODES*99, PCODE*3

      DATA PCODES(01:40) /'AZP TAN STG SIN ARC ZPN ZEA AIR CYP CEA '/
      DATA PCODES(41:80) /'CAR MER SFL PAR MOL AIT COP COE COD COO '/
      DATA PCODES(81:99) /'BON PCO TSC CSC QSC'/
      DATA IREF /2,2,2,2,2,2,2,2,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0/
d242 2
a243 1
*     Native latitude of the reference point.
d256 1
a256 1
      GO TO (101, 102, 103, 104, 105, 106, 107, 108,
d268 4
d273 1
a273 1
 102  CALL TANSET (PRJ, IERR)
d277 1
a277 1
 103  CALL STGSET (PRJ, IERR)
d281 1
a281 1
 104  CALL SINSET (PRJ, IERR)
d285 1
a285 1
 105  CALL ARCSET (PRJ, IERR)
d289 1
a289 1
 106  CALL ZPNSET (PRJ, IERR)
d293 1
a293 1
 107  CALL ZEASET (PRJ, IERR)
d297 1
a297 1
 108  CALL AIRSET (PRJ, IERR)
d373 1
a373 1
      INTEGER   I, IERR, J, K, NDX(7,8)
d376 7
a382 7
      DATA (NDX(1,J),J=1,8) / 1,  2,  3,  4,  5,  6,  7,  8/
      DATA (NDX(2,J),J=1,8) / 9, 10, 11, 12,  0,  0,  0,  0/
      DATA (NDX(3,J),J=1,8) /13, 14, 15,  0,  0,  0,  0,  0/
      DATA (NDX(4,J),J=1,8) /16,  0,  0,  0,  0,  0,  0,  0/
      DATA (NDX(5,J),J=1,8) /17, 18, 19, 20,  0,  0,  0,  0/
      DATA (NDX(6,J),J=1,8) /21, 22,  0,  0,  0,  0,  0,  0/
      DATA (NDX(7,J),J=1,8) /23, 24, 25,  0,  0,  0,  0,  0/
d388 1
a388 1
      IF (I.LT.1 .OR. I.GT.7 .OR. J.LT.1 .OR. J.GT.8) THEN
d395 1
a395 1
      GO TO (101, 102, 103, 104, 105, 106, 107, 108,
d410 4
d415 1
a415 1
 102  CALL TANFWD (PHI, THETA, PRJ, X, Y, IERR)
d419 1
a419 1
 103  CALL STGFWD (PHI, THETA, PRJ, X, Y, IERR)
d423 1
a423 1
 104  CALL SINFWD (PHI, THETA, PRJ, X, Y, IERR)
d427 1
a427 1
 105  CALL ARCFWD (PHI, THETA, PRJ, X, Y, IERR)
d431 1
a431 1
 106  CALL ZPNFWD (PHI, THETA, PRJ, X, Y, IERR)
d435 1
a435 1
 107  CALL ZEAFWD (PHI, THETA, PRJ, X, Y, IERR)
d439 1
a439 1
 108  CALL AIRFWD (PHI, THETA, PRJ, X, Y, IERR)
d515 1
a515 1
      INTEGER   I, IERR, J, K, NDX(7,8)
d518 7
a524 7
      DATA (NDX(1,J),J=1,8) / 1,  2,  3,  4,  5,  6,  7,  8/
      DATA (NDX(2,J),J=1,8) / 9, 10, 11, 12,  0,  0,  0,  0/
      DATA (NDX(3,J),J=1,8) /13, 14, 15,  0,  0,  0,  0,  0/
      DATA (NDX(4,J),J=1,8) /16,  0,  0,  0,  0,  0,  0,  0/
      DATA (NDX(5,J),J=1,8) /17, 18, 19, 20,  0,  0,  0,  0/
      DATA (NDX(6,J),J=1,8) /21, 22,  0,  0,  0,  0,  0,  0/
      DATA (NDX(7,J),J=1,8) /23, 24, 25,  0,  0,  0,  0,  0/
d530 1
a530 1
      IF (I.LT.1 .OR. I.GT.7 .OR. J.LT.1 .OR. J.GT.8) THEN
d537 1
a537 1
      GO TO (101, 102, 103, 104, 105, 106, 107, 108,
d552 4
d557 1
a557 1
 102  CALL TANREV (X, Y, PRJ, PHI, THETA, IERR)
d561 1
a561 1
 103  CALL STGREV (X, Y, PRJ, PHI, THETA, IERR)
d565 1
a565 1
 104  CALL SINREV (X, Y, PRJ, PHI, THETA, IERR)
d569 1
a569 1
 105  CALL ARCREV (X, Y, PRJ, PHI, THETA, IERR)
d573 1
a573 1
 106  CALL ZPNREV (X, Y, PRJ, PHI, THETA, IERR)
d577 1
a577 1
 107  CALL ZEAREV (X, Y, PRJ, PHI, THETA, IERR)
d581 1
a581 1
 108  CALL AIRREV (X, Y, PRJ, PHI, THETA, IERR)
d658 2
a659 1
*      PRJ(1)    AZP distance parameter, mu in units of r0.
d663 3
a665 1
*      PRJ(11)   State flag.
d667 7
a673 3
*      PRJ(13)   1/PRJ(12)
*      PRJ(14)   Boundary parameter, -mu    for |mu| <= 1,
*                                    -1/mu  for |mu| >= 1.
d683 2
d686 2
d695 10
a704 1
      PRJ(13) = 1D0/PRJ(12)
d706 2
a707 2
      IF (ABS(PRJ(1)).LE.1D0) THEN
         PRJ(14) = -PRJ(1)
d709 1
a709 1
         PRJ(14) = -1D0/PRJ(1)
d712 4
a715 2
      IF (PRJ(11).EQ.-1D0) THEN
         PRJ(11) = -101D0
d717 1
a717 1
         PRJ(11) =  101D0
d728 2
a729 1
      DOUBLE PRECISION PHI, PRJ(0:20), R, S, STHE, THETA, X, Y
d731 1
a731 1
      DOUBLE PRECISION COSD, SIND
d738 6
a743 3
      STHE = SIND(THETA)
      S = PRJ(1) + STHE
      IF (S.EQ.0D0) THEN
d748 1
a748 1
      R =  PRJ(12)*COSD(THETA)/S
d750 9
a758 1
      Y = -R*COSD(PHI)
d760 19
a778 3
      IF (PRJ(11).GT.0D0 .AND. STHE.LT.PRJ(14)) THEN
         IERR = 2
         RETURN
d789 2
a790 1
      DOUBLE PRECISION PHI, PRJ(0:20), R, RHO, S, THETA, TOL, X, Y
d801 3
a803 1
      R = SQRT(X*X + Y*Y)
d805 92
a896 1
         PHI = 0D0
d898 30
a927 1
         PHI = ATAN2D(X, -Y)
d930 7
a936 4
      RHO = R*PRJ(13)
      S = RHO*PRJ(1)/SQRT(RHO*RHO+1D0)
      IF (ABS(S).GT.1D0) THEN
         IF (ABS(S).GT.1D0+TOL) THEN
d940 60
a999 1
         THETA = ATAN2D(1D0,RHO) - SIGN(90D0,S)
d1001 39
a1039 1
         THETA = ATAN2D(1D0,RHO) - ASIND(S)
d1042 2
d1053 1
a1053 1
*      PRJ(11)   State flag.
d1064 2
a1067 6
      IF (PRJ(11).EQ.-1D0) THEN
         PRJ(11) = -102D0
      ELSE
         PRJ(11) =  102D0
      END IF

d1080 1
a1080 1
      IF (ABS(PRJ(11)).NE.102D0) THEN
d1112 1
a1112 1
      IF (ABS(PRJ(11)).NE.102D0) THEN
d1134 3
a1136 1
*      PRJ(11)   State flag.
d1149 2
a1159 2
      PRJ(11) = 103D0

d1172 1
a1172 1
      IF (PRJ(11).NE.103D0) THEN
d1199 1
a1199 1
      IF (PRJ(11).NE.103D0) THEN
d1220 1
a1220 1
*      PRJ(1:2)  SIN obliqueness parameters, alpha and beta.
d1224 3
a1226 1
*      PRJ(11)   State flag.
d1228 3
a1230 4
*      PRJ(13)   alpha**2 + beta**2
*      PRJ(14)   2*(alpha**2 + beta**2)
*      PRJ(15)   2*(alpha**2 + beta**2 + 1)
*      PRJ(16)   alpha**2 + beta**2 - 1
d1241 2
d1247 2
a1248 9
      PRJ(14) = 2D0*PRJ(13)
      PRJ(15) = PRJ(14) + 2D0
      PRJ(16) = PRJ(13) - 1D0

      IF (PRJ(11).EQ.-1D0) THEN
         PRJ(11) = -104D0
      ELSE
         PRJ(11) =  104D0
      END IF
d1267 1
a1267 1
      IF (ABS(PRJ(11)).NE.104D0) THEN
d1275 1
a1275 1
            Z = -T*T/2D0
d1277 1
a1277 1
            Z = -2D0 + T*T/2D0
d1281 1
a1281 1
         Z = SIND(THETA) - 1D0
d1288 1
a1288 1
      Y = -PRJ(10)*(CTHE*CPHI + PRJ(2)*Z)
d1300 1
a1300 1
            T = ATAND(PRJ(1)*SPHI + PRJ(2)*CPHI)
d1316 2
a1317 2
      DOUBLE PRECISION A, B, C, D, PHI, PRJ(0:20), R2, STH, STH1, STH2,
     *          SXY, THETA, TOL, X, X0, XP, Y, Y0, YP, Z
d1327 1
a1327 1
      IF (ABS(PRJ(11)).NE.104D0) THEN
d1356 4
d1362 2
a1363 2
            Z = -R2/2D0
            THETA = 90D0 - R2D*SQRT(R2/(1D0 - X0*PRJ(1) + Y0*PRJ(2)))
d1366 4
a1369 1
            SXY = 2D0*(PRJ(1)*X0 - PRJ(2)*Y0)
a1370 5
            A = PRJ(15)
            B = -(SXY + PRJ(14))
            C = R2 + SXY + PRJ(16)
            D = B*B - 2D0*A*C

d1381 4
a1384 4
            STH = MAX(STH1,STH2)
            IF (STH.GT.1D0) THEN
               IF (STH-1D0.LT.TOL) THEN
                  STH = 1D0
d1386 1
a1386 1
                  STH = MIN(STH1,STH2)
d1389 5
a1393 2
            IF (STH.LT.-1D0) THEN
               IF (STH+1D0.GT.-TOL) STH = -1D0
d1395 2
a1396 1
            IF (STH.GT.1D0 .OR. STH.LT.-1D0) THEN
d1401 2
a1402 2
            THETA = ASIND(STH)
            Z = STH - 1D0
d1405 1
a1405 1
         XP = -Y0 - PRJ(2)*Z
d1423 3
a1425 1
*      PRJ(11)   State flag.
d1438 2
a1448 2
      PRJ(11) = 105D0

d1461 1
a1461 1
      IF (PRJ(11).NE.105D0) THEN
d1482 1
a1482 1
      IF (PRJ(11).NE.105D0) THEN
d1507 3
a1509 1
*      PRJ(11)   State flag.
d1525 2
a1591 6
      IF (PRJ(11).EQ.-1D0) THEN
         PRJ(11) = -106D0
      ELSE
         PRJ(11) =  106D0
      END IF

d1608 1
a1608 1
      IF (ABS(PRJ(11)).NE.106D0) THEN
d1648 1
a1648 1
      IF (ABS(PRJ(11)).NE.106D0) THEN
d1764 3
a1766 1
*      PRJ(11)   State flag.
d1779 2
a1789 2
      PRJ(11) = 107D0

d1802 1
a1802 1
      IF (PRJ(11).NE.107D0) THEN
d1825 1
a1825 1
      IF (PRJ(11).NE.107D0) THEN
d1862 3
a1864 1
*      PRJ(11)   State flag.
d1887 2
a1908 2
      PRJ(11) = 108D0

d1925 1
a1925 1
      IF (PRJ(11).NE.108D0) THEN
d1968 1
a1968 1
      IF (PRJ(11).NE.108D0) THEN
d2046 3
a2048 1
*      PRJ(11)   State flag.
d2063 2
a2100 2
      PRJ(11) = 201D0

d2161 3
a2163 1
*      PRJ(11)   State flag.
d2178 2
a2200 2
      PRJ(11) = 202D0

d2262 3
a2264 1
*      PRJ(11)   State flag.
d2277 2
a2287 2
      PRJ(11) = 203D0

d2333 3
a2335 1
*      PRJ(11)   State flag.
d2348 2
a2358 2
      PRJ(11) = 204D0

d2413 3
a2415 1
*      PRJ(11)   State flag.
d2428 2
a2438 2
      PRJ(11) = 301D0

d2491 3
a2493 1
*      PRJ(11)   State flag.
d2508 2
a2522 2
      PRJ(11) = 302D0

d2590 3
a2592 1
*      PRJ(11)   State flag.
d2611 2
a2620 2
      PRJ(11) = 303D0

d2629 1
a2629 1
      DOUBLE PRECISION ALPHA, PHI, PRJ(0:20), RESID, THETA, TOL, U,
d2668 3
a2670 3
 20      ALPHA = V/2D0
         X = PRJ(13)*PHI*COS(ALPHA)
         Y = PRJ(12)*SIN(ALPHA)
d2745 3
a2747 1
*      PRJ(11)   State flag.
d2762 2
a2770 2
      PRJ(11) = 401D0

d2779 1
a2779 1
      DOUBLE PRECISION COSTHE, PHI, PRJ(0:20), THETA, W, X, Y
d2788 3
a2790 3
      COSTHE = COSD(THETA)
      W = SQRT(PRJ(12)/(1D0 + COSTHE*COSD(PHI/2D0)))
      X = 2D0*W*COSTHE*SIND(PHI/2D0)
d2856 3
a2858 1
*      PRJ(11)   State flag.
d2877 2
a2899 6
      IF (PRJ(11).EQ.-1D0) THEN
         PRJ(11) = -501D0
      ELSE
         PRJ(11) =  501D0
      END IF

d2979 3
a2981 1
*      PRJ(11)   State flag.
d3003 2
a3026 2
      PRJ(11) = 502D0

d3116 3
a3118 1
*      PRJ(11)   State flag.
d3135 2
a3153 2
      PRJ(11) = 503D0

d3221 3
a3223 1
*      PRJ(11)   State flag.
d3243 2
a3274 2
      PRJ(11) = 504D0

d3357 3
a3359 1
*      PRJ(11)   State flag.
d3374 2
a3384 2
      PRJ(11) = 601D0

d3422 1
a3422 1
      DOUBLE PRECISION A, DY, COSTHE, PHI, PRJ(0:20), R, THETA, X, Y
d3448 2
a3449 2
      COSTHE = COSD(THETA)
      IF (COSTHE.EQ.0D0) THEN
d3452 1
a3452 1
         PHI = A*(R/PRJ(10))/COSD(THETA)
d3464 3
a3466 1
*      PRJ(11)   State flag.
d3480 2
a3492 2
      PRJ(11) = 602D0

d3501 2
a3502 2
      DOUBLE PRECISION A, COSTHE, COTTHE, PHI, PRJ(0:20), SINTHE,
     *          THETA, X, Y
d3515 3
a3517 3
      COSTHE = COSD(THETA)
      SINTHE = SIND(THETA)
      A = PHI*SINTHE
d3519 1
a3519 1
      IF (SINTHE.EQ.0D0) THEN
d3523 1
a3523 1
         COTTHE = COSTHE/SINTHE
d3621 3
a3623 1
*      PRJ(11)   State flag.
d3636 2
a3646 2
      PRJ(11) = 701D0

d3655 1
a3655 1
      DOUBLE PRECISION COSTHE, L, M, N, PHI, PRJ(0:20), RHO, THETA, TOL,
d3667 3
a3669 3
      COSTHE = COSD(THETA)
      L = COSTHE*COSD(PHI)
      M = COSTHE*SIND(PHI)
d3835 3
a3837 1
*      PRJ(11)   State flag.
d3850 2
a3860 2
      PRJ(11) = 702D0

d3869 1
a3869 1
      DOUBLE PRECISION COSTHE, ETA, L, M, N, PHI, PRJ(0:20), RHO, THETA,
d3898 3
a3900 3
      COSTHE = COSD(THETA)
      L = COSTHE*COSD(PHI)
      M = COSTHE*SIND(PHI)
d4152 3
a4154 1
*      PRJ(11)   State flag.
d4167 2
a4177 2
      PRJ(11) = 703D0

d4186 2
a4187 2
      DOUBLE PRECISION CHI, COSTHE, ETA, L, M, N, P, PHI, PRJ(0:20),
     *          PSI, RHO, RHU, T, THETA, TOL, X, X0, XF, XI, Y, Y0, YF
d4209 3
a4211 3
      COSTHE = COSD(THETA)
      L = COSTHE*COSD(PHI)
      M = COSTHE*SIND(PHI)
d4318 4
a4321 4
         PSI = ETA/XI
         CHI = 1D0 + PSI*PSI
         XF  = -SQRT(RHU/(1D0-1D0/SQRT(1D0+CHI)))
         YF  = (XF/15D0)*(ATAND(PSI) - ASIND(PSI/SQRT(CHI+CHI)))
d4323 4
a4326 4
         PSI = ETA/XI
         CHI = 1D0 + PSI*PSI
         XF  =  SQRT(RHU/(1D0-1D0/SQRT(1D0+CHI)))
         YF  = (XF/15D0)*(ATAND(PSI) - ASIND(PSI/SQRT(CHI+CHI)))
d4328 4
a4331 4
         PSI = XI/ETA
         CHI = 1D0 + PSI*PSI
         YF  = -SQRT(RHU/(1D0-1D0/SQRT(1D0+CHI)))
         XF  = (YF/15D0)*(ATAND(PSI) - ASIND(PSI/SQRT(CHI+CHI)))
d4333 4
a4336 4
         PSI = XI/ETA
         CHI = 1D0 + PSI*PSI
         YF  =  SQRT(RHU/(1D0-1D0/SQRT(1D0+CHI)))
         XF  = (YF/15D0)*(ATAND(PSI) - ASIND(PSI/SQRT(CHI+CHI)))
d4367 1
a4367 1
      DOUBLE PRECISION CHI, L, M, N, PHI, PRJ(0:20), PSI, RHO, RHU,
d4424 2
a4425 2
            PSI = 0D0
            CHI = 1D0
d4430 3
a4432 3
            PSI = SIND(W)/(COSD(W) - C2)
            CHI = 1D0 + PSI*PSI
            RHU = XF*XF*(1D0 - 1D0/SQRT(1D0 + CHI))
d4437 2
a4438 2
            PSI = 0D0
            CHI = 1D0
d4443 3
a4445 3
            PSI = SIND(W)/(COSD(W) - C2)
            CHI = 1D0 + PSI*PSI
            RHU = YF*YF*(1D0 - 1D0/SQRT(1D0 + CHI))
d4460 1
a4460 1
         W = SQRT(RHU*(2D0-RHU)/CHI)
d4468 1
a4468 1
            L = -M*PSI
d4472 1
a4472 1
            M = -L*PSI
d4479 1
a4479 1
            N = M*PSI
d4483 1
a4483 1
            M = N*PSI
d4490 1
a4490 1
            N = -L*PSI
d4494 1
a4494 1
            L = -N*PSI
d4501 1
a4501 1
            N = -M*PSI
d4505 1
a4505 1
            M = -N*PSI
d4512 1
a4512 1
            N = L*PSI
d4516 1
a4516 1
            L = N*PSI
d4523 1
a4523 1
            L = M*PSI
d4527 1
a4527 1
            M = L*PSI
@


2.14
log
@Provide driver routines PRJSET, PRJFWD & PRJREV which key on the projection-
specific value of PRJ(11); renamed glsset(), glsfwd() & glsrev() to sflset(),
sflfwd() & sflrev(); reordered projections to match paper.
@
text
@d220 1
a220 1
*   $Id: proj.f,v 2.14 2000/12/05 04:00:30 mcalabre Exp mcalabre $
d231 1
a231 1
      DATA IREF /2,2,2,2,2,2,2,2,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0/
@


2.13
log
@Updated date in copyright notice.
@
text
@d42 2
d46 1
a47 1
*      STGSET STGFWD STGREV   STG: stereographic
d53 1
d56 4
a59 1
*      CEASET CEAFWD CEAREV   CEA: cylindrical equal area
d66 1
a66 4
*      GLSSET GLSFWD GLSREV   GLS: Sanson-Flamsteed (global sinusoidal)
*      PARSET PARFWD PARREV   PAR: parabolic
*      AITSET AITFWD AITREV   AIT: Hammer-Aitoff
*      MOLSET MOLFWD MOLREV   MOL: Molweide
d69 28
a96 1
*      TSCSET TSCFWD TSCREV   TSC: tangential spherical cube
d102 2
a103 1
*   invoked by PRJFWD and PRJREV if PRJ(11) is zero.
d170 2
a171 1
*         COP projections.
d220 420
a639 1
*   $Id: proj.f,v 2.12 2000/02/01 04:17:25 mcalabre Exp mcalabre $
d679 1
a679 1
         PRJ(11) = -137D0
d681 1
a681 1
         PRJ(11) = 137D0
d696 1
a696 1
      IF (ABS(PRJ(11)).NE.137D0) THEN
d712 1
a712 1
      IF (PRJ(11).EQ.137D0 .AND. STHE.LT.PRJ(14)) THEN
d731 1
a731 1
      IF (ABS(PRJ(11)).NE.137D0) THEN
d778 1
a778 1
         PRJ(11) = -137D0
d780 1
a780 1
         PRJ(11) = 137D0
d795 1
a795 1
      IF (ABS(PRJ(11)).NE.137D0) THEN
d810 1
a810 1
      IF (PRJ(11).EQ.137D0 .AND. S.LT.0D0) THEN
d827 1
a827 1
      IF (ABS(PRJ(11)).NE.137D0) THEN
d845 85
d962 1
a962 1
         PRJ(11) = -137D0
d964 1
a964 1
         PRJ(11) = 137D0
d984 1
a984 1
      IF (ABS(PRJ(11)).NE.137D0) THEN
d1008 1
a1008 1
      IF (PRJ(11).EQ.137D0) THEN
d1044 1
a1044 1
      IF (ABS(PRJ(11)).NE.137D0) THEN
a1129 85
*   STG: stereographic projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*      PRJ(11)   State flag.
*      PRJ(12)   2*r0
*      PRJ(13)   1/(2*r0)
*=======================================================================
      SUBROUTINE STGSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D
         PRJ(12) = 360D0/PI
         PRJ(13) = PI/360D0
      ELSE
         PRJ(12) = 2D0*PRJ(10)
         PRJ(13) = 1D0/PRJ(12)
      END IF

      PRJ(11) = 137D0

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE STGFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), R, S, THETA, X, Y

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.137D0) THEN
         CALL STGSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      S = 1D0 + SIND(THETA)
      IF (S.EQ.0D0) THEN
         IERR = 2
         RETURN
      END IF

      R =  PRJ(12)*COSD(THETA)/S
      X =  R*SIND(PHI)
      Y = -R*COSD(PHI)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE STGREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), R, THETA, X, Y

      DOUBLE PRECISION ATAN2D, ATAND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.137D0) THEN
         CALL STGSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      R = SQRT(X*X + Y*Y)
      IF (R.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = ATAN2D(X, -Y)
      END IF
      THETA = 90D0 - 2D0*ATAND(R*PRJ(13))

      IERR = 0
      RETURN
      END

*=======================================================================
d1156 1
a1156 1
      PRJ(11) = 137D0
d1170 1
a1170 1
      IF (PRJ(11).NE.137D0) THEN
d1191 1
a1191 1
      IF (PRJ(11).NE.137D0) THEN
d1298 1
a1298 1
         PRJ(11) = -137D0
d1300 1
a1300 1
         PRJ(11) = 137D0
d1319 1
a1319 1
      IF (ABS(PRJ(11)).NE.137D0) THEN
d1335 1
a1335 1
      IF (PRJ(11).EQ.137D0 .AND. S.GT.PRJ(13)) THEN
d1359 1
a1359 1
      IF (ABS(PRJ(11)).NE.137D0) THEN
d1497 1
a1497 1
      PRJ(11) = 137D0
d1511 1
a1511 1
      IF (PRJ(11).NE.137D0) THEN
d1534 1
a1534 1
      IF (PRJ(11).NE.137D0) THEN
d1614 1
a1614 1
      PRJ(11) = 137D0
d1632 1
a1632 1
      IF (PRJ(11).NE.137D0) THEN
d1675 1
a1675 1
      IF (PRJ(11).NE.137D0) THEN
d1804 1
a1804 1
      PRJ(11) = 137D0
d1818 1
a1818 1
      IF (PRJ(11).NE.137D0) THEN
d1844 1
a1844 1
      IF (PRJ(11).NE.137D0) THEN
d1858 103
d1987 1
a1987 1
      PRJ(11) = 137D0
d1999 1
a1999 1
      IF (PRJ(11).NE.137D0) THEN
d2017 1
a2017 1
      IF (PRJ(11).NE.137D0) THEN
d2056 1
a2056 1
      PRJ(11) = 137D0
d2070 1
a2070 1
      IF (PRJ(11).NE.137D0) THEN
d2095 1
a2095 1
      IF (PRJ(11).NE.137D0) THEN
d2108 1
a2108 1
*   CEA: cylindrical equal area projection.
d2110 75
a2184 3
*   Given:
*      PRJ(1)    Square of the cosine of the latitude at which the
*                projection is conformal, lambda.
d2191 2
a2192 2
*      PRJ(14)   r0/lambda
*      PRJ(15)   lambda/r0
d2194 1
a2194 1
      SUBROUTINE CEASET (PRJ, IERR)
d2207 205
a2411 2
         IF (PRJ(1).LE.0D0 .OR. PRJ(1).GT.1D0) THEN
            IERR = 1
d2414 1
a2414 2
         PRJ(14) = PRJ(10)/PRJ(1)
         PRJ(15) = PRJ(1)/PRJ(10)
d2416 6
a2421 4
         PRJ(12) = PRJ(10)*D2R
         PRJ(13) = R2D/PRJ(10)
         IF (PRJ(1).LE.0D0 .OR. PRJ(1).GT.1D0) THEN
            IERR = 1
d2424 1
a2424 2
         PRJ(14) = PRJ(10)/PRJ(1)
         PRJ(15) = PRJ(1)/PRJ(10)
d2427 34
a2460 1
      PRJ(11) = 137D0
d2467 1
a2467 1
      SUBROUTINE CEAFWD (PHI, THETA, PRJ, X, Y, IERR)
d2470 1
a2470 1
      DOUBLE PRECISION PHI, PRJ(0:20), THETA, X, Y
d2472 1
a2472 1
      DOUBLE PRECISION SIND
d2474 2
a2475 2
      IF (PRJ(11).NE.137D0) THEN
         CALL CEASET (PRJ, IERR)
d2479 4
a2482 2
      X = PRJ(12)*PHI
      Y = PRJ(14)*SIND(THETA)
d2489 1
a2489 1
      SUBROUTINE CEAREV (X, Y, PRJ, PHI, THETA, IERR)
d2492 1
a2492 1
      DOUBLE PRECISION PHI, PRJ(0:20), S, THETA, TOL, X, Y
d2496 1
a2496 1
      DOUBLE PRECISION ASIND
d2498 2
a2499 2
      IF (PRJ(11).NE.137D0) THEN
         CALL CEASET (PRJ, IERR)
d2503 12
a2514 1
      S = Y*PRJ(15)
d2524 7
a2530 1
      PHI   = X*PRJ(13)
d2588 1
a2588 1
         PRJ(11) = -137D0
d2590 1
a2590 1
         PRJ(11) = 137D0
d2605 1
a2605 1
      IF (ABS(PRJ(11)).NE.137D0) THEN
d2623 1
a2623 1
      IF (PRJ(11).EQ.137D0 .AND. R*PRJ(12).LT.0D0) THEN
d2640 1
a2640 1
      IF (ABS(PRJ(11)).NE.137D0) THEN
d2663 1
a2663 1
*   COD: conic equidistant projection.
d2673 1
a2673 1
*      PRJ(12)   C = r0*sin(sigma)*sin(delta)/delta
d2675 7
a2681 2
*      PRJ(14)   Y0 = delta*cot(delta)*cot(sigma)
*      PRJ(15)   Y0 + sigma
d2683 1
a2683 1
      SUBROUTINE CODSET (PRJ, IERR)
d2686 1
a2686 1
      DOUBLE PRECISION PRJ(0:20)
d2692 1
a2692 1
      DOUBLE PRECISION COSD, SIND
d2696 2
a2697 5
      IF (PRJ(2).EQ.0D0) THEN
         PRJ(12) = PRJ(10)*SIND(PRJ(1))*D2R
      ELSE
         PRJ(12) = PRJ(10)*SIND(PRJ(1))*SIND(PRJ(2))/PRJ(2)
      END IF
d2699 1
d2706 9
a2714 2
      PRJ(14) = PRJ(10)*COSD(PRJ(2))*COSD(PRJ(1))/PRJ(12)
      PRJ(15) = PRJ(14) + PRJ(1)
d2716 1
a2716 1
      PRJ(11) = 137D0
d2723 1
a2723 1
      SUBROUTINE CODFWD (PHI, THETA, PRJ, X, Y, IERR)
d2730 2
a2731 2
      IF (PRJ(11).NE.137D0) THEN
         CALL CODSET (PRJ, IERR)
d2735 6
a2740 2
      A = PRJ(12)*PHI
      R = PRJ(15) - THETA
d2750 1
a2750 1
      SUBROUTINE CODREV (X, Y, PRJ, PHI, THETA, IERR)
d2753 3
a2755 1
      DOUBLE PRECISION A, DY, PHI, PRJ(0:20), R, THETA, X, Y
d2757 1
a2757 1
      DOUBLE PRECISION ATAN2D
d2759 2
a2760 2
      IF (PRJ(11).NE.137D0) THEN
         CALL CODSET (PRJ, IERR)
d2774 18
a2791 2
      PHI   = A*PRJ(13)
      THETA = PRJ(15) - R
d2798 1
a2798 1
*   COE: conic equal area projection.
d2808 1
a2808 1
*      PRJ(12)   C = (sin(theta1) + sin(theta2))/2
d2810 2
a2811 7
*      PRJ(14)   Y0 = chi*sqrt(psi - 2C*sind(sigma))
*      PRJ(15)   chi = r0/C
*      PRJ(16)   psi = 1 + sin(theta1)*sin(theta2)
*      PRJ(17)   2C
*      PRJ(18)   (1 + sin(theta1)*sin(theta2))*(r0/C)**2
*      PRJ(19)   C/(2*r0**2)
*      PRJ(20)   Y0 = chi*sqrt(psi + 2C)
d2813 1
a2813 1
      SUBROUTINE COESET (PRJ, IERR)
d2816 1
a2816 1
      DOUBLE PRECISION PRJ(0:20), THETA1, THETA2
d2822 1
a2822 1
      DOUBLE PRECISION SIND
d2826 5
a2830 2
      THETA1 = PRJ(1) - PRJ(2)
      THETA2 = PRJ(1) + PRJ(2)
a2831 1
      PRJ(12) = (SIND(THETA1) + SIND(THETA2))/2D0
d2838 2
a2839 7

      PRJ(15) = PRJ(10)/PRJ(12)
      PRJ(16) = 1D0 + SIND(THETA1)*SIND(THETA2)
      PRJ(17) = 2D0*PRJ(12)
      PRJ(18) = PRJ(15)*PRJ(15)*PRJ(16)
      PRJ(19) = 1D0/(2D0*PRJ(10)*PRJ(15))
      PRJ(20) = PRJ(15)*SQRT(PRJ(16) + PRJ(17))
d2841 1
a2841 3
      PRJ(14) = PRJ(15)*SQRT(PRJ(16) - PRJ(17)*SIND(PRJ(1)))

      PRJ(11) = 137D0
d2848 1
a2848 1
      SUBROUTINE COEFWD (PHI, THETA, PRJ, X, Y, IERR)
d2855 2
a2856 2
      IF (PRJ(11).NE.137D0) THEN
         CALL COESET (PRJ, IERR)
d2860 2
a2861 6
      A = PHI*PRJ(12)
      IF (THETA.EQ.-90D0) THEN
         R = PRJ(20)
      ELSE
         R = PRJ(15)*SQRT(PRJ(16) - PRJ(17)*SIND(THETA))
      END IF
d2871 1
a2871 1
      SUBROUTINE COEREV (X, Y, PRJ, PHI, THETA, IERR)
d2874 1
a2874 1
      DOUBLE PRECISION A, DY, PHI, PRJ(0:20), R, THETA, TOL, W, X, Y
d2876 1
a2876 3
      PARAMETER (TOL = 1D-12)

      DOUBLE PRECISION ASIND, ATAN2D
d2878 2
a2879 2
      IF (PRJ(11).NE.137D0) THEN
         CALL COESET (PRJ, IERR)
d2893 2
a2894 18
      PHI = A*PRJ(13)
      IF (ABS(R - PRJ(20)).LT.TOL) THEN
         THETA = -90D0
      ELSE
         W = (PRJ(18) - R*R)*PRJ(19)
         IF (ABS(W).GT.1D0) THEN
            IF (ABS(W-1D0).LT.TOL) THEN
               THETA = 90D0
            ELSE IF (ABS(W+1D0).LT.TOL) THEN
               THETA = -90D0
            ELSE
               IERR = 2
               RETURN
            END IF
         ELSE
            THETA = ASIND(W)
         END IF
      END IF
d2960 1
a2960 1
      PRJ(11) = 137D0
d2974 1
a2974 1
      IF (PRJ(11).NE.137D0) THEN
d3006 1
a3006 1
      IF (PRJ(11).NE.137D0) THEN
d3068 1
a3068 1
      PRJ(11) = 137D0
d3084 1
a3084 1
         CALL GLSFWD (PHI, THETA, PRJ, X, Y, IERR)
d3088 1
a3088 1
      IF (PRJ(11).NE.137D0) THEN
d3113 1
a3113 1
         CALL GLSREV (X, Y, PRJ, PHI, THETA, IERR)
d3117 1
a3117 1
      IF (PRJ(11).NE.137D0) THEN
d3174 1
a3174 1
      PRJ(11) = 137D0
d3193 1
a3193 1
      IF (PRJ(11).NE.137D0) THEN
d3226 1
a3226 1
      IF (PRJ(11).NE.137D0) THEN
d3300 1
a3300 1
*   GLS: Sanson-Flamsteed ("global sinusoid") projection.
d3305 2
a3306 2
*      PRJ(12)   r0*(pi/180)
*      PRJ(13)   (180/pi)/r0
d3308 1
a3308 1
      SUBROUTINE GLSSET (PRJ, IERR)
d3319 2
a3320 2
         PRJ(12) = 1D0
         PRJ(13) = 1D0
d3322 1
a3322 1
         PRJ(12) = PRJ(10)*D2R
d3326 1
a3326 1
      PRJ(11) = 137D0
d3333 1
a3333 1
      SUBROUTINE GLSFWD (PHI, THETA, PRJ, X, Y, IERR)
d3335 5
a3339 2
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), THETA, X, Y
d3341 1
a3341 1
      DOUBLE PRECISION COSD
d3343 2
a3344 2
      IF (PRJ(11).NE.137D0) THEN
         CALL GLSSET (PRJ, IERR)
d3348 4
a3351 2
      X = PRJ(12)*PHI*COSD(THETA)
      Y = PRJ(12)*THETA
d3353 5
a3357 13
      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE GLSREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), THETA, W, X, Y
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.137D0) THEN
         CALL GLSSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
d3359 3
a3361 6

      W = COS(Y/PRJ(10))
      IF (W.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = X*PRJ(13)/COS(Y/PRJ(10))
d3363 3
a3365 37
      THETA = Y*PRJ(13)

      IERR = 0
      RETURN
      END

*=======================================================================
*   PAR: parabolic projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*      PRJ(11)   State flag.
*      PRJ(12)   r0*(pi/180)
*      PRJ(13)   (180/pi)/r0
*      PRJ(14)   pi*r0
*      PRJ(15)   1/(pi*r0)
*=======================================================================
      SUBROUTINE PARSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D
         PRJ(12) = 1D0
         PRJ(13) = 1D0
         PRJ(14) = 180D0
         PRJ(15) = 1D0/PRJ(14)
      ELSE
         PRJ(12) = PRJ(10)*D2R
         PRJ(13) = 1D0/PRJ(12)
         PRJ(14) = PI*PRJ(10)
         PRJ(15) = 1D0/PRJ(14)
d3367 3
a3369 18

      PRJ(11) = 137D0

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE PARFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), S, THETA, X, Y

      DOUBLE PRECISION SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.137D0) THEN
         CALL PARSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
d3371 3
a3373 20

      S = SIND(THETA/3D0)
      X = PRJ(12)*PHI*(1D0 - 4D0*S*S)
      Y = PRJ(14)*S

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE PARREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), S, T, THETA, X, Y

      DOUBLE PRECISION ASIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.137D0) THEN
         CALL PARSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
d3376 30
a3405 4
      S = Y*PRJ(15)
      IF (S.GT.1D0 .OR. S.LT.-1D0) THEN
         IERR = 2
         RETURN
d3408 2
a3409 5
      T = 1D0 - 4D0*S*S
      IF (T.EQ.0D0) THEN
         IF (X.EQ.0D0) THEN
            PHI = 0D0
         ELSE
d3413 1
a3413 2
      ELSE
         PHI = PRJ(13)*X/T
d3415 6
a3420 51

      THETA = 3D0*ASIND(S)

      IERR = 0
      RETURN
      END

*=======================================================================
*   AIT: Hammer-Aitoff projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*      PRJ(11)   State flag.
*      PRJ(12)   2*r0**2
*      PRJ(13)   1/(2*r0)**2
*      PRJ(14)   1/(4*r0)**2
*      PRJ(15)   1/(2*r0)
*=======================================================================
      SUBROUTINE AITSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      IF (PRJ(10).EQ.0D0) PRJ(10) = R2D

      PRJ(12) = 2D0*PRJ(10)*PRJ(10)
      PRJ(13) = 1D0/(2D0*PRJ(12))
      PRJ(14) = PRJ(13)/4D0
      PRJ(15) = 1D0/(2D0*PRJ(10))

      PRJ(11) = 137D0

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE AITFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION COSTHE, PHI, PRJ(0:20), THETA, W, X, Y

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.137D0) THEN
         CALL AITSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
d3423 2
a3424 4
      COSTHE = COSD(THETA)
      W = SQRT(PRJ(12)/(1D0 + COSTHE*COSD(PHI/2D0)))
      X = 2D0*W*COSTHE*SIND(PHI/2D0)
      Y = W*SIND(THETA)
d3431 1
a3431 1
      SUBROUTINE AITREV (X, Y, PRJ, PHI, THETA, IERR)
d3434 1
a3434 3
      DOUBLE PRECISION PHI, PRJ(0:20), S, THETA, TOL, U, X, XP, Y, YP, Z

      PARAMETER (TOL = 1D-13)
d3438 2
a3439 2
      IF (PRJ(11).NE.137D0) THEN
         CALL AITSET (PRJ, IERR)
d3443 6
a3448 3
      U = 1D0 - X*X*PRJ(14) - Y*Y*PRJ(13)
      IF (U.LT.0D0) THEN
         IF (U.LT.-TOL) THEN
d3452 2
a3453 8

         U = 0D0
      END IF

      Z = SQRT(U)
      S = Z*Y/PRJ(10)
      IF (ABS(S).GT.1D0) THEN
         IF (ABS(S).GT.1D0+TOL) THEN
a3456 2

         S = SIGN(1D0,S)
d3459 2
a3460 8
      XP = 2D0*Z*Z - 1D0
      YP = Z*X*PRJ(15)
      IF (XP.EQ.0D0 .AND. YP.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = 2D0*ATAN2D(YP, XP)
      END IF
      THETA = ASIND(S)
d3462 31
a3492 68
      IERR = 0
      RETURN
      END

*=======================================================================
*   MOL: Mollweide's projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*      PRJ(11)   State flag.
*      PRJ(12)   sqrt(2)*r0
*      PRJ(13)   sqrt(2)*r0/90
*      PRJ(14)   1/(sqrt(2)*r0)
*      PRJ(15)   90/r0
*      PRJ(16)   2/pi
*=======================================================================
      SUBROUTINE MOLSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      DOUBLE PRECISION SQRT2
      PARAMETER (SQRT2 = 1.4142135623730950488D0)
*-----------------------------------------------------------------------
      IF (PRJ(10).EQ.0D0) PRJ(10) = R2D

      PRJ(12) = SQRT2*PRJ(10)
      PRJ(13) = PRJ(12)/90D0
      PRJ(14) = 1D0/PRJ(12)
      PRJ(15) = 90D0/PRJ(10)
      PRJ(16) = 2D0/PI

      PRJ(11) = 137D0

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE MOLFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR, J
      DOUBLE PRECISION ALPHA, PHI, PRJ(0:20), RESID, THETA, TOL, U,
     *          V, V0, V1, X, Y

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      PARAMETER (TOL = 1D-13)

      DOUBLE PRECISION SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.137D0) THEN
         CALL MOLSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      IF (ABS(THETA).EQ.90D0) THEN
         X = 0D0
         Y = SIGN(PRJ(12),THETA)
      ELSE IF (THETA.EQ.0D0) THEN
         X = PRJ(13)*PHI
         Y = 0D0
d3494 4
a3497 19
         U  = PI*SIND(THETA)
         V0 = -PI
         V1 = +PI
         V  = U
         DO 10 J = 1, 100
            RESID = (V - U) + SIN(V)
            IF (RESID.LT.0D0) THEN
               IF (RESID.GT.-TOL) GO TO 20
               V0 = V
            ELSE
               IF (RESID.LT.TOL) GO TO 20
               V1 = V
            END IF
            V = (V0 + V1)/2D0
 10      CONTINUE

 20      ALPHA = V/2D0
         X = PRJ(13)*PHI*COS(ALPHA)
         Y = PRJ(12)*SIN(ALPHA)
d3500 1
a3500 36
      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE MOLREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PHI, PRJ(0:20), S, THETA, TOL, X, Y, Y0, Z

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)

      PARAMETER (TOL = 1D-12)

      DOUBLE PRECISION ASIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.137D0) THEN
         CALL MOLSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      Y0 = Y/PRJ(10)
      S  = 2D0 - Y0*Y0
      IF (S.LE.TOL) THEN
         IF (S.LT.-TOL) THEN
            IERR = 2
            RETURN
         END IF
         S = 0D0

         IF (ABS(X).GT.TOL) THEN
            IERR = 2
            RETURN
         END IF
d3503 1
a3503 2
         S = SQRT(S)
         PHI = PRJ(15)*X/S
d3505 1
a3505 21

      Z = Y*PRJ(14)
      IF (ABS(Z).GT.1D0) THEN
         IF (ABS(Z).GT.1D0+TOL) THEN
            IERR = 2
            RETURN
         END IF
         Z = SIGN(1D0,Z) + Y0*S/PI
      ELSE
         Z = ASIN(Z)*PRJ(16) + Y0*S/PI
      END IF

      IF (ABS(Z).GT.1D0) THEN
         IF (ABS(Z).GT.1D0+TOL) THEN
            IERR = 2
            RETURN
         END IF
         Z = SIGN(1D0,Z)
      END IF

      THETA = ASIND(Z)
d3538 1
a3538 1
      PRJ(11) = 137D0
d3572 1
a3572 1
      IF (PRJ(11).NE.137D0) THEN
d3724 1
a3724 1
      IF (PRJ(11).NE.137D0) THEN
d3853 1
a3853 1
      PRJ(11) = 137D0
d3874 1
a3874 1
      IF (PRJ(11).NE.137D0) THEN
d4054 1
a4054 1
      IF (PRJ(11).NE.137D0) THEN
a4205 212
      END IF

      IF (L.EQ.0D0 .AND. M.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = ATAN2D(M, L)
      END IF
      THETA = ASIND(N)

      IERR = 0
      RETURN
      END

*=======================================================================
*   TSC: tangential spherical cube projection.
*
*   Given and/or returned:
*      PRJ(10)   r0; reset to 180/pi if 0.
*      PRJ(11)   State flag.
*      PRJ(12)   r0*(pi/4)
*      PRJ(13)   (4/pi)/r0
*=======================================================================
      SUBROUTINE TSCSET (PRJ, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION PRJ(0:20)

      DOUBLE PRECISION D2R, PI, R2D
      PARAMETER (PI = 3.141592653589793238462643D0)
      PARAMETER (D2R = PI/180D0, R2D = 180D0/PI)
*-----------------------------------------------------------------------
      IF (PRJ(10).EQ.0D0) THEN
         PRJ(10) = R2D
         PRJ(12) = 45D0
         PRJ(13) = 1D0/45D0
      ELSE
         PRJ(12) = PRJ(10)*PI/4D0
         PRJ(13) = 1D0/PRJ(12)
      END IF

      PRJ(11) = 137D0

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE TSCFWD (PHI, THETA, PRJ, X, Y, IERR)
*-----------------------------------------------------------------------
      INTEGER   FACE, IERR
      DOUBLE PRECISION COSTHE, L, M, N, PHI, PRJ(0:20), RHO, THETA, TOL,
     *          X, X0, XF, Y, Y0, YF

      PARAMETER (TOL = 1D-12)

      DOUBLE PRECISION COSD, SIND
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.137D0) THEN
         CALL TSCSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      COSTHE = COSD(THETA)
      L = COSTHE*COSD(PHI)
      M = COSTHE*SIND(PHI)
      N = SIND(THETA)

      FACE = 0
      RHO  = N
      IF (L.GT.RHO) THEN
         FACE = 1
         RHO  = L
      END IF
      IF (M.GT.RHO) THEN
         FACE = 2
         RHO  = M
      END IF
      IF (-L.GT.RHO) THEN
         FACE = 3
         RHO  = -L
      END IF
      IF (-M.GT.RHO) THEN
         FACE = 4
         RHO  = -M
      END IF
      IF (-N.GT.RHO) THEN
         FACE = 5
         RHO  = -N
      END IF

      IF (FACE.EQ.0) THEN
         XF =  M/RHO
         YF = -L/RHO
         X0 =  0D0
         Y0 =  2D0
      ELSE IF (FACE.EQ.1) THEN
         XF =  M/RHO
         YF =  N/RHO
         X0 =  0D0
         Y0 =  0D0
      ELSE IF (FACE.EQ.2) THEN
         XF = -L/RHO
         YF =  N/RHO
         X0 =  2D0
         Y0 =  0D0
      ELSE IF (FACE.EQ.3) THEN
         XF = -M/RHO
         YF =  N/RHO
         X0 =  4D0
         Y0 =  0D0
      ELSE IF (FACE.EQ.4) THEN
         XF =  L/RHO
         YF =  N/RHO
         X0 =  6D0
         Y0 =  0D0
      ELSE IF (FACE.EQ.5) THEN
         XF =  M/RHO
         YF =  L/RHO
         X0 =  0D0
         Y0 = -2D0
      END IF

      IF (ABS(XF).GT.1D0) THEN
         IF (ABS(XF).GT.1D0+TOL) THEN
            IERR = 2
            RETURN
         END IF
         XF = SIGN(1D0,XF)
      END IF
      IF (ABS(YF).GT.1D0) THEN
         IF (ABS(YF).GT.1D0+TOL) THEN
            IERR = 2
            RETURN
         END IF
         YF = SIGN(1D0,YF)
      END IF

      X = PRJ(12)*(XF + X0)
      Y = PRJ(12)*(YF + Y0)

      IERR = 0
      RETURN
      END

*-----------------------------------------------------------------------
      SUBROUTINE TSCREV (X, Y, PRJ, PHI, THETA, IERR)
*-----------------------------------------------------------------------
      INTEGER   IERR
      DOUBLE PRECISION L, M, N, PHI, PRJ(0:20), THETA, X, XF, Y, YF

      DOUBLE PRECISION ASIND, ATAN2D
*-----------------------------------------------------------------------
      IF (PRJ(11).NE.137D0) THEN
         CALL TSCSET (PRJ, IERR)
         IF (IERR.NE.0) RETURN
      END IF

      XF = X*PRJ(13)
      YF = Y*PRJ(13)

*     Check bounds.
      IF (ABS(XF).LE.1D0) THEN
         IF (ABS(YF).GT.3D0) THEN
            IERR = 2
            RETURN
         END IF
      ELSE
         IF (ABS(XF).GT.7D0 .OR. ABS(YF).GT.1D0) THEN
            IERR = 2
            RETURN
         END IF
      END IF

*     Map negative faces to the other side.
      IF (XF.LT.-1D0) XF = XF + 8D0

*     Determine the face.
      IF (XF.GT.5D0) THEN
*        FACE = 4
         XF = XF - 6D0
         M  = -1D0/SQRT(1D0 + XF*XF + YF*YF)
         L  = -M*XF
         N  = -M*YF
      ELSE IF (XF.GT.3D0) THEN
*        FACE = 3
         XF = XF - 4D0
         L  = -1D0/SQRT(1D0 + XF*XF + YF*YF)
         M  =  L*XF
         N  = -L*YF
      ELSE IF (XF.GT.1D0) THEN
*        FACE = 2
         XF = XF - 2D0
         M  =  1D0/SQRT(1D0 + XF*XF + YF*YF)
         L  = -M*XF
         N  =  M*YF
      ELSE IF (YF.GT.1D0) THEN
*        FACE = 0
         YF = YF - 2D0
         N  = 1D0/SQRT(1D0 + XF*XF + YF*YF)
         L  = -N*YF
         M  =  N*XF
      ELSE IF (YF.LT.-1D0) THEN
*        FACE = 5
         YF = YF + 2D0
         N  = -1D0/SQRT(1D0 + XF*XF + YF*YF)
         L  = -N*YF
         M  = -N*XF
      ELSE
*        FACE = 1
         L  =  1D0/SQRT(1D0 + XF*XF + YF*YF)
         M  =  L*XF
         N  =  L*YF
@


2.12
log
@Check for invalid (X,Y) in ZEAREV.
@
text
@d4 1
a4 1
*   Copyright (C) 1995-1999, Mark Calabretta
d188 1
a188 1
*   $Id: proj.f,v 2.11 1999/12/13 03:05:34 mcalabre Exp mcalabre $
@


2.11
log
@Amended "RETURN 2" in COOREV.
@
text
@d188 1
a188 1
*   $Id: proj.f,v 2.10 1999/11/19 05:54:32 mcalabre Exp mcalabre $
d1077 1
a1077 1
      DOUBLE PRECISION PHI, PRJ(0:20), R, THETA, TOL, X, Y
d1095 8
a1102 2
      IF (ABS(R-PRJ(12)).LT.TOL) THEN
         THETA = -90D0
@


2.10
log
@Allow more general face layout in the inverse quadcube projections.
@
text
@d188 1
a188 1
*   $Id: proj.f,v 2.9 1999/11/19 01:17:18 mcalabre Exp mcalabre $
d2139 1
a2139 1
            RETURN 2
@


2.9
log
@Declare ACOSD in SINREV; mMade a small allowance for rounding error in
determining the validity of (X,Y) in CEAREV and AITREV; corrected a comment
in the BON prologue.
@
text
@d188 1
a188 1
*   $Id: proj.f,v 2.8 1999/04/09 05:09:27 mcalabre Exp mcalabre $
d3063 8
a3070 6
*     Determine the face.
      IF (XF.GT.7.0) THEN
         IERR = 2
         RETURN
      ELSE IF (XF.GT.5.0) THEN
         IF (ABS(YF).GT.1.0) THEN
d3074 7
a3083 4
         IF (ABS(YF).GT.1.0) THEN
            IERR = 2
            RETURN
         END IF
a3086 4
         IF (ABS(YF).GT.1.0) THEN
            IERR = 2
            RETURN
         END IF
a3088 3
      ELSE IF (XF.LT.-1.0) THEN
         IERR = 2
         RETURN
a3089 4
         IF (YF.GT.3.0) THEN
            IERR = 2
            RETURN
         END IF
a3092 4
         IF (YF.LT.-3.0) THEN
            IERR = 2
            RETURN
         END IF
d3393 8
a3400 6
*     Determine the face.
      IF (XF.GT.7D0) THEN
         IERR = 2
         RETURN
      ELSE IF (XF.GT.5D0) THEN
         IF (ABS(YF).GT.1D0) THEN
d3404 7
a3413 4
         IF (ABS(YF).GT.1D0) THEN
            IERR = 2
            RETURN
         END IF
a3416 4
         IF (ABS(YF).GT.1D0) THEN
            IERR = 2
            RETURN
         END IF
a3418 3
      ELSE IF (XF.LT.-1D0) THEN
         IERR = 2
         RETURN
a3419 4
         IF (YF.GT.3D0) THEN
            IERR = 2
            RETURN
         END IF
a3422 4
         IF (YF.LT.-3D0) THEN
            IERR = 2
            RETURN
         END IF
d3697 8
a3704 6
*     Determine the face.
      IF (XF.GT.7D0) THEN
         IERR = 2
         RETURN
      ELSE IF (XF.GT.5D0) THEN
         IF (ABS(YF).GT.1D0) THEN
d3708 7
a3720 4
         IF (ABS(YF).GT.1D0) THEN
            IERR = 2
            RETURN
         END IF
a3726 4
         IF (ABS(YF).GT.1D0) THEN
            IERR = 2
            RETURN
         END IF
a3731 3
      ELSE IF (XF.LT.-1D0) THEN
         IERR = 2
         RETURN
a3732 4
         IF (YF.GT.3D0) THEN
            IERR = 2
            RETURN
         END IF
a3738 4
         IF (YF.LT.-3D0) THEN
            IERR = 2
            RETURN
         END IF
@


2.8
log
@Fixed overlap test for COP, the signs of R and sigma must agree.
@
text
@d188 1
a188 1
*   $Id: proj.f,v 2.7 1999/04/06 07:58:16 mcalabre Exp mcalabre $
d506 1
a506 1
      DOUBLE PRECISION ASIND, ATAN2D
d1622 3
a1624 1
      DOUBLE PRECISION PHI, PRJ(0:20), S, THETA, X, Y
d1635 6
a1640 2
         IERR = 2
         RETURN
d2159 1
a2159 1
*      PRJ(13)   Y0 = r0*cot(theta1) + theta1
d2644 1
a2644 1
      DOUBLE PRECISION PHI, PRJ(0:20), S, THETA, U, X, XP, Y, YP, Z
d2646 2
d2657 6
a2662 2
         IERR = 2
         RETURN
d2664 1
d2667 7
a2673 3
      IF (S.LT.-1D0 .OR. S.GT.1D0) THEN
         IERR = 2
         RETURN
@


2.7
log
@Account for missing factor of 2 in the Airy projection.
@
text
@d4 1
a4 1
*   Copyright (C) 1995-1998, Mark Calabretta
d188 1
a188 1
*   $Id: proj.f,v 2.6 1998/02/25 07:02:23 mcalabre Exp mcalabre $
d1730 1
a1730 1
      IF (PRJ(11).EQ.137D0 .AND. R.LT.0D0) THEN
@


2.6
log
@Updated date in copyright notice.
@
text
@d188 1
a188 1
*   $Id: proj.f,v 2.5 1997/03/19 05:34:13 mcalabre Exp mcalabre $
d1115 5
a1119 4
*      PRJ(12)   ln(cos(xi_b))/tan(xi_b)**2, where xi_b = (90-theta_b)/2
*      PRJ(13)   1/2 - PRJ(12)
*      PRJ(14)   r0*PRJ(13)
*      PRJ(15)   tol, cutoff for using small angle approximation, in
d1121 2
a1122 2
*      PRJ(16)   PRJ(13)*tol
*      PRJ(17)   (180/pi)/PRJ(13)
d1139 1
d1141 2
a1142 2
         PRJ(12) = -0.5D0
         PRJ(13) = 1D0
d1145 2
a1146 2
         PRJ(12) = LOG(CXI)*(CXI*CXI)/(1D0-CXI*CXI)
         PRJ(13) = 0.5D0 - PRJ(12)
d1152 4
a1155 4
      PRJ(14) = PRJ(10)*PRJ(13)
      PRJ(15) = TOL
      PRJ(16) = PRJ(13)*TOL
      PRJ(17) = R2D/PRJ(13)
d1184 2
a1185 2
         IF (XI.LT.PRJ(15)) THEN
            R = XI*PRJ(14)
d1189 1
a1189 1
            R = -PRJ(10)*(LOG(CXI)/TXI + PRJ(12)*TXI)
d1223 1
a1223 1
      R = SQRT(X*X + Y*Y)/PRJ(10)
d1227 2
a1228 2
      ELSE IF (R.LT.PRJ(16)) THEN
         XI = R*PRJ(17)
d1236 1
a1236 1
            R2 = -(LOG(X2)/TXI + PRJ(12)*TXI)
d1256 1
a1256 1
            RT = -(LOG(CXI)/TXI + PRJ(12)*TXI)
@


2.5
log
@The projection routines for AZP, TAN, SIN, ZPN, and COP now return error 2
if (phi,theta) correspond to the overlapped (far) side of the projection.
This strict bounds checking can be switched off by setting PRJ(11) to -1
(rather than 0) when the projections are initialized.
@
text
@d4 1
a4 1
*   Copyright (C) 1995-1997 Mark Calabretta
d188 1
a188 1
*   $Id: proj.f,v 2.4 1996/09/10 06:46:15 mcalabre Exp mcalabre $
@


2.3
log
@Added additional notes on usage of the PRJ array.
@
text
@d4 1
a4 1
*   Copyright (C) 1995,1996 Mark Calabretta
d137 3
a139 1
*         recompute intermediaries.
d161 10
a170 6
*      Although many of the forward projections will accept values of
*      theta outside the range [-90,90] such latitudes are not
*      meaningful and should normally be marked as an error.  However,
*      in the interests of efficiency, the forward projection routines
*      do not check for this, although they do check for any invalid
*      values of theta within the range [-90,90].
a186 4
*   Notwithstanding this, absolutely no claim is made for the accuracy
*   or reliability of these routines.  They are supplied as is, with no
*   warranty of fitness for any purpose.
*
d188 1
a188 1
*   $Id: proj.f,v 2.2 1995/12/22 02:40:14 mcalabre Exp mcalabre $
d197 5
a201 2
*      PRJ(11)   r0*(mu+1)
*      PRJ(12)   1/PRJ(11)
d214 2
a215 2
      PRJ(11) = PRJ(10)*(PRJ(1) + 1D0)
      IF (PRJ(11).EQ.0D0) THEN
d219 7
d227 5
a231 1
      PRJ(12) = 1D0/PRJ(11)
d241 1
a241 1
      DOUBLE PRECISION PHI, PRJ(0:20), R, S, THETA, X, Y
d245 1
a245 1
      IF (PRJ(11).EQ.0D0) THEN
d250 2
a251 1
      S = PRJ(1) + SIND(THETA)
d257 1
a257 1
      R =  PRJ(11)*COSD(THETA)/S
d261 5
d280 1
a280 1
      IF (PRJ(11).EQ.0D0) THEN
d292 1
a292 1
      RHO = R*PRJ(12)
d313 1
d326 5
a330 1
      PRJ(11) = 1D0
d344 1
a344 1
      IF (PRJ(11).EQ.0D0) THEN
d359 5
d376 1
a376 1
      IF (PRJ(11).EQ.0D0) THEN
d401 6
a406 5
*      PRJ(11)   1/r0
*      PRJ(12)   alpha**2 + beta**2
*      PRJ(13)   2*(alpha**2 + beta**2)
*      PRJ(14)   2*(alpha**2 + beta**2 + 1)
*      PRJ(15)   alpha**2 + beta**2 - 1
d419 11
a429 5
      PRJ(11) = 1D0/PRJ(10)
      PRJ(12) = PRJ(1)*PRJ(1) + PRJ(2)*PRJ(2)
      PRJ(13) = 2D0*PRJ(12)
      PRJ(14) = PRJ(13) + 2D0
      PRJ(15) = PRJ(12) - 1D0
d439 2
a440 1
      DOUBLE PRECISION CTHE, PHI, PRJ(0:20), T, THETA, X, Y, Z
d446 1
a446 1
      DOUBLE PRECISION COSD, SIND
d448 1
a448 1
      IF (PRJ(11).EQ.0D0) THEN
d458 1
a458 1
            Z = 2D0 - T*T/2D0
d466 4
a469 2
      X =  PRJ(10)*(CTHE*SIND(PHI) + PRJ(1)*Z)
      Y = -PRJ(10)*(CTHE*COSD(PHI) + PRJ(2)*Z)
d471 18
d508 1
a508 1
      IF (PRJ(11).EQ.0D0) THEN
d514 2
a515 2
      X0 = X*PRJ(11)
      Y0 = Y*PRJ(11)
d518 2
a519 1
      IF (PRJ(12).EQ.0D0) THEN
d528 2
d531 2
a532 1
            THETA = ASIND(SQRT(1D0 - R2))
d536 1
d545 3
a547 3
            A = PRJ(14)
            B = -(SXY + PRJ(13))
            C = R2 + SXY + PRJ(15)
a578 1
      END IF
d580 7
a586 7
*     Compute native coordinates.
      XP = -Y0 - PRJ(2)*Z
      YP =  X0 - PRJ(1)*Z
      IF (XP.EQ.0D0 .AND. YP.EQ.0D0) THEN
         PHI = 0D0
      ELSE
         PHI = ATAN2D(YP,XP)
d598 3
a600 2
*      PRJ(11)   2*r0
*      PRJ(12)   1/(2*r0)
d613 2
a614 2
         PRJ(11) = 360D0/PI
         PRJ(12) = PI/360D0
d616 2
a617 2
         PRJ(11) = 2D0*PRJ(10)
         PRJ(12) = 1D0/PRJ(11)
d620 2
d634 1
a634 1
      IF (PRJ(11).EQ.0D0) THEN
d645 1
a645 1
      R =  PRJ(11)*COSD(THETA)/S
d661 1
a661 1
      IF (PRJ(11).EQ.0D0) THEN
d672 1
a672 1
      THETA = 90D0 - 2D0*ATAND(R*PRJ(12))
d683 3
a685 2
*      PRJ(11)   r0*(pi/180)
*      PRJ(12)   (180/pi)/r0
a697 1
         PRJ(11) = 1D0
d699 1
d701 2
a702 2
         PRJ(11) = PRJ(10)*D2R
         PRJ(12) = 1D0/PRJ(11)
d705 2
d719 1
a719 1
      IF (PRJ(11).EQ.0D0) THEN
d724 1
a724 1
      R =  PRJ(11)*(90D0 - THETA)
d740 1
a740 1
      IF (PRJ(11).EQ.0D0) THEN
d751 1
a751 1
      THETA = 90D0 - R*PRJ(12)
d765 4
a768 3
*      PRJ(11)   Degree of the polynomial, N.
*      PRJ(12)   Co-latitude of the first point of inflection (N > 2).
*      PRJ(13)   Radius of the first point of inflection (N > 2).
d790 1
a790 1
 20   PRJ(11) = DBLE(K)
d842 8
a849 2
         PRJ(12) = ZD
         PRJ(13) = R
d868 1
a868 1
      IF (PRJ(11).EQ.0D0) THEN
d874 1
d884 5
d908 1
a908 1
      IF (PRJ(11).EQ.0D0) THEN
d913 1
a913 1
      K = NINT(PRJ(11))
d959 2
a960 2
         ZD2 = PRJ(12)
         R2  = PRJ(13)
d1024 3
a1026 2
*      PRJ(11)   2*r0
*      PRJ(12)   1/(2*r0)
d1039 2
a1040 2
         PRJ(11) = 360D0/PI
         PRJ(12) = PI/360D0
d1042 2
a1043 2
         PRJ(11) = 2D0*PRJ(10)
         PRJ(12) = 1D0/PRJ(11)
d1046 2
d1060 1
a1060 1
      IF (PRJ(11).EQ.0D0) THEN
d1065 1
a1065 1
      R =  PRJ(11)*SIND((90D0 - THETA)/2D0)
d1083 1
a1083 1
      IF (PRJ(11).EQ.0D0) THEN
d1095 1
a1095 1
      IF (ABS(R-PRJ(11)).LT.TOL) THEN
d1098 1
a1098 1
         THETA = 90D0 - 2D0*ASIND(R*PRJ(12))
d1114 1
a1114 1
*      PRJ(11)   State flag for PRJ(12).
a1149 1
      PRJ(11) = 1D0
d1155 2
d1173 1
a1173 1
      IF (PRJ(11).EQ.0D0) THEN
d1216 1
a1216 1
      IF (PRJ(11).EQ.0D0) THEN
d1294 5
a1298 4
*      PRJ(11)   r0*lambda*(pi/180)
*      PRJ(12)   (180/pi)/(r0*lambda)
*      PRJ(13)   r0*(mu + lambda)
*      PRJ(14)   1/(r0*(mu + lambda))
d1312 2
a1313 2
         PRJ(11) = PRJ(2)
         IF (PRJ(11).EQ.0D0) THEN
d1318 1
a1318 1
         PRJ(12) = 1D0/PRJ(11)
d1320 2
a1321 2
         PRJ(13) = R2D*(PRJ(1) + PRJ(2))
         IF (PRJ(13).EQ.0D0) THEN
d1326 1
a1326 1
         PRJ(14) = 1D0/PRJ(13)
d1328 2
a1329 2
         PRJ(11) = PRJ(10)*PRJ(2)*D2R
         IF (PRJ(11).EQ.0D0) THEN
d1334 1
a1334 1
         PRJ(12) = 1D0/PRJ(11)
d1336 2
a1337 2
         PRJ(13) = PRJ(10)*(PRJ(1) + PRJ(2))
         IF (PRJ(13).EQ.0D0) THEN
d1342 1
a1342 1
         PRJ(14) = 1D0/PRJ(13)
d1345 2
d1359 1
a1359 1
      IF (PRJ(11).EQ.0D0) THEN
d1370 2
a1371 2
      X = PRJ(11)*PHI
      Y = PRJ(13)*SIND(THETA)/S
d1385 1
a1385 1
      IF (PRJ(11).EQ.0D0) THEN
d1390 2
a1391 2
      PHI   = X*PRJ(12)
      ETA   = Y*PRJ(14)
d1403 3
a1405 2
*      PRJ(11)   r0*(pi/180)
*      PRJ(12)   (180/pi)/r0
a1417 1
         PRJ(11) = 1D0
d1419 1
d1421 2
a1422 2
         PRJ(11) = PRJ(10)*D2R
         PRJ(12) = 1D0/PRJ(11)
d1425 2
d1437 1
a1437 1
      IF (PRJ(11).EQ.0D0) THEN
d1442 2
a1443 2
      X = PRJ(11)*PHI
      Y = PRJ(11)*THETA
d1455 1
a1455 1
      IF (PRJ(11).EQ.0D0) THEN
d1460 2
a1461 2
      PHI   = PRJ(12)*X
      THETA = PRJ(12)*Y
d1472 3
a1474 2
*      PRJ(11)   r0*(pi/180)
*      PRJ(12)   (180/pi)/r0
a1486 1
         PRJ(11) = 1D0
d1488 1
d1490 2
a1491 2
         PRJ(11) = PRJ(10)*D2R
         PRJ(12) = 1D0/PRJ(11)
d1494 2
d1508 1
a1508 1
      IF (PRJ(11).EQ.0D0) THEN
d1518 1
a1518 1
      X = PRJ(11)*PHI
d1533 1
a1533 1
      IF (PRJ(11).EQ.0D0) THEN
d1538 1
a1538 1
      PHI   = X*PRJ(12)
d1554 5
a1558 4
*      PRJ(11)   r0*(pi/180)
*      PRJ(12)   (180/pi)/r0
*      PRJ(13)   r0/lambda
*      PRJ(14)   lambda/r0
a1570 1
         PRJ(11) = 1D0
d1572 1
d1577 2
a1578 2
         PRJ(13) = PRJ(10)/PRJ(1)
         PRJ(14) = PRJ(1)/PRJ(10)
d1580 2
a1581 2
         PRJ(11) = PRJ(10)*D2R
         PRJ(12) = R2D/PRJ(10)
d1586 2
a1587 2
         PRJ(13) = PRJ(10)/PRJ(1)
         PRJ(14) = PRJ(1)/PRJ(10)
d1590 2
d1604 1
a1604 1
      IF (PRJ(11).EQ.0D0) THEN
d1609 2
a1610 2
      X = PRJ(11)*PHI
      Y = PRJ(13)*SIND(THETA)
d1624 1
a1624 1
      IF (PRJ(11).EQ.0D0) THEN
d1629 1
a1629 1
      S = Y*PRJ(14)
d1635 1
a1635 1
      PHI   = X*PRJ(12)
d1652 7
a1658 6
*      PRJ(11)   C  = sin(sigma)
*      PRJ(12)   1/C
*      PRJ(13)   Y0 = r0*cos(delta)*cot(sigma)
*      PRJ(14)   r0*cos(delta)
*      PRJ(15)   1/(r0*cos(delta))
*      PRJ(16)   cot(sigma)
d1673 2
a1674 2
      PRJ(11) = SIND(PRJ(1))
      IF (PRJ(11).EQ.0D0) THEN
d1679 1
a1679 1
      PRJ(12) = 1D0/PRJ(11)
d1681 2
a1682 2
      PRJ(14) = PRJ(10)*COSD(PRJ(2))
      IF (PRJ(14).EQ.0D0) THEN
d1687 4
a1690 2
      PRJ(15) = 1D0/PRJ(14)
      PRJ(16) = 1D0/TAND(PRJ(1))
d1692 5
a1696 1
      PRJ(13) = PRJ(14)*PRJ(16)
d1706 1
a1706 1
      DOUBLE PRECISION A, PHI, PRJ(0:20), R, THETA, X, Y
d1708 1
a1708 1
      DOUBLE PRECISION COSD, SIND, TAND
d1710 1
a1710 1
      IF (PRJ(11).EQ.0D0) THEN
d1715 9
a1723 2
      A = PRJ(11)*PHI
      R = PRJ(13) - PRJ(14)*TAND(THETA-PRJ(1))
d1726 6
a1731 1
      Y = PRJ(13) - R*COSD(A)
d1745 1
a1745 1
      IF (PRJ(11).EQ.0D0) THEN
d1750 1
a1750 1
      DY = PRJ(13) - Y
d1760 2
a1761 2
      PHI   = A*PRJ(12)
      THETA = PRJ(1) + ATAND(PRJ(16) - R*PRJ(15))
d1777 5
a1781 4
*      PRJ(11)   C = r0*sin(sigma)*sin(delta)/delta
*      PRJ(12)   1/C
*      PRJ(13)   Y0 = delta*cot(delta)*cot(sigma)
*      PRJ(14)   Y0 + sigma
d1797 1
a1797 1
         PRJ(11) = PRJ(10)*SIND(PRJ(1))*D2R
d1799 1
a1799 1
         PRJ(11) = PRJ(10)*SIND(PRJ(1))*SIND(PRJ(2))/PRJ(2)
d1802 1
a1802 1
      IF (PRJ(11).EQ.0D0) THEN
d1807 5
a1811 3
      PRJ(12) = 1D0/PRJ(11)
      PRJ(13) = PRJ(10)*COSD(PRJ(2))*COSD(PRJ(1))/PRJ(11)
      PRJ(14) = PRJ(13) + PRJ(1)
d1825 1
a1825 1
      IF (PRJ(11).EQ.0D0) THEN
d1830 2
a1831 2
      A = PRJ(11)*PHI
      R = PRJ(14) - THETA
d1834 1
a1834 1
      Y = PRJ(13) - R*COSD(A)
d1848 1
a1848 1
      IF (PRJ(11).EQ.0D0) THEN
d1853 1
a1853 1
      DY = PRJ(13) - Y
d1863 2
a1864 2
      PHI   = A*PRJ(12)
      THETA = PRJ(14) - R
d1880 10
a1889 9
*      PRJ(11)   C = (sin(theta1) + sin(theta2))/2
*      PRJ(12)   1/C
*      PRJ(13)   Y0 = chi*sqrt(psi - 2C*sind(sigma))
*      PRJ(14)   chi = r0/C
*      PRJ(15)   psi = 1 + sin(theta1)*sin(theta2)
*      PRJ(16)   2C
*      PRJ(17)   (1 + sin(theta1)*sin(theta2))*(r0/C)**2
*      PRJ(18)   C/(2*r0**2)
*      PRJ(19)   Y0 = chi*sqrt(psi + 2C)
d1907 2
a1908 2
      PRJ(11) = (SIND(THETA1) + SIND(THETA2))/2D0
      IF (PRJ(11).EQ.0D0) THEN
d1913 8
a1920 1
      PRJ(12) = 1D0/PRJ(11)
d1922 1
a1922 6
      PRJ(14) = PRJ(10)/PRJ(11)
      PRJ(15) = 1D0 + SIND(THETA1)*SIND(THETA2)
      PRJ(16) = 2D0*PRJ(11)
      PRJ(17) = PRJ(14)*PRJ(14)*PRJ(15)
      PRJ(18) = 1D0/(2D0*PRJ(10)*PRJ(14))
      PRJ(19) = prj(14)*SQRT(PRJ(15) + PRJ(16))
d1924 1
a1924 1
      PRJ(13) = PRJ(14)*SQRT(PRJ(15) - PRJ(16)*SIND(PRJ(1)))
d1938 1
a1938 1
      IF (PRJ(11).EQ.0D0) THEN
d1943 1
a1943 1
      A = PHI*PRJ(11)
d1945 1
a1945 1
         R = PRJ(19)
d1947 1
a1947 1
         R = PRJ(14)*SQRT(PRJ(15) - PRJ(16)*SIND(THETA))
d1951 1
a1951 1
      Y = PRJ(13) - R*COSD(A)
d1967 1
a1967 1
      IF (PRJ(11).EQ.0D0) THEN
d1972 1
a1972 1
      DY = PRJ(13) - Y
d1982 2
a1983 2
      PHI = A*PRJ(12)
      IF (ABS(R - PRJ(19)).LT.TOL) THEN
d1986 1
a1986 1
         W = (PRJ(17) - R*R)*PRJ(18)
d2015 2
a2016 1
*      PRJ(11)   C = ln(cos(theta2)/cos(theta1))/ln(tan(tau2)/tan(tau1))
d2019 4
a2022 4
*      PRJ(12)   1/C
*      PRJ(13)   Y0 = psi*tan((90-sigma)/2)**C
*      PRJ(14)   psi = (r0*cos(theta1)/C)/tan(tau1)**C
*      PRJ(15)   1/psi
d2044 1
a2044 1
         PRJ(11) = SIND(THETA1)
d2048 1
a2048 1
         PRJ(11) = LOG(COS2/COS1)/LOG(TAN2/TAN1)
d2050 1
a2050 1
      IF (PRJ(11).EQ.0D0) THEN
d2055 1
a2055 1
      PRJ(12) = 1D0/PRJ(11)
d2057 2
a2058 2
      PRJ(14) = PRJ(10)*(COS1/PRJ(11))/(TAN1)**PRJ(11)
      IF (PRJ(14).EQ.0D0) THEN
d2062 4
a2065 2
      PRJ(15) = 1D0/PRJ(14)
      PRJ(13) = PRJ(14)*TAND((90D0-PRJ(1))/2D0)**PRJ(11)
d2079 1
a2079 1
      IF (PRJ(11).EQ.0D0) THEN
d2084 1
a2084 1
      A = PRJ(11)*PHI
d2086 1
a2086 1
         IF (PRJ(11).LT.0D0) THEN
d2093 1
a2093 1
         R = PRJ(14)*TAND((90D0 - THETA)/2D0)**PRJ(11)
d2097 1
a2097 1
      Y = PRJ(13) - R*COSD(A)
d2111 1
a2111 1
      IF (PRJ(11).EQ.0D0) THEN
d2116 1
a2116 1
      DY = PRJ(13) - Y
d2126 1
a2126 1
      PHI = A*PRJ(12)
d2128 1
a2128 1
         IF (PRJ(11).LT.0D0) THEN
d2134 1
a2134 1
         THETA = 90D0 - 2D0*ATAND((R*PRJ(15))**PRJ(12))
d2149 1
d2173 2
d2193 1
a2193 1
      IF (PRJ(13).EQ.0D0) THEN
d2222 1
a2222 1
      IF (PRJ(13).EQ.0D0) THEN
d2254 4
a2257 3
*      PRJ(11)   r0*(pi/180)
*      PRJ(12)   1/r0
*      PRJ(13)   2*r0
a2269 1
         PRJ(11) = 1D0
d2271 2
a2272 1
         PRJ(13) = 360D0/PI
d2274 3
a2276 3
         PRJ(11) = PRJ(10)*D2R
         PRJ(12) = 1D0/PRJ(11)
         PRJ(13) = 2D0*PRJ(10)
d2279 2
d2298 1
a2298 1
      IF (PRJ(11).EQ.0D0) THEN
d2308 1
a2308 1
         X = PRJ(11)*PHI
d2331 1
a2331 1
      IF (PRJ(11).EQ.0D0) THEN
d2336 1
a2336 1
      W = ABS(Y*PRJ(12))
d2338 1
a2338 1
         PHI = X*PRJ(12)
d2353 1
a2353 1
         YMTHE = Y - PRJ(11)*THEPOS
d2373 1
a2373 1
            YMTHE = Y - PRJ(11)*THETA
d2375 1
a2375 1
            F = XX + YMTHE*(YMTHE - PRJ(13)/TANTHE)
d2409 3
a2411 2
*      PRJ(11)   r0*(pi/180)
*      PRJ(12)   (180/pi)/r0
a2423 1
         PRJ(11) = 1D0
d2425 1
d2427 2
a2428 2
         PRJ(11) = PRJ(10)*D2R
         PRJ(12) = 1D0/PRJ(11)
d2431 2
d2445 1
a2445 1
      IF (PRJ(11).EQ.0D0) THEN
d2450 2
a2451 2
      X = PRJ(11)*PHI*COSD(THETA)
      Y = PRJ(11)*THETA
d2463 1
a2463 1
      IF (PRJ(11).EQ.0D0) THEN
d2472 1
a2472 1
         PHI = X*PRJ(12)/COS(Y/PRJ(10))
d2474 1
a2474 1
      THETA = Y*PRJ(12)
d2485 5
a2489 4
*      PRJ(11)   r0*(pi/180)
*      PRJ(12)   (180/pi)/r0
*      PRJ(13)   pi*r0
*      PRJ(14)   1/(pi*r0)
a2501 1
         PRJ(11) = 1D0
d2503 3
a2505 2
         PRJ(13) = 180D0
         PRJ(14) = 1D0/PRJ(13)
d2507 4
a2510 4
         PRJ(11) = PRJ(10)*D2R
         PRJ(12) = 1D0/PRJ(11)
         PRJ(13) = PI*PRJ(10)
         PRJ(14) = 1D0/PRJ(13)
d2513 2
d2527 1
a2527 1
      IF (PRJ(11).EQ.0D0) THEN
d2533 2
a2534 2
      X = PRJ(11)*PHI*(1D0 - 4D0*S*S)
      Y = PRJ(13)*S
d2548 1
a2548 1
      IF (PRJ(11).EQ.0D0) THEN
d2553 1
a2553 1
      S = Y*PRJ(14)
d2568 1
a2568 1
         PHI = PRJ(12)*X/T
d2582 5
a2586 4
*      PRJ(11)   2*r0**2
*      PRJ(12)   1/(2*r0)**2
*      PRJ(13)   1/(4*r0)**2
*      PRJ(14)   1/(2*r0)
d2599 6
a2604 4
      PRJ(11) = 2D0*PRJ(10)*PRJ(10)
      PRJ(12) = 1D0/(2D0*PRJ(11))
      PRJ(13) = PRJ(12)/4D0
      PRJ(14) = 1D0/(2D0*PRJ(10))
d2618 1
a2618 1
      IF (PRJ(11).EQ.0D0) THEN
d2624 1
a2624 1
      W = SQRT(PRJ(11)/(1D0 + COSTHE*COSD(PHI/2D0)))
d2640 1
a2640 1
      IF (PRJ(11).EQ.0D0) THEN
d2645 1
a2645 1
      U = 1D0 - X*X*PRJ(13) - Y*Y*PRJ(12)
d2658 1
a2658 1
      YP = Z*X*PRJ(14)
d2675 6
a2680 4
*      PRJ(11)   sqrt(2)*r0
*      PRJ(12)   sqrt(2)*r0/90
*      PRJ(13)   1/(sqrt(2)*r0)
*      PRJ(14)   90/r0
d2696 7
a2702 5
      PRJ(11) = SQRT2*PRJ(10)
      PRJ(12) = PRJ(11)/90D0
      PRJ(13) = 1D0/PRJ(11)
      PRJ(14) = 90D0/PRJ(10)
      PRJ(15) = 2D0/PI
d2723 1
a2723 1
      IF (PRJ(11).EQ.0D0) THEN
d2730 1
a2730 1
         Y = SIGN(PRJ(11),THETA)
d2732 1
a2732 1
         X = PRJ(12)*PHI
d2752 2
a2753 2
         X = PRJ(12)*PHI*COS(ALPHA)
         Y = PRJ(11)*SIN(ALPHA)
d2774 1
a2774 1
      IF (PRJ(11).EQ.0D0) THEN
d2795 1
a2795 1
         PHI = PRJ(14)*X/S
d2798 1
a2798 1
      Z = Y*PRJ(13)
d2806 1
a2806 1
         Z = ASIN(Z)*PRJ(15) + Y0*S/PI
d2828 3
a2830 2
*      PRJ(11)   r0*(pi/4)
*      PRJ(12)   (4/pi)/r0
d2843 2
a2844 2
         PRJ(11) = 45D0
         PRJ(12) = 1D0/45D0
d2846 2
a2847 2
         PRJ(11) = PRJ(10)*PI/4D0
         PRJ(12) = 1D0/PRJ(11)
d2850 2
d2884 1
a2884 1
      IF (PRJ(11).EQ.0D0) THEN
d2890 2
a2891 2
      L = COSTHE*SIND(PHI)
      M = COSTHE*COSD(PHI)
d2896 1
a2896 1
      IF (M.GT.RHO) THEN
d2898 1
a2898 1
         RHO  = M
d2900 1
a2900 1
      IF (L.GT.RHO) THEN
d2902 1
a2902 1
         RHO  = L
d2904 1
a2904 1
      IF (-M.GT.RHO) THEN
d2906 1
a2906 1
         RHO  = -M
d2908 1
a2908 1
      IF (-L.GT.RHO) THEN
d2910 1
a2910 1
         RHO  = -L
d2918 2
a2919 2
         XI  =  L
         ETA = -M
d2923 1
a2923 1
         XI  =  L
d2928 1
a2928 1
         XI  = -M
d2933 1
a2933 1
         XI  = -L
d2938 1
a2938 1
         XI  =  M
d2943 2
a2944 2
         XI  =  L
         ETA =  M
d2988 2
a2989 2
      X = PRJ(11)*(X0 + XF)
      Y = PRJ(11)*(Y0 + YF)
d3036 1
a3036 1
      IF (PRJ(11).EQ.0D0) THEN
d3041 2
a3042 2
      XF = X*PRJ(12)
      YF = Y*PRJ(12)
d3044 9
a3052 1
      IF (XF.GT.5.0) THEN
d3056 4
d3063 4
d3069 3
d3073 4
d3080 4
d3113 2
a3114 2
         L =  A*N
         M = -B*N
d3116 4
d3121 1
a3121 1
         L =  A*M
a3122 4
      ELSE IF (FACE.EQ.2) THEN
         L =  1D0/SQRT(A*A + B*B + 1D0)
         M = -A*L
         N =  B*L
d3124 4
d3129 1
a3129 1
         L =  A*M
a3130 4
      ELSE IF (FACE.EQ.4) THEN
         L = -1D0/SQRT(A*A + B*B + 1D0)
         M = -A*L
         N = -B*L
d3133 2
a3134 2
         L = -A*N
         M = -B*N
d3140 1
a3140 1
         PHI = ATAN2D(L, M)
d3153 3
a3155 2
*      PRJ(11)   r0*(pi/4)
*      PRJ(12)   (4/pi)/r0
d3168 2
a3169 2
         PRJ(11) = 45D0
         PRJ(12) = 1D0/45D0
d3171 2
a3172 2
         PRJ(11) = PRJ(10)*PI/4D0
         PRJ(12) = 1D0/PRJ(11)
d3175 2
d3196 1
a3196 1
      IF (PRJ(11).EQ.0D0) THEN
d3203 1
a3203 1
         Y = SIGN(2D0*PRJ(11),THETA)
d3209 2
a3210 2
      L = COSTHE*SIND(PHI)
      M = COSTHE*COSD(PHI)
d3215 1
a3215 1
      IF (M.GT.RHO) THEN
d3217 1
a3217 1
         RHO  = M
d3219 1
a3219 1
      IF (L.GT.RHO) THEN
d3221 1
a3221 1
         RHO  = L
d3223 1
a3223 1
      IF (-M.GT.RHO) THEN
d3225 1
a3225 1
         RHO  = -M
d3227 1
a3227 1
      IF (-L.GT.RHO) THEN
d3229 1
a3229 1
         RHO  = -L
d3239 2
a3240 2
         XI  =  L
         ETA = -M
d3249 1
a3249 1
         XI  =  L
d3263 1
a3263 1
         XI  = -M
d3276 1
a3276 1
         XI  = -L
d3289 1
a3289 1
         XI  =  M
d3302 2
a3303 2
         XI  =  L
         ETA =  M
d3353 2
a3354 2
      X = PRJ(11)*(XF + X0)
      Y = PRJ(11)*(YF + Y0)
d3376 1
a3376 1
      IF (PRJ(11).EQ.0D0) THEN
d3381 2
a3382 2
      XF = X*PRJ(12)
      YF = Y*PRJ(12)
d3384 9
a3392 1
      IF (XF.GT.5D0) THEN
d3396 4
d3403 4
d3409 3
d3413 4
d3420 4
d3475 4
d3480 1
a3480 1
            IF (XF.LT.0D0) L = -L
d3482 7
d3490 3
a3492 3
            M = W
            IF (YF.GT.0D0) M = -M
            L = -M*PSI
d3494 2
a3495 2
      ELSE IF (FACE.EQ.1) THEN
         M = RHO
d3498 2
a3499 2
            IF (XF.LT.0D0) L = -L
            N = L*PSI
d3503 1
a3503 1
            L = N*PSI
d3505 2
a3506 2
      ELSE IF (FACE.EQ.2) THEN
         L =  RHO
d3516 1
a3516 1
      ELSE IF (FACE.EQ.3) THEN
d3520 2
a3521 2
            IF (XF.GT.0D0) L = -L
            N = -L*PSI
d3525 1
a3525 1
            L = -N*PSI
d3527 2
a3528 2
      ELSE IF (FACE.EQ.4) THEN
         L = -RHO
d3532 1
a3532 1
            N = M*PSI
a3533 7
            N = W
            IF (YF.LT.0D0) N = -N
            M = N*PSI
         END IF
      ELSE IF (FACE.EQ.5) THEN
         N = -RHO
         IF (DIRECT) THEN
d3535 1
a3535 1
            IF (XF.LT.0D0) L = -L
a3536 4
         ELSE
            M = W
            IF (YF.LT.0D0) M = -M
            L = M*PSI
d3543 1
a3543 1
         PHI = ATAN2D(L, M)
d3556 3
a3558 2
*      PRJ(11)   r0*(pi/4)
*      PRJ(12)   (4/pi)/r0
d3571 2
a3572 2
         PRJ(11) = 45D0
         PRJ(12) = 1D0/45D0
d3574 2
a3575 2
         PRJ(11) = PRJ(10)*PI/4D0
         PRJ(12) = 1D0/PRJ(11)
d3578 2
d3595 1
a3595 1
      IF (PRJ(11).EQ.0D0) THEN
d3601 2
a3602 2
      L = COSTHE*SIND(PHI)
      M = COSTHE*COSD(PHI)
d3607 1
a3607 1
      IF (M.GT.RHO) THEN
d3609 1
a3609 1
         RHO  = M
d3611 1
a3611 1
      IF (L.GT.RHO) THEN
d3613 1
a3613 1
         RHO  = L
d3615 1
a3615 1
      IF (-M.GT.RHO) THEN
d3617 1
a3617 1
         RHO  = -M
d3619 1
a3619 1
      IF (-L.GT.RHO) THEN
d3621 1
a3621 1
         RHO  = -L
d3629 2
a3630 2
         XF =  L/RHO
         YF = -M/RHO
d3634 1
a3634 1
         XF =  L/RHO
d3639 1
a3639 1
         XF = -M/RHO
d3644 1
a3644 1
         XF = -L/RHO
d3649 1
a3649 1
         XF =  M/RHO
d3654 2
a3655 2
         XF =  L/RHO
         YF =  M/RHO
d3675 2
a3676 2
      X = PRJ(11)*(XF + X0)
      Y = PRJ(11)*(YF + Y0)
d3685 1
a3685 1
      INTEGER   FACE, IERR
d3690 1
a3690 1
      IF (PRJ(11).EQ.0D0) THEN
d3695 2
a3696 2
      XF = X*PRJ(12)
      YF = Y*PRJ(12)
d3698 10
a3707 2
      IF (XF.GT.5D0) THEN
         FACE = 4
d3709 10
d3720 1
a3720 1
         M  = -L*XF
a3721 6
      ELSE IF (XF.GT.3D0) THEN
         FACE = 3
         XF = XF - 4D0
         M  = -1D0/SQRT(1D0 + XF*XF + YF*YF)
         L  =  M*XF
         N  = -M*YF
d3723 5
a3727 1
         FACE = 2
d3729 6
a3734 3
         L  =  1D0/SQRT(1D0 + XF*XF + YF*YF)
         M  = -L*XF
         N  =  L*YF
d3736 5
a3740 1
         FACE = 0
d3743 2
a3744 2
         L  =  N*XF
         M  = -N*YF
d3746 5
a3750 1
         FACE = 5
d3753 2
a3754 2
         L  = -N*XF
         M  = -N*YF
d3756 4
a3759 4
         FACE = 1
         M  =  1D0/SQRT(1D0 + XF*XF + YF*YF)
         L  =  M*XF
         N  =  M*YF
d3765 1
a3765 1
         PHI = ATAN2D(L, M)
@


2.2
log
@Amended the conic projections to allow for sigma < 0.  Likewise for
Bonne's projection.
@
text
@d4 1
a4 1
*   Copyright (C) 1995, Mark Calabretta
d126 4
a129 1
*         correspond to the PROJPn keywords in FITS.
d186 1
a186 1
*   $Id: proj.f,v 2.1 1995/11/09 03:14:18 mcalabre Exp mcalabre $
d689 3
@


2.1
log
@Added a special-case calculation to SINREV for simple orthographic
projections.
@
text
@d183 1
a183 1
*   $Id: proj.f,v 2.0 1995/09/11 08:35:22 mcalabre Exp mcalabre $
d1549 1
a1549 1
*      PRJ(13)   Y0 = 0
d1568 1
a1568 1
      IF (PRJ(11).LE.0D0) THEN
a1573 1
      PRJ(13) = 0D0
d1584 2
d1604 1
a1604 1
      R = PRJ(14)*(PRJ(16) - TAND(THETA-PRJ(1)))
d1606 2
a1607 2
      X =  R*SIND(A)
      Y = -R*COSD(A)
d1617 1
a1617 1
      DOUBLE PRECISION A, PHI, PRJ(0:20), R, THETA, X, Y
d1626 4
a1629 1
      R = SQRT(X*X + Y*Y)
d1633 1
a1633 1
         A = ATAN2D(X, -Y)
d1655 2
a1656 2
*      PRJ(13)   Y0 = sigma - 90 + delta*cot(delta)*cot(sigma)
*      PRJ(14)   Y0 + 90
d1677 1
a1677 1
      IF (PRJ(11).LE.0D0) THEN
d1681 1
d1683 2
a1684 3

      PRJ(14) = PRJ(1) + PRJ(10)*COSD(PRJ(2))*COSD(PRJ(1))/PRJ(11)
      PRJ(13) = PRJ(14) - 90D0
d1728 2
d1733 1
a1733 1
         A = ATAN2D(X, DY)
d1755 7
a1761 4
*      PRJ(13)   Y0 = r0*sqrt(1 + sin(theta1)*sin(theta2) - 2C)/C
*      PRJ(14)   (1 + sin(theta1)*sin(theta2))*(r0/C)**2
*      PRJ(15)   2*r0**2/C
*      PRJ(16)   C/(2*r0**2)
d1766 1
a1766 1
      DOUBLE PRECISION PRJ(0:20), THETA1, THETA2, W
d1780 1
a1780 1
      IF (PRJ(11).LE.0D0) THEN
d1787 8
a1794 7
      W = 1D0 + SIND(THETA1)*SIND(THETA2)
      PRJ(13) = PRJ(10)*SQRT(W - 2D0*PRJ(11))/PRJ(11)
      PRJ(14) = W*(PRJ(10)/PRJ(11))**2

      PRJ(15) = 2D0*PRJ(10)*PRJ(10)/PRJ(11)
      PRJ(16) = 1D0/PRJ(15)
      PRJ(17) = PRJ(10)*SQRT(W + PRJ(11)+PRJ(11))/PRJ(11)
d1815 1
a1815 1
         R = PRJ(17)
d1817 1
a1817 1
         R = SQRT(PRJ(14) - PRJ(15)*SIND(THETA))
d1842 1
a1842 1
      DY = PRJ(13)-Y
d1844 2
d1849 1
a1849 1
         A = ATAN2D(X, DY)
d1853 1
a1853 1
      IF (ABS(R-PRJ(17)).LT.TOL) THEN
d1856 1
a1856 1
         W = (PRJ(14) - R*R)*PRJ(16)
d1889 3
a1891 3
*      PRJ(13)   Y0 = 0
*      PRJ(14)   (r0*cos(theta1)/C)/tan(tau1)**C
*      PRJ(15)   C/(r0*cos(theta1))*tan(tau1)**C
d1919 1
a1919 1
      IF (PRJ(11).LE.0D0) THEN
a1924 1
      PRJ(13) = 0D0
d1932 1
d1953 3
d1959 3
a1961 1
      R = PRJ(14)*TAND((90D0 - THETA)/2D0)**PRJ(11)
d1963 2
a1964 2
      X =  R*SIND(A)
      Y = -R*COSD(A)
d1974 1
a1974 1
      DOUBLE PRECISION A, PHI, PRJ(0:20), R, THETA, X, Y
d1983 4
a1986 1
      R = SQRT(X*X + Y*Y)
d1990 1
a1990 1
         A = ATAN2D(X, -Y)
d1993 10
a2002 2
      PHI   = A*PRJ(12)
      THETA = 90D0 - 2D0*ATAND((R*PRJ(15))**PRJ(12))
d2016 2
a2017 1
*      PRJ(11)   Y0 = r0*cot(theta1) + theta1
a2031 1
         PRJ(11) = PRJ(10)*COSD(PRJ(1))/SIND(PRJ(1)) + PRJ(1)
d2033 1
a2034 1
         PRJ(11) = PRJ(10)*(COSD(PRJ(1))/SIND(PRJ(1)) + PRJ(1)*D2R)
d2036 1
d2057 1
a2057 1
      IF (PRJ(11).EQ.0D0) THEN
d2062 1
a2062 1
      R = PRJ(11) - THETA*PRJ(12)
d2066 1
a2066 1
      Y = PRJ(11) - R*COSD(A)
d2086 1
a2086 1
      IF (PRJ(11).EQ.0D0) THEN
d2091 1
a2091 1
      DY = PRJ(11)-Y
d2093 2
d2098 1
a2098 1
         A = ATAN2D(X, DY)
a2099 1
      IF (PRJ(1).LT.0D0) R = -R
d2101 1
a2101 1
      THETA = (PRJ(11) - R)/PRJ(12)
@


2.0
log
@WCSLIB 2.0
@
text
@d183 1
a183 1
*   $Id: proj.f,v 1.2 1995/09/11 08:34:50 mcalabre Exp mcalabre $
d455 12
a466 4
      IF (R2.LT.1D-10) THEN
*        Use small angle formula.
         Z = -R2/2D0
         THETA = 90D0 - R2D*SQRT(R2/(1D0 - X0*PRJ(1) + Y0*PRJ(2)))
d469 4
a472 1
         SXY = 2D0*(PRJ(1)*X0 - PRJ(2)*Y0)
d474 2
a475 4
         A = PRJ(14)
         B = -(SXY + PRJ(13))
         C = R2 + SXY + PRJ(15)
         D = B*B - 2D0*A*C
d477 11
a487 6
*        Check for a solution.
         IF (D.LT.0D0) THEN
            IERR = 2
            RETURN
         END IF
         D = SQRT(D)
d489 17
a505 9
*        Choose solution closest to pole.
         STH1 = (-B + D)/A
         STH2 = (-B - D)/A
         STH = MAX(STH1,STH2)
         IF (STH.GT.1D0) THEN
            IF (STH-1D0.LT.TOL) THEN
               STH = 1D0
            ELSE
               STH = MIN(STH1,STH2)
d507 3
a509 3
         END IF
         IF (STH.LT.-1D0) THEN
            IF (STH+1D0.GT.-TOL) STH = -1D0
a510 7
         IF (STH.GT.1D0 .OR. STH.LT.-1D0) THEN
            IERR = 2
            RETURN
         END IF

         THETA = ASIND(STH)
         Z = STH - 1D0
@


1.2
log
@1) Improve precision by using small angle formulae in SIN, AIR, and QSC
   and avoid floating underflows in CSC.
2) Revise the use of p[1] and p[2] for the conics, these now contain the
   half-sum and half-differences of the standard latitudes rather than the
   standard latitudes themselves.
@
text
@d183 1
a183 1
*   $Id: proj.f,v 1.1 1995/01/31 03:20:39 mcalabre Exp mcalabre $
@


1.1
log
@Initial revision
@
text
@d36 4
a39 3
*   forward (*FWD) and reverse (*REV) transformation.  A separate
*   routine (*SET) initializes intermediate values formed from the
*   projection parameters and need not be called explicitly - see the
d68 16
a83 2
*   Forward transformation
*   ----------------------
d101 2
a102 2
*   Reverse transformation
*   ----------------------
d120 2
a121 2
*   Projection parameters; PRJ
*   --------------------------
d183 1
a183 1
*   $Id$
d397 1
a397 1
      DOUBLE PRECISION PHI, PRJ(0:20), THETA, X, Y, Z
d399 4
d410 15
a424 4
*     Compute (X,Y) from (PHI,THETA).
      Z =  SIND(THETA) - 1D0
      X =  PRJ(10)*(COSD(THETA)*SIND(PHI) + PRJ(1)*Z)
      Y = -PRJ(10)*(COSD(THETA)*COSD(PHI) + PRJ(2)*Z)
d434 2
a435 2
      DOUBLE PRECISION A, B, C, D, PHI, PRJ(0:20), PXY, STH, STH1, STH2,
     *          THETA, X, X0, XP, Y, Y0, YP
d437 6
d453 14
a466 5
      PXY = 2D0*(PRJ(1)*X0 - PRJ(2)*Y0)
      A = PRJ(14)
      B = -(PXY + PRJ(13))
      C = (X0*X0 + Y0*Y0) + PXY + PRJ(15)
      D = B*B - 2D0*A*C
d468 25
a492 6
*     Check for a solution.
      IF (D.LT.0D0) THEN
         IERR = 2
         RETURN
      END IF
      D = SQRT(D)
d494 2
a495 8
*     Choose solution closest to pole.
      STH1 = (-B + D)/A
      STH2 = (-B - D)/A
      STH = MAX(STH1,STH2)
      IF (STH.GT.1D0) STH = MIN(STH1,STH2)
      IF (STH.GT.1D0 .OR. STH.LT.-1D0) THEN
         IERR = 2
         RETURN
d499 2
a500 2
      XP = -Y0 - PRJ(2)*(STH-1D0)
      YP =  X0 - PRJ(1)*(STH-1D0)
d504 1
a504 1
         PHI   = ATAN2D(YP,XP)
a505 1
      THETA = ASIND(STH)
d1009 6
d1019 1
a1019 1
      DOUBLE PRECISION CXI, PRJ(0:20)
d1025 2
d1033 1
d1037 1
d1044 4
d1057 1
a1057 1
      DOUBLE PRECISION CXI, PHI, PRJ(0:20), R, THETA, TXI, X, Y
d1073 8
a1080 3
         CXI = COSD((90D0 - THETA)/2D0)
         TXI = SQRT(1D0-CXI*CXI)/CXI
         R = -PRJ(10)*(LOG(CXI)/TXI + PRJ(12)*TXI)
d1117 2
d1149 3
a1151 3
                IF (R-RT.LT.TOL) GO TO 40
                R1 = RT
                X1 = CXI
d1153 3
a1155 3
                IF (RT-R.LT.TOL) GO TO 40
                R2 = RT
                X2 = CXI
d1526 3
a1528 2
*      PRJ(1:2)  Latitudes of the standard parallels, theta1 and theta2,
*                in degrees.
d1532 1
a1532 1
*      PRJ(11)   C  = sin((theta1 + theta2)/2)
d1535 3
a1537 4
*      PRJ(14)   (theta1 + theta2)/2
*      PRJ(15)   r0*cos((theta2 - theta1)/2)
*      PRJ(16)   1/(r0*cos((theta2 - theta1)/2))
*      PRJ(17)   cot((theta1 + theta2)/2)
d1552 1
a1552 1
      PRJ(11) = SIND((PRJ(2) + PRJ(1))/2D0)
a1559 1
      PRJ(14) = (PRJ(2) + PRJ(1))/2D0
d1561 2
a1562 2
      PRJ(15) = PRJ(10)*COSD((PRJ(2) - PRJ(1))/2D0)
      IF (PRJ(15).EQ.0D0) THEN
d1567 2
a1568 2
      PRJ(16) = 1D0/PRJ(15)
      PRJ(17) = 1D0/TAND(PRJ(14))
d1588 1
a1588 1
      R = PRJ(15)*(PRJ(17) - TAND(THETA-PRJ(14)))
d1618 1
a1618 1
      THETA = PRJ(14) + ATAND(PRJ(17) - R*PRJ(16))
d1628 3
a1630 2
*      PRJ(1:2)  Latitudes of the standard parallels, theta1 and theta2,
*                in degrees.
d1634 1
a1634 3
*      PRJ(11)   C = r0*sin(sigma)*sin(delta)/delta, where
*                    sigma = (theta1+theta2)/2
*                    delta = (theta2-theta1)/2
d1636 1
a1636 1
*      PRJ(13)   Y0 = sigma - 90 +alpha*cot(alpha)*cot(sigma)
d1642 1
a1642 1
      DOUBLE PRECISION PRJ(0:20), THEDIF, THESUM
d1652 2
a1653 5
      THESUM = (PRJ(2) + PRJ(1))/2D0
      THEDIF = (PRJ(2) - PRJ(1))/2D0

      IF (THEDIF.EQ.0D0) THEN
         PRJ(11) = PRJ(10)*SIND(THESUM)*D2R
d1655 1
a1655 1
         PRJ(11) = PRJ(10)*SIND(THESUM)*SIND(THEDIF)/THEDIF
d1664 1
a1664 1
      PRJ(14) = THESUM + PRJ(10)*COSD(THEDIF)*COSD(THESUM)/PRJ(11)
d1723 1
a1723 1
*   CEA: conic equal area projection.
d1726 3
a1728 2
*      PRJ(1:2)  Latitudes of the standard parallels, theta1 and theta2,
*                in degrees.
d1742 1
a1742 1
      DOUBLE PRECISION PRJ(0:20), W
d1752 4
a1755 1
      PRJ(11) = (SIND(PRJ(1)) + SIND(PRJ(2)))/2D0
d1763 1
a1763 1
      W = 1D0 + SIND(PRJ(1))*SIND(PRJ(2))
d1852 3
a1854 2
*      PRJ(1:2)  Latitudes of the standard parallels, theta1 and theta2,
*                in degrees.
d1869 1
a1869 1
      DOUBLE PRECISION COS1, COS2, PRJ(0:20), TAN1, TAN2
d1879 2
a1880 2
      TAN1 = TAND((90D0 - PRJ(1))/2D0)
      COS1 = COSD(PRJ(1))
d1882 5
a1886 2
      IF (PRJ(1).EQ.PRJ(2)) THEN
         PRJ(11) = SIND(PRJ(1))
d1888 2
a1889 2
         TAN2 = TAND((90D0 - PRJ(2))/2D0)
         COS2 = COSD(PRJ(2))
d2171 1
a2171 1
         FNEG = -999.0
d2663 2
a2664 2
      REAL      A, A2, A2B2, A4, B, B2, B4, CA2, CB2, TOL, X0, XF, Y0,
     *          YF
a2753 3
      A4 = A2*A2
      B4 = B2*B2
      A2B2 = A2*B2
d2757 9
d2955 6
a2960 2
      DOUBLE PRECISION CHI, COSTHE, ETA, L, M, N, PHI, PRJ(0:20), PSI,
     *          RHO, THETA, TOL, X, X0, XF, XI, Y, Y0, YF
d3006 2
d3011 5
d3021 9
d3035 8
d3048 8
d3061 8
d3074 5
d3089 1
a3089 1
         XF  = -SQRT((1D0-RHO)/(1D0-1D0/SQRT(1D0+CHI)))
d3094 1
a3094 1
         XF  =  SQRT((1D0-RHO)/(1D0-1D0/SQRT(1D0+CHI)))
d3099 1
a3099 1
         YF  = -SQRT((1D0-RHO)/(1D0-1D0/SQRT(1D0+CHI)))
d3104 1
a3104 1
         YF  =  SQRT((1D0-RHO)/(1D0-1D0/SQRT(1D0+CHI)))
d3136 2
a3137 2
      DOUBLE PRECISION CHI, L, M, N, PHI, PRJ(0:20), PSI, RHO, THETA,
     *          TOL, X, XF, Y, YF, W
d3179 1
d3184 2
a3185 1
            RHO = 1D0 - XF*XF*(1D0 - 1D0/SQRT(1D0 + CHI))
d3192 1
d3197 2
a3198 1
            RHO = 1D0 - YF*YF*(1D0 - 1D0/SQRT(1D0 + CHI))
d3202 2
a3203 2
      IF (ABS(RHO).GE.1D0) THEN
         IF (ABS(RHO).GE.1D0+TOL) THEN
d3207 4
a3210 2
         RHO = SIGN(1D0,RHO)
         W = 0D0
d3212 1
a3212 1
         W = SQRT((1D0-RHO*RHO)/CHI)
@
