solar_fraction Subroutine

public subroutine solar_fraction(d_s, rs, d_p, rp, fraction, info)

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.

Arguments

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

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

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

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

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

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

info string

Calls

proc~~solar_fraction~~CallsGraph proc~solar_fraction solar_fraction proc~unit unit proc~solar_fraction->proc~unit

Called by

proc~~solar_fraction~~CalledByGraph proc~solar_fraction solar_fraction proc~get_sun_fraction get_sun_fraction proc~get_sun_fraction->proc~solar_fraction proc~lighting_module_test lighting_module_test proc~lighting_module_test->proc~solar_fraction

Source Code

    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 ( (s<ds) .and. (ds-s<c) .and. (drho>theta) ) then
        ! total eclipse
        if (present(info)) info = 'umbra'
        fraction = zero
    else if ( (c<ds-s) .and. (drho<theta) ) then   ! JW : typo in original reference
        ! annular eclipse
        if (present(info)) info = 'antumbra'
        fraction = one - (one-crp) / (one-crs)
    else
        ! no eclipse
        if (present(info)) info = 'full sun'
        fraction = one
    end if

    end subroutine solar_fraction