IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 27, 2005, 9:55:16 AM (21 years ago)
Author:
desonia
Message:

* empty log message *

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r4613 r4620  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.81 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-07-27 00:36:01 $
     12*  @version $Revision: 1.82 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-07-27 19:55:15 $
    1414*
    1515*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    2929#include <math.h>
    3030#include <float.h>
    31 /******************************************************************************/
    32 /*  DEFINE STATEMENTS                                                         */
    33 /******************************************************************************/
    34 
    35 // Modified Julian Day 01/01/1900 00:00:00
    36 #define MJD_1900 15021.0
    37 
    38 // Days in Julian century
    39 #define JULIAN_CENTURY 36525.0
    40 
    41 /******************************************************************************/
    42 /*  TYPE DEFINITIONS                                                          */
    43 /******************************************************************************/
    44 
    45 // None
    46 
    47 /*****************************************************************************/
    48 /*  GLOBAL VARIABLES                                                         */
    49 /*****************************************************************************/
    50 
    51 // None
    52 
    53 /*****************************************************************************/
    54 /*  FILE STATIC VARIABLES                                                    */
    55 /*****************************************************************************/
    56 
    57 // None
    58 
    59 /*****************************************************************************/
    60 /*  FUNCTION IMPLEMENTATION - LOCAL                                          */
    61 /*****************************************************************************/
    62 
    63 /*****************************************************************************/
    64 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    65 /*****************************************************************************/
     31
     32
    6633static void planeFree(psPlane *p)
    6734{
    6835    // There are non dynamic allocated items
     36}
     37
     38static void sphereFree(psSphere* s)
     39{
     40    // There are non dynamic allocated items
     41}
     42
     43static void cubeFree(psCube* c)
     44{
     45    // There are non dynamic allocated items
     46}
     47
     48static void planeTransformFree(psPlaneTransform *pt)
     49{
     50    psFree(pt->x);
     51    psFree(pt->y);
    6952}
    7053
    … …  
    195178{
    196179    psPlane *p = psAlloc(sizeof(psPlane));
    197 
    198180    psMemSetDeallocator(p, (psFreeFunc) planeFree);
     181
    199182    return(p);
    200183}
    201184
    202185
    203 static void sphereFree(psSphere *s)
    204 {
    205     // There are non dynamic allocated items
    206 }
    207 
    208186psSphere* psSphereAlloc(void)
    209187{
    210188    psSphere *s = psAlloc(sizeof(psSphere));
    211 
    212189    psMemSetDeallocator(s, (psFreeFunc) sphereFree);
     190
    213191    return(s);
    214 }
    215 
    216 static void cubeFree(psCube *c)
    217 {
    218     // There are non dynamic allocated items
    219192}
    220193
    … …  
    225198    psMemSetDeallocator(c, (psFreeFunc) cubeFree);
    226199    return(c);
    227 }
    228 
    229 static void planeTransformFree(psPlaneTransform *pt)
    230 {
    231     psFree(pt->x);
    232     psFree(pt->y);
    233200}
    234201
    … …  
    330297}
    331298
    332 /******************************************************************************
    333 alpha is LONGITUDE
    334 delta is LATITUDE
    335  
    336     alphaP: Take the target pole in the source system; calculate its LONGITUDE
    337      in the target system.  That longitude is alphaP.
    338     DeltaP: Take the target pole in the source system; calculate its LATITUDE
    339      in the target system.  That longitude is deltaP.
    340     phiP:   This is the LONGITUDE of the ascending node in the target system.
    341  *****************************************************************************/
    342 psSphereTransform* psSphereTransformAlloc(psF64 alphaP,
    343         psF64 deltaP,
    344         psF64 phiP)
    345 {
    346     psSphereTransform* tmp = (psSphereTransform* ) psAlloc(sizeof(psSphereTransform));
    347 
    348     tmp->cosDeltaP = cos(deltaP);
    349     tmp->sinDeltaP = sin(deltaP);
    350     tmp->alphaP = alphaP;
    351     tmp->phiP = phiP;
    352 
    353     return (tmp);
    354 }
    355299
    356300/******************************************************************************
    … …  
    369313    }
    370314    return(angle);
    371 }
    372 
    373 /******************************************************************************
    374 XXX: We convert Right Ascension angles to the range 0:PI.  Is that acceptable?
    375 XXX: Should we do something for Declination as well?
    376  *****************************************************************************/
    377 psSphere* psSphereTransformApply(psSphere* out,
    378                                  const psSphereTransform* transform,
    379                                  const psSphere* coord)
    380 {
    381     PS_ASSERT_PTR_NON_NULL(transform, NULL);
    382     PS_ASSERT_PTR_NON_NULL(coord, NULL);
    383 
    384     if (out == NULL) {
    385         out = (psSphere* ) psAlloc(sizeof(psSphere));
    386     }
    387 
    388     psF64 alpha = coord->r;
    389     psF64 delta = coord->d;
    390     psF64 alphaMinusAlphaP = alpha - transform->alphaP;
    391 
    392     psF64 eq55 = (sin(delta) * transform->cosDeltaP) -
    393                  (cos(delta) * transform->sinDeltaP * sin(alphaMinusAlphaP));
    394     psF64 eq56 = (cos(delta) * transform->cosDeltaP * sin(alphaMinusAlphaP)) +
    395                  (sin(delta) * transform->sinDeltaP);
    396     psF64 eq57 = cos(delta) * cos(alphaMinusAlphaP);
    397 
    398     psF64 theta = asin(eq55);
    399     psF64 phi = atan2(eq56, eq57) + transform->phiP;
    400     out->r = piNormalize(phi);
    401     out->d = theta;
    402 
    403     return(out);
    404 }
    405 
    406 psSphereTransform* psSphereTransformICRSToEcliptic(psTime *time)
    407 {
    408     psF64 T;
    409 
    410     // Check for null parameter
    411     PS_ASSERT_PTR_NON_NULL(time, NULL);
    412 
    413     // Convert psTime to MJD
    414     psF64 MJD = psTimeToMJD(time);
    415 
    416     // Check the specified MJD is greater than 1900
    417     if ( MJD < MJD_1900 ) {
    418         psError(PS_ERR_BAD_PARAMETER_TYPE,true,PS_ERRORTEXT_psCoord_INVALID_MJD);
    419         return NULL;
    420     }
    421 
    422     // Calculate number of Julian centuries since 1900
    423     T = ( MJD - MJD_1900 ) / JULIAN_CENTURY;
    424 
    425     psF64 alphaP = 0.0;
    426     psF64 deltaP = DEG_TO_RAD(23.0) +
    427                    MIN_TO_RAD(27.0) +
    428                    SEC_TO_RAD(8.26) -
    429                    (SEC_TO_RAD(46.845) * T) -
    430                    (SEC_TO_RAD(0.0059) * T * T) +
    431                    (SEC_TO_RAD(0.00181) * T * T * T);
    432     psF64 phiP = 0.0;
    433 
    434     // Don't neglect the minus sign on deltaP (bug 244):
    435     return (psSphereTransformAlloc(alphaP, deltaP, phiP));
    436 }
    437 
    438 
    439 psSphereTransform* psSphereTransformEclipticToICRS(psTime *time)
    440 {
    441     psF64 T;
    442 
    443     // Check for null parameter
    444     PS_ASSERT_PTR_NON_NULL(time, NULL);
    445 
    446     // Convert psTime to MJD
    447     psF64 MJD = psTimeToMJD(time);
    448 
    449     // Check the specified MJD is greater than 1900
    450     if ( MJD < MJD_1900 ) {
    451         psError(PS_ERR_BAD_PARAMETER_TYPE,true,PS_ERRORTEXT_psCoord_INVALID_MJD);
    452         return NULL;
    453     }
    454 
    455     // Calculate number of Julian centuries since 1900
    456     T = ( MJD - MJD_1900 ) / JULIAN_CENTURY;
    457 
    458     psF64 alphaP = 0.0;
    459     psF64 deltaP = DEG_TO_RAD(23.0) +
    460                    MIN_TO_RAD(27.0) +
    461                    SEC_TO_RAD(8.26) -
    462                    (SEC_TO_RAD(46.845) * T) -
    463                    (SEC_TO_RAD(0.0059) * T * T) +
    464                    (SEC_TO_RAD(0.00181) * T * T * T);
    465     psF64 phiP = 0.0;
    466 
    467     return (psSphereTransformAlloc(alphaP, -deltaP, phiP));
    468 }
    469 
    470 // XXX: This is bug 245: alphaP swaps with phiP from psSphereTransformGalacticToICRS()
    471 psSphereTransform* psSphereTransformGalacticToICRS(void)
    472 {
    473     psF64 alphaP = DEG_TO_RAD(32.93192);
    474     psF64 deltaP = DEG_TO_RAD(-62.87175);
    475     psF64 phiP = DEG_TO_RAD(282.85948);
    476 
    477     return (psSphereTransformAlloc(alphaP, deltaP, phiP));
    478 }
    479 
    480 psSphereTransform* psSphereTransformICRSToGalactic(void)
    481 {
    482     psF64 alphaP = DEG_TO_RAD(282.85948);
    483     psF64 deltaP = DEG_TO_RAD(62.87175);
    484     psF64 phiP = DEG_TO_RAD(32.93192);
    485 
    486     return (psSphereTransformAlloc(alphaP, deltaP, phiP));
    487315}
    488316
    … …  
    617445}
    618446
    619 /******************************************************************************
    620 The basic idea is to project both positions onto the linear plane, with
    621 position1 at the center, then calculate the linear offset between those
    622 projections.
    623  
    624 XXX: Do I need to check for unacceptable transformation parameters?  Maybe,
    625      if the points are on the North/South Pole, etc?
    626  
    627 XXX: Do I need to somehow scale this projection?
    628  
    629 XXX: Does PS_LINEAR mode make sense?  The result must be returned in psSphere
    630      regardless of the mode.
    631  
    632 XXX: How to compound errors?
    633  *****************************************************************************/
    634 psSphere* psSphereGetOffset(const psSphere* position1,
    635                             const psSphere* position2,
    636                             psSphereOffsetMode mode,
    637                             psSphereOffsetUnit unit)
    638 {
    639     PS_ASSERT_PTR_NON_NULL(position1, NULL);
    640     PS_ASSERT_PTR_NON_NULL(position2, NULL);
    641 
    642     // Check positions near 90 degree and issue warnings if necessary
    643     if (position1->d >= DEG_TO_RAD(90.0)) {
    644         psLogMsg(__func__, PS_LOG_WARN,
    645                  "WARNING: psDeproject(): position1->d is larger than 90 degrees.  Returning NULL.");
    646         return NULL;
    647     }
    648     if (position2->d >= DEG_TO_RAD(90.0)) {
    649         psLogMsg(__func__, PS_LOG_WARN,
    650                  "WARNING: psDeproject(): position2->d is larger than 90 degrees.  Returning NULL.");
    651         return NULL;
    652     }
    653 
    654     // Allocate return structure
    655     psSphere* tmp = psSphereAlloc();
    656 
    657     // Mode is LINEAR - Use first position as projection center and project second point
    658     // onto tangent plane, set point projected into psSphere structure x->r y->d
    659     if (mode == PS_LINEAR) {
    660         psProjection* proj = psProjectionAlloc(position1->r,
    661                                                position1->d,
    662                                                1.0,
    663                                                1.0,
    664                                                PS_PROJ_TAN);
    665 
    666         // Perform projection onto tangent plane
    667         psPlane* lin = psProject(position2, proj);
    668 
    669         // Set return values
    670         tmp->r = lin->x;
    671         tmp->d = lin->y;
    672 
    673         // Free data structures allocated
    674         psFree(proj);
    675         psFree(lin);
    676 
    677         // Mode is SPHERICAL - Get difference between positiion 1 and position 2 and convert
    678         // offset value from radians to desired units and return
    679     } else if (mode == PS_SPHERICAL) {
    680         tmp->r = position2->r - position1->r;
    681         tmp->d = position2->d - position1->d;
    682 
    683         // Wrap these to an acceptable range.  This assumes that all
    684         // angles are in radians.
    685         tmp->r = fmod(tmp->r, 2*M_PI);
    686         tmp->d = fmod(tmp->d, 2*M_PI);
    687         tmp->rErr = 0.0;
    688         tmp->dErr = 0.0;
    689 
    690         // Convert to desired units
    691         if (unit == PS_ARCSEC) {
    692             tmp->r = RAD_TO_SEC(tmp->r);
    693             tmp->d = RAD_TO_SEC(tmp->d);
    694         } else if (unit == PS_ARCMIN) {
    695             tmp->r = RAD_TO_MIN(tmp->r);
    696             tmp->d = RAD_TO_MIN(tmp->d);
    697         } else if (unit == PS_DEGREE) {
    698             tmp->r = RAD_TO_DEG(tmp->r);
    699             tmp->d = RAD_TO_DEG(tmp->d);
    700         } else if (unit == PS_RADIAN) {}
    701         else {
    702             psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    703                     PS_ERRORTEXT_psCoord_UNITS_UNKNOWN,
    704                     unit);
    705             psFree(tmp);
    706             return NULL;
    707         }
    708         // Invalid mode
    709     } else {
    710 
    711         psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    712                 PS_ERRORTEXT_psCoord_OFFSET_MODE_UNKNOWN,
    713                 mode);
    714         psFree(tmp);
    715         return NULL;
    716     }
    717 
    718     // Return value
    719     return tmp;
    720 }
    721 
    722 /******************************************************************************
    723 XXX: Do we need to check for unacceptable transformation parameters?  Maybe,
    724      if the points are on the North/South Pole, etc?
    725  
    726 XXX: Do we need to somehow scale this projection?
    727  
    728 XXX: I copied the algorithm from the ADD exactly.
    729  
    730 XXX: Should we compound errors?
    731  *****************************************************************************/
    732 
    733 psSphere* psSphereSetOffset(const psSphere* position,
    734                             const psSphere* offset,
    735                             psSphereOffsetMode mode,
    736                             psSphereOffsetUnit unit)
    737 {
    738     PS_ASSERT_PTR_NON_NULL(position, NULL);
    739     PS_ASSERT_PTR_NON_NULL(offset, NULL);
    740 
    741     psSphere* tmp;
    742     psF64 tmpR = 0.0;
    743     psF64 tmpD = 0.0;
    744 
    745     // If mode is linear then set position to projection center
    746     // and offset to linear coordinate then deproject to obtain
    747     // new sphere coordinate
    748     if (mode == PS_LINEAR) {
    749 
    750         // Allocate plane coordinate and set coordinate
    751         psPlane*  lin = psPlaneAlloc();
    752         lin->x = offset->r;
    753         lin->y = offset->d;
    754 
    755         // Allocate and set projection structure
    756         psProjection* proj = psProjectionAlloc(position->r,
    757                                                position->d,
    758                                                1.0,
    759                                                1.0,
    760                                                PS_PROJ_TAN);
    761 
    762         // Project tangent plane coord to spherical coord
    763         tmp = psDeproject(lin, proj);
    764 
    765         // Free data structures used
    766         psFree(proj);
    767         psFree(lin);
    768 
    769         // If mode is spherical then convert offset to radians, add the offset
    770         // to the position and wrap to 0 to 2pi
    771     } else if (mode == PS_SPHERICAL) {
    772 
    773         // Convert offset unit to radians
    774         if (unit == PS_ARCSEC) {
    775             tmpR = SEC_TO_RAD(offset->r);
    776             tmpD = SEC_TO_RAD(offset->d);
    777         } else if (unit == PS_ARCMIN) {
    778             tmpR = MIN_TO_RAD(offset->r);
    779             tmpD = MIN_TO_RAD(offset->d);
    780         } else if (unit == PS_DEGREE) {
    781             tmpR = DEG_TO_RAD(offset->r);
    782             tmpD = DEG_TO_RAD(offset->d);
    783         } else if (unit == PS_RADIAN) {
    784             tmpR = offset->r;
    785             tmpD = offset->d;
    786         } else {
    787             psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    788                     PS_ERRORTEXT_psCoord_UNITS_UNKNOWN,
    789                     unit);
    790             return NULL;
    791         }
    792 
    793         // Allocate sphere structure to return
    794         tmp = psSphereAlloc();
    795 
    796         // Add offset and wrap to 0 to 2PI if necessary
    797         tmp->r = position->r + tmpR;
    798         tmp->r = fmod(tmp->r, 2.0*M_PI);
    799         tmp->d = position->d + tmpD;
    800         tmp->d = fmod(tmp->d, 2.0*M_PI);
    801         tmp->rErr = 0.0;
    802         tmp->dErr = 0.0;
    803 
    804         // Invalid mode report error
    805     } else {
    806 
    807         psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    808                 PS_ERRORTEXT_psCoord_OFFSET_MODE_UNKNOWN,
    809                 mode);
    810         return NULL;
    811     }
    812 
    813     return tmp;
    814 }
    815 
    816 
    817 
    818 /******************************************************************************
    819 psSpherePrecess(coords, fromTime, toTime):
    820  
    821 XXX: Use static memory for tmpST.
    822  *****************************************************************************/
    823 psSphere *psSpherePrecess(psSphere *coords,
    824                           const psTime *fromTime,
    825                           const psTime *toTime)
    826 {
    827     // Check input for NULL pointers
    828     PS_ASSERT_PTR_NON_NULL(coords, NULL);
    829     PS_ASSERT_PTR_NON_NULL(fromTime, NULL);
    830     PS_ASSERT_PTR_NON_NULL(toTime, NULL);
    831 
    832     // Calculate Julian centuries
    833     psF64 fromMJD = psTimeToMJD(fromTime);
    834     psF64 toMJD = psTimeToMJD(toTime);
    835     psF64 T = (toMJD - fromMJD) / JULIAN_CENTURY;
    836 
    837     // Calculate conversion constants
    838     psF64 alphaP = DEG_TO_RAD(90.0) - ((DEG_TO_RAD(0.6406161) * T) +
    839                                        (DEG_TO_RAD(0.0000839) * T * T) +
    840                                        (DEG_TO_RAD(0.000005) * T * T * T));
    841 
    842     psF64 deltaP = (DEG_TO_RAD(0.5567530) * T) -
    843                    (DEG_TO_RAD(0.0001185) * T * T) -
    844                    (DEG_TO_RAD(0.0000116) * T * T * T);
    845 
    846     psF64 phiP = DEG_TO_RAD(90.0) + ((DEG_TO_RAD(0.6406161) * T) +
    847                                      (DEG_TO_RAD(0.0003041) * T * T) +
    848                                      (DEG_TO_RAD(0.0000051) * T * T * T));
    849 
    850     // Create transform with proper constants
    851     psSphereTransform *tmpST = psSphereTransformAlloc(alphaP, deltaP, phiP);
    852 
    853     // Apply transform to coordinates
    854     psSphere *out = psSphereTransformApply(NULL, tmpST, coords);
    855 
    856     psFree(tmpST);
    857 
    858     return(out);
    859 }
    860 
    861447/*****************************************************************************
    862448multiplyDPoly2D(trans1, trans2): Takes two 2-D polynomials as input and
Note: See TracChangeset for help on using the changeset viewer.