Geometric Event Finding Hands-On Lesson (FORTRAN)





March 30, 2005

Updated to correct minor bugs in Fortran utilities: umbrella functions MODRIS and MODOCP, as well as entries in these functions SETRIS and SETOCP, now return the value .FALSE. Fortran callers of SETRIS and SETOCP were updated to use function call syntax.

June 1, 2004

Updated to include code to test routines written by student: SRFAZL and OCCPT.

May 7, 2004



Overview




This lesson illustrates how the SPICE Toolkit can be used to find time intervals when specified geometric conditions are satisfied. The student is asked to construct a program that finds the times within a specified time range when a point target is visible from a specified surface point on an extended body. As a concrete example, the student is asked to find the complete set of times within the interval

   2004 MAY 1 TDB
   2004 MAY 5 TDB
when the Mars Express Orbiter is visible from the DSN station DSS-14. The program is to consider occultation of the spacecraft by Mars and by the earth. We'll consider the spacecraft occulted by the earth if its elevation is below six degrees when measured in the topocentric frame DSS-14_TOPO, which is centered at DSS-14.

Although we've identified a specific set of objects participating in our observation geometry problem, the solution program should be written in such a way that these selections could be changed simply by modifying a meta-kernel read by the program.

In this lesson, we use the terms ``event finding'' or ``root finding'' to designate the process of solving for the set of times when a specified geometric condition is met. Rather than ask the student to devise event finding algorithms, we provide the student with subroutines to perform this task. See the section titled ``Utility routines'' below for details.

The utility routines provided here are not part of the SPICE Toolkit. NAIF plans to include event finding software in a future Toolkit release. The planned SPICE Toolkit subsystem dealing with event finding is called the ``Geometry Finder.'' The Geometry Finder's subroutine API will not closely resemble the event finding interface presented here. However, the event finding utilities comprise a small part of the overall set of routines the student will use in this lesson. Aside from these few utilities, the other routines presented in this lesson are in common use among SPICELIB programmers.

One set of SPICELIB routines that are particularly relevant to event finding are the SPICELIB window manipulation routines. In addition to this lesson, the ``Other Stuff'' hands-on programming lesson contains an exercise dealing with SPICELIB window routines. The student may wish to complete that exercise in addition to this lesson.

This lesson consists of a series of programming exercises. In each exercise, the student writes a portion of the program constituting the complete problem solution. The exercise topics are:

This document contains a chapter for each of the listed topics.



References




This section provides a list of SPICE documents that are referred to in this lesson.

Of these documents, the ``Tutorials'' contains the highest level descriptions with the fewest details while the ``Required Reading'' documents contain much more detailed specifications. The most complete specifications are provided in the ``Headers'' -- the comments in the top section of each SPICELIB source file.



Tutorials



The following SPICE tutorials may be particularly useful as references for this lesson:

   Name                    Lesson steps/routines that it describes
   ---------------         -----------------------------------------
   Time                    Conversion between UTC and ET
   Intro to Kernel Files   Loading SPICE kernels
   SPK                     Computing positions and velocities
   Frames                  Computing transformations between frames
These tutorials are available in printed form and as MS Office or PDF files from the NAIF server at JPL:

   ftp://naif.jpl.nasa.gov/pub/naif/toolkit_docs/Tutorials


Required Reading Documents



The Required Reading documents are provided with the Toolkit and are located under ``toolkit/doc'' or ``cspice/doc'' directory in the FORTRAN and C Toolkit installation trees.

   Name             Functions/routines that it describes
   ---------------  -----------------------------------------
   cells.req        The SPICE cell data type
   error.req        SPICE error handling
   frames.req       Using reference frames
   kernel.req       Loading SPICE kernels
   naif_ids.req     Body and reference frame names
   sets.req         The SPICE set data type
   spk.req          Computing positions and velocities
   time.req         UTC to ET time conversion
   windows.req      The SPICE window data type
Another very useful document, also distributed with the Toolkit, is the ``Permuted Index'', called ``spicelib.idx''; the Permuted index is located under ``doc'' directory in the Toolkit installation tree.

This text document provides an easy way to find which SPICE routine(s) performs a particular function in which you are interested. It also provides the name of the source file that contains this routine. This is especially useful for FORTRAN because some of the routines are entry points and, therefore, their names differ those of the source files in which they're defined.



Source Code Header Comments



The most detailed specification of a given SPICE FORTRAN or C routine is contained in the header section of its source code. The source code is distributed with the Toolkit and is located under ``toolkit/src/spicelib'' in FORTRAN and under ``cspice/src/cspice'' in C Toolkits.

For example the source code of the STR2ET/str2et_c routine is

   toolkit/src/spicelib/str2et.for
in the FORTRAN Toolkit and in

   cspice/src/cspice/str2et_c.c
in the C Toolkit.

Since some of the FORTRAN routines are entry points they are usually 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
   -----------------------         ----  --------------------------
   de405s.bsp                      SPK   Planet Ephemeris SPK
   dsnstns.bsp                     SPK   DSN station SPK
   DSN_topo.frm                    FK    DSN station frame definitions
   earth_031228_231229_predict.bpc PCK   Binary predict PCK for earth
   ORMM__040501000000_00069.BSP    SPK   MEX Orbiter Trajectory SPK
   naif0007.tls                    LSK   Generic LSK
   pck00007.tpc                    PCK   Generic PCK
These SPICE kernels are available from the NAIF server at JPL:

   ftp://naif.jpl.nasa.gov/pub/naif/misc/lessons/events


SPICE Routines Used




The following SPICE routines are used in the solution programs to be written by the student, or in the supporting utility routines provided in this lesson:

   Name        Function that it performs
   ----------  ---------------------------------------------------
   APPNDD      Appends an item to a d.p. cell
   BADKPV      Checks a kernel variable definition
   BRCKTD      Bracket a d.p. value
   CARDD       Return cardinality of d.p. cell
   CHKIN       SPICE call tree tracing check-in
   CHKOUT      SPICE call tree tracing check-out
   CNMFRM      Map a body name to an associated frame
   RPD         Return number of degrees per radian
   ERRCH       Substitute a character string into an error message
   ERRDP       Substitute a d.p. number into an error message
   FURNSH      Loads kernels, individually or listed in meta-kernel
   GCPOOL      Fetch string values from the kernel pool
   GDPOOL      Fetch d.p. values from the kernel pool
   KDATA       Return info on specified loaded kernel
   KTOTAL      Return the count of loaded kernels of a given type
   RECLAT      Transform from rectangular to latitudinal coordinates
   REPMC       Substitute a substring for a marker in a string
   REPMD       Substitute double precision value for marker in string
   RETURN      Test whether SPICELIB routines should return on entry
   RPD         Return number of radians per degree
   SETMSG      Set the SPICE long error message
   SIGERR      Signal a SPICE error
   SPKPOS      Computes positions of ephemeris objects
   SSIZED      Set the size of a d.p. cell
   STR2ET      Converts a time string to ET seconds past J2000
   SRFXPT      Find intersection of ray and ellipsoidal target body
   TIMOUT      Format a time string for output
   TOSTDO      Write a string to standard output
   VNORM       Find the norm of a 3-vector
   WNDIFD      Find the difference of two d.p. windows
   WNFETD      Fetch a specified interval from a d.p. window
   WNINSD      Insert an interval into a d.p. window
The most detailed documentation source for these routines is their headers.



Utility routines not included in SPICELIB




The following utility routines are used in the solution programs to be written by the student, or in the supporting utility routines provided in this lesson:

   Name        Function that it performs
   ----------  ---------------------------------------------------
   FNDEVT      Find set of times of state transitions
   MAKWIN      Convert state transition set to SPICELIB window
   MODOCP      Umbrella routine for occultation utilities
   MODRIS      Umbrella routine for rise-set utilities
   SETOCP      Set parameter values for occultation utilities
   SETRIS      Set parameter values for rise-set utilities
   STOCP       Indicate whether target is occulted
   STRIS       Indicate whether target elevation is above limit
The subroutine FNDEVT is the key root-finding utility. Given a LOGICAL function of time F(T), and a time interval, FNDEVT finds the set of times within the interval at which the value returned by F(T) changes.

In this context, we call F(T) a ``state function'' and the epochs at which the value of F(T) changes ``epochs of state transitions'' or simply ``state transitions.''

FNDEVT stores the state transition epochs it finds in a SPICELIB set, which is an output argument.

The utility MAKWIN converts the set returned by FNDEVT into a SPICELIB window whose intervals represent the time periods when the state returned by F(T) is .TRUE.

The purpose of the suite of routines contained in the module MODRIS is to enable a user to define and use a ``state function'' that indicates when a target is above a specified elevation limit, as seen by an observer at a surface location on an extended body. The module contains an initialization function SETRIS that allows a user to set parameters defining the geometry of interest and a state function STRIS that uses values established by SETRIS to decide whether the target is above the elevation limit at a specified time. The function STRIS computes target elevation using the subroutine SRFAZL, which is supplied by the student.

In the same vein, the purpose of the suite of routines contained in the module MODOCP is to enable a user to define and use a ``state function'' that indicates when a point target is occulted by a specified extended body, as seen by a specified observer. The module contains an initialization function SETOCP that allows a user to set parameters defining the geometry of interest and a state function STOCP that uses values established by SETOCP to decide whether the target is occulted at a specified time. The function STOCP tests for occultation using the subroutine OCCPT, which is supplied by the student.



Input Data







Task Statement




Write a program that loads SPICE kernels and reads input values required to solve the point target visibility problem.

Create a SPICE text kernel (``meta-kernel'') containing information required to initialize the program. The program is to read both command inputs and SPICE kernel names from this SPICE text kernel. The following items are to be read from the kernel:

The program should perform error checking on the inputs obtained from the kernel file. Each item, other than the kernels to be loaded, should be verified to be present and to have the correct data type. To test the program, use the following inputs:

Finally, the program should write to standard output the kernel variable inputs read from the text kernel. Times should be displayed in both TDB and UTC calendar format. The items displayed should include the names and types of the loaded SPICE kernels.



Learning Goals




The student is introduced to flexible, robust techniques for introducing data into SPICELIB based programs.

The student is exposed to the following SPICE components:



Approach






Solution steps



A possible solution would consist of the following steps:

Preparation:

Next, write a program that performs the following steps:

You may find it useful to consult the permuted index, the headers of various source modules, and the ``Time Required Reading'' and ``Kernel Required Reading'' documents.



Solution






Solution Meta-Kernel



The meta-kernel we created for the solution to this exercise is named 'geomevnt.mk'. Its contents follow:

 
 
      Example meta-kernel for geometric event finding hands-on
      coding lesson.
 
         Version 1.0.0 07-MAY-2004 (NJB)
 
 
      Identify names of kernels to load:
 
   \begindata
 
      KERNELS_TO_LOAD = (
 
              'kernels/spk/de405s.bsp'
              'kernels/spk/dsnstns.bsp'
              'kernels/fk/DSN_topo.frm'
              'kernels/pck/earth_031228_231229_predict.bpc'
              'kernels/lsk/naif0007.tls'
              'kernels/spk/ORMM__040501000000_00069.BSP'
              'kernels/pck/pck00007.tpc'
                        )
 
   \begintext
 
       Set values of other inputs:
 
           surface point
           target
           occulting body
           aberration correction
           search interval start time
           search interval stop time
           elevation limit
 
      The kernel variable names used here are examples; the
      programmer may choose other names.  The names must be no
      longer than 32 characters and are case-sensitive.  These
      names must be referenced in calls to the kernel pool fetch
      routines.
 
   \begindata
 
      SCVIEW_SRFPT  = 'DSS-14'
      SCVIEW_TARGET = 'MEX'
      SCVIEW_OCCBDY = 'MARS'
      SCVIEW_ABCORR = 'CN+S'
      SCVIEW_START  = '2004 MAY 1 TDB'
      SCVIEW_STOP   = '2004 MAY 5 TDB'
      SCVIEW_ELVLIM =  6.0
 
   \begintext


Solution Code



The example program below shows one possible solution.

         PROGRAM SCVIEW
         IMPLICIT NONE
   C
   C     Load SPICE Kernels and read additional inputs required to
   C     solve the spacecraft visibility problem posed by the
   C     geometric event finding hands-on code lesson.
   C
   C     The following items are loaded from a SPICE text kernel:
   C
   C        - Observer (surface point) name
   C        - Target name
   C        - Occulting body name
   C        - Aberration correction specification
   C        - Start time string
   C        - Stop time string
   C        - Elevation limit (in degrees)
   C
   C     This program assumes all necessary kernels and kernel
   C     variable definitions are specified in a SPICE meta-kernel.
   C     The meta-kernel name is assumed to be 'geomevnt.mk'.
   C
 
   C
   C     SPICELIB functions
   C
         DOUBLE PRECISION      RPD
         LOGICAL               BADKPV
 
   C
   C     Other functions
   C
 
   C
   C     Local parameters
   C
 
   C
   C     The meta-kernel:
   C
         CHARACTER*(*)         META
         PARAMETER           ( META   = 'geomevnt.mk' )
 
   C
   C     Labels of kernel pool items:
   C
         CHARACTER*(*)         OCCLBL
         PARAMETER           ( OCCLBL = 'SCVIEW_OCCBDY' )
 
         CHARACTER*(*)         SRFLBL
         PARAMETER           ( SRFLBL = 'SCVIEW_SRFPT' )
 
         CHARACTER*(*)         TRGLBL
         PARAMETER           ( TRGLBL = 'SCVIEW_TARGET' )
 
         CHARACTER*(*)         CORLBL
         PARAMETER           ( CORLBL = 'SCVIEW_ABCORR' )
 
         CHARACTER*(*)         STRLBL
         PARAMETER           ( STRLBL = 'SCVIEW_START' )
 
         CHARACTER*(*)         STPLBL
         PARAMETER           ( STPLBL = 'SCVIEW_STOP' )
 
         CHARACTER*(*)         ELVLBL
         PARAMETER           ( ELVLBL = 'SCVIEW_ELVLIM' )
 
   C
   C     String length parameters:
   C
         INTEGER               CORLEN
         PARAMETER           ( CORLEN = 10 )
 
         INTEGER               FILSIZ
         PARAMETER           ( FILSIZ = 255 )
 
         INTEGER               LNSIZE
         PARAMETER           ( LNSIZE = 200 )
 
         INTEGER               NAMLEN
         PARAMETER           ( NAMLEN = 32 )
 
         INTEGER               TIMLEN
         PARAMETER           ( TIMLEN = 50 )
 
         INTEGER               TYPLEN
         PARAMETER           ( TYPLEN = 50 )
 
   C
   C     Local variables
   C
         CHARACTER*(CORLEN)    ABCORR
         CHARACTER*(FILSIZ)    FILE
         CHARACTER*(TYPLEN)    FILTYP
         CHARACTER*(LNSIZE)    LINE
         CHARACTER*(NAMLEN)    OCCBDY
         CHARACTER*(FILSIZ)    SOURCE
         CHARACTER*(NAMLEN)    SRFPT
         CHARACTER*(TIMLEN)    START
         CHARACTER*(TIMLEN)    STOP
         CHARACTER*(NAMLEN)    TARGET
         CHARACTER*(TIMLEN)    TIMSTR
         CHARACTER*(LNSIZE)    TITLE
 
 
         DOUBLE PRECISION      ELVLIM
         DOUBLE PRECISION      ETBEG
         DOUBLE PRECISION      ETEND
         DOUBLE PRECISION      REVLIM
 
         INTEGER               HANDLE
         INTEGER               I
         INTEGER               N
         INTEGER               NKER
 
         LOGICAL               BAD
         LOGICAL               FOUND
 
   C
   C     Load the meta-kernel.
   C
         CALL FURNSH ( META )
 
   C
   C     Look up observation location, target,  name of the
   C     occulting body, aberration correction, start
   C     and stop times, and elevation limit.
   C
   C     We simplify error checking by using BADKPV.  BADKPV
   C     tests whether a kernel variable having specified
   C     attributes is present in the kernel pool. The fourth
   C     argument of BADKPV is the expected dimension; the
   C     last indicates the expected data type.  Read the
   C     header of BADKPV for details.
   C
         BAD = BADKPV ( 'SCVIEW', SRFLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', TRGLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', OCCLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', CORLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', STRLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', STPLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', ELVLBL, '=',  1,      1, 'N' )
 
   C
   C     Now we know the kernel variables of interest have been
   C     defined; look up the values.  We don't need to check
   C     the FOUND flag because BADKPV has ensured that each item
   C     will be found.
   C
   C     Note that the elevation limit has units of degrees on input,
   C     so we store the input in the variable ELVLIM.  We convert
   C     this value to radians for computation using SPICE routines;
   C     we'll store the equivalent value in radians in REVLIM.
   C
         CALL GCPOOL ( SRFLBL, 1, 1,      N,    SRFPT,  FOUND )
         CALL GCPOOL ( TRGLBL, 1, 1,      N,    TARGET, FOUND )
         CALL GCPOOL ( OCCLBL, 1, 1,      N,    OCCBDY, FOUND )
         CALL GCPOOL ( CORLBL, 1, 1,      N,    ABCORR, FOUND )
         CALL GCPOOL ( STRLBL, 1, 1,      N,    START,  FOUND )
         CALL GCPOOL ( STPLBL, 1, 1,      N,    STOP,   FOUND )
         CALL GDPOOL ( ELVLBL, 1, 1,      N,    ELVLIM, FOUND )
 
         REVLIM = RPD() * ELVLIM
 
   C
   C     Display to standard output a banner for the output report:
   C
         TITLE  = 'Inputs for geometric event finding program'
         CALL TOSTDO ( ' '   )
         CALL TOSTDO ( TITLE )
         CALL TOSTDO ( ' '   )
 
         LINE   = '   Target                       = #'
         CALL REPMC  ( LINE, '#', TARGET, LINE )
         CALL TOSTDO ( LINE )
 
         LINE   = '   Observation surface location = #'
         CALL REPMC  ( LINE, '#', SRFPT, LINE )
         CALL TOSTDO ( LINE )
 
         LINE   = '   Aberration correction        = #'
         CALL REPMC  ( LINE, '#', ABCORR, LINE )
         CALL TOSTDO ( LINE )
 
   C
   C     Convert the start and stop times to ET.
   C
         CALL STR2ET ( START, ETBEG )
         CALL STR2ET ( STOP,  ETEND )
 
   C
   C     Display the start time as both calendar UTC and
   C     calendar ET using the formats shown below.
   C
   C        2004 MAY 06 20:15:00.000 (UTC)
   C        2004 MAY 06 20:15:00.000 (TDB)
   C
         CALL TIMOUT ( ETBEG,
        .              'YYYY MON DD HR:MN:SC.### (UTC)',
        .              TIMSTR                           )
 
         LINE   = '   Start time                   = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
 
         CALL TIMOUT ( ETBEG,
        .              'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .              TIMSTR                           )
 
         LINE   = '   Start time                   = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
 
         CALL TIMOUT ( ETEND,
        .              'YYYY MON DD HR:MN:SC.### (UTC)',
        .              TIMSTR                           )
 
         LINE   = '   Stop time                    = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
 
         CALL TIMOUT ( ETEND,
        .              'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .              TIMSTR                                   )
 
         LINE   = '   Stop time                    = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
   C
   C     Display the elevation limit in degrees.
   C
         LINE   = '   Elevation limit (degrees)    = #'
         CALL REPMF  ( LINE, '#', ELVLIM, 7, 'F', LINE )
         CALL TOSTDO ( LINE )
 
   C
   C     Display the names of the SPICE kernels we've loaded.
   C
         CALL TOSTDO ( ' '   )
         CALL TOSTDO ( 'Loaded SPICE Kernels: ' )
 
         CALL KTOTAL ( 'ALL', NKER )
 
         DO I = 1, NKER
 
            CALL KDATA ( I,      'ALL',  FILE, FILTYP,
        .                SOURCE, HANDLE, FOUND         )
   C
   C        Due to the way we've constructed the loop, there's
   C        no need to check the FOUND flag.
   C
            CALL TOSTDO ( ' '   )
 
            LINE = '   Kernel name:    #'
            CALL REPMC  ( LINE, '#', FILE, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '   Kernel type:    #'
            CALL REPMC  ( LINE, '#', FILTYP, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '   Kernel source:  #'
            CALL REPMC  ( LINE, '#', SOURCE, LINE )
            CALL TOSTDO ( LINE )
 
         END DO
 
         END


Solution Sample Output



Numerical results shown for this example may differ across platforms since the results depend on the SPICE kernels used as input and on the host platform's arithmetic implementation.

After compiling the program, execute it. The output is:

 
   Inputs for geometric event finding program
 
      Target                       = MEX
      Observation surface location = DSS-14
      Occulting body               = MARS
      Aberration correction        = CN+S
      Start time                   = 2004 APR 30 23:58:55.814 (UTC)
      Start time                   = 2004 MAY 01 00:00:00.000 (TDB)
      Stop time                    = 2004 MAY 04 23:58:55.814 (UTC)
      Stop time                    = 2004 MAY 05 00:00:00.000 (TDB)
      Elevation limit (degrees)    = 6.000000
 
   Loaded SPICE Kernels:
 
      Kernel name:    geomevnt.mk
      Kernel type:    META
      Kernel source:
 
      Kernel name:    de405s.bsp
      Kernel type:    SPK
      Kernel source:  geomevnt.mk
 
      Kernel name:    dsnstns.bsp
      Kernel type:    SPK
      Kernel source:  geomevnt.mk
 
      Kernel name:    DSN_topo.frm
      Kernel type:    TEXT
      Kernel source:  geomevnt.mk
 
      Kernel name:    earth_031228_231229_predict.bpc
      Kernel type:    PCK
      Kernel source:  geomevnt.mk
 
      Kernel name:    naif0007.tls
      Kernel type:    TEXT
      Kernel source:  geomevnt.mk
 
      Kernel name:    ORMM__040501000000_00069.BSP
      Kernel type:    SPK
      Kernel source:  geomevnt.mk
 
      Kernel name:    pck00007.tpc
      Kernel type:    TEXT
      Kernel source:  geomevnt.mk


Find Azimuth and Elevation of Target







Task Statement




Write a subroutine to find the azimuth and elevation of a target as seen from a specified location on the surface of an extended body. Also return the range to the target.

The output angles are measured relative to a topocentric frame centered at the surface location. Angular units are radians.

In the topocentric frame, elevation is equivalent to latitude, and azimuth is the negative of longitude.

The subroutine is to have the exact signature shown below:

   SUBROUTINE SRFAZL ( SRFPT, ET, ABCORR, TARG, AZ, EL, R )
 
   CHARACTER*(*)         SRFPT
   DOUBLE PRECISION      ET
   CHARACTER*(*)         ABCORR
   CHARACTER*(*)         TARG
   DOUBLE PRECISION      AZ
   DOUBLE PRECISION      EL
   DOUBLE PRECISION      R
The meanings of the arguments are as follows:

   SRFPT      Name of surface point, e.g. 'DSS-14'
   ET         Computation epoch, seconds past J2000 TDB.
   ABCORR     Aberration correction flag: 'NONE', 'LT+S', etc.
   TARG       Target name.
   AZ         Azimuth in radians.
   EL         Elevation in radians.
   R          Range in kilometers.
This subroutine will be used by the supplied event finding utilities to find times when the target is above the specified elevation limit.

Test the routine by calling it twice and displaying the outputs. The first call should use the input values:

   SRFPT      'DSS-14'
   ET         Double precision count of seconds past J2000 TDB
              corresponding to 2004 MAY 01 16:10:26.686 TDB
   ABCORR     'CN+S'
   TARG       'MEX'
The second call should use the same inputs as the first, except that the time should be the double precision count of seconds past J2000 TDB corresponding to 2004 MAY 01 16:10:26.688 TDB.

The two times are intended to bracket a spacecraft ``rise'' event: the target spacecraft should have elevation below the specified limit at the first time and above the limit at the second.



Learning Goals




This exercise introduces the student to elementary geometry computations using SPICELIB. The student uses SPKPOS to look up position vectors and RECLAT to convert from rectangular to latitudinal coordinates.

When the SPICELIB event finding software is delivered, computations of this type usually will not have to be performed by users in the course of solving event-finding problems: the SPICELIB event finding software will do the job behind the scenes. However, this exercise is meant to provide the student with some insight into the geometry computations involved in an event finding system.



Approach






Solution steps



A possible solution would consist of the following steps:

Write the SRFAZL routine:

Add test code to the main program:



Solution






SRFAZL Solution Code



The example program below shows one possible implementation of the routine SRFAZL.

 
         SUBROUTINE SRFAZL ( SRFPT, ET, ABCORR, TARG, AZ, EL, R )
         IMPLICIT NONE
   C
   C     Find the azimuth and elevation of a target as seen from
   C     a specified location on the surface of an extended body.
   C     Also return the range to the target.
   C
   C     The output angles are measured relative to a topocentric
   C     frame centered at the surface location.  Angular units are
   C     radians.
   C
   C
   C     Subroutine arguments:
   C
   C
   C        SRFPT         is a string giving the name of the surface
   C                      point of interest.  A topocentric frame is
   C                      defined at this point.
   C
   C        ET            is the epoch at which the location of the
   C                      target is to be found.  ET is given in
   C                      seconds past J2000 (TDB).
   C
   C        ABCORR        indicates the aberration correction.
   C                      Values and meanings are:
   C
   C                         'NONE'   No correction.  Use geometric
   C                                  states.
   C
   C                        The following are for the "reception"
   C                        case where radiation travels from the
   C                        target to the surface point and arrives
   C                        at ET:
   C
   C                         'LT'     Correct for one-way light time.
   C
   C                         'LT+S'   Correct for one-way light time
   C                                  and stellar aberration.
   C
   C                         'CN'     Converged Newtonian light time
   C                                  correction.
   C
   C                         'CN+S'   Converged Newtonian light time
   C                                  and stellar aberration
   C                                  corrections.
   C
   C
   C                        The following are for the "transmission"
   C                        case where radiation travels from the
   C                        surface point to the target, departing the
   C                        surface location at ET:
   C
   C                         'XLT'    Correct for one-way light time.
   C
   C                         'XLT+S'  Correct for one-way light time
   C                                  and stellar aberration.
   C
   C                         'XCN'    Converged Newtonian light time
   C                                  correction.
   C
   C                         'XCN+S'  Converged Newtonian light time
   C                                  and stellar aberration
   C                                  corrections.
   C
   C
   C        TARG          is a string containing the target name.
   C
   C        AZ,
   C        EL            are the azimuth and elevation of the target
   C                      vector at ET.  Azimuth and elevation are
   C                      given with respect to a right-handed frame
   C                      with z-axis equal to the local surface
   C                      normal, x-axis orthogonal to the z-axis and
   C                      pointing north, and the y-axis equal to
   C                      z cross x.
   C
   C                      Azimuth is the angle between the x-axis and
   C                      the projection of the target vector onto the
   C                      topocentric x-y plane.  Azimuth increases in
   C                      the clockwise direction.
   C
   C                      Elevation is the angle between the target
   C                      vector and the x-y plane.
   C
   C                      Both angles are given in radians.
   C
   C        R             is the range to the target, in km.
   C
   C
   C- Version
   C
   C   Geometric Event Finding Code Lesson 1.0.0 07-MAY-2004 (NJB)
   C
 
         CHARACTER*(*)         SRFPT
         DOUBLE PRECISION      ET
         CHARACTER*(*)         ABCORR
         CHARACTER*(*)         TARG
 
         DOUBLE PRECISION      AZ
         DOUBLE PRECISION      EL
         DOUBLE PRECISION      R
 
   C
   C     SPICELIB functions
   C
         LOGICAL               RETURN
 
   C
   C     Local parameters
   C
         INTEGER               FRNMLN
         PARAMETER           ( FRNMLN = 32 )
 
   C
   C     Local variables
   C
         CHARACTER*(FRNMLN)    FRNAME
 
         DOUBLE PRECISION      LON
         DOUBLE PRECISION      LT
         DOUBLE PRECISION      TRGPOS ( 3 )
 
         INTEGER               FRCODE
 
         LOGICAL               FOUND
 
 
         IF ( RETURN() ) THEN
            RETURN
         END IF
 
         CALL CHKIN ( 'SRFAZL' )
 
 
   C
   C     Look up the topocentric reference frame associated with the
   C     surface point.  The definition of this frame is supplied in
   C     a frame kernel loaded by the calling program.
   C
         CALL CNMFRM ( SRFPT, FRCODE, FRNAME, FOUND )
 
         IF ( .NOT. FOUND ) THEN
 
            CALL SETMSG ( 'No frame is associated with body #. ' //
        .                 'This can be due to a required frame ' //
        .                 'kernel not having been loaded.'       )
            CALL ERRCH  ( '#',  SRFPT                            )
            CALL SIGERR ( 'SPICE(NOTDEFINED)'                    )
            CALL CHKOUT ( 'SRFAZL'                               )
            RETURN
 
         END IF
 
   C
   C     Find the position of the target as seen from the observer
   C     at ET.  Use the specified aberration corrections.  Ask for
   C     the position vector in the topocentric frame associated
   C     with the surface point.
   C
         CALL SPKPOS ( TARG, ET, FRNAME, ABCORR, SRFPT, TRGPOS, LT )
 
   C
   C     Elevation is the angular separation of the surface
   C     point-target vector from the x-y plane of the
   C     topocentric frame; this is simply latitude.  Azimuth is the
   C     negative of the longitude angle.
   C
         CALL RECLAT ( TRGPOS, R, LON, EL )
 
         AZ  =  - LON
 
 
         CALL CHKOUT ( 'SRFAZL' )
         END


Solution: Add Test Code to the Main Program



The example program below shows one possible solution.

         PROGRAM SCVIEW
         IMPLICIT NONE
   C
   C     Load SPICE Kernels and read additional inputs required to
   C     solve the spacecraft visibility problem posed by the
   C     geometric event finding hands-on code lesson.
   C
   C     The following items are loaded from a SPICE text kernel:
   C
   C        - Observer (surface point) name
   C        - Target name
   C        - Occulting body name
   C        - Aberration correction specification
   C        - Start time string
   C        - Stop time string
   C        - Elevation limit (in degrees)
   C
   C     This program assumes all necessary kernels and kernel
   C     variable definitions are specified in a SPICE meta-kernel.
   C     The meta-kernel name is assumed to be 'geomevnt.mk'.
   C
 
   C
   C     SPICELIB functions
   C
         DOUBLE PRECISION      DPR
         DOUBLE PRECISION      RPD
         LOGICAL               BADKPV
 
   C
   C     Other functions
   C
 
   C
   C     Local parameters
   C
 
   C
   C     The meta-kernel:
   C
         CHARACTER*(*)         META
         PARAMETER           ( META   = 'geomevnt.mk' )
 
   C
   C     Labels of kernel pool items:
   C
         CHARACTER*(*)         OCCLBL
         PARAMETER           ( OCCLBL = 'SCVIEW_OCCBDY' )
 
         CHARACTER*(*)         SRFLBL
         PARAMETER           ( SRFLBL = 'SCVIEW_SRFPT' )
 
         CHARACTER*(*)         TRGLBL
         PARAMETER           ( TRGLBL = 'SCVIEW_TARGET' )
 
         CHARACTER*(*)         CORLBL
         PARAMETER           ( CORLBL = 'SCVIEW_ABCORR' )
 
         CHARACTER*(*)         STRLBL
         PARAMETER           ( STRLBL = 'SCVIEW_START' )
 
         CHARACTER*(*)         STPLBL
         PARAMETER           ( STPLBL = 'SCVIEW_STOP' )
 
         CHARACTER*(*)         ELVLBL
         PARAMETER           ( ELVLBL = 'SCVIEW_ELVLIM' )
 
   C
   C     String length parameters:
   C
         INTEGER               CORLEN
         PARAMETER           ( CORLEN = 10 )
 
         INTEGER               FILSIZ
         PARAMETER           ( FILSIZ = 255 )
 
         INTEGER               LNSIZE
         PARAMETER           ( LNSIZE = 200 )
 
         INTEGER               NAMLEN
         PARAMETER           ( NAMLEN = 32 )
 
         INTEGER               TIMLEN
         PARAMETER           ( TIMLEN = 50 )
 
         INTEGER               TYPLEN
         PARAMETER           ( TYPLEN = 50 )
 
   C
   C     Local variables
   C
         CHARACTER*(CORLEN)    ABCORR
         CHARACTER*(FILSIZ)    FILE
         CHARACTER*(TYPLEN)    FILTYP
         CHARACTER*(LNSIZE)    LINE
         CHARACTER*(NAMLEN)    OCCBDY
         CHARACTER*(FILSIZ)    SOURCE
         CHARACTER*(NAMLEN)    SRFPT
         CHARACTER*(TIMLEN)    START
         CHARACTER*(TIMLEN)    STOP
         CHARACTER*(NAMLEN)    TARGET
         CHARACTER*(TIMLEN)    TIMSTR
         CHARACTER*(LNSIZE)    TITLE
         CHARACTER*(TIMLEN)    TSTTIM ( 2 )
 
         DOUBLE PRECISION      AZ
         DOUBLE PRECISION      EL
         DOUBLE PRECISION      ELVLIM
         DOUBLE PRECISION      ET
         DOUBLE PRECISION      ETBEG
         DOUBLE PRECISION      ETEND
         DOUBLE PRECISION      R
         DOUBLE PRECISION      REVLIM
 
         INTEGER               HANDLE
         INTEGER               I
         INTEGER               N
         INTEGER               NKER
 
         LOGICAL               BAD
         LOGICAL               FOUND
 
   C
   C     Load the meta-kernel.
   C
         CALL FURNSH ( META )
 
   C
   C     Look up observation location, target,  name of the
   C     occulting body, aberration correction, start
   C     and stop times, and elevation limit.
   C
   C     We simplify error checking by using BADKPV.  BADKPV
   C     tests whether a kernel variable having specified
   C     attributes is present in the kernel pool. The fourth
   C     argument of BADKPV is the expected dimension; the
   C     last indicates the expected data type.  Read the
   C     header of BADKPV for details.
   C
         BAD = BADKPV ( 'SCVIEW', SRFLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', TRGLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', OCCLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', CORLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', STRLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', STPLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', ELVLBL, '=',  1,      1, 'N' )
 
   C
   C     Now we know the kernel variables of interest have been
   C     defined; look up the values.  We don't need to check
   C     the FOUND flag because BADKPV has ensured that each item
   C     will be found.
   C
   C     Note that the elevation limit has units of degrees on input,
   C     so we store the input in the variable ELVLIM.  We convert
   C     this value to radians for computation using SPICE routines;
   C     we'll store the equivalent value in radians in REVLIM.
   C
         CALL GCPOOL ( SRFLBL, 1, 1,      N,    SRFPT,  FOUND )
         CALL GCPOOL ( TRGLBL, 1, 1,      N,    TARGET, FOUND )
         CALL GCPOOL ( OCCLBL, 1, 1,      N,    OCCBDY, FOUND )
         CALL GCPOOL ( CORLBL, 1, 1,      N,    ABCORR, FOUND )
         CALL GCPOOL ( STRLBL, 1, 1,      N,    START,  FOUND )
         CALL GCPOOL ( STPLBL, 1, 1,      N,    STOP,   FOUND )
         CALL GDPOOL ( ELVLBL, 1, 1,      N,    ELVLIM, FOUND )
 
         REVLIM = RPD() * ELVLIM
 
   C
   C     Display to standard output a banner for the output report:
   C
         TITLE  = 'Inputs for geometric event finding program'
         CALL TOSTDO ( ' '   )
         CALL TOSTDO ( TITLE )
         CALL TOSTDO ( ' '   )
 
         LINE   = '   Target                       = #'
         CALL REPMC  ( LINE, '#', TARGET, LINE )
         CALL TOSTDO ( LINE )
 
         LINE   = '   Observation surface location = #'
         CALL REPMC  ( LINE, '#', SRFPT, LINE )
         CALL TOSTDO ( LINE )
 
         LINE   = '   Aberration correction        = #'
         CALL REPMC  ( LINE, '#', ABCORR, LINE )
         CALL TOSTDO ( LINE )
 
   C
   C     Convert the start and stop times to ET.
   C
         CALL STR2ET ( START, ETBEG )
         CALL STR2ET ( STOP,  ETEND )
 
   C
   C     Display the start time as both calendar UTC and
   C     calendar ET using the formats shown below.
   C
   C        2004 MAY 06 20:15:00.000 (UTC)
   C        2004 MAY 06 20:15:00.000 (TDB)
   C
         CALL TIMOUT ( ETBEG,
        .              'YYYY MON DD HR:MN:SC.### (UTC)',
        .              TIMSTR                           )
 
         LINE   = '   Start time                   = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
 
         CALL TIMOUT ( ETBEG,
        .              'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .              TIMSTR                           )
 
         LINE   = '   Start time                   = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
 
         CALL TIMOUT ( ETEND,
        .              'YYYY MON DD HR:MN:SC.### (UTC)',
        .              TIMSTR                           )
 
         LINE   = '   Stop time                    = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
 
         CALL TIMOUT ( ETEND,
        .              'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .              TIMSTR                                   )
 
         LINE   = '   Stop time                    = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
   C
   C     Display the elevation limit in degrees.
   C
         LINE   = '   Elevation limit (degrees)    = #'
         CALL REPMF  ( LINE, '#', ELVLIM, 7, 'F', LINE )
         CALL TOSTDO ( LINE )
 
   C
   C     Display the names of the SPICE kernels we've loaded.
   C
         CALL TOSTDO ( ' '   )
         CALL TOSTDO ( 'Loaded SPICE Kernels: ' )
 
         CALL KTOTAL ( 'ALL', NKER )
 
         DO I = 1, NKER
 
            CALL KDATA ( I,      'ALL',  FILE, FILTYP,
        .                SOURCE, HANDLE, FOUND         )
   C
   C        Due to the way we've constructed the loop, there's
   C        no need to check the FOUND flag.
   C
            CALL TOSTDO ( ' '   )
 
            LINE = '   Kernel name:    #'
            CALL REPMC  ( LINE, '#', FILE, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '   Kernel type:    #'
            CALL REPMC  ( LINE, '#', FILTYP, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '   Kernel source:  #'
            CALL REPMC  ( LINE, '#', SOURCE, LINE )
            CALL TOSTDO ( LINE )
 
         END DO
 
   C
   C     Test SRFAZL:  make two calls at times bracketing
   C     a spacecraft "rise" event.
   C
         CALL TOSTDO ( ' ' )
         CALL TOSTDO ( 'SRFAZL test results:' )
         CALL TOSTDO ( ' ' )
 
         TSTTIM(1) = '2004 MAY 01 16:10:26.686 TDB'
         TSTTIM(2) = '2004 MAY 01 16:10:26.688 TDB'
 
         DO I = 1, 2
   C
   C        Convert the TDB calendar time to seconds past J2000.
   C        Find the azimuth, elevation, and range at this time.
   C
            CALL STR2ET ( TSTTIM(I), ET )
            CALL SRFAZL ( SRFPT,     ET, ABCORR, TARGET, AZ, EL, R )
   C
   C        Display the time and SRFAZL outputs.  Convert radians to
   C        degrees for output using the SPICELIB function DPR.
   C
            LINE = '   ET = #'
            CALL REPMC  ( LINE, '#', TSTTIM(I),   LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '      Elevation (degrees):  #'
            CALL REPMD  ( LINE, '#', DPR()*EL,  16, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '      Azimuth   (degrees):  #'
            CALL REPMD  ( LINE, '#', DPR()*AZ,  16, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '      Range     (km):       #'
            CALL REPMD  ( LINE, '#', R,         16, LINE )
            CALL TOSTDO ( LINE )
            CALL TOSTDO ( ' '  )
 
         END DO
 
         END


Solution Sample Output



Numerical results shown for this example may differ across platforms since the results depend on the SPICE kernels used as input and on the host platform's arithmetic implementation.

 
   Inputs for geometric event finding program
 
      Target                       = MEX
      Observation surface location = DSS-14
      Aberration correction        = CN+S
      Start time                   = 2004 APR 30 23:58:55.814 (UTC)
      Start time                   = 2004 MAY 01 00:00:00.000 (TDB)
      Stop time                    = 2004 MAY 04 23:58:55.814 (UTC)
      Stop time                    = 2004 MAY 05 00:00:00.000 (TDB)
      Elevation limit (degrees)    = 6.000000
 
   Loaded SPICE Kernels:
 
      Kernel name:    geomevnt.mk
      Kernel type:    META
      Kernel source:
 
      Kernel name:    de405s.bsp
      Kernel type:    SPK
      Kernel source:  geomevnt.mk
 
      Kernel name:    dsnstns.bsp
      Kernel type:    SPK
      Kernel source:  geomevnt.mk
 
      Kernel name:    DSN_topo.frm
      Kernel type:    TEXT
      Kernel source:  geomevnt.mk
 
      Kernel name:    earth_031228_231229_predict.bpc
      Kernel type:    PCK
      Kernel source:  geomevnt.mk
 
      Kernel name:    naif0007.tls
      Kernel type:    TEXT
      Kernel source:  geomevnt.mk
 
      Kernel name:    ORMM__040501000000_00069.BSP
      Kernel type:    SPK
      Kernel source:  geomevnt.mk
 
      Kernel name:    pck00007.tpc
      Kernel type:    TEXT
      Kernel source:  geomevnt.mk
 
   SRFAZL test results:
 
      ET = 2004 MAY 01 16:10:26.686 TDB
         Elevation (degrees):  5.9999951213034E+00
         Azimuth   (degrees):  6.3942106540558E+01
         Range     (km):       3.2480333721584E+08
 
      ET = 2004 MAY 01 16:10:26.688 TDB
         Elevation (degrees):  6.0000012274140E+00
         Azimuth   (degrees):  6.3942111060397E+01
         Range     (km):       3.2480333724230E+08
 


Test for Occultation of Point Target







Task Statement




Write a LOGICAL function to test when a point target is occulted by a specified extended body, as seen from a specified observing body.

The function is to have the exact signature shown below:

   LOGICAL FUNCTION OCCPT ( TARGET, OCCBDY, ET, ABCORR, OBSRVR )
The meanings of the arguments are as follows:

   TARGET     Target name.
   OCCBDY     Occulting body name.
   ET         Computation epoch, seconds past J2000 TDB.
   ABCORR     Aberration correction flag: 'NONE', 'LT+S', etc.
   OBSRVR     Observing body name.
   OCCPT      LOGICAL function value:  .TRUE. if the target
              is occulted; .FALSE. otherwise.
This subroutine will be used by the supplied event finding utilities to find times when the target is occulted by the specified occulting body.

Test the routine by calling it twice and displaying the outputs. The first call should use the input values:

   TARGET     'MEX'
   OCCBDY     'MARS'
   ET         Double precision count of seconds past J2000 TDB
              corresponding to 2004 MAY 01 13:38:36.740 TDB
   ABCORR     'CN+S'
   OBSRVR     'DSS-14'
The second call should use the same inputs as the first, except that the time should be the double precision count of seconds past J2000 TDB corresponding to 2004 MAY 01 13:38:36.742 TDB.

The two times are intended to bracket a spacecraft occultation ingress event: the aberration-corrected target spacecraft location should be visible at the first time and occulted at the second.



Learning Goals




This exercise introduces the student to further elementary geometry computations using SPICELIB. The student uses SPKPOS to look up position vectors and SRFXPT to test whether a specified ray intersects an extended body modeled as a triaxial ellipsoid.

This exercise is meant to provide the student with further understanding of the geometry computations involved in an event finding system.



Approach






Solution steps



A possible solution would consist of the following steps:

Write the OCCPT routine:

Add test code to the main program:



Solution






Solution Code



The example program below shows one possible solution.

 
         LOGICAL FUNCTION OCCPT ( TARGET, OCCBDY, ET, ABCORR, OBSRVR )
         IMPLICIT NONE
   C
   C     Determine whether a point target is occulted by extended body,
   C     as seen from a specified observation location.
   C
   C
   C     Subroutine arguments:
   C
   C
   C        TARGET        is a string giving the name of the target
   C                      body of interest
   C
   C        OCCBDY        is a string giving the name of a potentially
   C                      occulting extended body.
   C
   C        ET            is the epoch at which the visibility
   C                      computation is to be performed.  ET is
   C                      given in seconds past J2000 (TDB).
   C
   C        ABCORR        indicates the aberration correction.
   C                      Values and meanings are:
   C
   C                         'NONE'   No correction.  Use geometric
   C                                  states.
   C
   C                        The following are for the "reception"
   C                        case where radiation travels from the
   C                        target body and arrives at the observer
   C                        at ET:
   C
   C                         'LT'     Correct for one-way light time.
   C
   C                         'LT+S'   Correct for one-way light time
   C                                  and stellar aberration.
   C
   C                         'CN'     Converged Newtonian light time
   C                                  correction.
   C
   C                         'CN+S'   Converged Newtonian light time
   C                                  and stellar aberration
   C                                  corrections.
   C
   C
   C                        The following are for the "transmission"
   C                        case where radiation travels from the
   C                        observer to the target body, departing the
   C                        observer's location at ET:
   C
   C                         'XLT'    Correct for one-way light time.
   C
   C                         'XLT+S'  Correct for one-way light time
   C                                  and stellar aberration.
   C
   C                         'XCN'    Converged Newtonian light time
   C                                  correction.
   C
   C                         'XCN+S'  Converged Newtonian light time
   C                                  and stellar aberration
   C                                  corrections.
   C
   C        OBSRVR        is a string containing the observer's name.
   C
   C
   C     The function returns .TRUE. if the target is occulted and
   C     .FALSE. otherwise.
   C
   C
 
         CHARACTER*(*)         TARGET
         CHARACTER*(*)         OCCBDY
         DOUBLE PRECISION      ET
         CHARACTER*(*)         ABCORR
         CHARACTER*(*)         OBSRVR
 
   C
   C     SPICELIB functions
   C
         DOUBLE PRECISION      VNORM
 
         LOGICAL               RETURN
 
   C
   C     Local variables
   C
         DOUBLE PRECISION      DIST
         DOUBLE PRECISION      OBSPOS ( 3 )
         DOUBLE PRECISION      TRGEPC
         DOUBLE PRECISION      TRGLT
         DOUBLE PRECISION      TRGPOS ( 3 )
         DOUBLE PRECISION      XPT    ( 3 )
 
         LOGICAL               FOUND
 
         IF ( RETURN() ) THEN
            OCCPT = .FALSE.
            RETURN
         END IF
 
         CALL CHKIN ( 'OCCPT' )
 
   C
   C     Use SRFXPT to find the intersection of the
   C     observer-target vector with the occulting
   C     body, if any.  Before calling SRFXPT, we must
   C     find the observer-target vector, using the
   C     specified aberration corrections.  We'll look up
   C     the vector in the J2000 reference frame.
   C
         CALL SPKPOS ( TARGET, ET,    'J2000', ABCORR,
        .              OBSRVR, TRGPOS, TRGLT          )
 
         CALL SRFXPT ( 'Ellipsoid', OCCBDY,    ET,
        .              ABCORR,      OBSRVR,    'J2000',
        .              TRGPOS,      XPT,       DIST,
        .              TRGEPC,      OBSPOS,    FOUND    )
 
         IF ( .NOT. FOUND ) THEN
   C
   C        If there's no intercept, the target is not occulted.
   C        We're done.
   C
            OCCPT = .FALSE.
 
            CALL CHKOUT ( 'OCCPT' )
            RETURN
 
         END IF
 
   C
   C     At this point, we know there's an intercept.  If the target
   C     lies between the observer and the intercept point, the target
   C     is still visible.  Otherwise, the target is occulted.
   C
         OCCPT  =  VNORM( TRGPOS ) .GT. DIST
 
         CALL CHKOUT ( 'OCCPT' )
         END


Solution: Add Test Code to the Main Program



The example program below shows one possible solution.

         PROGRAM SCVIEW
         IMPLICIT NONE
   C
   C     Load SPICE Kernels and read additional inputs required to
   C     solve the spacecraft visibility problem posed by the
   C     geometric event finding hands-on code lesson.
   C
   C     The following items are loaded from a SPICE text kernel:
   C
   C        - Observer (surface point) name
   C        - Target name
   C        - Occulting body name
   C        - Aberration correction specification
   C        - Start time string
   C        - Stop time string
   C        - Elevation limit (in degrees)
   C
   C     This program assumes all necessary kernels and kernel
   C     variable definitions are specified in a SPICE meta-kernel.
   C     The meta-kernel name is assumed to be 'geomevnt.mk'.
   C
 
   C
   C     SPICELIB functions
   C
         DOUBLE PRECISION      DPR
         DOUBLE PRECISION      RPD
         LOGICAL               BADKPV
 
   C
   C     Other functions
   C
         LOGICAL               OCCPT
 
   C
   C     Local parameters
   C
 
   C
   C     The meta-kernel:
   C
         CHARACTER*(*)         META
         PARAMETER           ( META   = 'geomevnt.mk' )
 
   C
   C     Labels of kernel pool items:
   C
         CHARACTER*(*)         OCCLBL
         PARAMETER           ( OCCLBL = 'SCVIEW_OCCBDY' )
 
         CHARACTER*(*)         SRFLBL
         PARAMETER           ( SRFLBL = 'SCVIEW_SRFPT' )
 
         CHARACTER*(*)         TRGLBL
         PARAMETER           ( TRGLBL = 'SCVIEW_TARGET' )
 
         CHARACTER*(*)         CORLBL
         PARAMETER           ( CORLBL = 'SCVIEW_ABCORR' )
 
         CHARACTER*(*)         STRLBL
         PARAMETER           ( STRLBL = 'SCVIEW_START' )
 
         CHARACTER*(*)         STPLBL
         PARAMETER           ( STPLBL = 'SCVIEW_STOP' )
 
         CHARACTER*(*)         ELVLBL
         PARAMETER           ( ELVLBL = 'SCVIEW_ELVLIM' )
 
   C
   C     String length parameters:
   C
         INTEGER               CORLEN
         PARAMETER           ( CORLEN = 10 )
 
         INTEGER               FILSIZ
         PARAMETER           ( FILSIZ = 255 )
 
         INTEGER               LNSIZE
         PARAMETER           ( LNSIZE = 200 )
 
         INTEGER               NAMLEN
         PARAMETER           ( NAMLEN = 32 )
 
         INTEGER               TIMLEN
         PARAMETER           ( TIMLEN = 50 )
 
         INTEGER               TYPLEN
         PARAMETER           ( TYPLEN = 50 )
 
   C
   C     Local variables
   C
         CHARACTER*(CORLEN)    ABCORR
         CHARACTER*(FILSIZ)    FILE
         CHARACTER*(TYPLEN)    FILTYP
         CHARACTER*(LNSIZE)    LINE
         CHARACTER*(NAMLEN)    OCCBDY
         CHARACTER*(FILSIZ)    SOURCE
         CHARACTER*(NAMLEN)    SRFPT
         CHARACTER*(TIMLEN)    START
         CHARACTER*(TIMLEN)    STOP
         CHARACTER*(NAMLEN)    TARGET
         CHARACTER*(TIMLEN)    TIMSTR
         CHARACTER*(LNSIZE)    TITLE
         CHARACTER*(TIMLEN)    TSTTIM ( 2 )
 
         DOUBLE PRECISION      AZ
         DOUBLE PRECISION      EL
         DOUBLE PRECISION      ELVLIM
         DOUBLE PRECISION      ET
         DOUBLE PRECISION      ETBEG
         DOUBLE PRECISION      ETEND
         DOUBLE PRECISION      R
         DOUBLE PRECISION      REVLIM
 
         INTEGER               HANDLE
         INTEGER               I
         INTEGER               N
         INTEGER               NKER
 
         LOGICAL               BAD
         LOGICAL               FOUND
         LOGICAL               HIDDEN
 
   C
   C     Load the meta-kernel.
   C
         CALL FURNSH ( META )
 
   C
   C     Look up observation location, target,  name of the
   C     occulting body, aberration correction, start
   C     and stop times, and elevation limit.
   C
   C     We simplify error checking by using BADKPV.  BADKPV
   C     tests whether a kernel variable having specified
   C     attributes is present in the kernel pool. The fourth
   C     argument of BADKPV is the expected dimension; the
   C     last indicates the expected data type.  Read the
   C     header of BADKPV for details.
   C
         BAD = BADKPV ( 'SCVIEW', SRFLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', TRGLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', OCCLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', CORLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', STRLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', STPLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', ELVLBL, '=',  1,      1, 'N' )
 
   C
   C     Now we know the kernel variables of interest have been
   C     defined; look up the values.  We don't need to check
   C     the FOUND flag because BADKPV has ensured that each item
   C     will be found.
   C
   C     Note that the elevation limit has units of degrees on input,
   C     so we store the input in the variable ELVLIM.  We convert
   C     this value to radians for computation using SPICE routines;
   C     we'll store the equivalent value in radians in REVLIM.
   C
         CALL GCPOOL ( SRFLBL, 1, 1,      N,    SRFPT,  FOUND )
         CALL GCPOOL ( TRGLBL, 1, 1,      N,    TARGET, FOUND )
         CALL GCPOOL ( OCCLBL, 1, 1,      N,    OCCBDY, FOUND )
         CALL GCPOOL ( CORLBL, 1, 1,      N,    ABCORR, FOUND )
         CALL GCPOOL ( STRLBL, 1, 1,      N,    START,  FOUND )
         CALL GCPOOL ( STPLBL, 1, 1,      N,    STOP,   FOUND )
         CALL GDPOOL ( ELVLBL, 1, 1,      N,    ELVLIM, FOUND )
 
         REVLIM = RPD() * ELVLIM
 
   C
   C     Display to standard output a banner for the output report:
   C
         TITLE  = 'Inputs for geometric event finding program'
         CALL TOSTDO ( ' '   )
         CALL TOSTDO ( TITLE )
         CALL TOSTDO ( ' '   )
 
         LINE   = '   Target                       = #'
         CALL REPMC  ( LINE, '#', TARGET, LINE )
         CALL TOSTDO ( LINE )
 
         LINE   = '   Observation surface location = #'
         CALL REPMC  ( LINE, '#', SRFPT, LINE )
         CALL TOSTDO ( LINE )
 
         LINE   = '   Aberration correction        = #'
         CALL REPMC  ( LINE, '#', ABCORR, LINE )
         CALL TOSTDO ( LINE )
 
   C
   C     Convert the start and stop times to ET.
   C
         CALL STR2ET ( START, ETBEG )
         CALL STR2ET ( STOP,  ETEND )
 
   C
   C     Display the start time as both calendar UTC and
   C     calendar ET using the formats shown below.
   C
   C        2004 MAY 06 20:15:00.000 (UTC)
   C        2004 MAY 06 20:15:00.000 (TDB)
   C
         CALL TIMOUT ( ETBEG,
        .              'YYYY MON DD HR:MN:SC.### (UTC)',
        .              TIMSTR                           )
 
         LINE   = '   Start time                   = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
 
         CALL TIMOUT ( ETBEG,
        .              'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .              TIMSTR                           )
 
         LINE   = '   Start time                   = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
 
         CALL TIMOUT ( ETEND,
        .              'YYYY MON DD HR:MN:SC.### (UTC)',
        .              TIMSTR                           )
 
         LINE   = '   Stop time                    = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
 
         CALL TIMOUT ( ETEND,
        .              'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .              TIMSTR                                   )
 
         LINE   = '   Stop time                    = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
   C
   C     Display the elevation limit in degrees.
   C
         LINE   = '   Elevation limit (degrees)    = #'
         CALL REPMF  ( LINE, '#', ELVLIM, 7, 'F', LINE )
         CALL TOSTDO ( LINE )
 
   C
   C     Display the names of the SPICE kernels we've loaded.
   C
         CALL TOSTDO ( ' '   )
         CALL TOSTDO ( 'Loaded SPICE Kernels: ' )
 
         CALL KTOTAL ( 'ALL', NKER )
 
         DO I = 1, NKER
 
            CALL KDATA ( I,      'ALL',  FILE, FILTYP,
        .                SOURCE, HANDLE, FOUND         )
   C
   C        Due to the way we've constructed the loop, there's
   C        no need to check the FOUND flag.
   C
            CALL TOSTDO ( ' '   )
 
            LINE = '   Kernel name:    #'
            CALL REPMC  ( LINE, '#', FILE, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '   Kernel type:    #'
            CALL REPMC  ( LINE, '#', FILTYP, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '   Kernel source:  #'
            CALL REPMC  ( LINE, '#', SOURCE, LINE )
            CALL TOSTDO ( LINE )
 
         END DO
 
   C
   C     Test SRFAZL:  make two calls at times bracketing
   C     a spacecraft "rise" event.
   C
         CALL TOSTDO ( ' ' )
         CALL TOSTDO ( 'SRFAZL test results:' )
         CALL TOSTDO ( ' ' )
 
         TSTTIM(1) = '2004 MAY 01 16:10:26.686 TDB'
         TSTTIM(2) = '2004 MAY 01 16:10:26.688 TDB'
 
         DO I = 1, 2
   C
   C        Convert the TDB calendar time to seconds past J2000.
   C        Find the azimuth, elevation, and range at this time.
   C
            CALL STR2ET ( TSTTIM(I), ET )
            CALL SRFAZL ( SRFPT,     ET, ABCORR, TARGET, AZ, EL, R )
   C
   C        Display the time and elevation.  Convert radians to
   C        degrees for output using the SPICELIB function DPR.
   C
            LINE = '   ET = #'
            CALL REPMC  ( LINE, '#', TSTTIM(I),   LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '      Elevation (degrees):  #'
            CALL REPMD  ( LINE, '#', DPR()*EL,  16, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '      Azimuth   (degrees):  #'
            CALL REPMD  ( LINE, '#', DPR()*AZ,  16, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '      Range     (km):       #'
            CALL REPMD  ( LINE, '#', R,         16, LINE )
            CALL TOSTDO ( LINE )
            CALL TOSTDO ( ' '  )
 
         END DO
 
 
   C
   C     Test OCCPT:  make two calls at times bracketing
   C     a spacecraft occultation ingress event.
   C
         CALL TOSTDO ( 'OCCPT test results:' )
         CALL TOSTDO ( ' ' )
 
         TSTTIM(1) =  '2004 MAY 01 13:38:36.740 TDB'
         TSTTIM(2) =  '2004 MAY 01 13:38:36.742 TDB'
 
         DO I = 1, 2
   C
   C        Convert the TDB calendar time to seconds past J2000.
   C
            CALL STR2ET ( TSTTIM(I), ET )
            HIDDEN = OCCPT ( TARGET, OCCBDY, ET, ABCORR, SRFPT )
   C
   C        Display the time and occultation state.
   C
            LINE = '   ET = #.  #.'
            CALL REPMC  ( LINE, '#', TSTTIM(I),   LINE )
 
            IF ( HIDDEN ) THEN
               CALL REPMC ( LINE, '#', 'Occulted', LINE )
            ELSE
               CALL REPMC ( LINE, '#', 'Not occulted', LINE )
            END IF
 
            CALL TOSTDO ( LINE )
 
         END DO
 
         CALL TOSTDO ( ' ' )
 
         END


Solution Sample Output



Numerical results shown for this example may differ across platforms since the results depend on the SPICE kernels used as input and on the host platform's arithmetic implementation.

 
   Inputs for geometric event finding program
 
      Target                       = MEX
      Observation surface location = DSS-14
      Aberration correction        = CN+S
      Start time                   = 2004 APR 30 23:58:55.814 (UTC)
      Start time                   = 2004 MAY 01 00:00:00.000 (TDB)
      Stop time                    = 2004 MAY 04 23:58:55.814 (UTC)
      Stop time                    = 2004 MAY 05 00:00:00.000 (TDB)
      Elevation limit (degrees)    = 6.000000
 
   Loaded SPICE Kernels:
 
      Kernel name:    geomevnt.mk
      Kernel type:    META
      Kernel source:
 
      Kernel name:    de405s.bsp
      Kernel type:    SPK
      Kernel source:  geomevnt.mk
 
      Kernel name:    dsnstns.bsp
      Kernel type:    SPK
      Kernel source:  geomevnt.mk
 
      Kernel name:    DSN_topo.frm
      Kernel type:    TEXT
      Kernel source:  geomevnt.mk
 
      Kernel name:    earth_031228_231229_predict.bpc
      Kernel type:    PCK
      Kernel source:  geomevnt.mk
 
      Kernel name:    naif0007.tls
      Kernel type:    TEXT
      Kernel source:  geomevnt.mk
 
      Kernel name:    ORMM__040501000000_00069.BSP
      Kernel type:    SPK
      Kernel source:  geomevnt.mk
 
      Kernel name:    pck00007.tpc
      Kernel type:    TEXT
      Kernel source:  geomevnt.mk
 
   SRFAZL test results:
 
      ET = 2004 MAY 01 16:10:26.686 TDB
         Elevation (degrees):  5.9999951213034E+00
         Azimuth   (degrees):  6.3942106540558E+01
         Range     (km):       3.2480333721584E+08
 
      ET = 2004 MAY 01 16:10:26.688 TDB
         Elevation (degrees):  6.0000012274140E+00
         Azimuth   (degrees):  6.3942111060397E+01
         Range     (km):       3.2480333724230E+08
 
   OCCPT test results:
 
      ET = 2004 MAY 01 13:38:36.740 TDB.  Not occulted.
      ET = 2004 MAY 01 13:38:36.742 TDB.  Occulted.
 


Find Times when Target is Above the Elevation Limit







Task Statement




Extend the program of the previous chapter to find times when the target object is above the 6 degree elevation limit in the topocentric frame associated with the specified surface point.

Store the set of time intervals when the target is visible in a SPICELIB window. We'll call this the ``result window.''

Display each of the intervals in the result window as a pair of start and stop times. Express each time as a TDB calendar date using the same format as in the previous program.



Learning Goals




The student is introduced to topocentric frames and the SPICE calls required to compute azimuth and elevation of a vector expressed in a topocentric frame.

The student is introduced to a few of the most used SPICELIB cell, set and window routines.

The student is introduced to some elementary root-finding techniques. Some aspects of root-finding are widely applicable, such as selecting a search step size or a root-finding convergence tolerance. In addition, manipulation of SPICELIB cells, sets, and windows is quite useful for solving more complex root-finding problems, such as finding times when Boolean combinations of geometric criteria are satisfied.

The non-SPICELIB utility routines presented here are not emphasized, in part because these utilities unofficial, and more importantly because the planned SPICELIB geometric event finding software will replace the utilities presented here.



Approach






Solution steps



A possible solution would consist of the following steps:

            CALL FNDEVT ( ETBEG, ETEND, STEPSZ, STRIS, EVTSET )
            CALL MAKWIN ( ETBEG, ETEND, STRIS, EVTSET, RISWIN )


Solution






Solution Code



The example program below shows one possible solution.

 
         PROGRAM SCVIEW
         IMPLICIT NONE
   C
   C     Load SPICE Kernels and read additional inputs required to
   C     solve the spacecraft visibility problem posed by the
   C     geometric event finding hands-on code lesson.
   C
   C     The following items are loaded from a SPICE text kernel:
   C
   C        - Observer (surface point) name
   C        - Target name
   C        - Occulting body name
   C        - Aberration correction specification
   C        - Start time string
   C        - Stop time string
   C        - Elevation limit (in degrees)
   C
   C     This program assumes all necessary kernels and kernel
   C     variable definitions are specified in a SPICE meta-kernel.
   C     The meta-kernel name is assumed to be 'geomevnt.mk'.
   C
   C     Find and display the window of times when the target is
   C     above the elevation limit in the observer's topocentric
   C     reference frame.
   C
 
   C
   C     SPICELIB functions
   C
         DOUBLE PRECISION      DPR
         DOUBLE PRECISION      RPD
         INTEGER               CARDD
         LOGICAL               BADKPV
 
   C
   C     Other functions
   C
         EXTERNAL              STRIS
         LOGICAL               OCCPT
         LOGICAL               SETRIS
 
   C
   C     Local parameters
   C
 
   C
   C     The meta-kernel:
   C
         CHARACTER*(*)         META
         PARAMETER           ( META   = 'geomevnt.mk' )
 
   C
   C     Labels of kernel pool items:
   C
         CHARACTER*(*)         OCCLBL
         PARAMETER           ( OCCLBL = 'SCVIEW_OCCBDY' )
 
         CHARACTER*(*)         SRFLBL
         PARAMETER           ( SRFLBL = 'SCVIEW_SRFPT' )
 
         CHARACTER*(*)         TRGLBL
         PARAMETER           ( TRGLBL = 'SCVIEW_TARGET' )
 
         CHARACTER*(*)         CORLBL
         PARAMETER           ( CORLBL = 'SCVIEW_ABCORR' )
 
         CHARACTER*(*)         STRLBL
         PARAMETER           ( STRLBL = 'SCVIEW_START' )
 
         CHARACTER*(*)         STPLBL
         PARAMETER           ( STPLBL = 'SCVIEW_STOP' )
 
         CHARACTER*(*)         ELVLBL
         PARAMETER           ( ELVLBL = 'SCVIEW_ELVLIM' )
 
 
   C
   C     STEPSZ is the step size, measured in seconds, used to search
   C     for times bracketing a state transition.
   C
         DOUBLE PRECISION      STEPSZ
         PARAMETER           ( STEPSZ = 300.D0 )
 
   C
   C     Maximum number of events we can handle in our event set:
   C
         INTEGER               MAXEVT
         PARAMETER           ( MAXEVT = 1000 )
 
   C
   C     Maximum result window size:
   C
         INTEGER               MAXWIN
         PARAMETER           ( MAXWIN = 2 * MAXEVT )
 
   C
   C     SPICELIB cell bound:
   C
         INTEGER               LBCELL
         PARAMETER           ( LBCELL = -5 )
 
   C
   C     String length parameters:
   C
         INTEGER               CORLEN
         PARAMETER           ( CORLEN = 10 )
 
         INTEGER               FILSIZ
         PARAMETER           ( FILSIZ = 255 )
 
         INTEGER               LNSIZE
         PARAMETER           ( LNSIZE = 200 )
 
         INTEGER               NAMLEN
         PARAMETER           ( NAMLEN = 32 )
 
         INTEGER               TIMLEN
         PARAMETER           ( TIMLEN = 50 )
 
         INTEGER               TYPLEN
         PARAMETER           ( TYPLEN = 50 )
 
   C
   C     Local variables
   C
         CHARACTER*(CORLEN)    ABCORR
         CHARACTER*(FILSIZ)    FILE
         CHARACTER*(TYPLEN)    FILTYP
         CHARACTER*(LNSIZE)    LINE
         CHARACTER*(NAMLEN)    OCCBDY
         CHARACTER*(FILSIZ)    SOURCE
         CHARACTER*(NAMLEN)    SRFPT
         CHARACTER*(TIMLEN)    START
         CHARACTER*(TIMLEN)    STOP
         CHARACTER*(NAMLEN)    TARGET
         CHARACTER*(TIMLEN)    TIMSTR
         CHARACTER*(LNSIZE)    TITLE
         CHARACTER*(TIMLEN)    TSTTIM ( 2 )
 
 
         DOUBLE PRECISION      AZ
         DOUBLE PRECISION      EL
         DOUBLE PRECISION      ELVLIM
         DOUBLE PRECISION      ET
         DOUBLE PRECISION      ETBEG
         DOUBLE PRECISION      ETEND
         DOUBLE PRECISION      EVTSET ( LBCELL : MAXEVT )
         DOUBLE PRECISION      INTBEG
         DOUBLE PRECISION      INTEND
         DOUBLE PRECISION      R
         DOUBLE PRECISION      REVLIM
         DOUBLE PRECISION      RISWIN ( LBCELL : MAXWIN )
 
         INTEGER               HANDLE
         INTEGER               I
         INTEGER               N
         INTEGER               NKER
         INTEGER               WINSIZ
 
         LOGICAL               BAD
         LOGICAL               FOUND
         LOGICAL               HIDDEN
         LOGICAL               LVAL
 
   C
   C     Load the meta-kernel.
   C
         CALL FURNSH ( META )
 
   C
   C     Look up observation location, target,  name of the
   C     occulting body, aberration correction, start
   C     and stop times, and elevation limit.
   C
   C     We simplify error checking by using BADKPV.  BADKPV
   C     tests whether a kernel variable having specified
   C     attributes is present in the kernel pool. The fourth
   C     argument of BADKPV is the expected dimension; the
   C     last indicates the expected data type.  Read the
   C     header of BADKPV for details.
   C
         BAD = BADKPV ( 'SCVIEW', SRFLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', TRGLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', OCCLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', CORLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', STRLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', STPLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', ELVLBL, '=',  1,      1, 'N' )
 
   C
   C     Now we know the kernel variables of interest have been
   C     defined; look up the values.  We don't need to check
   C     the FOUND flag because BADKPV has ensured each item
   C     will be found.
   C
   C     Note that the elevation limit has units of degrees on input,
   C     so we store the input in the variable ELVLIM.  We convert
   C     this value to radians for computation using SPICE routines;
   C     we'll store the equivalent value in radians in REVLIM.
   C
         CALL GCPOOL ( SRFLBL, 1, 1,      N,    SRFPT,  FOUND )
         CALL GCPOOL ( TRGLBL, 1, 1,      N,    TARGET, FOUND )
         CALL GCPOOL ( OCCLBL, 1, 1,      N,    OCCBDY, FOUND )
         CALL GCPOOL ( CORLBL, 1, 1,      N,    ABCORR, FOUND )
         CALL GCPOOL ( STRLBL, 1, 1,      N,    START,  FOUND )
         CALL GCPOOL ( STPLBL, 1, 1,      N,    STOP,   FOUND )
         CALL GDPOOL ( ELVLBL, 1, 1,      N,    ELVLIM, FOUND )
 
         REVLIM = RPD() * ELVLIM
 
   C
   C     Display to standard output a banner for the output report:
   C
         TITLE  = 'Inputs for geometric event finding program'
         CALL TOSTDO ( ' '   )
         CALL TOSTDO ( TITLE )
         CALL TOSTDO ( ' '   )
 
         LINE   = '   Target                       = #'
         CALL REPMC  ( LINE, '#', TARGET, LINE )
         CALL TOSTDO ( LINE )
 
         LINE   = '   Observation surface location = #'
         CALL REPMC  ( LINE, '#', SRFPT, LINE )
         CALL TOSTDO ( LINE )
 
         LINE   = '   Aberration correction        = #'
         CALL REPMC  ( LINE, '#', ABCORR, LINE )
         CALL TOSTDO ( LINE )
 
   C
   C     Convert the start and stop times to ET.
   C
         CALL STR2ET ( START, ETBEG )
         CALL STR2ET ( STOP,  ETEND )
 
   C
   C     Display the start time as both calendar ET and
   C     calendar UTC using the formats shown below.
   C
   C        2004 MAY 06 20:15:00.000 (UTC)
   C        2004 MAY 06 20:15:00.000 (TDB)
   C
         CALL TIMOUT ( ETBEG,
        .              'YYYY MON DD HR:MN:SC.### (UTC)',
        .              TIMSTR                           )
 
         LINE   = '   Start time                   = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
 
         CALL TIMOUT ( ETBEG,
        .              'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .              TIMSTR                           )
 
         LINE   = '   Start time                   = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
 
         CALL TIMOUT ( ETEND,
        .              'YYYY MON DD HR:MN:SC.### (UTC)',
        .              TIMSTR                           )
 
         LINE   = '   Stop time                    = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
 
         CALL TIMOUT ( ETEND,
        .              'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .              TIMSTR                                   )
 
         LINE   = '   Stop time                    = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
   C
   C     Display the elevation limit in degrees.
   C
         LINE   = '   Elevation limit (degrees)    = #'
         CALL REPMF  ( LINE, '#', ELVLIM, 7, 'F', LINE )
         CALL TOSTDO ( LINE )
 
   C
   C     Display the names of the SPICE kernels we've loaded.
   C
         CALL TOSTDO ( ' '   )
         CALL TOSTDO ( 'Loaded SPICE Kernels: ' )
 
         CALL KTOTAL ( 'ALL', NKER )
 
         DO I = 1, NKER
 
            CALL KDATA ( I,      'ALL',  FILE, FILTYP,
        .                SOURCE, HANDLE, FOUND         )
   C
   C        Due to the way we've constructed the loop, there's
   C        no need to check the FOUND flag.
   C
            CALL TOSTDO ( ' '   )
 
            LINE = '   Kernel name:    #'
            CALL REPMC  ( LINE, '#', FILE, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '   Kernel type:    #'
            CALL REPMC  ( LINE, '#', FILTYP, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '   Kernel source:  #'
            CALL REPMC  ( LINE, '#', SOURCE, LINE )
            CALL TOSTDO ( LINE )
 
         END DO
 
   C
   C     Test SRFAZL:  make two calls at times bracketing
   C     a spacecraft "rise" event.
   C
         CALL TOSTDO ( ' ' )
         CALL TOSTDO ( 'SRFAZL test results:' )
         CALL TOSTDO ( ' ' )
 
         TSTTIM(1) = '2004 MAY 01 16:10:26.686 TDB'
         TSTTIM(2) = '2004 MAY 01 16:10:26.688 TDB'
 
         DO I = 1, 2
   C
   C        Convert the TDB calendar time to seconds past J2000.
   C        Find the azimuth, elevation, and range at this time.
   C
            CALL STR2ET ( TSTTIM(I), ET )
            CALL SRFAZL ( SRFPT,     ET, ABCORR, TARGET, AZ, EL, R )
   C
   C        Display the time and elevation.  Convert radians to
   C        degrees for output using the SPICELIB function DPR.
   C
            LINE = '   ET = #'
            CALL REPMC  ( LINE, '#', TSTTIM(I),   LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '      Elevation (degrees):  #'
            CALL REPMD  ( LINE, '#', DPR()*EL,  16, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '      Azimuth   (degrees):  #'
            CALL REPMD  ( LINE, '#', DPR()*AZ,  16, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '      Range     (km):       #'
            CALL REPMD  ( LINE, '#', R,         16, LINE )
            CALL TOSTDO ( LINE )
            CALL TOSTDO ( ' '  )
 
         END DO
 
 
   C
   C     Test OCCPT:  make two calls at times bracketing
   C     a spacecraft occultation ingress event.
   C
         CALL TOSTDO ( 'OCCPT test results:' )
         CALL TOSTDO ( ' ' )
 
         TSTTIM(1) =  '2004 MAY 01 13:38:36.740 TDB'
         TSTTIM(2) =  '2004 MAY 01 13:38:36.742 TDB'
 
         DO I = 1, 2
   C
   C        Convert the TDB calendar time to seconds past J2000.
   C
            CALL STR2ET ( TSTTIM(I), ET )
            HIDDEN = OCCPT ( TARGET, OCCBDY, ET, ABCORR, SRFPT )
   C
   C        Display the time and occultation state.
   C
            LINE = '   ET = #.  #.'
            CALL REPMC  ( LINE, '#', TSTTIM(I),   LINE )
 
            IF ( HIDDEN ) THEN
               CALL REPMC ( LINE, '#', 'Occulted', LINE )
            ELSE
               CALL REPMC ( LINE, '#', 'Not occulted', LINE )
            END IF
 
            CALL TOSTDO ( LINE )
 
         END DO
 
         CALL TOSTDO ( ' ' )
 
   C
   C     Get our initializations out of the way now.  Every SPICELIB
   C     set or window must have its size initialized.
   C
   C     Initialize the event set.
   C
         CALL SSIZED ( MAXEVT, EVTSET )
 
   C
   C     Initialize the rise/set window.
   C
         CALL SSIZED ( MAXWIN, RISWIN )
 
   C
   C     Set the parameters required by the rise/set state routine.
   C     This function works by side effects; the returned value
   C     LVAL is not used.
   C
         LVAL = SETRIS ( SRFPT, ABCORR, TARGET, REVLIM )
 
   C
   C     Locate the epochs where the target passes through the
   C     elevation limit. The logical function STRIS determines
   C     the "state" of the target:  .TRUE. when the target's
   C     elevation is positive, .FALSE. otherwise.
   C     The step size is given in seconds.
   C
         CALL FNDEVT ( ETBEG, ETEND, STEPSZ, STRIS, EVTSET )
 
   C
   C     Use the transition set to create a window of time intervals
   C     when the target is above the elevation limit at the
   C     surface point.
   C
         CALL MAKWIN ( ETBEG, ETEND, STRIS, EVTSET, RISWIN )
 
   C
   C     Display the rise and set times.
   C
         LINE = 'Rise and set times of # as seen from #:'
         CALL REPMC ( LINE, '#', TARGET, LINE )
         CALL REPMC ( LINE, '#', SRFPT,  LINE )
 
         CALL TOSTDO ( LINE )
         CALL TOSTDO ( ' '  )
 
         WINSIZ = CARDD(RISWIN) / 2
 
         DO I = 1, WINSIZ
   C
   C        Fetch the Ith interval from the window.
   C
            CALL WNFETD ( RISWIN, I, INTBEG, INTEND )
 
   C
   C        Convert the transition time to a TDB calendar string.
   C
            CALL TIMOUT ( INTBEG,
        .                 'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .                 TIMSTR                                    )
 
   C
   C        Write the string to standard output.
   C
            IF ( I .EQ. 1 ) THEN
               LINE = 'Target rise time (or window start):  #'
            ELSE
               LINE = 'Target rise time:                    #'
            END IF
 
            CALL REPMC ( LINE, '#', TIMSTR, LINE )
            CALL TOSTDO ( LINE )
 
 
            CALL TIMOUT ( INTEND,
        .                 'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .                 TIMSTR                                    )
 
   C
   C        Write the string to standard output.
   C
            IF ( I .EQ. WINSIZ ) THEN
               LINE = 'Target set time (or window end):     #'
            ELSE
               LINE = 'Target set time:                     #'
            END IF
 
            CALL REPMC ( LINE, '#', TIMSTR, LINE )
            CALL TOSTDO ( LINE )
 
            CALL TOSTDO ( ' ' )
 
         END DO
 
         END


Solution Sample Output



Numerical results shown for this example may differ across platforms since the results depend on the SPICE kernels used as input and on the host platform's arithmetic implementation.

After compiling the program, execute it. The output is:

 
   Inputs for geometric event finding program
 
      Target                       = MEX
      Observation surface location = DSS-14
      Occulting body               = MARS
      Aberration correction        = CN+S
      Start time                   = 2004 APR 30 23:58:55.814 (UTC)
      Start time                   = 2004 MAY 01 00:00:00.000 (TDB)
      Stop time                    = 2004 MAY 04 23:58:55.814 (UTC)
      Stop time                    = 2004 MAY 05 00:00:00.000 (TDB)
      Elevation limit (degrees)    = 6.000000
 
   Loaded SPICE Kernels:
 
      Kernel name:    geomevnt.mk
      Kernel type:    META
      Kernel source:
 
      Kernel name:    de405s.bsp
      Kernel type:    SPK
      Kernel source:  geomevnt.mk
 
      Kernel name:    dsnstns.bsp
      Kernel type:    SPK
      Kernel source:  geomevnt.mk
 
      Kernel name:    DSN_topo.frm
      Kernel type:    TEXT
      Kernel source:  geomevnt.mk
 
      Kernel name:    earth_031228_231229_predict.bpc
      Kernel type:    PCK
      Kernel source:  geomevnt.mk
 
      Kernel name:    naif0007.tls
      Kernel type:    TEXT
      Kernel source:  geomevnt.mk
 
      Kernel name:    ORMM__040501000000_00069.BSP
      Kernel type:    SPK
      Kernel source:  geomevnt.mk
 
      Kernel name:    pck00007.tpc
      Kernel type:    TEXT
      Kernel source:  geomevnt.mk
 
   SRFAZL test results:
 
      ET = 2004 MAY 01 16:10:26.686 TDB
         Elevation (degrees):  5.9999951213034E+00
         Azimuth   (degrees):  6.3942106540558E+01
         Range     (km):       3.2480333721584E+08
 
      ET = 2004 MAY 01 16:10:26.688 TDB
         Elevation (degrees):  6.0000012274140E+00
         Azimuth   (degrees):  6.3942111060397E+01
         Range     (km):       3.2480333724230E+08
 
   OCCPT test results:
 
      ET = 2004 MAY 01 13:38:36.740 TDB.  Not occulted.
      ET = 2004 MAY 01 13:38:36.742 TDB.  Occulted.
 
   Rise and set times of MEX as seen from DSS-14:
 
   Target rise time (or window start):  2004 MAY 01 00:00:00.000 (TDB)
   Target set time:                     2004 MAY 01 05:36:08.338 (TDB)
 
   Target rise time:                    2004 MAY 01 16:10:26.687 (TDB)
   Target set time:                     2004 MAY 02 05:35:03.020 (TDB)
 
   Target rise time:                    2004 MAY 02 16:09:14.024 (TDB)
   Target set time:                     2004 MAY 03 05:33:57.180 (TDB)
 
   Target rise time:                    2004 MAY 03 16:08:02.224 (TDB)
   Target set time:                     2004 MAY 04 05:32:50.686 (TDB)
 
   Target rise time:                    2004 MAY 04 16:06:51.203 (TDB)
   Target set time (or window end):     2004 MAY 05 00:00:00.000 (TDB)
 


Find Times when Target is Visible







Task Statement




Extend the program of the previous chapter to find times when the target object is:

Store the set of time intervals when the target is visible in a SPICELIB window. We'll call this the ``result window.''

Display each of the intervals in the result window as a pair of start and stop times. Express each time as a TDB calendar date using the same format as in the previous program.



Learning Goals




The student gains further experience with the SPICELIB cell, set and window routines, and with root finding techniques.



Approach






Solution steps



A possible solution would consist of the following steps:

            CALL FNDEVT ( ETBEG, ETEND, STEPSZ, STOCP, EVTSET )
            CALL MAKWIN ( ETBEG, ETEND, STOCP, EVTSET, OCCWIN )
This completes the assignment.



Solution






Solution Code



 
         PROGRAM SCVIEW
         IMPLICIT NONE
   C
   C     Load SPICE Kernels and read additional inputs required to
   C     solve the spacecraft visibility problem posed by the
   C     geometric event finding hands-on code lesson.
   C
   C     The following items are loaded from a SPICE text kernel:
   C
   C        - Observer (surface point) name
   C        - Target name
   C        - Occulting body name
   C        - Aberration correction specification
   C        - Start time string
   C        - Stop time string
   C        - Elevation limit (in degrees)
   C
   C     This program assumes all necessary kernels and kernel
   C     variable definitions are specified in a SPICE meta-kernel.
   C     The meta-kernel name is assumed to be 'geomevnt.mk'.
   C
   C     Find and display the window of times when the target is
   C     above the elevation limit in the observer's topocentric
   C     reference frame.
   C
   C     Find and display the window of times when the target is
   C     occulted as seen from the observer's location.
   C
   C     Find and display the window of times when the target is
   C     above the elevation limit and not occulted.
   C
 
   C
   C     SPICELIB functions
   C
         DOUBLE PRECISION      DPR
         DOUBLE PRECISION      RPD
 
         INTEGER               CARDD
 
         LOGICAL               BADKPV
 
   C
   C     Other functions
   C
         EXTERNAL              STOCP
         EXTERNAL              STRIS
 
         LOGICAL               OCCPT
         LOGICAL               SETOCP
         LOGICAL               SETRIS
 
   C
   C     Local parameters
   C
 
   C
   C     The meta-kernel:
   C
         CHARACTER*(*)         META
         PARAMETER           ( META   = 'geomevnt.mk' )
 
   C
   C     Labels of kernel pool items:
   C
         CHARACTER*(*)         OCCLBL
         PARAMETER           ( OCCLBL = 'SCVIEW_OCCBDY' )
 
         CHARACTER*(*)         SRFLBL
         PARAMETER           ( SRFLBL = 'SCVIEW_SRFPT' )
 
         CHARACTER*(*)         TRGLBL
         PARAMETER           ( TRGLBL = 'SCVIEW_TARGET' )
 
         CHARACTER*(*)         CORLBL
         PARAMETER           ( CORLBL = 'SCVIEW_ABCORR' )
 
         CHARACTER*(*)         STRLBL
         PARAMETER           ( STRLBL = 'SCVIEW_START' )
 
         CHARACTER*(*)         STPLBL
         PARAMETER           ( STPLBL = 'SCVIEW_STOP' )
 
         CHARACTER*(*)         ELVLBL
         PARAMETER           ( ELVLBL = 'SCVIEW_ELVLIM' )
 
 
   C
   C     STEPSZ is the step size, measured in seconds, used to search
   C     for times bracketing a state transition.
   C
         DOUBLE PRECISION      STEPSZ
         PARAMETER           ( STEPSZ = 300.D0 )
 
   C
   C     Maximum number of events we can handle in our event set:
   C
         INTEGER               MAXEVT
         PARAMETER           ( MAXEVT = 1000 )
 
   C
   C     Maximum result window size:
   C
         INTEGER               MAXWIN
         PARAMETER           ( MAXWIN = 2 * MAXEVT )
 
   C
   C     SPICELIB cell bound:
   C
         INTEGER               LBCELL
         PARAMETER           ( LBCELL = -5 )
 
   C
   C     String length parameters:
   C
         INTEGER               CORLEN
         PARAMETER           ( CORLEN = 10 )
 
         INTEGER               FILSIZ
         PARAMETER           ( FILSIZ = 255 )
 
         INTEGER               LNSIZE
         PARAMETER           ( LNSIZE = 200 )
 
         INTEGER               NAMLEN
         PARAMETER           ( NAMLEN = 32 )
 
         INTEGER               TIMLEN
         PARAMETER           ( TIMLEN = 50 )
 
         INTEGER               TYPLEN
         PARAMETER           ( TYPLEN = 50 )
 
   C
   C     Local variables
   C
         CHARACTER*(CORLEN)    ABCORR
         CHARACTER*(FILSIZ)    FILE
         CHARACTER*(TYPLEN)    FILTYP
         CHARACTER*(LNSIZE)    LINE
         CHARACTER*(NAMLEN)    OCCBDY
         CHARACTER*(FILSIZ)    SOURCE
         CHARACTER*(NAMLEN)    SRFPT
         CHARACTER*(TIMLEN)    START
         CHARACTER*(TIMLEN)    STOP
         CHARACTER*(NAMLEN)    TARGET
         CHARACTER*(TIMLEN)    TIMSTR
         CHARACTER*(LNSIZE)    TITLE
         CHARACTER*(TIMLEN)    TSTTIM ( 2 )
 
 
         DOUBLE PRECISION      AZ
         DOUBLE PRECISION      EL
         DOUBLE PRECISION      ELVLIM
         DOUBLE PRECISION      ET
         DOUBLE PRECISION      ETBEG
         DOUBLE PRECISION      ETEND
         DOUBLE PRECISION      EVTSET ( LBCELL : MAXEVT )
         DOUBLE PRECISION      INTBEG
         DOUBLE PRECISION      INTEND
         DOUBLE PRECISION      OCCWIN ( LBCELL : MAXWIN )
         DOUBLE PRECISION      R
         DOUBLE PRECISION      REVLIM
         DOUBLE PRECISION      RISWIN ( LBCELL : MAXWIN )
         DOUBLE PRECISION      VISWIN ( LBCELL : MAXWIN )
 
         INTEGER               HANDLE
         INTEGER               I
         INTEGER               N
         INTEGER               NKER
         INTEGER               WINSIZ
 
         LOGICAL               BAD
         LOGICAL               FOUND
         LOGICAL               HIDDEN
         LOGICAL               LVAL
 
   C
   C     Load the meta-kernel.
   C
         CALL FURNSH ( META )
 
   C
   C     Look up observation location, target,  name of the
   C     occulting body, aberration correction, start
   C     and stop times, and elevation limit.
   C
   C     We simplify error checking by using BADKPV.  BADKPV
   C     tests whether a kernel variable having specified
   C     attributes is present in the kernel pool. The fourth
   C     argument of BADKPV is the expected dimension; the
   C     last indicates the expected data type.  Read the
   C     header of BADKPV for details.
   C
         BAD = BADKPV ( 'SCVIEW', SRFLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', TRGLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', OCCLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', CORLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', STRLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', STPLBL, '=',  1,      1, 'C' )
         BAD = BADKPV ( 'SCVIEW', ELVLBL, '=',  1,      1, 'N' )
 
   C
   C     Now we know the kernel variables of interest have been
   C     defined; look up the values.  We don't need to check
   C     the FOUND flag because BADKPV has ensured the item
   C     will be found.
   C
   C     Note that the elevation limit has units of degrees on input,
   C     so we store the input in the variable ELVLIM.  We convert
   C     this value to radians for computation using SPICE routines;
   C     we'll store the equivalent value in radians in REVLIM.
   C
         CALL GCPOOL ( SRFLBL, 1, 1,      N,    SRFPT,  FOUND )
         CALL GCPOOL ( TRGLBL, 1, 1,      N,    TARGET, FOUND )
         CALL GCPOOL ( OCCLBL, 1, 1,      N,    OCCBDY, FOUND )
         CALL GCPOOL ( CORLBL, 1, 1,      N,    ABCORR, FOUND )
         CALL GCPOOL ( STRLBL, 1, 1,      N,    START,  FOUND )
         CALL GCPOOL ( STPLBL, 1, 1,      N,    STOP,   FOUND )
         CALL GDPOOL ( ELVLBL, 1, 1,      N,    ELVLIM, FOUND )
 
         REVLIM = RPD() * ELVLIM
 
   C
   C     Display to standard output a banner for the output report:
   C
         TITLE  = 'Inputs for geometric event finding program'
         CALL TOSTDO ( ' '   )
         CALL TOSTDO ( TITLE )
         CALL TOSTDO ( ' '   )
 
         LINE   = '   Target                       = #'
         CALL REPMC  ( LINE, '#', TARGET, LINE )
         CALL TOSTDO ( LINE )
 
         LINE   = '   Observation surface location = #'
         CALL REPMC  ( LINE, '#', SRFPT, LINE )
         CALL TOSTDO ( LINE )
 
         LINE   = '   Aberration correction        = #'
         CALL REPMC  ( LINE, '#', ABCORR, LINE )
         CALL TOSTDO ( LINE )
 
   C
   C     Convert the start and stop times to ET.
   C
         CALL STR2ET ( START, ETBEG )
         CALL STR2ET ( STOP,  ETEND )
 
   C
   C     Display the start time as both calendar ET and
   C     calendar UTC using the formats shown below.
   C
   C        2004 MAY 06 20:15:00.000 (UTC)
   C        2004 MAY 06 20:15:00.000 (TDB)
   C
         CALL TIMOUT ( ETBEG,
        .              'YYYY MON DD HR:MN:SC.### (UTC)',
        .              TIMSTR                           )
 
         LINE   = '   Start time                   = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
 
         CALL TIMOUT ( ETBEG,
        .              'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .              TIMSTR                           )
 
         LINE   = '   Start time                   = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
 
         CALL TIMOUT ( ETEND,
        .              'YYYY MON DD HR:MN:SC.### (UTC)',
        .              TIMSTR                           )
 
         LINE   = '   Stop time                    = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
 
         CALL TIMOUT ( ETEND,
        .              'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .              TIMSTR                                   )
 
         LINE   = '   Stop time                    = #'
         CALL REPMC  ( LINE, '#', TIMSTR, LINE )
         CALL TOSTDO ( LINE )
 
   C
   C     Display the elevation limit in degrees.
   C
         LINE   = '   Elevation limit (degrees)    = #'
         CALL REPMF  ( LINE, '#', ELVLIM, 7, 'F', LINE )
         CALL TOSTDO ( LINE )
 
   C
   C     Display the names of the SPICE kernels we've loaded.
   C
         CALL TOSTDO ( ' '   )
         CALL TOSTDO ( 'Loaded SPICE Kernels: ' )
 
         CALL KTOTAL ( 'ALL', NKER )
 
         DO I = 1, NKER
 
            CALL KDATA ( I,      'ALL',  FILE, FILTYP,
        .                SOURCE, HANDLE, FOUND         )
   C
   C        Due to the way we've constructed the loop, there's
   C        no need to check the FOUND flag.
   C
            CALL TOSTDO ( ' '   )
 
            LINE = '   Kernel name:    #'
            CALL REPMC  ( LINE, '#', FILE, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '   Kernel type:    #'
            CALL REPMC  ( LINE, '#', FILTYP, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '   Kernel source:  #'
            CALL REPMC  ( LINE, '#', SOURCE, LINE )
            CALL TOSTDO ( LINE )
 
         END DO
 
   C
   C     Test SRFAZL:  make two calls at times bracketing
   C     a spacecraft "rise" event.
   C
         CALL TOSTDO ( ' ' )
         CALL TOSTDO ( 'SRFAZL test results:' )
         CALL TOSTDO ( ' ' )
 
         TSTTIM(1) = '2004 MAY 01 16:10:26.686 TDB'
         TSTTIM(2) = '2004 MAY 01 16:10:26.688 TDB'
 
         DO I = 1, 2
   C
   C        Convert the TDB calendar time to seconds past J2000.
   C        Find the azimuth, elevation, and range at this time.
   C
            CALL STR2ET ( TSTTIM(I), ET )
            CALL SRFAZL ( SRFPT,     ET, ABCORR, TARGET, AZ, EL, R )
   C
   C        Display the time and elevation.  Convert radians to
   C        degrees for output using the SPICELIB function DPR.
   C
            LINE = '   ET = #'
            CALL REPMC  ( LINE, '#', TSTTIM(I),   LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '      Elevation (degrees):  #'
            CALL REPMD  ( LINE, '#', DPR()*EL,  16, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '      Azimuth   (degrees):  #'
            CALL REPMD  ( LINE, '#', DPR()*AZ,  16, LINE )
            CALL TOSTDO ( LINE )
 
            LINE = '      Range     (km):       #'
            CALL REPMD  ( LINE, '#', R,         16, LINE )
            CALL TOSTDO ( LINE )
            CALL TOSTDO ( ' '  )
 
         END DO
 
 
   C
   C     Test OCCPT:  make two calls at times bracketing
   C     a spacecraft occultation ingress event.
   C
         CALL TOSTDO ( 'OCCPT test results:' )
         CALL TOSTDO ( ' ' )
 
         TSTTIM(1) =  '2004 MAY 01 13:38:36.740 TDB'
         TSTTIM(2) =  '2004 MAY 01 13:38:36.742 TDB'
 
         DO I = 1, 2
   C
   C        Convert the TDB calendar time to seconds past J2000.
   C
            CALL STR2ET ( TSTTIM(I), ET )
            HIDDEN = OCCPT ( TARGET, OCCBDY, ET, ABCORR, SRFPT )
   C
   C        Display the time and occultation state.
   C
            LINE = '   ET = #.  #.'
            CALL REPMC  ( LINE, '#', TSTTIM(I),   LINE )
 
            IF ( HIDDEN ) THEN
               CALL REPMC ( LINE, '#', 'Occulted', LINE )
            ELSE
               CALL REPMC ( LINE, '#', 'Not occulted', LINE )
            END IF
 
            CALL TOSTDO ( LINE )
 
         END DO
 
         CALL TOSTDO ( ' ' )
 
   C
   C     Get our initializations out of the way now.  Every SPICELIB
   C     set or window must have its initial size set.
   C
   C     Initialize the event set.
   C
         CALL SSIZED ( MAXEVT, EVTSET )
 
   C
   C     Initialize the rise/set window.
   C
         CALL SSIZED ( MAXWIN, RISWIN )
 
   C
   C     Set the parameters required by the rise/set state routine.
   C     This function works by side effects; the returned value
   C     LVAL is not used.
   C
         LVAL = SETRIS ( SRFPT, ABCORR, TARGET, REVLIM )
 
   C
   C     Locate the epochs where the target passes through the
   C     elevation limit. The logical function STRIS determines
   C     the "state" of the target:  .TRUE. when the target's
   C     elevation is positive, .FALSE. otherwise.
   C     The step size is given in seconds.
   C
         CALL FNDEVT ( ETBEG, ETEND, STEPSZ, STRIS, EVTSET )
 
   C
   C     Use the transition set to create a window of time intervals
   C     when the target is above the elevation limit at the
   C     surface point.
   C
         CALL MAKWIN ( ETBEG, ETEND, STRIS, EVTSET, RISWIN )
 
   C
   C     Display the rise and set times.
   C
         CALL TOSTDO ( ' ' )
         CALL TOSTDO ( ' ' )
 
         LINE = 'Rise and set times of # as seen from #:'
         CALL REPMC ( LINE, '#', TARGET, LINE )
         CALL REPMC ( LINE, '#', SRFPT,  LINE )
 
         CALL TOSTDO ( LINE )
         CALL TOSTDO ( ' '  )
 
         WINSIZ = CARDD(RISWIN) / 2
 
         DO I = 1, WINSIZ
   C
   C        Fetch the Ith interval from the window.
   C
            CALL WNFETD ( RISWIN, I, INTBEG, INTEND )
 
   C
   C        Convert the transition time to a TDB calendar string.
   C
            CALL TIMOUT ( INTBEG,
        .                 'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .                 TIMSTR                                    )
 
   C
   C        Write the string to standard output.
   C
            IF ( I .EQ. 1 ) THEN
               LINE = 'Target rise time (or window start):  #'
            ELSE
               LINE = 'Target rise time:                    #'
            END IF
 
            CALL REPMC ( LINE, '#', TIMSTR, LINE )
            CALL TOSTDO ( LINE )
 
 
            CALL TIMOUT ( INTEND,
        .                 'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .                 TIMSTR                                    )
 
   C
   C        Write the string to standard output.
   C
            IF ( I .EQ. WINSIZ ) THEN
               LINE = 'Target set time (or window end):     #'
            ELSE
               LINE = 'Target set time:                     #'
            END IF
 
            CALL REPMC ( LINE, '#', TIMSTR, LINE )
            CALL TOSTDO ( LINE )
 
            CALL TOSTDO ( ' ' )
 
         END DO
 
   C
   C     Now find the window of times when the target is occulted
   C     by the designated blocking body.
   C
   C     Initialize the occultation window.
   C
         CALL SSIZED ( MAXWIN, OCCWIN )
 
   C
   C     Initialize the occultation finding utilities. This function
   C     works by side effects; the returned value LVAL is not used.
   C
         LVAL = SETOCP ( TARGET, OCCBDY, ABCORR, SRFPT )
 
   C
   C     Locate the epochs where the target passes into or out
   C     of occultation by the occulting body.  These epochs
   C     will be returned as a SPICELIB set.
   C
         CALL FNDEVT ( ETBEG, ETEND, STEPSZ, STOCP, EVTSET )
 
   C
   C     Make the set into a window of times when the target
   C     is occulted.
   C
         CALL MAKWIN ( ETBEG, ETEND, STOCP, EVTSET, OCCWIN )
 
   C
   C     Display the set of times when the target is occulted.
   C
         CALL TOSTDO ( ' ' )
         CALL TOSTDO ( ' ' )
 
         LINE = '# occultation of # ---'
         CALL REPMC ( LINE, '#', OCCBDY, LINE )
         CALL REPMC ( LINE, '#', TARGET, LINE )
         CALL TOSTDO ( LINE )
 
         LINE = 'Ingress and egress times as seen from #:'
         CALL REPMC ( LINE, '#', SRFPT,  LINE )
         CALL TOSTDO ( LINE )
         CALL TOSTDO ( ' '  )
 
         WINSIZ = CARDD(OCCWIN) / 2
 
         DO I = 1, WINSIZ
   C
   C        Fetch the Ith interval from the window.
   C
            CALL WNFETD ( OCCWIN, I, INTBEG, INTEND )
 
   C
   C        Convert the transition time to a TDB calendar string.
   C
            CALL TIMOUT ( INTBEG,
        .                 'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .                 TIMSTR                                    )
 
   C
   C        Write the string to standard output.
   C
            IF ( I .EQ. 1 ) THEN
               LINE = 'Occ. ingress time (or window start): #'
            ELSE
               LINE = 'Occultation ingress time:            #'
            END IF
 
            CALL REPMC ( LINE, '#', TIMSTR, LINE )
            CALL TOSTDO ( LINE )
 
 
            CALL TIMOUT ( INTEND,
        .                 'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .                 TIMSTR                                    )
 
   C
   C        Write the string to standard output.
   C
            IF ( I .EQ. WINSIZ ) THEN
               LINE = 'Occ. egress time (or window end):    #'
            ELSE
               LINE = 'Occultation egress time:             #'
            END IF
 
            CALL REPMC ( LINE, '#', TIMSTR, LINE )
            CALL TOSTDO ( LINE )
 
            CALL TOSTDO ( ' ' )
 
         END DO
 
   C
   C     Find the visibility window by subtracting the occultation
   C     window from the rise-set window.  Display the result.
   C
   C     Initialize the visibility window.
   C
         CALL SSIZED ( MAXWIN, VISWIN )
 
   C
   C     Subtract the occultation window from the rise-set window.
   C
         CALL WNDIFD ( RISWIN, OCCWIN, VISWIN )
 
   C
   C     Display the set of times when the target is visible.
   C
         CALL TOSTDO ( ' ' )
         CALL TOSTDO ( ' ' )
         LINE = '# visibility start and stop times ' //
        .       'as seen from #:'
         CALL REPMC ( LINE, '#', TARGET, LINE )
         CALL REPMC ( LINE, '#', SRFPT,  LINE )
 
         CALL TOSTDO ( LINE )
         CALL TOSTDO ( ' '  )
 
         WINSIZ = CARDD(VISWIN) / 2
 
         DO I = 1, WINSIZ
   C
   C        Fetch the Ith interval from the window.
   C
            CALL WNFETD ( VISWIN, I, INTBEG, INTEND )
 
   C
   C        Convert the transition time to a TDB calendar string.
   C
            CALL TIMOUT ( INTBEG,
        .                 'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .                 TIMSTR                                    )
 
   C
   C        Write the string to standard output.
   C
            IF ( I .EQ. 1 ) THEN
               LINE = 'Vis. start time (or window start):   #'
            ELSE
               LINE = 'Visibility start time:               #'
            END IF
 
            CALL REPMC ( LINE, '#', TIMSTR, LINE )
            CALL TOSTDO ( LINE )
 
 
            CALL TIMOUT ( INTEND,
        .                 'YYYY MON DD HR:MN:SC.### (TDB) ::TDB',
        .                 TIMSTR                                    )
 
   C
   C        Write the string to standard output.
   C
            IF ( I .EQ. WINSIZ ) THEN
               LINE = 'Vis. end time (or window end):       #'
            ELSE
               LINE = 'Visibility end time:                 #'
            END IF
 
            CALL REPMC ( LINE, '#', TIMSTR, LINE )
            CALL TOSTDO ( LINE )
 
            CALL TOSTDO ( ' ' )
 
         END DO
 
         END


Solution Sample Output



Numerical results shown for this example may differ across platforms since the results depend on the SPICE kernels used as input and on the host platform's arithmetic implementation.

After compiling the program, execute it. The output is:

   Inputs for geometric event finding program
 
      Target                       = MEX
      Observation surface location = DSS-14
      Occulting body               = MARS
      Aberration correction        = CN+S
      Start time                   = 2004 APR 30 23:58:55.814 (UTC)
      Start time                   = 2004 MAY 01 00:00:00.000 (TDB)
      Stop time                    = 2004 MAY 04 23:58:55.814 (UTC)
      Stop time                    = 2004 MAY 05 00:00:00.000 (TDB)
      Elevation limit (degrees)    = 6.000000
 
   Loaded SPICE Kernels:
 
      Kernel name:    geomevnt.mk
      Kernel type:    META
      Kernel source:
 
      Kernel name:    de405s.bsp
      Kernel type:    SPK
      Kernel source:  geomevnt.mk
 
      Kernel name:    dsnstns.bsp
      Kernel type:    SPK
      Kernel source:  geomevnt.mk
 
      Kernel name:    DSN_topo.frm
      Kernel type:    TEXT
      Kernel source:  geomevnt.mk
 
      Kernel name:    earth_031228_231229_predict.bpc
      Kernel type:    PCK
      Kernel source:  geomevnt.mk
 
      Kernel name:    naif0007.tls
      Kernel type:    TEXT
      Kernel source:  geomevnt.mk
 
      Kernel name:    ORMM__040501000000_00069.BSP
      Kernel type:    SPK
      Kernel source:  geomevnt.mk
 
      Kernel name:    pck00007.tpc
      Kernel type:    TEXT
      Kernel source:  geomevnt.mk
 
   SRFAZL test results:
 
      ET = 2004 MAY 01 16:10:26.686 TDB
         Elevation (degrees):  5.9999951213034E+00
         Azimuth   (degrees):  6.3942106540558E+01
         Range     (km):       3.2480333721584E+08
 
      ET = 2004 MAY 01 16:10:26.688 TDB
         Elevation (degrees):  6.0000012274140E+00
         Azimuth   (degrees):  6.3942111060397E+01
         Range     (km):       3.2480333724230E+08
 
   OCCPT test results:
 
      ET = 2004 MAY 01 13:38:36.740 TDB.  Not occulted.
      ET = 2004 MAY 01 13:38:36.742 TDB.  Occulted.
 
   Rise and set times of MEX as seen from DSS-14:
 
   Target rise time (or window start):  2004 MAY 01 00:00:00.000 (TDB)
   Target set time:                     2004 MAY 01 05:36:08.338 (TDB)
 
   Target rise time:                    2004 MAY 01 16:10:26.687 (TDB)
   Target set time:                     2004 MAY 02 05:35:03.020 (TDB)
 
   Target rise time:                    2004 MAY 02 16:09:14.024 (TDB)
   Target set time:                     2004 MAY 03 05:33:57.180 (TDB)
 
   Target rise time:                    2004 MAY 03 16:08:02.224 (TDB)
   Target set time:                     2004 MAY 04 05:32:50.686 (TDB)
 
   Target rise time:                    2004 MAY 04 16:06:51.203 (TDB)
   Target set time (or window end):     2004 MAY 05 00:00:00.000 (TDB)
 
 
 
   MARS occultation of MEX ---
   Ingress and egress times as seen from DSS-14:
 
   Occ. ingress time (or window start): 2004 MAY 01 06:03:15.313 (TDB)
   Occultation egress time:             2004 MAY 01 07:05:45.105 (TDB)
 
   Occultation ingress time:            2004 MAY 01 13:38:36.741 (TDB)
   Occultation egress time:             2004 MAY 01 14:40:51.279 (TDB)
 
   Occultation ingress time:            2004 MAY 01 21:14:05.206 (TDB)
   Occultation egress time:             2004 MAY 01 22:16:04.347 (TDB)
 
   Occultation ingress time:            2004 MAY 02 04:49:26.607 (TDB)
   Occultation egress time:             2004 MAY 02 05:51:12.175 (TDB)
 
   Occultation ingress time:            2004 MAY 02 12:24:51.007 (TDB)
   Occultation egress time:             2004 MAY 02 13:26:20.244 (TDB)
 
   Occultation ingress time:            2004 MAY 02 20:00:17.624 (TDB)
   Occultation egress time:             2004 MAY 02 21:01:31.988 (TDB)
 
   Occultation ingress time:            2004 MAY 03 03:35:37.097 (TDB)
   Occultation egress time:             2004 MAY 03 04:36:36.001 (TDB)
 
   Occultation ingress time:            2004 MAY 03 11:11:00.963 (TDB)
   Occultation egress time:             2004 MAY 03 12:11:43.066 (TDB)
 
   Occultation ingress time:            2004 MAY 03 18:46:20.320 (TDB)
   Occultation egress time:             2004 MAY 03 19:46:47.607 (TDB)
 
   Occultation ingress time:            2004 MAY 04 02:21:38.679 (TDB)
   Occultation egress time:             2004 MAY 04 03:21:49.169 (TDB)
 
   Occultation ingress time:            2004 MAY 04 09:57:02.737 (TDB)
   Occultation egress time:             2004 MAY 04 10:56:56.827 (TDB)
 
   Occultation ingress time:            2004 MAY 04 17:32:19.374 (TDB)
   Occ. egress time (or window end):    2004 MAY 04 18:31:58.252 (TDB)
 
 
 
   MEX visibility start and stop times as seen from DSS-14:
 
   Vis. start time (or window start):   2004 MAY 01 00:00:00.000 (TDB)
   Visibility end time:                 2004 MAY 01 05:36:08.338 (TDB)
 
   Visibility start time:               2004 MAY 01 16:10:26.687 (TDB)
   Visibility end time:                 2004 MAY 01 21:14:05.206 (TDB)
 
   Visibility start time:               2004 MAY 01 22:16:04.347 (TDB)
   Visibility end time:                 2004 MAY 02 04:49:26.607 (TDB)
 
   Visibility start time:               2004 MAY 02 16:09:14.024 (TDB)
   Visibility end time:                 2004 MAY 02 20:00:17.624 (TDB)
 
   Visibility start time:               2004 MAY 02 21:01:31.988 (TDB)
   Visibility end time:                 2004 MAY 03 03:35:37.097 (TDB)
 
   Visibility start time:               2004 MAY 03 04:36:36.001 (TDB)
   Visibility end time:                 2004 MAY 03 05:33:57.180 (TDB)
 
   Visibility start time:               2004 MAY 03 16:08:02.224 (TDB)
   Visibility end time:                 2004 MAY 03 18:46:20.320 (TDB)
 
   Visibility start time:               2004 MAY 03 19:46:47.607 (TDB)
   Visibility end time:                 2004 MAY 04 02:21:38.679 (TDB)
 
   Visibility start time:               2004 MAY 04 03:21:49.169 (TDB)
   Visibility end time:                 2004 MAY 04 05:32:50.686 (TDB)
 
   Visibility start time:               2004 MAY 04 16:06:51.203 (TDB)
   Visibility end time:                 2004 MAY 04 17:32:19.374 (TDB)
 
   Visibility start time:               2004 MAY 04 18:31:58.252 (TDB)
   Vis. end time (or window end):       2004 MAY 05 00:00:00.000 (TDB)
 


Utilities







Rise/Set State Function Module MODRIS




 
         LOGICAL FUNCTION MODRIS ( SRFPT, ET, ABCORR, TARGET, LIMIT )
         IMPLICIT NONE
   C
   C     Umbrella routine for rise-set function.
   C
         CHARACTER*(*)         SRFPT
         DOUBLE PRECISION      ET
         CHARACTER*(*)         ABCORR
         CHARACTER*(*)         TARGET
         DOUBLE PRECISION      LIMIT
 
   C
   C     Entry points
   C
         LOGICAL               SETRIS
         LOGICAL               STRIS
 
   C
   C     Local parameters
   C
         INTEGER               NAMLEN
         PARAMETER           ( NAMLEN = 32 )
 
         INTEGER               CORLEN
         PARAMETER           ( CORLEN = 10 )
 
   C
   C     Local variables
   C
         CHARACTER*(CORLEN)    SVCORR
         CHARACTER*(NAMLEN)    SVSFPT
         CHARACTER*(NAMLEN)    SVTARG
 
         DOUBLE PRECISION      AZ
         DOUBLE PRECISION      EL
         DOUBLE PRECISION      R
         DOUBLE PRECISION      SVLIM
 
   C
   C     Saved variables
   C
         SAVE                  SVCORR
         SAVE                  SVLIM
         SAVE                  SVSFPT
         SAVE                  SVTARG
 
 
   C
   C     A return value must be supplied because this is
   C     an entry point of a function.  The value has no
   C     significance.
   C
         MODRIS = .FALSE.
 
         CALL SIGERR ( 'SPICE(BOGUSENTRY)' )
         RETURN
 
 
 
 
         ENTRY SETRIS ( SRFPT, ABCORR, TARGET, LIMIT )
   C
   C     Save function arguments to be used while searching.
   C
         SVCORR = ABCORR
         SVSFPT = SRFPT
         SVTARG = TARGET
         SVLIM  = LIMIT
   C
   C     A return value must be supplied because this is
   C     an entry point of a function.  The value has no
   C     significance.
   C
         SETRIS = .FALSE.
 
         RETURN
 
 
 
 
         ENTRY STRIS ( ET )
   C
   C     Logical function to indicate "state" of target:  the state
   C     is .TRUE. when the target is above the elevation limit
   C     in the topocentric frame at the saved surface point;
   C     otherwise the state is .FALSE.
   C
         CALL SRFAZL ( SVSFPT, ET, SVCORR, SVTARG, AZ, EL, R )
 
         STRIS = EL .GT. SVLIM
 
         RETURN
 
         END
 


Point Target Occultation State Function Module MODOCP




 
         LOGICAL FUNCTION MODOCP ( TARGET, OCCBDY, ET,
        .                          ABCORR, OBSRVR     )
         IMPLICIT NONE
   C
   C     Umbrella routine for point target occultation function.
   C
         CHARACTER*(*)         TARGET
         CHARACTER*(*)         OCCBDY
         DOUBLE PRECISION      ET
         CHARACTER*(*)         ABCORR
         CHARACTER*(*)         OBSRVR
 
   C
   C     Non-SPICELIB functions
   C
         LOGICAL               OCCPT
 
   C
   C     Entry points
   C
         LOGICAL               SETOCP
         LOGICAL               STOCP
 
   C
   C     Local parameters
   C
         INTEGER               NAMLEN
         PARAMETER           ( NAMLEN = 32 )
 
         INTEGER               CORLEN
         PARAMETER           ( CORLEN = 10 )
 
   C
   C     Local variables
   C
         CHARACTER*(CORLEN)    SVCORR
         CHARACTER*(NAMLEN)    SVOBSV
         CHARACTER*(NAMLEN)    SVOCBD
         CHARACTER*(NAMLEN)    SVTARG
 
   C
   C     Saved variables
   C
         SAVE                  SVCORR
         SAVE                  SVOBSV
         SAVE                  SVOCBD
         SAVE                  SVTARG
 
 
   C
   C     A return value must be supplied because this is
   C     an entry point of a function.  The value has no
   C     significance.
   C
         MODOCP = .FALSE.
 
         CALL SIGERR ( 'SPICE(BOGUSENTRY)' )
         RETURN
 
 
 
 
         ENTRY SETOCP ( TARGET, OCCBDY, ABCORR, OBSRVR )
   C
   C     Save function arguments to be used while searching.
   C
         SVCORR = ABCORR
         SVOCBD = OCCBDY
         SVTARG = TARGET
         SVOBSV = OBSRVR
 
   C
   C     A return value must be supplied because this is
   C     an entry point of a function.  The value has no
   C     significance.
   C
         SETOCP = .FALSE.
 
         RETURN
 
 
 
 
         ENTRY STOCP ( ET )
   C
   C     Logical function to indicate "state" of target:  the state
   C     is .TRUE. when the target is occulted as seen from the
   C     observer's location; otherwise the state is .FALSE.
   C
         STOCP = OCCPT ( SVTARG, SVOCBD, ET, SVCORR, SVOBSV )
 
         RETURN
 
         END
 
 


Event finder FNDEVT




 
         SUBROUTINE FNDEVT ( ETBEG, ETEND, STEP, STATEF, XSET )
         IMPLICIT NONE
   C
   C     Find locations of state changes for a boolean function.
   C
   C       12-SEP-2002 (NJB)
   C       06-JUN-2002 (NJB)
   C
         INTEGER               LBCELL
         PARAMETER           ( LBCELL = -5 )
 
         DOUBLE PRECISION      ETBEG
         DOUBLE PRECISION      ETEND
         DOUBLE PRECISION      STEP
         LOGICAL               STATEF
         DOUBLE PRECISION      XSET  ( LBCELL : * )
 
         EXTERNAL              STATEF
 
 
   C
   C     SPICELIB functions
   C
         DOUBLE PRECISION      BRCKTD
 
   C
   C     Local parameters
   C
         DOUBLE PRECISION      CNVLIM
         PARAMETER           ( CNVLIM = 1.D-6 )
 
         INTEGER               MAXITR
         PARAMETER           ( MAXITR = 100 )
 
   C
   C     Local variables
   C
         DOUBLE PRECISION      DELTA
         DOUBLE PRECISION      LOWER
         DOUBLE PRECISION      LPT
         DOUBLE PRECISION      MIDPT
         DOUBLE PRECISION      UPPER
         DOUBLE PRECISION      UPT
 
         INTEGER               NITR
 
         LOGICAL               LSTATE
         LOGICAL               MSTATE
         LOGICAL               USTATE
 
 
         CALL CHKIN ( 'FNDEVT' )
   C
   C     Set the cardinality of the input cell to zero.
   C
         CALL SCARDD  ( 0, XSET )
 
   C
   C     If the input interval is empty, we're done.
   C
         IF ( ETEND .LE. ETBEG ) THEN
            CALL CHKOUT ( 'FNDEVT' )
            RETURN
         END IF
 
   C
   C     The step size must be at least as large as the convergence
   C     limit.
   C
         IF ( STEP .LE. CNVLIM ) THEN
 
            CALL SETMSG ( 'STEP must be > the convergence limit #; '//
        .                 'actual value was #.'                      )
            CALL ERRDP  ( '#',  CNVLIM                               )
            CALL ERRDP  ( '#',  STEP                                 )
            CALL SIGERR ( 'SPICE(INVALIDVALUE)'                      )
            CALL CHKOUT ( 'FNDEVT'                                   )
            RETURN
 
         END IF
 
   C
   C     Obtain the initial state; save this as the "lower" state.
   C
         LSTATE = STATEF ( ETBEG )
 
   C
   C     Step along the input interval, looking for state changes.
   C
         LPT    =  ETBEG
         UPT    =  BRCKTD ( ETBEG + STEP,  ETBEG,  ETEND )
 
         DO WHILE ( LPT .LT. ETEND )
   C
   C        Find the state at the epoch UPT.  First, check that UPT
   C        is actually greater than LPT.
   C
            IF ( UPT .LE. LPT ) THEN
   C
   C           We're not getting anywhere; the step is too small.
   C
               CALL SETMSG ( 'Upper bracketing epoch UPT = #; '     //
        .                    'lower epoch LPT = #. This condition ' //
        .                    'arises when the step size is too '    //
        .                    'small.'                               )
               CALL ERRDP  ( '#', UPT                               )
               CALL ERRDP  ( '#', LPT                               )
               CALL SIGERR ( 'SPICE(STEPTOOSMALL)'                  )
               CALL CHKOUT ( 'FNDEVT'                               )
               RETURN
 
            END IF
 
   C
   C        Find the state at the upper bound of our step interval.
   C
            USTATE = STATEF ( UPT )
 
 
            IF ( USTATE .NEQV. LSTATE ) THEN
   C
   C           There's a state change between the right endpoint UPT
   C           and the left endpoint LPT.  Do a binary search to
   C           locate the epoch at which the state changes.  Note:
   C           uniqueness of the root is not guaranteed; it's up
   C           to the caller to choose STEP small enough to ensure
   C           that only one root can occur in a time interval of
   C           length STEP.
   C
               LOWER  =  BRCKTD ( LPT, ETBEG, ETEND )
               UPPER  =  BRCKTD ( UPT, ETBEG, ETEND )
               MIDPT  =  ( LOWER + UPPER ) / 2
               MSTATE =  STATEF ( MIDPT )
 
               NITR   =  0
               DELTA  =  UPPER - LOWER
 
               DO WHILE (         ( DELTA .GT. CNVLIM )
        .                  .AND.  ( NITR  .LE. MAXITR )  )
   C
   C              Adjust our search interval so the length is reduced
   C              by a factor of two and so that the location of the
   C              state change remains between LOWER and UPPER.
   C
                  IF ( MSTATE .EQV. USTATE ) THEN
   C
   C                 The state is the same at the midpoint and the
   C                 upper bound, so a state change must occur
   C                 between the lower bound and the midpoint.
   C                 Make the midpoint the new upper bound.
   C
                     UPPER = MIDPT
 
                  ELSE
   C
   C                 There is a state change between the midpoint
   C                 and the upper bound; the state at the midpoint
   C                 matches the state at the lower bound.  Make
   C                 the midpoint the new lower bound.
   C
                     LOWER  =  MIDPT
 
                  END IF
 
                  MIDPT  = ( LOWER + UPPER ) / 2
                  MSTATE =   STATEF ( MIDPT )
 
                  DELTA  =   UPPER - LOWER
                  NITR   =   NITR  + 1
 
   C
   C              At this point, the state at LOWER matches the state
   C              at LPT (LSTATE), and the state at UPPER matches the
   C              state at UPT (USTATE).
   C
               END DO
 
   C
   C           If we dropped out of the loop because we hit the
   C           iteration limit, we have a problem.
   C
               IF ( DELTA .GT. CNVLIM ) THEN
 
                  CALL SETMSG ( 'Binary search failed to converge. '//
        .                       'ETBEG = #, LOWER = #, MIDPT = #, ' //
        .                       'UPPER = #, ETEND = #'               )
                  CALL ERRDP  ( '#',  ETBEG                          )
                  CALL ERRDP  ( '#',  LOWER                          )
                  CALL ERRDP  ( '#',  MIDPT                          )
                  CALL ERRDP  ( '#',  UPPER                          )
                  CALL ERRDP  ( '#',  ETEND                          )
                  CALL SIGERR ( 'SPICE(NOCONVERGENCE)'               )
                  CALL CHKOUT ( 'FNDEVT'                             )
                  RETURN
 
               END IF
 
   C
   C           The epoch of the state transition has been determined
   C           to within CNVLIM.  We'll use UPPER as the epoch of
   C           the transition.  This ensures that the state at the
   C           transition epoch is USTATE.
   C
               CALL APPNDD ( UPPER, XSET )
 
   C
   C           UPPER becomes the left endpoint of the next interval.
   C           The state at the right interval endpoint UPT becomes
   C           the state at the left endpoint of the next interval.
   C           Calculate UPT at the next step.  We look up
   C           USTATE at the top of the loop.
   C
               LPT     =  UPPER
               LSTATE  =  USTATE
               UPT     =  BRCKTD ( LPT + STEP,  ETBEG,  ETEND )
 
            ELSE
   C
   C           No state change was found on this step.  UPT becomes
   C           the lower bound of the next search interval.  LSTATE
   C           remains unchanged.  We look up USTATE at the top of
   C           the loop.
   C
               LPT     =  UPT
               UPT     =  BRCKTD ( LPT + STEP,  ETBEG,  ETEND )
 
            END IF
 
 
         END DO
   C
   C     The left endpoint equals the right endpoint of our search
   C     interval, so there are no more state changes to be found.
   C
 
         CALL CHKOUT ( 'FNDEVT' )
         RETURN
         END


Set to Window Mapping Utility MAKWIN




 
         SUBROUTINE MAKWIN ( ETBEG, ETEND, FSTATE, EVTSET, EVTWIN )
         IMPLICIT NONE
   C
   C     Utility program for creating a window representing
   C     time intervals when a specified binary state function
   C     returns .TRUE.
   C
   C
   C     Subroutine arguments:
   C
   C        ETBEG         is the start time of a time interval on
   C                      which the binary state function
   C                      FSTATE is defined.
   C
   C        ETEND         is the end time corresponding to ETBEG.
   C
   C        FSTATE        is a binary, state function of ET that
   C                      returns a value of type LOGICAL.
   C
   C        EVTSET        is a SPICELIB set containing epochs of
   C                      state transitions.  The members of EVTSET
   C                      must lie in the interval
   C
   C                         [ETBEG, ETEND]
   C
   C        EVTWIN        is a SPICELIB window constructed by this
   C                      routine.  EVTWIN contains subintervals
   C                      of the interval
   C
   C                         [ETBEG, ETEND]
   C
   C                      on which FSTATE(ET) is .TRUE.  EVTWIN
   C                      must be initialized by the caller of this
   C                      routine.
   C
   C
 
         INTEGER               LBCELL
         PARAMETER           ( LBCELL = -5 )
 
         DOUBLE PRECISION      ETBEG
         DOUBLE PRECISION      ETEND
         LOGICAL               FSTATE
         DOUBLE PRECISION      EVTSET ( LBCELL : * )
         DOUBLE PRECISION      EVTWIN ( LBCELL : * )
 
   C
   C
   C     SPICELIB functions
   C
         INTEGER               CARDD
 
   C
   C     External functions
   C
         EXTERNAL              FSTATE
 
   C
   C     Local variables
   C
         INTEGER               CARD
         INTEGER               I
 
 
         CALL CHKIN ( 'MAKWIN' )
 
   C
   C     Empty the output window.
   C
         CALL SCARDD ( 0, EVTWIN )
 
   C
   C     Get a local copy of the cardinality of the event set.
   C
         CARD = CARDD(EVTSET)
 
   C
   C     If FSTATE is .TRUE. at ETBEG, then the first state transition
   C     is the right endpoint of the first window interval, and every
   C     subsequent even-numbered transition marks the start of a
   C     window interval. Otherwise, every odd-numbered transition
   C     marks the start of a window interval.
   C
         IF ( FSTATE(ETBEG) ) THEN
   C
   C        The first window interval starts at ETBEG.
   C
            IF ( CARD .EQ. 0 ) THEN
   C
   C           The function is .TRUE. the whole time.  The
   C           window consists of the single interval
   C
   C              [ETBEG, ETEND]
   C
               CALL WNINSD ( ETBEG, ETEND, EVTWIN )
 
            ELSE
   C
   C           Insert the first interval into the window.
   C
               CALL WNINSD ( ETBEG, EVTSET(1), EVTWIN )
 
   C
   C           Insert into our window any remaining intervals bounded
   C           by transitions.
   C
               I = 2
 
               DO WHILE ( I+1  .LE. CARD )
 
                  CALL WNINSD ( EVTSET(I), EVTSET(I+1), EVTWIN )
 
                  I = I + 2
 
               END DO
 
               IF ( I .EQ. CARD ) THEN
   C
   C              The last interval begins at EVTSET(CARD)
   C              and ends at ETEND.
   C
                  CALL WNINSD ( EVTSET(CARD), ETEND, EVTWIN )
 
               END IF
 
            END IF
   C
   C        We've handled the case where FSTATE is .TRUE. at ETBEG.
   C
         ELSE
   C
   C        The first interval starts at EVTSET(1).
   C
            I = 1
 
            DO WHILE ( I+1  .LE. CARD )
 
               CALL WNINSD ( EVTSET(I), EVTSET(I+1), EVTWIN )
 
               I = I + 2
 
            END DO
 
            IF ( I .EQ. CARD ) THEN
   C
   C           The last interval begins at EVTSET(CARD)
   C           and ends at ETEND.
   C
               CALL WNINSD ( EVTSET(CARD), ETEND, EVTWIN )
 
            END IF
   C
   C        We've handled the case where FSTATE is .FALSE. at ETBEG.
   C
         END IF
 
   C
   C     EVTWIN is ready for use.
   C
         CALL CHKOUT ( 'MAKWIN' )
         END