c SUBROUTINE ClarkBert
c Integrated Lambert Transfers
c Copyright 2000, William H. Clark II
c All Rights Reserved c Last modified 6/3/00 c
CALL gui STOP END
C C
SUBROUTINE TailHookm (x,t,deltaV,Vrel) implicit real*8 (a-h,o-z)
common/output/tj common/const0/pi,tol,iflag,itarget,idim
common/constm2/xm, gmmars, radM,rm,soiM,omm,phM,rpM real*8 x(6),xm(6),x12(6),xg(6),tj(30,3)
dt = +100d0
do 75 n = 1,1500
CALL RK78T (t,x,dt,tol,6)
range = range2 (x,xm,x12,sense,Vrel) rhoM = dsqrt( xm(1)**2+xm(2)**2+xm(3)**2) rhoX = dsqrt( x(1)**2+ x(2)**2+ x(3)**2) cable1 = rhoM – soiM*50d0
cable2 = rhoM – soiM*3d00
IF (rhoX .GT. cable1) dt = +86400d0/1d0 IF (rhoX .GT. cable2) dt = +86400d0/4d0 IF (rhoX .GT. rhoM -soiM) THEN
deltaV=dmod(datan2(x(2),x(1))-datan2(xm(2),xm(1)),2*pi) GOTO 78
ENDIF
range1 = helioTOgeo (x,t,4,xg) axis = datan2(xg(2),xg(1)) zone = datan2(xg(5),xg(4))
IF (t/86400d0 .GT. 250d0) THEN GOTO 78
ENDIF
75 continue
78 return end
c
SUBROUTINE catapult (beta1,beta2,deltaVi,deltaVf,tsip) implicit real*8 (a-h,o-z)
common/output/tj common/constopt/alpha,beta,t1,t2 common/const0/pi,tol,iflag,itarget,idim dimension tj(30,3)
c
pi = dacos(-1.d0) tol = 1.d-13 CALL INRK78
idim = 2
c
CALL constants
CALL HawkeyM (beta2,deltaVi,dVmars) CALL constants
CALL HawkeyE (beta1,dVearth,deltaf) tsip = (t2-t1)/86400d0
deltaVf = ABS( dVearth) + ABS(deltaf) + ABS( deltaVi) + ABS(dVmars) tj(22,1) = tsip
tj(22,2) = deltaVf
RETURN END
C
SUBROUTINE HawkeyM (beta2,deltaVi,dVmin) implicit real*8 (a-h, o-z)
dimension x(6),xg(6),tj(30,3) common/output/tj common/constopt/alpha,beta,t1,t2 common/const0/pi,tol,iflag,itarget,idim
c
iflag = 1 beta = beta2 alpha = 0.d0
c
do 2 i = 1,10
CALL initial(deltaVi,x,t)
CALL TailHookm (x,t,delta,Vrel) alpha = alpha + delta
IF (ABS(delta) .LT. 1d-6) GOTO 3
2 continue
c
3 tmars = t tj(13,1) = x(1)
tj(13,2) = x(2)
tj(14,1) = x(4)
tj(14,2) = x(5)
tj(15,1) = tmars/86400.d0 range2 = helioTOgeo (x,t,4,xg) tj(16,1) = xg(1)
tj(16,2) = xg(2)
tj(17,1) = xg(4)
tj(17,2) = xg(5)
CALL phoenix (t,Vrel,angle,dVmin) t2 = tmars + t
RETURN END
c
SUBROUTINE phoenix (t,Vrel,angle,dVmin)
c level THREE diagnostics: seek free return trajectory implicit real*8 (a-h,o-z)
dimension xe(6),xm(6),xg(6),xtable(40,14),tj(30,3) common/output/tj common/const0/pi,tol,iflag,itarget,idim
common/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE common/constm2/xm, gmmars, radM,rm,soiM,omm,phM,rpM
c
Rf = radM + rpM
gamma = -pi/8.d0 + .0003d0 do 4 j = 1,40
gamma = gamma + pi/220.d0 xtable (j,1) = -soiM*dsin(gamma) xtable (j,2) = -soiM*dcos(gamma) xtable (j,3) = 0.d0
xtable (j,4) = 0.d0 xtable (j,5) = Vrel xtable (j,6) = 0.d0 xtable (j,7) = 0.d0 xtable (j,8) = 0.d0 xtable (j,9) = 0.d0 xtable (j,10)= 0.d0 xtable (j,11)= 0.d0 xtable (j,12)= 0.d0 xtable (j,13)= 0.d0 xtable (j,14)= 0.d0
4 continue
gmsun = gmmars gmearth = 0.d0 gmmars = 0.d0 gmjupitr = 0.d0
tmax = 10.d0*soiM/Vrel
c
do 78 m = 1,40
do 8 k = 1,6
8 xg(k) = xtable (m,k) t = 0.d0
dt = +1000.d0
dVmin = 100.d0 Vrel = xtable(m,5) do 75 n=1,1500
CALL RK78T(t,xg,dt,tol,6)
range = sqrt(xg(1)**2 + xg(2)**2 + xg(3)**2)
IF (range .LT. radM) THEN deltaV = deltaV*2.d0 GOTO 78
ENDIF
CALL velocity (xg,gmsun,range,Vrel,Rf,deltaV,V1,V2,V3)
IF (deltaV .LT. dVmin) THEN do 10 k = 1,6
10 xtable (m,k) = xg(k) xtable (m,10) = deltaV xtable (m,9) = gamma xtable (m,8) = range dVmin = deltaV xtable(m,11) = V1 xtable(m,12) = V2 xtable(m,13) = V3 xtable(m,14) = t/86400d0
ENDIF
IF (range .GT. soiM*1.1d0) GOTO 78 IF (t .GT. tmax) GOTO 78
75 continue
78 continue
do 80 m = 1,40
IF (xtable (m,10) .LT. dVmin .AND. xtable(m,8) .GT. radM) THEN tj(18,1) = xtable(m,1)
tj(18,2) = xtable(m,2) tj(19,1) = xtable(m,4) tj(19,2) = xtable(m,5) tj(20,1) = xtable(m,14) tj(21,1) = xtable(m,12) tj(21,2) = xtable(m,11) tj(21,3) = xtable(m,13) dVmin = xtable(m,10) angle = xtable(m,9)
ENDIF
80 continue CALL constants RETURN
END
C
SUBROUTINE trapE (ax,step,dX,ttol,adelta) implicit real*8 (a-h,o-z)
C beta angle is fixed >> find min. deltaV to Rpark CALL SG (ax,adelta)
c print *, ax,adelta 2 DO 3 k = 1,50
IF (step .LT. ttol) GOTO 5
bx = ax + step*dX CALL SG (bx,bdelta)
c print *, bx,bdelta
IF (bdelta .GT. adelta) GOTO 4 ax = bx
adelta = bdelta
-
- bx = ax – step*dX step = step/2.d0 GOTO 2
C
SUBROUTINE SG (dV,delta)
implicit real*8 (a-h,o-z) real*8 x(6)
common/const0/pi,tol,iflag,itarget,idim
common/constopt/alpha,beta,t1,t2 SELECT CASE (itarget)
CASE (1)
CALL initial (dV,x,t)
CALL TailHooke (X,t,0.d0,delta) t1 = t
CASE (2)
CALL initial (dV,x,t)
CALL TailHookm (X,t,0.d0,delta) t2= t
END SELECT
RETURN END
C
SUBROUTINE HawkeyE (beta1,deltaVf,deltaf)
c level ONE diagnostics: target quadrant of SOI implicit real*8 (a-h,o-z)
dimension x(6),xe(6),tj(30,3)
common/output/tj common/constopt/alpha,beta,t1,t2 common/const0/pi,tol,iflag,itarget,idim
common/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE itarget = 1 ! Earth
iflag = 2 beta = beta1
deltaV = 0.d0
t = 0d0
c >>> first be sure s/c crosses outside of Earth cable2 = radE + rpE
CALL initial (deltaV,x,t)
CALL TailHooke (X,t,0d0,delta)
IF (x(1)-xe(1) .LT. cable2) THEN dX = +1.d0
ELSE
dX = -1.d0 ENDIF
IF (x(1) – xe(1) .GT. cable2*1.5d0) GOTO 5
do 4 i = 1,10
deltaV = deltaV + 0.005d0*dX
CALL initial (deltaV,x,t)
CALL TailHooke (X,t,0d0,delta)
IF (x(1) – xe(1) .GT. cable2*1.5d0) GOTO 5
4 continue
5 iflag = 1 delta0 = deltaV dX = -1.0*dX
step = 0.0001d0 ttol = 1.0d-6
CALL trapE (delta0,step,dX,ttol,deltaVf) deltaVf = delta0
CALL initial (deltaVf,x,t)
CALL TailHooke (X,t,0d0,deltaf) RETURN
END
C
SUBROUTINE earth (x,t,deltaVf) implicit real*8 (a-h,o-z)
real*8 x(6),x2(6),xg(6),x12(6),xe(6),xm(6),xj(6),tj(30,3) common/output/tj common/const0/pi,tol,iflag,itarget,idim common/constopt/alpha,beta,t1,t2
common/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE common/constm2/xm, gmmars, radM,rm,soiM,omm,phM,rpM common/constj2/xj, gmjuptr,rj,omj,phJ
c
Rf = radE + rpE dt = -1000d0
range1 = helioTOgeo (x,t,3,xg) tj(7,1) = xg(1)
tj(7,2) = xg(2)
tj(8,1) = xg(4)
tj(8,2) = xg(5)
tj(9,1) = x(1)
tj(9,2) = x(2)
tj(10,1) = x(4)
tj(10,2) = x(5) tj(6,2) = t/86400d0
gm = gmsun gmsun = gmearth gmearth = gm
SELECT CASE (idim) CASE (1)
phE = pi – (t*ome + phE) CASE (2)
do 2, i = 1,6
2 xe(i) = -1d0*xe(i) END SELECT
gmmars = 0d0
gmjupitr = 0d0 deltaVf = 1000d0
c
do 75 n=1,250
CALL RK78T(t,xg,dt,tol,6)
range = range2 (xg,x2,x12,sense,Vrel)
IF (range .LT. radE) THEN deltaVf = deltaVf*2d0 GOTO 78
ENDIF
CALL velocity (xg,gmsun,range,Vrel,Rf,deltaV,V1,V2,V3)
IF (deltaV .LT. deltaVf) THEN deltaVf = deltaV
rangef = range xgf = xg(1) tj(3,1) = V1 tj(3,2) = V2 tj(3,3) = V3
tj(4,1) = xg(1)
tj(4,2) = xg(2)
tj(5,1) = xg(4)
tj(5,2) = xg(5) tj(6,1) = t/86400d0
ENDIF
IF (range .GT. soiE) GOTO 78 75 continue
78 t1 = t
CALL constants RETURN
END
SUBROUTINE TailHooke (x,t,final,deltaV) implicit real*8 (a-h,o-z)
common/output/tj common/const0/pi,tol,iflag,itarget,idim
common/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE real*8 x(6),xe(6),x12(6),tj(30,3)
dt = -100d0
delta = soiE*2d0 do 75 n = 1,1500
CALL RK78T (t,x,dt,tol,6)
rhoE = dsqrt(xe(1)**2+xe(2)**2+xe(3)**2) rhoX = dsqrt( x(1)**2+ x(2)**2+ x(3)**2) cable1 = rhoE + delta
cable2 = rhoE + soiE
IF (rhoX .LT. cable1) dt = -86400d0/3d0 IF (rhoX .LT. cable2) dt = -86400d0/6d0 SELECT CASE (iflag)
CASE (1)
range = range2(x,xe,x12,sense,Vrel)
IF (sense .GT. 0.d0 .OR. x12(1) .GT. soiE) THEN
Vrel = Vrel + range/soiE
CALL velocity (x,gmearth,soiE,Vrel,radE+rpE,deltaV,V1,V2,V3) deltaV = deltaV + 0.5d0
GOTO 78 ENDIF
IF (range .LT. soiE) THEN
IF (final .EQ. 1.d0) THEN
c CALL Phoenix
GOTO 78 ENDIF
CALL earth (x,t,deltaV)
GOTO 78 ENDIF CASE (2)
alphaX = datan2(x(2),x(1)) alphaE = datan2(xe(2),xe(1)) gamma = ABS(alphaX – alphaE)
IF (gamma .LT. 0.15d0) dt = -4*60*60
IF (gamma .LT. 0.02d0 .AND. t/86400d0 .LT. -40d0) GOTO 78 END SELECT
75 continue
78 return end
c======== sanitized routines =================
SUBROUTINE velocity (xg,gmsun,range,Vrel,Rf,deltaV,V1,V2,V3) implicit real*8 (a-h,o-z)
dimension xg(6)
Rp = range
phi1 = datan(xg(4)/xg(5)) phi2 = datan(xg(2)/xg(1)) gamma = ABS(phi1 – phi2) Vcirc = dsqrt(gmsun/range)
V3 =dsqrt( Vrel**2+Vcirc**2-2.d0*Vrel*Vcirc*dcos(gamma)) V2 = dsqrt(gmsun/Rp)*(1.d0-dsqrt(2.d0/(1.d0+Rp/Rf)))
V1 = dsqrt(gmsun/Rf)*( dsqrt((2.d0*Rp/Rf)/(1+Rp/Rf))-1.d0)
deltaV = ABS(V1) + ABS(V2) + ABS(V3) RETURN
END
c
SUBROUTINE thirdbc(t,n,x3) implicit real*8 (a-h,o-z)
common/const0/pi,tol,iflag,itarget,idim
common/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE common/constm2/xm, gmmars, radM,rm,soiM,omm,phM,rpM common/ephemers/ej
real*8 x3(6),xe(6),xm(6),ej(4,6) SELECT CASE (idim)
CASE (1)
SELECT CASE (n) CASE (3)
r3 = re
theta = t*ome + phE CASE (4)
r3 = rm
theta = t*omm + phM END SELECT
x3(1) = r3*dcos(theta)
x3(2) = r3*dsin(theta) x3(3) = 0.0d0
CASE (2) SELECT CASE (n) CASE (3)
pe = ej(3,1)*(1d0 – ej(3,2)**2) xeM = phE + ome*t
CALL randv (gmsun,pe,ej(3,2),ej(3,3),ej(3,4),ej(3,5),xeM,xe) CASE (4)
pm = ej(4,1)*(1d0 – ej(3,2)**2) xmM = phM + omm*t
CALL randv (gmsun,pm,ej(4,2),ej(4,3),ej(4,4),ej(4,5),xmM,xm) END SELECT
END SELECT RETURN END
FUNCTION helioTOgeo (xh,t,n,xg) implicit real*8 (a-h,o-z)
real*8 xh(6),x3(6),xg(6),xe(6),xm(6)
common/const0/pi,tol,iflag,itarget,idim common/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE common/constm2/xm, gmmars, radM,rm,soiM,omm,phM,rpM SELECT CASE (idim)
CASE (1)
SELECT CASE (n) CASE (1)
r3 = rm
theta = pi/2.d0 + datan(ABS(xh(1))/xh(2)) x3(1) = r3*dcos(theta)
x3(2) = r3*dsin(theta) x3(3) = 0.0d0
CASE (3)
r3 = re
theta = t*ome + phE CASE (4)
r3 = rm
theta = t*omm + phM END SELECT
x3(1) = r3*dcos(theta) x3(2) = r3*dsin(theta) x3(3) = 0.0d0
iquad = INT ((theta + pi/2.d0)*2.d0/pi) SELECT CASE (iquad)
CASE (1)
signVx = -1.d0 signVy = +1.d0 CASE (2)
signVx = -1.d0 signVy = -1.d0 CASE (3)
signVx = +1.d0 signVy = -1.d0 CASE (4)
signVx = +1.d0 signVy = +1.d0
END SELECT
Vcirc = dsqrt(gmsun/r3)
x3(4) = signVx*ABS(dsin(theta))*Vcirc x3(5) = signVy*ABS(dcos(theta))*Vcirc x3(6) = 0.d0
range = range2(xh,x3,xg,sense,Vrel) CASE (2)
SELECT CASE (n) CASE (3)
range = range2(xh,xe,xg,sense,Vrel) CASE (4)
range = range2(xh,xm,xg,sense,Vrel) END SELECT
helioTOgeo = range END SELECT RETURN
END
c
FUNCTION range2 (x1,x2,x12,sense,Vrel) implicit real*8 (a-h,o-z) common/const0/pi,tol,iflag,itarget,idim dimension x1(6),x2(6),x12(6)
c relative position & velocity of 1 w.r.t. 2 sense = 1.d0
range = 0.d0 Vrel = 0.d0 rdotV = 0.d0 do 4 i = 1,3
j = i+3
x12(i) = x1(i) – x2(i)
x12(j) = x1(j) – x2(j) range = range + x12(i)**2 Vrel = Vrel + x12(j)**2
rdotV = rdotV + x12(i)*x12(j) 4 continue
range2 = dsqrt(range) Vrel = dsqrt(Vrel)
beta = dacos(rdotV/(range2*Vrel))
IF (beta .LT. pi/2) sense = -1.d0
c range is openinging for beta < 90 deg RETURN
END
SUBROUTINE constants implicit real*8 (a-h,o-z) real*8 xe(6),xm(6),xj(6)
common/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE common/constm2/xm, gmmars, radM,rm,soiM,omm,phM,rpM common/constj2/xj, gmjuptr,rj,omj,phJ
gmsun = 1.32712428d11 gmearth = 3.986004415d5
re = 149598023d0 ! = a
ome = dsqrt((gmsun+gmearth)/re)/re soiE = 924647d0
radE = 6378.1363d0 rpE = 200.d0 gmmars = 4.305d4
rm = 227939186d0 ! = a
omm = dsqrt(( gmsun+gmmars)/rm)/rm soiM = 577213d0
radM = 3397.2d0
rpM = 80.d0 gmjuptr = 1.268d8
rj = 778298361d0 ! = a
omj = dsqrt((gmsun+gmjuptr)/rj)/rj RETURN
END
C
SUBROUTINE Hohmann (x,phase) implicit real*8 (a-h,o-z)
dimension x(6),xe(6),xm(6),ej(4,6),tj(30,3) common/output/tj
common/ephemers/ej common/const0/pi,tol,iflag,itarget,idim
common/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE common/constm2/xm, gmmars, radM,rm,soiM,omm,phM,rpM CALL ephemerides
r1 = ej(3,1) r2 = ej(4,1)
a = (r1 + r2)/2.d0
e = (r2 – r1)/(r1 + r2) p = a*(1.d0 – e**2)
period = pi*(sqrt(a)**3)/dsqrt(gmsun) tconj = (pi – period*omm)/(ome-omm)
c
phase = tconj*ome
IF (idim .EQ. 2) THEN
CALL conjunction (t,f0,fe,fm,xMe,xMm) phE = xMe
phM = xMm
f1 = dmod(fe – phase,2d0*pi)
f2 = dmod(fm + pi – phase,2d0*pi) pe = ej(3,1)*(1d0 – ej(3,2)**2)
pm = ej(4,1)*(1d0 – ej(4,2)**2) r1 = pe/(1d0 + ej(3,2)*dcos(f1)) r2 = pm/(1d0 + ej(4,2)*dcos(f2)) a = (r1 + r2)/2.d0
e = (r2 – r1)/(r1 + r2) p = a*(1.d0 – e**2)
omsc = dsqrt((gmsun)/a)/a
E0e=2d0*datan(dsqrt((1d0-ej(3,2))/(1d0+ej(3,2)))*tan(f1/2d0)) e0M= E0e – ej(3,2)*dsin(E0e)
tconj = (xMe – e0M)/ome
scM = omsc*tconj
CALL Ekepler (e,scM,scE) scf=2d0*datan(dsqrt((1d0+ej(3,2))/(1d0-ej(3,2)))*tan(scE/2d0)) omega = f0 – phase
CALL randv (gmsun,p,e,0.d0,omega,0.d0,scM,x)
GOTO 3 ENDIF
omsc = dsqrt((gmsun)/a)/a
scM = omsc*tconj omega = 0.d0
-
- CALL randv (gmsun,p,e,0.d0,omega,0.d0,scM,x)
C
SUBROUTINE initial (deltaV,X,t) implicit real*8 (a-h,o-z)
real*8 X(6),xe(6),xm(6),xj(6),tj(30,3) common/output/tj common/constopt/alpha,beta,t1,t2 common/const0/pi,tol,iflag,itarget,idim
common/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE common/constm2/xm, gmmars, radM,rm,soiM,omm,phM,rpM common/constj2/xj, gmjuptr,rj,omj,phJ
c
CALL Hohmann(x,phase) IF (idim .EQ. 1) THEN
c phase = 1.65265383684363d0
phE = phase phM = phase
phJ = phase ! pi – dasin(x(2)/rj) ENDIF
IF (beta .NE. 0.d0 .OR. deltaV .NE. 0.d0) THEN
gamma = ABS(datan(X(5)/X(4)))
X(4) = X(4) – deltaV*dcos(gamma + beta) IF (X(5) .GE. 0.d0) THEN
X(5) = X(5) + deltaV*dsin(gamma + beta) ELSE
X(5) = X(5) – deltaV*dsin(gamma – beta)
ENDIF ENDIF tj(11,1) = x(1)
tj(11,2) = x(2)
tj(12,1) = x(4)
tj(12,2) = x(5)
t = 0.d0 RETURN END
SUBROUTINE DERIV(t,x,f)
implicit real*8 (a-h,o-z)
real*8 x(7),f(6),xe(6),xm(6),xj(6) common/const0/pi,tol,iflag,itarget,idim common/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE common/constm2/xm, gmmars, radM,rm,soiM,omm,phM,rpM common/constj2/xj, gmjuptr,rj,omj,phJ
f(1) = x(4)
f(2) = x(5)
f(3) = x(6)
r2 = x(1)**2 + x(2)**2 + x(3)**2
r = dsqrt(r2)
f(4) = -gmsun*(x(1)/r)/r2
f(5) = -gmsun*(x(2)/r)/r2
f(6) = -gmsun*(x(3)/r)/r2 CALL thirdbc(t,3,xe)
CALL pert3b(gmearth,xe,x,f) CALL thirdbc(t,4,xm)
CALL pert3b(gmmars,xm,x,f) IF (idim .EQ. 1) THEN
CALL thirdbc(t,5,xj)
CALL pert3b(gmjuptr,xj,x,f) ENDIF
RETURN END
SUBROUTINE ephemerides ! J2000 (Vallado) implicit real*8 (a-h,o-z)
common/ephemers/ej
common/const0/pi,tol,iflag,itarget,idim real*8 ej (4,6)
factor = pi/180.d0
c
ej(3,1) = 149598023.d0 ! a
ej(3,2) = .016708617d0 ! e ej(3,3) = 0.d0 ! i
ej(3,4) = 0.d0 ! LONGITUDE of ascending node
ej(3,5) = 102.93734808d0*factor ! longongitude of perihelion ej(3,6) = 100.46644851d0*factor ! true longitude @ epoch
c
ej(4,1) = 227939186.d0 ej(4,2) = .093400620d0 IF (idim .EQ. 3) THEN
ej(4,3) = 1.84972648d0*factor ELSE
ej(4,3) = 0.d0 ENDIF
ej(4,4) = 49.55809321d0*factor
ej(4,5) = 336.06023398d0*factor ej(4,6) = 355.43327463d0*factor RETURN
END
SUBROUTINE energye(x,t,E) implicit real*8 (a-h,o-z) real*8 x(9),xb(6),xe(6)
common/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE r = dsqrt(x(1)**2 + x(2)**2 + x(3)**2)
E2body = 0.5d0*(x(4)**2 + x(5)**2 + x(6)**2) – gmsun/r
C
CALL thirdbc(t,re,ome,phE,xb)
r23 = dsqrt((xb(1)-x(1))**2 + (xb(2)-x(2))**2 + (xb(3)-x(3))**2) Um1 = gmearth/r23
Um2 = ome*(x(1)*x(5)-x(2)*x(4))
Um3 = (gmearth/re**3)*(X(1)*xb(1) + X(2)*xb(2))
C
E = E2body -Um1 – Um2 + Um3 RETURN
END
C
SUBROUTINE pert3b(GMP,x3,x,f)
real*8 f(6),x(3),x3(3),gmp,r12,r23
r12 = dsqrt(x3(1)**2 + x3(2)**2 + x3(3)**2)
r23 = dsqrt((x3(1)-x(1))**2+(x3(2)-x(2))**2+(x3(3)-x(3))**2) f(4) = f(4) – gmp*((x(1)-x3(1))/r23**3 + x3(1)/r12**3)
f(5) = f(5) – gmp*((x(2)-x3(2))/r23**3 + x3(2)/r12**3)
f(6) = f(6) – gmp*((x(3)-x3(3))/r23**3 + x3(3)/r12**3) RETURN
END
SUBROUTINE RK78T(T,X,DT,TOL,N)
IMPLICIT double precision (A-H,O-Z) common/RKCOM/CH(13),AL(13),B(13,12)
DIMENSION XD(21),F(21,13),X(N),F1(21),F2(21),F3(21),F4(21),F5(21)
c CALL INRK78 TO INITIALIZE
c IF(DABS(DT).LT.1.D-20)RETURN TM=T
DT1=DT
C
DO 40 I=1,N
GO TO 1020
C
1010 T=TM DT1=DT
DO 41 I=1,N
C
1020 CONTINUE
CALL DERIV(T,X,F1) C K=2
DO 602 I=1,N
TP=B(2,1)*F1(I)
602 X(I)=XD(I)+DT*TP T=TM+AL(2)*DT CALL DERIV(T,X,F2)
C K=3
DO 603 I=1,N TP=B(3,1)*F1(I)+B(3,2)*F2(I)
603 X(I)=XD(I)+DT*TP
T=TM+AL(3)*DT CALL DERIV(T,X,F3)
C K=4
DO 604 I=1,N TP=B(4,1)*F1(I)+B(4,3)*F3(I)
604 X(I)=XD(I)+DT*TP T=TM+AL(4)*DT CALL DERIV(T,X,F4)
C K=5
DO 605 I=1,N TP=B(5,1)*F1(I)+B(5,3)*F3(I)+B(5,4)*F4(I)
605 X(I)=XD(I)+DT*TP
T=TM+AL(5)*DT CALL DERIV(T,X,F5)
C
DO 50 K=6,13 KK=K-1
DO 71 I=1,N TP=B(K,1)*F1(I)+B(K,4)*F4(I)+B(K,5)*F5(I) IF(KK.LT.6) GOTO 71
DO 70 J=6,KK
70 TP=TP+B(K,J)*F(I,J)
71 X(I)=XD(I)+DT*TP T=TM+AL(K)*DT
50 CALL DERIV(T,X,F(1,K))
C
DO 101 I=1,N TP=0.D0
C TP=CH(1)*F1(I)+CH(2)*F2(I)+CH(3)*F3(I)+CH(4)*F4(I)+CH(5)*F5(I) DO 100 L=6,13
100 TP=TP+CH(L)*F(I,L)
101 X(I)=XD(I)+DT*TP C
IF(TOL.EQ.0.D0) GOTO 900 ER=0.D0
DO 112 I=1,N A=DABS(X(I)) IF(A.LT.1.D-6) A=1.D-6
TEI=DABS(F1(I)+F(I,11)-F(I,12)-F(I,13))*CH(12) / A IF(TEI.GT.ER) ER=TEI
112 CONTINUE ER=ER*DABS(DT)+1.D-16 DT=DT*(TOL/ER)**.125D0 IF(ER.GT.TOL) GOTO 1010
C
900 CONTINUE T=TM+DT1 RETURN END
SUBROUTINE INRK78
implicit double precision (a-h,o-z) common/RKCOM/CH(13),AL(13),B(13,12) DO 1 I=1,13
CH(I)=0.D0 AL(I)=0.D0 DO 1 J=1,12
1 B(I,J)=0.D0 CH(6)=34.D0/105.D0 CH(7)=9.D0/35.D0 CH(8)=CH(7) CH(9)=9.D0/280.D0 CH(10)=CH(9)
CH(12)=41.D0/840.D0 CH(13)=CH(12) AL(2)=2.D0/27.D0 AL(3)=1.D0/9.D0 AL(4)=5.D0/30.D0 AL(5)=5.D0/12.D0 AL(6)=1.D0/2.D0 AL(7)=5.D0/6.D0 AL(8)=5.D0/30.D0 AL(9)=2.D0/3.D0 AL(10)=1.D0/3.D0 AL(11)=1.D0 AL(13)=1.D0 B(2,1)=2.D0/27.D0 B(3,1)=1.D0/36.D0 B(4,1)=5.D0/120.D0 B(5,1)=5.D0/12.D0 B(6,1)=1.D0/20.D0 B(7,1)=-25.D0/108.D0 B(8,1)=31.D0/300.D0 B(9,1)=2.D0
B(10,1)=-91.D0/108.D0 B(11,1)=2383.D0/4100.D0 B(12,1)=3.D0/205.D0 B(13,1)=-1777.D0/4100.D0 B(3,2)=1.D0/12.D0 B(4,3)=1.D0/8.D0
B(5,3)=-25.D0/16.D0 B(5,4)=25.D0/16.D0 B(6,4)=1.D0/4.D0 B(7,4)=125.D0/108.D0 B(9,4)=-53.D0/6.D0 B(10,4)=23.D0/108.D0 B(11,4)=-341.D0/164.D0 B(13,4)=-341.D0/164.D0 B(6,5)=1.D0/5.D0
B(7,5)=-65.D0/27.D0 B(8,5)=61.D0/225.D0 B(9,5)=704.D0/45.D0 B(10,5)=-976.D0/135.D0 B(11,5)=4496.D0/1025.D0 B(13,5)=4496.D0/1025.D0 B(7,6)=125.D0/54.D0 B(8,6)=-2.D0/9.D0
B(9,6)=-107.D0/9.D0 B(10,6)=311.D0/54.D0 B(11,6)=-301.D0/82.D0 B(12,6)=-6.D0/41.D0 B(13,6)=-289.D0/82.D0 B(8,7)=13.D0/900.D0 B(9,7)=67.D0/90.D0 B(10,7)=-19.D0/60.D0 B(11,7)=2133.D0/4100.D0 B(12,7)=-3.D0/205.D0
B(13,7)=2193.D0/4100.D0 B(9,8)=3.D0 B(10,8)=17.D0/6.D0 B(11,8)=45.D0/82.D0 B(12,8)=-3.D0/41.D0 B(13,8)=51.D0/82.D0 B(10,9)=-5.D0/60.D0 B(11,9)=45.D0/164.D0 B(12,9)=3.D0/41.D0 B(13,9)=33.D0/164.D0 B(11,10)=18.D0/41.D0 B(12,10)=6.D0/41.D0 B(13,10)=12.D0/41.D0 B(13,12)=1.D0
RETURN END
*———————————————————–
SUBROUTINE vadd (v,t,n,r)
* scalar addition of arrays, v(n)+t(n)=r(n) real*8 v(n),t(n),r(n)
DO I=1,n r(I)=v(I)+t(I)
END DO RETURN END
SUBROUTINE cross (v,t,r)
* cross multiply v(3) x t(3) = r(3) real*8 v(3),t(3),r(3)
r(1) = v(2)*t(3) – v(3)*t(2)
r(2) = v(3)*t(1) – v(1)*t(3)
r(3) = v(1)*t(2) – v(2)*t(1) RETURN
END
SUBROUTINE dot (v,w,z)
* multiply v(3) dot w(3) = z real*8 z,v(3), w(3)
z = v(1)*w(1) + v(2)*w(2) + v(3)*w(3)
RETURN END
SUBROUTINE ds (s,v,n,r)
* scalar division of an array with n elements real*8 s,v(n),r(n)
DO I=1,n
r(I)=v(I)/s END DO RETURN END
SUBROUTINE exp (v,p,n,r)
* raising a vector v(n) to the power p real*8 p,v(n), r(n)
DO I=1,n
r(I) = v(I)**p END DO RETURN
END
SUBROUTINE mult (v,t,n,r)
* vector multiplication of array v(n) with array t(n) real*8 v(n),t(n),sum,r
sum=0.0D0 DO I=1,n
sum = sum + v(I)*t(I) END DO
r = sum RETURN END
SUBROUTINE mxm (e,f,n,r)
* matrix multiplication e(3,3) * f(3,3) = r(3,3) implicit real*8 (a-h,o-z)
real*8 e(3,3), f(3,3), r(3,3)
DO I=1,n
sum1=0.0d0 sum2=0.0d0 sum3=0.0d0 DO J=1,n
sum1 = sum1 + e(I,J)*f(J,1) sum2 = sum2 + e(I,J)*f(J,2) sum3 = sum3 + e(I,J)*f(J,3)
END DO
r(I,1) = sum1 r(I,2) = sum2 r(I,3) = sum3
END DO RETURN END
SUBROUTINE norm (v,n,s)
* calculation of the length of an n-dimensional array v(n) real*8 s,v(n),sum
sum = 0.0d0 DO I=1,n
sum = sum + v(I)**2
END DO
s = dsqrt(sum) RETURN END
SUBROUTINE vr3 (v,m,r)
* multiplication of 3-element row vector v(3) by matrix m(3,3) real*8 v(3),m(3,3),r(3),sum
sum = 0.0d0
DO I=1,3
DO J=1,3
sum = sum + v(I)*m(J,I) END DO
r(I) = sum sum = 0.0D0
END DO
RETURN END
SUBROUTINE xs (s,v,n,r)
* scalar multiplication of an array with n elements real*8 s,v(n),r(n)
DO I=1,n
r(I)=v(I)*s END DO RETURN END
*———celestial mechanics subroutines——————-SUBROUTINE eccen1 (gmu,r,v,e)
* calculate eccenticity vector from r,v,mu
implicit real*8 (a-h,l-z)
real*8 lrl,lvl,r(3),v(3),rte(3),rvv(3),rvv1(3),ex(3),e(3) CALL norm (r,3,lrl)
CALL norm (v,3,lvl) CALL dot (r,v,rv) CALL xs (rv,v,3,rvv)
CALL xs (-1.0D0,rvv,3,rvv1) te = lvl**2 – gmu/lrl
CALL xs (te,r,3,rte) CALL vadd (rte,rvv1,3,ex) CALL ds (gmu,ex,3,e) RETURN
END
SUBROUTINE elorb (x,gmu,p,a,lel,i,omega,w,f)
* calculate orbital elements from position and velocity vectors implicit real*8 (a-i,k-z)
real*8 r(3),v(3),h(3),k(3),n(3),e(3),x(6)
pi = dacos(-1.d0) do 2 j=1,3
k(j) = 0.d0
r(j) = x(j)
2 v(j) = x(j+3)
CALL norm (r,3,lrl) ! magnitude of r CALL norm (v,3,lvl) ! magnitude of v CALL dot (r,v,z)
CALL cross (r,v,h) ! angular momentum vector CALL norm (h,3,lhl) ! magnitude of h
k(3) = 1.0D0 ! create [0 0 1] vector
CALL cross (k,h,n) ! eccentricity vector CALL norm (n,3,lnl) ! magnitude of n
CALL eccen1 (gmu,r,v,e) CALL norm (e,3,lel)
energy = 0.5D0*lvl**2 – gmu/lrl
a=-0.5D0*gmu/energy p=lhl**2/gmu i=DACOS(h(3)/dabs(lhl))
IF (lnl .EQ. 0.d0) THEN
omega = 0.d0 ELSE
omega=DACOS(n(1)/dabs(lnl))
IF (n(2) .LT. 0.0D0) omega = 2.d0*pi – omega ENDIF
CALL dot (n,e,ne)
IF (i .eq. 0.d0) THEN
w = dacos(e(1)/lel) ELSE
w=DACOS(ne/(dabs(lnl*lel))) ENDIF
IF (e(3) .LT. 0.0D0) w = 2.d0*pi – w
CALL dot (e,r,er) f=DACOS(er/(dabs(lel*lrl)))
IF (z .LT. 0.0D0) f = 2.d0*pi – f
RETURN END
SUBROUTINE randv (gmu,p,e,i,omega,w,xM,x)
* calculate position and velocity vectors from orbital elements implicit real*8 (a-i,o-z)
real*8 rx1(3),rx(3),vx(3),vx1(3),ROT(3,3),x(6)
c
CALL Ekepler (e,xM,xE)
f = 2.d0*datan(dsqrt((1d0+e)/(1d0-e))*tan(xE/2.d0))
c
rx1(1)=p*DCOS(f)/(1.0D0+e*DCOS(f)) rx1(2)=p*DSIN(f)/(1.0D0+e*DCOS(f)) rx1(3)=0.0D0
vx1(1)=-1.0D0*DSIN(f)*DSQRT(gmu/p)
vx1(2)=(e+DCOS(f))*DSQRT(gmu/p) vx1(3)=0.0D0
ROT(1,1)=DCOS(omega)*DCOS(w)-DSIN(omega)*DSIN(w)*DCOS(i) ROT(1,2)=-1.0D0*DCOS(omega)*DSIN(w)-DSIN(omega)*DCOS(w)*DCOS(i) ROT(1,3)=DSIN(omega)*DSIN(i) ROT(2,1)=DSIN(omega)*DCOS(w)+DCOS(omega)*DSIN(w)*COS(i) ROT(2,2)=-1.0D0*DSIN(omega)*DSIN(w)+DCOS(omega)*DCOS(w)*DCOS(i) ROT(2,3)=-1.0D0*DCOS(omega)*DSIN(i)
ROT(3,1)=DSIN(w)*DSIN(i) ROT(3,2)=DCOS(w)*DSIN(i) ROT(3,3)=DCOS(i)
rx(1)= rx1(1)*ROT(1,1) + rx1(2)*ROT(1,2) + rx1(3)*ROT(1,3)
rx(2)= rx1(1)*ROT(2,1) + rx1(2)*ROT(2,2) + rx1(3)*ROT(2,3)
rx(3)= rx1(1)*ROT(3,1) + rx1(2)*ROT(3,2) + rx1(3)*ROT(3,3)
vx(1)= vx1(1)*ROT(1,1) + vx1(2)*ROT(1,2) + vx1(3)*ROT(1,3)
vx(2)= vx1(1)*ROT(2,1) + vx1(2)*ROT(2,2) + vx1(3)*ROT(2,3)
vx(3)= vx1(1)*ROT(3,1) + vx1(2)*ROT(3,2) + vx1(3)*ROT(3,3)
do 2 j = 1,3 x(j) = rx(j)
2 x(j+3) = vx(j)
RETURN END
SUBROUTINE Ekepler (e,aM,aE) implicit real*8 (a-h,o-z) common/const0/pi,tol,iflag,itarget,idim ttol = 0.d-8
IF (-pi .LT. aM .AND. aM .LT. 0.d0) aE = aM -e IF (0.d0 .LT. aM .AND. aM .LT. pi) aE = aM + e
IF (aM .GT. pi) aE = aM – e vary = 100.d0
do 2 i = 1,50
bE = aE +(aM – aE + e*dsin(aE))/(1.d0+e*dcos( aE)) aE = bE
IF (ABS(aE-bE) .LT. ttol) GOTO 3
SUBROUTINE conjunction (t,f0e,fe,fm,xMe,xMm) implicit real*8 (a-h,o-z)
common/output/tj common/ephemers/ej common/constopt/alpha,beta,t1,t2 common/const0/pi,tol,iflag,itarget,idim
common/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE common/constm2/xm, gmmars, radM,rm,soiM,omm,phM,rpM dimension xe(6),xm(6),ej(4,6),tj(30,3)
c
pi2 = 2.d0*pi
xM0e = ej(3,6) – ej(3,4) – ej(3,5)
xM0m = ej(4,6) – ej(4,4) – ej(4,5)
period = 1.d0 t = 0.d0
test = 2
do 2 i = 1,1000
t = t + 86400d0*period xMe = xM0e + ome*t xMm = xM0m + omm*t
CALL Ekepler (ej(3,2),xMe,Ee) fe=2d0*datan(dsqrt((1d0+ej(3,2))/(1d0-ej(3,2)))*tan(Ee/2d0)) CALL Ekepler (ej(4,2),xMm,Em) fm=2d0*datan(dsqrt((1d0+ej(4,2))/(1d0-ej(4,2)))*tan(Em/2d0)) f0e = dmod(fe + ej(3,4) + ej(3,5),pi2)
f0m = dmod(fm + ej(4,4) + ej(4,5),pi2)
IF (ABS(f0e – f0m) .LT. 1.d0 .AND. f0e .GT. f0m) THEN
t = t – 86400.d0*period period = period/2.d0
IF (ABS(f0e – f0m) .LT. 1.d-12) GOTO 3
ENDIF
-
- fe = dmod(fe,pi2) xMe = dmod(xMe,pi2)
pe = ej(3,1)*(1d0 – ej(3,2)**2)
CALL randv (gmsun,pe,ej(3,2),ej(3,3),ej(3,4),ej(3,5),xMe,xe)
IF (alpha .NE. 0.d0) THEN
fm = dmod(fm + alpha,pi2)
Em=2d0*datan(dsqrt((1d0-ej(4,2))/(1d0+ej(4,2)))*tan(fm/2d0)) xMm = Em – ej(4,2)*dsin(Em)
xMm = dmod(xMm,pi2) ELSE
fm = dmod(fm,pi2)
xMm = dmod(xMm,pi2) ENDIF
pm = ej(4,1)*(1d0 – ej(4,2)**2)
CALL randv (gmsun,pm,ej(4,2),ej(4,3),ej(4,4),ej(4,5),xMm,xm) tj(1,1) = t/86400.d0
tj(2,1) = f0e*180d0/pi tj(2,2) = fe*180d0/pi tj(2,3) = fm*180d0/pi RETURN
END
*——— acknowledgements ————————–
*
* The following subroutines were provided by
* Dr. Rodger Broucke, Professor Emeritus, Univ. of Texas, Austin
* used with his permission
*
* INRK78, RK78T, perturb3b, energye, deriv
*
*—————————————————-
SUBROUTINE gui
implicit real*8 (a-h,o-z)
print *, ‘ Mars Pathfinder ‘ print *, ‘ by’
print *, ‘ Bar X Software – Austin, Texas’ print *, ‘ Est. 1993’
print *, ‘——————————————————— ‘
print *, ‘.’
print *, ‘ Your mission is to reach Mars by adjusting values’ print *, ‘ of the trajectory parameters to (1) minimize’
print *, ‘ fuel usage, in km/sec of thrust, and (2) to’
print *, ‘ reach Mars as quickly as possible……… ‘
print *, ‘ You can change the trajectory in one place, ‘
print *, ‘ about half way to Mars. You select MAGNITUDE and’ print *, ‘ DIRECTION for this thrust. The program will apply’ print *, ‘ the thrust at the best time and place, and’
print *, ‘ optimize all the other four maneuvers of the’ print *, ‘ flight path based upon your criteria.’
print *, ‘.’
print *, ‘ Now enter the THRUST MAGNITUDE in km/sec.’
print *, ‘ Values should be between zero’
print *, ‘ (a standard Hohmann Transfer) and 2.5.’ print *, ‘.’
read *, dVmars print *, ‘.’
print *, ‘ Now the THRUST DIRECTION, in degrees above/below’ print *, ‘ the current velocity direction. Be careful! A value’
print *, ‘ below zero, with a low thrust and you will never reach’ print *, ‘ Mars because you will pass inside its orbit!’
print *, ‘.’ read *, beta1 print *, ‘.’
print *, ‘ The program will now compute the minimum energy’ print *, ‘ trajectory to get you to Mars in the least time,’
print *, ‘ using the values you have specified.’
print *, ‘ HEY!’
print *, ‘ If the program crashes – dont blame it!!’
print *, ‘ Try another scenario – more thrust, 0 angle… ‘
c
dVmars = -1.d0*dVmars
CALL catapult (beta1,beta2,dVmars,deltaVf,tsip) call XOUT
c
print *, ‘——————————————————- ‘
print *, ‘ Your total thrust requirement in km/sec was’ print *, deltaVf
print *, ‘ The total time required for this mission in days was,’ print *, tsip
print *, ‘ Please refer to the above data printout for a more’ print *, ‘ complete mission profile.’
print *, ‘——————————————————- ‘
print *, ‘ Copyright 2000, William H. Clark II’ print *, ‘ All Rights Reserved.’
print *, ‘.’
print *, ‘——————————————————- ‘
print *, ‘——————————————————- ‘
print *, ‘ Please enter 1 if you want to try another mission.’ print *, ‘.’
read *, repeat
IF (repeat .EQ. 1) THEN
-
- print *, ‘ New thrust,’ read *, dVmars
print *, ‘ New angle,’ read *, beta1
print *, ‘……. computing!’
dVmars = -1.d0*dVmars
CALL catapult (beta1,beta2,dVmars,deltaVf,tsip) call XOUT(LP)
print *, ‘.’
print *, ‘ Total thrust:’ print *, deltaVf
print *, ‘ Total time:’ print *, tsip
print *, ‘.’
print *, ‘ Enter 1 to try again!’ print *, ‘.’
read *, repeat
IF (repeat .EQ. 1) THEN GOTO 8
ELSE
GOTO 9 ENDIF ELSE GOTO 10 ENDIF
SUBROUTINE xout
implicit real*8 (a-h,o-z) dimension tj(30,3) common/output/tj
c
print *, ‘.’
print *, ‘——————————————————- ‘
print *, ‘ Mission Profile’
print *, ‘——————————————————- ‘
print *, ‘.’
print *, ‘ The program uses the J2000 Ephemerides. The first’ print *, ‘ Earth-Mars conjunction occurs, in days after’
print *, ‘ 1/1/2000 (midnight is at 0.5 days) at’ print *, tj(1,1)
print *, ‘ The position of Earth and Mars at this time is’ print *, ‘ (units, degrees)’
print *, tj(2,1)
print *, ‘ Where zero degrees is the Vernal Equinox by’ print *, ‘ convention. All subsequent coordinates will be in’
print *, ‘ this coordinate system, with the z-axis perpendicular’ print *, ‘ to the orbit plane. The true anomalies of Earth and’ print *, ‘ Mars at mid-course correction (about conjunction) are’ print *, tj(2,2)
print *, tj(2,3)
print *, ‘ The motion is assumed to be co-planar.’ print *, ‘——————————– ‘
print *, ‘.’
print *, ‘ The mission begins with the spacecraft (s/c) in a’
print *, ‘ 200 km parking orbit around Earth. The s/c maneuvers’ print *, ‘ to slightly higher circular orbit using a minimum’
print *, ‘ energy Hohmann transfer with a thrust of’ print *, ‘ (units, km/sec)’
print *, tj(3,1)
print *, ‘ The interplanetary trajectory begins with a thrust of’ print *, tj(3,3) + tj(3,2)
print *, ‘ at the following geocentric (earth-centered)’
print *, ‘ coordinates, initiatiating a hyperbolic escape’ print *, ‘ trajectory (perturbed by the sun and mars):’
print *, tj(4,1) print *, tj(4,2)
print *, ‘ The initial velocity coordinates are:’
print *, tj(5,1) print *, tj(5,2)
print *, ‘ The time (relative to the mid-course thrust) is’ print *, tj(6,1)
print *, ‘ The s/c leaves Earths sphere of influence (soi)’ print *, ‘ at the following geocentric coordinates and time:’ print *, tj(7,1)
print *, tj(7,2) print *, tj(8,1) print *, tj(8,2) print *, tj(6,2)
print *, ‘ Its heliocentric (sun-centered) coordinates’ print *, ‘ at this time are’
print *, tj(9,1) print *, tj(9,2) print *, tj(10,1) print *, tj(10,2)
print *, ‘——————————————————- ‘
print *, ‘.’
print *, ‘ The s/c is now on an elliptical heliocentric course’ print *, ‘ (perturbed by Earth, Mars, and Jupiter). The mid-‘ print *, ‘ course correction occurs at time zero, and at the’ print *, ‘ following heliocentric coordinates’
print *, tj(11,1) print *, tj(11,2) print *, tj(12,1) print *, tj(12,2)
print *, ‘ (This velocity is after the thrust is applied.)’ print *, ‘ The s/c reaches Mars soi at the heliocentric’ print *, ‘ coordinates and time of,’
print *, tj(13,1) print *, tj(13,2) print *, tj(14,1) print *, tj(14,2) print *, tj(15,1)
print *, ‘ (Time, measured after the mid-course corrrection.)’ print *, ‘ The mars-centered (geocentric) coordinates at this’ print *, ‘ point are’
print *, tj(16,1) print *, tj(16,2) print *, tj(17,1) print *, tj(17,2)
print *, ‘——————————————————- ‘
print *, ‘.’
print *, ‘ The s/c follows a hyperbolic capture orbit to Mars’ print *, ‘ (perturbed by the sun), and applies a thrust of’ print *, tj(21,3) + tj(21,1)
print *, ‘ to circularize the orbit at coordinates, and time,’ print *, tj(18,1)
print *, tj(18,2) print *, tj(19,1)
print *, tj(19,2) print *, tj(20,1)
print *, ‘ A final thrust, used to target a specific latitude’
print *, ‘ and longitude on the surface of Mars, of’ print *, tj(21,2)
print *, ‘ brings the s/c via Hohmann Transfer into a 100 km’ print *, ‘ parking orbit around Mars. The total time and’ print *, ‘ thrust required for this scenario are’
print *, tj(22,1) print *, tj(22,2) print *,’.’ RETURN END
