IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1532


Ignore:
Timestamp:
Aug 13, 2004, 1:43:29 PM (22 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/astro/psCoord.c

    r1531 r1532  
    1010*  @author George Gusciora, MHPCC
    1111*
    12 *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2004-08-13 23:33:13 $
     12*  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2004-08-13 23:43:29 $
    1414*
    1515*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1616*/
    17 
     17/******************************************************************************/
     18/*  INCLUDE FILES                                                             */
     19/******************************************************************************/
    1820#include "psType.h"
    1921#include "psImage.h"
    … …  
    2628#include <math.h>
    2729#include <float.h>
    28 
    29 static float cot(float x);
    30 static float arg(float x, float y);
    31 
    32 psPlane* psPlaneTransformApply(psPlane* out,
    33                                const psPlaneTransform* transform,
    34                                const psPlane* coords)
    35 {
    36     if (out == NULL) {
    37         out = (psPlane* ) psAlloc(sizeof(psPlane));
    38     }
    39     out->x = transform->x->coeff[0][0] +
    40              (transform->x->coeff[1][0] * coords->x) +
    41              (transform->x->coeff[0][1] * coords->y);
    42 
    43     out->y = transform->y->coeff[0][0] +
    44              (transform->y->coeff[1][0] * coords->x) +
    45              (transform->y->coeff[0][1] * coords->y);
    46 
    47     return (out);
    48 }
    49 
    50 // This transformation takes into account parameters beyond an objects
    51 // spatial coordinates: term3 and term4 (magnitude and color).
    52 psPlane* psPlaneDistortApply(psPlane* out,
    53                              const psPlaneDistort* transform,
    54                              const psPlane* coords,
    55                              float term3,
    56                              float term4)
    57 {
    58     if (out == NULL) {
    59         out = (psPlane* ) psAlloc(sizeof(psPlane));
    60     }
    61 
    62     out->x = transform->x->coeff[0][0][0][0] +
    63              (transform->x->coeff[1][0][0][0] * coords->x) +
    64              (transform->x->coeff[0][1][0][0] * coords->y) +
    65              (transform->x->coeff[0][0][1][0] * term3) +
    66              (transform->x->coeff[0][0][0][1] * term4);
    67 
    68     out->y = transform->y->coeff[0][0][0][0] +
    69              (transform->y->coeff[1][0][0][0] * coords->x) +
    70              (transform->y->coeff[0][1][0][0] * coords->y) +
    71              (transform->y->coeff[0][0][1][0] * term3) +
    72              (transform->y->coeff[0][0][0][1] * term4);
    73 
    74     return (out);
    75 }
    76 
    77 // This function prototype has been modified since the SDRS.
    78 psSphereTransform* psSphereTransformAlloc(double NPlat,
    79         double Xo,
    80         double xo)
    81 {
    82     psSphereTransform* tmp = (psSphereTransform* ) psAlloc(sizeof(psSphereTransform));
    83 
    84     tmp->sinPhi = sin(NPlat);
    85     tmp->cosPhi = cos(NPlat);
    86     tmp->Xo = Xo;
    87     tmp->xo = xo;
    88 
    89     return (tmp);
    90 }
    91 
    92 /******************************************************************************
    93 This algorithm comes from an email from Gene.  I assume (x,y) corresponds to
    94 (r,d) in the sphere coordinates.
    95  
    96 XXX: In Gene's email, there are different variables with similar names (y
    97 and Y) and in the code, there's cos_y, cos_Y, and cos(y).  Verify that there
    98 are no typo's.
    99  *****************************************************************************/
    100 psSphere* psSphereTransformApply(psSphere* out,
    101                                  const psSphereTransform* transform,
    102                                  const psSphere* coord)
    103 {
    104     double sinY = 0.0;
    105     double cosY = 0.0;
    106     double sinX = 0.0;
    107     double cosX = 0.0;
    108     double x = 0.0;
    109     double y = 0.0;
    110     double dx = 0.0;
    111 
    112     if (out == NULL) {
    113         out = (psSphere* ) psAlloc(sizeof(psSphere));
    114     }
    115 
    116     x = coord->r;
    117     y = coord->d;
    118     dx = x - transform->xo;
    119     sinY = cos(y) * sin(dx) * transform->sinPhi + sin(y) * transform->cosPhi;
    120     cosY = sqrt(1.0 - sinY * sinY);
    121     sinX = (cos(y) * sin(dx) * transform->cosPhi - sin(y) * transform->sinPhi) / cos(y);
    122     cosX = cos(y) * cos(dx) / cos(y);
    123 
    124     out->r = atan2(sinX, cosX) + transform->Xo;
    125     out->d = atan2(sinY, cosY);
    126 
    127     return (out);
    128 }
    129 
    130 psSphereTransform* psSphereTransformICRStoEcliptic(psTime time)
    131 {
    132     struct tm *tmTime = psTimeToTM(time);
    133     double year = (double)(1900 + tmTime->tm_year);
    134     double T = year / 100.0;
    135     double phi = -23.452294 + 0.013013 * T + 0.000001639 * T * T - 0.000000503 * T * T * T;
    136     double Xo = 0.0;
    137     double xo = 0.0;
    138 
    139     return (psSphereTransformAlloc(phi, Xo, xo));
    140 }
    141 
    142 psSphereTransform* psSphereTransformEcliptictoICRS(psTime time)
    143 {
    144     struct tm *tmTime = psTimeToTM(time);
    145     double year = (double)(1900 + tmTime->tm_year);
    146     double T = year / 100.0;
    147     double phi = +23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T;
    148     double Xo = 0.0;
    149     double xo = 0.0;
    150 
    151     return (psSphereTransformAlloc(phi, Xo, xo));
    152 }
    153 
    154 psSphereTransform* psSphereTransformICRStoGalatic(void)
    155 {
    156     return (psSphereTransformAlloc(62.6, 282.25, 33.0));
    157 }
    158 
    159 psSphereTransform* psSphereTransformGalatictoICRS(void)
    160 {
    161     return (psSphereTransformAlloc(-62.6, 33.0, 282.25));
    162 }
    163 
    164 float cot(float x)
     30/******************************************************************************/
     31/*  DEFINE STATEMENTS                                                         */
     32/******************************************************************************/
     33
     34// None
     35
     36/******************************************************************************/
     37/*  TYPE DEFINITIONS                                                          */
     38/******************************************************************************/
     39
     40// None
     41
     42/*****************************************************************************/
     43/*  GLOBAL VARIABLES                                                         */
     44/*****************************************************************************/
     45
     46// None
     47
     48/*****************************************************************************/
     49/*  FILE STATIC VARIABLES                                                    */
     50/*****************************************************************************/
     51
     52// None
     53
     54/*****************************************************************************/
     55/*  FUNCTION IMPLEMENTATION - LOCAL                                          */
     56/*****************************************************************************/
     57
     58static float p_psCot(float x);
     59static float p_psArg(float x, float y);
     60
     61/******************************************************************************
     62XXX: Do this with a macro.
     63 *****************************************************************************/
     64float p_psCot(float x)
    16565{
    16666    return (1.0 / atan(x));
    … …  
    17070XXX: Verify this arc tan function.
    17171 *****************************************************************************/
    172 float arg(float x,
    173           float y)
     72float p_psArg(float x,
     73              float y)
    17474{
    17575    if (x > 0) {
    … …  
    18585    }
    18686
    187     psAbort(__func__, "Unacceptable range for (arg(%f, %f).\n", x, y);
     87    psAbort(__func__, "Unacceptable range for (p_psArg(%f, %f).\n", x, y);
    18888    return (0.0);
     89}
     90
     91/*****************************************************************************/
     92/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
     93/*****************************************************************************/
     94psPlane* psPlaneTransformApply(psPlane* out,
     95                               const psPlaneTransform* transform,
     96                               const psPlane* coords)
     97{
     98    if (out == NULL) {
     99        out = (psPlane* ) psAlloc(sizeof(psPlane));
     100    }
     101    out->x = transform->x->coeff[0][0] +
     102             (transform->x->coeff[1][0] * coords->x) +
     103             (transform->x->coeff[0][1] * coords->y);
     104
     105    out->y = transform->y->coeff[0][0] +
     106             (transform->y->coeff[1][0] * coords->x) +
     107             (transform->y->coeff[0][1] * coords->y);
     108
     109    return (out);
     110}
     111
     112/******************************************************************************
     113This transformation takes into account parameters beyond an objects spatial
     114coordinates: term3 and term4 (magnitude and color).
     115 *****************************************************************************/
     116psPlane* psPlaneDistortApply(psPlane* out,
     117                             const psPlaneDistort* transform,
     118                             const psPlane* coords,
     119                             float color,
     120                             float magnitude)
     121{
     122    if (out == NULL) {
     123        out = (psPlane* ) psAlloc(sizeof(psPlane));
     124    }
     125
     126    out->x = transform->x->coeff[0][0][0][0] +
     127             (transform->x->coeff[1][0][0][0] * coords->x) +
     128             (transform->x->coeff[0][1][0][0] * coords->y) +
     129             (transform->x->coeff[0][0][1][0] * color) +
     130             (transform->x->coeff[0][0][0][1] * magnitude);
     131
     132    out->y = transform->y->coeff[0][0][0][0] +
     133             (transform->y->coeff[1][0][0][0] * coords->x) +
     134             (transform->y->coeff[0][1][0][0] * coords->y) +
     135             (transform->y->coeff[0][0][1][0] * color) +
     136             (transform->y->coeff[0][0][0][1] * magnitude);
     137
     138    return (out);
     139}
     140
     141/******************************************************************************
     142This function prototype has been modified since the SDRS.
     143 *****************************************************************************/
     144psSphereTransform* psSphereTransformAlloc(double NPlat,
     145        double Xo,
     146        double xo)
     147{
     148    psSphereTransform* tmp = (psSphereTransform* ) psAlloc(sizeof(psSphereTransform));
     149
     150    tmp->sinPhi = sin(NPlat);
     151    tmp->cosPhi = cos(NPlat);
     152    tmp->Xo = Xo;
     153    tmp->xo = xo;
     154
     155    return (tmp);
     156}
     157
     158/******************************************************************************
     159This algorithm comes from an email from Gene.  I assume (x,y) corresponds to
     160(r,d) in the sphere coordinates.
     161 
     162XXX: In Gene's email, there are different variables with similar names (y
     163and Y) and in the code, there's cos_y, cos_Y, and cos(y).  Verify that there
     164are no typo's.
     165 *****************************************************************************/
     166psSphere* psSphereTransformApply(psSphere* out,
     167                                 const psSphereTransform* transform,
     168                                 const psSphere* coord)
     169{
     170    double sinY = 0.0;
     171    double cosY = 0.0;
     172    double sinX = 0.0;
     173    double cosX = 0.0;
     174    double x = 0.0;
     175    double y = 0.0;
     176    double dx = 0.0;
     177
     178    if (out == NULL) {
     179        out = (psSphere* ) psAlloc(sizeof(psSphere));
     180    }
     181
     182    x = coord->r;
     183    y = coord->d;
     184    dx = x - transform->xo;
     185    sinY = cos(y) * sin(dx) * transform->sinPhi + sin(y) * transform->cosPhi;
     186    cosY = sqrt(1.0 - sinY * sinY);
     187    sinX = (cos(y) * sin(dx) * transform->cosPhi - sin(y) * transform->sinPhi) / cos(y);
     188    cosX = cos(y) * cos(dx) / cos(y);
     189
     190    out->r = atan2(sinX, cosX) + transform->Xo;
     191    out->d = atan2(sinY, cosY);
     192
     193    return (out);
     194}
     195
     196psSphereTransform* psSphereTransformICRStoEcliptic(psTime time)
     197{
     198    struct tm *tmTime = psTimeToTM(time);
     199    double year = (double)(1900 + tmTime->tm_year);
     200    double T = year / 100.0;
     201    double phi = -23.452294 + 0.013013 * T + 0.000001639 * T * T - 0.000000503 * T * T * T;
     202    double Xo = 0.0;
     203    double xo = 0.0;
     204
     205    return (psSphereTransformAlloc(phi, Xo, xo));
     206}
     207
     208psSphereTransform* psSphereTransformEcliptictoICRS(psTime time)
     209{
     210    struct tm *tmTime = psTimeToTM(time);
     211    double year = (double)(1900 + tmTime->tm_year);
     212    double T = year / 100.0;
     213    double phi = +23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T;
     214    double Xo = 0.0;
     215    double xo = 0.0;
     216
     217    return (psSphereTransformAlloc(phi, Xo, xo));
     218}
     219
     220psSphereTransform* psSphereTransformICRStoGalatic(void)
     221{
     222    return (psSphereTransformAlloc(62.6, 282.25, 33.0));
     223}
     224
     225psSphereTransform* psSphereTransformGalatictoICRS(void)
     226{
     227    return (psSphereTransformAlloc(-62.6, 33.0, 282.25));
    189228}
    190229
    … …  
    201240
    202241    if (projection->type == PS_PROJ_TAN) {
    203         R = cot(coord->r) * (180.0 / M_PI);
     242        R = p_psCot(coord->r) * (180.0 / M_PI);
    204243        tmp->x = R * sin(coord->d);
    205244        tmp->y = R * cos(coord->d);
    … …  
    249288    if (projection->type == PS_PROJ_TAN) {
    250289        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
    251         tmp->d = arg(-coord->y, coord->x);
     290        tmp->d = p_psArg(-coord->y, coord->x);
    252291        tmp->r = atan(180.0 / (R * M_PI));
    253292
    254293    } else if (projection->type == PS_PROJ_SIN) {
    255294        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
    256         tmp->d = arg(-coord->y, coord->x);
     295        tmp->d = p_psArg(-coord->y, coord->x);
    257296        tmp->r = acos((R * M_PI) / 180.0);
    258297
    … …  
    271310        chu2 *= chu2;
    272311        chu = sqrt(1.0 - chu1 - chu2);
    273         tmp->d = 2.0 * arg((2.0 * chu * chu) - 1.0, (coord->x * chu * M_PI) / 360.0);
     312        tmp->d = 2.0 * p_psArg((2.0 * chu * chu) - 1.0, (coord->x * chu * M_PI) / 360.0);
    274313        tmp->r = asin((coord->y * chu * M_PI) / 180.0);
    275314
    … …  
    314353        lin = psProject(position2, &proj);
    315354        tmp = psDeproject(lin, &proj);
     355        psFree(lin);
    316356
    317357        // XXX: Do we need to convert units in tmp?
    … …  
    349389
    350390/******************************************************************************
    351 XXX: Do I need to check for unacceptable transformation parameters?  Maybe,
     391XXX: Do we need to check for unacceptable transformation parameters?  Maybe,
    352392     if the points are on the North/South Pole, etc?
    353393 
    354 XXX: Do I need to somehow scale this projection?
     394XXX: Do we need to somehow scale this projection?
    355395 
    356396XXX: I copied the algorithm from the ADD exactly.

  • trunk/psLib/src/astronomy/psAstrometry.c

    r1530 r1532  
    88*  @author George Gusciora, MHPCC
    99*
    10 *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
    11 *  @date $Date: 2004-08-13 22:41:24 $
     10*  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2004-08-13 23:43:29 $
    1212*
    1313*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1414*/
    15 
     15/******************************************************************************/
     16/*  INCLUDE FILES                                                             */
     17/******************************************************************************/
    1618#include<math.h>
    17 
    1819#include "psFunctions.h"
    1920#include "psAstrometry.h"
    … …  
    2122#include "psAbort.h"
    2223#include "slalib.h"
    23 
    24 psExposure* psExposureAlloc(double ra, double dec, double hourAngle,
    25                             double zenith, double azimuth, double localTime, float date,
    26                             float rotAngle, float temperature, float pressure, float humidity,
    27                             float exposureTime)
    28 {
    29     psExposure* exp = psAlloc(sizeof(psExposure));
    30 
    31     *(double *)&exp->ra = ra;
    32     *(double *)&exp->dec = dec;
    33     *(double *)&exp->hourAngle = hourAngle;
    34     *(double *)&exp->zenith = zenith;
    35     *(double *)&exp->azimuth = azimuth;
    36     *(double *)&exp->localTime = localTime;
    37     *(float *)&exp->date = date;
    38     *(float *)&exp->rotAngle = rotAngle;
    39     *(float *)&exp->temperature = temperature;
    40     *(float *)&exp->pressure = pressure;
    41     *(float *)&exp->humidity = humidity;
    42     *(float *)&exp->exposureTime = exposureTime;
    43 
    44     return exp;
    45 }
    46 #define TBD 0.0
    47 /*****************************************************************************
    48 XXX: Several of the input params to sla_aoppa() are currently undefined.
    49 We are awaiting futher direction from IfA on this.
    50  *****************************************************************************/
    51 psGrommit* psGrommitAlloc(const psExposure* exp)
    52 {
    53     if (exp == NULL) {
    54         psAbort(__func__, "the 'exp' parameter is NULL\n");
    55     }
    56 
    57     double date = TBD;  // "mjd" in psExposure will become a psTime
    58     // from which it will be possible to get UTC.
    59     double dut = 0.0;
    60     double elongm = TBD;
    61     double phim = TBD;
    62     double hm = TBD;
    63     double xp = 0.0;
    64     double yp = 0.0;
    65     double tdk = exp->temperature;
    66     double pmb = exp->pressure;
    67     double rh = exp->humidity;
    68     double wl = exp->wavelength;
    69     double tlr = TBD;
    70     double *AOPRMS = NULL;
    71 
    72     slaAoppa(date, dut, elongm, phim, hm, xp, yp,
    73              tdk, pmb, rh, wl, tlr, AOPRMS);
    74 
    75     psGrommit* grommit = (psGrommit* ) psAlloc(sizeof(psGrommit));
    76     *(double *)&grommit->latitude = AOPRMS[0];
    77     *(double *)&grommit->sinLat = AOPRMS[1];
    78     *(double *)&grommit->cosLat = AOPRMS[2];
    79     *(double *)&grommit->height = AOPRMS[3];
    80     *(double *)&grommit->abberationMag = AOPRMS[4];
    81     *(double *)&grommit->temperature = AOPRMS[5];
    82     *(double *)&grommit->pressure = AOPRMS[6];
    83     *(double *)&grommit->humidity = AOPRMS[7];
    84     *(double *)&grommit->wavelength = AOPRMS[8];
    85     *(double *)&grommit->lapseRate = AOPRMS[9];
    86     *(double *)&grommit->refractA = AOPRMS[10];
    87     *(double *)&grommit->refractB = AOPRMS[11];
    88     *(double *)&grommit->longitudeOffset = AOPRMS[12];
    89     *(double *)&grommit->siderealTime = AOPRMS[13];
    90 
    91     return (grommit);
    92 }
     24/******************************************************************************/
     25/*  DEFINE STATEMENTS                                                         */
     26/******************************************************************************/
     27
     28// None
     29
     30/******************************************************************************/
     31/*  TYPE DEFINITIONS                                                          */
     32/******************************************************************************/
     33
     34// None
     35
     36/*****************************************************************************/
     37/*  GLOBAL VARIABLES                                                         */
     38/*****************************************************************************/
     39
     40// None
     41
     42/*****************************************************************************/
     43/*  FILE STATIC VARIABLES                                                    */
     44/*****************************************************************************/
     45
     46// None
     47
     48/*****************************************************************************/
     49/*  FUNCTION IMPLEMENTATION - LOCAL                                          */
     50/*****************************************************************************/
    9351
    9452void p_psGrommitFree(psGrommit* grommit)
    … …  
    9755}
    9856
    99 psCell* psCellInFPA(const psPlane* fpaCoord,
    100                     const psFPA* FPA)
    101 {
    102     psChip* tmpChip = NULL;
    103     psPlane* chipCoord = NULL;
    104     psCell* outCell = NULL;
    105 
    106     if (fpaCoord == NULL) {
    107         psAbort(__func__, "input parameter fpaCoord is NULL.");
    108     }
    109     if (FPA == NULL) {
    110         psAbort(__func__, "input parameter FPA is NULL.");
    111     }
    112 
    113     // Determine which chip contains the fpaCoords.
    114     tmpChip = psChipInFPA(fpaCoord, FPA);
    115 
    116     // Convert to those chip coordinates.
    117     chipCoord = psCoordFPAToChip(chipCoord, fpaCoord, tmpChip);
    118 
    119     // Determine which cell contains those chip coordinates.
    120     outCell = psCellInChip(chipCoord, tmpChip);
    121 
    122     psFree(tmpChip);
    123     psFree(chipCoord);
    124 
    125     return (outCell);
    126 }
    127 
     57/*****************************************************************************
     58p_psCheckValidImageCoords(): this is a private function which simply
     59determines if the supplied x,y coordinates are in the range for the supplied
     60psImage.
     61 *****************************************************************************/
    12862int p_psCheckValidImageCoords(double x,
    12963                              double y,
    … …  
    14276
    14377    return (1);
    144 }
    145 
    146 psChip* psChipInFPA(const psPlane* fpaCoord,
    147                     const psFPA* FPA)
    148 {
    149     psArray* chips = FPA->chips;
    150     int nChips = chips->n;
    151     psPlane* chipCoord = NULL;
    152     psCell *tmpCell = NULL;
    153 
    154     if (fpaCoord == NULL) {
    155         psAbort(__func__, "input parameter fpaCoord is NULL.");
    156     }
    157     if (FPA == NULL) {
    158         psAbort(__func__, "input parameter FPA is NULL.");
    159     }
    160 
    161     // Loop through every chip in this FPA.  Convert the original
    162     // FPA coordinates to chip coordinates for that chip.  Then,
    163     // determine if any cells in that chip contain those chip
    164     // coordinates.
    165     for (int i = 0; i < nChips; i++) {
    166         psChip* tmpChip = chips->data[i];
    167         chipCoord = psPlaneTransformApply(chipCoord, tmpChip->fromFPA,
    168                                           fpaCoord);
    169 
    170         tmpCell = psCellInChip(chipCoord, tmpChip);
    171         if (tmpCell != NULL) {
    172             psFree(chipCoord);
    173             return(tmpChip);
    174         }
    175         psFree(chipCoord);
    176     }
    177 
    178     return (NULL);
    179 }
    180 
    181 psCell* psCellInChip(const psPlane* chipCoord,
    182                      const psChip* chip)
    183 {
    184     psPlane* cellCoord = NULL;
    185     psArray* cells;
    186 
    187     // We return NULL if either of the input parameters is NULL.
    188     if (chipCoord == NULL) {
    189         psAbort(__func__, "the 'chipCoord' parameter is NULL\n");
    190     }
    191     if (chip == NULL) {
    192         psAbort(__func__, "the 'chip' parameter is NULL\n");
    193     }
    194 
    195     cells = chip->cells;
    196     if (cells == NULL) {
    197         return NULL;
    198     }
    199 
    200     // We loop over each cell in the chip.  We transform the chipCoord into
    201     // a cellCoord for that cell and determine if that cellCoord is valid.
    202     // If so, then we return that cell.
    203     for (int i = 0; i < cells->n; i++) {
    204         psCell* tmpCell = (psCell* ) cells->data[i];
    205         psArray* readouts = tmpCell->readouts;
    206 
    207         if (readouts != NULL) {
    208             for (int j = 0; j < readouts->n; j++) {
    209                 psReadout* tmpReadout = readouts->data[j];
    210 
    211                 cellCoord = psPlaneTransformApply(cellCoord,
    212                                                   tmpCell->fromChip,
    213                                                   chipCoord);
    214 
    215                 if (p_psCheckValidImageCoords(cellCoord->x,
    216                                               cellCoord->y,
    217                                               tmpReadout->image)) {
    218                     psFree(cellCoord);
    219                     return (tmpCell);
    220                 }
    221             }
    222         }
    223     }
    224 
    225     psFree(cellCoord);
    226     return (NULL);
    227 }
    228 
    229 psPlane* psCoordCellToChip(psPlane* outCoord,
    230                            const psPlane* inCoord,
    231                            const psCell* cell)
    232 {
    233     if (inCoord == NULL) {
    234         psAbort(__func__, "input parameter inCoord is NULL.");
    235     }
    236     if (cell == NULL) {
    237         psAbort(__func__, "input parameter cell is NULL.");
    238     }
    239 
    240     return (psPlaneTransformApply(outCoord, cell->toChip, inCoord));
    241 }
    242 
    243 psPlane* psCoordChipToFPA(psPlane* outCoord,
    244                           const psPlane* inCoord,
    245                           const psChip* chip)
    246 {
    247     if (inCoord == NULL) {
    248         psAbort(__func__, "input parameter inCoord is NULL.");
    249     }
    250     if (chip == NULL) {
    251         psAbort(__func__, "input parameter chip is NULL.");
    252     }
    253 
    254     return (psPlaneTransformApply(outCoord, chip->toFPA, inCoord));
    255 }
    256 
    257 psPlane* psCoordFPAToTP(psPlane* outCoord,
    258                         const psPlane* inCoord,
    259                         double color,
    260                         double magnitude,
    261                         const psFPA* fpa)
    262 {
    263     if (inCoord == NULL) {
    264         psAbort(__func__, "input parameter inCoord is NULL.");
    265     }
    266     if (fpa == NULL) {
    267         psAbort(__func__, "input parameter fpa is NULL.");
    268     }
    269 
    270     return(psPlaneDistortApply(outCoord, fpa->toTangentPlane, inCoord,
    271                                color, magnitude));
    272 }
    273 
    274 /*****************************************************************************
    275 XXX: What about units for the (x,y) coords?
    276  *****************************************************************************/
    277 psSphere* psCoordTPToSky(psSphere* outSphere,
    278                          const psPlane* tpCoord,
    279                          const psGrommit* grommit)
    280 {
    281     double AOB = 0.0;
    282     double ZOB = 0.0;
    283     double HOB = 0.0;
    284 
    285     if (tpCoord == NULL) {
    286         psAbort(__func__, "input parameter tpCoord is NULL.");
    287     }
    288     if (grommit == NULL) {
    289         psAbort(__func__, "input parameter grommit is NULL.");
    290     }
    291     if (outSphere == NULL) {
    292         outSphere = (psSphere* ) psAlloc(sizeof(psSphere));
    293     }
    294 
    295     slaAopqk(tpCoord->x, tpCoord->y, (double *) grommit,
    296              &AOB, &ZOB, &HOB, &outSphere->r, &outSphere->d);
    297 
    298     return (outSphere);
    299 }
    300 
    301 psPlane* psCoordCellToFPA(psPlane* fpaCoord,
    302                           const psPlane* cellCoord,
    303                           const psCell* cell)
    304 {
    305     if (cellCoord == NULL) {
    306         psAbort(__func__, "input parameter cellCoord is NULL.");
    307     }
    308     if (cell == NULL) {
    309         psAbort(__func__, "input parameter cell is NULL.");
    310     }
    311 
    312     return (psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord));
    313 }
    314 
    315 psSphere* psCoordCellToSky(psSphere* skyCoord,
    316                            const psPlane* cellCoord,
    317                            double color,
    318                            double magnitude,
    319                            const psCell* cell)
    320 {
    321     if (cellCoord == NULL) {
    322         psAbort(__func__, "input parameter cellCoord is NULL.");
    323     }
    324     if (cell == NULL) {
    325         psAbort(__func__, "input parameter cell is NULL.");
    326     }
    327 
    328     psPlane* fpaCoord = NULL;
    329     psPlane* tpCoord = NULL;
    330     psFPA* parFPA = (cell->parent)->parent;
    331     psGrommit* tmpGrommit = NULL;
    332 
    333     // Convert the input cell coordinates to FPA coordinates.
    334     fpaCoord = psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord);
    335 
    336     // Convert the FPA coordinates to tangent plane Coordinates.
    337     tpCoord = psPlaneDistortApply(tpCoord, parFPA->toTangentPlane,
    338                                   fpaCoord, color, magnitude);
    339 
    340     // Generate a grommit for this FPA.
    341     tmpGrommit = psGrommitAlloc(parFPA->exposure);
    342 
    343     // Convert the tangent plane Coordinates to sky coordinates.
    344     skyCoord = psCoordTPToSky(skyCoord, tpCoord, tmpGrommit);
    345 
    346     psFree(fpaCoord);
    347     psFree(tpCoord);
    348     psFree(tmpGrommit);
    349 
    350     return(skyCoord);
    351 }
    352 
    353 psSphere* psCoordCellToSkyQuick(psSphere* outSphere,
    354                                 const psPlane* cellCoord,
    355                                 const psCell* cell)
    356 {
    357     if (cellCoord == NULL) {
    358         psAbort(__func__, "input parameter cellCoord is NULL.");
    359     }
    360     if (cell == NULL) {
    361         psAbort(__func__, "input parameter cell is NULL.");
    362     }
    363 
    364     psPlane *tpCoord = NULL;
    365     psChip *chip = cell->parent;
    366     psFPA *FPA = chip->parent;
    367     psProjectionType oldProjectionType;
    368 
    369     if (outSphere == NULL) {
    370         outSphere = (psSphere* ) psAlloc(sizeof(psSphere));
    371     }
    372 
    373     // Determine the tangent plane coordinates.
    374     tpCoord = psPlaneTransformApply(tpCoord, cell->toTP, cellCoord);
    375 
    376     // Save the old projection type and set the new projection type to TAN.
    377     oldProjectionType = FPA->projection->type;
    378     FPA->projection->type = PS_PROJ_TAN;
    379 
    380     // Deproject the tangent plane coordinates a sphere.
    381     outSphere = psDeproject(tpCoord, FPA->projection);
    382 
    383     // Restore old projection type.  Free memory.
    384     FPA->projection->type = oldProjectionType;
    385     psFree(tpCoord);
    386 
    387     return (outSphere);
    388 }
    389 
    390 /*****************************************************************************
    391 XXX: What about units for the (x,y) coords?
    392  *****************************************************************************/
    393 psPlane* psCoordSkyToTP(psPlane* tpCoord,
    394                         const psSphere* in,
    395                         const psGrommit* grommit)
    396 {
    397     if (in == NULL) {
    398         psAbort(__func__, "input parameter in is NULL.");
    399     }
    400     if (grommit == NULL) {
    401         psAbort(__func__, "input parameter grommit is NULL.");
    402     }
    403     if (tpCoord == NULL) {
    404         tpCoord = (psPlane* ) psAlloc(sizeof(psPlane));
    405     }
    406 
    407     slaOapqk("RA", in->r, in->d, (double *) grommit, &tpCoord->x, &tpCoord->y);
    408 
    409     return(tpCoord);
    410 }
    411 
    412 
    413 psPlane* psCoordTPToFPA(psPlane* fpaCoord,
    414                         const psPlane* tpCoord,
    415                         double color,
    416                         double magnitude,
    417                         const psFPA* fpa)
    418 {
    419     if (tpCoord == NULL) {
    420         psAbort(__func__, "input parameter tpCoord is NULL.");
    421     }
    422     if (fpa == NULL) {
    423         psAbort(__func__, "input parameter fpa is NULL.");
    424     }
    425 
    426     return (psPlaneDistortApply(fpaCoord, fpa->fromTangentPlane,
    427                                 tpCoord, color, magnitude));
    428 }
    429 
    430 psPlane* psCoordFPAToChip(psPlane* chipCoord,
    431                           const psPlane* fpaCoord,
    432                           const psChip* chip)
    433 {
    434     if (fpaCoord == NULL) {
    435         psAbort(__func__, "input parameter fpaCoord is NULL.");
    436     }
    437     if (chip == NULL) {
    438         psAbort(__func__, "input parameter chip is NULL.");
    439     }
    440 
    441     psFPA *FPA = chip->parent;
    442     if (FPA == NULL) {
    443         psAbort(__func__, "chip->parent is NULL");
    444     }
    445 
    446     // Determine which chip contains these FPA coordinates.
    447     psChip *newChip = psChipInFPA(fpaCoord, FPA);
    448     if (newChip == NULL) {
    449         return(NULL);
    450     }
    451 
    452     chipCoord = psPlaneTransformApply(chipCoord, newChip->fromFPA, fpaCoord);
    453     psFree(newChip);
    454     return(chipCoord);
    455 }
    456 
    457 psPlane* psCoordChipToCell(psPlane* cellCoord,
    458                            const psPlane* chipCoord,
    459                            const psCell* cell)
    460 {
    461     if (chipCoord == NULL) {
    462         psAbort(__func__, "input parameter chipCoord is NULL.");
    463     }
    464     if (cell == NULL) {
    465         psAbort(__func__, "input parameter cell is NULL.");
    466     }
    467 
    468     psChip *chip = cell->parent;
    469     if (chip == NULL) {
    470         psAbort(__func__, "cell->parent is NULL");
    471     }
    472 
    473     // Determine which cell contains these FPA coordinates.
    474     psCell *newCell = psCellInChip(chipCoord, chip);
    475     if (newCell == NULL) {
    476         return(NULL);
    477     }
    478 
    479     cellCoord = psPlaneTransformApply(cellCoord, newCell->fromChip, chipCoord);
    480     psFree(newCell);
    481     return(cellCoord);
    482 }
    483 
    484 /*****************************************************************************
    485 XXX: Should we do anything to determine which cell contains the coordinates?
    486  *****************************************************************************/
    487 psPlane* psCoordSkyToCell(psPlane* cellCoord,
    488                           const psSphere* skyCoord,
    489                           double color,
    490                           double magnitude,
    491                           const psCell* cell)
    492 {
    493     if (skyCoord == NULL) {
    494         psAbort(__func__, "input parameter skyCoord is NULL.");
    495     }
    496     if (cell == NULL) {
    497         psAbort(__func__, "input parameter cell is NULL.");
    498     }
    499 
    500     psChip *parChip = cell->parent;
    501     psFPA *parFPA = parChip->parent;
    502     psGrommit* grommit = parFPA->grommit;
    503 
    504     // Convert the skyCoords to tangent plane coords.
    505     psPlane *tpCoord = psCoordSkyToTP(tpCoord, skyCoord, grommit);
    506 
    507     // Convert the tangent plane coords to FPA coords.
    508     psPlane *fpaCoord = psCoordTPToFPA(fpaCoord, tpCoord, color,
    509                                        magnitude, parFPA);
    510 
    511     // Convert the FPA coords to chip coords.
    512     psPlane *chipCoord = psCoordFPAToChip(chipCoord, fpaCoord, parChip);
    513 
    514     // Convert the chip coords to cell coords.
    515     cellCoord = psCoordChipToCell(cellCoord, chipCoord, cell);
    516 
    517     psFree(tpCoord);
    518     psFree(fpaCoord);
    519     psFree(chipCoord);
    520 
    521     return (cellCoord);
    52278}
    52379
    … …  
    602158}
    603159
     160/*****************************************************************************/
     161/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
     162/*****************************************************************************/
     163psExposure* psExposureAlloc(double ra, double dec, double hourAngle,
     164                            double zenith, double azimuth, double localTime, float date,
     165                            float rotAngle, float temperature, float pressure, float humidity,
     166                            float exposureTime)
     167{
     168    psExposure* exp = psAlloc(sizeof(psExposure));
     169
     170    *(double *)&exp->ra = ra;
     171    *(double *)&exp->dec = dec;
     172    *(double *)&exp->hourAngle = hourAngle;
     173    *(double *)&exp->zenith = zenith;
     174    *(double *)&exp->azimuth = azimuth;
     175    *(double *)&exp->localTime = localTime;
     176    *(float *)&exp->date = date;
     177    *(float *)&exp->rotAngle = rotAngle;
     178    *(float *)&exp->temperature = temperature;
     179    *(float *)&exp->pressure = pressure;
     180    *(float *)&exp->humidity = humidity;
     181    *(float *)&exp->exposureTime = exposureTime;
     182
     183    return exp;
     184}
     185#define TBD 0.0
     186/*****************************************************************************
     187XXX: Several of the input params to sla_aoppa() are currently undefined.
     188We are awaiting futher direction from IfA on this.
     189 *****************************************************************************/
     190psGrommit* psGrommitAlloc(const psExposure* exp)
     191{
     192    if (exp == NULL) {
     193        psAbort(__func__, "the 'exp' parameter is NULL\n");
     194    }
     195
     196    double date = TBD;  // "mjd" in psExposure will become a psTime
     197    // from which it will be possible to get UTC.
     198    double dut = 0.0;
     199    double elongm = TBD;
     200    double phim = TBD;
     201    double hm = TBD;
     202    double xp = 0.0;
     203    double yp = 0.0;
     204    double tdk = exp->temperature;
     205    double pmb = exp->pressure;
     206    double rh = exp->humidity;
     207    double wl = exp->wavelength;
     208    double tlr = TBD;
     209    double *AOPRMS = NULL;
     210
     211    slaAoppa(date, dut, elongm, phim, hm, xp, yp,
     212             tdk, pmb, rh, wl, tlr, AOPRMS);
     213
     214    psGrommit* grommit = (psGrommit* ) psAlloc(sizeof(psGrommit));
     215    *(double *)&grommit->latitude = AOPRMS[0];
     216    *(double *)&grommit->sinLat = AOPRMS[1];
     217    *(double *)&grommit->cosLat = AOPRMS[2];
     218    *(double *)&grommit->height = AOPRMS[3];
     219    *(double *)&grommit->abberationMag = AOPRMS[4];
     220    *(double *)&grommit->temperature = AOPRMS[5];
     221    *(double *)&grommit->pressure = AOPRMS[6];
     222    *(double *)&grommit->humidity = AOPRMS[7];
     223    *(double *)&grommit->wavelength = AOPRMS[8];
     224    *(double *)&grommit->lapseRate = AOPRMS[9];
     225    *(double *)&grommit->refractA = AOPRMS[10];
     226    *(double *)&grommit->refractB = AOPRMS[11];
     227    *(double *)&grommit->longitudeOffset = AOPRMS[12];
     228    *(double *)&grommit->siderealTime = AOPRMS[13];
     229
     230    return (grommit);
     231}
     232
     233psCell* psCellInFPA(const psPlane* fpaCoord,
     234                    const psFPA* FPA)
     235{
     236    psChip* tmpChip = NULL;
     237    psPlane* chipCoord = NULL;
     238    psCell* outCell = NULL;
     239
     240    if (fpaCoord == NULL) {
     241        psAbort(__func__, "input parameter fpaCoord is NULL.");
     242    }
     243    if (FPA == NULL) {
     244        psAbort(__func__, "input parameter FPA is NULL.");
     245    }
     246
     247    // Determine which chip contains the fpaCoords.
     248    tmpChip = psChipInFPA(fpaCoord, FPA);
     249
     250    // Convert to those chip coordinates.
     251    chipCoord = psCoordFPAToChip(chipCoord, fpaCoord, tmpChip);
     252
     253    // Determine which cell contains those chip coordinates.
     254    outCell = psCellInChip(chipCoord, tmpChip);
     255
     256    psFree(tmpChip);
     257    psFree(chipCoord);
     258
     259    return (outCell);
     260}
     261
     262psChip* psChipInFPA(const psPlane* fpaCoord,
     263                    const psFPA* FPA)
     264{
     265    psArray* chips = FPA->chips;
     266    int nChips = chips->n;
     267    psPlane* chipCoord = NULL;
     268    psCell *tmpCell = NULL;
     269
     270    if (fpaCoord == NULL) {
     271        psAbort(__func__, "input parameter fpaCoord is NULL.");
     272    }
     273    if (FPA == NULL) {
     274        psAbort(__func__, "input parameter FPA is NULL.");
     275    }
     276
     277    // Loop through every chip in this FPA.  Convert the original
     278    // FPA coordinates to chip coordinates for that chip.  Then,
     279    // determine if any cells in that chip contain those chip
     280    // coordinates.
     281    for (int i = 0; i < nChips; i++) {
     282        psChip* tmpChip = chips->data[i];
     283        chipCoord = psPlaneTransformApply(chipCoord, tmpChip->fromFPA,
     284                                          fpaCoord);
     285
     286        tmpCell = psCellInChip(chipCoord, tmpChip);
     287        if (tmpCell != NULL) {
     288            psFree(chipCoord);
     289            return(tmpChip);
     290        }
     291        psFree(chipCoord);
     292    }
     293
     294    return (NULL);
     295}
     296
     297psCell* psCellInChip(const psPlane* chipCoord,
     298                     const psChip* chip)
     299{
     300    psPlane* cellCoord = NULL;
     301    psArray* cells;
     302
     303    // We return NULL if either of the input parameters is NULL.
     304    if (chipCoord == NULL) {
     305        psAbort(__func__, "the 'chipCoord' parameter is NULL\n");
     306    }
     307    if (chip == NULL) {
     308        psAbort(__func__, "the 'chip' parameter is NULL\n");
     309    }
     310
     311    cells = chip->cells;
     312    if (cells == NULL) {
     313        return NULL;
     314    }
     315
     316    // We loop over each cell in the chip.  We transform the chipCoord into
     317    // a cellCoord for that cell and determine if that cellCoord is valid.
     318    // If so, then we return that cell.
     319    for (int i = 0; i < cells->n; i++) {
     320        psCell* tmpCell = (psCell* ) cells->data[i];
     321        psArray* readouts = tmpCell->readouts;
     322
     323        if (readouts != NULL) {
     324            for (int j = 0; j < readouts->n; j++) {
     325                psReadout* tmpReadout = readouts->data[j];
     326
     327                cellCoord = psPlaneTransformApply(cellCoord,
     328                                                  tmpCell->fromChip,
     329                                                  chipCoord);
     330
     331                if (p_psCheckValidImageCoords(cellCoord->x,
     332                                              cellCoord->y,
     333                                              tmpReadout->image)) {
     334                    psFree(cellCoord);
     335                    return (tmpCell);
     336                }
     337            }
     338        }
     339    }
     340
     341    psFree(cellCoord);
     342    return (NULL);
     343}
     344
     345psPlane* psCoordCellToChip(psPlane* outCoord,
     346                           const psPlane* inCoord,
     347                           const psCell* cell)
     348{
     349    if (inCoord == NULL) {
     350        psAbort(__func__, "input parameter inCoord is NULL.");
     351    }
     352    if (cell == NULL) {
     353        psAbort(__func__, "input parameter cell is NULL.");
     354    }
     355
     356    return (psPlaneTransformApply(outCoord, cell->toChip, inCoord));
     357}
     358
     359psPlane* psCoordChipToFPA(psPlane* outCoord,
     360                          const psPlane* inCoord,
     361                          const psChip* chip)
     362{
     363    if (inCoord == NULL) {
     364        psAbort(__func__, "input parameter inCoord is NULL.");
     365    }
     366    if (chip == NULL) {
     367        psAbort(__func__, "input parameter chip is NULL.");
     368    }
     369
     370    return (psPlaneTransformApply(outCoord, chip->toFPA, inCoord));
     371}
     372
     373psPlane* psCoordFPAToTP(psPlane* outCoord,
     374                        const psPlane* inCoord,
     375                        double color,
     376                        double magnitude,
     377                        const psFPA* fpa)
     378{
     379    if (inCoord == NULL) {
     380        psAbort(__func__, "input parameter inCoord is NULL.");
     381    }
     382    if (fpa == NULL) {
     383        psAbort(__func__, "input parameter fpa is NULL.");
     384    }
     385
     386    return(psPlaneDistortApply(outCoord, fpa->toTangentPlane, inCoord,
     387                               color, magnitude));
     388}
     389
     390/*****************************************************************************
     391XXX: What about units for the (x,y) coords?
     392 *****************************************************************************/
     393psSphere* psCoordTPToSky(psSphere* outSphere,
     394                         const psPlane* tpCoord,
     395                         const psGrommit* grommit)
     396{
     397    double AOB = 0.0;
     398    double ZOB = 0.0;
     399    double HOB = 0.0;
     400
     401    if (tpCoord == NULL) {
     402        psAbort(__func__, "input parameter tpCoord is NULL.");
     403    }
     404    if (grommit == NULL) {
     405        psAbort(__func__, "input parameter grommit is NULL.");
     406    }
     407    if (outSphere == NULL) {
     408        outSphere = (psSphere* ) psAlloc(sizeof(psSphere));
     409    }
     410
     411    slaAopqk(tpCoord->x, tpCoord->y, (double *) grommit,
     412             &AOB, &ZOB, &HOB, &outSphere->r, &outSphere->d);
     413
     414    return (outSphere);
     415}
     416
     417psPlane* psCoordCellToFPA(psPlane* fpaCoord,
     418                          const psPlane* cellCoord,
     419                          const psCell* cell)
     420{
     421    if (cellCoord == NULL) {
     422        psAbort(__func__, "input parameter cellCoord is NULL.");
     423    }
     424    if (cell == NULL) {
     425        psAbort(__func__, "input parameter cell is NULL.");
     426    }
     427
     428    return (psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord));
     429}
     430
     431psSphere* psCoordCellToSky(psSphere* skyCoord,
     432                           const psPlane* cellCoord,
     433                           double color,
     434                           double magnitude,
     435                           const psCell* cell)
     436{
     437    if (cellCoord == NULL) {
     438        psAbort(__func__, "input parameter cellCoord is NULL.");
     439    }
     440    if (cell == NULL) {
     441        psAbort(__func__, "input parameter cell is NULL.");
     442    }
     443
     444    psPlane* fpaCoord = NULL;
     445    psPlane* tpCoord = NULL;
     446    psFPA* parFPA = (cell->parent)->parent;
     447    psGrommit* tmpGrommit = NULL;
     448
     449    // Convert the input cell coordinates to FPA coordinates.
     450    fpaCoord = psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord);
     451
     452    // Convert the FPA coordinates to tangent plane Coordinates.
     453    tpCoord = psPlaneDistortApply(tpCoord, parFPA->toTangentPlane,
     454                                  fpaCoord, color, magnitude);
     455
     456    // Generate a grommit for this FPA.
     457    tmpGrommit = psGrommitAlloc(parFPA->exposure);
     458
     459    // Convert the tangent plane Coordinates to sky coordinates.
     460    skyCoord = psCoordTPToSky(skyCoord, tpCoord, tmpGrommit);
     461
     462    psFree(fpaCoord);
     463    psFree(tpCoord);
     464    psFree(tmpGrommit);
     465
     466    return(skyCoord);
     467}
     468
     469psSphere* psCoordCellToSkyQuick(psSphere* outSphere,
     470                                const psPlane* cellCoord,
     471                                const psCell* cell)
     472{
     473    if (cellCoord == NULL) {
     474        psAbort(__func__, "input parameter cellCoord is NULL.");
     475    }
     476    if (cell == NULL) {
     477        psAbort(__func__, "input parameter cell is NULL.");
     478    }
     479
     480    psPlane *tpCoord = NULL;
     481    psChip *chip = cell->parent;
     482    psFPA *FPA = chip->parent;
     483    psProjectionType oldProjectionType;
     484
     485    if (outSphere == NULL) {
     486        outSphere = (psSphere* ) psAlloc(sizeof(psSphere));
     487    }
     488
     489    // Determine the tangent plane coordinates.
     490    tpCoord = psPlaneTransformApply(tpCoord, cell->toTP, cellCoord);
     491
     492    // Save the old projection type and set the new projection type to TAN.
     493    oldProjectionType = FPA->projection->type;
     494    FPA->projection->type = PS_PROJ_TAN;
     495
     496    // Deproject the tangent plane coordinates a sphere.
     497    outSphere = psDeproject(tpCoord, FPA->projection);
     498
     499    // Restore old projection type.  Free memory.
     500    FPA->projection->type = oldProjectionType;
     501    psFree(tpCoord);
     502
     503    return (outSphere);
     504}
     505
     506/*****************************************************************************
     507XXX: What about units for the (x,y) coords?
     508 *****************************************************************************/
     509psPlane* psCoordSkyToTP(psPlane* tpCoord,
     510                        const psSphere* in,
     511                        const psGrommit* grommit)
     512{
     513    if (in == NULL) {
     514        psAbort(__func__, "input parameter in is NULL.");
     515    }
     516    if (grommit == NULL) {
     517        psAbort(__func__, "input parameter grommit is NULL.");
     518    }
     519    if (tpCoord == NULL) {
     520        tpCoord = (psPlane* ) psAlloc(sizeof(psPlane));
     521    }
     522
     523    slaOapqk("RA", in->r, in->d, (double *) grommit, &tpCoord->x, &tpCoord->y);
     524
     525    return(tpCoord);
     526}
     527
     528
     529psPlane* psCoordTPToFPA(psPlane* fpaCoord,
     530                        const psPlane* tpCoord,
     531                        double color,
     532                        double magnitude,
     533                        const psFPA* fpa)
     534{
     535    if (tpCoord == NULL) {
     536        psAbort(__func__, "input parameter tpCoord is NULL.");
     537    }
     538    if (fpa == NULL) {
     539        psAbort(__func__, "input parameter fpa is NULL.");
     540    }
     541
     542    return (psPlaneDistortApply(fpaCoord, fpa->fromTangentPlane,
     543                                tpCoord, color, magnitude));
     544}
     545
     546psPlane* psCoordFPAToChip(psPlane* chipCoord,
     547                          const psPlane* fpaCoord,
     548                          const psChip* chip)
     549{
     550    if (fpaCoord == NULL) {
     551        psAbort(__func__, "input parameter fpaCoord is NULL.");
     552    }
     553    if (chip == NULL) {
     554        psAbort(__func__, "input parameter chip is NULL.");
     555    }
     556
     557    psFPA *FPA = chip->parent;
     558    if (FPA == NULL) {
     559        psAbort(__func__, "chip->parent is NULL");
     560    }
     561
     562    // Determine which chip contains these FPA coordinates.
     563    psChip *newChip = psChipInFPA(fpaCoord, FPA);
     564    if (newChip == NULL) {
     565        return(NULL);
     566    }
     567
     568    chipCoord = psPlaneTransformApply(chipCoord, newChip->fromFPA, fpaCoord);
     569    psFree(newChip);
     570    return(chipCoord);
     571}
     572
     573psPlane* psCoordChipToCell(psPlane* cellCoord,
     574                           const psPlane* chipCoord,
     575                           const psCell* cell)
     576{
     577    if (chipCoord == NULL) {
     578        psAbort(__func__, "input parameter chipCoord is NULL.");
     579    }
     580    if (cell == NULL) {
     581        psAbort(__func__, "input parameter cell is NULL.");
     582    }
     583
     584    psChip *chip = cell->parent;
     585    if (chip == NULL) {
     586        psAbort(__func__, "cell->parent is NULL");
     587    }
     588
     589    // Determine which cell contains these FPA coordinates.
     590    psCell *newCell = psCellInChip(chipCoord, chip);
     591    if (newCell == NULL) {
     592        return(NULL);
     593    }
     594
     595    cellCoord = psPlaneTransformApply(cellCoord, newCell->fromChip, chipCoord);
     596    psFree(newCell);
     597    return(cellCoord);
     598}
     599
     600/*****************************************************************************
     601XXX: Should we do anything to determine which cell contains the coordinates?
     602 *****************************************************************************/
     603psPlane* psCoordSkyToCell(psPlane* cellCoord,
     604                          const psSphere* skyCoord,
     605                          double color,
     606                          double magnitude,
     607                          const psCell* cell)
     608{
     609    if (skyCoord == NULL) {
     610        psAbort(__func__, "input parameter skyCoord is NULL.");
     611    }
     612    if (cell == NULL) {
     613        psAbort(__func__, "input parameter cell is NULL.");
     614    }
     615
     616    psChip *parChip = cell->parent;
     617    psFPA *parFPA = parChip->parent;
     618    psGrommit* grommit = parFPA->grommit;
     619
     620    // Convert the skyCoords to tangent plane coords.
     621    psPlane *tpCoord = psCoordSkyToTP(tpCoord, skyCoord, grommit);
     622
     623    // Convert the tangent plane coords to FPA coords.
     624    psPlane *fpaCoord = psCoordTPToFPA(fpaCoord, tpCoord, color,
     625                                       magnitude, parFPA);
     626
     627    // Convert the FPA coords to chip coords.
     628    psPlane *chipCoord = psCoordFPAToChip(chipCoord, fpaCoord, parChip);
     629
     630    // Convert the chip coords to cell coords.
     631    cellCoord = psCoordChipToCell(cellCoord, chipCoord, cell);
     632
     633    psFree(tpCoord);
     634    psFree(fpaCoord);
     635    psFree(chipCoord);
     636
     637    return (cellCoord);
     638}
     639
    604640/*****************************************************************************
    605641XXX: Should we do anything to determine which cell contains the coordinates?

  • trunk/psLib/src/astronomy/psCoord.c

    r1531 r1532  
    1010*  @author George Gusciora, MHPCC
    1111*
    12 *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2004-08-13 23:33:13 $
     12*  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2004-08-13 23:43:29 $
    1414*
    1515*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    1616*/
    17 
     17/******************************************************************************/
     18/*  INCLUDE FILES                                                             */
     19/******************************************************************************/
    1820#include "psType.h"
    1921#include "psImage.h"
    … …  
    2628#include <math.h>
    2729#include <float.h>
    28 
    29 static float cot(float x);
    30 static float arg(float x, float y);
    31 
    32 psPlane* psPlaneTransformApply(psPlane* out,
    33                                const psPlaneTransform* transform,
    34                                const psPlane* coords)
    35 {
    36     if (out == NULL) {
    37         out = (psPlane* ) psAlloc(sizeof(psPlane));
    38     }
    39     out->x = transform->x->coeff[0][0] +
    40              (transform->x->coeff[1][0] * coords->x) +
    41              (transform->x->coeff[0][1] * coords->y);
    42 
    43     out->y = transform->y->coeff[0][0] +
    44              (transform->y->coeff[1][0] * coords->x) +
    45              (transform->y->coeff[0][1] * coords->y);
    46 
    47     return (out);
    48 }
    49 
    50 // This transformation takes into account parameters beyond an objects
    51 // spatial coordinates: term3 and term4 (magnitude and color).
    52 psPlane* psPlaneDistortApply(psPlane* out,
    53                              const psPlaneDistort* transform,
    54                              const psPlane* coords,
    55                              float term3,
    56                              float term4)
    57 {
    58     if (out == NULL) {
    59         out = (psPlane* ) psAlloc(sizeof(psPlane));
    60     }
    61 
    62     out->x = transform->x->coeff[0][0][0][0] +
    63              (transform->x->coeff[1][0][0][0] * coords->x) +
    64              (transform->x->coeff[0][1][0][0] * coords->y) +
    65              (transform->x->coeff[0][0][1][0] * term3) +
    66              (transform->x->coeff[0][0][0][1] * term4);
    67 
    68     out->y = transform->y->coeff[0][0][0][0] +
    69              (transform->y->coeff[1][0][0][0] * coords->x) +
    70              (transform->y->coeff[0][1][0][0] * coords->y) +
    71              (transform->y->coeff[0][0][1][0] * term3) +
    72              (transform->y->coeff[0][0][0][1] * term4);
    73 
    74     return (out);
    75 }
    76 
    77 // This function prototype has been modified since the SDRS.
    78 psSphereTransform* psSphereTransformAlloc(double NPlat,
    79         double Xo,
    80         double xo)
    81 {
    82     psSphereTransform* tmp = (psSphereTransform* ) psAlloc(sizeof(psSphereTransform));
    83 
    84     tmp->sinPhi = sin(NPlat);
    85     tmp->cosPhi = cos(NPlat);
    86     tmp->Xo = Xo;
    87     tmp->xo = xo;
    88 
    89     return (tmp);
    90 }
    91 
    92 /******************************************************************************
    93 This algorithm comes from an email from Gene.  I assume (x,y) corresponds to
    94 (r,d) in the sphere coordinates.
    95  
    96 XXX: In Gene's email, there are different variables with similar names (y
    97 and Y) and in the code, there's cos_y, cos_Y, and cos(y).  Verify that there
    98 are no typo's.
    99  *****************************************************************************/
    100 psSphere* psSphereTransformApply(psSphere* out,
    101                                  const psSphereTransform* transform,
    102                                  const psSphere* coord)
    103 {
    104     double sinY = 0.0;
    105     double cosY = 0.0;
    106     double sinX = 0.0;
    107     double cosX = 0.0;
    108     double x = 0.0;
    109     double y = 0.0;
    110     double dx = 0.0;
    111 
    112     if (out == NULL) {
    113         out = (psSphere* ) psAlloc(sizeof(psSphere));
    114     }
    115 
    116     x = coord->r;
    117     y = coord->d;
    118     dx = x - transform->xo;
    119     sinY = cos(y) * sin(dx) * transform->sinPhi + sin(y) * transform->cosPhi;
    120     cosY = sqrt(1.0 - sinY * sinY);
    121     sinX = (cos(y) * sin(dx) * transform->cosPhi - sin(y) * transform->sinPhi) / cos(y);
    122     cosX = cos(y) * cos(dx) / cos(y);
    123 
    124     out->r = atan2(sinX, cosX) + transform->Xo;
    125     out->d = atan2(sinY, cosY);
    126 
    127     return (out);
    128 }
    129 
    130 psSphereTransform* psSphereTransformICRStoEcliptic(psTime time)
    131 {
    132     struct tm *tmTime = psTimeToTM(time);
    133     double year = (double)(1900 + tmTime->tm_year);
    134     double T = year / 100.0;
    135     double phi = -23.452294 + 0.013013 * T + 0.000001639 * T * T - 0.000000503 * T * T * T;
    136     double Xo = 0.0;
    137     double xo = 0.0;
    138 
    139     return (psSphereTransformAlloc(phi, Xo, xo));
    140 }
    141 
    142 psSphereTransform* psSphereTransformEcliptictoICRS(psTime time)
    143 {
    144     struct tm *tmTime = psTimeToTM(time);
    145     double year = (double)(1900 + tmTime->tm_year);
    146     double T = year / 100.0;
    147     double phi = +23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T;
    148     double Xo = 0.0;
    149     double xo = 0.0;
    150 
    151     return (psSphereTransformAlloc(phi, Xo, xo));
    152 }
    153 
    154 psSphereTransform* psSphereTransformICRStoGalatic(void)
    155 {
    156     return (psSphereTransformAlloc(62.6, 282.25, 33.0));
    157 }
    158 
    159 psSphereTransform* psSphereTransformGalatictoICRS(void)
    160 {
    161     return (psSphereTransformAlloc(-62.6, 33.0, 282.25));
    162 }
    163 
    164 float cot(float x)
     30/******************************************************************************/
     31/*  DEFINE STATEMENTS                                                         */
     32/******************************************************************************/
     33
     34// None
     35
     36/******************************************************************************/
     37/*  TYPE DEFINITIONS                                                          */
     38/******************************************************************************/
     39
     40// None
     41
     42/*****************************************************************************/
     43/*  GLOBAL VARIABLES                                                         */
     44/*****************************************************************************/
     45
     46// None
     47
     48/*****************************************************************************/
     49/*  FILE STATIC VARIABLES                                                    */
     50/*****************************************************************************/
     51
     52// None
     53
     54/*****************************************************************************/
     55/*  FUNCTION IMPLEMENTATION - LOCAL                                          */
     56/*****************************************************************************/
     57
     58static float p_psCot(float x);
     59static float p_psArg(float x, float y);
     60
     61/******************************************************************************
     62XXX: Do this with a macro.
     63 *****************************************************************************/
     64float p_psCot(float x)
    16565{
    16666    return (1.0 / atan(x));
    … …  
    17070XXX: Verify this arc tan function.
    17171 *****************************************************************************/
    172 float arg(float x,
    173           float y)
     72float p_psArg(float x,
     73              float y)
    17474{
    17575    if (x > 0) {
    … …  
    18585    }
    18686
    187     psAbort(__func__, "Unacceptable range for (arg(%f, %f).\n", x, y);
     87    psAbort(__func__, "Unacceptable range for (p_psArg(%f, %f).\n", x, y);
    18888    return (0.0);
     89}
     90
     91/*****************************************************************************/
     92/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
     93/*****************************************************************************/
     94psPlane* psPlaneTransformApply(psPlane* out,
     95                               const psPlaneTransform* transform,
     96                               const psPlane* coords)
     97{
     98    if (out == NULL) {
     99        out = (psPlane* ) psAlloc(sizeof(psPlane));
     100    }
     101    out->x = transform->x->coeff[0][0] +
     102             (transform->x->coeff[1][0] * coords->x) +
     103             (transform->x->coeff[0][1] * coords->y);
     104
     105    out->y = transform->y->coeff[0][0] +
     106             (transform->y->coeff[1][0] * coords->x) +
     107             (transform->y->coeff[0][1] * coords->y);
     108
     109    return (out);
     110}
     111
     112/******************************************************************************
     113This transformation takes into account parameters beyond an objects spatial
     114coordinates: term3 and term4 (magnitude and color).
     115 *****************************************************************************/
     116psPlane* psPlaneDistortApply(psPlane* out,
     117                             const psPlaneDistort* transform,
     118                             const psPlane* coords,
     119                             float color,
     120                             float magnitude)
     121{
     122    if (out == NULL) {
     123        out = (psPlane* ) psAlloc(sizeof(psPlane));
     124    }
     125
     126    out->x = transform->x->coeff[0][0][0][0] +
     127             (transform->x->coeff[1][0][0][0] * coords->x) +
     128             (transform->x->coeff[0][1][0][0] * coords->y) +
     129             (transform->x->coeff[0][0][1][0] * color) +
     130             (transform->x->coeff[0][0][0][1] * magnitude);
     131
     132    out->y = transform->y->coeff[0][0][0][0] +
     133             (transform->y->coeff[1][0][0][0] * coords->x) +
     134             (transform->y->coeff[0][1][0][0] * coords->y) +
     135             (transform->y->coeff[0][0][1][0] * color) +
     136             (transform->y->coeff[0][0][0][1] * magnitude);
     137
     138    return (out);
     139}
     140
     141/******************************************************************************
     142This function prototype has been modified since the SDRS.
     143 *****************************************************************************/
     144psSphereTransform* psSphereTransformAlloc(double NPlat,
     145        double Xo,
     146        double xo)
     147{
     148    psSphereTransform* tmp = (psSphereTransform* ) psAlloc(sizeof(psSphereTransform));
     149
     150    tmp->sinPhi = sin(NPlat);
     151    tmp->cosPhi = cos(NPlat);
     152    tmp->Xo = Xo;
     153    tmp->xo = xo;
     154
     155    return (tmp);
     156}
     157
     158/******************************************************************************
     159This algorithm comes from an email from Gene.  I assume (x,y) corresponds to
     160(r,d) in the sphere coordinates.
     161 
     162XXX: In Gene's email, there are different variables with similar names (y
     163and Y) and in the code, there's cos_y, cos_Y, and cos(y).  Verify that there
     164are no typo's.
     165 *****************************************************************************/
     166psSphere* psSphereTransformApply(psSphere* out,
     167                                 const psSphereTransform* transform,
     168                                 const psSphere* coord)
     169{
     170    double sinY = 0.0;
     171    double cosY = 0.0;
     172    double sinX = 0.0;
     173    double cosX = 0.0;
     174    double x = 0.0;
     175    double y = 0.0;
     176    double dx = 0.0;
     177
     178    if (out == NULL) {
     179        out = (psSphere* ) psAlloc(sizeof(psSphere));
     180    }
     181
     182    x = coord->r;
     183    y = coord->d;
     184    dx = x - transform->xo;
     185    sinY = cos(y) * sin(dx) * transform->sinPhi + sin(y) * transform->cosPhi;
     186    cosY = sqrt(1.0 - sinY * sinY);
     187    sinX = (cos(y) * sin(dx) * transform->cosPhi - sin(y) * transform->sinPhi) / cos(y);
     188    cosX = cos(y) * cos(dx) / cos(y);
     189
     190    out->r = atan2(sinX, cosX) + transform->Xo;
     191    out->d = atan2(sinY, cosY);
     192
     193    return (out);
     194}
     195
     196psSphereTransform* psSphereTransformICRStoEcliptic(psTime time)
     197{
     198    struct tm *tmTime = psTimeToTM(time);
     199    double year = (double)(1900 + tmTime->tm_year);
     200    double T = year / 100.0;
     201    double phi = -23.452294 + 0.013013 * T + 0.000001639 * T * T - 0.000000503 * T * T * T;
     202    double Xo = 0.0;
     203    double xo = 0.0;
     204
     205    return (psSphereTransformAlloc(phi, Xo, xo));
     206}
     207
     208psSphereTransform* psSphereTransformEcliptictoICRS(psTime time)
     209{
     210    struct tm *tmTime = psTimeToTM(time);
     211    double year = (double)(1900 + tmTime->tm_year);
     212    double T = year / 100.0;
     213    double phi = +23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T;
     214    double Xo = 0.0;
     215    double xo = 0.0;
     216
     217    return (psSphereTransformAlloc(phi, Xo, xo));
     218}
     219
     220psSphereTransform* psSphereTransformICRStoGalatic(void)
     221{
     222    return (psSphereTransformAlloc(62.6, 282.25, 33.0));
     223}
     224
     225psSphereTransform* psSphereTransformGalatictoICRS(void)
     226{
     227    return (psSphereTransformAlloc(-62.6, 33.0, 282.25));
    189228}
    190229
    … …  
    201240
    202241    if (projection->type == PS_PROJ_TAN) {
    203         R = cot(coord->r) * (180.0 / M_PI);
     242        R = p_psCot(coord->r) * (180.0 / M_PI);
    204243        tmp->x = R * sin(coord->d);
    205244        tmp->y = R * cos(coord->d);
    … …  
    249288    if (projection->type == PS_PROJ_TAN) {
    250289        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
    251         tmp->d = arg(-coord->y, coord->x);
     290        tmp->d = p_psArg(-coord->y, coord->x);
    252291        tmp->r = atan(180.0 / (R * M_PI));
    253292
    254293    } else if (projection->type == PS_PROJ_SIN) {
    255294        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
    256         tmp->d = arg(-coord->y, coord->x);
     295        tmp->d = p_psArg(-coord->y, coord->x);
    257296        tmp->r = acos((R * M_PI) / 180.0);
    258297
    … …  
    271310        chu2 *= chu2;
    272311        chu = sqrt(1.0 - chu1 - chu2);
    273         tmp->d = 2.0 * arg((2.0 * chu * chu) - 1.0, (coord->x * chu * M_PI) / 360.0);
     312        tmp->d = 2.0 * p_psArg((2.0 * chu * chu) - 1.0, (coord->x * chu * M_PI) / 360.0);
    274313        tmp->r = asin((coord->y * chu * M_PI) / 180.0);
    275314
    … …  
    314353        lin = psProject(position2, &proj);
    315354        tmp = psDeproject(lin, &proj);
     355        psFree(lin);
    316356
    317357        // XXX: Do we need to convert units in tmp?
    … …  
    349389
    350390/******************************************************************************
    351 XXX: Do I need to check for unacceptable transformation parameters?  Maybe,
     391XXX: Do we need to check for unacceptable transformation parameters?  Maybe,
    352392     if the points are on the North/South Pole, etc?
    353393 
    354 XXX: Do I need to somehow scale this projection?
     394XXX: Do we need to somehow scale this projection?
    355395 
    356396XXX: I copied the algorithm from the ADD exactly.
Note: See TracChangeset for help on using the changeset viewer.