2
$\begingroup$

I want to figure out lunar eclipses as seen from a point on the Moon using SPICE. To generate ephemeris for the point on the Moon, I use the code provided in this answer. Instead of IAU_MOON I use MOON_ME; the adjusted code provided by Alfonso Gonzalez is below.

The relevant code in Matlab to look for occultations is the second piece of code below.

But I keep getting:

Error using mice
SPICE(SPKINSUFFDATA): [gfoclt_c->GFOCLT->GFOCCE->ZZGFSOLV->ZZGFOCST->SPKEZP->SPKGPS] Insufficient ephemeris data has been loaded to compute the position of 399
(EARTH) relative to 301100 at the ephemeris epoch 2004 JAN 01 09:33:04.183.

How should i know which kernel am i missing. 301100 is given wrt J2000, so there should not be anything else needed. If I put Moon instead of 301100, it works fine. The loaded metakernel is the third piece of code below.


Code to create ephemeris:

'''
Create SPK kernel from Lunar lat-lon coordinates
'''

from sys import path
import spiceypy as spice
import numpy as np

r2d = np.pi / 180.0

if __name__ == '__main__':
    spice.furnsh( '/home/me/kernels/naif0011.tls' )
    spice.furnsh( '/home/me/kernels/pck00010.tpc' )
    spice.furnsh( '/home/me/kernels/de421.bsp' )
    spice.furnsh( '/home/me/kernels/moon_pa_de421_1900-2050.bpc' )
    spice.furnsh( '/home/me/kernels/moon_080317.tf' )

    lat     = 75.0 * r2d
    lon     = 135.0 * r2d
    r       = 1737.4
    et0     = spice.str2et( '2003-01-01' )
    etf     = spice.str2et( '2006-01-01' )
    delta_t = 1000
    ets     = np.arange( et0, etf, delta_t )
    states  = np.zeros( ( len( ets ), 6 ) )
    pos_iau_moon = spice.latrec( r, lon, lat )
    # this converts the geocentric coordinates of the surface point from to the J2000 frame for each time one wants the ephemeris calculated
    for n in range( len( ets ) ):
        states[ n, :3 ] = np.dot(
            spice.pxform( 'MOON_ME', 'J2000', ets[ n ] ),
            pos_iau_moon )
    # here one makes the differences between the subsequent positions of the surface point wrt J2000 frame and divides by time difference leading to velocities
    for n in range( len( ets ) - 1 ):
        states[ n, 3: ] = ( states[ n + 1, :3 ] - states[ n, :3 ] ) / delta_t
    # this takes care of the last time instant
    pos_np1 = np.dot(
            spice.pxform( 'MOON_ME', 'J2000', ets[ -1 ] + delta_t ),
            pos_iau_moon )
    states[ -1, 3: ] = ( pos_np1 - states[ -1, :3 ] ) / delta_t
    # here one creates a new SPK file
    handle   = spice.spkopn( '/home/me/kernels/lunar_point.bsp', 'SPK_file', 0 )
    point_id = 301100
    center   = 301
    frame    = 'J2000'
    degree   = 5
    # this function Write a type 9 segment to an SPK file. whatever that is
    #degree = degree of the Lagrange polynomials used to interpolate the states.
    spice.spkw09( handle, point_id, center, frame,
        ets[ 0 ], ets[ -1 ], '0', 5, len( ets ),
        states.tolist(), ets.tolist() )
    # closes the SPK file
    spice.spkcls( handle )

Matlab code

cspice_furnsh(C.fn_spice_ker);
MAXWIN  =  1000;
TIMFMT  = 'YYYY-MON-DD HR:MN:SC.###### (TDB) ::TDB ::RND';
t_start  = cspice_str2et(C.t_utc_start);
t_end    = t_start + C.t_duration;
t_v      = t_start:C.t_density:t_end;
t_n      = length(t_v);
cnfine = cspice_wninsd( t_start, t_end );

step    = 180.;

occtyp  = 'any';
front   = 'earth';
fshape  = 'ellipsoid';
fframe  = 'ITRF93';
back    = 'sun';
bshape  = 'ellipsoid';
bframe  = 'iau_sun';
obsrvr  = '301100';
abcorr  = 'none';
result = cspice_gfoclt( occtyp, front, fshape, fframe, ...
                              back, bshape, bframe,          ...
                              abcorr, obsrvr, step, cnfine,  ...
                              MAXWIN);
for i=1:numel(result)/2
     [left, right] = cspice_wnfetd( result, i );
     output = cspice_timout( [left,right], TIMFMT );
     if( isequal( left, right) )
        disp( ['Event time: ' output(1,:)] )
     else
        disp( ['From : ' output(1,:)] )
        disp( ['To   : ' output(2,:)] )
        disp( ' ')
     end
end

Metakernel spice.ker (provided by the line cspice_furnsh(C.fn_spice_ker);)

KPL/MK
\begindata
    PATH_VALUES     = ( '/home/me/kernels' )
    PATH_SYMBOLS    = ( 'KRL' )
    KERNELS_TO_LOAD = ( '$KRL/moon_080317.tf',
						'$KRL/naif0011.tls',
                        '$KRL/pck00010.tpc',
						'$KRL/earth_000101_210920_210630.bpc',
                        '$KRL/de421.bsp',
						'$KRL/moon_pa_de421_1900-2050.bpc',
                        '$KRL/lunar_point.bsp',
						'$KRL/moon_topo.fk'
                      )

\begintext
$\endgroup$
1
  • 1
    $\begingroup$ This may have to do with the time coverage (maybe its like a millisecond off), since you've loaded de421.bsp and lunar_point.bsp. You can check the time coverage of lunar_point.bsp either by using the brief utility (naif.jpl.nasa.gov/naif/utilities.html) or using the spkcov and wnfetd functions $\endgroup$ Jul 2 at 15:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.