From 0cf2854b94bdd7a620a8c113e76857445115c447 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Mon, 29 Sep 2025 12:56:29 -0500 Subject: [PATCH 01/16] prototype lighting routine --- src/fortran_astrodynamics_toolkit.f90 | 1 + src/lighting_module.f90 | 99 +++++++++++++++++++++++++++ 2 files changed, 100 insertions(+) create mode 100644 src/lighting_module.f90 diff --git a/src/fortran_astrodynamics_toolkit.f90 b/src/fortran_astrodynamics_toolkit.f90 index 5ef7effc..9ac5c837 100644 --- a/src/fortran_astrodynamics_toolkit.f90 +++ b/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 diff --git a/src/lighting_module.f90 b/src/lighting_module.f90 new file mode 100644 index 00000000..37d9113b --- /dev/null +++ b/src/lighting_module.f90 @@ -0,0 +1,99 @@ +!***************************************************************************************** +!> author: Jacob Williams +! +! Routines for computing solar fraction, lighting, eclipses, etc. + + module lighting_module + + use kind_module, only: wp + use numbers_module, only: pi, zero, one + use vector_module, only: unit + + implicit none + + private + + public :: solar_fraction + + contains +!***************************************************************************************** + +!***************************************************************************************** +!> +! Compute the solar fraction visible due to an eclipse by another body. +! +!### Reference +! * J. Wertz, "Spacecraft Attitude Determination and Control", 1978. +! See Chapter 3 and Appendix A. + + function solar_fraction(rp, rs, d_s, d_p) result(fraction) + + real(wp),intent(in) :: rp !! radius of the planet + real(wp),intent(in) :: rs !! radius of the Sun + real(wp),dimension(3),intent(in) :: d_s !! vector from the spacecraft to the Sun + real(wp),dimension(3),intent(in) :: d_p !! vector from the spacecraft to the planet + real(wp) :: fraction !! fraction of the Sun visible [0=total eclipse, 1=no eclipse] + + real(wp) :: s !! distance from the planet to the Sun + real(wp) :: c !! distance from the center of the planet to the apex of the shadow cone + real(wp) :: rho_c !! angular radius of the shadow cone + real(wp) :: rho_s !! angular radius of the Sun + real(wp) :: rho_p !! angular radius of the planet + real(wp) :: theta !! angular separation of the sun and planet as viewed by the spacecraft + real(wp) :: ds !! distance from the spacecraft to the Sun + real(wp) :: dp !! distance from the spacecraft to the planet + real(wp) :: drho !! difference in angular radii of the planet and Sun + real(wp) :: crp, crs, srp, srs, cth, sth, t1, t2, t3 !! temp variables + + if (rp<=zero) then ! no eclipse possible if the planet has no radius + fraction = one + return + end if + + ds = norm2(d_s) + dp = norm2(d_p) + + if (ds<=rs) then ! inside the Sun + fraction = one + return + else if (dp<=rp) then ! inside the planet + fraction = zero + return + end if + + s = norm2(d_s - d_p) + c = (rp*s) / (rs - rp) + rho_c = asin((rs - rp) / s) ! appx = asin(rs/s) + rho_s = asin(rs/ds) + rho_p = asin(rp/dp) + theta = acos(dot_product(unit(d_s), unit(d_p))) + drho = rho_p - rho_s + crp = cos(rho_p) + crs = cos(rho_s) + srp = sin(rho_p) + srs = sin(rho_s) + cth = cos(theta) + sth = sin(theta) + + if ( (ds>s) .and. (rho_p+rho_s>theta) .and. (theta>abs(drho)) ) then + ! partial eclipse + t1 = pi - crs * acos( (crp-crs*cth)/(srs*sth) ) + t2 = -crp * acos( (crs-crp*cth)/(srp*sth) ) + t3 = -acos( (cth-crs*crp)/(srs*srp) ) + fraction = (t1 + t2 + t3) / (pi*(one-crs)) + else if ( (stheta) ) then + ! total eclipse + fraction = 0.0_wp + else if ( (s+ctheta) ) then + ! annular eclipse + fraction = (one-crp) / (one-crs) + else + ! no eclipse + fraction = 1.0_wp + end if + + end function solar_fraction +!***************************************************************************************** + + end module lighting_module +!***************************************************************************************** \ No newline at end of file From 27bf1de1346fdc84e3c09077f2eab069bcce17b5 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Fri, 3 Oct 2025 21:49:46 -0500 Subject: [PATCH 02/16] high-level lighting routine also cubic model from halo library. untested! --- src/celestial_body_module.f90 | 5 +- src/lighting_module.f90 | 291 +++++++++++++++++++++++++++++++--- 2 files changed, 273 insertions(+), 23 deletions(-) diff --git a/src/celestial_body_module.f90 b/src/celestial_body_module.f90 index 187485ba..751dfd8f 100644 --- a/src/celestial_body_module.f90 +++ b/src/celestial_body_module.f90 @@ -17,9 +17,12 @@ module celestial_body_module real(wp) :: mu = zero !! gravitational parameter \( \mu \) [\(km^3/s^2\)] end type celestial_body + + type(celestial_body),parameter,public :: body_ssb = & + celestial_body(0, 'SSB', 0.0_wp ) !! solar-system barycenter + !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 = & diff --git a/src/lighting_module.f90 b/src/lighting_module.f90 index 37d9113b..aa6a72c1 100644 --- a/src/lighting_module.f90 +++ b/src/lighting_module.f90 @@ -5,34 +5,113 @@ module lighting_module - use kind_module, only: wp - use numbers_module, only: pi, zero, one - use vector_module, only: unit + use kind_module, only: wp + use numbers_module, only: pi, zero, one + use vector_module, only: unit, cross, axis_angle_rotation + use ephemeris_module, only: ephemeris_class + use transformation_module, only: icrf_frame + use celestial_body_module, only: celestial_body, body_sun, body_ssb + use conversion_module, only: deg2rad, rad2deg + use math_module, only: wrap_angle implicit none private - public :: solar_fraction + real(wp),parameter :: c = 299792.458_wp !! speed of light in km/s - contains + public :: from_j2000body_to_j2000ssb + public :: apparent_position + public :: get_sun_fraction ! high-level routine + public :: solar_fraction ! low-leve routine + public :: cubic_shadow_model ! low-level routine + + contains !***************************************************************************************** +!******************************************************************************** +!> +! Compute the "sun fraction" using the selected shadow model. + + function get_sun_fraction(b, rad_body, rad_sun, eph, et, rv, model, rbubble, use_geometric) result (phi) + + type(celestial_body),intent(in) :: b !! eclipsing body + real(wp),intent(in) :: rad_body !! radius of the eclipsing body [km] + real(wp),intent(in) :: rad_sun !! radius of the Sun [km] + class(ephemeris_class),intent(inout) :: eph !! the ephemeris to use for sun and ssb (if necessary) + real(wp),intent(in) :: et !! observer ephemeris time (sec) + real(wp),dimension(6),intent(in) :: rv !! state of the spacecraft (j2000-body frame) + integer,intent(in) :: model !! 1=circular cubic shadow model + !! 2=solar fraction model + real(wp),intent(in) :: rbubble !! eclipse bubble [km]. see the reference. + !! if rbubble=0, then no bubble is used. + !! only used if model=1 + logical,intent(in),optional :: use_geometric !! if true, use geometric positions + !! (no light time or stellar aberration correction) + !! default = false + real(wp) :: phi !! if `model=1`, circular cubic sun frac value: + !! + !! * >0 no eclipse, + !! * <0 eclipse, + !! * =0 on the eclipse line + !! + !! if `model=2`, true solar fraction value [0=total eclipse, 1=no eclipse], + !! with model of umbra/penumbra/antumbra (Wertz, 1978) + + logical :: status_ok !! true if no problems + real(wp),dimension(3) :: r_sun !! apparent state of the sun (j2000-ssb frame) + real(wp),dimension(3) :: r_body !! apparent state of the eclipsing body (j2000-ssb frame) + real(wp),dimension(6) :: rv_ssb !! state of the spacecraft !! (j2000-ssb frame) + logical :: use_apparent + + if (present(use_geometric)) then + use_apparent = .not. use_geometric + else + use_apparent = .true. + end if + + if (use_apparent) then + call from_j2000body_to_j2000ssb(b, eph, et, rv, rv_ssb) ! state of spacecraft in j2000-ssb + call apparent_position(eph, body_sun, et, rv_ssb, r_sun, status_ok) ! apparent position of sun in j2000 + if (.not. status_ok) error stop 'error getting apparent sun position' + call apparent_position(eph, b, et, rv_ssb, r_body, status_ok) ! apparent position of body in j2000 + if (.not. status_ok) error stop 'error getting apparent body position' + else + ! use geometric positions + r_body = -rv(1:3) ! geometric position of body wrt spacecraft in j2000 + call eph%get_r(et, body_sun, b, r_sun, status_ok) ! geometric position of sun wrt body in j2000 + if (.not. status_ok) error stop 'error getting geometric sun position' + r_sun = r_body + r_sun ! geometric position of sun wrt spacecraft in j2000 + end if + + ! compute sun fraction value + select case(model) + case(1) + call cubic_shadow_model(r_sun, rad_sun, r_body, rad_body, phi, rbubble) + case(2) + call solar_fraction(r_sun, rad_sun, r_body, rad_body, phi) + case default + error stop 'invalid sun fraction model' + end select + + end function get_sun_fraction +!******************************************************************************** + !***************************************************************************************** !> ! Compute the solar fraction visible due to an eclipse by another body. ! !### Reference -! * J. Wertz, "Spacecraft Attitude Determination and Control", 1978. +! * J. Wertz, "Spacecraft Attitude Determination and Control", 1978. ! See Chapter 3 and Appendix A. - function solar_fraction(rp, rs, d_s, d_p) result(fraction) + subroutine solar_fraction(d_s, rs, d_p, rp, fraction) - real(wp),intent(in) :: rp !! radius of the planet - real(wp),intent(in) :: rs !! radius of the Sun - real(wp),dimension(3),intent(in) :: d_s !! vector from the spacecraft to the Sun - real(wp),dimension(3),intent(in) :: d_p !! vector from the spacecraft to the planet - real(wp) :: fraction !! fraction of the Sun visible [0=total eclipse, 1=no eclipse] + real(wp),dimension(3),intent(in) :: d_s !! vector from the spacecraft to the Sun + real(wp),intent(in) :: rs !! radius of the Sun + real(wp),dimension(3),intent(in) :: d_p !! vector from the spacecraft to the planet + real(wp),intent(in) :: rp !! radius of the planet + real(wp),intent(out) :: fraction !! fraction of the Sun visible [0=total eclipse, 1=no eclipse] real(wp) :: s !! distance from the planet to the Sun real(wp) :: c !! distance from the center of the planet to the apex of the shadow cone @@ -50,7 +129,7 @@ function solar_fraction(rp, rs, d_s, d_p) result(fraction) return end if - ds = norm2(d_s) + ds = norm2(d_s) dp = norm2(d_p) if (ds<=rs) then ! inside the Sun @@ -80,20 +159,188 @@ function solar_fraction(rp, rs, d_s, d_p) result(fraction) t1 = pi - crs * acos( (crp-crs*cth)/(srs*sth) ) t2 = -crp * acos( (crs-crp*cth)/(srp*sth) ) t3 = -acos( (cth-crs*crp)/(srs*srp) ) - fraction = (t1 + t2 + t3) / (pi*(one-crs)) - else if ( (stheta) ) then - ! total eclipse - fraction = 0.0_wp - else if ( (s+ctheta) ) then + fraction = (t1 + t2 + t3) / (pi*(one-crs)) + else if ( (stheta) ) then + ! total eclipse + fraction = zero + else if ( (s+ctheta) ) then ! annular eclipse fraction = (one-crp) / (one-crs) - else - ! no eclipse - fraction = 1.0_wp + else + ! no eclipse + fraction = one end if - end function solar_fraction + end subroutine solar_fraction !***************************************************************************************** +!******************************************************************************** + subroutine from_j2000body_to_j2000ssb(b, eph, et, rv, rv_ssb) + + !! convert from a j2000-body frame to a j2000-ssb frame. + + type(celestial_body),intent(in) :: b !! eclipsing body + class(ephemeris_class),intent(inout) :: eph + real(wp),intent(in) :: et !! ephemeris time (sec) + real(wp),dimension(6),intent(in) :: rv !! j2000-body state (km, km/s) + real(wp),dimension(6),intent(out) :: rv_ssb !! j2000-ssb state (km, km/s) + + type(icrf_frame) :: f1, f2 + logical :: status_ok + + f1 = icrf_frame(b=b) + f2 = icrf_frame(b=body_ssb) + call f1%transform(rv,f2,et,eph,rv_ssb,status_ok) ! from f1 to f2 + if (.not. status_ok) error stop 'transformation error in from_j2000body_to_j2000ssb' + + end subroutine from_j2000body_to_j2000ssb +!******************************************************************************** + +!******************************************************************************** + subroutine apparent_position(eph, b_target, et, rv_obs_ssb, r_target, status_ok) + + !! Return the position of a target body relative to an observer, + !! corrected for light time and stellar aberration. + !! + !! see the SPICELIB routine `spkapo` (with 'lt+s') + + class(ephemeris_class),intent(inout) :: eph !! the ephemeris + type(celestial_body),intent(in) :: b_target !! target body + real(wp),dimension(6),intent(in) :: rv_obs_ssb !! state of the observer + !! (j2000 frame w.r.t. solar system barycenter) + real(wp),intent(in) :: et !! observer ephemeris time (sec) + real(wp),dimension(3),intent(out) :: r_target !! apparant state of the target (j2000 frame) + !! Corrected for one-way light time and stellar aberration + logical,intent(out) :: status_ok !! true if no problems + + real(wp),dimension(3) :: r_targ_ssb !! target body r wrt. ssb + real(wp),dimension(6) :: rv_targ_ssb !! target body rv wrt. ssb + real(wp) :: lt !! one-way light time [sec] + + ! Find the geometric position of the target body with respect to the + ! solar system barycenter. Subtract the position of the observer + ! to get the relative position. Use this to compute the one-way + ! light time. + call eph%get_r(et,b_target,body_ssb,r_targ_ssb,status_ok) + if (.not. status_ok) return + r_targ_ssb = r_targ_ssb - rv_obs_ssb(1:3) ! relative pos of target + lt = norm2(r_targ_ssb) / c ! light time + + ! To correct for light time, find the position of the target body + ! at the current epoch minus the one-way light time. Note that + ! the observer remains where he is. + call eph%get_r(et-lt,b_target,body_ssb,r_targ_ssb,status_ok) + if (.not. status_ok) return + r_targ_ssb = r_targ_ssb - rv_obs_ssb(1:3) + + ! At this point, r_targ_ssb contains the geometric or light-time + ! corrected position of the target relative to the observer + + ! stellar aberration correction + r_target = stellar_aberration(r_targ_ssb,rv_obs_ssb(4:6)) + + contains + + function stellar_aberration ( pobj, vobs ) result(appobj) + !! Correct the apparent position of an object for stellar aberration. + !! see SPICELIB routine `STELAB` + + real(wp),dimension(3),intent(in) :: pobj + real(wp),dimension(3),intent(in) :: vobs + real(wp),dimension(3) :: appobj + + real(wp),dimension(3) :: u, vbyc,h + real(wp) :: lensqr, sinphi, phi + real(wp),parameter :: zero_tol = epsilon(1.0_wp) !! tolerance for zero + + u = unit(pobj) + vbyc = vobs / c + lensqr = dot_product ( vbyc, vbyc ) + if ( lensqr >= 1.0_wp) error stop 'velocity > speed of light' + h = cross(u, vbyc) + sinphi = norm2 ( h ) + if ( abs(sinphi) > zero_tol ) then ! if (sinphi /= 0) + ! rotate the position of the object by phi + ! radians about h to obtain the apparent position. + phi = asin ( sinphi ) + call axis_angle_rotation ( pobj, h, phi, appobj ) + else + ! observer is moving along the line of sight to the object, + ! and no correction is required + appobj = pobj + end if + + end function stellar_aberration + + end subroutine apparent_position +!******************************************************************************** + +!******************************************************************************** +!> +! The "circular cubic" shadow model. +! +!### Reference +! * J. Williams, et. al, "A new eclipse algorithm for use in +! spacecraft trajectory optimization", 2023, AAS 23-243 + + subroutine cubic_shadow_model(rsun, radsun, rplanet, radplanet, sunfrac, rbubble) + + real(wp),dimension(3),intent(in) :: rsun !! apparent position vector of sun wrt spacecraft [km] + real(wp), intent(in) :: radsun !! radius of sun [km] + real(wp),dimension(3), intent(in) :: rplanet !! apparent position vector of eclipsing body wrt spacecraft [km] + real(wp), intent(in) :: radplanet !! radius of the eclipsing body [km] + real(wp), intent(out) :: sunfrac !! value of the function (>0 no eclipse, + !! <0 eclipse, =0 on the shadow line) + real(wp),intent(in),optional :: rbubble !! eclipse bubble radius. if present, then `sunfrac` is + !! the value along an arc length of `rbubble` + !! in the direction of the max eclipse line. + + real(wp),dimension(3) :: r !! radius vector from eclipsing body to spacecraft + real(wp),dimension(3) :: rsb !! radius vector from the sun to the eclipsing body + real(wp) :: tmp !! temp value + real(wp) :: alpha !! [deg] + real(wp) :: alpha0 !! [deg] + real(wp) :: sin_alpha0 !! `sin(alpha0)` + real(wp) :: rsbmag !! magnitude of radius vector from the sun to the eclipsing body + real(wp) :: rmag !! magnitude of `r` + logical :: compute_bubble !! use the `rbubble` inputs to adjust `alpha` + + compute_bubble = present(rbubble) + if (compute_bubble) compute_bubble = rbubble > zero + + r = -rplanet + rmag = norm2(r) + if (rmag Date: Fri, 3 Oct 2025 22:05:39 -0500 Subject: [PATCH 03/16] minor edits --- src/lighting_module.f90 | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/lighting_module.f90 b/src/lighting_module.f90 index aa6a72c1..acd90bf8 100644 --- a/src/lighting_module.f90 +++ b/src/lighting_module.f90 @@ -71,6 +71,7 @@ function get_sun_fraction(b, rad_body, rad_sun, eph, et, rv, model, rbubble, use end if if (use_apparent) then + ! apparent position of sun and body wrt to the spacecraft call from_j2000body_to_j2000ssb(b, eph, et, rv, rv_ssb) ! state of spacecraft in j2000-ssb call apparent_position(eph, body_sun, et, rv_ssb, r_sun, status_ok) ! apparent position of sun in j2000 if (.not. status_ok) error stop 'error getting apparent sun position' @@ -175,12 +176,13 @@ end subroutine solar_fraction !***************************************************************************************** !******************************************************************************** - subroutine from_j2000body_to_j2000ssb(b, eph, et, rv, rv_ssb) +!> +! convert from a j2000-body frame to a j2000-ssb frame. - !! convert from a j2000-body frame to a j2000-ssb frame. + subroutine from_j2000body_to_j2000ssb(b, eph, et, rv, rv_ssb) type(celestial_body),intent(in) :: b !! eclipsing body - class(ephemeris_class),intent(inout) :: eph + class(ephemeris_class),intent(inout) :: eph !! the ephemeris to use for body and ssb real(wp),intent(in) :: et !! ephemeris time (sec) real(wp),dimension(6),intent(in) :: rv !! j2000-body state (km, km/s) real(wp),dimension(6),intent(out) :: rv_ssb !! j2000-ssb state (km, km/s) @@ -197,12 +199,13 @@ end subroutine from_j2000body_to_j2000ssb !******************************************************************************** !******************************************************************************** - subroutine apparent_position(eph, b_target, et, rv_obs_ssb, r_target, status_ok) +!> +! Return the position of a target body relative to an observer, +! corrected for light time and stellar aberration. +! +! see the SPICELIB routine `spkapo` (with 'lt+s') - !! Return the position of a target body relative to an observer, - !! corrected for light time and stellar aberration. - !! - !! see the SPICELIB routine `spkapo` (with 'lt+s') + subroutine apparent_position(eph, b_target, et, rv_obs_ssb, r_target, status_ok) class(ephemeris_class),intent(inout) :: eph !! the ephemeris type(celestial_body),intent(in) :: b_target !! target body @@ -214,7 +217,6 @@ subroutine apparent_position(eph, b_target, et, rv_obs_ssb, r_target, status_ok) logical,intent(out) :: status_ok !! true if no problems real(wp),dimension(3) :: r_targ_ssb !! target body r wrt. ssb - real(wp),dimension(6) :: rv_targ_ssb !! target body rv wrt. ssb real(wp) :: lt !! one-way light time [sec] ! Find the geometric position of the target body with respect to the From 9a0c9851661763acb807329524e96d4c71b85b1b Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Fri, 3 Oct 2025 23:02:20 -0500 Subject: [PATCH 04/16] added another eclipse routine --- src/lighting_module.f90 | 90 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 79 insertions(+), 11 deletions(-) diff --git a/src/lighting_module.f90 b/src/lighting_module.f90 index acd90bf8..be183afe 100644 --- a/src/lighting_module.f90 +++ b/src/lighting_module.f90 @@ -6,7 +6,7 @@ module lighting_module use kind_module, only: wp - use numbers_module, only: pi, zero, one + use numbers_module, only: pi, zero, one, two use vector_module, only: unit, cross, axis_angle_rotation use ephemeris_module, only: ephemeris_class use transformation_module, only: icrf_frame @@ -18,12 +18,13 @@ module lighting_module private - real(wp),parameter :: c = 299792.458_wp !! speed of light in km/s + real(wp),parameter :: c_light = 299792.458_wp !! speed of light in km/s public :: from_j2000body_to_j2000ssb public :: apparent_position public :: get_sun_fraction ! high-level routine - public :: solar_fraction ! low-leve routine + public :: solar_fraction ! low-level routine + public :: solar_fraction_alt ! low-level routine public :: cubic_shadow_model ! low-level routine contains @@ -57,6 +58,7 @@ function get_sun_fraction(b, rad_body, rad_sun, eph, et, rv, model, rbubble, use !! !! if `model=2`, true solar fraction value [0=total eclipse, 1=no eclipse], !! with model of umbra/penumbra/antumbra (Wertz, 1978) + !! if `model=3`, alternate version of solar fraction (Montenbruck and Gill) logical :: status_ok !! true if no problems real(wp),dimension(3) :: r_sun !! apparent state of the sun (j2000-ssb frame) @@ -87,10 +89,9 @@ function get_sun_fraction(b, rad_body, rad_sun, eph, et, rv, model, rbubble, use ! compute sun fraction value select case(model) - case(1) - call cubic_shadow_model(r_sun, rad_sun, r_body, rad_body, phi, rbubble) - case(2) - call solar_fraction(r_sun, rad_sun, r_body, rad_body, phi) + case(1); call cubic_shadow_model(r_sun, rad_sun, r_body, rad_body, phi, rbubble) + case(2); call solar_fraction( r_sun, rad_sun, r_body, rad_body, phi) + case(3); call solar_fraction_alt(r_sun, rad_sun, r_body, rad_body, phi) case default error stop 'invalid sun fraction model' end select @@ -210,10 +211,10 @@ subroutine apparent_position(eph, b_target, et, rv_obs_ssb, r_target, status_ok) class(ephemeris_class),intent(inout) :: eph !! the ephemeris type(celestial_body),intent(in) :: b_target !! target body real(wp),dimension(6),intent(in) :: rv_obs_ssb !! state of the observer - !! (j2000 frame w.r.t. solar system barycenter) + !! (j2000 frame w.r.t. solar system barycenter) real(wp),intent(in) :: et !! observer ephemeris time (sec) real(wp),dimension(3),intent(out) :: r_target !! apparant state of the target (j2000 frame) - !! Corrected for one-way light time and stellar aberration + !! Corrected for one-way light time and stellar aberration logical,intent(out) :: status_ok !! true if no problems real(wp),dimension(3) :: r_targ_ssb !! target body r wrt. ssb @@ -226,7 +227,7 @@ subroutine apparent_position(eph, b_target, et, rv_obs_ssb, r_target, status_ok) call eph%get_r(et,b_target,body_ssb,r_targ_ssb,status_ok) if (.not. status_ok) return r_targ_ssb = r_targ_ssb - rv_obs_ssb(1:3) ! relative pos of target - lt = norm2(r_targ_ssb) / c ! light time + lt = norm2(r_targ_ssb) / c_light ! light time ! To correct for light time, find the position of the target body ! at the current epoch minus the one-way light time. Note that @@ -256,7 +257,7 @@ function stellar_aberration ( pobj, vobs ) result(appobj) real(wp),parameter :: zero_tol = epsilon(1.0_wp) !! tolerance for zero u = unit(pobj) - vbyc = vobs / c + vbyc = vobs / c_light lensqr = dot_product ( vbyc, vbyc ) if ( lensqr >= 1.0_wp) error stop 'velocity > speed of light' h = cross(u, vbyc) @@ -343,6 +344,73 @@ end function safe_acosd end subroutine cubic_shadow_model !***************************************************************************************** +!***************************************************************************************** +!> +! Another eclipse model, using flat plate-assumptions. +! +!### References +! * Montenbruck and Gill, "Satellite Orbits". +! * The GMAT routine `ShadowState::FindShadowState`. + + subroutine solar_fraction_alt(d_s, rs, d_p, rp, percentsun) + + real(wp),dimension(3),intent(in) :: d_s !! vector from the spacecraft to the Sun + real(wp),intent(in) :: rs !! radius of the Sun + real(wp),dimension(3),intent(in) :: d_p !! vector from the spacecraft to the planet + real(wp),intent(in) :: rp !! radius of the planet + real(wp),intent(out) :: percentsun !! fraction of the Sun visible [0=total eclipse, 1=no eclipse] + + real(wp),dimension(3) :: unitsun, d_s_hat, d_p_hat + real(wp) :: rho_s, rho_p, theta, rdotsun, d_s_mag, d_p_mag, c, a2, b2, x, y, area + + ! [sc] + ! / \ + ! d_s d_p + ! / \ + ! [sun] ---------- [body] + + unitsun = unit(d_s - d_p) ! body to sun unit vector + rdotsun = dot_product(-d_p,unitsun) + + if (rdotsun > zero) then ! sunny side of central body is always fully lit + percentsun = one + else + + d_s_mag = norm2(d_s) + d_p_mag = norm2(d_p) + if (rs >= d_s_mag) then ! inside the Sun + percentsun = one + else if (rp >= d_p_mag) then ! inside the planet + percentsun = zero + else + rho_s = asin(rs/d_s_mag) + rho_p = asin(rp/d_p_mag) + d_p_hat = unit(d_p) + d_s_hat = unit(d_s) + theta = acos(dot_product(d_p_hat,d_s_hat)) ! apparant distance from sun to body + + if (rho_s + rho_p <= theta) then ! full sunlight + percentsun = one + else if (theta <= rho_p-rho_s) then ! umbra + percentsun = zero + else if ( (abs(rho_s-rho_p) Date: Sat, 4 Oct 2025 16:22:45 -0500 Subject: [PATCH 05/16] fix typo and added tests --- src/lighting_module.f90 | 106 ++++++++++++++++++++++++++++++--- tests/unit_tests/fat_tests.f90 | 1 + 2 files changed, 100 insertions(+), 7 deletions(-) diff --git a/src/lighting_module.f90 b/src/lighting_module.f90 index be183afe..ad3dadb8 100644 --- a/src/lighting_module.f90 +++ b/src/lighting_module.f90 @@ -27,6 +27,8 @@ module lighting_module public :: solar_fraction_alt ! low-level routine public :: cubic_shadow_model ! low-level routine + public :: lighting_module_test + contains !***************************************************************************************** @@ -34,7 +36,7 @@ module lighting_module !> ! Compute the "sun fraction" using the selected shadow model. - function get_sun_fraction(b, rad_body, rad_sun, eph, et, rv, model, rbubble, use_geometric) result (phi) + function get_sun_fraction(b, rad_body, rad_sun, eph, et, rv, model, rbubble, use_geometric, info) result (phi) type(celestial_body),intent(in) :: b !! eclipsing body real(wp),intent(in) :: rad_body !! radius of the eclipsing body [km] @@ -50,6 +52,7 @@ function get_sun_fraction(b, rad_body, rad_sun, eph, et, rv, model, rbubble, use logical,intent(in),optional :: use_geometric !! if true, use geometric positions !! (no light time or stellar aberration correction) !! default = false + character(len=:),allocatable,intent(out),optional :: info !! info string real(wp) :: phi !! if `model=1`, circular cubic sun frac value: !! !! * >0 no eclipse, @@ -90,8 +93,8 @@ function get_sun_fraction(b, rad_body, rad_sun, eph, et, rv, model, rbubble, use ! compute sun fraction value select case(model) case(1); call cubic_shadow_model(r_sun, rad_sun, r_body, rad_body, phi, rbubble) - case(2); call solar_fraction( r_sun, rad_sun, r_body, rad_body, phi) - case(3); call solar_fraction_alt(r_sun, rad_sun, r_body, rad_body, phi) + case(2); call solar_fraction( r_sun, rad_sun, r_body, rad_body, phi, info) + case(3); call solar_fraction_alt(r_sun, rad_sun, r_body, rad_body, phi, info) case default error stop 'invalid sun fraction model' end select @@ -107,13 +110,14 @@ end function get_sun_fraction ! * J. Wertz, "Spacecraft Attitude Determination and Control", 1978. ! See Chapter 3 and Appendix A. - subroutine solar_fraction(d_s, rs, d_p, rp, fraction) + subroutine solar_fraction(d_s, rs, d_p, rp, fraction, info) real(wp),dimension(3),intent(in) :: d_s !! vector from the spacecraft to the Sun real(wp),intent(in) :: rs !! radius of the Sun real(wp),dimension(3),intent(in) :: d_p !! vector from the spacecraft to the planet real(wp),intent(in) :: rp !! radius of the planet real(wp),intent(out) :: fraction !! fraction of the Sun visible [0=total eclipse, 1=no eclipse] + character(len=:),allocatable,intent(out),optional :: info !! info string real(wp) :: s !! distance from the planet to the Sun real(wp) :: c !! distance from the center of the planet to the apex of the shadow cone @@ -127,6 +131,7 @@ subroutine solar_fraction(d_s, rs, d_p, rp, fraction) real(wp) :: crp, crs, srp, srs, cth, sth, t1, t2, t3 !! temp variables if (rp<=zero) then ! no eclipse possible if the planet has no radius + if (present(info)) info = 'no eclipse: planet radius <= 0' fraction = one return end if @@ -135,9 +140,11 @@ subroutine solar_fraction(d_s, rs, d_p, rp, fraction) dp = norm2(d_p) if (ds<=rs) then ! inside the Sun + if (present(info)) info = 'inside Sun' fraction = one return else if (dp<=rp) then ! inside the planet + if (present(info)) info = 'inside Planet' fraction = zero return end if @@ -158,18 +165,22 @@ subroutine solar_fraction(d_s, rs, d_p, rp, fraction) if ( (ds>s) .and. (rho_p+rho_s>theta) .and. (theta>abs(drho)) ) then ! partial eclipse + if (present(info)) info = 'partial eclipse' t1 = pi - crs * acos( (crp-crs*cth)/(srs*sth) ) t2 = -crp * acos( (crs-crp*cth)/(srp*sth) ) t3 = -acos( (cth-crs*crp)/(srs*srp) ) - fraction = (t1 + t2 + t3) / (pi*(one-crs)) + fraction = one - (t1 + t2 + t3) / (pi*(one-crs)) else if ( (stheta) ) then ! total eclipse + if (present(info)) info = 'total eclipse' fraction = zero else if ( (s+ctheta) ) then ! annular eclipse - fraction = (one-crp) / (one-crs) + if (present(info)) info = 'annular eclipse' + fraction = one - (one-crp) / (one-crs) else ! no eclipse + if (present(info)) info = 'no eclipse' fraction = one end if @@ -352,13 +363,14 @@ end subroutine cubic_shadow_model ! * Montenbruck and Gill, "Satellite Orbits". ! * The GMAT routine `ShadowState::FindShadowState`. - subroutine solar_fraction_alt(d_s, rs, d_p, rp, percentsun) + subroutine solar_fraction_alt(d_s, rs, d_p, rp, percentsun, info) real(wp),dimension(3),intent(in) :: d_s !! vector from the spacecraft to the Sun real(wp),intent(in) :: rs !! radius of the Sun real(wp),dimension(3),intent(in) :: d_p !! vector from the spacecraft to the planet real(wp),intent(in) :: rp !! radius of the planet real(wp),intent(out) :: percentsun !! fraction of the Sun visible [0=total eclipse, 1=no eclipse] + character(len=:),allocatable,intent(out),optional :: info !! info string real(wp),dimension(3) :: unitsun, d_s_hat, d_p_hat real(wp) :: rho_s, rho_p, theta, rdotsun, d_s_mag, d_p_mag, c, a2, b2, x, y, area @@ -373,14 +385,17 @@ subroutine solar_fraction_alt(d_s, rs, d_p, rp, percentsun) rdotsun = dot_product(-d_p,unitsun) if (rdotsun > zero) then ! sunny side of central body is always fully lit + if (present(info)) info = 'sunny side of body' percentsun = one else d_s_mag = norm2(d_s) d_p_mag = norm2(d_p) if (rs >= d_s_mag) then ! inside the Sun + if (present(info)) info = 'inside Sun' percentsun = one else if (rp >= d_p_mag) then ! inside the planet + if (present(info)) info = 'inside Planet' percentsun = zero else rho_s = asin(rs/d_s_mag) @@ -390,10 +405,13 @@ subroutine solar_fraction_alt(d_s, rs, d_p, rp, percentsun) theta = acos(dot_product(d_p_hat,d_s_hat)) ! apparant distance from sun to body if (rho_s + rho_p <= theta) then ! full sunlight + if (present(info)) info = 'full sunlight' percentsun = one else if (theta <= rho_p-rho_s) then ! umbra + if (present(info)) info = 'umbra' percentsun = zero else if ( (abs(rho_s-rho_p) +! Unit tests for the listing module. + + subroutine lighting_module_test() + + real(wp) :: rs, rp + real(wp),dimension(3) :: d_s, d_p + real(wp) :: phi1, phi2 + character(len=:),allocatable :: info1, info2 + Real(wp) :: RSun(3) !! Position vector of Sun + Real(wp) :: RPlanet(3) !! Position vector of planet + Real(wp) :: RSpcraft(3) !! Position vector of spacecraft + + rs = 1.0_wp + rp = 1.0_wp + + ! sun -- body -- sc -> 0.0 + d_s = [-100.0_wp, 0.0_wp, 0.0_wp] + d_p = [-10.0_wp, 0.0_wp, 0.0_wp] + call solar_fraction( d_s, rs, d_p, rp, phi1, info1) ! this one is wrong + call solar_fraction_alt(d_s, rs, d_p, rp, phi2, info2) + write(*,*) '' + write(*,*) 'phi1 = ', phi1, info1 + write(*,*) 'phi2 = ', phi2, info2 + write(*,*) 'diff = ', abs(phi1-phi2) + + ! sc -- sun -- body -> 1.0 + d_s = [10.0_wp, 0.0_wp, 0.0_wp] + d_p = [100.0_wp, 0.0_wp, 0.0_wp] + call solar_fraction( d_s, rs, d_p, rp, phi1, info1) + call solar_fraction_alt(d_s, rs, d_p, rp, phi2, info2) + write(*,*) '' + write(*,*) 'phi1 = ', phi1, info1 + write(*,*) 'phi2 = ', phi2, info2 + write(*,*) 'diff = ', abs(phi1-phi2) + + ! sc -- body -- sun -> 0.0 + d_s = [100.0_wp, 0.0_wp, 0.0_wp] + d_p = [10.0_wp, 0.0_wp, 0.0_wp] + call solar_fraction( d_s, rs, d_p, rp, phi1, info1) + call solar_fraction_alt(d_s, rs, d_p, rp, phi2, info2) + write(*,*) '' + write(*,*) 'phi1 = ', phi1, info1 + write(*,*) 'phi2 = ', phi2, info2 + write(*,*) 'diff = ', abs(phi1-phi2) + + ! sc -- body -- sun -> penumbra + d_s = [100.0_wp, 0.0_wp, 0.0_wp] + d_p = [10.0_wp, 1.0_wp, 0.0_wp] + call solar_fraction( d_s, rs, d_p, rp, phi1, info1) + call solar_fraction_alt(d_s, rs, d_p, rp, phi2, info2) + write(*,*) '' + write(*,*) 'phi1 = ', phi1, info1 + write(*,*) 'phi2 = ', phi2, info2 + write(*,*) 'diff = ', abs(phi1-phi2) + + ! realistic sun/earth case: + ! sun -- earth -- sc + rs = 696000.0_wp + rp = 6378.0_wp + d_s = [-149597870.7_wp, 0.0_wp, 0.0_wp] + d_p = [-6778.0_wp, 6400.0_wp, 0.0_wp] + call solar_fraction( d_s, rs, d_p, rp, phi1, info1) + call solar_fraction_alt(d_s, rs, d_p, rp, phi2, info2) + write(*,*) '' + write(*,*) 'phi1 = ', phi1, info1 + write(*,*) 'phi2 = ', phi2, info2 + write(*,*) 'diff = ', abs(phi1-phi2) + + end subroutine lighting_module_test +!***************************************************************************************** + !***************************************************************************************** end module lighting_module !***************************************************************************************** \ No newline at end of file diff --git a/tests/unit_tests/fat_tests.f90 b/tests/unit_tests/fat_tests.f90 index e141190d..70bba4d6 100644 --- a/tests/unit_tests/fat_tests.f90 +++ b/tests/unit_tests/fat_tests.f90 @@ -29,6 +29,7 @@ program fat_tests call standish_module_test() call halo_orbit_test() call eispack_test() + call lighting_module_test() end program fat_tests !***************************************************************************************** From f7c7a8e152b7620ac472c79ce17088f39900aaad Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sat, 4 Oct 2025 22:23:18 -0500 Subject: [PATCH 06/16] fixed issues in solar_fraction there was a typo in the original reference causing the antumbra case to be wrong. also protect for a divide by zero when the sun and body are the same size. updates to tests and info strings. --- src/lighting_module.f90 | 102 ++++++++++++++++++++++------------------ 1 file changed, 56 insertions(+), 46 deletions(-) diff --git a/src/lighting_module.f90 b/src/lighting_module.f90 index ad3dadb8..cbe1e2d6 100644 --- a/src/lighting_module.f90 +++ b/src/lighting_module.f90 @@ -108,7 +108,8 @@ end function get_sun_fraction ! !### Reference ! * J. Wertz, "Spacecraft Attitude Determination and Control", 1978. -! See Chapter 3 and Appendix A. +! See Chapter 3 and Appendix A. Note that the implementation here corrects +! a typo in this reference, and also protects for a division by zero. subroutine solar_fraction(d_s, rs, d_p, rp, fraction, info) @@ -128,7 +129,7 @@ subroutine solar_fraction(d_s, rs, d_p, rp, fraction, info) real(wp) :: ds !! distance from the spacecraft to the Sun real(wp) :: dp !! distance from the spacecraft to the planet real(wp) :: drho !! difference in angular radii of the planet and Sun - real(wp) :: crp, crs, srp, srs, cth, sth, t1, t2, t3 !! temp variables + real(wp) :: crp, crs, srp, srs, cth, sth, t1, t2, t3, delr !! temp variables if (rp<=zero) then ! no eclipse possible if the planet has no radius if (present(info)) info = 'no eclipse: planet radius <= 0' @@ -149,9 +150,16 @@ subroutine solar_fraction(d_s, rs, d_p, rp, fraction, info) return end if - s = norm2(d_s - d_p) - c = (rp*s) / (rs - rp) - rho_c = asin((rs - rp) / s) ! appx = asin(rs/s) + s = norm2(d_s - d_p) + delr = rs - rp + if (delr==zero) then + ! special case when the bodies are the same size, + ! to avoid division by zero + c = huge(1.0_wp) + else + c = (rp*s) / delr + end if + rho_c = asin(delr / s) ! appx = asin(rs/s) rho_s = asin(rs/ds) rho_p = asin(rp/dp) theta = acos(dot_product(unit(d_s), unit(d_p))) @@ -165,22 +173,22 @@ subroutine solar_fraction(d_s, rs, d_p, rp, fraction, info) if ( (ds>s) .and. (rho_p+rho_s>theta) .and. (theta>abs(drho)) ) then ! partial eclipse - if (present(info)) info = 'partial eclipse' + if (present(info)) info = 'penumbra' t1 = pi - crs * acos( (crp-crs*cth)/(srs*sth) ) t2 = -crp * acos( (crs-crp*cth)/(srp*sth) ) t3 = -acos( (cth-crs*crp)/(srs*srp) ) fraction = one - (t1 + t2 + t3) / (pi*(one-crs)) else if ( (stheta) ) then ! total eclipse - if (present(info)) info = 'total eclipse' + if (present(info)) info = 'umbra' fraction = zero - else if ( (s+ctheta) ) then + else if ( (s+c zero) then ! sunny side of central body is always fully lit - if (present(info)) info = 'sunny side of body' + if (present(info)) info = 'full sun' percentsun = one else @@ -440,52 +448,46 @@ subroutine lighting_module_test() real(wp),dimension(3) :: d_s, d_p real(wp) :: phi1, phi2 character(len=:),allocatable :: info1, info2 - Real(wp) :: RSun(3) !! Position vector of Sun - Real(wp) :: RPlanet(3) !! Position vector of planet - Real(wp) :: RSpcraft(3) !! Position vector of spacecraft - rs = 1.0_wp - rp = 1.0_wp + rs = 1.0_wp ! sun radius + rp = 1.0_wp ! planet radius ! sun -- body -- sc -> 0.0 d_s = [-100.0_wp, 0.0_wp, 0.0_wp] d_p = [-10.0_wp, 0.0_wp, 0.0_wp] - call solar_fraction( d_s, rs, d_p, rp, phi1, info1) ! this one is wrong - call solar_fraction_alt(d_s, rs, d_p, rp, phi2, info2) - write(*,*) '' - write(*,*) 'phi1 = ', phi1, info1 - write(*,*) 'phi2 = ', phi2, info2 - write(*,*) 'diff = ', abs(phi1-phi2) + call go() ! sc -- sun -- body -> 1.0 d_s = [10.0_wp, 0.0_wp, 0.0_wp] d_p = [100.0_wp, 0.0_wp, 0.0_wp] - call solar_fraction( d_s, rs, d_p, rp, phi1, info1) - call solar_fraction_alt(d_s, rs, d_p, rp, phi2, info2) - write(*,*) '' - write(*,*) 'phi1 = ', phi1, info1 - write(*,*) 'phi2 = ', phi2, info2 - write(*,*) 'diff = ', abs(phi1-phi2) + call go() ! sc -- body -- sun -> 0.0 d_s = [100.0_wp, 0.0_wp, 0.0_wp] d_p = [10.0_wp, 0.0_wp, 0.0_wp] - call solar_fraction( d_s, rs, d_p, rp, phi1, info1) - call solar_fraction_alt(d_s, rs, d_p, rp, phi2, info2) - write(*,*) '' - write(*,*) 'phi1 = ', phi1, info1 - write(*,*) 'phi2 = ', phi2, info2 - write(*,*) 'diff = ', abs(phi1-phi2) + call go() ! sc -- body -- sun -> penumbra d_s = [100.0_wp, 0.0_wp, 0.0_wp] d_p = [10.0_wp, 1.0_wp, 0.0_wp] - call solar_fraction( d_s, rs, d_p, rp, phi1, info1) - call solar_fraction_alt(d_s, rs, d_p, rp, phi2, info2) - write(*,*) '' - write(*,*) 'phi1 = ', phi1, info1 - write(*,*) 'phi2 = ', phi2, info2 - write(*,*) 'diff = ', abs(phi1-phi2) + call go() + + ! body -- sc -- sun + d_s = [-100.0_wp, 0.0_wp, 0.0_wp] + d_p = [100.0_wp, 0.0_wp, 0.0_wp] + call go() + + !.................................... + ! sc -- body -- sun -> antumbra + rs = 100.0_wp + d_s = [20000.0_wp, 0.0_wp, 0.0_wp] + d_p = [400.0_wp, 0.0_wp, 0.0_wp] + call go() + + rs = 100.0_wp ! umbra + d_s = [20000.0_wp, 0.0_wp, 0.0_wp] + d_p = [100.0_wp, 0.0_wp, 0.0_wp] + call go() ! realistic sun/earth case: ! sun -- earth -- sc @@ -493,13 +495,21 @@ subroutine lighting_module_test() rp = 6378.0_wp d_s = [-149597870.7_wp, 0.0_wp, 0.0_wp] d_p = [-6778.0_wp, 6400.0_wp, 0.0_wp] - call solar_fraction( d_s, rs, d_p, rp, phi1, info1) - call solar_fraction_alt(d_s, rs, d_p, rp, phi2, info2) - write(*,*) '' - write(*,*) 'phi1 = ', phi1, info1 - write(*,*) 'phi2 = ', phi2, info2 - write(*,*) 'diff = ', abs(phi1-phi2) + call go() + + contains + subroutine go() + print*, '----------------------------------' + call solar_fraction( d_s, rs, d_p, rp, phi1, info1) + call solar_fraction_alt(d_s, rs, d_p, rp, phi2, info2) + write(*,*) '' + write(*,*) 'phi1 = ', phi1, info1 + write(*,*) 'phi2 = ', phi2, info2 + write(*,*) 'diff = ', abs(phi1-phi2) + if (abs(phi1-phi2)>1.0e-4_wp) error stop 'WARNING: large difference between models' + print*, '' + end subroutine go end subroutine lighting_module_test !***************************************************************************************** From 04fa92193c4b99165e0427a6062af75379963ba1 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 5 Oct 2025 11:30:05 -0500 Subject: [PATCH 07/16] added another circular area method --- src/lighting_module.f90 | 114 ++++++++++++++++++++++++++++++++++++---- 1 file changed, 104 insertions(+), 10 deletions(-) diff --git a/src/lighting_module.f90 b/src/lighting_module.f90 index cbe1e2d6..73cbcf52 100644 --- a/src/lighting_module.f90 +++ b/src/lighting_module.f90 @@ -62,6 +62,7 @@ function get_sun_fraction(b, rad_body, rad_sun, eph, et, rv, model, rbubble, use !! if `model=2`, true solar fraction value [0=total eclipse, 1=no eclipse], !! with model of umbra/penumbra/antumbra (Wertz, 1978) !! if `model=3`, alternate version of solar fraction (Montenbruck and Gill) + !! if `model=4`, alternate version of solar fraction (nyxspace) logical :: status_ok !! true if no problems real(wp),dimension(3) :: r_sun !! apparent state of the sun (j2000-ssb frame) @@ -92,9 +93,10 @@ function get_sun_fraction(b, rad_body, rad_sun, eph, et, rv, model, rbubble, use ! compute sun fraction value select case(model) - case(1); call cubic_shadow_model(r_sun, rad_sun, r_body, rad_body, phi, rbubble) - case(2); call solar_fraction( r_sun, rad_sun, r_body, rad_body, phi, info) - case(3); call solar_fraction_alt(r_sun, rad_sun, r_body, rad_body, phi, info) + case(1); call cubic_shadow_model( r_sun, rad_sun, r_body, rad_body, phi, rbubble) + case(2); call solar_fraction( r_sun, rad_sun, r_body, rad_body, phi, info) + case(3); call solar_fraction_alt( r_sun, rad_sun, r_body, rad_body, phi, info) + case(4); call solar_fraction_alt2(r_sun, rad_sun, r_body, rad_body, phi) case default error stop 'invalid sun fraction model' end select @@ -365,7 +367,7 @@ end subroutine cubic_shadow_model !***************************************************************************************** !> -! Another eclipse model, using flat plate-assumptions. +! Another eclipse model, using circular area assumptions. ! !### References ! * Montenbruck and Gill, "Satellite Orbits". @@ -392,7 +394,9 @@ subroutine solar_fraction_alt(d_s, rs, d_p, rp, percentsun, info) unitsun = unit(d_s - d_p) ! body to sun unit vector rdotsun = dot_product(-d_p,unitsun) - if (rdotsun > zero) then ! sunny side of central body is always fully lit + if (rdotsun > zero) then + ! sunny side of central body is always fully lit + ! [the assumption here is the sun is always bigger than the body?] if (present(info)) info = 'full sun' percentsun = one else @@ -438,6 +442,85 @@ subroutine solar_fraction_alt(d_s, rs, d_p, rp, percentsun, info) end subroutine solar_fraction_alt !***************************************************************************************** +!***************************************************************************************** +!> +! Another eclipse model, using circular area assumptions, +! coded up based on the nixspace documentation. +! The results are very similar to `solar_fraction_alt`. +! +!### References +! * https://nyxspace.com/nyxspace/MathSpec/celestial/eclipse/#nomenclature + + subroutine solar_fraction_alt2(r_l, Rl, r_e, Re, percentsun, info) + + real(wp),dimension(3),intent(in) :: r_l !! vector from the spacecraft to the Sun + real(wp),intent(in) :: Rl !! radius of the Sun + real(wp),dimension(3),intent(in) :: r_e !! vector from the spacecraft to the planet + real(wp),intent(in) :: Re !! radius of the planet + real(wp),intent(out) :: percentsun !! fraction of the Sun visible [0=total eclipse, 1=no eclipse] + character(len=:),allocatable,intent(out),optional :: info !! info string + + real(wp) :: rlp, rep, dp, r_l_mag, r_e_mag, & + d1, d2, dp2, rlp2, rep2, At, Astar + + ! this check isn't mentioned in the reference, but needed + ! for sc -- sun -- body case + if (dot_product(-r_e,unit(r_l-r_e)) > zero) then + ! sunny side of body is always fully lit + ! [the assumption here is the sun is always bigger than the body?] + if (present(info)) info = 'full sun' + percentsun = one + return + end if + + r_l_mag = norm2(r_l) + r_e_mag = norm2(r_e) + + ! these checks also aren't in the writeup: + if (Rl >= r_l_mag) then ! inside the Sun + if (present(info)) info = 'inside Sun' + percentsun = one; return + else if (Re >= r_e_mag) then ! inside the planet + if (present(info)) info = 'inside Planet' + percentsun = zero; return + end if + + rlp = asin(Rl/r_l_mag) + rep = asin(Re/r_e_mag) + dp = acos(dot_product(r_l,r_e)/(r_l_mag*r_e_mag)) + + ! modified these two checks: + !if (dp-rlpdp+rlp) then ! original + else if (dp<=rep-rlp) then ! corrected + if (present(info)) info = 'umbra' + percentsun = zero ! umbra + else if (rlp-rep>=dp .or. dp>=rlp+rep) then ! antumbra + if (present(info)) info = 'antumbra' + percentsun = one - rep*rep/(rlp*rlp) + else ! penumbra + if (present(info)) info = 'penumbra' + dp2 = dp*dp + rlp2 = rlp*rlp + rep2 = rep*rep + d1 = (dp2 - rlp2 + rep2) / (two*dp) + d2 = (dp2 + rlp2 - rep2) / (two*dp) + At = A(rep, rep2, d1) + A(rlp, rlp2, d2) + Astar = pi*rlp2 + percentsun = (Astar - At) / Astar + end if + + contains + pure real(wp) function A(r,r2,d) + real(wp),intent(in) :: r, r2, d + A = r2*acos(d/r) - d*sqrt(r2 - d*d) + end function A + end subroutine solar_fraction_alt2 +!***************************************************************************************** + !***************************************************************************************** !> ! Unit tests for the listing module. @@ -446,8 +529,6 @@ subroutine lighting_module_test() real(wp) :: rs, rp real(wp),dimension(3) :: d_s, d_p - real(wp) :: phi1, phi2 - character(len=:),allocatable :: info1, info2 rs = 1.0_wp ! sun radius rp = 1.0_wp ! planet radius @@ -497,16 +578,29 @@ subroutine lighting_module_test() d_p = [-6778.0_wp, 6400.0_wp, 0.0_wp] call go() + ! ! an edge case, a very small sun very close to the body on x-axis, + ! ! sc on y-axis very close to body .. i don't think any properly handle this .. .double check... + ! rs = 0.0001_wp ! sun radius + ! rp = 10.0_wp ! planet radius + ! d_p = [0.0001_wp, -rp-0.01_wp, 0.0_wp] + ! d_s = d_p + [-rp-0.01_wp, 0.0_wp, 0.0_wp] + ! call go() + contains subroutine go() + real(wp) :: phi1, phi2, phi3 + character(len=:),allocatable :: info1, info2, info3 print*, '----------------------------------' - call solar_fraction( d_s, rs, d_p, rp, phi1, info1) - call solar_fraction_alt(d_s, rs, d_p, rp, phi2, info2) write(*,*) '' + call solar_fraction( d_s, rs, d_p, rp, phi1, info1) + call solar_fraction_alt( d_s, rs, d_p, rp, phi2, info2) + call solar_fraction_alt2(d_s, rs, d_p, rp, phi3, info3) write(*,*) 'phi1 = ', phi1, info1 write(*,*) 'phi2 = ', phi2, info2 - write(*,*) 'diff = ', abs(phi1-phi2) + write(*,*) 'phi3 = ', phi3, info3 + write(*,*) 'diff 1= ', abs(phi1-phi2) ! spherical vs circular + write(*,*) 'diff 2= ', abs(phi2-phi3) ! two circular models if (abs(phi1-phi2)>1.0e-4_wp) error stop 'WARNING: large difference between models' print*, '' end subroutine go From a94aa9efd497e4ebc2298a4a553d38cd5783a0b8 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 5 Oct 2025 11:38:55 -0500 Subject: [PATCH 08/16] this one not needed --- src/lighting_module.f90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/lighting_module.f90 b/src/lighting_module.f90 index 73cbcf52..a0b97a6d 100644 --- a/src/lighting_module.f90 +++ b/src/lighting_module.f90 @@ -489,13 +489,12 @@ subroutine solar_fraction_alt2(r_l, Rl, r_e, Re, percentsun, info) rep = asin(Re/r_e_mag) dp = acos(dot_product(r_l,r_e)/(r_l_mag*r_e_mag)) - ! modified these two checks: + ! modified this check: !if (dp-rlpdp+rlp) then ! original - else if (dp<=rep-rlp) then ! corrected + else if (rep>dp+rlp) then if (present(info)) info = 'umbra' percentsun = zero ! umbra else if (rlp-rep>=dp .or. dp>=rlp+rep) then ! antumbra From 332c270363de040ccb2a785bf81c63b0da62432a Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 5 Oct 2025 14:00:05 -0500 Subject: [PATCH 09/16] update --- src/lighting_module.f90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/lighting_module.f90 b/src/lighting_module.f90 index a0b97a6d..67975673 100644 --- a/src/lighting_module.f90 +++ b/src/lighting_module.f90 @@ -23,9 +23,10 @@ module lighting_module public :: from_j2000body_to_j2000ssb public :: apparent_position public :: get_sun_fraction ! high-level routine - public :: solar_fraction ! low-level routine - public :: solar_fraction_alt ! low-level routine - public :: cubic_shadow_model ! low-level routine + public :: solar_fraction ! low-level routine + public :: solar_fraction_alt ! low-level routine + public :: solar_fraction_alt2 ! low-level routine + public :: cubic_shadow_model ! low-level routine public :: lighting_module_test From c49bd99dfe44c78fa175278e9c93f7b69a7c4b75 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 5 Oct 2025 14:37:23 -0500 Subject: [PATCH 10/16] minor update --- src/lighting_module.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lighting_module.f90 b/src/lighting_module.f90 index 67975673..3f0d3372 100644 --- a/src/lighting_module.f90 +++ b/src/lighting_module.f90 @@ -495,7 +495,7 @@ subroutine solar_fraction_alt2(r_l, Rl, r_e, Re, percentsun, info) if (rlp+rep<=dp) then ! corrected if (present(info)) info = 'full sun' percentsun = one ! full sun - else if (rep>dp+rlp) then + else if (rep>=dp+rlp) then if (present(info)) info = 'umbra' percentsun = zero ! umbra else if (rlp-rep>=dp .or. dp>=rlp+rep) then ! antumbra From 942974732d924bc15db363707d05fac970962a1e Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 5 Oct 2025 17:33:11 -0500 Subject: [PATCH 11/16] fixed test compile for REAL128 numbers --- src/jpl_ephemeris_module.f90 | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/src/jpl_ephemeris_module.f90 b/src/jpl_ephemeris_module.f90 index 6e387a11..71e98e43 100644 --- a/src/jpl_ephemeris_module.f90 +++ b/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 From adbdbcebcc7cf35d306ec39163b047f8c5063c1c Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 5 Oct 2025 17:44:38 -0500 Subject: [PATCH 12/16] docstrings --- src/lighting_module.f90 | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/src/lighting_module.f90 b/src/lighting_module.f90 index 3f0d3372..63500102 100644 --- a/src/lighting_module.f90 +++ b/src/lighting_module.f90 @@ -45,8 +45,10 @@ function get_sun_fraction(b, rad_body, rad_sun, eph, et, rv, model, rbubble, use class(ephemeris_class),intent(inout) :: eph !! the ephemeris to use for sun and ssb (if necessary) real(wp),intent(in) :: et !! observer ephemeris time (sec) real(wp),dimension(6),intent(in) :: rv !! state of the spacecraft (j2000-body frame) - integer,intent(in) :: model !! 1=circular cubic shadow model - !! 2=solar fraction model + integer,intent(in) :: model !! algorithm to use: + !! + !! * 1: circular cubic shadow model + !! * 2-4: solar fraction model real(wp),intent(in) :: rbubble !! eclipse bubble [km]. see the reference. !! if rbubble=0, then no bubble is used. !! only used if model=1 @@ -54,16 +56,16 @@ function get_sun_fraction(b, rad_body, rad_sun, eph, et, rv, model, rbubble, use !! (no light time or stellar aberration correction) !! default = false character(len=:),allocatable,intent(out),optional :: info !! info string - real(wp) :: phi !! if `model=1`, circular cubic sun frac value: + real(wp) :: phi !! solar fraction returned: !! - !! * >0 no eclipse, - !! * <0 eclipse, - !! * =0 on the eclipse line - !! - !! if `model=2`, true solar fraction value [0=total eclipse, 1=no eclipse], - !! with model of umbra/penumbra/antumbra (Wertz, 1978) - !! if `model=3`, alternate version of solar fraction (Montenbruck and Gill) - !! if `model=4`, alternate version of solar fraction (nyxspace) + !! * if `model=1`, circular cubic sun frac value: + !! * `>0` no eclipse + !! * `<0` eclipse + !! * `=0` on the eclipse line + !! * if `model=2`, true solar fraction value [0=total eclipse, 1=no eclipse], + !! with model of umbra/penumbra/antumbra (Wertz, 1978) + !! * if `model=3`, alternate version of solar fraction (Montenbruck and Gill) + !! * if `model=4`, alternate version of solar fraction (nyxspace) logical :: status_ok !! true if no problems real(wp),dimension(3) :: r_sun !! apparent state of the sun (j2000-ssb frame) From 03f0e4bc0e98634ebe4f37eb1914e771da511b81 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 5 Oct 2025 22:21:49 -0500 Subject: [PATCH 13/16] minor change to avoid potential overflow if sun and body are the same size --- src/lighting_module.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/lighting_module.f90 b/src/lighting_module.f90 index 63500102..52017440 100644 --- a/src/lighting_module.f90 +++ b/src/lighting_module.f90 @@ -183,11 +183,11 @@ subroutine solar_fraction(d_s, rs, d_p, rp, fraction, info) t2 = -crp * acos( (crs-crp*cth)/(srp*sth) ) t3 = -acos( (cth-crs*crp)/(srs*srp) ) fraction = one - (t1 + t2 + t3) / (pi*(one-crs)) - else if ( (stheta) ) then + else if ( (stheta) ) then ! total eclipse if (present(info)) info = 'umbra' fraction = zero - else if ( (s+c Date: Sun, 5 Oct 2025 22:27:25 -0500 Subject: [PATCH 14/16] move a constant some formatting and notes --- src/celestial_body_module.f90 | 41 ++++++++++++----------------------- src/lighting_module.f90 | 4 +--- src/numbers_module.f90 | 2 ++ 3 files changed, 17 insertions(+), 30 deletions(-) diff --git a/src/celestial_body_module.f90 b/src/celestial_body_module.f90 index 751dfd8f..6f3424b6 100644 --- a/src/celestial_body_module.f90 +++ b/src/celestial_body_module.f90 @@ -15,38 +15,25 @@ 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 + 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 = & - celestial_body(199,'Mercury',2.2031780000000021E+04_wp ) - type(celestial_body),parameter,public :: body_venus = & - celestial_body(299,'Venus',3.2485859200000006E+05_wp ) - type(celestial_body),parameter,public :: body_earth = & - celestial_body(399,'Earth',3.9860043543609598E+05_wp ) - type(celestial_body),parameter,public :: body_earth_moon_barycenter = & - celestial_body(3,'Earth-Moon Barycenter',4.0350323550225981E+05_wp ) - type(celestial_body),parameter,public :: body_moon = & - celestial_body(301,'Moon',4.9028000661637961E+03_wp ) - type(celestial_body),parameter,public :: body_mars = & - celestial_body(499,'Mars',4.282837362069909E+04_wp ) - type(celestial_body),parameter,public :: body_jupiter = & - celestial_body(599,'Jupiter',1.266865349218008E+08_wp ) - type(celestial_body),parameter,public :: body_saturn = & - celestial_body(699,'Saturn',3.793120749865224E+07_wp ) - type(celestial_body),parameter,public :: body_uranus = & - celestial_body(799,'Uranus',5.793951322279009E+06_wp ) - type(celestial_body),parameter,public :: body_neptune = & - celestial_body(899,'Neptune',6.835099502439672E+06_wp ) - type(celestial_body),parameter,public :: body_pluto = & - celestial_body(999,'Pluto',8.696138177608748E+02_wp ) + type(celestial_body),parameter,public :: body_sun = celestial_body(10, 'Sun', 1.3271244004193938E+11_wp) + type(celestial_body),parameter,public :: body_mercury = celestial_body(199,'Mercury', 2.2031780000000021E+04_wp) + type(celestial_body),parameter,public :: body_venus = celestial_body(299,'Venus', 3.2485859200000006E+05_wp) + type(celestial_body),parameter,public :: body_earth = celestial_body(399,'Earth', 3.9860043543609598E+05_wp) + type(celestial_body),parameter,public :: body_earth_moon_barycenter = celestial_body(3, 'Earth-Moon Barycenter',4.0350323550225981E+05_wp) + type(celestial_body),parameter,public :: body_moon = celestial_body(301,'Moon', 4.9028000661637961E+03_wp) + type(celestial_body),parameter,public :: body_mars = celestial_body(499,'Mars', 4.282837362069909E+04_wp ) + type(celestial_body),parameter,public :: body_jupiter = celestial_body(599,'Jupiter', 1.266865349218008E+08_wp ) + type(celestial_body),parameter,public :: body_saturn = celestial_body(699,'Saturn', 3.793120749865224E+07_wp ) + type(celestial_body),parameter,public :: body_uranus = celestial_body(799,'Uranus', 5.793951322279009E+06_wp ) + type(celestial_body),parameter,public :: body_neptune = celestial_body(899,'Neptune', 6.835099502439672E+06_wp ) + type(celestial_body),parameter,public :: body_pluto = celestial_body(999,'Pluto', 8.696138177608748E+02_wp ) !***************************************************************************************** end module celestial_body_module diff --git a/src/lighting_module.f90 b/src/lighting_module.f90 index 52017440..68b354ae 100644 --- a/src/lighting_module.f90 +++ b/src/lighting_module.f90 @@ -6,7 +6,7 @@ module lighting_module use kind_module, only: wp - use numbers_module, only: pi, zero, one, two + use numbers_module, only: pi, zero, one, two, c_light use vector_module, only: unit, cross, axis_angle_rotation use ephemeris_module, only: ephemeris_class use transformation_module, only: icrf_frame @@ -18,8 +18,6 @@ module lighting_module private - real(wp),parameter :: c_light = 299792.458_wp !! speed of light in km/s - public :: from_j2000body_to_j2000ssb public :: apparent_position public :: get_sun_fraction ! high-level routine diff --git a/src/numbers_module.f90 b/src/numbers_module.f90 index fc0a81d7..2f11306c 100644 --- a/src/numbers_module.f90 +++ b/src/numbers_module.f90 @@ -29,6 +29,8 @@ module numbers_module real(wp),parameter,public :: universal_grav_constant = 6.67408e-20_wp !! CODATA-recommended universal gravitational !! constant \( km^3/kg-s^2 \) + real(wp),parameter,public :: c_light = 299792.458_wp !! speed of light in km/s + !> 3x3 identity matrix: real(wp),dimension(3,3),parameter,public :: identity_3x3 = reshape(& [[one,zero,zero],& From d490e6a0678a9a79e674616ab0e4234d6877c0d4 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 5 Oct 2025 22:28:10 -0500 Subject: [PATCH 15/16] update license --- LICENSE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LICENSE b/LICENSE index 7044ac49..ef9c0036 100644 --- a/LICENSE +++ b/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, From 7497d255f2e0f383cecf2ba2acf0f01f441113c5 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 5 Oct 2025 22:34:11 -0500 Subject: [PATCH 16/16] avoid line truncation issues --- src/celestial_body_module.f90 | 39 +++++++++++++++++++++++------------ 1 file changed, 26 insertions(+), 13 deletions(-) diff --git a/src/celestial_body_module.f90 b/src/celestial_body_module.f90 index 6f3424b6..b6bbc430 100644 --- a/src/celestial_body_module.f90 +++ b/src/celestial_body_module.f90 @@ -18,22 +18,35 @@ module celestial_body_module !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] + 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 = celestial_body(199,'Mercury', 2.2031780000000021E+04_wp) - type(celestial_body),parameter,public :: body_venus = celestial_body(299,'Venus', 3.2485859200000006E+05_wp) - type(celestial_body),parameter,public :: body_earth = celestial_body(399,'Earth', 3.9860043543609598E+05_wp) - type(celestial_body),parameter,public :: body_earth_moon_barycenter = celestial_body(3, 'Earth-Moon Barycenter',4.0350323550225981E+05_wp) - type(celestial_body),parameter,public :: body_moon = celestial_body(301,'Moon', 4.9028000661637961E+03_wp) - type(celestial_body),parameter,public :: body_mars = celestial_body(499,'Mars', 4.282837362069909E+04_wp ) - type(celestial_body),parameter,public :: body_jupiter = celestial_body(599,'Jupiter', 1.266865349218008E+08_wp ) - type(celestial_body),parameter,public :: body_saturn = celestial_body(699,'Saturn', 3.793120749865224E+07_wp ) - type(celestial_body),parameter,public :: body_uranus = celestial_body(799,'Uranus', 5.793951322279009E+06_wp ) - type(celestial_body),parameter,public :: body_neptune = celestial_body(899,'Neptune', 6.835099502439672E+06_wp ) - type(celestial_body),parameter,public :: body_pluto = celestial_body(999,'Pluto', 8.696138177608748E+02_wp ) + type(celestial_body),parameter,public :: body_sun = & + celestial_body(10, 'Sun',1.3271244004193938E+11_wp ) + type(celestial_body),parameter,public :: body_mercury = & + celestial_body(199,'Mercury',2.2031780000000021E+04_wp ) + type(celestial_body),parameter,public :: body_venus = & + celestial_body(299,'Venus',3.2485859200000006E+05_wp ) + type(celestial_body),parameter,public :: body_earth = & + celestial_body(399,'Earth',3.9860043543609598E+05_wp ) + type(celestial_body),parameter,public :: body_earth_moon_barycenter = & + celestial_body(3,'Earth-Moon Barycenter',4.0350323550225981E+05_wp ) + type(celestial_body),parameter,public :: body_moon = & + celestial_body(301,'Moon',4.9028000661637961E+03_wp ) + type(celestial_body),parameter,public :: body_mars = & + celestial_body(499,'Mars',4.282837362069909E+04_wp ) + type(celestial_body),parameter,public :: body_jupiter = & + celestial_body(599,'Jupiter',1.266865349218008E+08_wp ) + type(celestial_body),parameter,public :: body_saturn = & + celestial_body(699,'Saturn',3.793120749865224E+07_wp ) + type(celestial_body),parameter,public :: body_uranus = & + celestial_body(799,'Uranus',5.793951322279009E+06_wp ) + type(celestial_body),parameter,public :: body_neptune = & + celestial_body(899,'Neptune',6.835099502439672E+06_wp ) + type(celestial_body),parameter,public :: body_pluto = & + celestial_body(999,'Pluto',8.696138177608748E+02_wp ) !***************************************************************************************** end module celestial_body_module