solar_fraction_alt Subroutine

public subroutine solar_fraction_alt(d_s, rs, d_p, rp, percentsun, info)

Another eclipse model, using circular area assumptions.

References

  • Montenbruck and Gill, "Satellite Orbits".
  • The GMAT routine ShadowState::FindShadowState.

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) :: 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_alt~~CallsGraph proc~solar_fraction_alt solar_fraction_alt proc~unit unit proc~solar_fraction_alt->proc~unit

Called by

proc~~solar_fraction_alt~~CalledByGraph proc~solar_fraction_alt solar_fraction_alt proc~get_sun_fraction get_sun_fraction proc~get_sun_fraction->proc~solar_fraction_alt proc~lighting_module_test lighting_module_test proc~lighting_module_test->proc~solar_fraction_alt

Source Code

    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)<theta) .and. (theta < rho_s + rho_p) ) then ! penumbra
                if (present(info)) info = 'penumbra'
                ! see montenbruck and gill, eq. 3.87-3.94
                c          = acos(dot_product(d_p_hat,d_s_hat))
                a2         = rho_s*rho_s
                b2         = rho_p*rho_p
                x          = (c*c + a2 - b2) / (two * c)
                y          = sqrt(a2 - x*x)
                area       = a2*acos(x/rho_s) + b2*acos((c-x)/rho_p) - c*y
                percentsun = one - area / (pi * a2)
            else ! antumbra
                if (present(info)) info = 'antumbra'
                percentsun =  one - rho_p*rho_p/(rho_s*rho_s)
            end if
        end if
    end if

    end subroutine solar_fraction_alt