# Lab 5

## Get the code

Find the Lab 5 code in the class repository.

 ls
$cd labs/lab5 README cpp/ julia/ python/ rust/$ ls labs/lab5/cpp
makefile   orbit.cpp  root.cpp   view1.gp   view3.gp


The instructions here focus on C++, but you may carry out the comparable steps in Julia, Python, or Rust.

## Elliptic orbit

Consider a satellite orbiting the earth with period $$\mathsf{P}$$ whose distances of closest and farthest approach are $$\mathsf{r_1}$$ and $$\mathsf{r_2}$$. The satellite’s path traces out an ellipse,

$\mathsf{r}(\theta) = \mathsf{\frac{a(1-e^2)}{1+ e\cos(\theta)}} = \mathsf{a}(\mathsf{1}-\mathsf{e}\,\textsf{cos}\,\mathsf{E} ).$

The formula above assumes polar coordinates $$(\mathsf{r},\theta)$$ in the lane of the orbit with the earth centred on one of the ellipse’s foci. The ratio

$\mathsf{e} = \mathsf{\frac{r_2-r_1}{r_1+r_2}} = \mathsf{\frac{r_2 - r_1}{2a}}$

defines the eccentricity of the orbit; $$\mathsf{a}$$ is the semi-major axis, and $$\mathsf{E}$$ is the so-called eccentric anomaly.

Kepler worked out the following procedure for determining the location of the satellite at time $$\mathsf{t}$$, as measured from the moment of closest approach:

• Define the mean anomaly $$\mathsf{M} = \mathsf{2\pi t} / \mathsf{P}$$.

• Determine the eccentric anomaly $$\mathsf{E}$$ by solving Kepler’s equation, $$\mathsf{M} = \mathsf{E} - \mathsf{e}\,\textsf{sin}\,\mathsf{E}$$.

• Compute the true anomaly $$\theta$$ from the equation

$\textsf{tan}\,\mathsf{\frac{\theta}{2}} = \sqrt{\mathsf{\frac{1+e}{1-e}}}\,\textsf{tan}\,\mathsf{\frac{E}{2}}.$

Kepler’s equation is transcendental and has no closed-form solution. It has to be inverted numerically. In the file orbit.cpp, solve Kepler’s equation by finding the root of the function $$\mathsf{f(E)} = \mathsf{E}-\mathsf{M} - \mathsf{e\,\textsf{sin}\,E}$$ via Newton-Raphson iteration:

$\mathsf{E_{n+1}} := \mathsf{E_n} - \mathsf{\frac{f(E_n)}{f'(E_n)}} = \mathsf{E_n} - \frac{\mathsf{E_n} - \mathsf{M} - \mathsf{e}\, \textsf{sin}\,\mathsf{E_n}}{\mathsf{1}-\mathsf{e}\,\textsf{cos}\, \mathsf{E_n}}.$

Use the initial guess $$\mathsf{E_0} = \mathsf{0}$$ at time $$\mathsf{t} = \mathsf{0}$$, and for each subsequent time step use the previous step’s converged $$\mathsf{E}$$ value as the initial guess. Make sure that a minimum number of iterations are always performed. The view1.gp script should give you the following plots.

## More eccentric still

Write a new program orbit_vel.cpp that functions identically to orbit.cpp except that it outputs five columns of data: time $$\mathsf{t}$$, radius $$\mathsf{r}$$, angle $$\theta$$, radial velocity $$\dot{\mathsf{r}}$$, and angular velocity $$\dot{\theta}$$. The time derivatives $$\dot{\mathsf{r}}$$ and $$\dot{\theta}$$ should be approximated as symmetric finite difference. Some care must be taken in computing $$\dot{\theta}$$ since $$\theta$$ is compact on $$[\mathsf{0,2\pi}]$$.

In orbit.cpp the closest- and farthest-approach values were set to $$\mathsf{r_1} = \mathsf{R_{\textsf{e}}} + \mathsf{500}\,\textsf{km}$$ and $$\mathsf{r_2} = \mathsf{R_{\textsf{e}}} + \mathsf{3000}\,\textsf{km}$$, where $$\mathsf{R_\textsf{e}}$$ is the radius of the earth. For orbit_vel.cpp consider a more eccentric orbit with $$\mathsf{r_1} = \mathsf{R_{\textsf{e}}} + \mathsf{1000}\,\textsf{km}$$ and $$\mathsf{r_2} = \mathsf{R_{\textsf{e}}} + \mathsf{8500}\,\textsf{km}$$. Compose a gnuplot script view2.gp that plots the satellite’s spatial trajectory, its radial velocity versus time, its angular velocity versus time, and its total speed versus $$\textsf{sin}\,\theta$$, assuming that the data is in a file ov.dat. (Arrange these as four plots in succession, separated by a pause -1 command.) In the last plot, be sure to compute the speed as the magnitude of the vector $$\vec{\mathsf{v}} = \dot{\mathsf{r}}\,\hat{\mathsf{r}} + \mathsf{r}\dot{\mathsf{\theta}} \,\hat{\mathsf{\theta}}$$.

$make orbit_vel$ ./orbit_vel > ov.dat
$gnuplot -persist view3.gp  ## Square root by series expansion Recall that Newton’s method is an iterative scheme for finding the zeros of an arbitrary function $$\mathsf{f(x)}$$. It involves making an initial guess $$\mathsf{x_0}$$ and then generating a sequence of improved estimates according to $$\mathsf{x_{n+1}} := \mathsf{x_n} - \mathsf{f(x_n)/f'(x_n)}$$. If we choose $$\mathsf{f(x)} = \mathsf{x^2 - a}$$, then finding the zeros of $$\mathsf{f(x)}$$ is equivalent to computing the square root of $$\mathsf{a}$$. The correct recurrence relation is $\mathsf{x_{n+1}} := \mathsf{\frac{1}{2}}\biggl( \mathsf{x_n}+\mathsf{\frac{a}{x_n}}\biggr).$ The program root.cpp implements the recurrence relation shown above, starting from x = 1. The loop terminates when the next value in the sequence is sufficiently close to the old one. $ make root
g++ -o root root.cpp -O2 -ansi -pedantic -Wall -lm
$./root Returns the square root of the provided argument: Usage: root # [--verbose] ./root 2 Newton's method value: 1.414213562373095 C Math library value: 1.414213562373095$ ./root 9
Newton's method value: 3
C Math library value: 3
./root 101010 Newton's method value: 317.8207041713928 C Math library value: 317.8207041713928  In general, computing square roots via series expansion is much less reliable, but let’s give it a try. The square root $$\sqrt{\mathsf{a^2+b}}$$ can be expanded in powers of $$\mathsf{b/4a^2}$$ as follows: \begin{aligned} \sqrt{\mathsf{a^2 + b}} &= \mathsf{a} + \mathsf{\frac{1}{2}\frac{b}{a}} - \mathsf{\frac{1}{8}\frac{b^2}{a^3}} + \mathsf{\frac{1}{16}\frac{b^3}{a^5}} - \mathsf{\frac{5}{128}\frac{b^4}{a^7}} + \cdots\\ &= \mathsf{a} + \mathsf{\frac{b}{2a}} + \sum_{\mathsf{n=1}}^\infty \mathsf{(-1)^{n}} \mathcal{C}_n \mathsf{\frac{b^{n+1}}{(2a)^{2n+1}}}\\ &=\mathsf{a} + \mathsf{\frac{b}{2a}}\biggl(\mathsf{1} + \sum_{\mathsf{n=1}}^\infty \mathsf{(-1)^{n}} \mathcal{C}_n \mathsf{\frac{b^{n}}{(2a)^{2n}}}\biggr). \end{aligned} Here, $$(\mathcal{C}_n) = (\mathsf{1}, \mathsf{2}, \mathsf{5}, \mathsf{14}, \mathsf{42}, \mathsf{132}, \mathsf{429}, \mathsf{1430}, \ldots)$$ are the Catalan numbers. They are defined by $\mathcal{C}_\mathsf{n} = \mathsf{\frac{1}{n+1} {2n \choose n}} = \mathsf{\frac{(2n)!}{n!(n+1)!}}$ and describe the number of ways a polygon with $$\mathsf{n+2}$$ sides can be cut into $$\mathsf{n}$$ triangles. For large values of $$\mathsf{n}$$, the factorials are too large to compute, so we should use the trick of computing each term from the previous one. The ratio of two consecutive terms is $\frac{\mathsf{(-1)^{n+1}}\mathcal{C}_{\mathsf{n+1}}\mathsf{b^{n+1}}} {\mathsf{(2a)^{2n+2}}} \frac{\mathsf{(2a)^{2n}}}{\mathsf{(-1)^n} \mathcal{C}_\mathsf{n} \mathsf{b^n}} = -\mathsf{\frac{b(2n+2)(2n+1)}{4a^2(n+1)(n+2)}}.$ Write a program seriesroot.cpp that computes the truncated $$N$$-term series expansion for a given list of $$N$$ values. (In other words, argc can have any value greater than 3, and the program should loop over all N assigned from argv[3], argv[4], …, argv[argc-1].) You should be able to generate the one- through ten-term approximations to $$\sqrt{\mathsf{2}} = \sqrt{\mathsf{1^2+1}}$$ and $$\sqrt{\mathsf{3}} = \sqrt{\mathsf{2^2-1}}$$ as follows. The script view3.gp illustrates the convergence rates for Heron’s method and three different series approximations to $$\sqrt{\mathsf{3}}$$.  make seriesroot
g++ -o seriesroot seriesroot.cpp -O2 -ansi -pedantic -Wall -lm
$./seriesroot Computes the N-term series expansion of sqrt(a^2+b): Usage: seriesroot a b N1 [N2 N3 N4 ...]$ ./seriesroot 1 1 $(seq 10) 1 1 2 1.5 3 1.375 4 1.4375 5 1.3984375 6 1.42578125 7 1.4052734375 8 1.42138671875 9 1.408294677734375 10 1.419204711914062$ ./seriesroot 2 -1 $(seq 10) 1 2 2 1.75 3 1.734375 4 1.732421875 5 1.73211669921875 6 1.732063293457031 7 1.732053279876709 8 1.732051312923431 9 1.732050913386047 10 1.732050830149092$ gnuplot -persist view3.gp