All angular quantities in this document are in radians, following the psLib convention.


1. Start with RA, Dec:

    r = 2.14527700450406
    d = 0.84758370322182

These get stuffed into a psSphere (obj): psSphereAlloc, followed by
explicitly setting the values in the struct.


2. Convert obj to a psCube (psSphereToCube), and check that:

    x=-0.3596195125758298
    y=0.5555613903455866
    z=0.7496834983724809


3. Stuff the Sun position into a psCube (sunCube):

    x=1.467797790127511E11
    y=2.5880956908748722E10
    z=1.1220046291457653E10

Convert to a psSphere (psCubeToSphere *sun) for later use.


4. Gravitational deflection:

    psGravityDeflection(obj, obj, sun);

Convert obj to a psCube, and check that:

    x=-0.35961949760293604
    y=0.5555613950298085
    z=0.7496835020836093


5. Earth's direction of motion is:

Barycentric velocity of the Earth in m/s:
    x=5148.713262821658
    y=-26945.04752348012
    z=-11682.787302030947
Diurnal velocity of the observer in m/s:
    x=-357.6031690489248
    y=248.46429758174693
    z=0.09694774143797581

These two vectors need to be summed (standard vector addition).
Converting to a psSphere gives the direction of motion (psSphere
*direction), while the magnitude gives the speed, after dividing by
the speed of light (double speed).


6. Aberration:

    psAberration(obj, obj, direction, speed);

Convert obj to a psCube, check that:

    x=-0.35963388069046304
    y=0.5555192509816625
    z=0.7497078321908413


7. Generate a psTime (time) for 2003-04-01T01:30:00 UTC.  Check that:

time->sec == 1049160600
time->nsec == 0
time->leapsecond = false

Copy time and convert to UT1:

     timeUT1 = psTimeAlloc(PS_TIME_UTC);
     timeUT1->sec = time->sec;
     timeUT1->nsec = time->nsec;
     timeUT1->leapsecond = time->leapsecond;
     timeUT1 = psTimeConvert(timeUT1, PS_TIME_UT1);

Then check timeUT1:

     timeUT1->sec == 1049160599
     timeUT1->nsec == 657017200 (+/- hundreds of nsec at least)


8. Precession:

(a) psEarthPole *precession = psEOC_PrecessionModel(time);

Check that:

    X=2.857175590089105E-4
    Y=2.3968739377734732E-5
    S=-1.3970066457904322E-8

(b) psEarthPole *precessionCorr = psEOC_PrecessionCorr(time, PS_IERS_B);

Check that:

    X = 3.05224300720406e-10
    Y = -1.39441339235822e-10
    S = 0

(c) Generate the transformation:

    precession->x += precessionCorr->x;
    precession->y += precessionCorr->y;
    precession->s += precessionCorr->s;
    psSphereRot *precessionNutationInv = psSphereRot_CEOtoGCRS(precession); // This is CEO->GCRS
    psSphereRot precessionNutation = psSphereRotInvert(precessionNutationInv); // This is GCRS->CEO

Check that the rotation quaternion of precessionNutationInv is:

    -1.1984522406756289E-5
    1.4285893358610674E-4
    1.2191193518914336E-10
    -0.9999999897238481

(d) Apply the rotation to the position:

    obj = psSphereRotApply(obj, precessionNutation, obj);

Convert obj to a psCube and check that:

    x=-0.3598480726985338
    y=0.5555012823608123
    z=0.7496183628158023


9. Earth rotation

(a) Get the tidal correction:

    psEarthPole *tidalCorr = psEOC_PolarTideCorr(time);

(b) Generate the transformation:

    psSphereRot *earthRot = psSphereRot_TEOtoCEO(time, tidalCorr); // This is TEO->CEO

Check that the rotation quaternion of earthRot is:

    0.0
    0.0
    0.9625401009002903
    -0.2711393629830588

Within psSphereRot_TEOtoCEO, the Earth rotation angle should be 428251.4641536639 degrees

(c) Apply the rotation to the position:

    psSphereRot *earthRotInv = psSphereRotInvert(earthRot); // This is CEO->TEO
    obj = psSphereRotApply(obj, earthRotInv, obj);

Convert obj to a psCube and check that:

    x=0.01698625430807123
    y=-0.6616523084626379
    z=0.7496183628158023


10. Polar Motion:

(a) psEarthPole *earthPM = psEOC_GetPolarMotion(time, PS_IERS_B);
(b) psEarthPole *nutationCorr = psEOC_NutationCorr(time);

Add the above components:

    earthPM->x += nutationCorr->x;
    earthPM->y += nutationCorr->y;
    earthPM->s += nutationCorr->s;

Add the tidal correction:

    earthPM->x += tidalCorr->x;
    earthPM->y += tidalCorr->y;

Check that:

    x=-6.43607313124045e-07
    y=2.11351436973568e-06
    s=-7.39617581324646e-12

(c) Generate the transformation:

   psSphereRot *polarMotion = psSphereRot_ITRStoTEO(earthPM); // This is ITRS->TEO

Check that the rotation quaternion of polarMotion is:

    -1.0567571848664005E-6
    3.218036562931509E-7
    -3.3580195807204483E-12
    -0.9999999999993899
 
(d) Apply these to the position:

   psSphereRot *polarMotionInv = psSphereRotInvert(polarMotion); // This is TEO->ITRS
   obj = psSphereRotApply(obj, polarMotionInv, obj);

Convert obj to a psCube, should get:
   
    x=0.01698577185310146
    y=-0.6616538927902393
    z=0.7496169753347885


11. The psLib SDRS does not define functions to go beyond this point
(specifically, atmospheric refraction correction, and conversion to
the observer's horizon coordinates), so we stop here.  But check that
the above sequence can be executed faithfully in the reverse order as
well.
