solar_fraction_alt2 Subroutine

public subroutine solar_fraction_alt2(r_l, Rl, r_e, Re, percentsun, info)

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

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in), dimension(3) :: r_l

vector from the spacecraft to the Sun
real(kind=wp), intent(in) :: Rl

radius of the Sun
real(kind=wp), intent(in), dimension(3) :: r_e

vector from the spacecraft to the planet
real(kind=wp), intent(in) :: Re

radius of the planet
real(kind=wp), intent(out) :: percentsun

fraction of the Sun visible [0=total eclipse, 1=no eclipse]
character(len=:), intent(out), optional, allocatable :: info

info string

Calls

proc~~solar_fraction_alt2~~CallsGraph proc~solar_fraction_alt2 solar_fraction_alt2 proc~unit unit proc~solar_fraction_alt2->proc~unit

Called by

proc~~solar_fraction_alt2~~CalledByGraph proc~solar_fraction_alt2 solar_fraction_alt2 proc~get_sun_fraction get_sun_fraction proc~get_sun_fraction->proc~solar_fraction_alt2 proc~lighting_module_test lighting_module_test proc~lighting_module_test->proc~solar_fraction_alt2

Source Code

    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<rep) then  ! original
    if (rlp+rep<=dp) then  ! corrected
        if (present(info)) info = 'full sun'
        percentsun = one ! full sun
    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
        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