June 2000 Mars FORTRAN (copyrighted) – Grad Student Advocate

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

continue
bx = ax – step*dX step = step/2.d0 GOTO 2
RETURN END

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)
RETURN END

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

XD(I)=X(I)

GO TO 1020 C 1010 T=TM DT1=DT DO 41 I=1,N

X(I)=XD(I)

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

continue
RETURN END

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

continue

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

close (LP)
RETURN END

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