/** @file  psAstrometry.c
 *
 *  @brief This file defines the basic types for astronomical coordinate
 *  transformation
 *
 *  @ingroup AstroImage
 *
 *  @author GLG, MHPCC
 *
 *  @version $Revision: 1.58 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-01-17 23:25:09 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 */

/******************************************************************************/
/*  INCLUDE FILES                                                             */
/******************************************************************************/
#include <string.h>
#include <math.h>

#include "psFunctions.h"
#include "psAstrometry.h"
#include "psMemory.h"
#include "psError.h"
#include "psConstants.h"
#include "psAstronomyErrors.h"
#include "psTrace.h"
#include "psLogMsg.h"

#include "slalib.h"

/*****************************************************************************
checkValidImageCoords(): this is a private function which simply
determines if the supplied x,y coordinates are in the range for the supplied
psImage.
 *****************************************************************************/
static psS32 checkValidImageCoords(double x,
                                   double y,
                                   psImage* tmpImage)
{
    PS_IMAGE_CHECK_NULL(tmpImage, 0);

    if ((x < 0.0) || (x > (double)tmpImage->numCols) ||
            (y < 0.0) || (y > (double)tmpImage->numRows)) {
        return (0);
    }

    return (1);
}

/*****************************************************************************
isProjectionLinear(): this is a private function which simply determines
if the supplied psPlaneTransform transform is linear: if any of the
cooefficients of order 2 are higher are non-zero, then it is not linear.
 *****************************************************************************/
static psS32 isProjectionLinear(psPlaneTransform *transform)
{
    PS_PTR_CHECK_NULL(transform, 0);
    PS_PTR_CHECK_NULL(transform->x, 0);
    PS_PTR_CHECK_NULL(transform->y, 0);

    int i;

    for (i=2;i<(transform->x->nX);i++) {
        if (transform->x->coeff[i][0] != 0.0) {
            return(0);
        }
    }
    for (i=2;i<(transform->x->nY);i++) {
        if (transform->x->coeff[0][i] != 0.0) {
            return(0);
        }
    }

    for (i=2;i<(transform->y->nX);i++) {
        if (transform->y->coeff[i][0] != 0.0) {
            return(0);
        }
    }
    for (i=2;i<(transform->y->nY);i++) {
        if (transform->y->coeff[0][i] != 0.0) {
            return(0);
        }
    }

    return(1);
}

/*****************************************************************************
invertPlaneTransform(transform): : this is a private function which
simply inverts the supplied psPlaneTransform transform.  It assumes that
"transform" is linear.
 
This program assumes that the inverse of the following linear equations:
        X2 = A + (B * X1) + (C * Y1);
        Y2 = D + (E * X1) + (F * Y1);
is
        Y1 = (Y2 - ((E/B) * X2) - D + ((E*A)/B)) / (F - ((C*E)/B));
        X1 = (Y2 - ((F/C) * X2) - D + ((F*A)/C)) / (E - ((F*B)/C));
or
 X1 = (-D + ((F*A)/C)) / (E - ((F*B)/C)) +
      (X2 * -((F/C) / (E - ((F*B)/C)))) +
             (Y2 * (1.0 / (E - ((F*B)/C))));
        Y1 = (-D + ((E*A)/B))/(F - ((C*E)/B)) +
             (X2 * -((E/B) / (F - ((C*E)/B)))) +
             (Y2 * (1.0 / (F - ((C*E)/B))));
 *****************************************************************************/
static psPlaneTransform *invertPlaneTransform(psPlaneTransform *transform)
{
    PS_PTR_CHECK_NULL(transform, 0);
    PS_PTR_CHECK_NULL(transform->x, 0);
    PS_PTR_CHECK_NULL(transform->y, 0);

    double A = transform->x->coeff[0][0];
    double B = transform->x->coeff[1][0];
    double C = transform->x->coeff[0][1];
    double D = transform->y->coeff[0][0];
    double E = transform->y->coeff[1][0];
    double F = transform->y->coeff[0][1];

    // XXX: Use the constructor here.
    psPlaneTransform *out = psAlloc(sizeof(psPlaneTransform));

    out->x->coeff[0][0] = -D + ((F*A)/C) / (E - ((F*B)/C));
    out->x->coeff[1][0] = -(F/C) / (E - ((F*B)/C));
    out->x->coeff[0][1] =  1.0 / (E - ((F*B)/C));
    out->y->coeff[0][0] = -D + ((E*A)/B) / (F - ((C*E)/B));
    out->y->coeff[1][0] = -(E/B) / (F - ((C*E)/B));
    out->y->coeff[0][1] =  1.0 / (F - ((C*E)/B));

    return(out);
}

static void FPAFree(psFPA* fpa)
{
    if (fpa != NULL) {
        psFree(fpa->chips);
        psFree(fpa->grommit);
        psFree((psExposure*)fpa->exposure);
        psFree(fpa->metadata);
        psFree(fpa->fromTangentPlane);
        psFree(fpa->toTangentPlane);
        psFree(fpa->pattern);
        psFree(fpa->colorPlus);
        psFree(fpa->colorMinus);
        psFree(fpa->projection);
    }
}

static void chipFree(psChip* chip)
{
    if (chip != NULL) {
        psFree(chip->cells);
        psFree(chip->metadata);
        psFree(chip->toFPA);
        psFree(chip->fromFPA);
    }
}

static void cellFree(psCell* cell)
{
    if (cell != NULL) {
        psFree(cell->readouts);
        psFree(cell->metadata);
        psFree(cell->toChip);
        psFree(cell->fromChip);
        psFree(cell->toFPA);
        psFree(cell->toTP);
        psFree(cell->toSky);
    }
}

static void readoutFree(psReadout* readout)
{
    if (readout != NULL) {
        psFree(readout->image);
        psFree(readout->mask);
        psFree(readout->objects);
        psFree(readout->metadata);
    }
}

static void observatoryFree(psObservatory* obs)
{
    if (obs != NULL) {
        psFree((psPtr)obs->name);
    }
}

static void exposureFree(psExposure* exp)
{
    if (exp != NULL) {
        psFree((psPtr)exp->time);
        psFree((psPtr)exp->observatory);
        psFree((psPtr)exp->cameraName);
        psFree((psPtr)exp->telescopeName);
    }
}

static void fixedPatternFree(psFixedPattern* fp)
{
    if (fp != NULL) {
        for (psS32 i = 0; i < fp->p_ps_xRows; i++) {
            psFree(fp->x[i]);
        }

        for (psS32 j = 0; j < fp->p_ps_yRows; j++) {
            psFree(fp->y[j]);
        }

        psFree(fp->x);
        psFree(fp->y);
    }
}

/*****************************************************************************/
/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
/*****************************************************************************/

/*
 * XXX: Verify that you interpreted the SDR correctly.
 *
 * XXX: This assumes that x,y must be of type F64
 */
psFixedPattern* psFixedPatternAlloc(double x0,
                                    double y0,
                                    double xScale,
                                    double yScale,
                                    const psImage *x,
                                    const psImage *y)
{
    psFixedPattern *tmp;
    psS32 i;
    psS32 j;

    PS_IMAGE_CHECK_NULL(x, NULL);
    PS_IMAGE_CHECK_NULL(y, NULL);
    PS_IMAGE_CHECK_TYPE(x, PS_TYPE_F64, NULL);
    PS_IMAGE_CHECK_TYPE(y, PS_TYPE_F64, NULL);

    tmp = (psFixedPattern *) psAlloc(sizeof(psFixedPattern));
    // XXX: Is this correct?
    tmp->nX = (x->numCols * x->numRows);
    tmp->nY = (y->numCols * y->numRows);
    tmp->x0 = x0;
    tmp->y0 = y0;
    tmp->xScale = xScale;
    tmp->yScale = yScale;
    tmp->p_ps_xRows = x->numRows;
    tmp->p_ps_xCols = x->numCols;
    tmp->p_ps_yRows = y->numRows;
    tmp->p_ps_yCols = y->numCols;
    tmp->x = (double **) psAlloc(x->numRows * sizeof(double *));
    for (i=0;i<x->numRows;i++) {
        (tmp->x)[i] = (double *) psAlloc(x->numCols * sizeof(double));
    }
    for (i=0;i<x->numRows;i++) {
        for (j=0;j<x->numCols;j++) {
            (tmp->x)[i][j] = x->data.F64[i][j];
        }
    }

    tmp->y = (double **) psAlloc(y->numRows * sizeof(double *));
    for (i=0;i<y->numRows;i++) {
        (tmp->y)[i] = (double *) psAlloc(y->numCols * sizeof(double));
    }
    for (i=0;i<y->numRows;i++) {
        for (j=0;j<y->numCols;j++) {
            (tmp->y)[i][j] = y->data.F64[i][j];
        }
    }

    p_psMemSetDeallocator(tmp,(psFreeFcn)fixedPatternFree);

    return(tmp);
}


psExposure* psExposureAlloc(double ra,
                            double dec,
                            double hourAngle,
                            double zenithDistance,
                            double azimuth,
                            const psTime* time,
                            float rotAngle,
                            float temperature,
                            float pressure,
                            float humidity,
                            float exposureTime,
                            float wavelength,
                            const psObservatory* observatory)
{
    PS_PTR_CHECK_NULL(observatory, NULL);

    psExposure* exp = psAlloc(sizeof(psExposure));
    *(double *)&exp->ra = ra;
    *(double *)&exp->dec = dec;
    *(double *)&exp->hourAngle = hourAngle;
    *(double *)&exp->zenithDistance = zenithDistance;
    *(double *)&exp->azimuth = azimuth;
    *(float *)&exp->rotAngle = rotAngle;
    *(float *)&exp->temperature = temperature;
    *(float *)&exp->pressure = pressure;
    *(float *)&exp->humidity = humidity;
    *(float *)&exp->exposureTime = exposureTime;
    *(float *)&exp->wavelength = wavelength;

    exp->time = psMemIncrRefCounter((psPtr)time);
    exp->observatory = psMemIncrRefCounter((psPtr)observatory);

    // XXX: how is this value derived?
    *(double *)&exp->lst = psTimeToLST((psTime*)time,observatory->longitude);
    *(float *)&exp->positionAngle = 0.0f; // XXX: need input, see Bug #207
    *(float *)&exp->parallacticAngle = 0.0f; // XXX: need input, see Bug #207
    *(float *)&exp->airmass = slaAirmas(zenithDistance);
    *(float *)&exp->parallacticFactor = 0.0f;
    exp->cameraName = NULL;
    exp->telescopeName = NULL;

    p_psMemSetDeallocator(exp,(psFreeFcn)exposureFree);

    return exp;
}

psObservatory* psObservatoryAlloc(const char* name,
                                  double latitude,
                                  double longitude,
                                  double height,
                                  double tlr)
{
    psObservatory* obs = psAlloc(sizeof(psObservatory));

    if (name == NULL) {
        obs->name = NULL;
    } else {
        obs->name = psAlloc(strlen(name)+1);
        strcpy((char*)obs->name, name);
    }

    *(double *)&obs->latitude = latitude;
    *(double *)&obs->longitude = longitude;
    *(double *)&obs->height = height;
    *(double *)&obs->tlr = tlr;

    p_psMemSetDeallocator(obs,(psFreeFcn)observatoryFree);

    return obs;
}

psFPA* psFPAAlloc(psS32 nChips,
                  const psExposure* exp)
{
    PS_INT_CHECK_NON_NEGATIVE(nChips, NULL);

    psFPA* newFPA = psAlloc(sizeof(psFPA));

    // create array of NULL chips of the size nChips
    newFPA->chips = psArrayAlloc(nChips);
    psPtr* chips = newFPA->chips->data;
    for (psS32 i=0;i<nChips;i++) {
        chips[i] = NULL;
    }
    newFPA->chips->n = 0; // per requirement

    newFPA->metadata = NULL;
    newFPA->fromTangentPlane = NULL;
    newFPA->toTangentPlane = NULL;
    newFPA->pattern = NULL;

    if (exp != NULL) {
        newFPA->exposure = psMemIncrRefCounter((psExposure*)exp);
        newFPA->grommit = psGrommitAlloc(exp);
    } else {
        newFPA->exposure = NULL;
        newFPA->grommit = NULL;
    }

    newFPA->colorPlus = NULL;
    newFPA->colorMinus = NULL;
    newFPA->projection = NULL;

    newFPA->rmsX = 0.0f;
    newFPA->rmsY = 0.0f;
    newFPA->chi2 = 0.0f;

    p_psMemSetDeallocator(newFPA,(psFreeFcn)FPAFree);

    return newFPA;
}

/*
 * psChip constructor
 */
psChip* psChipAlloc(psS32 nCells,
                    psFPA *parentFPA)
{
    PS_INT_CHECK_NON_NEGATIVE(nCells, NULL);

    psChip* chip = psAlloc(sizeof(psChip));

    // create array of NULL psCells
    int n = (nCells > 0) ? nCells : 1;
    chip->cells = psArrayAlloc(n);
    psPtr* cells = chip->cells->data;
    for (psS32 i=0;i<n;i++) {
        cells[i] = NULL;
    }
    chip->cells->n = 0; // per requirement

    *(int*)&chip->row0 = 0;
    *(int*)&chip->col0 = 0;

    chip->metadata = NULL;

    chip->toFPA = NULL;
    chip->fromFPA = NULL;

    chip->parent = parentFPA;

    p_psMemSetDeallocator(chip,(psFreeFcn)chipFree);

    return chip;

}

/*
 * psCell constructor
 */
psCell* psCellAlloc(psS32 nReadouts,
                    psChip* parentChip)
{
    PS_INT_CHECK_NON_NEGATIVE(nReadouts, NULL);

    psCell* cell = psAlloc(sizeof(psCell));

    // create array of NULL psReadouts
    int n = (nReadouts > 0) ? nReadouts : 1;
    cell->readouts = psArrayAlloc(n);
    psPtr* readouts = cell->readouts->data;
    for (psS32 i=0;i<n;i++) {
        readouts[i] = NULL;
    }
    cell->readouts->n = 0; // per requirement

    *(int*)&cell->row0 = 0;
    *(int*)&cell->col0 = 0;

    cell->metadata = NULL;

    cell->toChip = NULL;
    cell->fromChip = NULL;
    cell->toFPA = NULL;
    cell->toTP = NULL;
    cell->toSky = NULL;

    cell->parent = parentChip;

    p_psMemSetDeallocator(cell,(psFreeFcn)cellFree);

    return cell;


}

psReadout* psReadoutAlloc()
{
    psReadout* readout = psAlloc(sizeof(psReadout));

    *(psU32*)&readout->colBins = 1;
    *(psU32*)&readout->rowBins = 1;
    *(psU32*)&readout->rowParity = 0;
    *(psU32*)&readout->colParity = 0;
    *(psS32*)&readout->col0 = 0;
    *(psS32*)&readout->row0 = 0;

    readout->image = NULL;
    readout->mask = NULL;
    readout->objects = NULL;
    readout->metadata = NULL;

    p_psMemSetDeallocator(readout,(psFreeFcn)readoutFree);

    return readout;
}

psGrommit* psGrommitAlloc(const psExposure* exp)
{
    PS_PTR_CHECK_NULL(exp, NULL);

    double date = psTimeToMJD(exp->time);
    double dut = psTimeGetUT1Delta(exp->time);
    double elongm = exp->observatory->longitude;
    double phim = exp->observatory->latitude;
    double hm = exp->observatory->height;

    psSphere* polarMotion = psTimeGetPoleCoords(exp->time);
    double xp = polarMotion->r;
    double yp = polarMotion->d;
    psFree(polarMotion);

    double tdk = exp->temperature;
    double pmb = exp->pressure;
    double rh = exp->humidity;
    double wl = exp->wavelength;
    double tlr = exp->observatory->tlr;

    psGrommit* grommit = (psGrommit* ) psAlloc(sizeof(psGrommit));

    slaAoppa(date, dut, elongm, phim, hm, xp, yp,
             tdk, pmb, rh, wl, tlr, (double*)grommit);

    return (grommit);
}

psCell* psCellInFPA(const psPlane* fpaCoord,
                    const psFPA* FPA)
{
    PS_PTR_CHECK_NULL(fpaCoord, NULL);
    PS_PTR_CHECK_NULL(FPA, NULL);

    psChip* tmpChip = NULL;
    psPlane chipCoord;
    psCell* outCell = NULL;

    // Determine which chip contains the fpaCoords.
    tmpChip = psChipInFPA(fpaCoord, FPA);
    if (tmpChip == NULL) {
        return(NULL);
    }

    // Convert to those chip coordinates.
    psCoordFPAToChip(&chipCoord, fpaCoord, tmpChip);

    // Determine which cell contains those chip coordinates.
    outCell = psCellInChip(&chipCoord, tmpChip);

    return (outCell);
}

psChip* psChipInFPA(const psPlane* fpaCoord,
                    const psFPA* FPA)
{
    PS_PTR_CHECK_NULL(fpaCoord, NULL);
    PS_PTR_CHECK_NULL(FPA, NULL);
    PS_PTR_CHECK_NULL(FPA->chips, NULL);

    psArray* chips = FPA->chips;
    psS32 nChips = chips->n;
    psPlane chipCoord;
    psCell *tmpCell = NULL;

    // Loop through every chip in this FPA.  Convert the original FPA
    // coordinates to chip coordinates for that chip.  Then, determine if any
    // cells in that chip contain those chip coordinates.

    for (psS32 i = 0; i < nChips; i++) {
        psChip* tmpChip = chips->data[i];
        PS_PTR_CHECK_NULL(tmpChip, NULL);
        PS_PTR_CHECK_NULL(tmpChip->fromFPA, NULL);

        psPlaneTransformApply(&chipCoord, tmpChip->fromFPA, fpaCoord);

        tmpCell = psCellInChip(&chipCoord, tmpChip);
        if (tmpCell != NULL) {
            return(tmpChip);
        }
    }

    // XXX: Print warning here?
    return (NULL);
}

psCell* psCellInChip(const psPlane* chipCoord,
                     const psChip* chip)
{
    PS_PTR_CHECK_NULL(chipCoord, NULL);
    PS_PTR_CHECK_NULL(chip, NULL);

    psPlane cellCoord;
    psArray* cells;

    cells = chip->cells;
    if (cells == NULL) {
        return NULL;
    }

    // We loop over each cell in the chip.  We transform the chipCoord into
    // a cellCoord for that cell and determine if that cellCoord is valid.
    // If so, then we return that cell.

    for (psS32 i = 0; i < cells->n; i++) {
        psCell* tmpCell = (psCell* ) cells->data[i];
        PS_PTR_CHECK_NULL(tmpCell, NULL);
        PS_PTR_CHECK_NULL(tmpCell->fromChip, NULL);
        psArray* readouts = tmpCell->readouts;

        if (readouts != NULL) {
            for (psS32 j = 0; j < readouts->n; j++) {
                psReadout* tmpReadout = readouts->data[j];
                PS_READOUT_CHECK_NULL(tmpReadout, NULL);

                psPlaneTransformApply(&cellCoord,
                                      tmpCell->fromChip,
                                      chipCoord);

                if (checkValidImageCoords(cellCoord.x,
                                          cellCoord.y,
                                          tmpReadout->image)) {
                    return (tmpCell);
                }
            }
        }
    }

    return (NULL);
}

psPlane* psCoordCellToChip(psPlane* outCoord,
                           const psPlane* inCoord,
                           const psCell* cell)
{
    PS_PTR_CHECK_NULL(inCoord, NULL);
    PS_PTR_CHECK_NULL(cell, NULL);

    return (psPlaneTransformApply(outCoord, cell->toChip, inCoord));
}

psPlane* psCoordChipToFPA(psPlane* outCoord,
                          const psPlane* inCoord,
                          const psChip* chip)
{
    PS_PTR_CHECK_NULL(inCoord, NULL);
    PS_PTR_CHECK_NULL(chip, NULL);

    return (psPlaneTransformApply(outCoord, chip->toFPA, inCoord));
}

psPlane* psCoordFPAToTP(psPlane* outCoord,
                        const psPlane* inCoord,
                        double color,
                        double magnitude,
                        const psFPA* fpa)
{
    PS_PTR_CHECK_NULL(inCoord, NULL);
    PS_PTR_CHECK_NULL(fpa, NULL);

    return(psPlaneDistortApply(outCoord, fpa->toTangentPlane, inCoord,
                               color, magnitude));
}

/*****************************************************************************
XXX: What about units for the (x,y) coords?
 *****************************************************************************/
psSphere* psCoordTPToSky(psSphere* outSphere,
                         const psPlane* tpCoord,
                         const psGrommit* grommit)
{
    PS_PTR_CHECK_NULL(tpCoord, NULL);
    PS_PTR_CHECK_NULL(grommit, NULL);

    double AOB = 0.0;
    double ZOB = 0.0;
    double HOB = 0.0;
    if (outSphere == NULL) {
        outSphere = (psSphere* ) psAlloc(sizeof(psSphere));
    }

    slaAopqk(tpCoord->x, tpCoord->y, (double*)grommit,
             &AOB, &ZOB, &HOB, &outSphere->r, &outSphere->d);

    return (outSphere);
}

psPlane* psCoordCellToFPA(psPlane* fpaCoord,
                          const psPlane* cellCoord,
                          const psCell* cell)
{
    PS_PTR_CHECK_NULL(cellCoord, NULL);
    PS_PTR_CHECK_NULL(cell, NULL);

    return (psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord));
}

psSphere* psCoordCellToSky(psSphere* skyCoord,
                           const psPlane* cellCoord,
                           double color,
                           double magnitude,
                           const psCell* cell)
{
    PS_PTR_CHECK_NULL(cellCoord, NULL);
    PS_PTR_CHECK_NULL(cell, NULL);
    PS_PTR_CHECK_NULL(cell->toFPA, NULL);
    PS_PTR_CHECK_NULL(cell->parent, NULL);
    PS_PTR_CHECK_NULL(cell->parent->parent, NULL);
    PS_PTR_CHECK_NULL(cell->parent->parent->toTangentPlane, NULL);
    PS_PTR_CHECK_NULL(cell->parent->parent->exposure, NULL);

    psPlane* fpaCoord = NULL;
    psPlane* tpCoord = NULL;
    psFPA* parFPA = (cell->parent)->parent;
    psGrommit* tmpGrommit = NULL;

    // Convert the input cell coordinates to FPA coordinates.
    fpaCoord = psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord);

    // Convert the FPA coordinates to tangent plane Coordinates.
    tpCoord = psPlaneDistortApply(tpCoord, parFPA->toTangentPlane,
                                  fpaCoord, color, magnitude);

    // Generate a grommit for this FPA.
    tmpGrommit = psGrommitAlloc(parFPA->exposure);

    // Convert the tangent plane Coordinates to sky coordinates.
    skyCoord = psCoordTPToSky(skyCoord, tpCoord, tmpGrommit);

    psFree(fpaCoord);
    psFree(tpCoord);
    psFree(tmpGrommit);

    return(skyCoord);
}

psSphere* psCoordCellToSkyQuick(psSphere* outSphere,
                                const psPlane* cellCoord,
                                const psCell* cell)
{
    PS_PTR_CHECK_NULL(cellCoord, NULL);
    PS_PTR_CHECK_NULL(cell, NULL);
    PS_PTR_CHECK_NULL(cell->toSky, NULL);
    PS_PTR_CHECK_NULL(cell->parent, NULL);
    PS_PTR_CHECK_NULL(cell->parent->parent, NULL);
    PS_PTR_CHECK_NULL(cell->parent->parent->projection, NULL);
    if (cell->toSky) {
        // XXX: Should we use toTP or toSky?
        psLogMsg(__func__, PS_LOG_WARN,
                 "WARNING: psCoordCellToSkyQuick(): The cell->toSky transform is ignored.  The cell->toTP transform is being used.");
    }

    psPlane *tpCoord = NULL;
    psChip *chip = cell->parent;
    psFPA *FPA = chip->parent;
    psProjectionType oldProjectionType;

    if (outSphere == NULL) {
        outSphere = (psSphere* ) psAlloc(sizeof(psSphere));
    }

    // Determine the tangent plane coordinates.
    tpCoord = psPlaneTransformApply(NULL, cell->toTP, cellCoord);

    // Save the old projection type and set the new projection type to TAN.
    oldProjectionType = FPA->projection->type;
    FPA->projection->type = PS_PROJ_TAN;

    // Deproject the tangent plane coordinates a sphere.
    outSphere = psDeproject(tpCoord, FPA->projection);

    // Restore old projection type.  Free memory.
    FPA->projection->type = oldProjectionType;
    psFree(tpCoord);

    return (outSphere);
}

/*****************************************************************************
XXX: What about units for the (x,y) coords?
 *****************************************************************************/
psPlane* psCoordSkyToTP(psPlane* tpCoord,
                        const psSphere* in,
                        const psGrommit* grommit)
{
    PS_PTR_CHECK_NULL(in, NULL);
    PS_PTR_CHECK_NULL(grommit, NULL);

    char* type = "RA";

    if (tpCoord == NULL) {
        tpCoord = (psPlane* ) psAlloc(sizeof(psPlane));
    }

    slaOapqk(type, in->r, in->d, (double*)grommit, &tpCoord->x, &tpCoord->y);

    return(tpCoord);
}


psPlane* psCoordTPToFPA(psPlane* fpaCoord,
                        const psPlane* tpCoord,
                        double color,
                        double magnitude,
                        const psFPA* fpa)
{
    PS_PTR_CHECK_NULL(tpCoord, NULL);
    PS_PTR_CHECK_NULL(fpa, NULL);
    PS_PTR_CHECK_NULL(fpa->fromTangentPlane, NULL);

    return (psPlaneDistortApply(fpaCoord, fpa->fromTangentPlane,
                                tpCoord, color, magnitude));
}

psPlane* psCoordFPAToChip(psPlane* chipCoord,
                          const psPlane* fpaCoord,
                          const psChip* chip)
{
    PS_PTR_CHECK_NULL(fpaCoord, NULL);
    PS_PTR_CHECK_NULL(chip, NULL);
    PS_PTR_CHECK_NULL(chip->fromFPA, NULL);

    chipCoord = psPlaneTransformApply(chipCoord, chip->fromFPA, fpaCoord);
    return(chipCoord);
}

psPlane* psCoordChipToCell(psPlane* cellCoord,
                           const psPlane* chipCoord,
                           const psCell* cell)
{
    PS_PTR_CHECK_NULL(chipCoord, NULL);
    PS_PTR_CHECK_NULL(cell, NULL);
    PS_PTR_CHECK_NULL(cell->fromChip, NULL);

    cellCoord = psPlaneTransformApply(cellCoord, cell->fromChip, chipCoord);
    return(cellCoord);
}

psPlane* psCoordSkyToCell(psPlane* cellCoord,
                          const psSphere* skyCoord,
                          double color,
                          double magnitude,
                          const psCell* cell)
{
    PS_PTR_CHECK_NULL(skyCoord, NULL);
    PS_PTR_CHECK_NULL(cell, NULL);
    PS_PTR_CHECK_NULL(cell->parent, NULL);
    PS_PTR_CHECK_NULL(cell->parent->parent, NULL);
    PS_PTR_CHECK_NULL(cell->parent->parent->grommit, NULL);

    psChip *parChip = cell->parent;
    psFPA *parFPA = parChip->parent;
    psGrommit* grommit = parFPA->grommit;

    // Convert the skyCoords to tangent plane coords.
    psPlane *tpCoord = psCoordSkyToTP(tpCoord, skyCoord, grommit);

    // Convert the tangent plane coords to FPA coords.
    psPlane *fpaCoord = psCoordTPToFPA(fpaCoord, tpCoord, color,
                                       magnitude, parFPA);

    // Convert the FPA coords to chip coords.
    psPlane *chipCoord = psCoordFPAToChip(chipCoord, fpaCoord, parChip);

    // Convert the chip coords to cell coords.
    cellCoord = psCoordChipToCell(cellCoord, chipCoord, cell);

    psFree(tpCoord);
    psFree(fpaCoord);
    psFree(chipCoord);

    return (cellCoord);
}

psPlane* psCoordSkyToCellQuick(psPlane* cellCoord,
                               const psSphere* skyCoord,
                               const psCell* cell)
{
    PS_PTR_CHECK_NULL(skyCoord, NULL);
    PS_PTR_CHECK_NULL(cell, NULL);
    PS_PTR_CHECK_NULL(cell->parent, NULL);
    PS_PTR_CHECK_NULL(cell->parent->parent, NULL);
    PS_PTR_CHECK_NULL(cell->parent->parent->projection, NULL);
    PS_PTR_CHECK_NULL(cell->toTP, NULL);
    if (cell->toSky) {
        // XXX: Should we use toTP or toSky?
        psLogMsg(__func__, PS_LOG_WARN,
                 "WARNING: psCoordSkyToCellQuick(): The cell->toSky transform is ignored.  The cell->toTP transform is being used.");
    }

    psPlane *tpCoord = NULL;
    psChip *whichChip = cell->parent;
    psFPA *whichFPA = whichChip->parent;
    psProjectionType oldProjectionType;
    psPlaneTransform *TPtoCell = NULL;

    // Save the old projection type and set the new projection type to TAN.
    oldProjectionType = whichFPA->projection->type;
    whichFPA->projection->type = PS_PROJ_TAN;

    if (cellCoord == NULL) {
        cellCoord = (psPlane* ) psAlloc(sizeof(psPlane));
    }

    tpCoord = psProject(skyCoord, whichFPA->projection);

    // generate an error if cell->toTP is not linear.
    if (0 == isProjectionLinear(cell->toTP)) {
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                PS_ERRORTEXT_psAstrometry_NONLINEAR_TRANSFORM,
                "cell to tangent plane");
    }

    TPtoCell = invertPlaneTransform(cell->toTP);
    cellCoord = psPlaneTransformApply(cellCoord, TPtoCell, tpCoord);

    // Restore old projection type.  Free memory.
    whichFPA->projection->type = oldProjectionType;
    psFree(tpCoord);
    return (cellCoord);
}

