IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 2983


Ignore:
Timestamp:
Jan 13, 2005, 12:45:40 PM (22 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r2941 r2983  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-01-10 19:47:10 $
     12*  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-01-13 22:45:28 $
    1414*
    1515*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2424#include "psConstants.h"
    2525#include "psError.h"
     26#include "psLogMsg.h"
    2627#include "psAstronomyErrors.h"
    2728#include <math.h>
    … …  
    183184
    184185/******************************************************************************
     186XXX: Private Function.
     187 
     188piNormalize(): take an input angle in radians and convert it to the range 0:PI.
     189 *****************************************************************************/
     190psF32 piNormalize(psF32 angle)
     191{
     192    while (angle < FLT_EPSILON) {
     193        angle+=M_PI;
     194    }
     195
     196    while (angle >= M_PI) {
     197        angle-=M_PI;
     198    }
     199    return(angle);
     200}
     201
     202/******************************************************************************
     203XXX: We convert Right Ascension angles to the range 0:PI.  Is that acceptable?
     204XXX: Should we do something for Declination as well?
    185205 *****************************************************************************/
    186206psSphere* psSphereTransformApply(psSphere* out,
    … …  
    210230    psF64 phi = atan2(eq56, eq57) + transform->alphaP;
    211231
    212     out->r = phi;
     232
     233    out->r = piNormalize(phi);
    213234    out->d = theta;
    214235
    … …  
    281302XXX: Waiting for the definition of the PS_PROJ_PAR projection.
    282303XXX: Waiting for the definition of the PS_PROJ_GLS projection.
     304XXX: Must apply scaling at the end.
    283305 *****************************************************************************/
    284306psPlane* psProject(const psSphere* coord,
    … …  
    288310    PS_PTR_CHECK_NULL(projection, NULL);
    289311
    290     //    float R = 0.0;
    291     //    float alpha = 0.0;
    292312    psPlane* tmp = (psPlane* ) psAlloc(sizeof(psPlane));
    293 
    294313    if (projection->type == PS_PROJ_TAN) {
    295         tmp->x = cos(coord->r) * sin(coord->d) / sin(coord->r);
    296         tmp->y = -cos(coord->r) * cos(coord->d) / sin(coord->r);
     314
     315        /*
     316                // ********************************************
     317                // From the ADD
     318                // ********************************************
     319                // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work.");
     320                // return(NULL);
     321                psF32 sinTheta;
     322                psF32 cosThetaCosPhi;
     323                psF32 cosThetaSinPhi;
     324         
     325                sinTheta = (sin(coord->d) * sin(projection->D)) +
     326                           (cos(coord->d) * cos(projection->D) * (cos(coord->r - projection->R)));
     327                cosThetaCosPhi = (sin(coord->d) * cos(projection->D)) -
     328                                 (cos(coord->d) * sin(projection->D) * (cos(coord->r - projection->R)));
     329                cosThetaSinPhi = - cos(coord->d) * sin(coord->r - projection->R);
     330         
     331                tmp->x =  -cosThetaSinPhi / sinTheta;
     332                tmp->y = cosThetaCosPhi / sinTheta;
     333         
     334        */
     335
     336        // ********************************************
     337        // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html
     338        // delta_0 and phi_1 are the projection centers, not the point being projected.
     339        // (delta_0, phi_1) == (projection->R, projection->D)
     340        // ********************************************
     341        psF32 cosC = (sin(projection->D) * sin(coord->d)) +
     342                     (cos(projection->D) * cos(coord->d) * cos(coord->r - projection->R));
     343        tmp->x = (cos(coord->d) * sin(coord->r - projection->R)) / cosC;
     344        tmp->y = ((cos(projection->D) * sin(coord->d)) -
     345                  (sin(projection->D) * cos(coord->d) * cos(coord->r - projection->R))) /
     346                 cosC;
    297347
    298348    } else if (projection->type == PS_PROJ_SIN) {
     349        // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work.");
     350        // return(NULL);
     351
    299352        tmp->x = cos(coord->d) * sin(coord->d);
    300353        tmp->y = -cos(coord->d) * cos(coord->d);
    … …  
    332385
    333386/******************************************************************************
     387XXX: Waiting for algorithms.
    334388XXX: Waiting for the definition of the PS_PROJ_PAR projection.
    335389XXX: Waiting for the definition of the PS_PROJ_GLS projection.
     390XXX: Must apply scaling at the beginning.
    336391 *****************************************************************************/
    337392psSphere* psDeproject(const psPlane* coord,
    … …  
    340395    PS_PTR_CHECK_NULL(coord, NULL);
    341396    PS_PTR_CHECK_NULL(projection, NULL);
    342 
     397    //    psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): These projections don't work.");
     398    //    return(NULL);
    343399    float R = 0.0;
    344400    float chu = 0.0;
    … …  
    348404
    349405    if (projection->type == PS_PROJ_TAN) {
    350         tmp->r = atan2(-coord->y, coord->x);
    351         R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
    352         tmp->d = atan(1.0 / R);
     406        /*
     407                // ********************************************
     408                // From the ADD
     409                // ********************************************
     410                psF32 phi;
     411                psF32 theta;
     412                psF32 sinDelta;
     413                psF32 cosDeltaCosAlphaMinusAlphaP;
     414                psF32 cosDeltaSinAlphaMinusAlphaP;
     415         
     416                phi = atan2(-coord->y, coord->x);
     417                R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
     418                theta = atan(1.0 / R);
     419         
     420                sinDelta = (sin(theta) * sin(projection->D)) +
     421                           (cos(theta) * cos(projection->D) * cos(phi));
     422                cosDeltaCosAlphaMinusAlphaP =
     423                           (sin(theta) * cos(projection->D)) -
     424                           (cos(theta) * sin(projection->D) * cos(phi));
     425                cosDeltaSinAlphaMinusAlphaP = -cos(theta) * sin(phi);
     426                tmp->d = asin(sinDelta);
     427                tmp->r = atan2(cosDeltaSinAlphaMinusAlphaP, cosDeltaCosAlphaMinusAlphaP) + projection->R;
     428        */
     429        // ********************************************
     430        // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html
     431        // delta_0 and phi_1 are the projection centers, not the point being projected.
     432        // (delta_0, phi_1) == (projection->R, projection->D)
     433        // ********************************************
     434        psF32 row = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
     435        psF32 C = atan(row);
     436        tmp->d = asin((cos(C) * sin(projection->D)) + ((coord->y * sin(C) * cos(projection->D)) / row));
     437        psF32 tmpAtan = atan((coord->x * sin(C)) /
     438                             ((row * cos(projection->D) * cos(C)) - (coord->y * sin(projection->D) * sin(C))));
     439        tmp->r = projection->R + tmpAtan;
     440
    353441
    354442    } else if (projection->type == PS_PROJ_SIN) {
    … …  
    415503    PS_PTR_CHECK_NULL(position2, NULL);
    416504
     505    if (position1->d >= DEG_TO_RAD(90.0)) {
     506        psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): position1->d is larger than 90 degrees.  Returning NULL.\n");
     507        return(NULL);
     508    }
     509    if (position2->d >= DEG_TO_RAD(90.0)) {
     510        psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): position2->d is larger than 90 degrees.  Returning NULL.\n");
     511        return(NULL);
     512    }
     513
    417514    psPlane* lin;
    418515    psProjection proj;
    … …  
    426523        proj.type = PS_PROJ_TAN;
    427524
    428         printf("------------ position 1 is (%f, %f)\n", position1->r, position1->d);
    429         printf("------------ position 2 is (%f, %f)\n", position2->r, position2->d);
    430525        lin = psProject(position1, &proj);
    431         printf("------------ projected position 1 is (%f, %f)\n", lin->y, lin->x);
    432526        lin = psProject(position2, &proj);
    433         printf("------------ projected position 2 is (%f, %f)\n", lin->y, lin->x);
    434527        tmp = psDeproject(lin, &proj);
    435         printf("------------ deprojected position 2 is (%f, %f)\n", tmp->r, tmp->d);
    436 
    437528        psFree(lin);
    438 
    439529        // XXX: Do we need to convert units in tmp?
    440530        return (tmp);

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

    r2931 r2983  
    88 *  @author GLG, MHPCC
    99 *
    10  *  @version $Revision: 1.54 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2005-01-07 23:52:15 $
     10 *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2005-01-13 22:45:28 $
    1212 *
    1313 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    566566        }
    567567    }
    568 
     568    printf("Hmmm, returning NULL\n");
    569569    return (NULL);
    570570}

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

    r2941 r2983  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2005-01-10 19:47:10 $
     12*  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2005-01-13 22:45:28 $
    1414*
    1515*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2424#include "psConstants.h"
    2525#include "psError.h"
     26#include "psLogMsg.h"
    2627#include "psAstronomyErrors.h"
    2728#include <math.h>
    … …  
    183184
    184185/******************************************************************************
     186XXX: Private Function.
     187 
     188piNormalize(): take an input angle in radians and convert it to the range 0:PI.
     189 *****************************************************************************/
     190psF32 piNormalize(psF32 angle)
     191{
     192    while (angle < FLT_EPSILON) {
     193        angle+=M_PI;
     194    }
     195
     196    while (angle >= M_PI) {
     197        angle-=M_PI;
     198    }
     199    return(angle);
     200}
     201
     202/******************************************************************************
     203XXX: We convert Right Ascension angles to the range 0:PI.  Is that acceptable?
     204XXX: Should we do something for Declination as well?
    185205 *****************************************************************************/
    186206psSphere* psSphereTransformApply(psSphere* out,
    … …  
    210230    psF64 phi = atan2(eq56, eq57) + transform->alphaP;
    211231
    212     out->r = phi;
     232
     233    out->r = piNormalize(phi);
    213234    out->d = theta;
    214235
    … …  
    281302XXX: Waiting for the definition of the PS_PROJ_PAR projection.
    282303XXX: Waiting for the definition of the PS_PROJ_GLS projection.
     304XXX: Must apply scaling at the end.
    283305 *****************************************************************************/
    284306psPlane* psProject(const psSphere* coord,
    … …  
    288310    PS_PTR_CHECK_NULL(projection, NULL);
    289311
    290     //    float R = 0.0;
    291     //    float alpha = 0.0;
    292312    psPlane* tmp = (psPlane* ) psAlloc(sizeof(psPlane));
    293 
    294313    if (projection->type == PS_PROJ_TAN) {
    295         tmp->x = cos(coord->r) * sin(coord->d) / sin(coord->r);
    296         tmp->y = -cos(coord->r) * cos(coord->d) / sin(coord->r);
     314
     315        /*
     316                // ********************************************
     317                // From the ADD
     318                // ********************************************
     319                // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work.");
     320                // return(NULL);
     321                psF32 sinTheta;
     322                psF32 cosThetaCosPhi;
     323                psF32 cosThetaSinPhi;
     324         
     325                sinTheta = (sin(coord->d) * sin(projection->D)) +
     326                           (cos(coord->d) * cos(projection->D) * (cos(coord->r - projection->R)));
     327                cosThetaCosPhi = (sin(coord->d) * cos(projection->D)) -
     328                                 (cos(coord->d) * sin(projection->D) * (cos(coord->r - projection->R)));
     329                cosThetaSinPhi = - cos(coord->d) * sin(coord->r - projection->R);
     330         
     331                tmp->x =  -cosThetaSinPhi / sinTheta;
     332                tmp->y = cosThetaCosPhi / sinTheta;
     333         
     334        */
     335
     336        // ********************************************
     337        // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html
     338        // delta_0 and phi_1 are the projection centers, not the point being projected.
     339        // (delta_0, phi_1) == (projection->R, projection->D)
     340        // ********************************************
     341        psF32 cosC = (sin(projection->D) * sin(coord->d)) +
     342                     (cos(projection->D) * cos(coord->d) * cos(coord->r - projection->R));
     343        tmp->x = (cos(coord->d) * sin(coord->r - projection->R)) / cosC;
     344        tmp->y = ((cos(projection->D) * sin(coord->d)) -
     345                  (sin(projection->D) * cos(coord->d) * cos(coord->r - projection->R))) /
     346                 cosC;
    297347
    298348    } else if (projection->type == PS_PROJ_SIN) {
     349        // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work.");
     350        // return(NULL);
     351
    299352        tmp->x = cos(coord->d) * sin(coord->d);
    300353        tmp->y = -cos(coord->d) * cos(coord->d);
    … …  
    332385
    333386/******************************************************************************
     387XXX: Waiting for algorithms.
    334388XXX: Waiting for the definition of the PS_PROJ_PAR projection.
    335389XXX: Waiting for the definition of the PS_PROJ_GLS projection.
     390XXX: Must apply scaling at the beginning.
    336391 *****************************************************************************/
    337392psSphere* psDeproject(const psPlane* coord,
    … …  
    340395    PS_PTR_CHECK_NULL(coord, NULL);
    341396    PS_PTR_CHECK_NULL(projection, NULL);
    342 
     397    //    psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): These projections don't work.");
     398    //    return(NULL);
    343399    float R = 0.0;
    344400    float chu = 0.0;
    … …  
    348404
    349405    if (projection->type == PS_PROJ_TAN) {
    350         tmp->r = atan2(-coord->y, coord->x);
    351         R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
    352         tmp->d = atan(1.0 / R);
     406        /*
     407                // ********************************************
     408                // From the ADD
     409                // ********************************************
     410                psF32 phi;
     411                psF32 theta;
     412                psF32 sinDelta;
     413                psF32 cosDeltaCosAlphaMinusAlphaP;
     414                psF32 cosDeltaSinAlphaMinusAlphaP;
     415         
     416                phi = atan2(-coord->y, coord->x);
     417                R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
     418                theta = atan(1.0 / R);
     419         
     420                sinDelta = (sin(theta) * sin(projection->D)) +
     421                           (cos(theta) * cos(projection->D) * cos(phi));
     422                cosDeltaCosAlphaMinusAlphaP =
     423                           (sin(theta) * cos(projection->D)) -
     424                           (cos(theta) * sin(projection->D) * cos(phi));
     425                cosDeltaSinAlphaMinusAlphaP = -cos(theta) * sin(phi);
     426                tmp->d = asin(sinDelta);
     427                tmp->r = atan2(cosDeltaSinAlphaMinusAlphaP, cosDeltaCosAlphaMinusAlphaP) + projection->R;
     428        */
     429        // ********************************************
     430        // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html
     431        // delta_0 and phi_1 are the projection centers, not the point being projected.
     432        // (delta_0, phi_1) == (projection->R, projection->D)
     433        // ********************************************
     434        psF32 row = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
     435        psF32 C = atan(row);
     436        tmp->d = asin((cos(C) * sin(projection->D)) + ((coord->y * sin(C) * cos(projection->D)) / row));
     437        psF32 tmpAtan = atan((coord->x * sin(C)) /
     438                             ((row * cos(projection->D) * cos(C)) - (coord->y * sin(projection->D) * sin(C))));
     439        tmp->r = projection->R + tmpAtan;
     440
    353441
    354442    } else if (projection->type == PS_PROJ_SIN) {
    … …  
    415503    PS_PTR_CHECK_NULL(position2, NULL);
    416504
     505    if (position1->d >= DEG_TO_RAD(90.0)) {
     506        psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): position1->d is larger than 90 degrees.  Returning NULL.\n");
     507        return(NULL);
     508    }
     509    if (position2->d >= DEG_TO_RAD(90.0)) {
     510        psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): position2->d is larger than 90 degrees.  Returning NULL.\n");
     511        return(NULL);
     512    }
     513
    417514    psPlane* lin;
    418515    psProjection proj;
    … …  
    426523        proj.type = PS_PROJ_TAN;
    427524
    428         printf("------------ position 1 is (%f, %f)\n", position1->r, position1->d);
    429         printf("------------ position 2 is (%f, %f)\n", position2->r, position2->d);
    430525        lin = psProject(position1, &proj);
    431         printf("------------ projected position 1 is (%f, %f)\n", lin->y, lin->x);
    432526        lin = psProject(position2, &proj);
    433         printf("------------ projected position 2 is (%f, %f)\n", lin->y, lin->x);
    434527        tmp = psDeproject(lin, &proj);
    435         printf("------------ deprojected position 2 is (%f, %f)\n", tmp->r, tmp->d);
    436 
    437528        psFree(lin);
    438 
    439529        // XXX: Do we need to convert units in tmp?
    440530        return (tmp);

  • trunk/psLib/test/astronomy/tst_psAstrometry01.c

    r2537 r2983  
    55*  @author GLG, MHPCC
    66*
    7 *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
    8 *  @date $Date: 2004-11-30 21:56:20 $
     7*  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
     8*  @date $Date: 2005-01-13 22:45:40 $
    99*
    1010*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2121                              {test1, 0000, "Tests psFixedPatternAlloc()", 0, false},
    2222                              {test2, 0000, "Tests psGrommitAlloc()", 0, false},
    23                               {test3, 0000, "psFunctionBar", 0, false},
     23                              {test3, 0000, "MISC", 0, false},
    2424                              {NULL}
    2525                          };
    … …  
    189189    psLogSetLevel( PS_LOG_INFO );
    190190
    191     return ( ! runTestSuite( stderr, "psImage", tests, argc, argv ) );
     191    return ( ! runTestSuite( stderr, "psAstrometry", tests, argc, argv ) );
    192192}
    193193
    … …  
    440440
    441441                    if ((myCell == NULL) || (tmpChip == NULL)) {
    442                         printf("ERROR: NULL\n");
     442                        if (myCell == NULL) {
     443                            printf("ERROR: myCell is NULL\n");
     444                        }
     445                        if (tmpChip == NULL) {
     446                            printf("ERROR: tmpChip is NULL\n");
     447                        }
    443448                    } else {
    444449                        psCoordFPAToChip(&chipCoord, &fpaCoord, tmpChip);
    … …  
    517522    }
    518523
     524    printf("HERE 00\n");
    519525    psFree(myFPA);
     526    printf("HERE 01\n");
    520527
    521528    return(0);

  • trunk/psLib/test/astronomy/tst_psCoord.c

    r2932 r2983  
    66*  @author GLG, MHPCC
    77*
    8 *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
    9 *  @date $Date: 2005-01-08 00:12:44 $
     8*  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
     9*  @date $Date: 2005-01-13 22:45:40 $
    1010*
    1111*  XXX: Must test psSpherePrecess.
    … …  
    1818#include "psTest.h"
    1919#include "pslib.h"
     20/*
    2021static psS32 test1( void );
    2122static psS32 test1b( void );
    … …  
    2425static psS32 test3( void );
    2526static psS32 test4( void );
     27static psS32 test4b( void );
    2628static psS32 test5( void );
    2729static psS32 test6( void );
    … …  
    3234static psS32 test40( void );
    3335static psS32 test41( void );
     36static psS32 test43( void );
     37*/
     38static psS32 test42( void );
    3439testDescription tests[] = {
    35                               {test1, 0000, "psSphereTransformAlloc()", 0, false},
    36                               {test1b, 0000, "psPlaneTransformAlloc()", 0, false},
    37                               {test1c, 0000, "psPlaneDistortAlloc()", 0, false},
    38                               {test2, 0000, "psPlaneTransformApply()", 0, false},
    39                               {test3, 0000, "psPlaneDistortApply()", 0, false},
    40                               {test4, 0000, "psPSphereTransformApply()", 0, false},
    41                               {test5, 0000, "psSphereTransformICRSToEcliptic()", 0, false},
    42                               {test6, 0000, "psSphereTransformEclipticToICRS()", 0, false},
    43                               {test7, 0000, "psSphereTransformICRSToGalatic()", 0, false},
    44                               {test8, 0000, "psSphereTransformGalaticToICRS()", 0, false},
    45                               {test20, 0000, "psProject()", 0, false},
    46                               {test21, 0000, "psDeProject()", 0, false},
    47                               {test40, 0000, "psSphereGetOffset()", 0, false},
    48                               {test41, 0000, "psSphereSetOffset()", 0, false},
     40                              /*
     41                                                            {test1, 0000, "psSphereTransformAlloc()", 0, false},
     42                                                            {test1b, 0000, "psPlaneTransformAlloc()", 0, false},
     43                                                            {test1c, 0000, "psPlaneDistortAlloc()", 0, false},
     44                                                            {test2, 0000, "psPlaneTransformApply()", 0, false},
     45                                                            {test3, 0000, "psPlaneDistortApply()", 0, false},
     46                                                            {test4, 0000, "psPSphereTransformApply()", 0, false},
     47                                                            {test4b, 0000, "psPSphereTransformApply()", 0, false},
     48                                                            {test5, 0000, "psSphereTransformICRSToEcliptic()", 0, false},
     49                                                            {test6, 0000, "psSphereTransformEclipticToICRS()", 0, false},
     50                                                            {test7, 0000, "psSphereTransformICRSToGalatic()", 0, false},
     51                                                            {test8, 0000, "psSphereTransformGalaticToICRS()", 0, false},
     52                                                            {test20, 0000, "psProject()", 0, false},
     53                                                            {test21, 0000, "psDeProject()", 0, false},
     54                                                            {test40, 0000, "psSphereGetOffset()", 0, false},
     55                                                            {test41, 0000, "psSphereSetOffset()", 0, false},
     56                                                            {test43, 0000, "psProject(), psDeproject", 0, false},
     57                              */
     58                              {test42, 0000, "psProject(), psDeproject", 0, false},
    4959                              {NULL}
    5060                          };
    … …  
    436446    psSphereTransform *myST = psSphereTransformAlloc(0.0, 0.0, 0.0);
    437447
    438     in.r = 45.0 * MY_DEG_TO_RAD;
    439     in.d = 30.0 * MY_DEG_TO_RAD;
     448    in.r = DEG_TO_RAD(45.0);
     449    in.d = DEG_TO_RAD(30.0);
    440450    in.rErr = 0.0;
    441451    in.dErr = 0.0;
    … …  
    443453    for (float r=0.0;r<180.0;r+=DEG_INC) {
    444454        for (float d=0.0;d<90.0;d+=DEG_INC) {
    445             in.r = r * MY_DEG_TO_RAD;
    446             in.d = d * MY_DEG_TO_RAD;
     455            in.r = DEG_TO_RAD(r);
     456            in.d = DEG_TO_RAD(d);
    447457            in.rErr = 0.0;
    448458            in.dErr = 0.0;
    … …  
    483493    printf("-------------------------------------------------------------------\n");
    484494    psFree(myST);
     495    return(testStatus);
     496}
     497
     498/******************************************************************************
     499test4b(): This test verifies that psSphereTransformApply() works properly.  We
     500create two psSphereTransforms: a forward transform and a reverse transform
     501(which is the mathematical inverse of the forward transform).  We apply both
     502transforms to several spherical coordinates and ensure that the original input
     503coordinate is obtained after applying both transforms.
     504 
     505XXX: We currently test the alpha and delta offsets independently.  Attempts to
     506test them both concurrently failed.  Determine why this is.  Are the following
     507spherical transforms not mathematical inverses?
     508    psSphereTransformAlloc(X, Y, 0.0)
     509    psSphereTransformAlloc(-X, -Y, 0.0)
     510 *****************************************************************************/
     511#define ERROR_PERCENT 0.01
     512psS32 test4b( void )
     513{
     514    psS32 testStatus = 0;
     515    psSphere in;
     516    psSphere out;
     517    psSphere out2;
     518    psSphereTransform *mySphereTransformForward = NULL;
     519    psSphereTransform *mySphereTransformReverse = NULL;
     520
     521
     522    mySphereTransformForward = psSphereTransformAlloc(DEG_TO_RAD(22.0),
     523                               0.0,
     524                               0.0);
     525    mySphereTransformReverse = psSphereTransformAlloc(DEG_TO_RAD(-22.0),
     526                               0.0,
     527                               0.0);
     528
     529    for (float r=0.1;r<180.0;r+=(DEG_INC/5.0)) {
     530        for (float d=0.1;d<90.0;d+=(DEG_INC/5.0)) {
     531            in.r = DEG_TO_RAD(r);
     532            in.d = DEG_TO_RAD(d);
     533            in.rErr = 0.0;
     534            in.dErr = 0.0;
     535
     536            psSphereTransformApply(&out, mySphereTransformForward, &in);
     537            psSphereTransformApply(&out2, mySphereTransformReverse, &out);
     538
     539            if ((fabs((in.r - out2.r) / in.r) > ERROR_PERCENT) ||
     540                    (fabs((in.d - out2.d) / in.d) > ERROR_PERCENT)) {
     541                printf("ERROR: \n");
     542                printf("Input  coords (R, D) are (%f, %f)\n", in.r, in.d);
     543                printf("Output coords (R, D) are (%f, %f)\n", out2.r, out2.d);
     544                testStatus = 4;
     545            }
     546        }
     547    }
     548    psFree(mySphereTransformForward);
     549    psFree(mySphereTransformReverse);
     550
     551    mySphereTransformForward = psSphereTransformAlloc(0.0,
     552                               DEG_TO_RAD(33.0),
     553                               0.0);
     554    mySphereTransformReverse = psSphereTransformAlloc(0.0,
     555                               DEG_TO_RAD(-33.0),
     556                               0.0);
     557    for (float r=0.1;r<180.0;r+=(DEG_INC/5.0)) {
     558        for (float d=0.1;d<90.0;d+=(DEG_INC/5.0)) {
     559            in.r = DEG_TO_RAD(r);
     560            in.d = DEG_TO_RAD(d);
     561            in.rErr = 0.0;
     562            in.dErr = 0.0;
     563
     564            psSphereTransformApply(&out, mySphereTransformForward, &in);
     565            psSphereTransformApply(&out2, mySphereTransformReverse, &out);
     566
     567            if ((fabs((in.r - out2.r) / in.r) > ERROR_PERCENT) ||
     568                    (fabs((in.d - out2.d) / in.d) > ERROR_PERCENT)) {
     569                printf("ERROR: \n");
     570                printf("Input  coords (R, D) are (%f, %f)\n", in.r, in.d);
     571                printf("Output coords (R, D) are (%f, %f)\n", out2.r, out2.d);
     572                testStatus = 4;
     573            }
     574        }
     575    }
     576    psFree(mySphereTransformForward);
     577    psFree(mySphereTransformReverse);
     578
    485579    return(testStatus);
    486580}
    … …  
    677771    psPlane *out;
    678772    psProjection myProjection;
    679     myProjection.R = 20.0 * MY_DEG_TO_RAD;
    680     myProjection.D = 10.0 * MY_DEG_TO_RAD;
     773    myProjection.R = DEG_TO_RAD(20.0);
     774    myProjection.D = DEG_TO_RAD(10.0);
    681775    myProjection.Xs = 1.0;
    682776    myProjection.Ys = 1.0;
    … …  
    684778    for (float r=0.0;r<180.0;r+=DEG_INC) {
    685779        for (float d=0.0;d<90.0;d+=DEG_INC) {
    686             in.r = r * MY_DEG_TO_RAD;
    687             in.d = d * MY_DEG_TO_RAD;
     780            in.r = DEG_TO_RAD(r);
     781            in.d = DEG_TO_RAD(d);
    688782            in.rErr = 0.0;
    689783            in.dErr = 0.0;
    … …  
    789883    psSphere *out;
    790884    psProjection myProjection;
    791     myProjection.R = 20.0 * MY_DEG_TO_RAD;
    792     myProjection.D = 10.0 * MY_DEG_TO_RAD;
     885    myProjection.R = DEG_TO_RAD(20.0);
     886    myProjection.D = DEG_TO_RAD(10.0);
    793887    myProjection.Xs = 1.0;
    794888    myProjection.Ys = 1.0;
    … …  
    904998    psSphere *offset = NULL;
    905999
    906     position1.r = 90.0 * MY_DEG_TO_RAD;
    907     position1.d = 45.0 * MY_DEG_TO_RAD;
     1000    position1.r = DEG_TO_RAD(90.0);
     1001    position1.d = DEG_TO_RAD(45.0);
    9081002    position1.rErr = 0.0;
    9091003    position1.dErr = 0.0;
    … …  
    9111005    for (float r=0.0;r<180.0;r+=DEG_INC) {
    9121006        for (float d=0.0;d<90.0;d+=DEG_INC) {
    913             position2.r = r * MY_DEG_TO_RAD;
    914             position2.d = d * MY_DEG_TO_RAD;
     1007            position2.r = DEG_TO_RAD(r);
     1008            position2.d = DEG_TO_RAD(d);
    9151009            position2.rErr = 0.0;
    9161010            position2.dErr = 0.0;
    … …  
    9241018            offset = psSphereGetOffset(&position1, &position2,
    9251019                                       PS_SPHERICAL, PS_DEGREE);
    926             offset->r = offset->r * MY_DEG_TO_RAD;
    927             offset->d = offset->d * MY_DEG_TO_RAD;
     1020            offset->r = DEG_TO_RAD(offset->r);
     1021            offset->d = DEG_TO_RAD(offset->d);
    9281022            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 2);
    9291023            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 3);
    … …  
    9461040            psFree(offset);
    9471041
    948             /* XXX: This code does not work correctly.
    949                         offset = psSphereGetOffset(&position1, &position2,
    950                                                    PS_LINEAR, 0);
    951                         printf("--------------- (%f, %f) to (%f, %f) is (%f, %f) ---------------\n",
    952                                 position1.r, position1.d, position2.r, position2.d,
    953                                 offset->r, offset->d);
     1042            //HEY            /* XXX: This code does not work correctly.
     1043            printf("--------------- (%f, %f) to (%f, %f) is (%f, %f) ---------------\n",
     1044                   position1.r, position1.d, position2.r, position2.d,
     1045                   offset->r, offset->d);
     1046            /*
     1047                        offset = psSphereGetOffset(&position1, &position2, PS_LINEAR, 0);
    9541048                        PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 1);
    9551049                        PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 1);
    956                         psFree(offset);
    9571050            */
     1051            psFree(offset);
     1052            //            */
    9581053        }
    9591054    }
    … …  
    10071102    psSphere tmpOffset;
    10081103
    1009     position1.r = 90.0 * MY_DEG_TO_RAD;
    1010     position1.d = 45.0 * MY_DEG_TO_RAD;
     1104    position1.r = DEG_TO_RAD(90.0);
     1105    position1.d = DEG_TO_RAD(45.0);
    10111106    position1.rErr = 0.0;
    10121107    position1.dErr = 0.0;
    … …  
    10141109    for (float r=0.0;r<180.0;r+=DEG_INC) {
    10151110        for (float d=0.0;d<90.0;d+=DEG_INC) {
    1016             offset.r = r * MY_DEG_TO_RAD;
    1017             offset.d = d * MY_DEG_TO_RAD;
     1111            offset.r = DEG_TO_RAD(offset.r);
     1112            offset.d = DEG_TO_RAD(offset.d);
    10181113            offset.rErr = 0.0;
    10191114            offset.dErr = 0.0;
    … …  
    10251120            psFree(position2);
    10261121
    1027             tmpOffset.r = offset.r * MY_RAD_TO_DEG;
    1028             tmpOffset.d = offset.d * MY_RAD_TO_DEG;
     1122            tmpOffset.r = DEG_TO_RAD(offset.r);
     1123            tmpOffset.d = DEG_TO_RAD(offset.d);
    10291124            tmpOffset.rErr = 0.0;
    10301125            tmpOffset.dErr = 0.0;
    … …  
    11081203}
    11091204
     1205/******************************************************************************
     1206test43(): Not really a test.  We simply project many points and print the
     1207results.
     1208 *****************************************************************************/
     1209psS32 test43( void )
     1210{
     1211    psS32 testStatus = 0;
     1212    psSphere projCenter;
     1213    psSphere startCoords;
     1214    psPlane *planeCoords;
     1215
     1216    projCenter.r = DEG_TO_RAD(33.0);
     1217    projCenter.d = DEG_TO_RAD(22.0);
     1218    projCenter.r = 0.0;
     1219    projCenter.d = 0.0;
     1220    // Create the psProjection plane:
     1221    psProjection myProj;
     1222    myProj.R = projCenter.r;
     1223    myProj.D = projCenter.d;
     1224    myProj.Xs = 1.0;
     1225    myProj.Ys = 1.0;
     1226    myProj.type = PS_PROJ_TAN;
     1227
     1228    printf("    projCenter is (%f, %f)\n", projCenter.r, projCenter.d);
     1229    for (float r=-90.0;r<90.0;r+=DEG_INC/5.0) {
     1230        for (float d=-90.0;d<90.0;d+=DEG_INC/5.0) {
     1231            printf("----------------------------------------------------------\n");
     1232            printf("Angles (R, D) are (%.2f, %.2f)\n", r, d);
     1233            startCoords.r = DEG_TO_RAD(r);
     1234            startCoords.d = DEG_TO_RAD(d);
     1235
     1236            planeCoords = psProject(&startCoords, &myProj);
     1237            printf("    startCoords (R, D) is (%f, %f)\n", startCoords.r, startCoords.d);
     1238            printf("        plane Coords (X, Y) is (%f, %f)\n", planeCoords->x, planeCoords->y);
     1239            psFree(planeCoords);
     1240        }
     1241    }
     1242    printf("----------------------------------------------------------\n");
     1243    return(testStatus);
     1244}
     1245
     1246#define NUM_DEGREES 5
     1247/******************************************************************************
     1248test42(): This test attempts to verify the psProject() and psDeproject()
     1249functions.  We create an arbitrary psProjection.  We then project several
     1250spherical coordinates onto the plane and then immediately deproject them; we
     1251ensure that the resulting spherical coordinates are identical to the start
     1252coordinates.
     1253 *****************************************************************************/
     1254psS32 test42( void )
     1255{
     1256    psS32 testStatus = 0;
     1257    psSphere projCenter;
     1258    psSphere startCoords;
     1259    psSphere *endCoords;
     1260    psPlane *planeCoords;
     1261
     1262    #define PROG_CENTER_R_DEG 20.0
     1263    #define PROG_CENTER_D_DEG 20.0
     1264
     1265    projCenter.r = DEG_TO_RAD(PROG_CENTER_R_DEG);
     1266    projCenter.d = DEG_TO_RAD(PROG_CENTER_D_DEG);
     1267    // Create the psProjection plane:
     1268    psProjection myProj;
     1269    myProj.R = projCenter.r;
     1270    myProj.D = projCenter.d;
     1271    myProj.Xs = 1.0;
     1272    myProj.Ys = 1.0;
     1273    myProj.type = PS_PROJ_TAN;
     1274
     1275    #define RADIUS_DEG 30.0
     1276    #define NUM_ANGLES 20
     1277    #define MY_DEG_INC ((2.0 * RADIUS_DEG) / NUM_ANGLES)
     1278
     1279    for (float r= (PROG_CENTER_R_DEG - RADIUS_DEG);
     1280            r<=(PROG_CENTER_R_DEG + RADIUS_DEG);
     1281            r+=MY_DEG_INC) {
     1282        printf("------------------------ r is %.1f ------------------------\n", r);
     1283        for (float d= (PROG_CENTER_D_DEG - RADIUS_DEG);
     1284                d<=(PROG_CENTER_D_DEG + RADIUS_DEG);
     1285                d+=MY_DEG_INC) {
     1286            startCoords.r = DEG_TO_RAD(r);
     1287            startCoords.d = DEG_TO_RAD(d);
     1288            printf("----- startCoords (R, D) is (%f, %f) is deg (%.1f, %.1f)  -----\n", startCoords.r, startCoords.d, r, d);
     1289
     1290            planeCoords = psProject(&startCoords, &myProj);
     1291            endCoords = psDeproject(planeCoords, &myProj);
     1292
     1293            if ((fabs((endCoords->r - startCoords.r) / startCoords.r) > ERROR_PERCENT) ||
     1294                    (fabs((endCoords->d - startCoords.d) / startCoords.d) > ERROR_PERCENT)) {
     1295                printf("ERROR: \n");
     1296                printf("    startCoords (R, D) is (%f, %f)\n", startCoords.r, startCoords.d);
     1297                //                printf("        plane Coords is (%f, %f)\n", planeCoords->x, planeCoords->y);
     1298                printf("    endCoords (R, D) is (%f, %f)\n", endCoords->r, endCoords->d);
     1299            }
     1300
     1301            psFree(planeCoords);
     1302            psFree(endCoords);
     1303        }
     1304    }
     1305    return(testStatus);
     1306}
     1307
    11101308// This code will ...
Note: See TracChangeset for help on using the changeset viewer.