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, diff --git a/src/celestial_body_module.f90 b/src/celestial_body_module.f90 index 187485ba..b6bbc430 100644 --- a/src/celestial_body_module.f90 +++ b/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 = & 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/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 diff --git a/src/lighting_module.f90 b/src/lighting_module.f90 new file mode 100644 index 00000000..68b354ae --- /dev/null +++ b/src/lighting_module.f90 @@ -0,0 +1,612 @@ +!***************************************************************************************** +!> 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, two, c_light + 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 :: 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 :: solar_fraction_alt2 ! low-level routine + public :: cubic_shadow_model ! low-level routine + + public :: lighting_module_test + + 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, info) 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 !! 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 + 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 !! solar fraction returned: + !! + !! * 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) + 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 + ! 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' + 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, 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 + + 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. +! 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) + + 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 + 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, 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' + fraction = one + return + end if + + ds = norm2(d_s) + 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 + + 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))) + 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 + 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 = 'umbra' + fraction = zero + else if ( (c +! 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 !! 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) + + 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 +!******************************************************************************** + +!******************************************************************************** +!> +! 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 + 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) :: 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 ! 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_light + 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 +! Another eclipse model, using circular area assumptions. +! +!### References +! * Montenbruck and Gill, "Satellite Orbits". +! * The GMAT routine `ShadowState::FindShadowState`. + + 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 + + ! [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 + ! [the assumption here is the sun is always bigger than the body?] + if (present(info)) info = 'full sun' + 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) + 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 + 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) +! 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 this check: + !if (dp-rlp=dp+rlp) then + 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. + + subroutine lighting_module_test() + + real(wp) :: rs, rp + real(wp),dimension(3) :: d_s, d_p + + 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 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 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 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 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 + 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 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*, '----------------------------------' + 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(*,*) '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 + end subroutine lighting_module_test +!***************************************************************************************** + +!***************************************************************************************** + end module lighting_module +!***************************************************************************************** \ No newline at end of file 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],& 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 !*****************************************************************************************