## Chapter 3Astrodynamics

### 3.2 Table of common equations

The following table contains the common relations to use for elliptic motion. Equation of ellipse is $$\frac {x^{2}}{a^{2}}+\frac {y^{2}}{b^{2}}=1$$

 term to ﬁnd relation conversion between $$E$$ and $$\theta$$ \begin {aligned} \tan \left (\frac {\theta }{2}\right ) &=\sqrt {\frac {1+e}{1-e}}\tan \left ( \frac {E}{2}\right ) \\ \cos E &= \frac {e+\cos \theta }{1+e\cos \theta }\\ \cos \theta &= \frac {e-\cos E}{e\cos E-1} \end {aligned} position of satellite at time $$t$$ Solve for $$E$$, then ﬁnd $$\theta$$. $$\tau$$ here is time at $$perigee$$ and $$n$$ is mean satellite speed. $$E-e\sin E=n\left ( t-\tau \right )$$ eccentricity $$e$$ $$e=\frac {c}{a}=\frac {r_{a}-r_{p}}{r_{a}+r_{p}} =\sqrt {1+\frac {2\mathcal {E}h^{2}}{\mu ^{2}}}$$ Major axes $$a$$ \begin {aligned} a &= \frac {r_p (1+e)}{1-e^2} \\ &= \frac {r_a (1-e)}{1-e^2} \\ &= -\frac {\mu }{2\mathcal {E}}\\ &= \sqrt {b^{2}+c^{2}}\\ &= \frac {p}{1-e^{2}} \end {aligned} Minor axes $$b$$ \begin {aligned} b &=a\sqrt {1-e^{2}} \end {aligned} $$r_p$$ \begin {aligned} r_p &= \frac {a (1-e^2)}{1+e}\\ &= a(1-e) \end {aligned} $$r_a$$ \begin {aligned} r_a &= \frac {a (1-e^2)}{1-e}\\ &= a\left ( 1+e\right ) \\ &=\frac {h^{2}}{\mu }\frac {1}{1-e} \end {aligned} $$\frac {r_{p}}{r_{a}} =\frac {1-e}{1+e}$$ $$p$$ $$p=a\left ( 1-e^{2}\right ) =\frac {h^{2}}{\mu }=r_{p}\left ( 1+e\right ) =r_{a}\left ( 1-e\right )$$ speciﬁc angular momentum $$h$$ \begin {aligned} h &= r_{p}v_{p}=r_{a}v_{a}=\vec {r}\times \vec {v}=\sqrt {p\mu }\\ h &=\sqrt {\mu r} \quad \text {(circular orbit)} \end {aligned} Total Energy $$\mathcal {E}$$ $$\mathcal {E}=\frac {v^{2}}{2}-\frac {\mu }{r} =-\frac {\mu }{2a}$$ velocity $$v$$ \begin {aligned} v &= \sqrt {\mu \left ( \frac {2}{r}-\frac {1}{a}\right ) } \quad \text {(vis-viva)}\\ v_{escape} &=\sqrt {\frac {2\mu }{r}} \quad \text {(escape velocity for parabola)}\\ v_{radial} &=\sqrt {\frac {\mu }{p}}e\sin \theta \\ v_{normal} &=\sqrt {\frac {\mu }{p}}\left ( 1+e\cos \theta \right ) \end {aligned} $$v_{perigee}$$ (closest) \begin {aligned} v_p &= \sqrt {\frac {\mu }{a}\left ( \frac {1+e}{1-e}\right ) }\\ &= \sqrt {\frac {\mu }{p}}(1+e)\\ &= \sqrt {\mu \left (\frac {2}{r_p} - \frac {1}{a} \right )} \end {aligned} $$v_{apogee}$$ (furthest) \begin {aligned} v_{a} &= \sqrt {\frac {\mu }{a}\left ( \frac {1-e}{1+e}\right ) }\\ &= \sqrt {\frac {\mu }{p}}(1-e) \\ &= \sqrt {\mu \left (\frac {2}{r_a} - \frac {1}{a} \right )} \end {aligned} magnitude of $$\vec {r}$$ \begin {aligned} r &= \frac {a\left ( 1-e^{2}\right ) }{1+e\cos \theta }=\frac {h^{2}}{\mu }\frac {1}{1+e\cos \theta }\\ r\cos \theta &= a\left ( \cos E-e\right ) \\ r &= a\left ( 1-e\cos E\right ) \quad \text {(eq 4.2-14 Bate book)} \end {aligned} period $$T$$ $$T = \frac {2}{h}\pi ab=2\pi \sqrt {\frac {a^{3}}{\mu }}$$ mean satellite speed $$n$$ $$n=\frac {2\pi }{T}=\sqrt {\frac {\mu }{a^{3}}}$$ eccentric anomaly $$E$$ $$\tan \frac {\theta }{2}=\sqrt {\frac {1+e}{1-e}}\tan \frac {E}{2}$$ area sweep rate $$\frac {dA}{dt}=\frac {h}{2}$$ equation of motion $$\ddot {\vec {r}}+\frac {\mu }{r^{3}}\vec {r}=0$$ spherical coordinates relation $$\cos (i) = \sin (A_z) \cos (\phi )$$ where $$i$$ is the inclination and $$A_z$$ is the azimuth and $$\phi$$ is latitude 1

Notice in the above, that the period $$T$$ of satellite depends only on $$a$$ (for same $$\mu$$)

In the above, $$\mu =GM$$ where $$M$$ is the mass of the body at the focus of the ellipse and $$G$$ is the gravitational constant. $$h$$ is the speciﬁc mass angular momentum (moment of linear momentum) of the satellite. Hence the units of $$\frac {h^{2}}{\mu }$$ is length.

To draw the locus of the satellite (the small body moving around the ellipse, all what we need is the eccentricity $$e$$ and $$a$$, the major axes length. Then by changing the angle $$\theta$$ the path of the satellite is drawn. I have a demo on this here

See http://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html for earth facts

This table below is from my class EMA 550 handouts (astrodynamics, spring 2014)

### 3.3 Flight path angle for ellipse $$\gamma$$

\begin {align*} v & =\left \vert \boldsymbol {\vec {v}}\right \vert =\sqrt {\left \vert \boldsymbol {\vec {v}}_{r}\right \vert ^{2}+\left \vert \boldsymbol {\vec {v}}_{n}\right \vert ^{2}}\\ & =\sqrt {\mu \left ( \frac {2}{r}-\frac {1}{a}\right ) } \end {align*}

To ﬁnd $$\gamma$$, if $$r$$ is given then use

$\cos \gamma =\sqrt {\frac {ap}{r\left ( 2a-r\right ) }}=\sqrt {\frac {a^{2}\left ( 1-e^{2}\right ) }{r\left ( 2a-r\right ) }}$

If $$\theta$$ is given, then use

$\tan \gamma =\frac {e\sin \theta }{1+e\cos \theta }$

### 3.4 Parabolic trajectory

This diagram shows the parabolic trajectory

### 3.5 Hyperbolic trajectory

This diagram shows the hyperbolic trajectory

This diagram below from Orbital mechanics for Engineering students, second edition, by Howard D. Curtis, page 109

### 3.6 showing that energy is constant

Showing that energy $$\mathcal {E}=\frac {v^{2}}{2}-\frac {\mu }{r}$$ is constant.

Most of such relations starts from the same place. The equation of motion of satellite under the assumption that its mass is much smaller than the mass of the large body (say earth) it is rotating around. Hence we can use $$\nu =GM$$ and the equation of motion reduces to $\ddot {\vec {r}} + \frac {\mu }{r^{3}} \vec {r} = 0$ In the above equation, the vector $$\vec {r}$$ is the relative vector from the center of the earth to the center of the satellite. The reason the center of earth is used as the origin of the inertial frame of reference is due to the assumption that $$M\gg {m}$$ where $$M$$ is the mass of earth (or the body at the focal of the ellipse) and $$m$$ is the mass of the satellite. Hence the median center of mass between the earth and the satellite is taken to be the center of earth. This is an approximation, but a very good approximation.

The ﬁrst step is to dot product the above equation with $$\dot {\vec {r}}$$ giving

\begin {align} \dot {\vec {r}} \cdot \ddot {\vec {r}} + \dot {\vec {r}} \cdot \frac {\mu }{r^{3}} \vec {r} & = 0 \tag {1} \end {align}

And there is the main trick. We look ahead and see that $$\dot {\vec {r}}\cdot \ddot {\vec {r}}=\dot {r}\ddot {r}$$ but $$\dot {r}\ddot {r}=\frac {d}{dt}\left ( \frac {\dot {r}^{2}}{2}\right )$$ and we also see that $$\dot {\vec {r}}\cdot \frac {\mu }{r^{3}}\vec {r}=\mu \frac {\dot {r}}{r^{2}}$$ but $$\mu \frac {\dot {r}}{r^{2}}=\frac {d}{dt}\left ( \frac {-\mu }{r}\right )$$ Hence equation 1 above can be written as $\frac {d}{dt}\left ( \frac {v^{2}}{2}-\frac {r}{\mu }\right ) =0$ Hence $\mathcal {E}=\frac {v^{2}}{2}-\frac {r}{\mu }$ Where $$\mathcal {E}$$ is a constant, which is the total energy of the satellite.

### 3.7 Earth satellite Transfer orbits

#### 3.7.1 Hohmann transfer

This diagram shows the Hohmann transfer

### 3.8 Rocket engines, Hohmann transfer, plane change at equator

two cases: Hohmann transfer, 2 burns, or semi-tangential. All burns at equator.

### 3.10 interplanetary transfer orbits

##### 3.10.0.1 interplanetary hohmann transfer orbit, case one

The following are the steps to accomplish the above. The ﬁrst stage is getting into the Hohmann orbit from planet 1, then reaching the sphere of inﬂuence of the second planet. Then we either do a ﬂy-by or do a parking orbit around the second planet. These steps below show how to reach the second planet and do a parking orbit around it.

The input is the following.

1. $$\mu _{1}$$ planet one standard gravitational parameter
2. $$\mu _{2}$$ planet two standard gravitational parameter
3. $$\mu _{sun}$$ standard gravitational parameter for the sun $$1.327\times 10^{8}$$ km
4. $$r_{1}$$ planet one radius
5. $$r_{2}$$ planet two radius
6. $$alt_{1}$$ original satellite altitude above planet one. For example, for LEO use 300 km
7. $$alt_{2}$$ satellite altitude above second planet. (since goal is to send satellite for circular orbit around second planet)
8. $$R_{1}$$ mean distance of center of ﬁrst planet from the sun. For earth use $$AU=1.495978\times 10^{8}$$ km
9. $$R_{2}$$ mean distance of center of second planet from the sun. For Mars use $$1.524$$ $$AU$$
10. $$SOI_{1}$$ sphere of inﬂuence for ﬁrst planet. For earth use $$9.24\times 10^{8}$$ km
11. $$SOI_{2}$$ sphere of inﬂuence for second planet.

Given the above input, there are the steps to achieve the above maneuver

1. Find the burn out distance of the satellite $$r_{bo}=r_{1}+alt_{1}$$
2. Find satellite speed around planet earth (relative to planet) $$V_{sat}=\sqrt {\frac {\mu _{1}}{r_{bo}}}$$
3. Find Hohmann ellipse $$a=\frac {R_{1}+R_{2}}{2}$$
4. Find speed of satellite at perigee relative to sun $$V_{perigee}=\sqrt {\mu _{sun}\left ( \frac {2}{R_{1}}-\frac {1}{a}\right ) }$$
5. Find speed of earth (ﬁrst planet) relative to sun $$V_{1}=\sqrt {\frac {\mu _{sun}}{R_{1}}}$$
6. Find escape velocity from ﬁrst planet $$V_{\infty ,out}=V_{perigee}-V_{1}$$
7. Find burn out speed at ﬁrst planet by solving the energy equation $$\frac {V_{bo}^{2}}{2}-\frac {\mu _{1}}{r_{bo}}=\frac {V_{\infty ,out}^{2}}{2}-\frac {\mu _{1}}{SOI_{1}}$$for $$V_{bo}$$
8. Find $$\Delta V_{1}$$ needed at planet one $$\Delta V_{1}=V_{bo}-V_{sat}$$
9. Find $$e$$ the eccentricity of the escape hyperbola $$e=\sqrt {1+\frac {V_{\infty }^{2}V_{bo}^{2}r_{bo}^{2}}{\mu _{1}^{2}}}$$
10. Find the angle with the path of planet one velocity vector $$\eta =\arccos \left ( -\frac {1}{e}\right )$$
11. Find the dusk-line angle $$\theta =180^{0}-\eta$$

The above completes the ﬁrst stage, now the satellite is in the Hohmann transfer orbit. Assuming it reached the orbit of the second planet ahead of it as shown in the diagram above. Now we start the second stage to land the satellite on a parking orbit around the second planet at altitude $$alt_{2}$$ above the surface of the second planet. These are the steps needed.

1. Find the apogee speed of the satellite $$V_{apogree}=\sqrt {\mu _{sun}\left ( \frac {2}{R_{2}}-\frac {1}{a}\right ) }$$
2. Find speed of second planet $$V_{2}=\sqrt {\frac {\mu _{sun}}{R_{2}}}$$
3. Find $$V_{\infty }$$ entering the second planet sphere of inﬂuence $$V_{\infty ,in}=V_{2}-V_{apogree}$$
4. Find burn in radius where the satellite will be closest to the second planet. $$r_{bo}=r_{1}+alt_{2}$$
5. Find burn out speed at second planet by solving the energy equation $$\frac {V_{bo}^{2}}{2}-\frac {\mu _{2}}{r_{bo}}=\frac {V_{\infty ,in}^{2}}{2}-\frac {\mu _{2}}{SOI_{2}}$$for $$V_{bo}$$
6. Find impact parameter $$b$$ on entry to second planet SOI $$b=\frac {r_{bo}V_{bo}}{V_{\infty ,in}}$$
7. Find the required satellite speed around the second planet $$V_{sat}=\sqrt {\frac {\mu _{2}}{r_{bo}}}$$
8. Find $$\Delta V_{2}$$ needed at planet two $$\Delta V_{2}=V_{sat}-V_{bo}$$
9. Find $$e$$ the eccentricity of the approaching hyperbola on second planet $$e=\sqrt {1+\frac {V_{\infty }^{2}V_{bo}^{2}r_{bo}^{2}}{\mu _{2}^{2}}}$$
10. Find the angle with the path of planet two velocity vector $$\eta =\arccos \left ( -\frac {1}{e}\right )$$
11. Find the dusk-line angle for second planet $$\theta =2\eta -180^{0}$$

#### 3.10.1 rendezvous orbits

##### 3.10.1.1 Two satellite, walking rendezvous using Hohmann transfer
1:   function hohmann_walking_rendezvous($$\theta _0,r,altitude,\mu$$)
2:    $$\theta _0 := \theta _0 \frac {\pi }{180}$$ $$\triangleright$$ convert from degrees to radian
3:    $$r_a := r+\text {altitude}$$
4:    $$T := 2\pi \sqrt {\frac {r_a^3}{\mu }}$$ $$\triangleright$$ period of circular orbit
5:    $$N:=1$$
6:    done:=false
7:    while not(done) do
8:    TOF := $$\left (N-\frac {\theta _0}{2\pi }\right )T$$
9:    $$a := \textbf {solve}\, \left (TOF = 2\pi \sqrt {\frac {a^3}{\mu }}\right ) \text {for}\, a$$
10:    $$r_p:=2 a-r_a$$
11:    if $$rp<r$$ then
12:    $$N=N+1$$
13:    else
14:    done:=true
15:    end if
16:    end while
17:    $$V_{befor}:=\sqrt {\frac {\mu }{h}}$$
18:    $$V_{after}:=\sqrt {\mu \left ( \frac {2}{h} - \frac {1}{a} \right )}$$
19:    $$\Delta V:= 2(V_{after}-V_{before})$$
20:    return $$(\text {TOF},\Delta V)$$
21:   end function

An example implementation is below

hohmannRendezvousSameOrbit[\[Theta]00_, r_, alt_, mu_] :=
Module[{\[Theta]0 = \[Theta]00*Pi/180, n = 1, delT, v1, v2, period, a,
rp, ra, done = False, vBefore, vAfter},
ra     = r + alt;
period = 2 Pi Sqrt[ra^3/mu];
While[Not[done],
delT = (n - \[Theta]0 /(2 Pi)) period ;
a = First@Select[a /. NSolve[delT == 2 Pi Sqrt[a^3/mu], a],
Element[#, Reals] &];
rp = 2 a - ra;
If[rp < r,(*we hit the earth, try again*)
n = n + 1,
done = True
]
];
vBefore = Sqrt[mu/h];
vAfter  = Sqrt[mu (2/h - 1/a)];
{delT, 2 (vAfter - vBefore)}  (*return value*)
]



And calling the above

mu = 324859;
alt = 1475.776;
r = 6052;
\[Theta]0 = 3.80562; (*degree*)
hohmannRendezvousSameOrbit[\[Theta]0, r, alt, mu]



gives

{7123.89, -0.0467913}


##### 3.10.1.2 Two satellite, separate orbits, rendezvous using Hohmann transfer, coplaner
1:   function hohmann_rendezvous_1($$\theta _0,r_a,r_b,\mu$$)
2:    $$\theta _0 := \theta _0 \frac {\pi }{180}$$ $$\triangleright$$ convert from degrees to radian
3:    $$a := \frac {r_a+r_b}{2}$$ $$\triangleright$$ Hohmann orbit semi-major axes
4:    $$\text {TOF} := \pi \sqrt {\frac {a^3}{\mu }}$$ $$\triangleright$$ time of ﬂight on Hohmann orbit
5:    $$\theta _H := \pi \left ( 1- \left ( \frac {r_a+r_b}{2 r_b} \right )^{3/2} \right )$$ $$\triangleright$$ required phase angle before starting Hohmann transfer
6:    $$\omega _a := \sqrt {\frac {\mu }{r_a^3}}$$ $$\triangleright$$ angular speed of lower rad/sec
7:    $$\omega _b := \sqrt {\frac {\mu }{r_b^3}}$$ $$\triangleright$$ angular speed of higher satellite rad/sec
8:    if $$\theta _0\leq \theta _H$$ then $$\triangleright$$ adjust initial angle if needed
9:    $$\theta _0 := \theta _0+2 \pi$$
10:    end if
11:    $$\text {wait\_time} := \frac {\theta _0-\theta _H}{\omega _a-\omega _b}$$ $$\triangleright$$ how long to wait before starting Hohmann transfer
12:    $$\text {wait\_time} := \text {wait\_time} + \text {TOF}$$ $$\triangleright$$ now ready to go, add Hohmann transfer time
13:    return wait_time
14:   end function

An example implementation is below (in Maple)

hohmann_rendezvous_1:= proc({
theta::numeric:=0,
r1::numeric:=0,
r2::numeric:=0,
mu::numeric:=3.986*10^5})
local theta0,thetaH,TOF,a,omega1,omega2,wait_time;
theta0  := evalf(theta*Pi/180);
a       := (r1+r2)/2;
TOF     := Pi*(sqrt(a^3/mu));
omega1  := sqrt(mu/r1^3);
omega2  := sqrt(mu/r2^3);
thetaH  := evalf(Pi*(1-((r1+r2)/(2*r2))^(3/2)));
if theta0 <= thetaH then
theta0 := theta0+2*Pi;
fi;
wait_time := TOF+(theta0-thetaH)/(omega1-omega2);
eval(wait_time);
end proc:



And calling the above for two diﬀerent cases gives (times in hrs)

TOF:=hohmann_rendezvous_1(r1=6678,r2=6878,theta=0):
evalf(TOF/(60*60));
35.23480353



And

TOF:=hohmann_rendezvous_1(r1=6678,r2=6878,theta=280):
evalf(TOF/(60*60));
27.49212919


##### 3.10.1.3 Two satellite, separate orbits, rendezvous using bi-elliptic transfer, coplaner

In this transfer, the lower (fast satellite) does not have to wait for phase lock as in the case with Hohmann transfer. The transfer can starts immediately. There is a free parameter $$N$$ that one select depending on fuel cost requiments or any limitiation on the ﬁrst transfer orbit semi-major axes distance required. One can start with $$N=0$$ and adjust as needed.

1:   function hohmann_rendezvous_2($$\theta _0,r_1,r_2,N,\mu$$)
2:    $$\theta _0 := \theta _0 \frac {\pi }{180}$$ $$\triangleright$$ convert from degrees to radian
3:    $$\theta _H := \pi \left ( 1- \left ( \frac {r_1+r_2}{2 r_2} \right )^{3/2} \right )$$ $$\triangleright$$ Find Hohmann ideal phase angle before transfer
4:    if $$\theta _0 = \theta _H\, \textbf {and}\, N=0$$ then $$\triangleright$$ adjust for special case
5:    $$a := \frac {r_2+r_1}{2}$$
6:    $$\text {TOF} := \pi \left ( \sqrt {\frac {a^3}{\mu } } \right )$$
7:    else
8:    $$\omega _2 := \sqrt {\frac {\mu }{r_2^3}}$$ $$\triangleright$$ angular speed of slower satellite in rad/sec
9:    $$t_2 :=\frac {(2 \pi -\theta _0)+ 2\pi N}{\omega _2}$$ $$\triangleright$$ ﬁnd time of light of the slower satellite
10:    $$a_1 := \frac {r_t+r_1}{2}$$
11:    $$a_2 := \frac {r_t+r_2}{2}$$
12:    $$\text {TOF} := \pi \left ( \sqrt {\frac {a_1^3}{\mu } + \frac {a_2^3}{\mu }} \right )$$ $$\triangleright$$ time of ﬂight for the fast satellite
13:    $$r_t := \textbf {solve}\,\left (t_2 = TOF \right ) \text {for}\, r_t$$ $$\triangleright$$ Solve numerically for $$r_t$$
14:    end if
16:   end function

An example implementation is below in Maple

hohmann_rendezvous_2:= proc({
theta::numeric:=0,
r1::numeric:=0,
r2::numeric:=0,
N::nonnegint:=0,
mu::numeric:=3.986*10^5})
local theta0,thetaH,TOF;
theta0  := theta*Pi/180;
thetaH  := Pi*(1-((r1+r2)/(2*r2))^(3/2));
if theta0 = thetaH and N = 0 then
proc()
local a:=(r1+r2)/2;
TOF:= Pi*(sqrt(a^3/mu));
end proc()
else
proc()
local t2,a1,a2,rt,omega2;
omega2 := sqrt(mu/r2^3);
t2     := ((2*Pi-theta0)+2*Pi*N)/omega2;
a1     := (rt+r1)/2;
a2     := (rt+r2)/2;
TOF    := Pi*(sqrt(a1^3/mu)+sqrt(a2^3/mu));
rt     := op(select(is, [solve(t2=TOF,rt)], real));
end proc()
fi;
eval(TOF);
end proc:



And calling the above for two diﬀerent cases gives

TOF:=hohmann_rendezvous_2(theta=0,r1=6678,r2=6878,N=0):
evalf(TOF/(60*60)); #in hrs
1.576892101
TOF:=hohmann_rendezvous_2(theta=160,r1=6678,r2=6878,N=1):
evalf(TOF/(60*60));  #in hrs
2.452943266



#### 3.10.5 References

1. Orbital mechanics for Engineering students, second edition, by Howard D. Curtis
2. Orbital Mechanics, Vladimir Chobotov, second edition, AIAA
3. Fundamentals of Astrodynamics, Bates, Muller and White. Dover1971

1Image is from my class notes, page 7-9, EMA 550, Univ. Of Wisconsin, Madison by professor S. Sandrik