Lab 5 ===== Get the code ------------ Use the ``curl`` command to download from the class website everything you'll need for the lab. .. code-block:: console $ curl https://www.phy.olemiss.edu/~kbeach/phys540/src/lab5.tgz -O $ tar xzf lab5.tgz $ cd lab5 cpp/ julia/ python/ rust/ $ ls lab5/cpp makefile orbit.cpp root.cpp view1.gp view3.gp $ cd lab5/cpp *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 :math:`\mathsf{P}` whose distances of closest and farthest approach are :math:`\mathsf{r_1}` and :math:`\mathsf{r_2}`. The satellite's path traces out an ellipse, .. math:: \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 :math:`(\mathsf{r},\theta)` in the lane of the orbit with the earth centred on one of the ellipse's foci. The ratio .. math:: \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; :math:`\mathsf{a}` is the semi-major axis, and :math:`\mathsf{E}` is the so-called *eccentric anomaly*. Kepler worked out the following procedure for determining the location of the satellite at time :math:`\mathsf{t}`, as measured from the moment of closest approach: * Define the mean anomaly :math:`\mathsf{M} = \mathsf{2\pi t} / \mathsf{P}`. * Determine the eccentric anomaly :math:`\mathsf{E}` by solving Kepler's equation, :math:`\mathsf{M} = \mathsf{E} - \mathsf{e}\,\textsf{sin}\,\mathsf{E}`. * Compute the true anomaly :math:`\theta` from the equation .. math:: \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 :math:`\mathsf{f(E)} = \mathsf{E}-\mathsf{M} - \mathsf{e\,\textsf{sin}\,E}` via Newton-Raphson iteration: .. math:: \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 :math:`\mathsf{E_0} = \mathsf{0}` at time :math:`\mathsf{t} = \mathsf{0}`, and for each subsequent time step use the previous step's converged :math:`\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. .. figure:: figures/orbit1.* .. figure:: figures/orbit2.* .. figure:: figures/orbit3.* 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 :math:`\mathsf{t}`, radius :math:`\mathsf{r}`, angle :math:`\theta`, radial velocity :math:`\dot{\mathsf{r}}`, and angular velocity :math:`\dot{\theta}`. The time derivatives :math:`\dot{\mathsf{r}}` and :math:`\dot{\theta}` should be approximated as symmetric finite difference. Some care must be taken in computing :math:`\dot{\theta}` since :math:`\theta` is compact on :math:`[\mathsf{0,2\pi}]`. In ``orbit.cpp`` the closest- and farthest-approach values were set to :math:`\mathsf{r_1} = \mathsf{R_{\textsf{e}}} + \mathsf{500}\,\textsf{km}` and :math:`\mathsf{r_2} = \mathsf{R_{\textsf{e}}} + \mathsf{3000}\,\textsf{km}`, where :math:`\mathsf{R_\textsf{e}}` is the radius of the earth. For ``orbit_vel.cpp`` consider a more eccentric orbit with :math:`\mathsf{r_1} = \mathsf{R_{\textsf{e}}} + \mathsf{1000}\,\textsf{km}` and :math:`\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 :math:`\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 :math:`\vec{\mathsf{v}} = \dot{\mathsf{r}}\,\hat{\mathsf{r}} + \mathsf{r}\dot{\mathsf{\theta}} \,\hat{\mathsf{\theta}}`. .. code-block:: console $ make orbit_vel $ ./orbit_vel > ov.dat $ gnuplot -persist view3.gp .. figure:: figures/orbit_vel1.* .. figure:: figures/orbit_vel2.* .. figure:: figures/orbit_vel3.* .. figure:: figures/orbit_vel4.* Square root by series expansion ------------------------------- Recall that Newton's method is an iterative scheme for finding the zeros of an arbitrary function :math:`\mathsf{f(x)}`. It involves making an initial guess :math:`\mathsf{x_0}` and then generating a sequence of improved estimates according to :math:`\mathsf{x_{n+1}} := \mathsf{x_n} - \mathsf{f(x_n)/f'(x_n)}`. If we choose :math:`\mathsf{f(x)} = \mathsf{x^2 - a}`, then finding the zeros of :math:`\mathsf{f(x)}` is equivalent to computing the square root of :math:`\mathsf{a}`. The correct recurrence relation is .. math:: \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. .. code-block:: console $ 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 :math:`\sqrt{\mathsf{a^2+b}}` can be expanded in powers of :math:`\mathsf{b/4a^2}` as follows: .. math:: \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, :math:`(\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 .. math:: \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 :math:`\mathsf{n+2}` sides can be cut into :math:`\mathsf{n}` triangles. For large values of :math:`\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 .. math:: \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 :math:`N`-term series expansion for a given list of :math:`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 :math:`\sqrt{\mathsf{2}} = \sqrt{\mathsf{1^2+1}}` and :math:`\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 :math:`\sqrt{\mathsf{3}}`. .. code-block:: console $ 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