sub qnopt (n%, tol, maxiter%, x(), f, niter%, ier%, nsig%)

' quasi-newton multi-dimensional minimization subroutine

' input
' n% = number of variables
' tol = convergence criteria
' maxiter% = maximum number of function evaluations
' x() = current guess for solution
' (n% rows by 1 column)

' output
' x() = solution vector (n% rows by 1 column)
' f = current objective function value
' niter% = number of function evaluations
' ier% = error indicator
' nsig% = number of significant digits

' note: requires function log10 and subroutines
' utility and usr_func

sub usr_func (n%, x(), f)

' objective function subroutine

'''''''''''''''''''''''''''''''

local i%, jdate1, jdate2, tof

local ri(3), vi(3), rf(3), vf(3)

local sv1(6), sv2(6)

local dvm1, dvm2

jdate1 = x(1) + jdate0

jdate2 = x(2) + jdate0

' compute initial state vector

planet(ip1%, jdate1, oevi(), ri(), vi())

' compute final state vector

planet(ip2%, jdate2, oevf(), rf(), vf())

'svprint(rf(), vf())

tof = (jdate2 - jdate1) * 86400.0

' set up for Gooding's method

for i% = 1 to 3

sv1(i%) = ri(i%)

sv2(i%) = rf(i%)

sv1(i% + 3) = vi(i%)

sv2(i% + 3) = vf(i%)

next i%

' solve Lambert's TPBVP using Gooding's method

glambert(xmu, sv1(), sv2(), tof, vito(), vfto())

' calculate departure delta-v (kilometers/second)

dv1(1) = vito(1) - vi(1)
dv1(2) = vito(2) - vi(2)
dv1(3) = vito(3) - vi(3)

dvm1 = vecmag(dv1())

' calculate arrival delta-v (kilometers/second)

dv2(1) = vf(1) - vfto(1)
dv2(2) = vf(2) - vfto(2)
dv2(3) = vf(3) - vfto(3)

dvm2 = vecmag(dv2())

' load scalar objective function

select case otype%

case 1

' launch delta-v

f = dvm1

case 2

' arrival delta-v

f = dvm2

case 3

' launch + arrival delta-v

f = dvm1 + dvm2

end select

end sub

program ipto - interplanetary trajectory optimization
=====================================================

initial orbital elements and state vector
(heliocentric ecliptic and equinox of J2000)
--------------------------------------------

Earth

calendar date June 7 2003

UTC time 20 hours 36 minutes 49.43 seconds

semimajor axis 1.0000010180 AU
eccentricity 0.0167071774
orbital inclination 0.0004 degrees
argument of perigee 288.0835 degrees
RAAN 174.8649 degrees
true anomaly 153.6979 degrees

position vector x-component -3.506686626636e+07 kilometers
position vector y-component -1.477252709901e+08 kilometers
position vector z-component 1.174872460684e+03 kilometers

position vector magnitude 1.518303026416e+08 kilometers

velocity vector x-component 2.849840251725e+01 kilometers/second
velocity vector y-component -6.991575298582e+00 kilometers/second
velocity vector z-component 3.450069449798e-05 kilometers/second

velocity vector magnitude 2.934350134514e+01 kilometers/second


final orbital elements and state vector
(heliocentric ecliptic and equinox of J2000)
--------------------------------------------

Mars

calendar date December 28 2003

UTC time 12 hours 56 minutes 32.04 seconds

semimajor axis 1.5236793420 AU
eccentricity 0.0934042294
orbital inclination 1.8494 degrees
argument of perigee 286.5316 degrees
RAAN 49.5463 degrees
true anomaly 72.9019 degrees

position vector x-component 1.443324368932e+08 kilometers
position vector y-component 1.659191955418e+08 kilometers
position vector z-component -7.016935399128e+04 kilometers

position vector magnitude 2.199114292449e+08 kilometers

velocity vector x-component -1.735727441052e+01 kilometers/second
velocity vector y-component 1.796707126749e+01 kilometers/second
velocity vector z-component 8.028848351405e-01 kilometers/second

velocity vector magnitude 2.499470441817e+01 kilometers/second

initial delta-v vector and magnitude
(heliocentric ecliptic and equinox of J2000)
--------------------------------------------

delta-v vector x-component 2895.449599 meters/second
delta-v vector y-component -653.699580 meters/second
delta-v vector z-component -21.760090 meters/second

delta-v magnitude 2700.672241 meters/second

orbital energy 7.293631 km^2/sec^2)


final delta-v vector and magnitude
(heliocentric ecliptic and equinox of J2000)
--------------------------------------------

delta-v vector x-component -2019.677630 meters/second
delta-v vector y-component 1609.283132 meters/second
delta-v vector z-component 790.405355 meters/second

delta-v magnitude 2700.672241 meters/second

orbital energy 7.293631 km^2/sec^2)


total deltav magnitude 5401.344482 meters/second

total orbital energy 14.587261 km^2/sec^2)

total time of flight 203.680354 days


transfer orbital elements and state vector after first delta-v
(heliocentric ecliptic and equinox of J2000)
--------------------------------------------

semimajor axis 1.2598682075 AU
eccentricity 0.1945049751
orbital inclination 0.0385 degrees
argument of perigee 177.2951 degrees
RAAN 77.3056 degrees
true anomaly 2.0456 degrees


transfer orbital elements and state vector prior to second delta-v
(heliocentric ecliptic and equinox of J2000)
--------------------------------------------

semimajor axis 1.2598682075 AU
eccentricity 0.1945049751
orbital inclination 0.0385 degrees
argument of perigee 177.2951 degrees
RAAN 77.3056 degrees
true anomaly 154.3793 degrees