Orbit Propagation

Figure 5

proc twobdy(in: psi0, r0, v0, tau; out: psi, r, v)
    set initial psi based on input psi0
    determine acceptable bounds for psi

    do until "exit loop" condition:
      compute terms of "C series", all functions of psi
      compute terms of "S series", all functions of C series

      compute dtau
      if (dtau == 0 ) exit loop
      compute d(dtau)/d(psi )

      do a Newton step: psi = psi - dtau/[d(dtau)/d(psi)]
      reset acceptable bounds for psi
      if psi lies within acceptable bounds, go to loop beginning
        and continue iterating
      else
        try several alternative simple adjustments to psi
        if any of these yield a psi which lies within acceptable
          bounds, go to loop beginning and continue iterating
        else
          exit loop
        end if
      end if

  end do

  compute "f,g, fdot, and gdot expressions" as functions of
    the S terms, evaluated at the final psi from the loop above

  compute r and v as functions of f,g,fdot, and gdot

end proc