 
Binary PCK Hands-On Lesson (FORTRAN)
===========================================================================
 
   March 01, 2023
 
 
Overview
--------------------------------------------------------
 
   In this lesson you will develop two programs that demonstrate geometric
   computations using ``high-accuracy'' Earth and Moon binary PCKs. The
   programs also demonstrate use of frame kernels and SPK files normally
   used together with these high-accuracy PCKs.
 
 
Note About HTML Links
--------------------------------------------------------
 
   The HTML version of this lesson contains links pointing to various HTML
   documents provided with the Toolkit. All of these links are relative
   and, in order to function, require this document to be in a certain
   location in the Toolkit HTML documentation directory tree.
 
   In order for the links to be resolved, if not done already by installing
   the lessons package under the Toolkit's ``doc/html'' directory, create a
   subdirectory called ``lessons'' under the ``doc/html'' directory of the
   ``toolkit/'' tree and copy this document to that subdirectory before
   loading it into a Web browser.
 
 
References
--------------------------------------------------------
 
   This section lists SPICE documents referred to in this lesson.
 
   Of these documents, the ``Tutorials'' contains the highest level
   descriptions with the least number of details while the ``Required
   Reading'' documents contain much more detailed specifications. The most
   complete specifications are provided in the ``API Documentation''.
 
 
Tutorials
 
   The following SPICE tutorials serve as references for the discussions in
   this lesson:
 
 
      Name              Lesson steps/routines it describes
      ----------------  -----------------------------------------------
      Frames            Moon rotation, Earth rotation
      PCK               Moon rotation, Earth rotation
      ``High Accuracy
      Orientation and
      Body-Fixed
      frames for Moon
      and Earth''
      (backup)          Moon rotation, Earth rotation
 
 
   These tutorials are available from the NAIF server at JPL:
 
      http://naif.jpl.nasa.gov/naif/tutorials.html
 
 
Required Readings
 
   The Required Reading documents are provided with the Toolkit and are
   located under the ``toolkit/doc'' directory in the SPICE Toolkit
   installation tree.
 
      Name             Lesson steps/routines that it describes
      ---------------  -----------------------------------------
      frames.req       Using reference frames
      pck.req          Obtaining planetary constants data
      spk.req          Obtaining ephemeris data
      time.req         Time conversion
 
 
The Permuted Index
 
   Another useful document distributed with the Toolkit is the permuted
   index. It is located under the ``toolkit/doc'' directory in the FORTRAN
   installation tree.
 
   This text document provides a simple mechanism by which users can
   discover which SPICE routines perform functions of interest, as well as
   the names of the source files that contain these routines. It is
   particularly useful for FORTRAN programmers because some of the routines
   are entry points; the names of these routines do not translate directly
   into the name of the respective source files that contain them.
 
 
API Documentation
 
   The most detailed specification of a given SPICE FORTRAN routine is
   contained in the header section of its source code. The source code is
   distributed with the Toolkit and is located under the
   ``toolkit/src/spicelib'' path.
 
   For example the path of the source code of the STR2ET routine is
 
      toolkit/src/spicelib/str2et.for
 
   Since some of the FORTRAN routines are entry points they may be part of
   a source file that has different name. The ``Permuted Index'' document
   mentioned above can be used to locate the name of their source file.
 
 
Kernels Used
--------------------------------------------------------
 
   The following kernels are used in examples provided in this lesson:
 
      #  FILE NAME                      TYPE DESCRIPTION
      -- ------------------------------ ---- ------------------------------
      1  naif0008.tls                   LSK  Generic LSK
      2  de414_2000_2020.bsp            SPK  Solar System Ephemeris
      3  moon_060721.tf                 FK   Lunar FK
      4  pck00008.tpc                   PCK  NAIF text PCK
      5  moon_pa_de403_1950-2198.bpc    PCK  Moon binary PCK
      6  earthstns_itrf93_050714.bsp    SPK  DSN station Ephemeris
      7  earth_topo_050714.tf           FK   Earth topocentric FK
      8  earth_000101_070725_070503.bpc PCK  Earth binary PCK
 
   These SPICE kernels are included in the lesson package.
 
 
SPICE Modules Used
--------------------------------------------------------
 
   This section provides a complete list of the routines and kernels that
   are suggested for usage in each of the exercises in this lesson. (You
   may wish to not look at this list unless/until you ``get stuck'' while
   working on your own.)
 
      CHAPTER EXERCISE   ROUTINES   FUNCTIONS  KERNELS
      ------- ---------  ---------  ---------  ----------
         1    mrotat     FURNSH     DPR        1-5
                         STR2ET     VSEP
                         SPKPOS     VDIST
                         RECLAT
                         SUBPNT
 
         2    erotat     FURNSH     DPR        1-2,4,6-8
                         STR2ET     VSEP
                         SPKPOS     SPD
                         RECLAT     TWOPI
                         TIMOUT     VDIST
                         PXFORM
                         SUBSLR
 
   Refer to the headers of the various routines listed above, as detailed
   interface specifications are provided with the source code.
 
 
Moon rotation (mrotat)
===========================================================================
 
 
Task Statement
--------------------------------------------------------
 
   Write a program that performs the following computations:
 
       1.   Convert the time string 2007 JAN 1 00:00:00 UTC to a double
            precision number representing seconds past J2000 TDB.
 
            In the following instructions, we'll call the result of this
            computation ET.
 
       2.   Compute the apparent position of the Earth as seen from the
            Moon in the IAU_MOON reference frame at the epoch ET. Use light
            time and stellar aberration corrections. Use RECLAT to compute
            the planetocentric longitude and latitude of the Earth position
            vector; display these coordinates in degrees.
 
       3.   Repeat the computation of step 2 using the MOON_ME reference
            frame. Display the results as above.
 
       4.   Compute the angular separation of the position vectors found in
            steps 2 and 3. Display the result in degrees.
 
       5.   Repeat the computation of step 2 using the MOON_PA reference
            frame. Display the results as above.
 
       6.   Compute the angular separation of the position vectors found in
            steps 3 and 5 (these vectors are expressed in the MOON_ME and
            MOON_PA frames). Display the result in degrees.
 
       7.   Compute the apparent sub-Earth point on the Moon at ET,
            expressed in the MOON_ME reference frame and using light time
            and stellar aberration corrections. Convert the sub-Earth point
            to latitudinal coordinates using RECLAT. Display the longitude
            and latitude of the sub-Earth point in degrees.
 
       8.   Repeat step 7, now using the MOON_PA frame.
 
       9.   Compute the distance between the two sub-Earth points found
            above in steps 7 and 8. Display the result in kilometers.
 
 
Learning Goals
--------------------------------------------------------
 
   Familiarity with SPICE kernels required to obtain high-accuracy
   orientation of the Moon. Understanding the differences between results
   obtained using low and high-accuracy Moon orientation data.
   Understanding the difference between the MOON_ME and MOON_PA frames.
 
 
Approach
--------------------------------------------------------
 
   The following ``tips'' may simplify the solution process.
 
       --   Examine the SPICE kernels provided with this lesson. Use BRIEF
            to find coverage periods of SPK kernels and binary PCKs. Use
            COMMNT to view the comment areas of binary PCKs. Examine text
            kernels, in particular text kernel comments, using a text
            editor or browser.
 
       --   Decide which SPICE kernels are necessary. Prepare a meta-kernel
            listing the kernels and load it into the program.
 
       --   Consult the above list titled ``SPICE Modules Used'' to see
            which routines are needed.
 
       --   The computational steps listed above should be followed in the
            order shown.
 
   You may find it useful to consult the permuted index, the headers of
   various source modules, and the tutorials titled ``PCK'' and ``High
   Accuracy Orientation and Body-Fixed frames for Moon and Earth.''
 
 
Solution
--------------------------------------------------------
 
 
Solution Meta-Kernel
 
   The meta-kernel we created for the solution to this exercise is named
   'mrotat.tm'. Its contents follow:
 
      KPL/MK
 
      Meta-kernel for the ``Moon Rotation'' task in the Binary PCK
      Hands On Lesson.
 
      The names and contents of the kernels referenced by this
      meta-kernel are as follows:
 
      File name                    Contents
      ---------------------------  ------------------------------------
      naif0008.tls                 Generic LSK
      de414_2000_2020.bsp          Solar System Ephemeris
      moon_060721.tf               Lunar FK
      pck00008.tpc                 NAIF text PCK
      moon_pa_de403_1950-2198.bpc  Moon binary PCK
 
      \begindata
 
         KERNELS_TO_LOAD = ( 'kernels/lsk/naif0008.tls'
                             'kernels/spk/de414_2000_2020.bsp'
                             'kernels/fk/moon_060721.tf'
                             'kernels/pck/pck00008.tpc'
                             'kernels/pck/moon_pa_de403_1950-2198.bpc' )
      \begintext
 
 
Solution Source Code
 
   A sample solution to the problem follows:
 
            PROGRAM MROTAT
            IMPLICIT NONE
      C
      C     SPICELIB functions
      C
            DOUBLE PRECISION      DPR
            DOUBLE PRECISION      VDIST
            DOUBLE PRECISION      VSEP
      C
      C     Local parameters
      C
            CHARACTER*(*)         METAKR
            PARAMETER           ( METAKR = 'mrotat.tm' )
 
            INTEGER               TIMLEN
            PARAMETER           ( TIMLEN = 35 )
      C
      C     Local variables
      C
            CHARACTER*(TIMLEN)    TIMSTR
 
            DOUBLE PRECISION      DIST
            DOUBLE PRECISION      ET
            DOUBLE PRECISION      IMOONV ( 3 )
            DOUBLE PRECISION      LAT
            DOUBLE PRECISION      LON
            DOUBLE PRECISION      LTIME
            DOUBLE PRECISION      MMOONV ( 3 )
            DOUBLE PRECISION      MSUB   ( 3 )
            DOUBLE PRECISION      PMOONV ( 3 )
            DOUBLE PRECISION      PSUB   ( 3 )
            DOUBLE PRECISION      R
            DOUBLE PRECISION      SEP
            DOUBLE PRECISION      SRFVEC ( 3 )
            DOUBLE PRECISION      TRGEPC
 
      C
      C     Load the kernels this program requires.
      C
            CALL FURNSH ( METAKR )
      C
      C     Convert our UTC string to seconds past J2000 TDB.
      C
            TIMSTR = '2007 JAN 1 00:00:00'
            CALL STR2ET ( TIMSTR, ET )
      C
      C     Look up the apparent position of the Earth relative
      C     to the Moon's center in the IAU_MOON frame at ET.
      C
            CALL SPKPOS ( 'EARTH', ET,     'IAU_MOON',
           .              'LT+S', 'MOON',  IMOONV,     LTIME )
      C
      C     Express the Earth direction in terms of longitude
      C     and latitude in the IAU_MOON frame.
      C
            CALL RECLAT ( IMOONV, R, LON, LAT )
 
            WRITE (*,'(A)') ' '
            WRITE (*,'(A)') 'Moon-Earth direction using low accuracy'
            WRITE (*,'(A)') 'PCK and IAU_MOON frame:'
            WRITE (*, '(A,F15.6)' ) 'Earth lon (deg):        ',
           .                           LON*DPR()
            WRITE (*, '(A,F15.6)' ) 'Earth lat (deg):        ',
           .                           LAT*DPR()
            WRITE (*,*) ' '
      C
      C     Look up the apparent position of the Earth relative
      C     to the Moon's center in the MOON_ME frame at ET.
      C
            CALL SPKPOS ( 'EARTH', ET,     'MOON_ME',
           .              'LT+S', 'MOON',  MMOONV,     LTIME )
      C
      C     Express the Earth direction in terms of longitude
      C     and latitude in the MOON_ME frame.
      C
            CALL RECLAT ( MMOONV, R, LON, LAT )
 
            WRITE (*,'(A)') 'Moon-Earth direction using high accuracy'
            WRITE (*,'(A)') 'PCK and MOON_ME frame:'
            WRITE (*, '(A,F15.6)' ) 'Earth lon (deg):        ',
           .                           LON*DPR()
            WRITE (*, '(A,F15.6)' ) 'Earth lat (deg):        ',
           .                           LAT*DPR()
            WRITE (*,*) ' '
      C
      C     Find the angular separation of the Earth position
      C     vectors in degrees.
      C
            SEP = DPR() * VSEP( MMOONV, IMOONV )
 
            WRITE (*, '(A)'       ) 'For IAU_MOON vs MOON_ME frames:'
            WRITE (*, '(A,F15.6)' )
           .'Moon-Earth vector separation angle (deg):     ', SEP
            WRITE (*,*) ' '
      C
      C     Look up the apparent position of the Earth relative
      C     to the Moon's center in the MOON_PA frame at ET.
      C
            CALL SPKPOS ( 'EARTH', ET,     'MOON_PA',
           .              'LT+S', 'MOON',  PMOONV,     LTIME )
      C
      C     Express the Earth direction in terms of longitude
      C     and latitude in the MOON_PA frame.
      C
            CALL RECLAT ( PMOONV, R, LON, LAT )
 
            WRITE (*,'(A)') 'Moon-Earth direction using high accuracy'
            WRITE (*,'(A)') 'PCK and MOON_PA frame:'
            WRITE (*, '(A,F15.6)' ) 'Earth lon (deg):        ',
           .                           LON*DPR()
            WRITE (*, '(A,F15.6)' ) 'Earth lat (deg):        ',
           .                           LAT*DPR()
            WRITE (*,*) ' '
      C
      C     Find the angular separation of the Earth position
      C     vectors in degrees.
      C
            SEP = DPR() * VSEP( MMOONV, PMOONV )
 
            WRITE (*, '(A)'       ) 'For MOON_PA vs MOON_ME frames:'
            WRITE (*, '(A,F15.6)' )
           .'Moon-Earth vector separation angle (deg):     ', SEP
            WRITE (*,*) ' '
      C
      C     Find the apparent sub-Earth point on the Moon at ET using
      C     the MOON_ME frame.
      C
            CALL SUBPNT ( 'NEAR POINT: ELLIPSOID',
           .              'MOON',  ET,   'MOON_ME', 'LT+S',
           .              'EARTH', MSUB, TRGEPC,    SRFVEC )
      C
      C     Display the sub-Earth point in latitudinal coordinates.
      C
            CALL RECLAT ( MSUB, R, LON, LAT )
 
            WRITE (*,'(A)')
           .'Sub-Earth point on Moon using high accuracy'
            WRITE (*,'(A)') 'PCK and MOON_ME frame:'
            WRITE (*, '(A,F15.6)' ) 'Sub-Earth lon (deg):   ',
           .                           LON *DPR()
            WRITE (*, '(A,F15.6)' ) 'Sub-Earth lat (deg):   ',
           .                           LAT*DPR()
            WRITE (*,*) ' '
      C
      C     Find the apparent sub-Earth point on the Moon at ET using
      C     the MOON_PA frame.
      C
            CALL SUBPNT ( 'NEAR POINT: ELLIPSOID',
           .              'MOON',  ET,   'MOON_PA', 'LT+S',
           .              'EARTH', PSUB, TRGEPC,    SRFVEC )
 
      C
      C     Display the sub-point in latitudinal coordinates.
      C
            CALL RECLAT ( PSUB, R, LON, LAT )
 
            WRITE (*,'(A)')
           .'Sub-Earth point on Moon using high accuracy'
            WRITE (*,'(A)') 'PCK and MOON_PA frame:'
            WRITE (*, '(A,F15.6)' ) 'Sub-Earth lon (deg):   ',
           .                           LON *DPR()
            WRITE (*, '(A,F15.6)' ) 'Sub-Earth lat (deg):   ',
           .                           LAT*DPR()
            WRITE (*,*) ' '
      C
      C     Find the distance between the sub-Earth points
      C     in km.
      C
            DIST = VDIST( MSUB, PSUB )
 
            WRITE (*, '(A,F15.6)' )
           .'Distance between sub-Earth points (km): ', DIST
 
            END
 
 
Solution Sample Output
 
   After compiling the program, execute it:
 
 
      Moon-Earth direction using low accuracy
      PCK and IAU_MOON frame:
      Earth lon (deg):               3.613102
      Earth lat (deg):              -6.438342
 
      Moon-Earth direction using high accuracy
      PCK and MOON_ME frame:
      Earth lon (deg):               3.611229
      Earth lat (deg):              -6.439501
 
      For IAU_MOON vs MOON_ME frames:
      Moon-Earth vector separation angle (deg):            0.002194
 
      Moon-Earth direction using high accuracy
      PCK and MOON_PA frame:
      Earth lon (deg):               3.593319
      Earth lat (deg):              -6.417582
 
      For MOON_PA vs MOON_ME frames:
      Moon-Earth vector separation angle (deg):            0.028235
 
      Sub-Earth point on Moon using high accuracy
      PCK and MOON_ME frame:
      Sub-Earth lon (deg):          3.611419
      Sub-Earth lat (deg):         -6.439501
 
      Sub-Earth point on Moon using high accuracy
      PCK and MOON_PA frame:
      Sub-Earth lon (deg):          3.593509
      Sub-Earth lat (deg):         -6.417582
 
      Distance between sub-Earth points (km):        0.856182
 
 
Earth rotation (erotat)
===========================================================================
 
 
Task Statement
--------------------------------------------------------
 
   Write a program that performs the following computations:
 
       1.   Convert the time string 2007 JAN 1 00:00:00 UTC to a double
            precision number representing seconds past J2000 TDB.
 
            In the following instructions, we'll call the result of this
            computation ET.
 
       2.   Compute the apparent position of the Moon as seen from the
            Earth in the IAU_EARTH reference frame at the epoch ET. Use
            light time and stellar aberration corrections. Display the
            planetocentric longitude and latitude of the Moon position
            vector in degrees.
 
       3.   Repeat the first computation using the ITRF93 reference frame.
            Display the results as above.
 
       4.   Compute the angular separation of the position vectors found
            the the previous two steps. Display the result in degrees.
 
   The following computations (steps 5-10) examine the cause of the angular
   offset found above, which is attributable to the rotation between the
   ITRF93 and IAU_EARTH frames. Steps 11 and up don't rely on the results
   of steps 5-10, so steps 5-10 may be safely skipped if they're not of
   interest to you.
 
   For each of the two epochs ET and ET + 100 days, examine the differences
   between the axes of the ITRF93 and IAU_EARTH frames using the following
   method:
 
       5.   Convert the epoch of interest to a string in the format style
            ``2007-MAY-16 02:29:00.000 (UTC).'' Display this string.
 
       6.   Look up the 3x3 position transformation matrix that converts
            vectors from the IAU_EARTH to the ITRF93 frame at the epoch of
            interest. We'll call the returned matrix RMAT.
 
       7.   Extract the first row of RMAT into a 3-vector, which we'll call
            ITRFX. This is the X-axis of the ITRF93 frame expressed
            relative to the IAU_EARTH frame.
 
       8.   Extract the third row of RMAT into a 3-vector, which we'll call
            ITRFZ. This is the Z-axis of the ITRF93 frame expressed
            relative to the IAU_EARTH frame.
 
       9.   Compute the angular separation between the vector ITRFX and the
            X-axis (1, 0, 0) of the IAU_EARTH frame. Display the result in
            degrees.
 
      10.   Compute the angular separation between the vector ITRFZ and the
            Z-axis (0, 0, 1) of the IAU_EARTH frame. Display the result in
            degrees.
 
   This is the end of the computations to be performed for the epochs ET
   and ET + 100 days. The following steps are part of a new computation.
 
   Find the azimuth and elevation of the apparent position of the Moon as
   seen from the DSN station DSS-13 by the following steps:
 
      11.   Find the apparent position vector of the Moon relative to the
            DSN station DSS-13 in the topocentric reference frame
            DSS-13_TOPO at epoch ET. Use light time and stellar aberration
            corrections.
 
            For this step, you'll need to have loaded a station SPK file
            providing geocentric station position vectors, as well as a
            frame kernel specifying topocentric reference frames centered
            at the respective DSN stations. (Other kernels will be needed
            as well; you must choose these.)
 
      12.   Convert the position vector to latitudinal coordinates. Use the
            routine RECLAT for this computation.
 
      13.   Compute the Moon's azimuth and elevation as follows: azimuth is
            the negative of topocentric longitude and lies within the range
            0-360 degrees; elevation is equal to the topocentric latitude.
            Display the results in degrees.
 
   The next computations demonstrate ``high-accuracy'' geometric
   computations using the Earth as the target body. These computations are
   *not* realistic; they are simply meant to demonstrate SPICE system
   features used for geometry computations involving the Earth as a target
   body. For example, the same basic techniques would be used to find the
   sub-solar point on the Earth as seen from an Earth-orbiting spacecraft.
 
      14.   Compute the apparent sub-solar point on the Earth at ET,
            expressed relative to the IAU_EARTH reference frame, using
            light time and stellar aberration corrections and using the Sun
            as the observer. Convert the sub-solar point to latitudinal
            coordinates using RECLAT. Display the longitude and latitude of
            the sub-solar point in degrees.
 
      15.   Repeat the sub-solar point computation described above, using
            the ITRF93 Earth body-fixed reference frame. Display the
            results as above.
 
      16.   Compute the distance between the two sub-solar points found
            above. Display the result in kilometers.
 
 
Learning Goals
--------------------------------------------------------
 
   Familiarity with SPICE kernels required to obtain high-accuracy
   orientation of the Earth. Understanding the differences between results
   obtained using low and high-accuracy Earth orientation data.
 
   Understanding of topocentric frames and computation of target geometry
   relative to a surface location on the Earth. Knowledge of SPICE kernels
   required to support such computations.
 
 
Approach
--------------------------------------------------------
 
   The following ``tips'' may simplify the solution process.
 
       --   Examine the SPICE kernels provided with this lesson. Use BRIEF
            to find coverage periods of SPK kernels and binary PCKs. Use
            COMMNT to view the comment areas of binary PCKs. Examine text
            kernels, in particular text kernel comments, using a text
            editor or browser.
 
       --   Decide which SPICE kernels are necessary. Prepare a meta-kernel
            listing the kernels and load it into the program.
 
       --   Consult the above list titled ``SPICE Modules Used'' to see
            which routines are needed. Note the functions used to provide
            the values ``seconds per day,'' ``degrees per radian,'' and ``2
            times Pi.''
 
       --   Examine the header of the routine RECLAT. Note that this
            routine may be used for coordinate conversions in situations
            where the input rectangular coordinates refer to any reference
            frame, not only a body-centered, body-fixed frame whose X-Y
            plane coincides with the body's equator.
 
       --   The computational steps listed above should be followed in the
            order shown, but steps 5-10 may be omitted.
 
   You may find it useful to consult the permuted index, the headers of
   various source modules, and the tutorials titled ``PCK'' and ``High
   Accuracy Orientation and Body-Fixed frames for Moon and Earth.''
 
 
Solution
--------------------------------------------------------
 
 
Solution Meta-Kernel
 
   The meta-kernel we created for the solution to this exercise is named
   'erotat.tm'. Its contents follow:
 
      KPL/MK
 
      Meta-kernel for the ``Earth Rotation'' task
      in the Binary PCK Hands On Lesson.
 
      The names and contents of the kernels referenced by this
      meta-kernel are as follows:
 
      File name                       Contents
      ------------------------------  ---------------------------------
      naif0008.tls                    Generic LSK
      de414_2000_2020.bsp             Solar System Ephemeris
      earthstns_itrf93_050714.bsp     DSN station Ephemeris
      earth_topo_050714.tf            Earth topocentric FK
      pck00008.tpc                    NAIF text PCK
      earth_000101_070725_070503.bpc  Earth binary PCK
 
 
      \begindata
 
      KERNELS_TO_LOAD = ( 'kernels/lsk/naif0008.tls'
                          'kernels/spk/de414_2000_2020.bsp'
                          'kernels/spk/earthstns_itrf93_050714.bsp'
                          'kernels/fk/earth_topo_050714.tf'
                          'kernels/pck/pck00008.tpc'
                          'kernels/pck/earth_000101_070725_070503.bpc' )
 
      \begintext
 
 
Solution Source Code
 
   A sample solution to the problem follows:
 
            PROGRAM EROTAT
            IMPLICIT NONE
 
      C
      C     SPICELIB functions
      C
            DOUBLE PRECISION      DPR
            DOUBLE PRECISION      SPD
            DOUBLE PRECISION      TWOPI
            DOUBLE PRECISION      VDIST
            DOUBLE PRECISION      VSEP
 
      C
      C     Local parameters
      C
            CHARACTER*(*)         METAKR
            PARAMETER           ( METAKR = 'erotat.tm' )
 
            INTEGER               TIMLEN
            PARAMETER           ( TIMLEN = 35 )
 
      C
      C     Local variables
      C
            CHARACTER*(TIMLEN)    TIMSTR
            CHARACTER*(TIMLEN)    OUTSTR
 
            DOUBLE PRECISION      AZ
            DOUBLE PRECISION      DIST
            DOUBLE PRECISION      EL
            DOUBLE PRECISION      ET
            DOUBLE PRECISION      HMOONV ( 3 )
            DOUBLE PRECISION      HSUB   ( 3 )
            DOUBLE PRECISION      ITRFX  ( 3 )
            DOUBLE PRECISION      ITRFZ  ( 3 )
            DOUBLE PRECISION      LMOONV ( 3 )
            DOUBLE PRECISION      LAT
            DOUBLE PRECISION      LON
            DOUBLE PRECISION      LSUB   ( 3 )
            DOUBLE PRECISION      LTIME
            DOUBLE PRECISION      R
            DOUBLE PRECISION      RMAT   ( 3, 3 )
            DOUBLE PRECISION      SEP
            DOUBLE PRECISION      SRFVEC ( 3 )
            DOUBLE PRECISION      T
            DOUBLE PRECISION      TOPOV  ( 3 )
            DOUBLE PRECISION      TRGEPC
            DOUBLE PRECISION      X      ( 3 )
            DOUBLE PRECISION      Z      ( 3 )
 
            INTEGER               I
            INTEGER               J
 
            DATA                  X / 1.D0, 0.D0, 0.D0 /
            DATA                  Z / 0.D0, 0.D0, 1.D0 /
 
      C
      C     Load the kernels this program requires.
      C
            CALL FURNSH ( METAKR )
 
      C
      C     Convert our UTC string to seconds past J2000 TDB.
      C
            TIMSTR = '2007 JAN 1 00:00:00'
            CALL STR2ET ( TIMSTR, ET )
 
      C
      C     Look up the apparent position of the Moon relative
      C     to the Earth's center in the IAU_EARTH frame at ET.
      C
            CALL SPKPOS ( 'MOON', ET,       'IAU_EARTH',
           .              'LT+S', 'EARTH',  LMOONV,     LTIME )
 
      C
      C     Express the Moon direction in terms of longitude
      C     and latitude in the IAU_EARTH frame.
      C
            CALL RECLAT ( LMOONV, R, LON, LAT )
 
            WRITE (*,'(A)') ' '
            WRITE (*,'(A)') 'Earth-Moon direction using low accuracy'
            WRITE (*,'(A)') 'PCK and IAU_EARTH frame:'
            WRITE (*, '(A,F15.6)' ) 'Moon lon (deg):        ',
           .                           LON*DPR()
            WRITE (*, '(A,F15.6)' ) 'Moon lat (deg):        ',
           .                           LAT*DPR()
            WRITE (*,*) ' '
 
      C
      C     Look up the apparent position of the Moon relative
      C     to the Earth's center in the ITRF93 frame at ET.
      C
            CALL SPKPOS ( 'MOON', ET,       'ITRF93',
           .              'LT+S', 'EARTH',  HMOONV,     LTIME )
 
      C
      C     Express the Moon direction in terms of longitude
      C     and latitude in the ITRF93 frame.
      C
            CALL RECLAT ( HMOONV, R, LON, LAT )
 
            WRITE (*,'(A)') 'Earth-Moon direction using high accuracy'
            WRITE (*,'(A)') 'PCK and ITRF93 frame:'
            WRITE (*, '(A,F15.6)' ) 'Moon lon (deg):        ',
           .                           LON*DPR()
            WRITE (*, '(A,F15.6)' ) 'Moon lat (deg):        ',
           .                           LAT*DPR()
            WRITE (*,*) ' '
 
      C
      C     Find the angular separation of the Moon position
      C     vectors in degrees.
      C
            SEP = DPR() * VSEP( LMOONV, HMOONV )
 
            WRITE (*, '(A,F15.6)' )
           .'Earth-Moon vector separation angle (deg):     ', SEP
            WRITE (*,*) ' '
 
      C
      C     Next, express the +Z and +X axes of the ITRF93 frame in
      C     the IAU_EARTH frame. We'll do this for two times: ET
      C     and ET + 100 days.
      C
            DO I = 1, 2
      C
      C        Set the time, expressing the time delta in
      C        seconds.
      C
               T = ET + (I-1)*SPD()*100
      C
      C        Convert the TDB time T to a string for output.
      C
               CALL TIMOUT ( T,
           .                 'YYYY-MON-DD HR:MN:SC.### (UTC)',
           .                 OUTSTR )
 
               WRITE (*,'(2A)') 'Epoch: ', OUTSTR
 
      C
      C        Find the rotation matrix for conversion of position
      C        vectors from the IAU_EARTH to the ITRF93 frame.
      C
               CALL PXFORM ( 'IAU_EARTH', 'ITRF93', T, RMAT )
 
               DO J = 1, 3
                  ITRFX(J) = RMAT(1,J)
                  ITRFZ(J) = RMAT(3,J)
               END DO
 
      C
      C        Display the angular offsets of the ITRF93 +X and +Z
      C        axes from their IAU_EARTH counterparts.
      C
               SEP = VSEP ( ITRFX, X )
               WRITE (*, '(A,F13.6)' )
           .         'ITRF93 - IAU_EARTH +X axis separation angle '
           .   //    '(deg): ', SEP * DPR()
 
               SEP = VSEP ( ITRFZ, Z )
               WRITE (*, '(A,F13.6)' )
           .         'ITRF93 - IAU_EARTH +Z axis separation angle '
           .   //    '(deg): ', SEP * DPR()
               WRITE (*,*) ' '
 
            END DO
 
      C
      C     Find the azimuth and elevation of the apparent
      C     position of the Moon in the local topocentric
      C     reference frame at the DSN station DSS-13.
      C     First look up the Moon's position relative to the
      C     station in that frame.
      C
            CALL SPKPOS ( 'MOON', ET,       'DSS-13_TOPO',
           .              'LT+S', 'DSS-13',  TOPOV,       LTIME )
 
      C
      C     Express the station-Moon direction in terms of
      C     longitude and latitude in the DSS-13_TOPO frame.
      C
            CALL RECLAT ( TOPOV, R, LON, LAT )
      C
      C     Convert to azimuth/elevation.
      C
            AZ = -LON
 
            IF ( AZ .LT. 0 ) THEN
               AZ = AZ + TWOPI()
            END IF
 
            EL = LAT
 
            WRITE (*,'(A)')
           .   'DSS-13-Moon az/el using high accuracy '
           .// 'PCK and DSS-13_TOPO frame:'
            WRITE (*, '(A,F15.6)' ) 'Moon Az (deg):         ',
           .                           AZ * DPR()
            WRITE (*, '(A,F15.6)' ) 'Moon El (deg):         ',
           .                           EL * DPR()
            WRITE (*,*) ' '
 
      C
      C     Find the sub-solar point on the Earth at ET using the
      C     Earth body-fixed frame IAU_EARTH. Treat the Sun as
      C     the observer.
      C
            CALL SUBSLR ( 'NEAR POINT: ELLIPSOID',
           .              'EARTH', ET,   'IAU_EARTH', 'LT+S',
           .              'SUN',   LSUB, TRGEPC,       SRFVEC )
 
      C
      C     Display the sub-point in latitudinal coordinates.
      C
            CALL RECLAT ( LSUB, R, LON, LAT )
 
            WRITE (*,'(A)')
           .'Sub-Solar point on Earth using low accuracy'
            WRITE (*,'(A)') 'PCK and IAU_EARTH frame:'
            WRITE (*, '(A,F15.6)' ) 'Sub-Solar lon (deg):   ',
           .                           LON *DPR()
            WRITE (*, '(A,F15.6)' ) 'Sub-Solar lat (deg):   ',
           .                           LAT*DPR()
            WRITE (*,*) ' '
 
      C
      C     Find the sub-solar point on the Earth at ET using the
      C     Earth body-fixed frame ITRF93. Treat the Sun as
      C     the observer.
      C
            CALL SUBSLR ( 'NEAR POINT: ELLIPSOID',
           .              'EARTH', ET,   'ITRF93', 'LT+S',
           .              'SUN',   HSUB, TRGEPC,   SRFVEC )
 
      C
      C     Display the sub-point in latitudinal coordinates.
      C
            CALL RECLAT ( HSUB, R, LON, LAT )
 
            WRITE (*,'(A)')
           .'Sub-Solar point on Earth using high accuracy'
            WRITE (*,'(A)') 'PCK and ITRF93 frame:'
            WRITE (*, '(A,F15.6)' ) 'Sub-Solar lon (deg):   ',
           .                           LON*DPR()
            WRITE (*, '(A,F15.6)' ) 'Sub-Solar lat (deg):   ',
           .                           LAT*DPR()
            WRITE (*,*) ' '
 
      C
      C     Find the distance between the sub-solar point
      C     vectors in km.
      C
            DIST = VDIST( LSUB, HSUB )
 
            WRITE (*, '(A,F15.6)' )
           .'Distance between sub-solar points (km): ', DIST
 
            END
 
 
Solution Sample Output
 
   After compiling the program, execute it:
 
 
      Earth-Moon direction using low accuracy
      PCK and IAU_EARTH frame:
      Moon lon (deg):             -35.496272
      Moon lat (deg):              26.416959
 
      Earth-Moon direction using high accuracy
      PCK and ITRF93 frame:
      Moon lon (deg):             -35.554286
      Moon lat (deg):              26.419156
 
      Earth-Moon vector separation angle (deg):            0.052002
 
      Epoch: 2007-JAN-01 00:00:00.000 (UTC)
      ITRF93 - IAU_EARTH +X axis separation angle (deg):      0.057677
      ITRF93 - IAU_EARTH +Z axis separation angle (deg):      0.002326
 
      Epoch: 2007-APR-10 23:59:59.998 (UTC)
      ITRF93 - IAU_EARTH +X axis separation angle (deg):      0.057787
      ITRF93 - IAU_EARTH +Z axis separation angle (deg):      0.002458
 
      DSS-13-Moon az/el using high accuracy PCK and DSS-13_TOPO frame:
      Moon Az (deg):               72.169006
      Moon El (deg):               20.689488
 
      Sub-Solar point on Earth using low accuracy
      PCK and IAU_EARTH frame:
      Sub-Solar lon (deg):       -177.100531
      Sub-Solar lat (deg):        -22.910377
 
      Sub-Solar point on Earth using high accuracy
      PCK and ITRF93 frame:
      Sub-Solar lon (deg):       -177.157874
      Sub-Solar lat (deg):        -22.912593
 
      Distance between sub-solar points (km):        5.881861
 
