[Spice_discussion] Solving for constraints on L_s: why not a loop?
Nat Bachman
nathaniel.bachman at jpl.nasa.gov
Mon Apr 2 01:37:04 PDT 2012
Hi Michael,
A bit more about the L_s constraint topic:
The GF solution is complicated.
Why not solve the problem by simply calling
LSPCN in a loop?
In fact calling LSPCN in a loop may be a perfect
solution for this problem and others.
The GF approach tends to be worthwhile when at least
some of the following requirements or conditions apply:
- The solution (roots) must be accurate.
There are (at least) three ways the GF subsystem enhances
accuracy of solutions:
- Where applicable, the GF subsystem finds
local maxima and minima using zero-crossings
of derivatives.
- The GF subsystem provides accurate algorithms
for geometric quantities. Sometimes these algorithms
require more development work than is practical for
individual SPICE users to expend on their own.
- The GF subsystem uses tolerance values independent
of the step sizes used for root bracketing. For example,
the LSPCN search was performed with a 30-day step size
for root bracketing and the default 1 microsecond
convergence tolerance.
- There are multiple roots.
- Bracketing the roots is problematic due to lack
of knowledge of even approximate locations of the
roots.
- The solution process must run quickly.
For example, the GF Required Reading contains an example
program "CASCADE" which demonstrates speeding up an
occultation search by a factor of 100.
- Multiple constraints are involved. For example,
one may wish to find times when a Mars orbiter
is not occulted by Mars and is above a certain
elevation limit as seen from a earth-based tracking
station.
- The time window over which the search is conducted
is not a single time interval. This situation can
arise when there are multiple constraints, or when
a time window representing data availability consists
of multiple, disjoint intervals.
As an example of a slightly non-obvious search performance
improvement provided by the GF subsystem, I extended the
previous example to search for times when L_s is greater
than 359 degrees. The solution intervals are roughly one
day long, but they're found using a search step of 30 days.
An explanation of why this works can be found in the backup
pages of the GF tutorial.
The source code is attached as ls_ex2.txt, since security
software doesn't seem to like the .for suffix.
Best regards,
-Nat
Nathaniel.Bachman at jpl.nasa.gov
-------------- next part --------------
Search constraint:
Coordinate system: LATITUDINAL
Coordinate: LONGITUDE
Relation: =
Reference value: 90.000000000000000
Number of result intervals found: 2
Result window
1 2012 JUN 20 23:02:07.083569 TDB 2012 JUN 20 23:02:07.083569 TDB
LS (deg) = 90.000000000000256
2 2013 JUN 21 05:09:00.484632 TDB 2013 JUN 21 05:09:00.484632 TDB
LS (deg) = 90.000000000000426
Search constraint:
Coordinate system: RA/DEC
Coordinate: RIGHT ASCENSION
Relation: >
Reference value: 359.00000000000000
Number of result intervals found: 2
Result window
1 2012 MAR 19 05:04:41.339143 TDB 2012 MAR 20 05:13:09.185274 TDB
LS @start (deg) = 359.00000000000063
LS @stop (deg) = 359.99999999942656
2 2013 MAR 19 10:51:06.842982 TDB 2013 MAR 20 11:00:09.046245 TDB
LS @start (deg) = 359.00000000000000
LS @stop (deg) = 359.99999999942719
-------------- next part --------------
C
C Program: LS_EX2
C Author: Nat Bachman (NAIF/JPL)
C Date: 31-MAR-2012
C
C Example program created for use by Michael Aye.
C
C This program searches for times when specified
C constraints on L_s are met.
C
C This program extends the example program LS_EX1
C by running a second search: the program finds
C the time interval when L_s is greater than 359
C degrees. Although the time intervals over which
C the condition is met are about one day long,
C the search is successfully carried out using a
C bracketing step size of 30 days.
C
PROGRAM LS_EX2
IMPLICIT NONE
C
C SPICELIB functions
C
DOUBLE PRECISION DPR
DOUBLE PRECISION LSPCN
DOUBLE PRECISION SPD
DOUBLE PRECISION RPD
INTEGER WNCARD
C
C Global parameters
C
INCLUDE 'gf.inc'
C
C Local parameters
C
C
C Meta-kernel:
C
CHARACTER*(*) META
PARAMETER ( META = 'lspcn.tm' )
C
C Output time string format:
C
CHARACTER*(*) TIMFMT
PARAMETER ( TIMFMT =
. 'YYYY MON DD HR:MN:SC.###### TDB::TDB' )
C
C SPICELIB cell lower bound:
C
INTEGER LBCELL
PARAMETER ( LBCELL = -5 )
C
C Window size parameters:
C
INTEGER MAXIVL
PARAMETER ( MAXIVL = 10000 )
INTEGER MAXWIN
PARAMETER ( MAXWIN = 2 * MAXIVL )
C
C Length of time strings:
C
INTEGER TIMLEN
PARAMETER ( TIMLEN = 40 )
INTEGER BDNMLN
PARAMETER ( BDNMLN = 36 )
INTEGER FRNMLN
PARAMETER ( FRNMLN = 32 )
INTEGER SYSLEN
PARAMETER ( SYSLEN = 50 )
INTEGER CRDLEN
PARAMETER ( CRDLEN = 50 )
INTEGER CORLEN
PARAMETER ( CORLEN = 10 )
INTEGER OPLEN
PARAMETER ( OPLEN = 50 )
C
C Local variables
C
CHARACTER*(CORLEN) ABCORR
CHARACTER*(TIMLEN) BEGSTR
CHARACTER*(CRDLEN) COORD
CHARACTER*(SYSLEN) CRDSYS
CHARACTER*(TIMLEN) ENDSTR
CHARACTER*(FRNMLN) FRAME
CHARACTER*(BDNMLN) OBSRVR
CHARACTER*(OPLEN) RELATE
CHARACTER*(BDNMLN) TARGET
CHARACTER*(TIMLEN) TIMSTR
C
C Confinement and result windows:
C
DOUBLE PRECISION CNFINE ( LBCELL : MAXWIN )
DOUBLE PRECISION RESULT ( LBCELL : MAXWIN )
C
C Workspace:
C
DOUBLE PRECISION WORK ( LBCELL : MAXWIN, NWMAX )
C
C Other local variables
C
DOUBLE PRECISION ADJUST
DOUBLE PRECISION ET
DOUBLE PRECISION ET0
DOUBLE PRECISION ET1
DOUBLE PRECISION REFVAL
DOUBLE PRECISION START
DOUBLE PRECISION STEP
DOUBLE PRECISION FINISH
DOUBLE PRECISION LS
INTEGER I
INTEGER N
C
C Executable code follows.
C
C
C Initialize SPICE windows. Note that the workspace
C need not be initialized; the GF routines will take
C care of that.
C
C Although our confinement window will contain at
C most a few time intervals, we'll leave room for
C a large number of them.
C
CALL SSIZED ( MAXWIN, CNFINE )
CALL SSIZED ( MAXWIN, RESULT )
C
C Load kernels.
C
CALL FURNSH ( META )
C
C Create the "confinement window": the time window
C over which we'll conduct the search.
C
C We can insert as many as MAXIVL disjoint time intervals
C into CNFINE. In this case we'll just insert one.
C
TIMSTR = '2012 MAR 1'
CALL STR2ET ( TIMSTR, ET0 )
TIMSTR = '2013 OCT 1'
CALL STR2ET ( TIMSTR, ET1 )
CALL WNINSD ( ET0, ET1, CNFINE )
C
C We'll search for times within the confinement window
C when L_s (relative to the earth) is 90 degrees.
C
C Set up our search parameters.
C
FRAME = 'EARTH_LSPCN'
TARGET = 'SUN'
ABCORR = 'NONE'
OBSRVR = 'EARTH'
CRDSYS = 'LATITUDINAL'
COORD = 'LONGITUDE'
RELATE = '='
REFVAL = 90.D0 * RPD()
ADJUST = 0.D0
C
C Use a 30-day time step for root bracketing.
C
STEP = 30 * SPD()
C
C Search for events.
C
CALL GFPOSC ( TARGET, FRAME, ABCORR, OBSRVR,
. CRDSYS, COORD, RELATE, REFVAL,
. ADJUST, STEP, CNFINE, MAXWIN,
. NWMAX, WORK, RESULT )
N = WNCARD( RESULT )
WRITE (*,*) ' '
WRITE (*,*) 'Search constraint:'
WRITE (*,*) ' '
WRITE (*,*) 'Coordinate system: ', CRDSYS
WRITE (*,*) 'Coordinate: ', COORD
WRITE (*,*) 'Relation: ', RELATE
WRITE (*,*) 'Reference value: ', REFVAL * DPR()
WRITE (*,*) ' '
WRITE (*,*) 'Number of result intervals found: ', N
IF ( N .GT. 0 ) THEN
WRITE (*,*) ' '
WRITE (*,*) 'Result window'
WRITE (*,*) ' '
DO I = 1, N
C
C Fetch the Ith interval from the result window,
C convert the time bounds to strings, and display
C the interval.
C
CALL WNFETD ( RESULT, I, START, FINISH )
CALL TIMOUT ( START, TIMFMT, BEGSTR )
CALL TIMOUT ( FINISH, TIMFMT, ENDSTR )
WRITE (*, '(1X,I5,3X,2A35)' ) I, BEGSTR, ENDSTR
C
C Check results: find L_s at the interval start time.
C
LS = LSPCN ( OBSRVR, START, ABCORR )
WRITE (*,*) ' LS (deg) = ', LS * DPR()
WRITE (*,*) ' '
END DO
END IF
C
C Now for the second search.
C
RELATE = '>'
REFVAL = 359.D0 * RPD()
C
C For compatibility with LSPCN, which outputs an angle
C in the range 0 : 2*pi, use RA/DEC and right ascension
C as the coordinate system and coordinate.
C
CRDSYS = 'RA/DEC'
COORD = 'RIGHT ASCENSION'
C
C Search for events.
C
CALL GFPOSC ( TARGET, FRAME, ABCORR, OBSRVR,
. CRDSYS, COORD, RELATE, REFVAL,
. ADJUST, STEP, CNFINE, MAXWIN,
. NWMAX, WORK, RESULT )
N = WNCARD( RESULT )
WRITE (*,*) ' '
WRITE (*,*) ' '
WRITE (*,*) 'Search constraint:'
WRITE (*,*) ' '
WRITE (*,*) 'Coordinate system: ', CRDSYS
WRITE (*,*) 'Coordinate: ', COORD
WRITE (*,*) 'Relation: ', RELATE
WRITE (*,*) 'Reference value: ', REFVAL * DPR()
WRITE (*,*) ' '
WRITE (*,*) ' '
WRITE (*,*) 'Number of result intervals found: ', N
IF ( N .GT. 0 ) THEN
WRITE (*,*) ' '
WRITE (*,*) 'Result window'
WRITE (*,*) ' '
DO I = 1, N
C
C Fetch the Ith interval from the result window,
C convert the time bounds to strings, and display
C the interval.
C
CALL WNFETD ( RESULT, I, START, FINISH )
CALL TIMOUT ( START, TIMFMT, BEGSTR )
CALL TIMOUT ( FINISH, TIMFMT, ENDSTR )
WRITE (*, '(1X,I5,3X,2A35)' ) I, BEGSTR, ENDSTR
C
C Check results: find L_s at the interval boundary
C times.
C
LS = LSPCN ( OBSRVR, START, ABCORR )
WRITE (*,*) ' LS @start (deg) = ', LS * DPR()
LS = LSPCN ( OBSRVR, FINISH, ABCORR )
WRITE (*,*) ' LS @stop (deg) = ', LS * DPR()
WRITE (*,*) ' '
END DO
END IF
WRITE (*,*) ' '
END
-------------- next part --------------
KPL/FK
File: lspcn.tf
Author: Nat Bachman (JPL/NAIF)
Date: 30-MAR-2012
Example dynamic frame kernel created by Nat Bachman for
Michael Aye.
This kernel can be used together with the SPICE GF
(geometry finder) subsystem to determine times when
solar longitude L_s with respect to the earth
satisfies given constraints.
This kernel can be edited to define new reference frames
for computing L_s with respect to bodies other than the
earth.
Frame Specifications
====================
GSE: Geocentric Solar Ecliptic
This is a reference frame based on the instantaneous orbital
angular velocity vector of the earth. Note that it's also
possible to define this frame using the mean ecliptic of
date as the X-Y plane.
The +X axis points from the earth to the sun.
The +Z axis is aligned with the cross product
of the +X axis and velocity of the sun with respect
to the earth.
The +Y axis is the cross product +Z x +X, completing
the right-handed frame.
\begindata
FRAME_GSE = 1400000
FRAME_1400000_NAME = 'GSE'
FRAME_1400000_CLASS = 5
FRAME_1400000_CLASS_ID = 1400000
FRAME_1400000_CENTER = 399
FRAME_1400000_RELATIVE = 'J2000'
FRAME_1400000_DEF_STYLE = 'PARAMETERIZED'
FRAME_1400000_FAMILY = 'TWO-VECTOR'
FRAME_1400000_PRI_AXIS = 'X'
FRAME_1400000_PRI_VECTOR_DEF = 'OBSERVER_TARGET_POSITION'
FRAME_1400000_PRI_OBSERVER = 'EARTH'
FRAME_1400000_PRI_TARGET = 'SUN'
FRAME_1400000_PRI_ABCORR = 'NONE'
FRAME_1400000_SEC_AXIS = 'Y'
FRAME_1400000_SEC_VECTOR_DEF = 'OBSERVER_TARGET_VELOCITY'
FRAME_1400000_SEC_OBSERVER = 'EARTH'
FRAME_1400000_SEC_TARGET = 'SUN'
FRAME_1400000_SEC_ABCORR = 'NONE'
FRAME_1400000_SEC_FRAME = 'J2000'
\begintext
The EARTH_LSPCN reference frame is the one in which solar
planetocentric longitude L_s is actually measured.
This frame is defined as specified in the header of
the SPICELIB routine EARTH_LSPCN:
This is the longitude of the body-sun vector in a
right-handed frame whose basis vectors are defined
as follows:
- The positive Z direction is given by the instantaneous
angular velocity vector of the orbit of the body about
the sun.
- The positive X direction is that of the cross product
of the instantaneous north spin axis of the body
with the positive Z direction.
- The positive Y direction is Z x X.
We can re-cast the specification of the X and Y axes
in a form usable by the dynamic frames subsystem:
- The positive Z direction is as above.
- The positive Y direction is aligned with the
component orthogonal to +Z of the north spin
axis.
- The positive X direction is Y x Z.
For simplicity, the earth rotational elements determining
the north pole are taken from the IAU rotation model.
Changing the reference frame for the pole from IAU_EARTH
to ITRF93 will yield a more accurate instantaneous pole
direction, if desired. The application using this kernel
will need to load a binary earth PCK file in order to use
the ITRF93 frame.
\begindata
FRAME_EARTH_LSPCN = 1400001
FRAME_1400001_NAME = 'EARTH_LSPCN'
FRAME_1400001_CLASS = 5
FRAME_1400001_CLASS_ID = 1400001
FRAME_1400001_CENTER = 399
FRAME_1400001_RELATIVE = 'J2000'
FRAME_1400001_DEF_STYLE = 'PARAMETERIZED'
FRAME_1400001_FAMILY = 'TWO-VECTOR'
FRAME_1400001_PRI_AXIS = 'Z'
FRAME_1400001_PRI_VECTOR_DEF = 'CONSTANT'
FRAME_1400001_PRI_FRAME = 'GSE'
FRAME_1400001_PRI_SPEC = 'RECTANGULAR'
FRAME_1400001_PRI_VECTOR = ( 0, 0, 1 )
FRAME_1400001_PRI_ABCORR = 'NONE'
FRAME_1400001_SEC_AXIS = 'Y'
FRAME_1400001_SEC_VECTOR_DEF = 'CONSTANT'
FRAME_1400001_SEC_FRAME = 'IAU_EARTH'
FRAME_1400001_SEC_SPEC = 'RECTANGULAR'
FRAME_1400001_SEC_VECTOR = ( 0, 0, 1 )
FRAME_1400001_SEC_ABCORR = 'NONE'
\begintext
[End of kernel]
-------------- next part --------------
KPL/MK
File: lspcn.tm
Author: Nat Bachman (JPL/NAIF)
Date: 30-MAR-2012
Example meta-kernel created by Nat Bachman for
Michael Aye.
This kernel can be used together with the SPICE GF
(geometry finder) subsystem to determine times when
solar longitude L_s with respect to the earth
satisfies given constraints.
The kernels shown here, other than the frame
kernel, are available from the NAIF server.
Kernels loaded by this meta-kernel:
naif0010.tls Leapseconds
pck00010.tpc Generic PCK
(provides earth orientation)
de421.bsp JPL planetary ephemeris
lspcn.tf Frame kernel specifying the
reference frame in which
L_s is measured
\begindata
KERNELS_TO_LOAD = ( 'naif0010.tls'
'pck00010.tpc'
'de421.bsp'
'lspcn.tf' )
\begintext
[End of kernel]
More information about the Spice_discussion
mailing list