W19 Assignment 3

Updated 2019-09-29

1. Orbital Changes 🪐

A small particle of mass $m$ is on a circular orbit of radius $R$ around a much larger mass $M$ . Suppose that the speed of mass $m$ is suddenly increased (in the same direction) by a factor $α > 1$ , that is, ${\vec{v}}_{final} = α {\vec{v}}_{initial}$ .

(a) Find expressions for the semi-major axis, the eccentricity, the semi-minor axis and the pericenter and apocentre distances, all in terms of just $R$ and $α$ . You may follow the logic used in the examples done in class.

We first start with the vis-visa equation from astrodynamics:

v^{2} = G M (\frac{2}{r} - \frac{1}{a})

By rearranging the equation, we can obtain the semi-major axis of the new orbit after the velocity change.

\frac{1}{a} = \frac{2}{R} - \frac{(α v_{i})^{2}}{G M}, v_{i} = \sqrt{\frac{G M}{R}}

Simplify:

\begin{aligned} \frac{1}{a} & = \frac{2}{R} - \frac{{(α \sqrt{\frac{G M}{R}})}^{2}}{G M} \\ = \frac{2}{R} - \frac{α^{2}}{R} \\ a & = \frac{R}{2 - α^{2}} \end{aligned}

Because $α > 1$ , then we the small particle mass $m$ is now at the periapsis. We can use the equation to get periapsis orbital altitude to find the new eccentricity of the orbit:

q = a (1 - e) = R e = 1 - \frac{R}{a} = 1 - \frac{R}{\frac{R}{2 - α^{2}}} = α^{2} - 1

The semi-minor axis is derived by Pythagorean theorem, since we already did that for assignment 1, I will jump into using the equation:

\begin{aligned} b^{2} & = a^{2} (1 - e^{2}) \\ b & = (\frac{R}{2 - α^{2}}) \sqrt{1 - (α^{2} - 1)^{2}} \\ = \frac{R}{2 - α^{2}} (α \sqrt{2 - α^{2}}) \\ = \frac{R α}{\sqrt{2 - α^{2}}} \end{aligned}

Since we increased the speed in the prograde direction, we increased our apoapsis. The new apoapsis is:

\begin{aligned} Q & = a (1 + e) \\ = \frac{R}{2 - α^{2}} (1 + (α^{2} - 1)) \\ = \frac{R α^{2}}{2 - α^{2}} \end{aligned}

(b) There is a maximum value of $α$ for which this is a sensible problem. Why is that, and what is the maximum value of $α$ ?

Because the geometry of the elliptical orbit is a conic section – a section cut of a cone shape. Given enough $α$ . the “tilt of the slice” will be a parabola and then a hyperbola (see figure below – source: voer.edu.vn). In physical sense, if we increase $α$ , our particle velocity would be so great that it would reach the escape velocity, and it will escape the orbit.

Precalculus: Rotation of Axes - VOER

The escape velocity, the factor we multiply our orbiting speed by $α$ is given as:

\begin{aligned} v_{e s c} = v_{f} & = \sqrt{\frac{2 G M}{a}} \\ α v_{i} & = \sqrt{\frac{2 G M}{R}}, v_{i} = \sqrt{\frac{G M}{R}} \\ α \sqrt{\frac{G M}{R}} & = \sqrt{\frac{2 G M}{R}} \\ α & = \sqrt{2} \end{aligned}

2. Kepler’s Equation

Kepler solved his own equation, $E - e \sin E = ϕ$ (where $ϕ = n t^{⋆}$ is the mean anomaly), by an iterative process; here you will explore that method. Them method relies on rewriting the equation as $E = ϕ + e \sin E$ . For an initial guess $E_{0} = ϕ$ , we can then proceed to a next guess for $E$ : $E_{1} = ϕ + e \sin E_{0}$ and then iterate until $E_{n} = ϕ + e \sin E_{n - 1} ≃ E_{n - 1}$ to within some tolerance.

Consider an ellipse with eccentricity 0.5. Use mean-anomaly points corresponding to (0, π/10, π/5, 3π/10, 2π/5, π/2, 3π/5, 7π/10, 4π/5, 9π/10, π).

(a) For ONE of these mean-anomaly points (which is not at 0 or π radians), produce a table showing each step in your iteration of Kepler’s equation until you reach a satisfactory value of $E$ to within a difference $| E_{n} - E_{n - 1} | < 0.01$ rad. Also show the calculation to compute the true anomaly.

I will choose 3π/10 (0.9425 rad) as the mean-anomaly point. I wrote a C program to help me compute the iterations. See the attached Appendix A for the C program code. Here are the results:

Iterations	Eccentric anomaly
0	0.9425
1	1.3470
2	1.4300
3	1.4375

The true anomaly ( $θ$ ) can be found using:

\begin{aligned} \tan (\frac{θ}{2}) & = \sqrt{\frac{1 + e}{1 - e}} \tan (\frac{E}{2}) \\ θ & = 2 \cdot \tan^{- 1} (\sqrt{\frac{1 + e}{1 - e}} \tan (\frac{E}{2})) \end{aligned}

(b) Make a second table with the full set of mean-anomaly points, the corresponding eccentric-anomaly $E$ values and the corresponding true-anomaly $θ$ values.

Again, I wrote a C program to compute the all the values for the table. See Appendix A for the C program code.

Mean-anomaly $ϕ$	Eccentric-anomaly $E$	True-anomaly $θ$
0.0 $π$	0.0000	0.0000
0.1 $π$	0.5900	0.9690
0.2 $π$	1.0635	1.5896
0.3 $π$	1.4375	1.9750
0.4 $π$	1.7486	2.2419
0.5 $π$	2.0204	2.4462
0.6 $π$	2.2692	2.6158
0.7 $π$	2.5016	2.7635
0.8 $π$	2.7168	2.8938
0.9 $π$	2.9345	3.0218
1.0 $π$	3.1416	3.1416

(c) Make a plot showing $E$ and $θ$ as functions of $ϕ$ and comment on whether it would be OK to neglect the orbit’s eccentricity when estimating the position of the object.

Using the data from Part (B), I used Excel to create the following plot.

1569798435590

As the graph clearly shows, there is a significant discrepancy between the eccentric anomaly, and the true anomaly of an elliptical orbit. Therefore it is not OK to neglect the eccentricity, in this case $e = 0.5$ , when estimating the position of the object.

3. Temperature Changes in an Eccentric Orbit 🌗

The planet HD 20782b orbits HD 20782, a G3V star with mass 1.05 times that of the Sun and luminosity 1.25 times that of the Sun. The planet’s orbit has a semi-major axis of 1.397 au, an orbital period of 597 days, and an eccentricity of 0.956. The planet’s mass is about 2 times that of Jupiter and you may assume its albedo is similar to that of Jupiter at about 0.5.

(a) (Review) Compute the periapsis and apoapsis distances of the planet from the star.

q = a (1 - e) = (1.397) (1 - 0.956) = 0.0615 au Q = a (1 + e) = (1.397) (1 + 0.956) = 2.7325 au

(b) Assume the planet is spinning rapidly. In class we derived the dependence of planetary blackbody temperature as a function of orbital distance in the Solar system: $T_{p} = 279 K (1 - A)^{1 / 4} \frac{1}{\sqrt{r_{au}}}$ . Adapt this formula to this new stellar system — that is, work out the new constant that should take the place of 279 K, and take into account the numerical value of the albedo.

First, we need to adapt the luminosity. The sun’s luminosity is the total power outputted by the sun: $L_\odot=4\pi R^2\odot \sigma T\odot^4=3.839\times 10^{26}$W. We multiply this by 1.25 and get:

L_{s} = 4.799 \times 10^{26} W

To find the temperature, we need to equate the power intercepted from the star HD 20782 and the power emitted (assuming that temperature on the planet is much lower than the star). Power intercepted is a function of distance from the star ( $r$ ) and is given by flux × surface area × absorption:

P_{i n} (r) = F_{s} \times A_{p_{intercept}} \times (1 - A)

Where $F_{s}$ is the flux [W/m²] which is the luminosity divided by surface area: $F_{s} (r) = L_{s} / 4 π r^{2}$ . The area intercepted by the planet is the “shadow” of the planet which is a disk, given by $A_{p_{intercept}} = π R_{p}^{2}$ where $R_{p}$ is the radius of the planet.

The power outputted by the planet is by:

P_{o u t} = A_{p} \times σ \times T_{p}^{4}

Where $A_{p}$ is the planet’s surface area, and $T_{p}$ is the surface temperature of the planet.

We equate power in, and power out and simplify:

\begin{aligned} P_{i n} & = P_{o u t} \\ F_{s} \cdot A_{p_{intercept}} \cdot (1 - A) & = A_{p} \cdot σ \cdot T_{p}^{4} \\ \frac{L_{s}}{4 π r^{2}} \cdot π R_{p}^{2} \cdot (1 - A) & = 4 π R_{p}^{2} \cdot σ \cdot T_{p}^{4} \\ \frac{L_{s}}{4 r^{2}} (1 - A) & = 4 π σ T_{p}^{4} \end{aligned}

Now we rearrange the equation to match the equation in question by isolating for $T_{p}$ :

\begin{aligned} \frac{L_{s}}{4 r^{2}} (1 - A) & = 4 π σ T_{p}^{4} \\ T_{p}^{4} & = \frac{L_{s}}{16 π r^{2} σ} (1 - A) \\ T_{p} & = {(\frac{L_{s}}{16 π σ})}^{1 / 4} \frac{1}{\sqrt{r}} (1 - A) \end{aligned}

We need to do one more thing. The equation we had at the top has the distance from the star $r$ in units of AU. So we need to account for it here. 149597900000m = 1 AU.

\begin{aligned} T_{p} & = {(\frac{L_{s}}{16 π σ})}^{1 / 4} (\frac{1}{\sqrt{r_{m}}} \times \sqrt{\frac{1 au}{1.496 \times 10^{11} m}}) (1 - A) \\ = [{(\frac{L_{s}}{16 π σ})}^{1 / 4} \frac{1 \sqrt{au}}{\sqrt{1.496 \times 10^{11} m}}] (1 - A) \frac{1}{\sqrt{r_{a u}}} \end{aligned}

Verifying the units: in this section, I’m only showing my work to see that the units for the above derivation and conversion makes sense.
$\begin{aligned} K & = [{(\frac{W}{{Wm}^{- 2} K^{- 4}})}^{1 / 4} \frac{{au}^{1 / 2}}{m^{1 / 2}}] \frac{1}{{au}^{1 / 2}} \\ = [{(m^{2} K^{4})}^{1 / 4} \frac{{au}^{1 / 2}}{m^{1 / 2}}] \frac{1}{{au}^{1 / 2}} \\ = m^{1 / 2} K \frac{{au}^{1 / 2}}{m^{1 / 2}} \frac{1}{{au}^{1 / 2}} \\ = K {au}^{1 / 2} \frac{1}{{au}^{1 / 2}} \\ = K \end{aligned}$
👌 Cool. We can proceed.

Plugging in all the values we calculated so far, where $σ$ is Stefan-Boltzmann consant. The constant that replaces the 279K is 294.5K:

T_{p} = [{(\frac{L_{s}}{16 π σ})}^{1 / 4} \frac{1 \sqrt{au}}{\sqrt{1.496 \times 10^{11} m}}] (1 - A) \frac{1}{\sqrt{r_{a u}}} T_{p} = 294.5 K (1 - A) \frac{1}{\sqrt{r_{a u}}}

(c) Assuming that there no signicant internal heating from the planet, determine its temperature and the wavelength of peak emission at both periapsis and apoapsis. In what part of the spectrum do these wavelengths fall?

At Periapsis: the planet’s distance from the star is only 0.0615au. The temperature is 594K. 🔥 $T_{p_{p e}} = 294.5 K (1 - (0.5)) \frac{1}{\sqrt{0.0615}} = 594 K$

The wavelength of maximum intensity or peak wavelength $λ_{p}$ is in accordance with Wien’s law. We plug in the numbers here: $λ_{p_{p e}} = \frac{2.9 \times 10^{3} μ m}{T} = \frac{2.9 \times 10^{3} μ m}{594} = 4.884 μ m$

At Apoapsis: the planet’s distance from the star is 2.7325au. The temperature is 89.1K. ⛄ $T_{p_{a p}} = 294.5 K (1 - (0.5)) \frac{1}{\sqrt{2.7325}} = 89.1 K$

Using similar calculations as above, we get peak wavelength at:

λ_{p_{a p}} = \frac{2.9 \times 10^{3} μ m}{89.1} = 32.6 μ m

So when the planet is at periapsis, the planet is very hot at 594K and emits peak wavelength of 4.9μm in the infrared spectrum. When the planet is at apoapsis, the planet is cold at only 89.1K and emits a peak wavelength of 32.6μm in the far-infrared spectrum.

Appendix A

Code Used for Problem 2 Part 1

#include <stdio.h>
#include <math.h>

#define THRES 0.01
#define PI    acos(-1)

void part1(void)
{
	printf("Part 1 #################\n");
	const double starting_mean_anomaly = 3.0 * PI / 10.0;
	const double eccentricity = 0.5;

	double eccentric_anomaly;
	double eccentric_anomaly_prev = starting_mean_anomaly;

	int iterations = 1;

	printf("| Iterations | Eccentric anomaly |\n");
	printf("|--:|--:|\n");
	printf("| 0 | %.4f |\n", starting_mean_anomaly);

	eccentric_anomaly = starting_mean_anomaly + eccentricity * sin(eccentric_anomaly_prev);
	printf("| 1 | %.4f |\n", eccentric_anomaly);

	while (fabs(eccentric_anomaly_prev - eccentric_anomaly) > THRES)
	{
		eccentric_anomaly_prev = eccentric_anomaly;
		eccentric_anomaly = starting_mean_anomaly + eccentricity * sin(eccentric_anomaly_prev);
		iterations++;
		printf("| %d | %.4f |\n", iterations, eccentric_anomaly);
	}

	printf("E_%d = %.4f rad\n", iterations, eccentric_anomaly);
}

int main(void)
{
	part1();
}

Code Used for Problem 2 Part 2

#include <stdio.h>
#include <math.h>

#define THRES 0.01
#define PI    acos(-1)

void part2(void)
{
	const double eccentricity = 0.5;
	const double mean_anomalies[] = {
        0.0, 0.1 * PI, 0.2 * PI, 0.3 * PI, 0.4 * PI,
        0.5 * PI, 0.6 * PI, 0.7 * PI, 0.8 * PI, 0.9 * PI, PI
	};
	const int mean_anomalies_length = 11;

	// Ouptut
	printf("Part 2 #################\n");
	printf("| Mean-anomaly | Eccentric-anomaly $E$ | True-anomaly $\\theta$ |\n");
	printf("|--:|--:|--:|\n");

	for (int i = 0; i < mean_anomalies_length; i++)
	{
		// Mean anomaly
		const double mult = mean_anomalies[i] / PI;

		// Get E
		double eccentric_anomaly;
		double eccentric_anomaly_prev = mean_anomalies[i];

		eccentric_anomaly = mean_anomalies[i] + eccentricity * sin(eccentric_anomaly_prev);
		while (fabs(eccentric_anomaly_prev - eccentric_anomaly) > THRES)
		{
			eccentric_anomaly_prev = eccentric_anomaly;
			eccentric_anomaly = mean_anomalies[i] + eccentricity * sin(eccentric_anomaly_prev);
		}
		
		// Get theta
		double true_anomaly = 2.0 * atan(
            sqrt((1 + eccentricity) / (1 - eccentricity)) * tan(eccentric_anomaly / 2)
         );

		// Print row
		printf("| %.1f$\\pi$ | %.4f | %.4f |\n", mult, eccentric_anomaly, true_anomaly);
	}
}

int main(void)
{
	part2();
}