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
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 TDBwhen 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 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.
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 framesThese 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
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 typeAnother 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.
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.forin the FORTRAN Toolkit and in
cspice/src/cspice/str2et_c.cin 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.
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 PCKThese SPICE kernels are available from the NAIF server at JPL:
ftp://naif.jpl.nasa.gov/pub/naif/misc/lessons/events
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. windowThe most detailed documentation source for these routines is their headers.
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 limitThe 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.
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 student is introduced to flexible, robust techniques for
introducing data into SPICELIB based programs.
The student is exposed to the following SPICE components:
A possible solution would consist of the following steps:
Preparation:
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
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
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
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 RThe 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.
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.
A possible solution would consist of the following steps:
Write the SRFAZL routine:
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
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
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
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.
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.
A possible solution would consist of the following steps:
Write the OCCPT routine:
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
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
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.
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.
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.
A possible solution would consist of the following steps:
CALL FNDEVT ( ETBEG, ETEND, STEPSZ, STRIS, EVTSET )
CALL MAKWIN ( ETBEG, ETEND, STRIS, EVTSET, RISWIN )
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
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)
Extend the program of the previous chapter to find times when the
target object is:
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.
The student gains further experience with the SPICELIB cell, set and
window routines, and with root finding techniques.
A possible solution would consist of the following steps:
CALL FNDEVT ( ETBEG, ETEND, STEPSZ, STOCP, EVTSET )
CALL MAKWIN ( ETBEG, ETEND, STOCP, EVTSET, OCCWIN )
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
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)
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
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
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
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