cubic_shadow_model Subroutine

private subroutine cubic_shadow_model(rsun, radsun, rplanet, radplanet, sunfrac, rbubble)

The "circular cubic" shadow model.

Reference

  • J. Williams, et. al, "A new eclipse algorithm for use in spacecraft trajectory optimization", 2023, AAS 23-243

Arguments

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

apparent position vector of sun wrt spacecraft [km]
real(kind=wp), intent(in) :: radsun

radius of sun [km]
real(kind=wp), intent(in), dimension(3) :: rplanet

apparent position vector of eclipsing body wrt spacecraft [km]
real(kind=wp), intent(in) :: radplanet

radius of the eclipsing body [km]
real(kind=wp), intent(out) :: sunfrac

value of the function (>0 no eclipse, <0 eclipse, =0 on the shadow line)
real(kind=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.

Calls

proc~~cubic_shadow_model~~CallsGraph proc~cubic_shadow_model cubic_shadow_model acosd acosd proc~cubic_shadow_model->acosd asind asind proc~cubic_shadow_model->asind unit unit proc~cubic_shadow_model->unit wrap_angle wrap_angle proc~cubic_shadow_model->wrap_angle

Called by

proc~~cubic_shadow_model~~CalledByGraph proc~cubic_shadow_model cubic_shadow_model proc~get_sun_fraction get_sun_fraction proc~get_sun_fraction->proc~cubic_shadow_model proc~generate_eclipse_data mission_type%generate_eclipse_data proc~generate_eclipse_data->proc~get_sun_fraction proc~halo_solver_main halo_solver_main proc~halo_solver_main->proc~generate_eclipse_data

Source Code

    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<radplanet) then
        ! if inside the body, just return value from the surface
        r    = radplanet * unit(r)
        rmag = radplanet
    end if
    rsb        = rplanet - rsun
    alpha      = safe_acosd(dot_product(unit(r),unit(rsb)))
    if (compute_bubble) alpha = rad2deg*abs(wrap_angle(alpha*deg2rad-abs(rbubble)/rmag))
    rsbmag     = norm2(rsb)
    tmp        = (radsun + radplanet) / rsbmag
    sin_alpha0 = (one/rmag)*(radplanet*sqrt((one/tmp)**2-one)+sqrt(rmag**2-radplanet**2))*tmp
    alpha0     = safe_asind(sin_alpha0)
    sunfrac    = (alpha**2/(alpha0**2-alpha0**3/270.0_wp))*(one-alpha/270.0_wp)-one

    contains

        pure real(wp) function safe_asind(x)
            !! `asind` with range checking
            real(wp),intent(in) :: x
            safe_asind = asind(min(one,max(-one,x)))
        end function safe_asind

        pure real(wp) function safe_acosd(x)
            !! `acosd` with range checking
            real(wp),intent(in) :: x
            safe_acosd = acosd(min(one,max(-one,x)))
        end function safe_acosd

    end subroutine cubic_shadow_model