IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 6030


Ignore:
Timestamp:
Jan 17, 2006, 2:41:29 PM (21 years ago)
Author:
drobbin
Message:

Updated Table parser, Fixed misc EOC, worked primarily on GravityDeflection

Location:
trunk/psLib
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/psTableParse.c

    r5970 r6030  
    77 *  @author Dave Robbins, MHPCC
    88 *
    9  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2006-01-12 21:18:48 $
     9 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2006-01-18 00:41:29 $
    1111 *
    1212 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    4444    input = fopen(argv[1], "r");
    4545    output = fopen(argv[2], "w");
     46    if (input == NULL) {
     47        printf("\nERROR.  Could not open input file.\n");
     48        return 0;
     49    }
     50    if (output == NULL) {
     51        printf("\nERROR.  Could not open or create output file.\n");
     52        return 0;
     53    }
    4654    char data2[200];
    4755    printHeader(argv[1]);

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

    r5969 r6030  
    88 *  @author Robert Daniel DeSonia, MHPCC
    99 *
    10  *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2006-01-12 20:40:12 $
     10 *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2006-01-18 00:41:29 $
    1212 *
    1313 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    3939#define PS_AU 149597890000.0 /* meters */
    4040// Modified Julian Day 1/1/2000 00:00:00
    41 #define MJD_2000  51544.0
     41#define MJD_2000  51544.5
    4242// Days in Julian century
    4343#define JULIAN_CENTURY 36525.0
    … …  
    299299    double a = 0.0;
    300300
    301     //mu = apparent * direction;  I believe this should be:
    302     //    mu = apparent->r * direction->r + apparent->d * direction->d;
    303301    psCube* directionVector = psSphereToCube(direction);
    304302    if (directionVector == NULL)
    305303        printf("directionVector is null\n");
    306304    psCube* actualVector = psSphereToCube(actual);
    307     if (directionVector == NULL)
     305    if (actualVector == NULL)
    308306        printf("actualVector is null\n");
    309307    //    mu = acos(directionVector->x*actualVector->x +
    … …  
    311309          directionVector->y*actualVector->y +
    312310          directionVector->z*actualVector->z);
    313 
    314     //    rp = apparent - mu * direction;
    315     //    rp->r = actual->r - mu * direction->r;
    316     //    rp->d = actual->d - mu * direction->d;
    317311    rp->x = actualVector->x - mu * directionVector->x;
    318312    rp->y = actualVector->y - mu * directionVector->y;
    319313    rp->z = actualVector->z - mu * directionVector->z;
    320314
    321     mu_p = mu + speed * ((mu * mu - 1.0) / (1.0 - speed * mu));
     315    mu_p = mu - speed * ((mu * mu - 1.0) / (1.0 - speed * mu));
    322316
    323317    //Not sure if this is right.  ADD gets a scalar from division of rp (a vector)
    … …  
    345339    //XXX: Must be a sign error somewhere above?  Magnitude of change is correct but wrong way...
    346340    //This will fix the problem but is somewhat of a hack.
    347     r_p->x = actualVector->x - (r_p->x - actualVector->x);
    348     r_p->y = actualVector->y - (r_p->y - actualVector->y);
    349     r_p->z = actualVector->z - (r_p->z - actualVector->z);
     341    //    r_p->x = actualVector->x - (r_p->x - actualVector->x);
     342    //    r_p->y = actualVector->y - (r_p->y - actualVector->y);
     343    //    r_p->z = actualVector->z - (r_p->z - actualVector->z);
    350344
    351345    psSphere *r_pSphere = psCubeToSphere(r_p);
    … …  
    370364    PS_ASSERT_PTR_NON_NULL(sun, NULL);
    371365
     366    psSphere *temp = psSphereAlloc();
    372367    // calculating the apparent angle from the actual angle and the sun position
    373368
    … …  
    378373    // N.B., assuming the psSphereToCube function returns a unit vector.
    379374    //    double theta = acos(sunVector->x*actualVector->x +
    380     double theta = (sunVector->x*actualVector->x +
    381                     sunVector->y*actualVector->y +
    382                     sunVector->z*actualVector->z);
    383     printf(" Theta = %lf\n", theta);
     375    double dotProd = (sunVector->x*actualVector->x +
     376                      sunVector->y*actualVector->y + sunVector->z*actualVector->z);
     377    double theta, sunMag, actMag;
     378    sunMag = sqrt(sunVector->x*sunVector->x + sunVector->y*sunVector->y +
     379                  sunVector->z*sunVector->z);
     380    actMag = sqrt(actualVector->x*actualVector->x + actualVector->y*actualVector->y +
     381                  actualVector->z*actualVector->z);
     382    dotProd = dotProd / (sunMag * actMag);
     383    theta = acos(dotProd);
     384
     385    //    theta = acos( cos(sun->d)*cos(actual->d)*cos(sun->r-actual->r) + sin(sun->r)*sin(actual->r) );
     386
     387    printf(" Theta = %.13g\n", theta);
    384388    //    theta = acos(-theta);
    385389    //    printf("\n Theta = %lf\n", theta);
    386390    double r0 = PS_AU * tan(theta);
     391    printf(" r0 = %.19e\n", r0);
    387392    double deflection = 4.0*PS_G*PS_M/(PS_C0*PS_C0*r0);
    388393
    … …  
    390395    double limit = SEC_TO_RAD(1.75);
    391396    printf(" deflection = %.13g\n", deflection);
    392     printf(" limit = %lf\n", limit);
     397    //printf(" limit = %lf\n", limit);
    393398    if (deflection > limit) {
    394399        //       deflection = limit;
    … …  
    400405    }
    401406
    402     if (apparent == NULL) {
    403         apparent = psSphereAlloc();
    404     } else {
    405         apparent->r = 0.0;
    406         apparent->d = 0.0;
    407         apparent->rErr = 0.0;
    408         apparent->dErr = 0.0;
    409     }
    410 
     407    /*    if (apparent == NULL) {
     408            apparent = psSphereAlloc();
     409        } else {
     410            apparent->r = 0.0;
     411            apparent->d = 0.0;
     412            apparent->rErr = 0.0;
     413            apparent->dErr = 0.0;
     414        }
     415    */
     416    if (apparent != NULL) {
     417        psFree(apparent);
     418    }
    411419    // bend the actual vector away from the sun vector by deflection angle.
    412420    // XXX: Not sure how to do this.  Dave thinks the formula should be:
    … …  
    414422    theta = 0.0;
    415423    double phi = 0.0;
    416     //    deflection = SEC_TO_RAD(deflection);
     424    //    deflection = SEC_TO_RAD(deflection) * 1e6;
     425    //    deflection *= M_PI * 1e-2;
    417426    theta = atan(r0/PS_AU) * tan(deflection);
    418     printf(" Theta = %.13g\n", theta);
    419     printf(" deflection = %.13g\n", deflection);
    420427    //    phi = sqrt( deflection*deflection - theta*theta );
    421     phi = deflection * cos(asin(theta/deflection)) * 3e-2;
     428    //    phi = deflection * cos(asin(theta/deflection));
     429
     430    //    phi = sqrt(theta*theta - deflection*deflection);
     431    //    phi = deflection * cos(asin(theta/deflection)) * 3e-2;
    422432    //    phi = cos(asin(theta/deflection));
    423433    //    phi = asin(theta/deflection);
    424     apparent->r = theta;
    425     apparent->d = phi;
     434
     435    //    apparent->r = theta;
     436    //    apparent->d = phi;
     437    /*
     438        actualVector->x += actualVector->x*deflection;
     439        actualVector->y += actualVector->y*deflection;
     440        actualVector->z += actualVector->z*deflection;
     441        apparent = psCubeToSphere(actualVector);
     442    */
     443    theta = tan(sun->r - actual->r) * deflection;
     444    phi = tan(sun->d - actual->d) * deflection;
     445
     446    printf(" Theta = %.13g\n", theta);
     447    printf(" phi = %.13g\n", phi);
     448
     449    temp->r = theta;
     450    temp->d = phi;
     451    apparent = psSphereSetOffset(actual, temp, PS_SPHERICAL, PS_RADIAN);
    426452
    427453    psFree(actualVector);
    428454    psFree(sunVector);
     455    psFree(temp);
     456
    429457    return apparent;
    430458}
    … …  
    807835}
    808836
    809 psSphereRot* psSphereRot_TEOtoCEO(const psTime *time)
     837psSphereRot* psSphereRot_TEOtoCEO(const psTime *time,
     838                                  psEarthPole *tidalCorr)
    810839{
    811840    PS_ASSERT_PTR_NON_NULL(time,NULL);
    … …  
    815844        in = psTimeConvert(in, PS_TIME_UT1);
    816845    }
    817     //    double T = (double)(in->sec) + (double)(in->nsec / 1e9);
     846    if (tidalCorr != NULL && tidalCorr->s != 0.0) {
     847        int nsec = in->nsec + (int)(tidalCorr->s * 1e9);
     848        if (nsec < 0.0) {
     849            in->sec += -1;
     850            in->nsec = (int)(1e9) - nsec;
     851        } else {
     852            in->nsec = nsec;
     853        }
     854    }
    818855    double T = psTimeToJD(in);
    819856    T += -2451545.0;
    820857    double theta = 2.0 * M_PI * (0.7790572732640 + 1.00273781191135448 * T);
    821858    psSphereRot *out = psSphereRotAlloc(theta, 0.0, 0.0);
    822     //    psSphereRot *out = psSphereRotInvert(theta, 0.0, 0.0);
    823859
    824860    psFree(in);
    … …  
    919955    out->x = SEC_TO_RAD(xOut);
    920956    out->y = SEC_TO_RAD(yOut);
     957    //    psEarthPole *polarTideCorr = psEOC_PolarTideCorr(time);
     958    //    out->x += polarTideCorr->x;
     959    //    out->y += polarTideCorr->y;
     960    //    psFree(polarTideCorr);
     961
    921962    //    out->s = SEC_TO_RAD(sOut);
    922963
    … …  
    10341075    CORX = SEC_TO_RAD(CORX);
    10351076    CORY = SEC_TO_RAD(CORY);
    1036     CORZ = SEC_TO_RAD(CORZ);
     1077    //    CORZ = SEC_TO_RAD(CORZ);
    10371078
    10381079    out->x = CORX;

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

    r5533 r6030  
    99*  @author Robert Daniel DeSonia, MHPCC
    1010*
    11 *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    12 *  @date $Date: 2005-11-17 03:59:05 $
     11*  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     12*  @date $Date: 2006-01-18 00:41:29 $
    1313*
    1414*  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    6262psEarthPole *psEarthPoleAlloc(void);
    6363
    64 /** Calculates the actual position of a star, given its apparent position and the
     64/** Calculates the apparent position of a star, given its actual position and the
    6565 *  velocity vector of the observer.
    6666 *
    6767 *  The actual and apparent positions are represented as psSphere entries, as is the
    6868 *  direction of motion.  The speed in that direction is given in units of the speed
    69  *  of light.  If the value of actual is NULL, a new psSphere is allocated, otherwise
    70  *  the point to actual is used for the result.
     69 *  of light.  If the value of apparent is NULL, a new psSphere is allocated, otherwise
     70 *  the point to apparent is used for the result.
    7171 *
    7272 *  @return psSphere*:      the actual position of a star.
    … …  
    7979);
    8080
    81 /** Calculates the actual position of a star, given its apparent position and the
     81/** Calculates the apparent position of a star, given its actual position and the
    8282 *  position of the sun.
    8383 *
    8484 *  The actual and apparent positions are represented as psSphere entries, as is
    85  *  position of the sun.  If the value of actual is NULL, a new psSphere is allocated,
    86  *  otherwise the point to actual is used for the result.
     85 *  position of the sun.  If the value of apparent is NULL, a new psSphere is allocated,
     86 *  otherwise the point to apparent is used for the result.
    8787 *
    8888 *  @return psSphere*:      the apparent position of a star.
    … …  
    133133 */
    134134psSphereRot *psSphereRot_TEOtoCEO(
    135     const psTime *time                 ///< specified time
     135    const psTime *time,                ///< specified time
     136    psEarthPole *tidalCorr             ///< UT1 polar tide correction or NULL
    136137);
    137138

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

    r5969 r6030  
    55*  @author d-Rob, MHPCC
    66*
    7 *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
    8 *  @date $Date: 2006-01-12 20:40:13 $
     7*  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
     8*  @date $Date: 2006-01-18 00:41:29 $
    99*
    1010*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    8585{
    8686    double r1, r2, d1, d2;
    87     r1 = in1->r;
    88     r2 = in2->r;
    89     d1 = in1->d;
    90     d2 = in2->d;
     87    d1 = in1->r;
     88    d2 = in2->r;
     89    r1 = in1->d;
     90    r2 = in2->d;
    9191    double out = 0.0;
    9292    double c1, c2, cd, s1, s2, sum, ac;
    … …  
    101101
    102102    ac = acos(sum);
    103     printf("\n  c1=%lf, c2=%lf, cd=%lf, s1=%lf, s2=%lf, sum=%.19g, ac=%g\n", c1,c2,cd,s1,s2,sum,ac);
     103    out = ac;
     104    //    printf("\n  c1=%lf, c2=%lf, cd=%lf, s1=%lf, s2=%lf, sum=%.19g, ac=%g\n", c1,c2,cd,s1,s2,sum,ac);
    104105    //    out = acos(cos(r1)*cos(r2)*cos(d1-d2) + sin(r1)*sin(r2));
    105106    //    if (out != 0.0) printf("\n   in greatCircle, out = %lf\n", out);
    … …  
    109110psS32 testAberration(void)
    110111{
    111 
    112112    psSphere *apparent = NULL;
    113     //    psSphere *actual = psSphereAlloc();
    114     //    psSphere *direction = psSphereAlloc();
    115113    psSphere *empty = NULL;
    116 
    117114    psCube *actualCube = psCubeAlloc();
    118     //    actual->r = DEG_TO_RAD(122.9153182445501);
    119     //    actual->d = DEG_TO_RAD(48.562968978679194);
    120     //    actualCube->x = -0.3596195125758298;
    121     //    actualCube->y = 0.5555613903455866;
    122     //    actualCube->z = 0.7496834983724809;
    123 
    124115    //values from After gravity deflection//
    125116    actualCube->x = -0.35961949760293604;
    … …  
    127118    actualCube->z = 0.7496835020836093;
    128119    psSphere *actual = psCubeToSphere(actualCube);
    129 
    130120    psCube *cubeDir = psCubeAlloc();
    131121    cubeDir->x = 5148.713262821658;
    … …  
    135125    cubeDir->y += 248.46429758174693;
    136126    cubeDir->z += 0.09694774143797581;
    137     double length = sqrt(cubeDir->x*cubeDir->x+cubeDir->y*cubeDir->y+cubeDir->z*cubeDir->z);
    138     cubeDir->x = cubeDir->x/length;
    139     cubeDir->y = cubeDir->y/length;
    140     cubeDir->z = cubeDir->z/length;
    141127    psSphere *direction = psCubeToSphere(cubeDir);
    142 
    143128    double speed = sqrt(cubeDir->x*cubeDir->x + cubeDir->y*cubeDir->y + cubeDir->z*cubeDir->z);
    144     // Speed of light in vacuum (src:NIST)
    145     double c = 299792458.0; /* m/s */
    146     //    speed = speed / c;
    147     speed = length / c;
     129    double c = 299792458.0; // Speed of light in vacuum (src:NIST)   /* m/s */
     130    speed /= c;
    148131
    149132    empty = psAberration(empty, apparent, direction, speed);
    … …  
    161144
    162145    apparent = psAberration(apparent, actual, direction, speed);
    163 
    164     //    printf("\nSphere Difference  =  r,d  =  %.13g, %.13g\n",
    165     //           (actual->r - apparent->r), (actual->d - apparent->d));
    166146    psCube *outCube = psSphereToCube(apparent);
    167     printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
     147    printf("\n -- resultCube = x,y,z =     %.13g,      %.13g,    %.13g  -- \n",
    168148           outCube->x, outCube->y, outCube->z);
    169149    //expected cube values
    … …  
    173153    z = 0.7497078321908413;
    174154
    175     printf(" -- expectedCube = x,y,z = %.13g, %.13g, %.13g  -- \n", x, y, z);
     155    printf(" -- expectedCube = x,y,z =   %.13g,     %.13g,    %.13g  -- \n", x, y, z);
    176156    printf("Cube Difference  =  x,y,z  = %.13g, %.13g, %.13g \n\n",
    177157           (x - outCube->x), (y - outCube->y), (z - outCube->z) );
    178158    psFree(actual);
    179     //        psFree(actualCube);
    180159    actualCube->x = x;
    181160    actualCube->y = y;
    … …  
    190169        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    191170                "psAberration returned incorrect values.\n");
    192         printf("\nMagnitude of expected change = x,y,z = %.13g, %.13g, %.13g \n",
    193                (actualCube->x - x), (actualCube->y - y), (actualCube->z - z) );
    194         printf("Magnitude of actual change = x,y,z = %.13g, %.13g, %.13g \n",
    195                (actualCube->x - outCube->x), (actualCube->y - outCube->y),
    196                (actualCube->z - outCube->z) );
     171        //        printf("\nMagnitude of expected change = x,y,z = %.13g, %.13g, %.13g \n",
     172        //               (actualCube->x - x), (actualCube->y - y), (actualCube->z - z) );
     173        //        printf("Magnitude of actual change = x,y,z = %.13g, %.13g, %.13g \n",
     174        //               (actualCube->x - outCube->x), (actualCube->y - outCube->y),
     175        //               (actualCube->z - outCube->z) );
    197176        psFree(actual);
    198         //        psFree(actualCube);
    199177        actualCube->x = x;
    200178        actualCube->y = y;
    … …  
    209187    psFree(outCube);
    210188    psFree(actualCube);
    211 
    212189    psFree(cubeDir);
    213190    psFree(apparent);
    … …  
    220197psS32 testGravityDeflect(void)
    221198{
    222 
    223     //    psSphere *actual = psSphereAlloc();
    224199    psSphere *apparent = NULL;
    225     //    psSphere *sun = psSphereAlloc();
    226200    psSphere *empty = NULL;
    227 
    228201    psCube *actualCube = psCubeAlloc();
    229202    actualCube->x = -0.3596195125758298;
    … …  
    235208    sunCube->y = 2.5880956908748722e10;
    236209    sunCube->z = 1.1220046291457653e10;
     210    double sunLength = sqrt(sunCube->x*sunCube->x + sunCube->y*sunCube->y + sunCube->z*sunCube->z);
     211    sunCube->x /= sunLength;
     212    sunCube->y /= sunLength;
     213    sunCube->z /= sunLength;
     214    printf("sunCube = x,y,z = %.13g, %.13g, %.13g\n", sunCube->x, sunCube->y, sunCube->z);
    237215    psSphere *sun = psCubeToSphere(sunCube);
     216    printf("sunSphere  = r, d = %.13g, %.13g\n", sun->r, sun->d);
     217    psCube *outCube = psCubeAlloc();
    238218
    239219    empty = psGravityDeflection(apparent, empty, sun);
    … …  
    251231
    252232    apparent = psGravityDeflection(NULL, actual, sun);
    253     psSphere *result2 = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
    254     apparent->r *= -1.0;
    255     apparent->d *= -1.0;
    256     psSphere *result = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
     233    //    apparent->r *= -1.0;
     234    //    apparent->d *= -1.0;
     235    psSphere *result2;
     236    //    psSphere *result2 = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
     237    //    psSphere *result = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
    257238    //    psSphere *result = psSphereGetOffset(apparent, actual, PS_SPHERICAL, PS_RADIAN);
    258239    printf(" -- actualCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
    259240           actualCube->x, actualCube->y, actualCube->z);
    260     psCube *outCube = psSphereToCube(result);
    261     printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
    262            outCube->x, outCube->y, outCube->z);
    263     psCube *outCube2 = psSphereToCube(result2);
     241    //    psCube *outCube = psSphereToCube(result);
     242    //    printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
     243    //           outCube->x, outCube->y, outCube->z);
     244    psCube *outCube2 = psSphereToCube(apparent);
    264245    printf(" -- resultCube2= x,y,z = %.13g, %.13g, %.13g  -- \n",
    265246           outCube2->x, outCube2->y, outCube2->z);
    … …  
    275256    //        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    276257    //                "psGravityDeflection returned incorrect values.\n");
    277     printf("expect-actual=  x,y,z  = %.11g, %.11g, %.11g \n",
     258    printf("expect-actual=  x,y,z  = %.13g, %.13g, %.13g \n",
    278259           (x - actualCube->x), (y - actualCube->y), (z - actualCube->z) );
    279     printf("result-actual=  x,y,z  = %.11g, %.11g, %.11g \n",
    280            (outCube->x - actualCube->x), (outCube->y - actualCube->y), (outCube->z - actualCube->z) );
    281     printf("expect-result=  x,y,z  = %.11g, %.11g, %.11g \n",
    282            (x - outCube->x), (y - outCube->y), (z - outCube->z) );
     260    printf("expect-result=  x,y,z  = %.13g, %.13g, %.13g \n",
     261           (x - outCube2->x), (y - outCube2->y), (z - outCube2->z) );
     262    printf("result-actual=  x,y,z  = %.13g, %.13g, %.13g \n",
     263           (outCube2->x - actualCube->x), (outCube2->y - actualCube->y),
     264           (outCube2->z - actualCube->z) );
    283265    //        return 1;
    284266    //    }
    285     psFree(result2);
     267    //    psFree(result2);
    286268    outCube2->x = x;
    287269    outCube2->y = y;
    288270    outCube2->z = z;
    289271    result2 = psCubeToSphere(outCube2);
    290     psFree(result);
    291     result = psSphereGetOffset(actual, result2, PS_SPHERICAL, PS_RADIAN);
    292     printf("The apparent output sphere = r,d = %.13g, %.13g\n", apparent->r, apparent->d);
     272    //    psFree(result);
     273    psSphere *result = psSphereGetOffset(actual, result2, PS_SPHERICAL, PS_RADIAN);
     274    psFree(result2);
     275    result2 = psSphereGetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
     276    printf("The apparent output sphere = r,d = %.13g, %.13g\n", result2->r, result2->d);
    293277    printf("The expected output sphere = r,d = %.13g, %.13g\n\n", result->r, result->d);
    294278    psFree(result2);
    … …  
    349333            return 4;
    350334        }
    351         printf("\n  Precession Model output = x,y,s = %.8g, %.8g, %.8g\n",
     335        printf("\n  PrecessionModel output = x,y,s = %.13g,  %.13g,  %.13g\n",
    352336               pmodel->x, pmodel->y, pmodel->s);
    353         printf("  Expected output = x,y,s = %.13g, %.13g, %.13g\n", x, y, s);
    354         printf("  A difference of:   %.13g, %.13g, %.13g\n\n",
     337        printf("  Expected output        = x,y,s = %.13g, %.13g, %.13g\n", x, y, s);
     338        printf("  A difference of:                 %.13g, %.13g, %.13g\n",
    355339               (pmodel->x - x), (pmodel->y - y), (pmodel->s - s) );
    356340    }
    … …  
    402386        return 7;
    403387    } else {
    404         printf("\nPrecession Correction output (IERSA) = x,y,s = %.13g, %.13g, %.13g\n\n",
     388        printf("\nPrecessionCorr output (IERSA) = x,y,s = %.13g,  %.13g, %.13g\n",
    405389               pcorr->x, pcorr->y, pcorr->s);
    406390    }
    407391    psFree(pcorr);
    408 
    409     //    printf("\n >>MJD value = %lf \n", psTimeToMJD(time2));
    410392
    411393    //Check values from IERS table B
    … …  
    428410        printf("PrecessionCorr output (IERSB) = x,y,s = %.13g, %.13g, %.13g\n",
    429411               pcorr->x, pcorr->y, pcorr->s);
    430         printf("Expected output               = x,y,s = %.13g, %.13g, %.13g\n\n", xx, yy, ss);
    431         printf("  A difference of:   %.13g, %.13g, %.13g\n\n",
     412        printf("Expected output               = x,y,s = %.13g, %.13g, %.13g\n", xx, yy, ss);
     413        printf("          A difference of:              %.13g, %.13g, %.13g\n\n",
    432414               (pcorr->x - xx), (pcorr->y - yy), (pcorr->s - ss) );
    433415        //            return 10;
    434416        //        }
    435417    }
    436 
    437 
     418    //precess is the *actual* output from PrecessionModel + PrecessionCorr
     419    psEarthPole *precess = psEOC_PrecessionModel(time2);
     420    precess->x += pcorr->x;
     421    precess->y += pcorr->y;
    438422    double xCorr, yCorr;
    439423    xCorr = 3.05224300720406e-10;
    … …  
    445429    pcorr->x += xCorr;
    446430    pcorr->y += yCorr;
    447     psEarthPole *precess = psEOC_PrecessionModel(time2);
    448     precess->x += xCorr;
    449     precess->y += yCorr;
    450431    psSphereRot *precessNutInv = psSphereRot_CEOtoGCRS(precess);
    451432    psSphereRot *precessNut = psSphereRotConjugate(NULL, precessNutInv);
    452     //    pcorr->x += 3.05224300720406e-7;
    453     //    pcorr->y += -1.39441339235822e-7;
    454433    double q0, q1, q2;
    455434    q0 = -1.1984522406756289e-5;
    … …  
    461440        printf("\n Error at CEOtoGCRS, output psSphereRot doesn't match expected.\n");
    462441    }
    463     printf("  Output sphere rotation   = %.13g,%.13g,%.13g,%.13g\n",
    464            pni->q0, pni->q1, pni->q2, pni->q3);
    465     printf("  Expected sphere rotation = %.13g,%.13g,%.13g\n", q0,q1,q2);
    466     printf("  MY sphere rotation       = %.13g,%.13g,%.13g\n",
     442    //    printf("  Output from CEOtoGCRS only  = %.13g,%.13g,%.13g,%.13g\n",
     443    //           pni->q0, pni->q1, pni->q2, pni->q3);
     444    printf("  Expected sphere rotation    = %.13g, %.13g, %.13g\n", q0,q1,q2);
     445    printf("  Result sphere rotation      = %.13g, %.13g, %.13g\n",
    467446           precessNutInv->q0, precessNutInv->q1, precessNutInv->q2);
     447    printf("     Difference         =        %.13g, %.13g, %.13g\n\n",
     448           precessNutInv->q0-q0, precessNutInv->q1-q1, precessNutInv->q2-q2);
    468449    psCube *objC = psCubeAlloc();
    469450    //    objC->x = -3.5963388069046304;
    … …  
    491472    y = 0.5555012823608123;
    492473    z = 0.7496183628158023;
    493     printf("\n<<Expected out      = x,y,z = %.13g, %.13g, %.13g\n", x, y, z);
    494     psFree(objC);
    495     objC = psSphereToCube(expect);
    496     printf("<<Expected out real = x,y,z = %.13g, %.13g, %.13g\n", objC->x, objC->y, objC->z);
    497     printf("     Difference     =      %.13g, %.13g, %.13g\n", objC->x-x, objC->y-y, objC->z-z);
    498     x = objC->x;
    499     y = objC->y;
    500     z = objC->z;
     474    printf("\n<<Expected out       = x,y,z = %.13g, %.13g, %.13g\n", x, y, z);
     475    //    psFree(objC);
     476    //    objC = psSphereToCube(expect);
     477    //printf("<<Expected out (CEO) = x,y,z = %.13g, %.13g, %.13g\n", objC->x, objC->y, objC->z);
     478    //printf("     Difference     =           %.13g, %.13g, %.13g\n", objC->x-x, objC->y-y, objC->z-z);
     479    //    x = objC->x;
     480    //    y = objC->y;
     481    //    z = objC->z;
    501482    psSphere *result = psSphereRotApply(NULL, precessNut, sphere);
    502483    psFree(objC);
    503484    objC = psSphereToCube(result);
    504     printf("<<Resulting out     = x,y,z = %.13g, %.13g, %.13g\n", objC->x, objC->y, objC->z);
    505     printf("     Difference     =      %.13g, %.13g, %.13g\n\n", objC->x-x, objC->y-y, objC->z-z);
     485    printf("<<Resulting out      = x,y,z = %.13g, %.13g, %.13g\n", objC->x, objC->y, objC->z);
     486    printf("     Difference         =      %.13g, %.13g, %.13g\n\n", objC->x-x, objC->y-y, objC->z-z);
    506487
    507488    double xx = greatCircle(result, expect);
    … …  
    562543        return 4;
    563544    }
    564     if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
    565             || fabs(polarMotion->s - s) > DBL_EPSILON) {
    566         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    567                 "psEOC_GetPolarMotion returned incorrect values.\n");
    568         printf("\n  <>PolarMotion output (IERSA)   = x,y,s = %.13g, %.13g, %.13g\n",
    569                polarMotion->x, polarMotion->y, polarMotion->s);
    570         printf("\n  <>PolarMotion expected (IERSA) = x,y,s = %.13g, %.13g, %.13g\n",
    571                x, y, s);
    572         //        return 5;
    573     }
     545    //    if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
     546    //            || fabs(polarMotion->s - s) > DBL_EPSILON) {
     547    //        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     548    //                "psEOC_GetPolarMotion returned incorrect values.\n");
     549    printf("  <>PolarMotion output (IERSA)   = x,y,s = %.13g, %.13g, %.13g\n",
     550           polarMotion->x, polarMotion->y, polarMotion->s);
     551    //        printf("  <>PolarMotion expected (IERSA) = x,y,s = %.13g, %.13g, %.13g\n",
     552    //               x, y, s);
     553    //        return 5;
     554    //    }
    574555    psFree(polarMotion);
    575 
    576556
    577557    //Return valid values for correct input time.  Test IERS Bulletin B.
    578558    polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B);
    579     x = -6.45381397904e-07;
    580     y = 2.112819726698e-06;
    581     s = 0.0;
    582 
     559    //    x = -6.45381397904e-07;
     560    //    y = 2.112819726698e-06;
     561    //    s = 0.0;
    583562    if (polarMotion == NULL) {
    584563        psError(PS_ERR_BAD_PARAMETER_NULL, false,
    … …  
    586565        return 2;
    587566    }
    588     if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
    589             || fabs(polarMotion->s - s) > DBL_EPSILON) {
    590         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    591                 "psEOC_GetPolarMotion returned incorrect values.\n");
    592         printf("\n  <>PolarMotion output (IERSB)   = x,y,s = %.13g, %.13g, %.13g\n",
    593                polarMotion->x, polarMotion->y, polarMotion->s);
    594         printf("\n  <>PolarMotion expected (IERSB) = x,y,s = %.13g, %.13g, %.13g\n",
    595                x, y, s);
    596         //        return 3;
    597     }
    598 
    599 
     567    //    if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
     568    //            || fabs(polarMotion->s - s) > DBL_EPSILON) {
     569    //        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     570    //                "psEOC_GetPolarMotion returned incorrect values.\n");
     571    printf("  <>PolarMotion output (IERSB)   = x,y,s = %.13g,  %.13g, %.13g\n",
     572           polarMotion->x, polarMotion->y, polarMotion->s);
     573    //        printf("  <>PolarMotion expected (IERSB) = x,y,s = %.13g, %.13g, %.13g\n",
     574    //               x, y, s);
     575    //        return 3;
     576    //    }
    600577    psEarthPole *nutationCorr = psEOC_NutationCorr(in);
    601578    if (nutationCorr == NULL) {
    … …  
    603580        return 6;
    604581    }
    605 
    606582    polarMotion->x += nutationCorr->x;
    607583    polarMotion->y += nutationCorr->y;
    608584    polarMotion->s += nutationCorr->s;
     585    psEarthPole *polarTide = psEOC_PolarTideCorr(in);
     586    polarMotion->x += polarTide->x;
     587    polarMotion->y += polarTide->y;
     588    psFree(polarTide);
    609589    x = -6.43607313124045e-7;
    610590    y = 2.11351436973568e-6;
    611591    s = -7.39617581324646e-12;
    612 
    613     if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
    614             || fabs(polarMotion->s - s) > DBL_EPSILON) {
    615         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    616                 "   psEOC_GetPolarMotion returned incorrect values.\n");
    617         printf("\n  PolarMotion + NutationCorr out = x,y,s = %.13g, %.13g, %.13g\n",
    618                polarMotion->x, polarMotion->y, polarMotion->s);
    619         printf("  Expected output = x,y,s = %.13g, %.13g, %.13g\n", x, y, s);
    620         printf("\n  Output Nutation = x,y,s = %.13g, %.13g, %.13g\n",
    621                nutationCorr->x, nutationCorr->y, nutationCorr->s);
    622         //        printf("  A difference of:   %.13g, %.13g, %.13g\n\n",
    623         //               (polarMotion->x - x), (polarMotion->y - y), (polarMotion->s - s) );
    624         //        return 10;
    625     }
    626 
    627     psEarthPole *polarTide = psEOC_PolarTideCorr(in);
    628     polarMotion->x += polarTide->x;
    629     polarMotion->y += polarTide->y;
    630     polarMotion->s += polarTide->s;
    631     psFree(polarTide);
    632     printf("\n  PolarMotion + PolarTideCorr out = x,y,s = %.13g, %.13g, %.13g\n",
     592    //    if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
     593    //            || fabs(polarMotion->s - s) > DBL_EPSILON) {
     594    //        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     595    //                "   psEOC_GetPolarMotion returned incorrect values.\n");
     596    printf("\n  PolarMotion + NutationCorr out = x,y,s = %.13g, %.13g, %.13g\n",
    633597           polarMotion->x, polarMotion->y, polarMotion->s);
     598    printf("  Expected output                = x,y,s = %.13g,  %.13g, %.13g\n", x, y, s);
     599    printf("                     Difference  = x,y,s = %.13g, %.13g,  %.13g\n",
     600           polarMotion->x - x, polarMotion->y - y, polarMotion->s - s);
     601    //    }
    634602
    635603    if (!p_psEOCFinalize() ) {
    … …  
    679647{
    680648    psTime *in = psTimeAlloc(PS_TIME_UTC);
    681     //    in->sec = 1131579114;
    682     //    in->nsec = 498489000;
    683649    in->sec = timesec;
    684650    in->nsec = 0;
    … …  
    728694    psSphereRot *rot = NULL;
    729695    psTime *empty = NULL;
    730     //    psTime *now = psTimeAlloc(PS_TIME_UT1);
    731     //    now->sec = 1128530000;
    732     //    now->nsec = 931154510;
    733696    psTime *time = psTimeAlloc(PS_TIME_UT1);
    734697    time->sec = timesec-1;
    … …  
    736699    time->leapsecond = false;
    737700
    738     //    psSphereRot *teoceo = psSphereRot_TEOtoCEO(now);
    739 
    740701    //return NULL for NULL input time
    741     rot = psSphereRot_TEOtoCEO(empty);
     702    rot = psSphereRot_TEOtoCEO(empty, NULL);
    742703    if (rot != NULL) {
    743704        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    … …  
    746707    }
    747708
    748     psSphereRot *teoceo = psSphereRot_TEOtoCEO(time);
     709    psEarthPole *polarTideCorr = psEOC_PolarTideCorr(time);
     710    psSphereRot *teoceo = psSphereRot_TEOtoCEO(time, polarTideCorr);
    749711    //Make sure values match for other psTime type
    750712    empty = psTimeAlloc(PS_TIME_UTC);
    … …  
    752714    empty->nsec = 0;
    753715    empty->leapsecond = false;
    754     rot = psSphereRot_TEOtoCEO(empty);
     716
     717    rot = psSphereRot_TEOtoCEO(empty, polarTideCorr);
    755718    if (fabs(rot->q0-teoceo->q0) > DBL_EPSILON || fabs(rot->q1-teoceo->q1) > DBL_EPSILON ||
    756719            fabs(rot->q2-teoceo->q2) > DBL_EPSILON || fabs(rot->q3-teoceo->q3) > DBL_EPSILON) {
    757720        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    758721                "psSphereRot_TEOtoCEO failed to return matching values for different time types.\n");
     722        printf("\n  Output Rotation1 = q0,q1,q2,q3 = %.13g, %.13g, %.13g, %.13g\n",
     723               teoceo->q0, teoceo->q1, teoceo->q2, teoceo->q3 );
     724        printf("\n  Output Rotation2 = q0,q1,q2,q3 = %.13g, %.13g, %.13g, %.13g\n",
     725               rot->q0, rot->q1, rot->q2, rot->q3 );
    759726        return 2;
    760727    }
    … …  
    783750        return 3;
    784751    }
    785     //    psFree(now);
     752
     753    psFree(polarTideCorr);
    786754    psFree(rot);
    787755    psFree(empty);
    … …  
    800768    psEarthPole *empty = NULL;
    801769    psSphereRot *rot = NULL;
    802 
    803     //    in->x = M_PI / 4.0;
    804     //    in->y = M_PI / 6.0;
    805     //    in->s = M_PI / 8.0;
    806 
    807770    in->x = 2.857175590089105e-4;
    808771    in->y = 2.3968739377734732e-5;
    809772    in->s = -1.3970066457904322e-8;
    810     //    in->x += 0.06295703125;
    811     //    in->y += -0.0287618408203125;
    812     //    in->s += 0.0;
    813     //    in->x += 3.05224300720406e-7;
    814     //    in->y += -1.39441339235822e-7;
    815     in->s += 0.0;
    816     //    in->x += -2.491942320903e-10;
    817     //    in->y += -1.648366515772e-11;
    818773
    819774    double q0,q1,q2,q3;
    … …  
    837792    }
    838793
    839     if (fabs(rot->q0-q0) > DBL_EPSILON || fabs(rot->q1-q1) > DBL_EPSILON ||
    840             fabs(rot->q2-q2) > DBL_EPSILON || fabs(rot->q3-q3) > DBL_EPSILON) {
     794    printf("\n  Output sphere rotation   = %.13g, %.13g, %.13g, %.13g\n",
     795           rot->q0, rot->q1, rot->q2, rot->q3);
     796    printf("  Expected sphere rotation = %.13g, %.13g, %.13g, %.13g\n", q0,q1,q2,q3);
     797    printf("  difference:                 %.13g, %.13g, %.13g \n",
     798           (rot->q0-q0), (rot->q1-q1), (rot->q2-q2) );
     799    if (fabs(rot->q0-q0) > FLT_EPSILON || fabs(rot->q1-q1) > FLT_EPSILON ||
     800            fabs(rot->q2-q2) > FLT_EPSILON || fabs(rot->q3+q3) > FLT_EPSILON) {
    841801        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    842802                "psSphereRot_CEOtoGCRS failed to return expected values.\n");
    843         printf("\n  Output sphere rotation   = %.13g,%.13g,%.13g,%.13g\n",
    844                rot->q0, rot->q1, rot->q2, rot->q3);
    845         printf("  Expected sphere rotation = %.13g,%.13g,%.13g,%.13g\n", q0,q1,q2,q3);
    846         printf("  a difference:   %.13g, %.13g, %.13g, %.13g \n", (rot->q0-q0), (rot->q1-q1),
    847                (rot->q2-q2), (rot->q3-q3) );
    848         printf("  abs difference:   %.13g, %.13g, %.13g, %.13g \n\n", (rot->q0+q0), (rot->q1+q1),
    849                (rot->q2+q2), (rot->q3+q3) );
    850         //return 3;
    851     }
    852 
    853 
     803        return 3;
     804    }
    854805    psCube *tempCube = psCubeAlloc();
    855806    tempCube->x = -0.35963388069046304;
    … …  
    858809    obj = psCubeToSphere(tempCube);
    859810    psFree(tempCube);
    860     //    obj = psSphereAlloc();
    861     //    obj->r = objR;
    862     //    obj->d = objD;
    863811    psSphereRot *precessionNutation = psSphereRotConjugate(NULL, rot);
    864812    psSphere *result = psSphereRotApply(NULL, precessionNutation, obj);
    … …  
    868816    y = 0.5555012823608123;
    869817    z = 0.7496183628158023;
    870     printf("\nOutput cube = x,y,z = %.13g, %.13g, %.13g\n", cube->x, cube->y, cube->z);
    871     printf("Expected cube = x,y,z = %.13g, %.13g, %.13g\n", x, y, z);
    872     printf("A difference of:   %.13g, %.13g, %.13g\n\n",
     818    printf("\n  Output cube = x,y,z = %.13g,  %.13g,    %.13g\n", cube->x, cube->y, cube->z);
     819    printf("Expected cube = x,y,z = %.13g,  %.13g,  %.13g\n", x, y, z);
     820    printf("  A difference of:    %.13g, %.13g, %.13g\n\n",
    873821           (x-cube->x), (y-cube->y), (z-cube->z));
    874 
    875     //great sphere difference between result and expected
    876     //Re = 6.38x10^6m
    877822    psCube *expect = psCubeAlloc();
    878823    expect->x = x;
    … …  
    880825    expect->z = z;
    881826    psSphere *expected = psCubeToSphere(expect);
     827    double d = greatCircle(result, expected);
     828    printf("   The great circle angular distance between expected & result = %.13g\n", d);
     829
    882830    psFree(expect);
    883     double d = acos(cos(result->r)*cos(expected->r)*cos(result->d-expected->d) + sin(result->r)*sin(expected->r));
    884     printf("Great circle difference of:  %.13g \n", d);
    885 
    886831    psFree(expected);
    887832    psFree(obj);
    888 
    889 
    890 
    891     //    psSphere *test = psSphereAlloc();
    892     //    test->r = 0.0;
    893     //    test->d = 0.0;
    894     //    test = psSphereRotApply(test, rot, test);
    895     //    printf("\n  Sphere -test- has values r,d = %.8g, %.8g \n", test->r, test->d);
    896 
    897833    psFree(precessionNutation);
    898834    psFree(result);
    899835    psFree(cube);
    900 
    901     //    psFree(test);
    902836    psFree(rot);
    903837    psFree(in);
     838
    904839    return 0;
    905840}
    … …  
    910845    psEarthPole *empty = NULL;
    911846    psSphereRot *rot = NULL;
    912 
    913     //    in->x = M_PI / 4.0;
    914     //    in->y = M_PI / 6.0;
    915     //    in->s = M_PI / 8.0;
    916847    in->x = -0.13275353774074533;
    917848    in->y = 0.4359436319739848;
    918849    in->s = -4.2376965863576153e-10;
    919 
    920850    in->x = SEC_TO_RAD(in->x);
    921851    in->y = SEC_TO_RAD(in->y);
    … …  
    928858        return 1;
    929859    }
    930 
    931860    rot = psSphereRot_ITRStoTEO(in);
    932861    if (rot == NULL) {
    … …  
    958887    printf("  Expected sphere rotation = %.13g, %.13g, %.13g, %.13g\n", q0,q1,q2,q3);
    959888
    960     //    psSphere *test = psSphereAlloc();
    961     //    test->r = 0.0;
    962     //    test->d = 0.0;
    963889    psCube *temp = psCubeAlloc();
    964890    //    temp->x = -0.3596195125758298;
    … …  
    969895    temp->z = 0.7496183628158023;
    970896
    971 
    972897    obj = psCubeToSphere(temp);
    973898    psSphere *test = NULL;
    … …  
    975900    psSphereRot *newRot = psSphereRotConjugate(NULL, rot);
    976901    test = psSphereRotApply(NULL, newRot, obj);
    977     printf("\n  Sphere -test- has values r,d = %.8g, %.8g \n", test->r, test->d);
    978902    temp = psSphereToCube(test);
    979903    printf("\n  Cube -test- has x,y,z =     %.13g, %.13g, %.13g \n", temp->x, temp->y, temp->z);
    … …  
    988912    psSphere *sphere = psCubeToSphere(temp);
    989913    double d = greatCircle(sphere, test);
    990 
    991914    printf("Great circle difference of:  %.13g \n", d);
     915
    992916    psFree(sphere);
    993 
    994917    psFree(newRot);
    995918    psFree(obj);
    … …  
    998921    psFree(rot);
    999922    psFree(in);
    1000     return 0;
    1001 }
    1002 
     923
     924    return 0;
     925}
     926
Note: See TracChangeset for help on using the changeset viewer.