/*=============================================================================
*
*   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 Lesser 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
*   Lesser General Public License for more details.
*   
*   You should have received a copy of the GNU Lesser General Public
*   License along with this library; if not, write to the Free Software
*   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*
*   Correspondence concerning WCSLIB may be directed to:
*      Internet email: mcalabre@atnf.csiro.au
*      Postal address: Dr. Mark Calabretta,
*                      Australia Telescope National Facility,
*                      P.O. Box 76,
*                      Epping, NSW, 2121,
*                      AUSTRALIA
*
*=============================================================================
*
*   C routines which implement the FITS World Coordinate System (WCS)
*   convention.
*
*   Summary of routines
*   -------------------
*   wcsfwd() and wcsrev() are high level driver routines for the WCS linear
*   transformation, spherical coordinate transformation, and spherical
*   projection routines.
*
*   Given either the celestial longitude or latitude plus an element of the
*   pixel coordinate a hybrid routine, wcsmix(), iteratively solves for the
*   unknown elements.
*
*   An initialization routine, wcsset(), computes indices from the ctype
*   array but need not be called explicitly - see the explanation of
*   wcs.flag below.
*
*
*   Initialization routine; wcsset()
*   --------------------------------
*   Initializes elements of a wcsprm data structure which holds indices into
*   the coordinate arrays.  Note that this routine need not be called directly;
*   it will be invoked by wcsfwd() and wcsrev() if the "flag" structure member
*   is anything other than a predefined magic value.
*
*   Given:
*      naxis    const int
*                        Number of image axes.
*      ctype[][9]
*               const char
*                        Coordinate axis types corresponding to the FITS
*                        CTYPEn header cards.
*
*   Returned:
*      wcs      wcsprm*  Indices for the celestial coordinates obtained
*                        by parsing the ctype[] array (see below).
*
*   Function return value:
*               int      Error status
*                           0: Success.
*                           1: Inconsistent or unrecognized coordinate axis
*                              types.
*
*
*   Forward transformation; wcsfwd()
*   --------------------------------
*   Compute the pixel coordinate for given world coordinates.
*
*   Given:
*      ctype[][9]
*               const char
*                        Coordinate axis types corresponding to the FITS
*                        CTYPEn header cards.
*
*   Given or returned:
*      wcs      wcsprm*  Indices for the celestial coordinates obtained
*                        by parsing the ctype[] array (see below).
*
*   Given:
*      world    const double[]
*                        World coordinates.  world[wcs->lng] and
*                        world[wcs->lat] are the celestial longitude and
*                        latitude, in degrees.
*
*   Given:
*      crval    const double[]
*                        Coordinate reference values corresponding to the FITS
*                        CRVALn header cards (see note 2).
*
*   Given and returned:
*      cel      celprm*  Spherical coordinate transformation parameters (usage
*                        is described in the prologue to "cel.c").
*
*   Returned:
*      phi,     double*  Longitude and latitude in the native coordinate
*      theta             system of the projection, in degrees.
*
*   Given and returned:
*      prj      prjprm*  Projection parameters (usage is described in the
*                        prologue to "proj.c").
*
*   Returned:
*      imgcrd   double[] Image coordinate.  imgcrd[wcs->lng] and
*                        imgcrd[wcs->lat] are the projected x-, and
*                        y-coordinates, in "degrees".  For quadcube
*                        projections with a CUBEFACE axis the face number is
*                        also returned in imgcrd[wcs->cubeface].
*
*   Given and returned:
*      lin      linprm*  Linear transformation parameters (usage is described
*                        in the prologue to "lin.c").
*
*   Returned:
*      pixcrd   double[] Pixel coordinate.
*
*   Function return value:
*               int      Error status
*                           0: Success.
*                           1: Invalid coordinate transformation parameters.
*                           2: Invalid projection parameters.
*                           3: Invalid world coordinate.
*                           4: Invalid linear transformation parameters.
*
*
*   Reverse transformation; wcsrev()
*   --------------------------------
*   Compute world coordinates for a given pixel coordinate.
*
*   Given:
*      ctype[][9]
*               const char
*                        Coordinate axis types corresponding to the FITS
*                        CTYPEn header cards.
*
*   Given or returned:
*      wcs      wcsprm*  Indices for the celestial coordinates obtained
*                        by parsing the ctype[] array (see below).
*
*   Given:
*      pixcrd   const double[]
*                        Pixel coordinate.
*
*   Given and returned:
*      lin      linprm*  Linear transformation parameters (usage is described
*                        in the prologue to "lin.c").
*
*   Returned:
*      imgcrd   double[] Image coordinate.  imgcrd[wcs->lng] and
*                        imgcrd[wcs->lat] are the projected x-, and
*                        y-coordinates, in "degrees".
*
*   Given and returned:
*      prj      prjprm*  Projection parameters (usage is described in the
*                        prologue to "proj.c").
*
*   Returned:
*      phi,     double*  Longitude and latitude in the native coordinate
*      theta             system of the projection, in degrees.
*
*   Given:
*      crval    const double[]
*                        Coordinate reference values corresponding to the FITS
*                        CRVALn header cards (see note 2).
*
*   Given and returned:
*      cel      celprm*  Spherical coordinate transformation parameters
*                        (usage is described in the prologue to "cel.c").
*
*   Returned:
*      world    double[] World coordinates.  world[wcs->lng] and
*                        world[wcs->lat] are the celestial longitude and
*                        latitude, in degrees.
*
*   Function return value:
*               int      Error status
*                           0: Success.
*                           1: Invalid coordinate transformation parameters.
*                           2: Invalid projection parameters.
*                           3: Invalid pixel coordinate.
*                           4: Invalid linear transformation parameters.
*
*
*   Hybrid transformation; wcsmix()
*   -------------------------------
*   Given either the celestial longitude or latitude plus an element of the
*   pixel coordinate solve for the remaining elements by iterating on the
*   unknown celestial coordinate element using wcsfwd().
*
*   Given:
*      ctype[][9]
*               const char
*                        Coordinate axis types corresponding to the FITS
*                        CTYPEn header cards.
*
*   Given or returned:
*      wcs      wcsprm*  Indices for the celestial coordinates obtained
*                        by parsing the ctype[] array (see below).
*
*   Given:
*      mixpix   const int
*                        Which element of the pixel coordinate is given.
*      mixcel   const int
*                        Which element of the celestial coordinate is
*                        given:
*                           1: Celestial longitude is given in
*                              world[wcs->lng], latitude returned in
*                              world[wcs->lat].
*                           2: Celestial latitude is given in
*                              world[wcs->lat], longitude returned in
*                              world[wcs->lng].
*      vspan[2] const double
*                        Solution interval for the celestial coordinate, in
*                        degrees.  The ordering of the two limits is
*                        irrelevant.  Longitude ranges may be specified with
*                        any convenient normalization, for example [-120,+120]
*                        is the same as [240,480], except that the solution
*                        will be returned with the same normalization, i.e.
*                        lie within the interval specified.
*      vstep    const double
*                        Step size for solution search, in degrees.  If zero,
*                        a sensible, although perhaps non-optimal default will
*                        be used.
*      viter    int
*                        If a solution is not found then the step size will be
*                        halved and the search recommenced.  viter controls
*                        how many times the step size is halved.  The allowed
*                        range is 5 - 10.
*
*   Given and returned:
*      world    double[] World coordinates.  world[wcs->lng] and
*                        world[wcs->lat] are the celestial longitude and
*                        latitude, in degrees.  Which is given and which
*                        returned depends on the value of mixcel.  All other
*                        elements are given.
*
*   Given:
*      crval    const double[]
*                        Coordinate reference values corresponding to the FITS
*                        CRVALn header cards (see note 2).
*
*   Given and returned:
*      cel      celprm*  Spherical coordinate transformation parameters
*                        (usage is described in the prologue to "cel.c").
*
*   Returned:
*      phi,     double*  Longitude and latitude in the native coordinate
*      theta             system of the projection, in degrees.
*
*   Given and returned:
*      prj      prjprm*  Projection parameters (usage is described in the
*                        prologue to "proj.c").
*
*   Returned:
*      imgcrd   double[] Image coordinate.  imgcrd[wcs->lng] and
*                        imgcrd[wcs->lat] are the projected x-, and
*                        y-coordinates, in "degrees".
*
*   Given and returned:
*      lin      linprm*  Linear transformation parameters (usage is described
*                        in the prologue to "lin.c").
*
*   Given and returned:
*      pixcrd   double[] Pixel coordinate.  The element indicated by mixpix is
*                        given and the remaining elements are returned.
*
*   Function return value:
*               int      Error status
*                           0: Success.
*                           1: Invalid coordinate transformation parameters.
*                           2: Invalid projection parameters.
*                           3: Coordinate transformation error.
*                           4: Invalid linear transformation parameters.
*                           5: No solution found in the specified interval.
*
*
*   Notes
*   -----
*    1) The CTYPEn must in be upper case and there must be 0 or 1 pair of
*       matched celestial axis types.  The ctype[][9] should be padded with
*       blanks on the right and null-terminated.
*
*    2) Elements of the crval[] array which correspond to celestial axes are
*       ignored, the reference coordinate values in cel->ref[0] and
*       cel->ref[1] are the ones used.
*
*    3) These functions recognize the NCP projection and convert it to the
*       equivalent SIN projection.
*
*       They also recognize GLS as a synonym for SFL.
*
*    4) The quadcube projections (TSC, CSC, QSC) may be represented in FITS in
*       either of two ways:
*
*          a) The six faces may be laid out in one plane and numbered as
*             follows:
*
*                                       0
*
*                              4  3  2  1  4  3  2
*
*                                       5
*
*             Faces 2, 3 and 4 may appear on one side or the other (or both).
*             The forward routines map faces 2, 3 and 4 to the left but the
*             inverse routines accept them on either side.
*
*          b) The "COBE" convention in which the six faces are stored in a
*             three-dimensional structure using a "CUBEFACE" axis indexed from
*             0 to 5 as above.
*
*       These routines support both methods; wcsset() determines which is
*       being used by the presence or absence of a CUBEFACE axis in ctype[].
*       wcsfwd() and wcsrev() translate the CUBEFACE axis representation to
*       the single plane representation understood by the lower-level WCSLIB
*       projection routines.
*
*
*   WCS indexing parameters
*   -----------------------
*   The wcsprm struct consists of the following:
*
*      int flag
*         The wcsprm struct contains indexes and other information derived
*         from the CTYPEn.  Whenever any of the ctype[] are set or changed
*         this flag must be set to zero to signal the initialization routine,
*         wcsset() to redetermine the indices.  The flag is set to 999 if
*         there is no celestial axis pair in the CTYPEn.
*
*      char pcode[4]
*         The WCS projection code.
*
*      char lngtyp[5], lattyp[5]
*         WCS celestial axis types.
*
*      int lng,lat
*         Indices into the imgcrd[], and world[] arrays as described above.
*         These may also serve as indices for the celestial longitude and
*         latitude axes in the pixcrd[] array provided that the PC matrix
*         does not transpose axes.
*
*      int cubeface
*         Index into the pixcrd[] array for the CUBEFACE axis.  This is
*         optionally used for the quadcube projections where each cube face is
*         stored on a separate axis.
*
*
*   wcsmix() algorithm
*   ------------------
*      Initially the specified solution interval is checked to see if it's a
*      "crossing" interval.  If it isn't, a search is made for a crossing
*      solution by iterating on the unknown celestial coordinate starting at
*      the upper limit of the solution interval and decrementing by the
*      specified step size.  A crossing is indicated if the trial value of the
*      pixel coordinate steps through the value specified.  If a crossing
*      interval is found then the solution is determined by a modified form of
*      "regula falsi" division of the crossing interval.  If no crossing
*      interval was found within the specified solution interval then a search
*      is made for a "non-crossing" solution as may arise from a point of
*      tangency.  The process is complicated by having to make allowance for
*      the discontinuities that occur in all map projections.
*
*      Once one solution has been determined others may be found by subsequent
*      invokations of wcsmix() with suitably restricted solution intervals.
*
*      Note the circumstance which arises when the solution point lies at a
*      native pole of a projection in which the pole is represented as a
*      finite curve, for example the zenithals and conics.  In such cases two
*      or more valid solutions may exist but WCSMIX only ever returns one.
*
*      Because of its generality wcsmix() is very compute-intensive.  For
*      compute-limited applications more efficient special-case solvers could
*      be written for simple projections, for example non-oblique cylindrical
*      projections.
*
*   Author: Mark Calabretta, Australia Telescope National Facility
*   $Id: wcs.c,v 2.23 2002/04/03 01:25:29 mcalabre Exp $
*===========================================================================*/

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdio.h>
#include "wcslib.h"

/* Map error number to error message for each function. */
const char *wcsset_errmsg[] = {
   0,
   "Inconsistent or unrecognized coordinate axis types"};

const char *wcsfwd_errmsg[] = {
   0,
   "Invalid coordinate transformation parameters",
   "Invalid projection parameters",
   "Invalid world coordinate",
   "Invalid linear transformation parameters"};

const char *wcsrev_errmsg[] = {
   0,
   "Invalid coordinate transformation parameters",
   "Invalid projection parameters",
   "Invalid pixel coordinate",
   "Invalid linear transformation parameters"};

const char *wcsmix_errmsg[] = {
   0,
   "Invalid coordinate transformation parameters",
   "Invalid projection parameters",
   "Coordinate transformation error",
   "Invalid linear transformation parameters",
   "No solution found in the specified interval"};

#define signb(X) ((X) < 0.0 ? 1 : 0)

int
wcsset (naxis, ctype, wcs)

const int naxis;
const char ctype[][9];
struct wcsprm *wcs;

{
   int  nalias = 2;
   char aliases [2][4] = {"NCP", "GLS"};

   int j, k;
   int *ndx = NULL;
   char requir[9];

   strcpy(wcs->pcode, "");
   strcpy(requir, "");
   wcs->lng = -1;
   wcs->lat = -1;
   wcs->cubeface = -1;

   for (j = 0; j < naxis; j++) {
      if (ctype[j][4] != '-') {
         if (strcmp(ctype[j], "CUBEFACE") == 0) {
            if (wcs->cubeface == -1) {
               wcs->cubeface = j;
            } else {
               /* Multiple CUBEFACE axes! */
               return 1;
            }
         }
         continue;
      }

      /* Got an axis qualifier, is it a recognized WCS projection? */
      for (k = 0; k < npcode; k++) {
         if (strncmp(&ctype[j][5], pcodes[k], 3) == 0) break;
      }

      if (k == npcode) {
         /* Maybe it's a projection alias. */
         for (k = 0; k < nalias; k++) {
            if (strncmp(&ctype[j][5], aliases[k], 3) == 0) break;
         }

         /* Not recognized. */
         if (k == nalias) {
            continue;
         }
      }

      /* Parse the celestial axis type. */
      if (strcmp(wcs->pcode, "") == 0) {
         sprintf(wcs->pcode, "%.3s", &ctype[j][5]);

         if (strncmp(ctype[j], "RA--", 4) == 0) {
            wcs->lng = j;
            strcpy(wcs->lngtyp, "RA");
            strcpy(wcs->lattyp, "DEC");
            ndx = &wcs->lat;
            sprintf(requir, "DEC--%s", wcs->pcode);
         } else if (strncmp(ctype[j], "DEC-", 4) == 0) {
            wcs->lat = j;
            strcpy(wcs->lngtyp, "RA");
            strcpy(wcs->lattyp, "DEC");
            ndx = &wcs->lng;
            sprintf(requir, "RA---%s", wcs->pcode);
         } else if (strncmp(&ctype[j][1], "LON", 3) == 0) {
            wcs->lng = j;
            sprintf(wcs->lngtyp, "%cLON", ctype[j][0]);
            sprintf(wcs->lattyp, "%cLAT", ctype[j][0]);
            ndx = &wcs->lat;
            sprintf(requir, "%s-%s", wcs->lattyp, wcs->pcode);
         } else if (strncmp(&ctype[j][1], "LAT", 3) == 0) {
            wcs->lat = j;
            sprintf(wcs->lngtyp, "%cLON", ctype[j][0]);
            sprintf(wcs->lattyp, "%cLAT", ctype[j][0]);
            ndx = &wcs->lng;
            sprintf(requir, "%s-%s", wcs->lngtyp, wcs->pcode);
         } else if (strncmp(&ctype[j][2], "LN", 2) == 0) {
            wcs->lng = j;
            sprintf(wcs->lngtyp, "%c%cLN", ctype[j][0], ctype[j][1]);
            sprintf(wcs->lattyp, "%c%cLT", ctype[j][0], ctype[j][1]);
            ndx = &wcs->lat;
            sprintf(requir, "%s-%s", wcs->lattyp, wcs->pcode);
         } else if (strncmp(&ctype[j][2], "LT", 2) == 0) {
            wcs->lat = j;
            sprintf(wcs->lngtyp, "%c%cLN", ctype[j][0], ctype[j][1]);
            sprintf(wcs->lattyp, "%c%cLT", ctype[j][0], ctype[j][1]);
            ndx = &wcs->lng;
            sprintf(requir, "%s-%s", wcs->lngtyp, wcs->pcode);
         } else {
            /* Unrecognized celestial type. */
            return 1;
         }
      } else {
         if (strncmp(ctype[j], requir, 8) != 0) {
            /* Inconsistent projection types. */
            return 1;
         }

	if (ndx == NULL)
	    return 1;
         *ndx = j;
         strcpy(requir, "");
      }
   }

   if (strcmp(requir, "")) {
      /* Unmatched celestial axis. */
      return 1;
   }

   /* Do simple alias translations. */
   if (strncmp(wcs->pcode, "GLS", 3) == 0) {
      strcpy(wcs->pcode, "SFL");
   }

   if (strcmp(wcs->pcode, "")) {
      wcs->flag = WCSSET;
   } else {
      /* Signal for no celestial axis pair. */
      wcs->flag = 999;
   }

   return 0;
}

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

int
wcsfwd(ctype, wcs, world, crval, cel, phi, theta, prj, imgcrd, lin, pixcrd)

const char ctype[][9];
struct wcsprm* wcs;
const double world[];
const double crval[];
struct celprm *cel;
double *phi, *theta;
struct prjprm *prj;
double imgcrd[];
struct linprm *lin;
double pixcrd[];

{
   int    err, j;
   double offset;

   /* Initialize if required. */
   if (wcs->flag != WCSSET) {
      if (wcsset(lin->naxis, ctype, wcs)) return 1;
   }

   /* Convert to relative physical coordinates. */
   for (j = 0; j < lin->naxis; j++) {
      if (j == wcs->lng) continue;
      if (j == wcs->lat) continue;
      imgcrd[j] = world[j] - crval[j];
   }

   if (wcs->flag != 999) {
      /* Compute projected coordinates. */
      if (strcmp(wcs->pcode, "NCP") == 0) {
         /* Convert NCP to SIN. */
         if (cel->ref[1] == 0.0) {
            return 2;
         }

         strcpy(wcs->pcode, "SIN");
         prj->p[1] = 0.0;
         prj->p[2] = cosdeg (cel->ref[1])/sindeg (cel->ref[1]);
         prj->flag = (prj->flag < 0) ? -1 : 0;
      }

      if ((err = celfwd(wcs->pcode, world[wcs->lng], world[wcs->lat], cel,
                   phi, theta, prj, &imgcrd[wcs->lng], &imgcrd[wcs->lat]))) {
         return err;
      }

      /* Do we have a CUBEFACE axis? */
      if (wcs->cubeface != -1) {
         /* Separation between faces. */
         if (prj->r0 == 0.0) {
            offset = 90.0;
         } else {
            offset = prj->r0*PI/2.0;
         }

         /* Stack faces in a cube. */
         if (imgcrd[wcs->lat] < -0.5*offset) {
            imgcrd[wcs->lat] += offset;
            imgcrd[wcs->cubeface] = 5.0;
         } else if (imgcrd[wcs->lat] > 0.5*offset) {
            imgcrd[wcs->lat] -= offset;
            imgcrd[wcs->cubeface] = 0.0;
         } else if (imgcrd[wcs->lng] > 2.5*offset) {
            imgcrd[wcs->lng] -= 3.0*offset;
            imgcrd[wcs->cubeface] = 4.0;
         } else if (imgcrd[wcs->lng] > 1.5*offset) {
            imgcrd[wcs->lng] -= 2.0*offset;
            imgcrd[wcs->cubeface] = 3.0;
         } else if (imgcrd[wcs->lng] > 0.5*offset) {
            imgcrd[wcs->lng] -= offset;
            imgcrd[wcs->cubeface] = 2.0;
         } else {
            imgcrd[wcs->cubeface] = 1.0;
         }
      }
   }

   /* Apply forward linear transformation. */
   if (linfwd(imgcrd, lin, pixcrd)) {
      return 4;
   }

   return 0;
}

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

int
wcsrev(ctype, wcs, pixcrd, lin, imgcrd, prj, phi, theta, crval, cel, world)

const char ctype[][9];
struct wcsprm *wcs;
const double pixcrd[];
struct linprm *lin;
double imgcrd[];
struct prjprm *prj;
double *phi, *theta;
const double crval[];
struct celprm *cel;
double world[];

{
   int    err, face, j;
   double offset;

   /* Initialize if required. */
   if (wcs->flag != WCSSET) {
      if (wcsset(lin->naxis, ctype, wcs)) return 1;
   }

   /* Apply reverse linear transformation. */
   if (linrev(pixcrd, lin, imgcrd)) {
      return 4;
   }

   /* Convert to world coordinates. */
   for (j = 0; j < lin->naxis; j++) {
      if (j == wcs->lng) continue;
      if (j == wcs->lat) continue;
      world[j] = imgcrd[j] + crval[j];
   }


   if (wcs->flag != 999) {
      /* Do we have a CUBEFACE axis? */
      if (wcs->cubeface != -1) {
         face = (int)(imgcrd[wcs->cubeface] + 0.5);
         if (fabs(imgcrd[wcs->cubeface]-face) > 1e-10) {
            return 3;
         }

         /* Separation between faces. */
         if (prj->r0 == 0.0) {
            offset = 90.0;
         } else {
            offset = prj->r0*PI/2.0;
         }

         /* Lay out faces in a plane. */
         switch (face) {
         case 0:
            imgcrd[wcs->lat] += offset;
            break;
         case 1:
            break;
         case 2:
            imgcrd[wcs->lng] += offset;
            break;
         case 3:
            imgcrd[wcs->lng] += offset*2;
            break;
         case 4:
            imgcrd[wcs->lng] += offset*3;
            break;
         case 5:
            imgcrd[wcs->lat] -= offset;
            break;
         default:
            return 3;
         }
      }

      /* Compute celestial coordinates. */
      if (strcmp(wcs->pcode, "NCP") == 0) {
         /* Convert NCP to SIN. */
         if (cel->ref[1] == 0.0) {
            return 2;
         }

         strcpy(wcs->pcode, "SIN");
         prj->p[1] = 0.0;
         prj->p[2] = cosdeg (cel->ref[1])/sindeg (cel->ref[1]);
         prj->flag = (prj->flag < 0) ? -1 : 0;
      }

      if ((err = celrev(wcs->pcode, imgcrd[wcs->lng], imgcrd[wcs->lat], prj,
                   phi, theta, cel, &world[wcs->lng], &world[wcs->lat]))) {
         return err;
      }
   }

   return 0;
}

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

int
wcsmix(ctype, wcs, mixpix, mixcel, vspan, vstep, viter, world, crval, cel,
           phi, theta, prj, imgcrd, lin, pixcrd)

const char ctype[][9];
struct wcsprm *wcs;
const int mixpix, mixcel;
const double vspan[2], vstep;
int viter;
double world[];
const double crval[];
struct celprm *cel;
double *phi, *theta;
struct prjprm *prj;
double imgcrd[];
struct linprm *lin;
double pixcrd[];

{
   const int niter = 60;
   int    crossed, err, istep, iter, j, k, nstep, retry;
   const double tol  = 1.0e-10;
   const double tol2 = 100.0*tol;
   double lambda, span[2], step;
   double pixmix;
   double dlng, lng, lng0, lng0m, lng1, lng1m;
   double dlat, lat, lat0, lat0m, lat1, lat1m;
   double d, d0, d0m, d1, d1m;
   double dx = 0.0;
   double dabs, dmin, lmin;
   double dphi, phi0, phi1;
   struct celprm cel0;

   /* Initialize if required. */
   if (wcs->flag != WCSSET) {
      if (wcsset(lin->naxis, ctype, wcs)) return 1;
   }

   /* Check vspan. */
   if (vspan[0] <= vspan[1]) {
      span[0] = vspan[0];
      span[1] = vspan[1];
   } else {
      /* Swap them. */
      span[0] = vspan[1];
      span[1] = vspan[0];
   }

   /* Check vstep. */
   step = fabs(vstep);
   if (step == 0.0) {
      step = (span[1] - span[0])/10.0;
      if (step > 1.0 || step == 0.0) step = 1.0;
   }

   /* Check viter. */
   nstep = viter;
   if (nstep < 5) {
      nstep = 5;
   } else if (nstep > 10) {
      nstep = 10;
   }

   /* Given pixel element. */
   pixmix = pixcrd[mixpix];

   /* Iterate on the step size. */
   for (istep = 0; istep <= nstep; istep++) {
      if (istep) step /= 2.0;

      /* Iterate on the sky coordinate between the specified range. */
      if (mixcel == 1) {
         /* Celestial longitude is given. */

         /* Check whether the solution interval is a crossing interval. */
         lat0 = span[0];
         world[wcs->lat] = lat0;
         if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, prj,
                          imgcrd, lin, pixcrd))) {
            return err;
         }
         d0 = pixcrd[mixpix] - pixmix;

         dabs = fabs(d0);
         if (dabs < tol) return 0;

         lat1 = span[1];
         world[wcs->lat] = lat1;
         if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, prj,
                          imgcrd, lin, pixcrd))) {
            return err;
         }
         d1 = pixcrd[mixpix] - pixmix;

         dabs = fabs(d1);
         if (dabs < tol) return 0;

         lmin = lat1;
         dmin = dabs;

         /* Check for a crossing point. */
         if (signb(d0) != signb(d1)) {
            crossed = 1;
            dx = d1;
         } else {
            crossed = 0;
            lat0 = span[1];
         }

         for (retry = 0; retry < 4; retry++) {
            /* Refine the solution interval. */
            while (lat0 > span[0]) {
               lat0 -= step;
               if (lat0 < span[0]) lat0 = span[0];
               world[wcs->lat] = lat0;
               if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                                prj, imgcrd, lin, pixcrd))) {
                  return err;
               }
               d0 = pixcrd[mixpix] - pixmix;

               /* Check for a solution. */
               dabs = fabs(d0);
               if (dabs < tol) return 0;

               /* Record the point of closest approach. */
               if (dabs < dmin) {
                  lmin = lat0;
                  dmin = dabs;
               }

               /* Check for a crossing point. */
               if (signb(d0) != signb(d1)) {
                  crossed = 2;
                  dx = d0;
                  break;
               }

               /* Advance to the next subinterval. */
               lat1 = lat0;
               d1 = d0;
            }

            if (crossed) {
               /* A crossing point was found. */
               for (iter = 0; iter < niter; iter++) {
                  /* Use regula falsi division of the interval. */
                  lambda = d0/(d0-d1);
                  if (lambda < 0.1) {
                     lambda = 0.1;
                  } else if (lambda > 0.9) {
                     lambda = 0.9;
                  }

                  dlat = lat1 - lat0;
                  lat = lat0 + lambda*dlat;
                  world[wcs->lat] = lat;
                  if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                                   prj, imgcrd, lin, pixcrd))) {
                     return err;
                  }

                  /* Check for a solution. */
                  d = pixcrd[mixpix] - pixmix;
                  dabs = fabs(d);
                  if (dabs < tol) return 0;

                  if (dlat < tol) {
                     /* An artifact of numerical imprecision. */
                     if (dabs < tol2) return 0;

                     /* Must be a discontinuity. */
                     break;
                  }

                  /* Record the point of closest approach. */
                  if (dabs < dmin) {
                     lmin = lat;
                     dmin = dabs;
                  }

                  if (signb(d0) == signb(d)) {
                     lat0 = lat;
                     d0 = d;
                  } else {
                     lat1 = lat;
                     d1 = d;
                  }
               }

               /* No convergence, must have been a discontinuity. */
               if (crossed == 1) lat0 = span[1];
               lat1 = lat0;
               d1 = dx;
               crossed = 0;

            } else {
               /* No crossing point; look for a tangent point. */
               if (lmin == span[0]) break;
               if (lmin == span[1]) break;

               lat = lmin;
               lat0 = lat - step;
               if (lat0 < span[0]) lat0 = span[0];
               lat1 = lat + step;
               if (lat1 > span[1]) lat1 = span[1];

               world[wcs->lat] = lat0;
               if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                                prj, imgcrd, lin, pixcrd))) {
                  return err;
               }
               d0 = fabs(pixcrd[mixpix] - pixmix);

               d  = dmin;

               world[wcs->lat] = lat1;
               if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                                prj, imgcrd, lin, pixcrd))) {
                  return err;
               }
               d1 = fabs(pixcrd[mixpix] - pixmix);

               for (iter = 0; iter < niter; iter++) {
                  lat0m = (lat0 + lat)/2.0;
                  world[wcs->lat] = lat0m;
                  if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                                   prj, imgcrd, lin, pixcrd))) {
                     return err;
                  }
                  d0m = fabs(pixcrd[mixpix] - pixmix);

                  if (d0m < tol) return 0;

                  lat1m = (lat1 + lat)/2.0;
                  world[wcs->lat] = lat1m;
                  if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                                   prj, imgcrd, lin, pixcrd))) {
                     return err;
                  }
                  d1m = fabs(pixcrd[mixpix] - pixmix);

                  if (d1m < tol) return 0;

                  if (d0m < d && d0m <= d1m) {
                     lat1 = lat;
                     d1   = d;
                     lat  = lat0m;
                     d    = d0m;
                  } else if (d1m < d) {
                     lat0 = lat;
                     d0   = d;
                     lat  = lat1m;
                     d    = d1m;
                  } else {
                     lat0 = lat0m;
                     d0   = d0m;
                     lat1 = lat1m;
                     d1   = d1m;
                  }
               }
            }
         }

      } else {
         /* Celestial latitude is given. */

         /* Check whether the solution interval is a crossing interval. */
         lng0 = span[0];
         world[wcs->lng] = lng0;
         if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, prj,
                          imgcrd, lin, pixcrd))) {
            return err;
         }
         d0 = pixcrd[mixpix] - pixmix;

         dabs = fabs(d0);
         if (dabs < tol) return 0;

         lng1 = span[1];
         world[wcs->lng] = lng1;
         if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta, prj,
                          imgcrd, lin, pixcrd))) {
            return err;
         }
         d1 = pixcrd[mixpix] - pixmix;

         dabs = fabs(d1);
         if (dabs < tol) return 0;
         lmin = lng1;
         dmin = dabs;

         /* Check for a crossing point. */
         if (signb(d0) != signb(d1)) {
            crossed = 1;
            dx = d1;
         } else {
            crossed = 0;
            lng0 = span[1];
         }

         for (retry = 0; retry < 4; retry++) {
            /* Refine the solution interval. */
            while (lng0 > span[0]) {
               lng0 -= step;
               if (lng0 < span[0]) lng0 = span[0];
               world[wcs->lng] = lng0;
               if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                          prj, imgcrd, lin, pixcrd))) {
                  return err;
               }
               d0 = pixcrd[mixpix] - pixmix;

               /* Check for a solution. */
               dabs = fabs(d0);
               if (dabs < tol) return 0;

               /* Record the point of closest approach. */
               if (dabs < dmin) {
                  lmin = lng0;
                  dmin = dabs;
               }

               /* Check for a crossing point. */
               if (signb(d0) != signb(d1)) {
                  crossed = 2;
                  dx = d0;
                  break;
               }

               /* Advance to the next subinterval. */
               lng1 = lng0;
               d1 = d0;
            }

            if (crossed) {
               /* A crossing point was found. */
               for (iter = 0; iter < niter; iter++) {
                  /* Use regula falsi division of the interval. */
                  lambda = d0/(d0-d1);
                  if (lambda < 0.1) {
                     lambda = 0.1;
                  } else if (lambda > 0.9) {
                     lambda = 0.9;
                  }

                  dlng = lng1 - lng0;
                  lng = lng0 + lambda*dlng;
                  world[wcs->lng] = lng;
                  if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                                   prj, imgcrd, lin, pixcrd))) {
                     return err;
                  }

                  /* Check for a solution. */
                  d = pixcrd[mixpix] - pixmix;
                  dabs = fabs(d);
                  if (dabs < tol) return 0;

                  if (dlng < tol) {
                     /* An artifact of numerical imprecision. */
                     if (dabs < tol2) return 0;

                     /* Must be a discontinuity. */
                     break;
                  }

                  /* Record the point of closest approach. */
                  if (dabs < dmin) {
                     lmin = lng;
                     dmin = dabs;
                  }

                  if (signb(d0) == signb(d)) {
                     lng0 = lng;
                     d0 = d;
                  } else {
                     lng1 = lng;
                     d1 = d;
                  }
               }

               /* No convergence, must have been a discontinuity. */
               if (crossed == 1) lng0 = span[1];
               lng1 = lng0;
               d1 = dx;
               crossed = 0;

            } else {
               /* No crossing point; look for a tangent point. */
               if (lmin == span[0]) break;
               if (lmin == span[1]) break;

               lng = lmin;
               lng0 = lng - step;
               if (lng0 < span[0]) lng0 = span[0];
               lng1 = lng + step;
               if (lng1 > span[1]) lng1 = span[1];

               world[wcs->lng] = lng0;
               if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                                prj, imgcrd, lin, pixcrd))) {
                  return err;
               }
               d0 = fabs(pixcrd[mixpix] - pixmix);

               d  = dmin;

               world[wcs->lng] = lng1;
               if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                                prj, imgcrd, lin, pixcrd))) {
                  return err;
               }
               d1 = fabs(pixcrd[mixpix] - pixmix);

               for (iter = 0; iter < niter; iter++) {
                  lng0m = (lng0 + lng)/2.0;
                  world[wcs->lng] = lng0m;
                  if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                                   prj, imgcrd, lin, pixcrd))) {
                     return err;
                  }
                  d0m = fabs(pixcrd[mixpix] - pixmix);

                  if (d0m < tol) return 0;

                  lng1m = (lng1 + lng)/2.0;
                  world[wcs->lng] = lng1m;
                  if ((err = wcsfwd(ctype, wcs, world, crval, cel, phi, theta,
                                   prj, imgcrd, lin, pixcrd))) {
                     return err;
                  }
                  d1m = fabs(pixcrd[mixpix] - pixmix);

                  if (d1m < tol) return 0;

                  if (d0m < d && d0m <= d1m) {
                     lng1 = lng;
                     d1   = d;
                     lng  = lng0m;
                     d    = d0m;
                  } else if (d1m < d) {
                     lng0 = lng;
                     d0   = d;
                     lng  = lng1m;
                     d    = d1m;
                  } else {
                     lng0 = lng0m;
                     d0   = d0m;
                     lng1 = lng1m;
                     d1   = d1m;
                  }
               }
            }
         }
      }
   }


   /* Set cel0 to the unity transformation. */
   cel0.flag = CELSET;
   cel0.ref[0] = cel->ref[0];
   cel0.ref[1] = cel->ref[1];
   cel0.ref[2] = cel->ref[2];
   cel0.ref[3] = cel->ref[3];
   cel0.euler[0] = -90.0;
   cel0.euler[1] =   0.0;
   cel0.euler[2] =  90.0;
   cel0.euler[3] =   1.0;
   cel0.euler[4] =   0.0;

   /* No convergence, check for aberrant behaviour at a native pole. */
   *theta = -90.0;
   for (j = 1; j <= 2; j++) {
      /* Could the celestial coordinate element map to a native pole? */
      *theta = -*theta;
      err = sphrev(0.0, *theta, cel->euler, &lng, &lat);

      if (mixcel == 1) {
         if (fabs(fmod(world[wcs->lng]-lng,360.0)) > tol) continue;
         if (lat < span[0]) continue;
         if (lat > span[1]) continue;
         world[wcs->lat] = lat;
      } else {
         if (fabs(world[wcs->lat]-lat) > tol) continue;
         if (lng < span[0]) lng += 360.0;
         if (lng > span[1]) lng -= 360.0;
         if (lng < span[0]) continue;
         if (lng > span[1]) continue;
         world[wcs->lng] = lng;
      }

      /* Is there a solution for the given pixel coordinate element? */
      lng = world[wcs->lng];
      lat = world[wcs->lat];

      /* Feed native coordinates to wcsfwd() with cel0 set to unity. */
      world[wcs->lng] = -180.0;
      world[wcs->lat] = *theta;
      if ((err = wcsfwd(ctype, wcs, world, crval, &cel0, phi, theta, prj,
                       imgcrd, lin, pixcrd))) {
         return err;
      }
      d0 = pixcrd[mixpix] - pixmix;

      /* Check for a solution. */
      if (fabs(d0) < tol) {
         /* Recall saved world coordinates. */
         world[wcs->lng] = lng;
         world[wcs->lat] = lat;
         return 0;
      }

      /* Search for a crossing interval. */
      phi0 = -180.0;
      for (k = -179; k <= 180; k++) {
         phi1 = (double) k;
         world[wcs->lng] = phi1;
         if ((err = wcsfwd(ctype, wcs, world, crval, &cel0, phi, theta, prj,
                          imgcrd, lin, pixcrd))) {
            return err;
         }
         d1 = pixcrd[mixpix] - pixmix;

         /* Check for a solution. */
         dabs = fabs(d1);
         if (dabs < tol) {
            /* Recall saved world coordinates. */
            world[wcs->lng] = lng;
            world[wcs->lat] = lat;
            return 0;
         }

         /* Is it a crossing interval? */
         if (signb(d0) != signb(d1)) break;

         phi0 = phi1;
         d0 = d1;
      }

      for (iter = 1; iter <= niter; iter++) {
         /* Use regula falsi division of the interval. */
         lambda = d0/(d0-d1);
         if (lambda < 0.1) {
            lambda = 0.1;
         } else if (lambda > 0.9) {
            lambda = 0.9;
         }

         dphi = phi1 - phi0;
         world[wcs->lng] = phi0 + lambda*dphi;
         if ((err = wcsfwd(ctype, wcs, world, crval, &cel0, phi, theta, prj,
                          imgcrd, lin, pixcrd))) {
            return err;
         }

         /* Check for a solution. */
         d = pixcrd[mixpix] - pixmix;
         dabs = fabs(d);
         if (dabs < tol || (dphi < tol && dabs < tol2)) {
            /* Recall saved world coordinates. */
            world[wcs->lng] = lng;
            world[wcs->lat] = lat;
            return 0;
         }

         if (signb(d0) == signb(d)) {
            phi0 = world[wcs->lng];
            d0 = d;
         } else {
            phi1 = world[wcs->lng];
            d1 = d;
         }
      }
   }


   /* No solution. */
   return 5;

}
/* Dec 20 1999  Doug Mink - Change signbit() to signb() and always define it
 * Dec 20 1999	Doug Mink - Include wcslib.h, which includes wcs.h, wcstrig.h
 *
 * Mar 20 2001	Doug Mink - Include stdio.h for sprintf()
 * Mar 20 2001	Doug Mink - Add () around err assignments in if statements
 * Sep 19 2001	Doug Mink - Add above changes to WCSLIB-2.7 version
 *
 * Mar 15 2002	Doug Mink - Add above changes to WCSLIB-2.8.2
 * Apr  3 2002	Mark Calabretta - Fix bug in code checking section
 *
 * Jun 20 2006	Doug Mink - Initialized uninitialized variables
 */