IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5657


Ignore:
Timestamp:
Dec 1, 2005, 6:40:00 PM (21 years ago)
Author:
drobbin
Message:

Debugging/Test code updates. Aberration, Gravity, and CEOtoTEO mainly

Location:
trunk/psLib
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r5642 r5657  
    88 *  @author Robert Daniel DeSonia, MHPCC
    99 *
    10  *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2005-11-30 23:50:40 $
     10 *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2005-12-02 04:40:00 $
    1212 *
    1313 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    302302    // calculating the apparent angle from the actual angle and the sun position
    303303
    304     // first, calculate the angle between the sun vector and the actual vector
    305 
    306     // Moving to cartesian first:  XXX -- is this required?
    307304    psCube* sunVector = psSphereToCube(sun);
    308305    psCube* actualVector = psSphereToCube(actual);
    309     printf("\n Actual sphere = %.13g, %.13g", actual->r, actual->d);
    310     printf("\n Actual cube = %.13g, %.13g, %.13g \n", actualVector->x, actualVector->y, actualVector->z);
    311     printf("\n SUN sphere = %.13g, %.13g", sun->r, sun->d);
    312     printf("\n SUN cube = %.13g, %.13g, %.13g \n", sunVector->x, sunVector->y, sunVector->z);
    313     psSphere *new = psCubeToSphere(actualVector);
    314     printf(" New sphere = %.13g, %.13g", new->r, new->d);
    315     psFree(new);
     306
    316307    // use dot product to calculate the angle of separation
    317308    // N.B., assuming the psSphereToCube function returns a unit vector.
    … …  
    320311                    sunVector->y*actualVector->y +
    321312                    sunVector->z*actualVector->z);
    322     printf("\n Theta = %.13g ", theta);
    323313    double r0 = PS_AU * tan(theta);
    324     printf("   r0 = %.13g", r0);
    325314    double deflection = 4.0*PS_G*PS_M/(PS_C0*PS_C0*r0);
    326     printf("   deflection = %.13g \n", deflection);
     315
    327316    // make sure the deflection is not greater than 1.75 arcsec
    328317    double limit = SEC_TO_RAD(1.75);
    … …  
    352341    deflection = SEC_TO_RAD(deflection);
    353342    theta = atan(r0/PS_AU) * tan(deflection);
    354     printf("        deflection = %.13g  Theta=%.13g ", deflection, theta);
    355343    //    phi = sqrt( deflection*deflection - theta*theta );
    356344    phi = deflection * cos(asin(theta/deflection));
    357     printf("   Phi = %.13g \n", phi);
    358345    apparent->r = theta;
    359346    apparent->d = phi;
    … …  
    726713        in = psTimeConvert(in, PS_TIME_UT1);
    727714    }
    728     double T = (double)(in->sec) + (double)(in->nsec / 1e9);
    729     printf("\nThe Value of T is = %.19g\n", T);
     715    //    double T = (double)(in->sec) + (double)(in->nsec / 1e9);
     716    double T = psTimeToJD(in);
    730717    T += -2451545.0;
    731718    double theta = 2.0 * M_PI * (0.7790572732640 + 1.00273781191135448 * T);
    732     printf("\nThe Value of theta is = %.13g\n", theta);
     719    psSphereRot *out = psSphereRotAlloc(theta, 0.0, 0.0);
    733720    //    psSphereRot *out = psSphereRotInvert(theta, 0.0, 0.0);
    734     psSphereRot *out = psSphereRotInvert(theta, 0.0, 0.0);
    735721
    736722    psFree(in);

  • trunk/psLib/test/astro/tst_psEarthOrientation.c

    r5642 r5657  
    55*  @author d-Rob, MHPCC
    66*
    7 *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
    8 *  @date $Date: 2005-11-30 23:50:40 $
     7*  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
     8*  @date $Date: 2005-12-02 04:40:00 $
    99*
    1010*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    5656
    5757#define timesec 1049160600
     58#define objR DEG_TO_RAD(122.9153182445501)
     59#define objD DEG_TO_RAD(48.562968978679194)
     60static psSphere *obj = NULL;
     61static void objSetup(void);
     62
     63static void objSetup(void)
     64{
     65    if (obj == NULL) {
     66        obj = psSphereAlloc();
     67        obj->r = objR;
     68        obj->d = objD;
     69    }
     70}
    5871
    5972psS32 testAberration(void)
    … …  
    8699    double speed = sqrt(cubeDir->x*cubeDir->x + cubeDir->y*cubeDir->y + cubeDir->z*cubeDir->z);
    87100    // Speed of light in vacuum (src:NIST)
    88     double c;
    89     c = 299792458.0; /* m/s */
     101    double c = 299792458.0; /* m/s */
    90102    speed = speed / c;
    91103
    … …  
    104116
    105117    apparent = psAberration(apparent, actual, direction, speed);
    106     if (apparent == NULL)
    107         printf("\nApparent is NULL\n");
    108     else {
    109         printf("\napparent = r,d  = %.8g, %.8g\n", apparent->r, apparent->d);
    110         psFree(cubeDir);
    111         cubeDir = psSphereToCube(apparent);
    112         printf("  ><><Apparent = x,y,z = %.13g, %.13g, %.13g \n", cubeDir->x, cubeDir->y, cubeDir->z);
    113     }
     118
     119    printf("\nSphere Difference  =  r,d  =  %.13g, %.13g\n",
     120           (actual->r - apparent->r), (actual->d - apparent->d));
     121    psCube *outCube = psSphereToCube(apparent);
     122    printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
     123           outCube->x, outCube->y, outCube->z);
     124
     125    double x, y, z;
     126    x = -0.35963388069046304;
     127    y = 0.5555192509816625;
     128    z = 0.7497078321908413;
     129
     130    printf(" -- expectedCube = x,y,z = %.13g, %.13g, %.13g  -- \n\n", x, y, z);
     131
     132    if ( fabs(x - outCube->x) > DBL_EPSILON || fabs(y - outCube->y) > DBL_EPSILON ||
     133            fabs(z - outCube->z) > DBL_EPSILON ) {
     134        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     135                "psAberration returned incorrect values.\n");
     136        printf("Cube Difference  =  x,y,z  = %.13g, %.13g, %.13g \n\n",
     137               (x - outCube->x), (y - outCube->y), (z - outCube->z) );
     138        return 3;
     139    }
     140
     141    psFree(outCube);
     142
    114143    psFree(cubeDir);
    115144    psFree(apparent);
    … …  
    138167    sunCube->y = 2.5880956908748722e10;
    139168    sunCube->z = 1.1220046291457653e10;
    140     //    double sum = sqrt(sunCube->x*sunCube->x + sunCube->y*sunCube->y + sunCube->z*sunCube->z);
    141     //    sunCube->x = sunCube->x/sum;
    142     //    sunCube->y = sunCube->y/sum;
    143     //    sunCube->z = sunCube->z/sum;
    144     //    printf("\n <<<<SUNvec = x,y,z = %.13g, %.13g %.13g", sunCube->x, sunCube->y, sunCube->z);
    145169    psSphere *sun = psCubeToSphere(sunCube);
    146     printf("\n <<<<SUN = r,d = %.13g, %.13g\n", sun->r, sun->d);
    147170
    148171    empty = psGravityDeflection(apparent, empty, sun);
    … …  
    161184    apparent = psGravityDeflection(NULL, actual, sun);
    162185    psSphere *result = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
     186
    163187    printf("\nActual r,d = %.13g,%.13g    Apparent r,d = %.16g, %.16g \n",
    164188           actual->r, actual->d, result->r, result->d);
     189    printf("Sphere Difference  =  r,d  =  %.13g, %.13g\n",
     190           (actual->r - result->r), (actual->d - result->d));
    165191    psCube *outCube = psSphereToCube(result);
    166     printf(" -- outCube = x,y,z = %.13g, %.13g, %.13g  -- \n", outCube->x, outCube->y, outCube->z);
     192    printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
     193           outCube->x, outCube->y, outCube->z);
     194    psCube *outCube2 = psSphereToCube(actual);
     195    printf(" -- actualCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
     196           outCube2->x, outCube2->y, outCube2->z);
     197
     198    double x, y, z;
     199    x = -0.35961949760293604;
     200    y = 0.5555613950298085;
     201    z = 0.7496835020836093;
     202
     203    printf(" -- expectedCube = x,y,z = %.13g, %.13g, %.13g  -- \n\n", x, y, z);
     204
     205    if ( fabs(x - outCube->x) > DBL_EPSILON || fabs(y - outCube->y) > DBL_EPSILON ||
     206            fabs(z - outCube->z) > DBL_EPSILON ) {
     207        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     208                "psGravityDeflection returned incorrect values.\n");
     209        printf("Cube Difference  =  x,y,z  = %.13g, %.13g, %.13g \n",
     210               (x - outCube->x), (y - outCube->y), (z - outCube->z) );
     211        return 1;
     212    }
     213
    167214    psFree(outCube);
    168     psCube *outCube2 = psSphereToCube(apparent);
    169     printf(" -- outCube = x,y,z = %.13g, %.13g, %.13g  -- \n", outCube2->x, outCube2->y, outCube2->z);
    170215    psFree(outCube2);
    171216    psFree(sunCube);
    … …  
    416461psS32 testSphereRot_TEOtoCEO(void)
    417462{
     463    psSphereRot *rot = NULL;
     464    psTime *empty = NULL;
    418465    //    psTime *now = psTimeAlloc(PS_TIME_UT1);
    419466    //    now->sec = 1128530000;
    420467    //    now->nsec = 931154510;
    421     psTime *time = psTimeAlloc(PS_TIME_UTC);
    422     time->sec = timesec;
    423     time->nsec = 0;
     468    psTime *time = psTimeAlloc(PS_TIME_UT1);
     469    time->sec = timesec-1;
     470    time->nsec = 657017200;
    424471    time->leapsecond = false;
     472
    425473    //    psSphereRot *teoceo = psSphereRot_TEOtoCEO(now);
     474
     475    //return NULL for NULL input time
     476    rot = psSphereRot_TEOtoCEO(empty);
     477    if (rot != NULL) {
     478        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     479                "psSphereRot_TEOtoCEO failed to return NULL for NULL input time.\n");
     480        return 1;
     481    }
     482
    426483    psSphereRot *teoceo = psSphereRot_TEOtoCEO(time);
    427     printf("Output sphere rotation = %lf,%lf,%lf,%lf\n",
     484    //Make sure values match for other psTime type
     485    empty = psTimeAlloc(PS_TIME_UTC);
     486    empty->sec = timesec;
     487    empty->nsec = 0;
     488    empty->leapsecond = false;
     489    rot = psSphereRot_TEOtoCEO(empty);
     490    if (fabs(rot->q0-teoceo->q0) > DBL_EPSILON || fabs(rot->q1-teoceo->q1) > DBL_EPSILON ||
     491            fabs(rot->q2-teoceo->q2) > DBL_EPSILON || fabs(rot->q3-teoceo->q3) > DBL_EPSILON) {
     492        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     493                "psSphereRot_TEOtoCEO failed to return matching values for different time types.\n");
     494        return 2;
     495    }
     496    printf("\nOutput sphere rotation = %.13g, %.13g, %.13g, %.13g\n\n",
    428497           teoceo->q0, teoceo->q1, teoceo->q2, teoceo->q3);
    429498
     499    objSetup();
     500    psSphereRot *earthRot = psSphereRotConjugate(NULL, teoceo);
     501    psSphere *result = psSphereRotApply(NULL, earthRot, obj);
     502    //    psCube *cube = psSphereToCube(obj);
     503    psCube *cube = psSphereToCube(result);
     504
     505
     506    double x, y, z;
     507    x = 0.01698625430807123;
     508    y = -0.6616523084626379;
     509    z = 0.7496183628158023;
     510    if ( fabs(x-cube->x) > DBL_EPSILON  || fabs(y-cube->y) > DBL_EPSILON ||
     511            fabs(z-cube->z) > DBL_EPSILON) {
     512        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     513                "psSphereRot_TEOtoCEO returned incorrect values.\n");
     514        printf("\nOutput cube = x,y,z = %.13g, %.13g, %.13g\n",
     515               cube->x, cube->y, cube->z);
     516        printf("Expected cube = x,y,z = %.13g, %.13g, %.13g\n", x, y, z);
     517        printf("A difference of:   %.13g, %.13g, %.13g\n\n",
     518               (x-cube->x), (y-cube->y), (z-cube->z));
     519        return 3;
     520    }
    430521    //    psFree(now);
     522    psFree(result);
     523    psFree(obj);
     524    psFree(earthRot);
     525    psFree(cube);
    431526    psFree(time);
    432527    psFree(teoceo);

  • trunk/psLib/test/astro/verified/tst_psEarthOrientation.stderr

    r5533 r5657  
    99<DATE><TIME>|<HOST>|E|psAberration (FILE:LINENO)
    1010    Unallowable operation: direction is NULL.
    11 
    12 ---> TESTPOINT PASSED (psEarthOrientation{psAberration()} | tst_psEarthOrientation.c)
     11<DATE><TIME>|<HOST>|E|testAberration (FILE:LINENO)
     12    psGravityDeflection returned incorrect values.
     13Return value mismatch: expected 0, got 3
     14---> TESTPOINT FAILED (psEarthOrientation{psAberration()} | tst_psEarthOrientation.c)
    1315
    1416/***************************** TESTPOINT ******************************************\
    … …  
    99101---> TESTPOINT PASSED (psEarthOrientation{psSphereRRot_ITRStoTEO()} | tst_psEarthOrientation.c)
    100102
     103<DATE><TIME>|<HOST>|E|p_runTestSuite (FILE:LINENO)
     104    One or more tests failed

  • trunk/psLib/test/astro/verified/tst_psEarthOrientation.stdout

    r5642 r5657  
    11
    2 apparent = r,d  = 2.1452242, 0.84754694
    3   ><><Apparent = x,y,z = -0.3596051297757, 0.5556035294057, 0.7496591683899
     2Sphere Difference  =  r,d  =  5.284231132663e-05, 3.676277340525e-05
     3 -- resultCube = x,y,z = -0.3596051297757, 0.5556035294057, 0.7496591683899  --
     4 -- expectedCube = x,y,z = -0.3596338806905, 0.5555192509817, 0.7497078321908  --
    45
    5  <<<<SUN = r,d = 0.1745310859793, 0.0751383461293
     6Cube Difference  =  x,y,z  = -2.87509148123e-05, -8.427842401337e-05, 4.866380096979e-05
    67
    7  Actual sphere = 2.145277004504, 0.8475837032218
    8  Actual cube = -0.3596195125758, 0.5555613903456, 0.7496834983725
    9 
    10  SUN sphere = 0.1745310859793, 0.0751383461293
    11  SUN cube = 0.9820293796219, 0.1731564131522, 0.0750676637684
    12  New sphere = 2.145277004504, 0.8475837032218
    13  Theta = -0.2006809204081    r0 = -30431058247.35   deflection = -1.941413573621e-07
    14         deflection = -9.412238611831e-13  Theta=1.888856707723e-13    Phi = -9.220762225745e-13
    158
    169Actual r,d = 2.145277004504,0.8475837032218    Apparent r,d = 2.145277004504245, 0.8475837032208979
    17  -- outCube = x,y,z = -0.3596195125763, 0.5555613903461, 0.7496834983719  --
    18  -- outCube = x,y,z = 1, 1.888856707723e-13, -9.220762225745e-13  --
     10Sphere Difference  =  r,d  =  -1.887379141863e-13, 9.220402219512e-13
     11 -- resultCube = x,y,z = -0.3596195125763, 0.5555613903461, 0.7496834983719  --
     12 -- actualCube = x,y,z = -0.3596195125758, 0.5555613903456, 0.7496834983725  --
     13 -- expectedCube = x,y,z = -0.3596194976029, 0.5555613950298, 0.7496835020836  --
     14
    1915Precession Model output = x,y,s = 58931469, 4942375.1, -1.4563071e+14
    2016Expected output = x,y,s = 0.0002857175590089, 2.396873937773e-05, -1.39700664579e-08
    … …  
    3531
    3632
    37 The Value of T is = 1049160599.657008886
     33The Value of T is = 2452730.562496030238
    3834
    39 The Value of theta is = 6594672635.129
    40 Output sphere rotation = 0.000000,0.000000,0.145165,0.989407
     35The Value of theta is = 7474.398075968
     36Output sphere rotation = 0, 0, 0.9625401009478, 0.2711393628144
     37Input cube = x,y,z = -0.3596195125758, 0.5555613903456, 0.7496834983725
     38Output cube = x,y,z = 0.5967271497347, -0.286166666538, 0.7496834983725
    4139Output sphere rotation = 1.198437080989e-05,-0.0001428587808888,-6.985033157173e-09,0.9999999897239
    4240Expected sphere rotation = -1.198452240676e-05,0.0001428589335861,1.219119351891e-10,-0.9999999897238
Note: See TracChangeset for help on using the changeset viewer.