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

  PGSBOX 4.8 - draw curvilinear coordinate axes for PGPLOT.
  Copyright (C) 1997-2011, 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.8.1.1 2011/08/15 08:07:07 cal103 Exp cal103 $
*===========================================================================*/

#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;
}