IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

...

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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.