Project 3 — Root Finding and Stellar Inference

Due Wednesday, October 7, 11:59 pm PT · pairs allowed, own code and own write-up

What you are building¶

How do you find an unknown that appears on both sides of an equation?

Build a root-finding library to recover Wien’s constant, compute orbital positions, and locate transit contacts. Then use algebra and a Star class to compare equilibrium-temperature zones for your fifteen Project 1 stars.

Answer each part’s interpretation questions in one or two short paragraphs after its calculations and figures in analysis.ipynb. The Deliverables callouts summarize the requirements; there is no separate report. They are completion reminders, not a points rubric. Your work is assessed using the course’s four project dimensions: correct model and validation, scientific reasoning, reproducibility and code design, and communication and figures.

Check each layer before combining them¶

A passing check supports that layer; it does not establish that every later calculation is correct. The required tests and comparisons below implement this sequence.

Structure: classes start here¶

project-3/
├── README.md
├── roots.py           # bisection, Newton, secant — plain functions
├── star.py            # the Star class
├── orbit.py           # Kepler's equation, z(t), contact points
├── transit.py         # your existing transit functions
├── integrate.py       # retained from Project 2
├── conftest.py        # retain the Project 2 test layout
├── tests/
└── analysis.ipynb

Use python -m pytest from the project directory. Reuse your tested Project 1–2 functions; you do not need to implement another quadrature method here.

Validate the star’s stored radius and temperature. Compute luminosity and peak wavelength from that state so the derived quantities stay consistent.

Part 1 — The library¶

Implement these functions in roots.py:

def bisect(f, a, b, atol=1e-12, rtol=0.0, max_iter=200): ...
def newton(f, fprime, x0, atol=1e-12, rtol=0.0, max_iter=200): ...
def secant(f, x0, x1, atol=1e-12, rtol=0.0, max_iter=200): ...

Each returns (root, n_iter, history):

• root: the final numerical estimate of the root.

• n_iter: the integer iteration count, defined for each method below.

• history: the ordered sequence of root estimates recorded by the method.

Keep the Lecture 8 bisection contract and use these recording conventions:

• Bisection: history contains evaluated midpoints; n_iter == len(history). Return the last evaluated midpoint. An endpoint root returns that endpoint, zero iterations, and an empty history.

• Newton: history begins with x0, followed by each updated iterate; n_iter == len(history) - 1.

• Secant: history begins with x0, x1, followed by each updated iterate; n_iter == len(history) - 2. Check both initial residuals before dividing. An initial root returns with zero updates and the two initial entries.

Newton and secant may return immediately at a computed zero residual; otherwise use the stopping rule below. Bisection counts midpoint evaluations; the open methods count updates. These are not total function-evaluation costs.

Use ValueError for invalid inputs, an invalid bracket, nonfinite function or derivative values, a zero Newton derivative at a nonroot, or a zero secant denominator when neither input is a root. Require finite real scalar inputs and function values, finite nonnegative tolerances with at least one positive, and a positive integer max_iter excluding Booleans. Bisection requires a < b; secant requires distinct initial inputs.

Use RuntimeError for iteration-budget exhaustion or a bisection interval whose tolerance is unmet but has no representable interior midpoint. Compare signs directly rather than multiplying function values. Continuity and reliable signs are mathematical preconditions; input validation cannot establish them.

What a tolerance means¶

The allowance at a candidate $x$ is

atol is the absolute tolerance, in the same units as the root. A positive value supplies an allowance near zero. rtol is the relative tolerance, a dimensionless fraction that scales the allowance with the candidate root’s magnitude. Near 10¹², an absolute tolerance of 10^-12 is below double-precision spacing.

For bisection, evaluate the midpoint $m$ of the current bracket and use $B_m=\max(m-a,b-m)$. Accept a computed zero residual or $B_m\le\tau_m$. For an exact midpoint, $B_m=(b-a)/2$. Use the bracket before updating it to bound that midpoint’s error, and record the midpoint before returning.

For Newton and secant, accept an updated iterate if its computed residual is zero or its last step is no larger than $\tau_x$. A small step does not certify a small residual or root error; check the defining equation in each application. A computed zero can also result from rounding or underflow.

Measure convergence¶

Use $f(x)=x-5(1-e^{-x})$, with bisection on $[3,10]$, Newton from $x_0=4$, and secant from $(x_0,x_1)=(3,6)$. Use $x_\star=4.965114231744276$ as the double-precision reference; do not round it before computing errors.

Plot $e_n=|x_n-x_\star|$ against history index on a log $y$ axis. State which entries are starting guesses. Report each solver’s returned iteration count and final error. Distinguish satisfying its stopping rule from the first recorded error below 10^-12. Omit zero errors from the logarithmic plot and identify them in the caption rather than replacing them with invented values.

• Newton and secant: estimate the order using

  $$p_n=\frac{\ln(e_{n+1}/e_n)}{\ln(e_n/e_{n-1})}.$$

  (2)

  Use and report the last consecutive triple with $10^{-13}<e_n<1$ for all three errors and a nonzero denominator. This band avoids the numerical floor here but does not establish an asymptotic regime. If no triple qualifies, report that. Compare with orders 2 and $(1+\sqrt{5})/2$.

• Bisection: reconstruct the first 20 brackets from the initial bracket and recorded midpoints without changing the return interface. Measure successive width ratios and the slope of $\ln(b_n-a_n)$ against $n$. Compare with $1/2$ and $-\ln 2$. Plot the bracket half-width alongside the actual errors. Actual midpoint errors need not decrease at every update.

Part 2 — Where Wien’s constant comes from¶

The wavelength form of Wien’s law is $\lambda_{\max}T=b$, with tabulated $b\approx0.2897772$ cm·K. Recover this value from the Planck spectrum:

Use the same CGS constants as Project 2:

Investigate root selection¶

Run your Newton method from $x_0=1$ and inspect its history.

1. Report the returned root and residual. Compute the implied $b$ and $\lambda_{\max}$ at 5772 K if the root is nonzero. If it is exactly zero, state why these quantities are undefined instead of dividing by zero.

2. Explain why $f(0)=0$ does not identify a finite-wavelength spectral peak.

3. Recover the desired root using a bracket or starting guess justified by the function and its physical domain.

4. Plot $f$ and $f'$ in two panels, marking both roots, the Newton history, and the derivative’s zero at $x=\ln5$. Run Newton from that point and report the outcome. The evaluated derivative may be zero or merely tiny in floating point; an enormous step or reported failure is evidence to explain.

The near-zero result depends on how you evaluate $f$; no particular sign or last digit is expected. A small derivative can send Newton outside the domain.

Part 3 — The orbit¶

Read the Science Background for diagrams and derivations of the orbital angles, sky projection, and transit contacts.

Project 1 used $z(t)=(a/R_\star)|\sin(2\pi t/P)|$ for a circular, edge-on orbit. Now compute the projected separation for an eccentric orbit.

For a bound Kepler orbit, $0\le e<1$. With $t=0$ at periastron,

Mean anomaly $M$ increases uniformly with time. Eccentric anomaly $E$ is a geometric parameter on the ellipse. Given $M$ and $e$, solve for $E$.

What to build¶

In orbit.py:

def solve_kepler(M, e, atol=1e-12, rtol=0.0, max_iter=100): ...  # (E, n_iter)
def true_anomaly(E, e): ...                                 # radians
def separation(t, period, a_over_rstar, e, inc): ...          # scalar or array

Use your own Newton from Part 1, with closures:

resid  = lambda E: E - e*np.sin(E) - M
dresid = lambda E: 1 - e*np.cos(E)

Use the quadrant-aware expression

For $\omega=\pi/2$, a stable equivalent projected separation is

This sum of squares avoids cancellation near an edge-on mid-transit. Explain its equivalence to the geometric expression above.

Conventions: angles in radians; t and period in the same units (days for the examples); positive period and a_over_rstar; $0\le e<1$ and $0\le i\le\pi/2$. solve_kepler takes scalar $M,e$ and returns the solution on the same unwrapped branch as $M$. If you reduce $M$ modulo $2\pi$ internally, restore that branch in the returned $E$. separation returns a scalar for a scalar input and preserves an array’s shape. Reject nonfinite or out-of-domain inputs with ValueError.

Compare starting guesses and check the circular limit¶

Use M = np.linspace(0, 2*np.pi, 400) at $e=0,0.3,0.6,0.9,0.99$. Compare $E_0=M$ with $E_0=\pi$ on this interval, using the same tolerances and max_iter=100. You may call newton directly for this comparison; no extra argument in solve_kepler is required. Report failures, mean and maximum iterations among successful runs, and their residuals. Use this evidence to choose and document a solve_kepler starting strategy. Failures on this grid do not establish a universal eccentricity threshold. The high-$e$ cases are solver tests, not viable transit models at $a/R_\star=10$.

For the circular-limit check, use $P=3.0$, $a/R_\star=10$, t = np.linspace(-0.25, 0.25, 501), $e=0$, and $i=\pi/2$. Compare with Project 1’s formula using rtol=0, atol=1e-11. Report the maximum absolute error, not a prescribed round-off pattern. Relative error is undefined where the reference separation is zero.

Part 4 — Contact points and duration¶

Let $p=R_p/R_\star$ be the planet-to-star radius ratio. For a nongrazing transit with $0<p<1$, four contacts mark tangency:

Solve $z(t)-(1\pm p)=0$. Return ordered times $(t_1,t_2,t_3,t_4)$ from contact_points(period, a_over_rstar, p, e, inc), and report $T_{14}=t_4-t_1$ and $T_{23}=t_3-t_2$. $T_{23}$ is the full-overlap duration; it is a flat-bottom duration for the uniform-brightness model used here.

Use $P=3.0$ days, $a/R_\star=10$, $p=0.1$, $i=\pi/2$, and $e=0,0.3,0.6$, unless a check specifies otherwise. These required configurations have four distinct contacts. Report failure if four distinct contacts are not found; a general grazing/no-transit model is not required.

Bracket first, then solve. For these configurations, sample $[-P/4,P/4]$, locate sign changes for each contact equation, and pass each bracket to bisection. Retain exact sampled roots once. Refine the sampling if needed and check that the four times are distinct and correctly ordered.

Circular check and eccentric light curves¶

At $e=0$, $i=\pi/2$, compare against your Project 1 grid estimates for $p=0.05$ and 0.25, and against

Compute the analytic reference; the table is rounded. Report the root-found duration error and the original grid spacing. Errors in both contact times and the nested Kepler solve affect duration: a 10^-12 stopping tolerance on one time does not promise 10^-13 duration accuracy.

Plot the three required eccentric light curves using your uniform-disk flux(p, z) from transit.py. Mark all four contacts on each curve and report $T_{14}$ and $T_{23}$. The period, scaled semimajor axis, inclination, and periastron orientation remain fixed as eccentricity changes.

Part 5 — Stars, and where the algebra wins¶

Implement the complete Star class in star.py with this interface:

class Star:
    def __init__(self, name, radius, teff): ...  # radius in R_sun; teff in K
    @property
    def luminosity(self): ...                  # L_sun
    @property
    def lambda_max(self): ...                  # cm, using your derived b
    def habitable_zone(self, albedo=0.3): ...   # (inner, outer) in AU

Store readable name, radius, and teff attributes; radius is in solar radii and temperature is in kelvin. Require finite, positive radius and temperature at construction. If those values can change later, validate each new value; a read-only design is also acceptable. Compute luminosity and lambda_max from the current state rather than storing separate copies.

Reuse your tested Project 1 Stefan–Boltzmann calculation. You may copy its luminosity(radius, teff, solar=False) function into star.py, import it from a submitted module, or refactor it while preserving its units and behavior. The Star.luminosity property must return $L_\odot$ from the object’s current radius and temperature. Check that the Sun gives the same result as in Project 1.

Use $1\ \mathrm{AU}=1.495979\times10^{13}$ cm and $R_\odot=6.957\times10^{10}$ cm.

For a planet with Bond albedo $A$, uniform redistribution of absorbed energy, unit thermal emissivity, and no internal heating, radiative balance gives

Use the existing simplified criterion $273<T_{\rm eq}<373$ K at $A=0.3$. Require $0\le A<1$. The inner boundary corresponds to 373 K and the outer boundary to 273 K. Rearrange this equation for $d$ by hand and implement the algebraic result. No root finder is needed.

Instantiate all fifteen Project 1 stars from the local data. Plot one horizontal interval per star, ordered and labelled by $T_{\rm eff}$. Use a log distance axis in AU; mark Earth at 1 AU and identify the Sun and Proxima Centauri.

One last question¶

End the notebook with a markdown cell headed Looking back:

Which single number in this notebook do you trust least, and why?

Answer in two or three sentences. Any honest answer earns full credit.

Required tests¶

Include these eight test cases in tests/, retaining your Lecture 8 solver contract tests. A case may use several assertions or parameterized inputs.

Analytic-limit checks and defining-equation residuals catch different errors; keep both. Use explicit tolerances rather than relying on comparison defaults.

Required figures¶

1. Convergence: the three error histories and bisection’s bracket half-width; log error axis, with starting guesses and omitted zero errors identified.

2. Wien root selection: $f$ and $f'$ in two panels, both roots, the $x_0=1$ history, and the derivative zero. Axis limits must contain the displayed iterates.

3. Transit timing: uniform-disk light curves at $e=0,0.3,0.6$, with all four contacts marked on each.

4. Fifteen-star zones: labelled temperature ordering, log distance in AU, Earth at 1 AU, and the Sun and Proxima identified.

Use Project 1’s figure standards: labelled units, readable mathematical labels, legends clear of the data, and figures saved to files.

What you turn in¶

• The repository with modules and tests, passing under python -m pytest.

• analysis.ipynb with calculations, the five interpretations, and Looking back.

• The four figures saved as files.

• README.md explaining the contents and how to run the work.

Pairs are allowed, but each student writes their own code and answers and must be able to explain every submitted line.