SPICE toolkit occultation finder doesn't work with spacecraft

Question

I've to know if the sun is occulted by earth from my spacecraft which is in low earth orbit with these parameters (a:6700 km, e:0.0, i:5°, ...). I run the following method with search window start at 2021 MAR 03 20:00 and stop at 2021 MAR 04 23:00 :

gfoclt_c("ANY", "earth", "ELLIPSOID","ITRF93", "sun", "ELLIPSOID", "IAU_sun", "LT", "MySpacecraft", 120.0, &cnfine, &results);

This method run without exception but result doesn't contain data.

I tried this method to detect sun eclipse by moon from earth center and it works perfectly, the problem seems to be my spacecraft but I don't understand what. Bellow my spk, ck and frame informations : Thank you for your help.

Summary for CK file: MySpacecraft.ck
Leapseconds File   : latest_leapseconds.tls
SCLK File          : MySpacecraft.tsc
Summary Type       : Entire File
 
--------------------------------------------------------------------------------
   Segment ID     : Seg1
   Instrument Code: -1789000
   Spacecraft     : Body -1789
   Reference Frame: Frame 1, J2000
   CK Data Type   : Type 5
      Description : Continuous Pointing: MEX/Rosetta Polynomial Interpolation
   Available Data : Pointing and Angular Velocity
   UTC Start Time : 2021 MAR 03 19:05:40.085
   UTC Stop Time  : 2021 MAR 05 00:00:00.841
   SCLK Start Time: 1/0005339209:17746
   SCLK Stop Time : 1/0005443270:01719

Summary for SPK file: MySpacecraft.spk
Leapseconds File    : latest_leapseconds.tls
Summary Type        : Entire File
 
--------------------------------------------------------------------------------
   Segment ID     : Seg1
   Target Body    : Body -1789
   Center Body    : Body 399, EARTH
   Reference frame: Frame 1, J2000
   SPK Data Type  : Type 9
      Description : Discrete States, Unevenly Spaced, Lagrange Interpolation
   UTC Start Time : 2021 MAR 03 19:05:40.085
   UTC Stop Time  : 2021 MAR 05 00:00:00.841
   ET Start Time  : 2021 MAR 03 19:06:49.270
   ET Stop time   : 2021 MAR 05 00:01:10.026

KPL/FK

\begindata
      
 
      FRAME_MySpacecraft_SPACECRAFT   = -1789000
      FRAME_-1789000_NAME      = 'MySpacecraft_SPACECRAFT'
      FRAME_-1789000_CLASS     =  3
      FRAME_-1789000_CLASS_ID  = -1789000
      FRAME_-1789000_CENTER    = -1789
 
      CK_-1789000_SCLK         = -1789
      CK_-1789000_SPK          = -1789
 
      OBJECT_-1789_FRAME       = 'MySpacecraft_SPACECRAFT'

      NAIF_BODY_NAME              += 'MySpacecraft_SPACECRAFT'
      NAIF_BODY_CODE              += -1789000
      
      NAIF_BODY_NAME              += 'MySpacecraft'
      NAIF_BODY_CODE              += -1789

\begintext

I also load these kernels :

• de440s.bsp
• earth_assoc_itrf93.tf
• earth_fixed.tf
• earth_latest_high_prec.bpc
• earthstns_itrf93_201023.bsp
• earth_topo_201023.tf
• geophysical.ker
• gm_de431.tpc
• latest_leapseconds.tls
• moon_080317.tf
• moon_assoc_me.tf
• moon_assoc_pa.tf
• moon_pa_de421_1900-2050.bpc
• pck00010.tpc

Answer 1

The error comes from bad formatted spk file. My application uses international unit system and spice uses km and km/s for state vector. When the occultation finder runs it see my spacecraft thousand times further than he really is, so earth apparent angular size is much smaller and occultation never occurs. Now it works perfectly. Thank you for your help you put me on the track!

Answer 2

As in your comments, plugging in: sma = 7480km, e = 0.09, i = 5.5 deg, RAAN = 0.2 deg, AOP = 5.95 deg, true anomaly = 6.26 (I used true anomaly because its close since e is small and I don't remember the conversion to mean anomaly off the top of my head, can fix if you'd like) at epoch 2021-03-03 22:10:40 TDB

Plugging in these numbers gives these eclipse results:

With custom function:

Using SPICE's gfoclt:

Visual reference for the eclipse entrance:

Eclipse exit:

For the custom function, I'm using SPICE here just for ephemeris data (DE432), not for the eclipse calculation. Eclipse calculations are done with a custom function that assumes spherical Earth and no atmospheric diffractions (purely geometrical, but with the one-way light time aberration correction for the Sun's position w.r.t Earth).

Here is what the geometry calculations look like for penumbra:

To summarize, the vector pointing from the Sun to Earth (with one way light time correction) is the 3D eclipse cone's center, and we use the magnitude of the projection of the spacecraft's position vector w.r.t Earth onto the sun vector (cyan colored vector $p\vec{ro}j$) to calculate the height of the cone. Then we use the cone's half angle to calculate its radius at that height. If the magnitude of the Spacecraft's rejection vector (lime colored vector $\vec{rej}$) is less than the radius of the cone, then spacecraft is in penumbra.

Similarly, here is the umbra geometry:

Here is a more detailed paper on this geometry from NASA Johnson from 1995: https://ntrs.nasa.gov/api/citations/19950023025/downloads/19950023025.pdf

The times of entrance and exit (from the custom function) are not "exact", since those times are just being found from the time and states arrays that the ODE solver used to propagate the orbit. There is not iterative (probably root) solver that refines that time. Note that this is an adaptive step size solver (which is why the dots aren't evenly spaced in this plot, where 1->penumbra, 2->umbra, -1->no eclipse):

I've been working on this this week, so it still may have small bugs, so I currently have these eclipse calculations in a separate branch from main in this GitHub repo: https://github.com/alfonsogonzalez/AWP/commit/a1c251b8942d492a40936cd759a1c676f2277d59

I'll be making a video going all over this analysis this weekend, I can share a link with you if you'd like (this is why I made these geometric diagrams)

And here is the main script for this example problem (orbit propagated for 1 period):

# AWP library
from Spacecraft import Spacecraft as SC
import spice_data  as sd
import spice_tools as st

# 3rd party libraries
import spiceypy as spice
    
if __name__ == '__main__':
    spice.furnsh( sd.leapseconds_kernel )
    spice.furnsh( sd.de432 )
    spice.furnsh( sd.pck00010 )

    sc = SC( {
        'date0': '2021-03-03 22:10:35 TDB',
        'coes' : [ 7480.0, 0.09, 5.5, 6.26, 5.95, 0.2 ],
        'tspan': '40',
        } )

    st.write_bsp( sc.ets, sc.states[ :, :6 ], { 'bsp_fn': 'leo.bsp' } )
    spice.furnsh( 'leo.bsp' )

    et0  = spice.str2et( '2021-03-03 22:10:40 TDB' )
    etf  = spice.str2et( '2021-03-04 TDB' )

    timecell = spice.utils.support_types.SPICEDOUBLE_CELL( 2 )
    spice.appndd( et0, timecell )
    spice.appndd( etf, timecell )

    cell = spice.gfoclt( 'ANY', '399', 'ELLIPSOID', 'IAU_EARTH',
        '10', 'ELLIPSOID', 'IAU_SUN', 'LT', '-999', 120.0, timecell )
    ets_SPICE = spice.wnfetd( cell, 0 )
    cal0      = spice.et2utc( ets_SPICE[ 0 ], 'C', 1 )
    cal1      = spice.et2utc( ets_SPICE[ 1 ], 'C', 1 )
    print( '\n*** SPICE RESULTS ***' )
    print( f'{cal0} --> {cal1}')

    sc.calc_eclipses( vv = True )
  
    sc.plot_eclipse_array( {
        'time_unit': 'seconds',
        'show'     : True
    } )