4
$\begingroup$

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
$\endgroup$
5
  • 1
    $\begingroup$ This might sound trivial, but could it be that there is no occultation? If you have access to the orbital parameters, you could simulate that same trajectory (albeit in lower fidelity) in GMAT and check if there are any occultations according to it. $\endgroup$
    – ChrisR
    Commented Nov 18, 2021 at 16:43
  • 1
    $\begingroup$ With these orbital parameters at this epoch in 24h I should have occultations but it's a good Idea to cross check my data with GMAT. $\endgroup$
    – sl20
    Commented Nov 18, 2021 at 17:50
  • $\begingroup$ @sl20 can you post all the orbital elements of your spacecraft? I have an umbra / penumbra finder I can run it through real quick. A 5 degree inclination LEO spacecraft should go into eclipse, but its possible for a high inclination orbit to orbit around the sun vector (terminator orbit) thus no eclipse $\endgroup$ Commented Nov 18, 2021 at 19:12
  • 1
    $\begingroup$ @AlfonsoGonzalez thank for your help. So we have a:7480 km e: 0.09 i:5.5° RAAN : 0.2° perigeeArgument:5.95° MeanAnomaly : 6.26° at epoch "2021-03-03 22:10:40 (TDB)" $\endgroup$
    – sl20
    Commented Nov 18, 2021 at 19:49
  • $\begingroup$ @sl20 Just posted an answer $\endgroup$ Commented Nov 18, 2021 at 20:45

2 Answers 2

4
$\begingroup$

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!

$\endgroup$
1
$\begingroup$

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: enter image description here

Using SPICE's gfoclt: enter image description here

Visual reference for the eclipse entrance: enter image description here

Eclipse exit: enter image description here

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:

enter image description here

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:

enter image description here

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): enter image description here

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
    } )
$\endgroup$
2
  • 2
    $\begingroup$ Thank you for your really interesting explanations but that confirm I've problem with the spice geometry finder :( $\endgroup$
    – sl20
    Commented Nov 18, 2021 at 21:04
  • 1
    $\begingroup$ @sl20 I just updated the answer to include the call to gfoclt which seemed to work for me. I furnshed pck00010.tpc for the IAU_EARTH, IAU_SUN frames and for the ellipsoid models of Earth and Sun. The SPICE results are right below the original screenshot of the terminal output $\endgroup$ Commented Nov 18, 2021 at 21:45

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.