Merge pull request #41 from jacobwilliams/develop

Develop

Author: jacobwilliams

LICENSE

@@ -1,7 +1,7 @@
 Fortran Astrodynamics Toolkit
 https://github.com/jacobwilliams/Fortran-Astrodynamics-Toolkit
 
-Copyright (c) 2014-2022, Jacob Williams
+Copyright (c) 2014-2025, Jacob Williams
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without modification,

src/celestial_body_module.f90

@@ -15,11 +15,14 @@ module celestial_body_module
         !! A celestial body (Planet, moon, etc.)
         !! The `ID` from the [[base_class]] is the NAIF SPICE ID code for the body
         real(wp) :: mu = zero   !! gravitational parameter \( \mu \) [\(km^3/s^2\)]
+        !note: also should add radius, etc.
     end type celestial_body
 
+    type(celestial_body),parameter,public :: body_ssb = &
+                    celestial_body(0, 'SSB', 0.0_wp ) !! solar-system barycenter [note: don't have mu defined here yet]
+
     !define some bodies:
     ! MU values from: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/gm_de431.tpc
-
     type(celestial_body),parameter,public :: body_sun = &
                     celestial_body(10, 'Sun',1.3271244004193938E+11_wp )
     type(celestial_body),parameter,public :: body_mercury = &

src/fortran_astrodynamics_toolkit.f90

@@ -30,6 +30,7 @@ module fortran_astrodynamics_toolkit
     use kepler_module
     use kind_module
     use lambert_module
+    use lighting_module
     use math_module
     use matrix_module
     use minpack_module

src/jpl_ephemeris_module.f90

@@ -982,18 +982,23 @@ end subroutine get_constants
     subroutine ephemeris_test()
 
     use time_module, only: jd_to_et
-    use celestial_body_module
+    use celestial_body_module, fat_wp => wp
 
     implicit none
 
+    ! note: the low-level functions use real64 variables.
     character(len=6),dimension(nmax) :: nams
-    real(wp) :: jd, et
-    real(wp),dimension(6) :: rv,rv1,rv2,diffrv
-    real(wp),dimension(3) :: ss, r
+    real(wp),dimension(6) :: diffrv
+    real(wp),dimension(3) :: ss
     real(wp),dimension(nmax) :: vals
     integer :: nvs,ntarg,nctr,i,j
     type(jpl_ephemeris) :: eph405, eph421
     logical :: status_ok_405,status_ok_421
+    real(wp) :: jd_64, rv_64(6), rv1_64(6), rv2_64(6)
+    real(fat_wp) :: et
+    real(fat_wp),dimension(3) :: r
+    real(fat_wp),dimension(6) :: rv,rv1,rv2
+    real(fat_wp) :: jd
 
     character(len=*),parameter :: ephemeris_file_405 = './eph/JPLEPH.405' !! JPL DE405 ephemeris file
     character(len=*),parameter :: ephemeris_file_421 = './eph/JPLEPH.421' !! JPL DE421 ephemeris file
@@ -1021,7 +1026,7 @@ subroutine ephemeris_test()
             write(*,'(A,1X,D25.16)') nams(i), vals(i)
         end do
 
-        jd = 2451536.5d0   ! julian date
+        jd = 2451536.5_wp   ! julian date
         et = jd_to_et(jd)  ! ephemeris time
 
         if (jd < ss(1) .or. jd > ss(2)) then
@@ -1051,7 +1056,9 @@ subroutine ephemeris_test()
                             '" wrt "'//trim(list_of_bodies(nctr))//'"'
 
                 do i=1,10
-                    call eph405%get_state( jd, ntarg, nctr, rv, status_ok_405)
+                    jd_64 = jd
+                    call eph405%get_state( jd_64, ntarg, nctr, rv_64, status_ok_405)
+                    rv = rv_64
                     write(*,'(F15.2,1X,*(E25.16,1X))') jd, norm2(rv(1:3)), rv
                     jd = jd + 10.0_wp
                 end do
@@ -1082,8 +1089,12 @@ subroutine ephemeris_test()
                         '" wrt "'//trim(list_of_bodies(nctr))//'"'
 
             do i=1,10
-                call eph405%get_state( jd, ntarg, nctr, rv1, status_ok_405)
-                call eph421%get_state( jd, ntarg, nctr, rv2, status_ok_421)
+                jd_64 = jd
+                call eph405%get_state( jd_64, ntarg, nctr, rv1_64, status_ok_405)
+                rv1 = rv1_64
+                jd_64 = jd
+                call eph421%get_state( jd_64, ntarg, nctr, rv2_64, status_ok_421)
+                rv2 = rv2_64
                 diffrv = rv2 - rv1
                 write(*,'(F15.2,1X,*(E25.16,1X))') jd, norm2(diffrv(1:3)), norm2(diffrv(4:6))
                 jd = jd + 10.0_wp