specify real kind via preprocessor directive

fpm.toml

@@ -7,7 +7,7 @@ homepage = "https://github.com/jacobwilliams/fortran-astrodynamics-toolkit"
 keywords = ["astrodynamics"]
 
 [dev-dependencies]
-pyplot-fortran= { git = "https://github.com/jacobwilliams/pyplot-fortran", rev = "3.1.0" }
+pyplot-fortran= { git = "https://github.com/jacobwilliams/pyplot-fortran", rev = "3.3.0" }
 
 [build]
 auto-executables = false

src/c_interface_module.f90

@@ -9,6 +9,7 @@
 
     use iso_c_binding
     use geopotential_module
+    use kind_module, only: wp
 
     implicit none
 
@@ -142,13 +143,19 @@
 
     type(container),pointer :: grav_container  !! Fortran version of `cp`
 
+    ! just in case wp /= c_double, we have to make a copy here
+    real(wp),dimension(3) :: rvec_f !! position vector
+    real(wp),dimension(3) :: acc_f !! acceleration vector
+
     ! convert cp to fortran:
     call c_f_pointer(cp,grav_container)
 
     if (associated(grav_container)) then
         select type (g => grav_container%data)
         class is (geopotential_model)
-            call g%get_acc(rvec,n,m,acc)
+            rvec_f = rvec
+            call g%get_acc(rvec_f,n,m,acc_f)
+            acc = acc_f
         end select
     else
         error stop 'error: pointer is not associated'

src/geopotential_module.f90

@@ -1133,8 +1133,6 @@ end subroutine kuga_carrara_geopotential
 
     pure function pinesnorm(mu,req,r_f,cnm,snm,nmax,mmax) result(accel)
 
-    use iso_fortran_env, only: wp => real64
-
     implicit none
 
     real(wp),intent(in)                 :: mu     !! gravitational constant [km^3/s^2]

src/halo_orbit_module.f90

@@ -10,7 +10,8 @@
 
     module halo_orbit_module
 
-    use iso_fortran_env, only: wp => real64, error_unit
+    use kind_module
+    use iso_fortran_env, only: error_unit
     use numbers_module
     use crtbp_module
 

src/jpl_ephemeris_module.f90

@@ -44,10 +44,15 @@
 !### History
 !  * [Original code from JPL](ftp://ssd.jpl.nasa.gov/pub/eph/planets/fortran/) Version : March 25, 2013
 !  * Extensive modifications by Jacob Williams for the Fortran Astrodynamics Toolkit.
+!
+!@note Warning: all calculations here are done with 64 bit reals. if using 128-bit reals
+!      the state is just copied into a 128-bit variable. Perhaps we should do the interpolations
+!      in 128-bit? (the data from the file must be read as 64-bit)
 
     module jpl_ephemeris_module
 
-    use, intrinsic :: iso_fortran_env, only: real64,error_unit
+    use kind_module, only: fat_wp => wp
+    use, intrinsic :: iso_fortran_env, only: real64, error_unit
     use ephemeris_module
 
     implicit none
@@ -177,16 +182,17 @@ subroutine get_rv_from_jpl_ephemeris(me,et,targ,obs,rv,status_ok)
 
     implicit none
 
-    class(jpl_ephemeris),intent(inout) :: me
-    real(wp),intent(in)                :: et         !! ephemeris time [sec]
-    type(celestial_body),intent(in)    :: targ       !! target body
-    type(celestial_body),intent(in)    :: obs        !! observer body
-    real(wp),dimension(6),intent(out)  :: rv         !! state of targ w.r.t. obs [km,km/s] in ICRF frame
-    logical,intent(out)                :: status_ok  !! true if there were no problems
+    class(jpl_ephemeris),intent(inout)    :: me
+    real(fat_wp),intent(in)               :: et         !! ephemeris time [sec]
+    type(celestial_body),intent(in)       :: targ       !! target body
+    type(celestial_body),intent(in)       :: obs        !! observer body
+    real(fat_wp),dimension(6),intent(out) :: rv         !! state of targ w.r.t. obs [km,km/s] in ICRF frame
+    logical,intent(out)                   :: status_ok  !! true if there were no problems
 
     real(wp) :: jd     !! julian date for input to [[get_state]].
     integer :: ntarg   !! id code for target body
     integer :: ncent   !! id code for observer body
+    real(wp),dimension(6) :: rv_ !! in case `wp /= fat_wp` we need a copy
 
     if (targ==obs) then
         !don't bother if target and observer are the same body
@@ -200,7 +206,8 @@ subroutine get_rv_from_jpl_ephemeris(me,et,targ,obs,rv,status_ok)
         ncent = spice_id_to_old_id(obs%id)
 
         if (ntarg>0 .and. ncent>0) then
-            call me%get_state(jd,ntarg,ncent,rv,status_ok)
+            call me%get_state(jd,ntarg,ncent,rv_,status_ok)
+            rv = rv_
             if (status_ok) then
                 if (.not. me%km) then
                     !we must return in units of km/s
@@ -776,7 +783,7 @@ subroutine state(me,et2,list,pv,pnut,status_ok)
     ! error return for epoch out of range
 
     if (pjd(1)+pjd(4)<me%ss(1) .or. pjd(1)+pjd(4)>me%ss(2)) then
-        write(error_unit,'(A,F12.2,A,2F12.2)') &
+        write(error_unit,'(A,F12.2,A,2F22.2)') &
                 'Error: requested jed,',&
                 et2(1)+et2(2),&
                 ' not within ephemeris limits,',&

src/kind_module.F90

@@ -0,0 +1,25 @@
+!*****************************************************************************************
+!> author: Jacob Williams
+!
+!  Define the numeric kinds.
+
+    module kind_module
+
+    use, intrinsic :: iso_fortran_env,    only: real32,real64,real128
+
+    implicit none
+
+    private
+
+#ifdef REAL32
+    integer,parameter,public :: wp = real32   !! real kind used by this module [4 bytes]
+#elif REAL64
+    integer,parameter,public :: wp = real64   !! real kind used by this module [8 bytes]
+#elif REAL128
+    integer,parameter,public :: wp = real128  !! real kind used by this module [16 bytes]
+#else
+    integer,parameter,public :: wp = real64   !! real kind used by this module [8 bytes]
+#endif
+
+    end module kind_module
+!*****************************************************************************************

src/minpack_module.f90

@@ -817,7 +817,7 @@ subroutine hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,diag,mode, &
           ncfail = 0
           ncsuc = ncsuc + 1
           if ( ratio>=p5 .or. ncsuc>1 ) delta = dmax1(delta,pnorm/p5)
-          if ( dabs(ratio-one)<=p1 ) delta = pnorm/p5
+          if ( abs(ratio-one)<=p1 ) delta = pnorm/p5
        else
           ncsuc = 0
           ncfail = ncfail + 1

tests/ephemeris/ephemeris_comparison.f90

@@ -5,12 +5,13 @@
     program ephemeris_comparison
 
     use fortran_astrodynamics_toolkit, wp => fat_wp
+    use iso_fortran_env, only: real64
 
     implicit none
 
     real(wp)              :: jd
     real(wp)              :: jd0
-    real(wp),dimension(6) :: rv
+    real(real64),dimension(6) :: rv
     real(wp),dimension(3) :: r_s,v_s
     real(wp)              :: max_err_s,err_s
     integer               :: ntarg,nctr
@@ -43,7 +44,7 @@ program ephemeris_comparison
 
             jd = real(i,wp)
 
-            call eph405%get_state(jd, ntarg, nctr, rv, status_ok)
+            call eph405%get_state(real(jd,real64), ntarg, nctr, rv, status_ok)
             if (.not.status_ok) error stop 'error computing state.'
 
             call simpson_lunar_ephemeris(jd, r_s, v_s)

tests/lambert/lambert_test.f90

@@ -27,9 +27,9 @@ program lambert_test_1
     real(wp),parameter :: tolerance = 1.0e-9_wp
     integer,parameter :: max_iterations = 1000
 
-    r1 = [20.0d3, 20.0d3, zero]
-    r2 = [-20.0d3, 10.0d3, zero]
-    tof = 1.0d0 * day2sec
+    r1 = [20.0e3_wp, 20.0e3_wp, zero]
+    r2 = [-20.0e3_wp, 10.0e3_wp, zero]
+    tof = 1.0_wp * day2sec
     mu = 398600.0_wp
 
     do icase = 1, 2 ! Gooding, Izzo

tests/lambert/lambert_test_2.f90

@@ -26,8 +26,8 @@ program lambert_test_2
     open(newunit=iunit_short,file='test2_short.csv',status='REPLACE')
     open(newunit=iunit_long, file='test2_long.csv', status='REPLACE')
 
-    r1 = [20.0d3, 20.0d3, zero]  ! initial position
-    r2 = [-20.0d3, 10.0d3, zero] ! final position
+    r1 = [20.0e3_wp, 20.0e3_wp, zero]  ! initial position
+    r2 = [-20.0e3_wp, 10.0e3_wp, zero] ! final position
     mu = 398600.0_wp             ! earth
     multi_revs = 5               ! try up to 5 revs
     n = 2*day2sec                ! 2 days