The Mars Pathfinder Source Code

C Last change: 15 Nov 2003 12:43 pm
c "Grabet" Orbits from Earth to Mars
c ~ Gravity Assisted Bi-Elliptic Transfer Orbits
c copyright 2003 sahmain space systems
c
implicit real*8 (a-h,o-z)
common/RKCOM/CH(13),AL(13),B(13,12)
CALL INRK78
CALL Earth (deltaV1)
thrust = 0.0d0
CALL Mars (thrust,deltaV2)
PRINT *, deltaV1 + deltaV2 ! = 0.24 km/sec < simple Hohmann
c Apply this thrust @ f=90 and arrive 30 days < simple Hohmann
STOP
END
c
SUBROUTINE Mars (thrust,dVmin)
implicit real*8 (a-h,o-z)
REAL*8 x(6),xe(6),xf(6),xm(6),xg(6),xtable(40,14)
common/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE
common/constm2/xm, gmmars, radM,rm,soiM,omm,phM,rpM
pi = DACOS(-1.0d0)
CALL initial (t,x)
theta =-DATAN(x(5)/x(4))
c x(4) = x(4) - thrust*DCOS(theta)
c x(5) = x(5) + thrust*DSIN(theta)
tol = 1.0d-13
dt = 1000d0
do 1 n = 1,1500
CALL RK78T (t,x,dt,tol,6)
rangeM = DSQRT((xm(1)-x(1))**2+(xm(2)-x(2))**2)
IF (rangeM .LT. soiM*2) dt =1000.0d0
IF (rangeM .LT. soiM) GOTO 2
IF (dsqrt(x(1)**2+x(2)**2) .GT. rm) GOTO 2
1 continue
c transfer to geocentric coordinates
2 xg(1) = x(1) - xm(1)
xg(2) = x(2) - xm(2)
xg(3) = 0.0d0
xg(4) = x(4) + DSQRT(gmsun/rm)*DSIN(t*omm + phM)
xg(5) = x(5) - DSQRT(gmsun/rm)*DCOS(t*omm + phM)
xg(6) = 0.d0
soi= soiM
soiM = dsqrt(xg(1)**2 + xg(2)**2)
Vrel = DSQRT(xg(4)**2+xg(5)**2)
Rf = radM + 80.0d0 ! 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
tmax = 10.d0*soiM/Vrel
c
c do 78 m = 1,40
m = 27 ! tried all 40 cases; best results with this one
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)
rangeM = Dsqrt(xg(1)**2 + xg(2)**2 + xg(3)**2)
IF (rangeM .LT. radM) THEN
deltaV = 100.0d0
GOTO 78
ENDIF
CALL velocityM (xg,gmsun,rangeM,Rf,Vrel,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) = rangeM
dVmin = deltaV
xtable(m,11) = V1
xtable(m,12) = V2
xtable(m,13) = V3
xtable(m,14) = t/86400d0
ENDIF
IF (rangeM .GT. soiM*1.1d0) GOTO 78
IF (t .GT. tmax) GOTO 78
75 continue
78 continue
RETURN
END
c
SUBROUTINE Earth (deltaV)
implicit real*8 (a-h,o-z)
REAL*8 x(6), xf(10)
step = 0.0001d0
ttol = 1.0d-7
ax = -.005d0
adelta = 100.0d0
CALL initial (t,x)
CALL TargetE (t,x,ax,adelta,xf)
bdelta = adelta
delta = adelta
2 DO 3 k = 1,3
IF (step .LT. ttol) GOTO 5
bx = ax + step
CALL initial (t,x)
CALL TargetE (t,x,bx,bdelta,xf)
IF (bdelta .LT. delta) delta = bdelta
IF (bdelta .GT. adelta) GOTO 4
ax = bx
adelta = bdelta
3 continue
4 bx = bx - step
step = step/2.d0
GOTO 2
5 deltaV = delta + bx
RETURN
END
c
SUBROUTINE TargetE (t,x,deltaV,deltaVf,xf)
implicit real*8 (a-h,o-z)
REAL*8 x(6),xe(6),xf(6),xm(6),xg(6)
common/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE
common/constm2/xm, gmmars, radM,rm,soiM,omm,phM,rpM
theta = DATAN(x(5)/x(4))
x(4) = x(4) + deltaV*DCOS(theta)
x(5) = x(5) + deltaV*DSIN(theta)
tol = 1.0d-13
dt = -1000d0
do 75 n = 1,1500
CALL RK78T (t,x,dt,tol,6)
rangeE = DSQRT((xe(1)-x(1))**2+(xe(2)-x(2))**2)
IF (rangeE .LT. soiE) GOTO 78
75 continue
c change to geocentric coordinates (xe + xf = x)
78 xg(1) = x(1) - xe(1)
xg(2) = x(2) - xe(2)
xg(3) = 0.0d0
xg(4) = x(4) + DSQRT(gmsun/re)*DSIN(t*ome)
xg(5) = x(5) - DSQRT(gmsun/re)*DCOS(t*ome)
xg(6) = 0.d0
gm = gmsun
gmm = gmmars
gmsun = gmearth
gmearth = gm
re = -re
gmmars = 0.0d0
Rf = radE + 200.0d0
dt = -100.0d0
deltaVf = 100.0d0
do 85 n = 1,85
CALL RK78T (t,xg,dt,tol,6)
rangeE = DSQRT(xg(1)**2+xg(2)**2)
IF (rangeE .LT. radE) GOTO 88
IF (rangeE .GT. soiE*1.2) GOTO 88
Vrel = DSQRT(xg(4)**2+xg(5)**2)
IF (xg(2) .LT. 0.0d0) THEN
CALL VelocityE (xg,gmsun,rangeE,Rf,Vrel,DeltaVx,V1,V2,V3)
IF (deltaVx .LT. deltaVf) deltaVf = deltaVx
END IF
85 continue
88 gmearth = gmsun
gmsun = gm
gmmars = gmm
re = -re
RETURN
END
c
SUBROUTINE velocityM (xg,gmsun,range,Rf,Vrel,deltaV,V1,V2,V3)
implicit real*8 (a-h,o-z)
REAL*8 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 velocityE (xg,gmsun,range,Rf,Vrel,deltaV,V1,V2,V3)
implicit real*8 (a-h,o-z)
REAL*8 xg(6)
Rp = range
Vcirc = dsqrt(gmsun/range)
V3 = dabs(Vcirc - Vrel)
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/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE
common/constm2/xm, gmmars, radM,rm,soiM,omm,phM,rpM
real*8 x3(6),xe(6),xm(6)
pi = DACOS(-1.0d0)
IF (n .EQ. 3) THEN ! position of Earth
x3(1) = re*DCOS(t*ome)
x3(2) = re*DSIN(t*ome)
x3(3) = 0.0d0
xe(1) = x3(1)
xe(2) = x3(2)
xe(3) = x3(3)
ENDIF
IF (n .EQ. 4) THEN ! position of Mars
x3(1) = rm*DCOS(t*omm + phM)
x3(2) = rm*DSIN(t*omm + phM)
x3(3) = 0.0d0
xm(1) = x3(1)
xm(2) = x3(2)
xm(3) = x3(3)
ENDIF
RETURN
END
c
SUBROUTINE initial (t,x)
implicit real*8 (a-h,o-z)
real*8 x(6),xe(6),xm(6),M,i
common/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE
common/constm2/xm, gmmars, radM,rm,soiM,omm,phM,rpM
c input the planetary constants
gmsun = 1.32712428d11
gmearth = 3.986004415d5
re = 149598023d0 ! = a
ome = dsqrt((gmsun+gmearth)/re)/re
soiE = 924647d0
radE = 6378.1363d0
rpE = 200.0d0
gmmars = .305d4
rm = 227939186d0 ! = a
omm = dsqrt((gmsun+gmmars)/rm)/rm
soiM = 577213d0
radM = 3397.2d0
rpM = 80.0d0
c determine the Hohmann Transfer orbit paramenters
pi = DACOS(-1.0d0)
a = (re + rm)/2.d0
e = (rm - re)/(re + rm)
p = a*(1.d0 - e**2)
period = pi*(dsqrt(a)**3)/dsqrt(gmsun)
omsc = dsqrt((gmsun)/a)/a
tconj = (pi - period*omm)/(ome-omm)
phM = pi - period*omm
c start when f=90 for the spacecraft
EA = 2.0d0*DATAN(DTAN(pi/4)*DSQRT((1.0d0-e)/(1.0d0+e)))
M = EA -e*DSIN(EA)
t = M/omsc
x(1) = 0.0d0
x(2) = p
x(3) = 0.0d0
x(4) = -1.0d0*DSQRT(gmsun/p)
x(5) = DSQRT(gmsun/p)*e
x(6) = 0.0d0
RETURN
END
c
SUBROUTINE DERIV(t,x,f)
implicit real*8 (a-h,o-z)
real*8 x(6),f(6),xe(6),xm(6)
common/conste2/xe,gmsun,gmearth,radE,re,soiE,ome,phE,rpE
common/constm2/xm, gmmars, radM,rm,soiM,omm,phM,rpM
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)
RETURN
END

SUBROUTINE energye(x,t,E)
implicit real*8 (a-h,o-z)
real*8 x(6),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))
E = E2body -Um1 - Um2 + Um3
RETURN
END
C
SUBROUTINE pert3b(GMP,x3,x,f)
implicit real*8 (a-h,o-z)
real*8 f(6),x(6),x3(6),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 REAL*8 (A-H,O-Z)
common/RKCOM/CH(13),AL(13),B(13,12)
REAL*8 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
40 XD(I)=X(I)
GO TO 1020
C
1010 T=TM
DT1=DT
DO 41 I=1,N
41 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 REAL*8 (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