1. Start with RA, Dec:

    r = 122.9153182445501
    d = 48.562968978679194

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

2. Convert obj to a psCube (psSphereToCube), should get:

    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).

4. Gravitational deflection:

    psGravityDeflection(obj, obj, sun);

Convert obj to a psCube, should get:

    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 need to be summed.  Converting to a psSphere gives the
direction of motion (psSphere *direction), while the magnitude gives
the speed (double speed).

6. Aberration:

    psAberration(obj, obj, direction, speed);

Convert obj to a psCube, should get:

    x=-0.35963388069046304
    y=0.5555192509816625
    z=0.7497078321908413

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

8. Precession

(a) psEarthPole *precession = psEOC_PrecessionModel(time);
(b) psEarthPole *precessionCorr = psEOC_PrecessionCorr(time, PS_IERS_A);

Add the above components:

    precession->x += precessionCorr->x;
    precession->y += precessionCorr->y;
    precession->s += precessionCorr->s;

Apply these to the position:

    psSphereRot *precessionNutationInv = psSphereRot_CEOtoGCRS(earth); // This is CEO->GCRS
    psSphereRot precessionNutation = psSphereRotInvert(precessionNutationInv); // This is GCRS->CEO
    obj = psSphereRotApply(obj, precessionNutation, obj);

Convert obj to a psCube, should get:

    x=-0.3598480726985338
    y=0.5555012823608123
    z=0.7496183628158023

9. Earth rotation:

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

Convert obj to a psCube, should get:

    x=0.01698625430807123
    y=-0.6616523084626379
    z=0.7496183628158023

10. Polar Motion:

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

Add the above components:

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

Apply these to the position:

   psSphereRot *polarMotionInv = psSphereRot_ITRStoTEO(earthPM); // This is ITRS->TEO
   psSphereRot *polarMotion = psSphereRotInvert(polarMotionInv); // This is TEO->ITRS
   obj = psSphereRotApply(obj, polarMotion, 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.
