/*============================================================================

    PGSBOX 4.3 - an implementation of the FITS WCS standard.
    Copyright (C) 1997-2007, Mark Calabretta

    This file is part of PGSBOX.

    PGSBOX 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 3 of the License, or (at your
    option) any later version.

    PGSBOX 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 PGSBOX.  If not, see <http://www.gnu.org/licenses/>.

    Correspondence concerning PGSBOX may be directed to:
       Internet email: mcalabre@atnf.csiro.au
       Postal address: Dr. Mark Calabretta
                       Australia Telescope National Facility, CSIRO
                       PO Box 76
                       Epping NSW 1710
                       AUSTRALIA

    Author: Mark Calabretta, Australia Telescope National Facility
    http://www.atnf.csiro.au/~mcalabre/index.html
    $Id: pgwcsl.c,v 4.3 2007/12/27 05:49:14 cal103 Exp $
*===========================================================================*/

#include <math.h>

#include <wcs.h>
#include <sph.h>

/* Fortran name mangling. */
#include <wcsconfig_f77.h>
#define pgwcsl_ F77_FUNC(pgwcsl, PGWCSL)

void pgwcsl_(opcode, nlc, nli, nld, nlcprm, wcs, nldprm, world, pixel, contrl,
            contxt, ierr)

const int *opcode, *nlc, *nli, *nld;
const char *nlcprm;
int *wcs;
double *nldprm;
double *world;
double *pixel;
int *contrl;
double contxt[20];
int *ierr;

{
   int i, outside, stat;
   double dp, imgcrd[9], lat, lng, ph, phi, sdummy, th, theta;
   static double wrld[9];
   struct wcsprm *wcsp;
   struct celprm *wcscel;

   *ierr = 0;

   wcsp = (struct wcsprm *)wcs;
   wcscel = &(wcsp->cel);

   if (*opcode == 2) {
      /* Compute pixel coordinates from world coordinates. */
      if (wcsp->lng < 0) {
         /* Simple linear coordinates. */
         if (wcss2p(wcsp, 1, 0, world, &phi, &theta, imgcrd, pixel, &stat)) {
            *ierr = 1;
         }
         return;
      }

      wrld[wcsp->lng] = world[0];
      if (world[1] > 90.0) {
         wrld[wcsp->lat] = 90.0;
         outside = -2;
      } else if (world[1] < -90.0) {
         wrld[wcsp->lat] = -90.0;
         outside = -2;
      } else {
         wrld[wcsp->lat] = world[1];
         outside = 0;
      }

      if (*contrl == 0) {
         if (wcss2p(wcsp, 1, 0, wrld, &phi, &theta, imgcrd, pixel, &stat)) {
            /* Translate status return values. */
            *ierr = stat ? 2 : 1;
            return;
         }

         if (fabs(phi-contxt[2]) > 180.0) {
            /* Hit a discontinuity at phi = +/- 180. */
            contxt[4] = pixel[0];
            contxt[5] = pixel[1];

            if (contxt[2] > phi) {
               ph = 179.9999;
               dp = (phi - contxt[2]) + 360.0;
            } else {
               ph = -179.9999;
               dp = (phi - contxt[2]) - 360.0;
            }

            /* First approximation for theta. */
            if (dp == 0.0) {
               th = contxt[3];
            } else {
               th = contxt[3] + (ph-contxt[2])*(theta-contxt[3])/dp;
            }

            /* Iterate once to refine the value of theta. */
            sphx2s(wcscel->euler, 1, 1, 1, 1, &ph, &th, &lng, &lat);
            if (wrld[wcsp->lng] == contxt[0]) {
               /* We are following a meridian of longitude. */
               lng = wrld[wcsp->lng];
            } else {
               /* We are following a parallel of latitude. */
               lat = wrld[wcsp->lat];
            }
            sphs2x(wcscel->euler, 1, 1, 1, 1, &lng, &lat, &sdummy, &th);

            contxt[0] = wrld[wcsp->lng];
            contxt[1] = wrld[wcsp->lat];
            contxt[2] = phi;
            contxt[3] = theta;

            /* Pixel coordinates crossing into the discontinuity. */
            sphx2s(wcscel->euler, 1, 1, 1, 1, &ph, &th, wrld+wcsp->lng,
               wrld+wcsp->lat);
            if (wcss2p(wcsp, 1, 0, wrld, &phi, &theta, imgcrd, pixel,
               &stat)) {
               /* Translate status return values. */
               *ierr = stat ? 2 : 1;
               return;
            }

            /* Pixel coordinates crossing out of the discontinuity. */
            ph *= -1.0;
            sphx2s(wcscel->euler, 1, 1, 1, 1, &ph, &th, wrld+wcsp->lng,
               wrld+wcsp->lat);
            if (wcss2p(wcsp, 1, 0, wrld, &phi, &theta, imgcrd, contxt+6,
               &stat)) {
               /* Translate status return values. */
               *ierr = stat ? 2 : 1;
               return;
            }

            *contrl = 1;
         } else {
            /* Normal mode, no discontinuity. */
            contxt[0] = wrld[wcsp->lng];
            contxt[1] = wrld[wcsp->lat];
            contxt[2] = phi;
            contxt[3] = theta;
         }
      } else {
         if (*contrl == 1) {
            /* Move to the other side of the discontinuity. */
            pixel[0] = contxt[6];
            pixel[1] = contxt[7];
            *contrl = 2;
         } else {
            /* Complete the traversal. */
            pixel[0] = contxt[4];
            pixel[1] = contxt[5];
            *contrl = 0;
         }
      }

      *ierr = outside;

   } else if (*opcode == 1) {
      /* Compute pixel coordinates from world coordinates. */
      if (wcsp->lng < 0) {
         /* Simple linear coordinates. */
         if (wcss2p(wcsp, 1, 0, world, &phi, &theta, imgcrd, pixel, &stat)) {
            *ierr = 1;
         }
         return;
      }

      wrld[wcsp->lng] = world[0];
      if (world[1] > 90.0) {
         wrld[wcsp->lat] = 90.0;
         outside = -2;
      } else if (world[1] < -90.0) {
         wrld[wcsp->lat] = -90.0;
         outside = -2;
      } else {
         wrld[wcsp->lat] = world[1];
         outside = 0;
      }

      if (wcss2p(wcsp, 1, 0, wrld, &phi, &theta, imgcrd, pixel, &stat)) {
         /* Translate status return values. */
         *ierr = stat ? 2 : 1;
         return;
      }

      contxt[0] = wrld[wcsp->lng];
      contxt[1] = wrld[wcsp->lat];
      contxt[2] = phi;
      contxt[3] = theta;

      *ierr = outside;

   } else if (*opcode == 0) {
      /* Initialize. */
      if (*nli < WCSLEN) {
         *ierr = 1;
         return;
      }

      if ((*ierr = wcsset(wcsp))) {
         *ierr = *ierr <= 2 ? 1 : 2;
      }

      for (i = 2; i < 9; i++) {
         wrld[i] = 0.0;
      }

      *contrl = 0;

   } else if (*opcode == -1) {
      /* Compute world coordinates from pixel coordinates. */
      if (wcsp2s(wcsp, 1, 0, pixel, imgcrd, &phi, &theta, wrld, &stat)) {
         /* Translate status return values. */
         *ierr = stat ? 3 : 1;
         return;
      }

      if (wcsp->lng < 0) {
         /* Simple linear coordinates. */
         world[0] = wrld[0];
         world[1] = wrld[1];
      } else {
         world[0] = wrld[wcsp->lng];
         world[1] = wrld[wcsp->lat];

         if (phi < -180.0 || phi > 180.0) {
            /* Pixel is outside the principle range of native longitude. */
            *ierr = 3;
            return;
         }
      }

   } else {
      *ierr = 1;
   }

   return;
}