In-situ Sensing Hands-On Lesson (FORTRAN) |
Table of ContentsIn-situ Sensing Hands-On Lesson (FORTRAN) Overview Note About HTML Links References Tutorials Required Reading Documents Source Code Header Comments Kernels Used SPICE Routines Used Step-1: ``UTC to ET'' ``UTC to ET'' Task Statement ``UTC to ET'' Hints ``UTC to ET'' Solution Steps ``UTC to ET'' Code Step-2: ``SCLK to ET'' ``SCLK to ET'' Task Statement ``SCLK to ET'' Hints ``SCLK to ET'' Solution Steps ``SCLK to ET'' Code Step-3: ``Spacecraft State'' ``Spacecraft State'' Task Statement ``Spacecraft State'' Hints ``Spacecraft State'' Solution Steps ``Spacecraft State'' Code Step-4: ``Sun Direction'' ``Sun Direction'' Task Statement ``Sun Direction'' Hints ``Sun Direction'' Solution Steps ``Sun Direction'' Code Step-5: ``Sub-Spacecraft Point'' ``Sub-Spacecraft Point'' Task Statement ``Sub-Spacecraft Point'' Hints ``Sub-Spacecraft Point'' Solution Steps ``Sub-Spacecraft Point'' Code Step-6: ``Spacecraft Velocity'' ``Spacecraft Velocity'' Task Statement ``Spacecraft Velocity'' Hints ``Spacecraft Velocity'' Solution Steps ``Spacecraft Velocity'' Code Program ``spice_example.f'': In-situ Sensing Hands-On Lesson (FORTRAN)
Overview
Note About HTML Links
In order for the links to be resolved, create a subdirectory called ``lessons'' under the ``doc/html'' directory of the Toolkit tree and copy this document to that subdirectory before loading it into a Web browser. References
Of these documents, the ``Tutorials'' contains the highest level descriptions with the least number of details while the ``Required Reading'' documents contain much more detailed specifications. The most complete specifications are provided in the ``Headers'' -- the comments in the top section of the source file. In some cases the lesson explanations also refer to the information provided in the meta-data area of the kernels used in the lesson examples. It is especially true in case of the FK and IK files, which often contain comprehensive descriptions of the frames, instrument FOVs, etc. Since both FK and IK are text kernels, the information provided in them can be viewed using any text editor, while the meta information provided in binary kernels -- SPKs and CKs -- can be viewed using ``commnt'' or ``spacit'' utility programs located in ``toolkit/exe'' of Toolkit installation tree. Tutorials
Name Lesson steps/routines that it describes --------------- ----------------------------------------- Time UTC to ET and SCLK to ET Loading Kernels Loading SPICE kernels SCLK SCLK to ET time conversion SPK Computing positions and velocities Frames Computing transformations between framesThese tutorials are available in printed form and as MS Office or PDF files from NAIF server at JPL:
http://naif.jpl.nasa.gov/naif/tutorials.html Required Reading Documents
Name Lesson steps/routines that it describes --------------- ----------------------------------------- time.req UTC to ET time conversion kernel.req Loading SPICE kernels sclk.req SCLK to ET time conversion naif_ids.req Body and reference frame names spk.req Computing positions and velocitiesAnother very useful document, also distributed with the Toolkit, is ``Permuted Index'', called ``spicelib.idx'' for FORTRAN or ``cspice.idx'' for C, IDL, and MATLAB also found in the ``doc'' directory. This text document provides an easy way to find what SPICE routine(s) performs a particular function of interest and the name of the source file that contains this function (this is especially useful for FORTRAN because some of the routines are entry points and, therefore, their name is different from the name of the source file in which they are located.) Source Code Header Comments
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. Kernels Used
File Name Type Description ------------------------- ---- -------------------------- naif0008.tls LSK Generic LSK cpck05Mar2004.tpc PCK Cassini project PCK cas00084.tsc SCLK Cassini SCLK 020514_SE_SAT105.bsp SPK Saturnian Satellite Ephemeris SPK 030201AP_SK_SM546_T45.bsp SPK Cassini Spacecraft SPK 981005_PLTEPH-DE405S.bsp SPK Planetary Ephemeris SPK sat128.bsp SPK Saturnian Satellite Ephemeris SPK 04135_04171pc_psiv2.bc CK Cassini Spacecraft CK cas_v37.tf FK Cassini FKThese SPICE kernels are included in the lesson package available from the NAIF server at JPL:
ftp://naif.jpl.nasa.gov/pub/naif/toolkit_docs/Lessons/ SPICE Routines Used
Name Function that it performs ---------- --------------------------------------------------- FURNSH Loads kernels, individually or listed in meta-kernel STR2ET Converts UTC to ET SCS2E Converts SCLK to ET SPKEZR Computes states (position & velocity) SPKPOS Computes positions VHAT Find unit vector along a 3d vector VHATIP Find unit vector along a 3d vector, in place SUBPNT Computes body-fixed coordinates of sub-observer point RECLAT Converts rectangular coordinated to latitudinal VSUB Subtracts 3d vectors PXFORM Computes 3x3 matrix rotating vectors between frames MXV Multiplies 3d vector by 3x3 matrix VPACK Packs 3 number into a 3d vectorThe most detailed documentation source for these routines are their headers. Step-1: ``UTC to ET''``UTC to ET'' Task Statement
``UTC to ET'' Hints
Find necessary kernel(s) on the NAIF's FTP site. Find out what routine should be called to load necessary kernel(s). Reference the ``kernel.req'' and/or ``Loading Kernels'' tutorial. Find the ``loader'' routine calling sequence specification. Look at the ``time.req'' and that routine's source code header. This routine may be an entry point, in which case there will be no source file with the same name. To find out in which source file this entry point is, search for its name in the ``Permuted Index''. Find the routine(s) used to convert time between UTC and ET. Look at the ``time.req'' and/or ``Time'' tutorial. Find the ``converter'' routine(s) calling sequence specification. Look in the ``time.req'' and the routine's source code header. Put all calls together in a program, add variable declarations (the routine header's ``Declarations'' and ``Examples'' sections are a good place to look for declaration specification and examples) and output print statements. Compile it and link it against SPICELIB. ``UTC to ET'' Solution Steps
As any other SPICE kernel this file can be loaded by the FURNSH routine (entry in ``toolkit/src/spicelib/keeper.f''). For that, the name of the file can put be provided as a sole argument of this routine:
CHARACTER*(256) LSKFLE ... LSKFLE = 'naif0008.tls' ... CALL FURNSH( LSKFLE )or it can be listed in a meta-kernel:
\begindata KERNELS_TO_LOAD = ( 'kernels/lsk/naif0008.tls' ) \begintextthe name of which, let's call it ``spice_example.tm'', can be then provided as a sole argument of the FURNSH routine:
CHARACTER*(256) MKFILE ... MKFILE = 'spice_example.tm' ... CALL FURNSH( MKFILE )While the second option seems to involve a bit more work -- it requires making an extra file -- it is a much better way to go if you plan to load more kernels as you extend the program. With the meta-kernel approach simply adding more kernels to the list in KERNEL_TO_LOAD without changing the program code will accomplish that. The highest level SPICE time routine converting UTC to ET is STR2ET (``toolkit/src/spicelib/str2et.f''). It has two arguments -- input time string representing UTC in a variety of formats (see STR2ET header's section ``Particulars'' for the complete description of input time formats) and output DP number of ET seconds past J2000. A call to STR2ET converting a given UTC to ET could look like this:
CHARACTER*(32) UTC DOUBLE PRECISION ET ... UTC = '2004-06-11T19:32:00' CALL STR2ET( UTC, ET )By combining FURNSH and STR2ET calls and required declarations and by adding a simple print statement, one would get a complete program that prints ET for the given UTC epoch. The program's source code then needs to be compiled and linked against SPICELIB. Assuming that the program was saved in a file called ``spice_example.f'', this can be done with the following command on a Sun workstation:
> f77 -C -u -o spice_example spice_example.f spicelib.aThe command assumes that the library file(s) ``spicelib.a'' is located in current directory, which may not be the case. When you run the program's executable, ``spice_example'', it produces the following output (the output below was generated by this program compiled with g77 on a PC running Linux; your output may differ slightly in its format and numeric precision):
> ./spice_example UTC = 2004-06-11T19:32:00 ET = 140254384. > ``UTC to ET'' Code
IMPLICIT NONE CHARACTER*(256) MKFILE CHARACTER*(32) UTC DOUBLE PRECISION ET MKFILE = 'spice_example.tm' CALL FURNSH( MKFILE ) UTC = '2004-06-11T19:32:00' CALL STR2ET( UTC, ET ) WRITE(*,*) 'UTC = ', UTC WRITE(*,*) 'ET = ', ET ENDMeta-kernel file ``spice_example.tm'':
\begindata KERNELS_TO_LOAD = ( 'kernels/lsk/naif0008.tls' ) \begintext Step-2: ``SCLK to ET''``SCLK to ET'' Task Statement
``SCLK to ET'' Hints
Find necessary kernel(s) on the NAIF's FTP site. Modify the program or meta-kernel to load this(these) kernels. Find the routine(s) needed to convert time between SCLK and ET. Look at the ``sclk.req'' and/or ``Time'' and ``SCLK'' tutorials. Find the ``converter'' routine's calling sequence specification. Look in the ``sclk.req'' and the routine's source code header. Look at ``naif_ids.req'' and the comments in the additional kernel(s) that you have loaded for information on proper values of input arguments of this routine. Add calls to the ``converter'' routine(s), necessary variable declarations (the routine header's ``Declarations'' and ``Examples'' sections are a good place to look for declaration specification and examples), and output print statements to the program. Re-compile and re-link it against SPICELIB. ``SCLK to ET'' Solution Steps
No code change is needed in the loading portion of the program if a meta-kernel approach was used in the Step-1. The program will load the file if it will be added to the list of kernels in the KERNELS_TO_LOAD variable:
\begindata KERNELS_TO_LOAD = ( 'kernels/lsk/naif0008.tls' 'kernels/sclk/cas00084.tsc' ) \begintextThe highest level SPICE routine converting SCLK to ET is SCS2E (``toolkit/src/spicelib/scs2e.f''). It has three arguments -- NAIF ID for CASSINI s/c (-82 as described by ``naif_ids.req'' document), input time string representing CASSINI SCLK, and output DP number of ET seconds past J2000. A call to STR2ET converting given SCLK to ET could look like this:
CHARACTER*(32) SCLK INTEGER SCID ... SCID = -82 SCLK = '1465674964.105' CALL SCS2E( SCID, SCLK, ET )By adding the SCS2E call, required declarations and a simple print statement, one would get a complete program that prints ET for the given SCLK epoch. The program's source code then needs to be re-compiled and re-linked against SPICELIB. It can be done using the same compile command as in Step-1:
> f77 -C -u -o spice_example spice_example.f spicelib.aWhen you run the program's executable, ``spice_example'', it produces the following output (the output below was generated by this program compiled with g77 on a PC running Linux; your output may differ slightly in its format and numeric precision):
> ./spice_example UTC = 2004-06-11T19:32:00 ET = 140254384. SCLK = 1465674964.105 ET = 140254384. > ``SCLK to ET'' Code
IMPLICIT NONE CHARACTER*(256) MKFILE CHARACTER*(32) UTC CHARACTER*(32) SCLK DOUBLE PRECISION ET INTEGER SCID MKFILE = 'spice_example.tm' CALL FURNSH( MKFILE ) UTC = '2004-06-11T19:32:00' CALL STR2ET( UTC, ET ) WRITE(*,*) 'UTC = ', UTC WRITE(*,*) 'ET = ', ET SCID = -82 SCLK = '1465674964.105' CALL SCS2E( SCID, SCLK, ET ) WRITE(*,*) 'SCLK = ', SCLK WRITE(*,*) 'ET = ', ET ENDMeta-kernel file ``spice_example.tm'':
\begindata KERNELS_TO_LOAD = ( 'kernels/lsk/naif0008.tls' 'kernels/sclk/cas00084.tsc' ) \begintext Step-3: ``Spacecraft State''``Spacecraft State'' Task Statement
``Spacecraft State'' Hints
Find necessary kernel(s) on the NAIF's FTP site. Verify that the kernels contain enough data to compute the state of interest. Use ``brief'' utility program located under ``toolkit/exe'' directory for that. Modify the meta-kernel to load this(these) kernels. Determine the routine(s) needed to compute states. Look at the ``spk.req'' and/or ``SPK'' tutorial presentation. Find the the routine(s) calling sequence specification. Look in the ``spk.req'' and the routine's source code header. Reference the ``naif_ids.req'' and ``frames.req'' and the routine(s) header ``Inputs'' and ``Particulars'' sections to determine proper values of the input arguments of this routine. Add calls to the routine(s), necessary variable declarations and output print statements to the program. Re-compile and re-link it against SPICELIB. ``Spacecraft State'' Solution Steps
The file names can be added to the meta-kernel to get them loaded into the program:
\begindata KERNELS_TO_LOAD = ( 'kernels/lsk/naif0008.tls' 'kernels/sclk/cas00084.tsc' 'kernels/spk/020514_SE_SAT105.bsp' 'kernels/spk/030201AP_SK_SM546_T45.bsp' 'kernels/spk/981005_PLTEPH-DE405S.bsp' 'kernels/spk/sat128.bsp' ) \begintextThe highest level SPICE routine computing states is SPKEZR (``toolkit/src/spicelib/spkezr.f''). We are interested in computing CASSINI position and velocity with respect to the Sun, therefore the target and observer names should be set to 'CASSINI' and 'Sun' (both names can be found in ``naif_ids.req''). The state should be in ecliptic frame, therefore the name of the frame in which the state should be computed is 'ECLIPJ2000' (see ``frames.req'' document.) Since we need only the geometric position, the ABCORR argument of the routine should be set to 'NONE' (see aberration correction discussion in the (``toolkit/src/spicelib/spkezr.f''). header). Putting it all together, we get:
CHARACTER*(32) TARGET CHARACTER*(32) FRAME CHARACTER*(32) CORRTN CHARACTER*(32) OBSERV DOUBLE PRECISION STATE ( 6 ) DOUBLE PRECISION LT ... TARGET = 'CASSINI' FRAME = 'ECLIPJ2000' CORRTN = 'NONE' OBSERV = 'SUN' CALL SPKEZR( TARGET, ET, FRAME, CORRTN, OBSERV, STATE, LT )The updated program with added calls, required declarations and simple print statements produces the following output (the output below was generated by this program compiled with g77 on a PC running Linux; your output may differ slightly in its format and numeric precision):
> ./spice_example UTC = 2004-06-11T19:32:00 ET = 140254384. SCLK = 1465674964.105 ET = 140254384. STATE = -376599062. 1.29448778E+09 -7064853.05 -5.16422619 0.801718909 0.0406030574 > ``Spacecraft State'' Code
IMPLICIT NONE CHARACTER*(256) MKFILE CHARACTER*(32) UTC CHARACTER*(32) SCLK CHARACTER*(32) TARGET CHARACTER*(32) FRAME CHARACTER*(32) CORRTN CHARACTER*(32) OBSERV DOUBLE PRECISION ET DOUBLE PRECISION STATE ( 6 ) DOUBLE PRECISION LT INTEGER SCID MKFILE = 'spice_example.tm' CALL FURNSH( MKFILE ) UTC = '2004-06-11T19:32:00' CALL STR2ET( UTC, ET ) WRITE(*,*) 'UTC = ', UTC WRITE(*,*) 'ET = ', ET SCID = -82 SCLK = '1465674964.105' CALL SCS2E( SCID, SCLK, ET ) WRITE(*,*) 'SCLK = ', SCLK WRITE(*,*) 'ET = ', ET TARGET = 'CASSINI' FRAME = 'ECLIPJ2000' CORRTN = 'NONE' OBSERV = 'SUN' CALL SPKEZR( TARGET, ET, FRAME, CORRTN, OBSERV, STATE, LT ) WRITE(*,*) 'STATE = ', STATE ENDMeta-kernel file ``spice_example.tm'':
\begindata KERNELS_TO_LOAD = ( 'kernels/lsk/naif0008.tls' 'kernels/sclk/cas00084.tsc' 'kernels/spk/020514_SE_SAT105.bsp' 'kernels/spk/030201AP_SK_SM546_T45.bsp' 'kernels/spk/981005_PLTEPH-DE405S.bsp' 'kernels/spk/sat128.bsp' ) \begintext Step-4: ``Sun Direction''``Sun Direction'' Task Statement
``Sun Direction'' Hints
Verify that the orientation data in the kernels have adequate coverage to support computation of the direction of interest. Use ``ckbrief'' utility program located under ``toolkit/exe'' directory for that. Modify the meta-kernel to load this(these) kernels. Determine the proper input arguments for the SPKPOS call to calculate the direction (which is the position portion of the output state). Look through the Frames Kernel find the name of the frame to used. Add calls to the routine(s), necessary variable declarations and output print statements to the program. Re-compile and re-link it against SPICELIB. ``Sun Direction'' Solution Steps
The file names can be added to the meta-kernel to get them loaded into the program:
\begindata KERNELS_TO_LOAD = ( 'kernels/lsk/naif0008.tls' 'kernels/sclk/cas00084.tsc' 'kernels/spk/020514_SE_SAT105.bsp' 'kernels/spk/030201AP_SK_SM546_T45.bsp' 'kernels/spk/981005_PLTEPH-DE405S.bsp' 'kernels/spk/sat128.bsp' 'kernels/ck/04135_04171pc_psiv2.bc' 'kernels/fk/cas_v37.tf' ) \begintextThe same highest level SPICE routine computing positions, SPKPOS, can be used to compute this direction. Since this is the direction of the Sun as seen from the s/c, the target argument should be set to 'Sun' and the observer argument should be set to 'CASSINI'. The name of the INMS frame is 'CASSINI_INMS', the definition and description of this frame are provided in the CASSINI FK file, ``cassini_v02.tf''. Since the apparent, or ``as seen'', position is sought for, the ABCORR argument of the routine should be set to 'LT+S' (see aberration correction discussion in the (``toolkit/src/spicelib/spkpos.f'') header). If desired, the position can then be turned into a unit vector using VHATIP routine (``toolkit/src/spicelib/vhatip.f''). Putting it all together, we get:
DOUBLE PRECISION SUNDIR ( 3 ) ... TARGET = 'SUN' FRAME = 'CASSINI_INMS' CORRTN = 'LT+S' OBSERV = 'CASSINI' CALL SPKPOS( TARGET, ET, FRAME, CORRTN, OBSERV, SUNDIR, LT ) CALL VHATIP( SUNDIR )The updated program with added calls, required declarations and simple print statements produces the following output (the output below was generated by this program compiled with g77 on a PC running Linux; your output may differ slightly in its format and numeric precision):
> ./spice_example UTC = 2004-06-11T19:32:00 ET = 140254384. SCLK = 1465674964.105 ET = 140254384. STATE = -376599062. 1.29448778E+09 -7064853.05 -5.16422619 0.801718909 0.0406030574 SUNDIR = -0.290204017 0.881631188 0.372166733 > ``Sun Direction'' Code
IMPLICIT NONE CHARACTER*(256) MKFILE CHARACTER*(32) UTC CHARACTER*(32) SCLK CHARACTER*(32) TARGET CHARACTER*(32) FRAME CHARACTER*(32) CORRTN CHARACTER*(32) OBSERV DOUBLE PRECISION ET DOUBLE PRECISION STATE ( 6 ) DOUBLE PRECISION SUNDIR ( 3 ) DOUBLE PRECISION LT INTEGER SCID MKFILE = 'spice_example.tm' CALL FURNSH( MKFILE ) UTC = '2004-06-11T19:32:00' CALL STR2ET( UTC, ET ) WRITE(*,*) 'UTC = ', UTC WRITE(*,*) 'ET = ', ET SCID = -82 SCLK = '1465674964.105' CALL SCS2E( SCID, SCLK, ET ) WRITE(*,*) 'SCLK = ', SCLK WRITE(*,*) 'ET = ', ET TARGET = 'CASSINI' FRAME = 'ECLIPJ2000' CORRTN = 'NONE' OBSERV = 'SUN' CALL SPKEZR( TARGET, ET, FRAME, CORRTN, OBSERV, STATE, LT ) WRITE(*,*) 'STATE = ', STATE TARGET = 'SUN' FRAME = 'CASSINI_INMS' CORRTN = 'LT+S' OBSERV = 'CASSINI' CALL SPKPOS( TARGET, ET, FRAME, CORRTN, OBSERV, SUNDIR, LT ) CALL VHATIP( SUNDIR ) WRITE(*,*) 'SUNDIR = ', SUNDIR ENDMeta-kernel file ``spice_example.tm'':
\begindata KERNELS_TO_LOAD = ( 'kernels/lsk/naif0008.tls' 'kernels/sclk/cas00084.tsc' 'kernels/spk/020514_SE_SAT105.bsp' 'kernels/spk/030201AP_SK_SM546_T45.bsp' 'kernels/spk/981005_PLTEPH-DE405S.bsp' 'kernels/spk/sat128.bsp' 'kernels/ck/04135_04171pc_psiv2.bc' 'kernels/fk/cas_v37.tf' ) \begintext Step-5: ``Sub-Spacecraft Point''``Sub-Spacecraft Point'' Task Statement
``Sub-Spacecraft Point'' Hints
Refer to the routine's header to determine the additional kernels needed for this direction computation. Get these kernels from the NAIF's FTP site. Modify the meta-kernel to load this(these) kernels. Determine the proper input arguments for the routine. Refer to the routine's header for that information. Convert the surface point Cartesian vector returned by this routine to latitudinal coordinates. Use ``Permuted Index'' to find the routine that does this conversion. Refer to the routine's header for input/output argument specifications. Since the Cartesian vector from the spacecraft to the sub-spacecraft point is computed in the Phoebe body-fixed frame, it should be transformed into the instrument frame get the direction we are looking for. Refer to ``frames.req'' and/or ``Frames'' tutorial to determine the name of the routine computing transformations and use it to compute transformation from Phoebe body-fixed to the INMS frame. Using ``Permuted Index'' find the routine that multiplies 3x3 matrix by 3d vector and use it to rotate the vector to the instrument frame. Add calls to the routine(s), necessary variable declarations and output print statements to the program. Re-compile and re-link it against SPICELIB. ``Sub-Spacecraft Point'' Solution Steps
Since the s/c is close to Phoebe, light time does not need to be taken into account and, therefore, the ABCORR argument can be set to 'NONE'. In order for SUBPNT to compute the nearest point location, a PCK file containing Phoebe radii has to be loaded into the program (see ``Files'' section of the routine's header.) All other files required for this computation are already being loaded by the program. With PCK file name added to it, the updated meta-kernel will look like this:
\begindata KERNELS_TO_LOAD = ( 'kernels/lsk/naif0008.tls' 'kernels/sclk/cas00084.tsc' 'kernels/spk/020514_SE_SAT105.bsp' 'kernels/spk/030201AP_SK_SM546_T45.bsp' 'kernels/spk/981005_PLTEPH-DE405S.bsp' 'kernels/spk/sat128.bsp' 'kernels/ck/04135_04171pc_psiv2.bc' 'kernels/fk/cas_v37.tf' 'kernels/pck/cpck05Mar2004.tpc' ) \begintextThe sub-spacecraft point Cartesian vector can be converted to planetocentric radius, longitude and latitude using the RECLAT routine (``toolkit/src/spicelib/reclat.f''). The vector from the spacecraft to the sub-spacecraft point returned by SUBPNT has to be rotated from the body-fixed frame to the instrument frame. The name of the routine that computes 3x3 matrices rotating vectors from one frame to another is PXFORM (``toolkit/src/spicelib/pxform.f''). In our case the FROM argument should be set to 'IAU_PHOEBE' and the TO argument should be set to 'CASSINI_INMS' The vector should be then multiplied by this matrix to rotate it to the instrument frame. The MXV routine performs that function (``toolkit/src/spicelib/mxv.f'') After applying the rotation, normalize the resultant vector using the VHATIP routine. For output the longitude and latitude angles returned by RECLAT in radians can be converted to degrees by multiplying by DPR function (``toolkit/src/spicelib/dpr.f''). Putting it all together, we get:
CHARACTER*(32) METHOD CHARACTER*(32) FROMFR CHARACTER*(32) TOFR DOUBLE PRECISION SPOINT ( 3 ) DOUBLE PRECISION TRGEPC DOUBLE PRECISION SRFVEC ( 3 ) DOUBLE PRECISION SRAD DOUBLE PRECISION SLON DOUBLE PRECISION SLAT DOUBLE PRECISION SBPDIR ( 3 ) DOUBLE PRECISION M2IMAT ( 3, 3 ) DOUBLE PRECISION DPR ... METHOD = 'NEAR POINT: ELLIPSOID' TARGET = 'PHOEBE' FRAME = 'IAU_PHOEBE' CORRTN = 'NONE' OBSERV = 'CASSINI' CALL SUBPNT ( METHOD, TARGET, ET, FRAME, CORRTN, OBSERV, SPOINT, TRGEPC, SRFVEC ) CALL RECLAT( SPOINT, SRAD, SLON, SLAT ) FROMFR = 'IAU_PHOEBE' TOFR = 'CASSINI_INMS' CALL PXFORM( FROMFR, TOFR, ET, M2IMAT ) CALL MXV ( M2IMAT, SRFVEC, SBPDIR ) CALL VHATIP( SBPDIR ) WRITE(*,*) 'LON = ', SLON*DPR() WRITE(*,*) 'LAT = ', SLAT*DPR() WRITE(*,*) 'SBPDIR = ', SBPDIRThe updated program with added calls, required declarations and simple print statements produces the following output (the output below was generated by this program compiled with g77 on a PC running Linux; your output may differ slightly in its format and numeric precision):
> ./spice_example UTC = 2004-06-11T19:32:00 ET = 140254384. SCLK = 1465674964.105 ET = 140254384. STATE = -376599062. 1.29448778E+09 -7064853.05 -5.16422619 0.801718909 0.0406030574 SUNDIR = -0.290204017 0.881631188 0.372166733 LON = 23.4231584 LAT = 3.70979716 SBPDIR = -0.000776207056 -0.999873199 -0.015905456 > ``Sub-Spacecraft Point'' Code
IMPLICIT NONE CHARACTER*(256) MKFILE CHARACTER*(32) UTC CHARACTER*(32) SCLK CHARACTER*(32) TARGET CHARACTER*(32) FRAME CHARACTER*(32) CORRTN CHARACTER*(32) OBSERV CHARACTER*(32) METHOD CHARACTER*(32) FROMFR CHARACTER*(32) TOFR DOUBLE PRECISION ET DOUBLE PRECISION STATE ( 6 ) DOUBLE PRECISION SUNDIR ( 3 ) DOUBLE PRECISION LT DOUBLE PRECISION SPOINT ( 3 ) DOUBLE PRECISION TRGEPC DOUBLE PRECISION SRFVEC ( 3 ) DOUBLE PRECISION SRAD DOUBLE PRECISION SLON DOUBLE PRECISION SLAT DOUBLE PRECISION SBPDIR ( 3 ) DOUBLE PRECISION M2IMAT ( 3, 3 ) INTEGER SCID DOUBLE PRECISION DPR MKFILE = 'spice_example.tm' CALL FURNSH( MKFILE ) UTC = '2004-06-11T19:32:00' CALL STR2ET( UTC, ET ) WRITE(*,*) 'UTC = ', UTC WRITE(*,*) 'ET = ', ET SCID = -82 SCLK = '1465674964.105' CALL SCS2E( SCID, SCLK, ET ) WRITE(*,*) 'SCLK = ', SCLK WRITE(*,*) 'ET = ', ET TARGET = 'CASSINI' FRAME = 'ECLIPJ2000' CORRTN = 'NONE' OBSERV = 'SUN' CALL SPKEZR( TARGET, ET, FRAME, CORRTN, OBSERV, STATE, LT ) WRITE(*,*) 'STATE = ', STATE TARGET = 'SUN' FRAME = 'CASSINI_INMS' CORRTN = 'LT+S' OBSERV = 'CASSINI' CALL SPKPOS( TARGET, ET, FRAME, CORRTN, OBSERV, SUNDIR, LT ) CALL VHATIP( SUNDIR ) WRITE(*,*) 'SUNDIR = ', SUNDIR METHOD = 'NEAR POINT: ELLIPSOID' TARGET = 'PHOEBE' FRAME = 'IAU_PHOEBE' CORRTN = 'NONE' OBSERV = 'CASSINI' CALL SUBPNT ( METHOD, TARGET, ET, FRAME, CORRTN, OBSERV, SPOINT, TRGEPC, SRFVEC ) CALL RECLAT( SPOINT, SRAD, SLON, SLAT ) FROMFR = 'IAU_PHOEBE' TOFR = 'CASSINI_INMS' CALL PXFORM( FROMFR, TOFR, ET, M2IMAT ) CALL MXV ( M2IMAT, SRFVEC, SBPDIR ) CALL VHATIP( SBPDIR ) WRITE(*,*) 'LON = ', SLON*DPR() WRITE(*,*) 'LAT = ', SLAT*DPR() WRITE(*,*) 'SBPDIR = ', SBPDIR ENDMeta-kernel file ``spice_example.tm'':
\begindata KERNELS_TO_LOAD = ( 'kernels/lsk/naif0008.tls' 'kernels/sclk/cas00084.tsc' 'kernels/spk/020514_SE_SAT105.bsp' 'kernels/spk/030201AP_SK_SM546_T45.bsp' 'kernels/spk/981005_PLTEPH-DE405S.bsp' 'kernels/spk/sat128.bsp' 'kernels/ck/04135_04171pc_psiv2.bc' 'kernels/fk/cas_v37.tf' 'kernels/pck/cpck05Mar2004.tpc' ) \begintext Step-6: ``Spacecraft Velocity''``Spacecraft Velocity'' Task Statement
``Spacecraft Velocity'' Hints
Since the velocity vector is computed in the inertial frame, it should be rotated to the instrument frame. Look at the previous step the routine that compute necessary rotation and rotate vectors. Add calls to the routine(s), necessary variable declarations and output print statements to the program. Re-compile and re-link it against SPICELIB. ``Spacecraft Velocity'' Solution Steps
The spacecraft velocity vector is the last three components of the state returned by SPKEZR. To compute velocity of CASSINI with respect to Phoebe in the J2000 inertial frame the SPKEZR arguments should be set to 'CASSINI' (TARG), 'PHOEBE' (OBS), 'J2000' (REF) and 'NONE' (ABCORR). For convenience the velocity can be copied from the output state in to a 3d vector using the VPACK routine (``toolkit/src/spicelib/vpack.f''). The computed velocity vector has to be rotated from the J2000 frame to the instrument frame. The PXFORM routine used in the previous step can be used to compute the rotation matrix needed for that. In this case the frame name arguments should be set to 'J2000' (FROM) and 'CASSINI_INMS' (TO). As in the previous step the difference vector should be then multiplied by this rotation matrix using the MXV routine. After applying the rotation, normalize the resultant vector using the VHAT routine. Putting it all together, we get:
DOUBLE PRECISION SCVDIR ( 3 ) DOUBLE PRECISION TMPDIR ( 3 ) DOUBLE PRECISION J2IMAT( 3, 3 ) ... TARGET = 'CASSINI' FRAME = 'J2000' CORRTN = 'NONE' OBSERV = 'PHOEBE' CALL SPKEZR( TARGET, ET, FRAME, CORRTN, OBSERV, STATE, LT ) CALL VPACK ( STATE(4), STATE(5), STATE(6), SCVDIR ) FROMFR = 'J2000' TOFR = 'CASSINI_INMS' CALL PXFORM( FROMFR, TOFR, ET, J2IMAT ) CALL MXV ( J2IMAT, SCVDIR, TMPDIR ) CALL VHAT ( TMPDIR, SCVDIR ) WRITE(*,*) 'SCVDIR = ', SCVDIRThe updated program with added calls, required declarations and simple print statements produces the following output (the output below was generated by this program compiled with g77 on a PC running Linux; your output may differ slightly in its format and numeric precision):
> ./spice_example UTC = 2004-06-11T19:32:00 ET = 140254384. SCLK = 1465674964.105 ET = 140254384. STATE = -376599062. 1.29448778E+09 -7064853.05 -5.16422619 0.801718909 0.0406030574 SUNDIR = -0.290204017 0.881631188 0.372166733 LON = 23.4231584 LAT = 3.70979716 SBPDIR = -0.000776207056 -0.999873199 -0.015905456 SCVDIR = 0.395784873 -0.292807673 0.870412546 >Note that computing the spacecraft velocity in the instrument frame by a single call to SPKEZR by specifying 'CASSINI_INMS' in the REF argument returns an incorrect result. Such computation will take into account the spacecraft angular velocity from the CK files, which should not be considered in this case. ``Spacecraft Velocity'' Code Program ``spice_example.f'':
IMPLICIT NONE CHARACTER*(256) MKFILE CHARACTER*(32) UTC CHARACTER*(32) SCLK CHARACTER*(32) TARGET CHARACTER*(32) FRAME CHARACTER*(32) CORRTN CHARACTER*(32) OBSERV CHARACTER*(32) METHOD CHARACTER*(32) FROMFR CHARACTER*(32) TOFR DOUBLE PRECISION ET DOUBLE PRECISION STATE ( 6 ) DOUBLE PRECISION SUNDIR ( 3 ) DOUBLE PRECISION LT DOUBLE PRECISION SPOINT ( 3 ) DOUBLE PRECISION TRGEPC DOUBLE PRECISION SRFVEC ( 3 ) DOUBLE PRECISION SRAD DOUBLE PRECISION SLON DOUBLE PRECISION SLAT DOUBLE PRECISION SBPDIR ( 3 ) DOUBLE PRECISION M2IMAT ( 3, 3 ) DOUBLE PRECISION SCVDIR ( 3 ) DOUBLE PRECISION TMPDIR ( 3 ) DOUBLE PRECISION J2IMAT ( 3, 3 ) INTEGER SCID DOUBLE PRECISION DPR MKFILE = 'spice_example.tm' CALL FURNSH( MKFILE ) UTC = '2004-06-11T19:32:00' CALL STR2ET( UTC, ET ) WRITE(*,*) 'UTC = ', UTC WRITE(*,*) 'ET = ', ET SCID = -82 SCLK = '1465674964.105' CALL SCS2E( SCID, SCLK, ET ) WRITE(*,*) 'SCLK = ', SCLK WRITE(*,*) 'ET = ', ET TARGET = 'CASSINI' FRAME = 'ECLIPJ2000' CORRTN = 'NONE' OBSERV = 'SUN' CALL SPKEZR( TARGET, ET, FRAME, CORRTN, OBSERV, STATE, LT ) WRITE(*,*) 'STATE = ', STATE TARGET = 'SUN' FRAME = 'CASSINI_INMS' CORRTN = 'LT+S' OBSERV = 'CASSINI' CALL SPKPOS( TARGET, ET, FRAME, CORRTN, OBSERV, SUNDIR, LT ) CALL VHATIP( SUNDIR ) WRITE(*,*) 'SUNDIR = ', SUNDIR METHOD = 'NEAR POINT: ELLIPSOID' TARGET = 'PHOEBE' FRAME = 'IAU_PHOEBE' CORRTN = 'NONE' OBSERV = 'CASSINI' CALL SUBPNT ( METHOD, TARGET, ET, FRAME, CORRTN, OBSERV, SPOINT, TRGEPC, SRFVEC ) CALL RECLAT( SPOINT, SRAD, SLON, SLAT ) FROMFR = 'IAU_PHOEBE' TOFR = 'CASSINI_INMS' CALL PXFORM( FROMFR, TOFR, ET, M2IMAT ) CALL MXV ( M2IMAT, SRFVEC, SBPDIR ) CALL VHATIP( SBPDIR ) WRITE(*,*) 'LON = ', SLON*DPR() WRITE(*,*) 'LAT = ', SLAT*DPR() WRITE(*,*) 'SBPDIR = ', SBPDIR TARGET = 'CASSINI' FRAME = 'J2000' CORRTN = 'NONE' OBSERV = 'PHOEBE' CALL SPKEZR( TARGET, ET, FRAME, CORRTN, OBSERV, STATE, LT ) CALL VPACK ( STATE(4), STATE(5), STATE(6), SCVDIR ) FROMFR = 'J2000' TOFR = 'CASSINI_INMS' CALL PXFORM( FROMFR, TOFR, ET, J2IMAT ) CALL MXV ( J2IMAT, SCVDIR, TMPDIR ) CALL VHAT ( TMPDIR, SCVDIR ) WRITE(*,*) 'SCVDIR = ', SCVDIR ENDMeta-kernel file ``spice_example.tm'':
\begindata KERNELS_TO_LOAD = ( 'kernels/lsk/naif0008.tls' 'kernels/sclk/cas00084.tsc' 'kernels/spk/020514_SE_SAT105.bsp' 'kernels/spk/030201AP_SK_SM546_T45.bsp' 'kernels/spk/981005_PLTEPH-DE405S.bsp' 'kernels/spk/sat128.bsp' 'kernels/ck/04135_04171pc_psiv2.bc' 'kernels/fk/cas_v37.tf' 'kernels/pck/cpck05Mar2004.tpc' ) \begintext |