diff --git a/fpm.toml b/fpm.toml index 61da785..7d2a9c5 100644 --- a/fpm.toml +++ b/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 diff --git a/src/c_interface_module.f90 b/src/c_interface_module.f90 index 7bea827..adc4ed9 100644 --- a/src/c_interface_module.f90 +++ b/src/c_interface_module.f90 @@ -9,6 +9,7 @@ module c_interface_module use iso_c_binding use geopotential_module + use kind_module, only: wp implicit none @@ -142,13 +143,19 @@ subroutine get_acceleration(cp,n,m,rvec,acc) bind(c,name='get_acceleration') 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' diff --git a/src/geopotential_module.f90 b/src/geopotential_module.f90 index e70f8d2..7ac3143 100644 --- a/src/geopotential_module.f90 +++ b/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] diff --git a/src/halo_orbit_module.f90 b/src/halo_orbit_module.f90 index 5a5d657..9233a7e 100755 --- a/src/halo_orbit_module.f90 +++ b/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 diff --git a/src/jpl_ephemeris_module.f90 b/src/jpl_ephemeris_module.f90 index 335b4cb..e85a253 100644 --- a/src/jpl_ephemeris_module.f90 +++ b/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(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,',& diff --git a/src/kind_module.F90 b/src/kind_module.F90 new file mode 100644 index 0000000..f33a5ec --- /dev/null +++ b/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 +!***************************************************************************************** diff --git a/src/kind_module.f90 b/src/kind_module.f90 deleted file mode 100644 index 8bc86a5..0000000 --- a/src/kind_module.f90 +++ /dev/null @@ -1,19 +0,0 @@ -!***************************************************************************************** -!> author: Jacob Williams -! -! Define the numeric kinds. - - module kind_module - - use, intrinsic :: iso_fortran_env, only: real32,real64,real128 - - implicit none - - private - - !integer,parameter,public :: wp = real32 !! single precision reals - integer,parameter,public :: wp = real64 !! double precision reals - !integer,parameter,public :: wp = real128 !! quad precision reals - - end module kind_module -!***************************************************************************************** diff --git a/src/minpack_module.f90 b/src/minpack_module.f90 index c56c65a..42eb51b 100644 --- a/src/minpack_module.f90 +++ b/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 diff --git a/tests/ephemeris/ephemeris_comparison.f90 b/tests/ephemeris/ephemeris_comparison.f90 index 168590a..1c332cb 100644 --- a/tests/ephemeris/ephemeris_comparison.f90 +++ b/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) diff --git a/tests/lambert/lambert_test.f90 b/tests/lambert/lambert_test.f90 index 51644f8..d97f3b0 100644 --- a/tests/lambert/lambert_test.f90 +++ b/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 diff --git a/tests/lambert/lambert_test_2.f90 b/tests/lambert/lambert_test_2.f90 index 8725026..fee5743 100644 --- a/tests/lambert/lambert_test_2.f90 +++ b/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