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$ Commented Jul 2, 2021 at 15:57

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.