matthewholman / assist

ASSIST is a software package for ephemeris-quality integrations of test particles.
https://assist.readthedocs.io/
GNU General Public License v3.0
24 stars 10 forks source link

Use DAF file record's FWARD integer to skip to first summary record #90

Closed moeyensj closed 1 year ago

moeyensj commented 1 year ago

Hi Hanno and Matt,

I don't know if its particularly helpful but DAF files have a specific attribute that indicates where the first summary record starts. This PR removes the while loop in spk.c in favor of a direct seek to the location of the first summary record using the FWARD attribute (https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/daf.html#The%20File%20Record).

I've tested this with the example in assist/examples/asteroid and get the same results. An off-by-one error in the interpretation of the location of the first record yields different state vectors and an error.

Results of example propagation using code on "main":

┌─(assist_py310)[moeyensj][prometheus][±][main ✓][~/repos/forked/assist/examples/asteroid]
└─▪ ./rebound 
-0.170145 0.895984 0.388437 
t = 2458869.5
particles[0]:   x = 3.390374780107      y = -0.765446062257     z = -0.443591042042
particles[1]:   x = 1.878434585839      y = 1.068680422732      z = 0.002968378656
particles[2]:   x = -0.520279833813     y = 0.779825882254      z = 0.341506750710
t = 2458889.5
particles[0]:   x = 3.432579176258      y = -0.611129773171     z = -0.382101923128
particles[1]:   x = 1.749132826943      y = 1.276964523800      z = 0.002955203013
particles[2]:   x = -0.817371985524     y = 0.567022464475      z = 0.253793997588
t = 2458909.5
particles[0]:   x = 3.465399787550      y = -0.455130472630     z = -0.319563145268
particles[1]:   x = 1.599405415013      y = 1.470451838759      z = 0.002941053960
particles[2]:   x = -1.025026707189     y = 0.291815842383      z = 0.138139433099
t = 2458929.5
particles[0]:   x = 3.488764039215      y = -0.297876494459     z = -0.256147101651
particles[1]:   x = 1.431153526114      y = 1.647035041297      z = 0.002925964498
particles[2]:   x = -1.134780116148     y = -0.011574348097     z = 0.009152298995
t = 2458949.5
particles[0]:   x = 3.502615350000      y = -0.139797159344     z = -0.192027205249
particles[1]:   x = 1.246512134000      y = 1.804897617429      z = 0.002909955990
particles[2]:   x = -1.152492124123     y = -0.314204896422     z = -0.120654880203

Results with this PR:

    lseek(fd, (record.file.fward - 1) * record_length, SEEK_SET);
┌─(assist_py310)[moeyensj][prometheus][±][record-start ✓][~/repos/forked/assist/examples/asteroid]
└─▪ ./rebound 
-0.170145 0.895984 0.388437 
t = 2458869.5
particles[0]:   x = 3.390374780107      y = -0.765446062257     z = -0.443591042042
particles[1]:   x = 1.878434585839      y = 1.068680422732      z = 0.002968378656
particles[2]:   x = -0.520279833813     y = 0.779825882254      z = 0.341506750710
t = 2458889.5
particles[0]:   x = 3.432579176258      y = -0.611129773171     z = -0.382101923128
particles[1]:   x = 1.749132826943      y = 1.276964523800      z = 0.002955203013
particles[2]:   x = -0.817371985524     y = 0.567022464475      z = 0.253793997588
t = 2458909.5
particles[0]:   x = 3.465399787550      y = -0.455130472630     z = -0.319563145268
particles[1]:   x = 1.599405415013      y = 1.470451838759      z = 0.002941053960
particles[2]:   x = -1.025026707189     y = 0.291815842383      z = 0.138139433099
t = 2458929.5
particles[0]:   x = 3.488764039215      y = -0.297876494459     z = -0.256147101651
particles[1]:   x = 1.431153526114      y = 1.647035041297      z = 0.002925964498
particles[2]:   x = -1.134780116148     y = -0.011574348097     z = 0.009152298995
t = 2458949.5
particles[0]:   x = 3.502615350000      y = -0.139797159344     z = -0.192027205249
particles[1]:   x = 1.246512134000      y = 1.804897617429      z = 0.002909955990
particles[2]:   x = -1.152492124123     y = -0.314204896422     z = -0.120654880203

Results with intentional off-by-one error:

    lseek(fd, record.file.fward * record_length, SEEK_SET);
┌─(assist_py310)[moeyensj][prometheus][±][record-start U:1 ✗][~/repos/forked/assist/examples/asteroid]
└─▪ ./rebound 
Error parsing DAF/SPL file. Cannot find summary block.
(ASSIST) The JPL asteroid ephemeris file has not been found. Asteroid forces have been disabled.
-0.170145 0.895984 0.388437 
t = 2458869.5
particles[0]:   x = 3.390374780111      y = -0.765446062256     z = -0.443591042042
particles[1]:   x = 1.878434585839      y = 1.068680422732      z = 0.002968378655
particles[2]:   x = -0.520279832264     y = 0.779825877555      z = 0.341506750642
t = 2458889.5
particles[0]:   x = 3.432579176275      y = -0.611129773168     z = -0.382101923127
particles[1]:   x = 1.749132826947      y = 1.276964523800      z = 0.002955203007
particles[2]:   x = -0.817371981027     y = 0.567022452699      z = 0.253793996940
t = 2458909.5
particles[0]:   x = 3.465399787589      y = -0.455130472623     z = -0.319563145264
particles[1]:   x = 1.599405415021      y = 1.470451838759      z = 0.002941053945
particles[2]:   x = -1.025026697579     y = 0.291815823092      z = 0.138139431156
t = 2458929.5
particles[0]:   x = 3.488764039283      y = -0.297876494444     z = -0.256147101644
particles[1]:   x = 1.431153526128      y = 1.647035041298      z = 0.002925964469
particles[2]:   x = -1.134780098408     y = -0.011574374140     z = 0.009152295424
t = 2458949.5
particles[0]:   x = 3.502615350104      y = -0.139797159319     z = -0.192027205236
particles[1]:   x = 1.246512134024      y = 1.804897617431      z = 0.002909955943
particles[2]:   x = -1.152492095509     y = -0.314204927054     z = -0.120654885146
moeyensj commented 1 year ago

We could also use the BWARD argument to replace the while-loop with a for-loop (see https://github.com/matthewholman/assist/blob/main/src/spk.c#L137) but I suspect that will have a negligible impact on performance

hannorein commented 1 year ago

Cool! Thank you. I have to admit, I don't understand what's going on with the offset? Was there a bug before?

moeyensj commented 1 year ago

There wasn't a bug before! I was just trying to test that the direct seek was working and that it had an effect on the test case, so I intentionally set the record number to be incorrect and then reran the test hoping to see different results. To be clear, as far as I can tell this PR correctly reads the development ephemeris and it matches the test case when compared to the code on the main branch.

It's definitely a weird way to do things but I wasn't sure how best to test the changes.

hannorein commented 1 year ago

Ok. Got it! Thanks. Well, I think this is safe to merge in then. I don't see any downside. Are you ok with this too, @matthewholman?

matthewholman commented 1 year ago

This sounds good to me. Many thanks, @moeyensj !