Phugoid Motion - Practical Numerical Methods

In the first module of Practical Numerical Methods we study problems modeled by ordinary differential equations.

Our motivating problem is the phugoid model of glider flight. In this first lesson, we discuss some background, explain the physics, and work out the mathematical model.

We first consider the idealized model with no drag, resulting in a simple harmonic motion. We will plot some interesting trajectories of a glider under phugoid oscillations under the ideal model. In the next lesson, you’ll learn to numerically integrate the differential equation using Euler’s method.

The lesson pattern

This is our first complete pass through a practice that will repeat throughout the course:

Derive — On paper: establish an expectation independently of the computer.
Reconstruct — In your notebook: rebuild the central calculation from the worked lesson.
Specify — Before using an agent: state what is being delegated and what evidence would make its result acceptable.
Audit — Review and verify: inspect the proposal and test whether the evidence can expose a plausible defect.
Explain — Your verdict: decide what the evidence supports and what remains uncertain.

The specification step may be unfamiliar. It is not merely the wording of a prompt: it records the outcome, model context, interface, constraints, permitted access, and acceptance evidence that you have decided before delegating work. This lesson provides a detailed scaffold; later lessons will ask you to make more of these decisions yourself.

The term phugoid in aeronautics refers to a motion pattern where an aircraft oscillates up and down —nose-up and climb, then nose-down and descend— around an equilibrium trajectory. The aircraft oscillates in altitude, speed and pitch, with only small (neglected) variations in the angle of attack, as it repeatedly exchanges kinetic and potential energy Sinha & Ananthkrishnan, 2013, chap. 5.

A low-amplitude phugoid motion can be just a nuisance, as the aircraft does not exceed the stall angle of attack and nothing bad happens. But the mode can also be unstable leading to a stall or even a loop.

Look at this video of a simulator showing a Cessna single-engine airplane in phugoid motion:

from IPython.display import YouTubeVideo
YouTubeVideo('ysdU4mnRYdM')

That doesn’t look too good! What’s happening?

It can get a lot worse when an aircraft enters one of these modes that is unstable. A famous example is that of NASA’s Helios Solar Powered Aircraft prototype, which broke up in mid air due to extreme phugoid oscillations.

Helios was a proof-of-concept solar electric-powered flying wing that broke the world altitude record for a non-rocket-powered aircraft in August 2001. But in June 26, 2003, it broke something else. The aircraft entered phugoid motion after encountering turbulence near the Hawaiian Island of Kauai. The high speed reached in the oscillations exceeded the design limits, and it ended up wrecked in the Pacific Ocean. Luckily, the Helios was remotely operated, and nobody got hurt.

The physics of phugoids¶

In a phugoid oscillation the aircraft pitches up and down, as it decelerates and accelerates. The trajectory might look like a sinusoid, as shown in Figure 1. The assumption is that the forward velocity of the aircraft, $v$ , varies in such a way that the angle of attack remains (nearly) constant, which means that we can assume a constant lift coefficient.

A sinusoidal trajectory showing an aircraft oscillating above and below a central flight path — Figure 1:Trajectory of an aircraft in phugoid motion.

In the descending portion of the trajectory, the aircraft’s velocity increases as it proceeds from a peak to the minimum height—gaining kinetic energy at the expense of potential energy. The contrary happens in the upward segment, as its velocity decreases there.

We measure the pitch angle (between the aircraft’s longitudinal axis and the horizontal) as positive when the aircraft’s nose is pointing up. In the portion of the trajectory below the center-line, where it curves upwards, the pitch angle $\theta$ is increasing: $\dot{\theta}>0$ . And where the trajectory curves down, the pitch angle is decreasing: $\dot{\theta}<0$ , as shown in Figure 1.

Let’s review the forces affecting an aircraft in a downward glide. In Figure 2, we show the flight path, the forces on the glider (no thrust), and the glide angle or flight path angle, $\gamma$ , between the flight path and the horizontal.

Lift, drag, and weight acting on a glider along a descending flight path — Figure 2:Forces on a glider.

The force of lift, $L$ —created by the airflow around the wings— is perpendicular to the trajectory, and the force of drag, $D$ , is parallel to the trajectory. Both forces are expressed in terms of coefficients of lift and drag, $C_L$ and $C_D$ , respectively, that depend on the wing design and angle of attack: the angle between the wing chord and the flight path.

If you are not familiar with airplane aerodynamics, you might be getting confused with some terms here ... and all those angles! But be patient and look things up, if you need to. We’re giving you a quick summary here.

Lift and drag are proportional to a surface area, $S$ , and to the dynamic pressure: $1/2 \rho v^2$ , where $\rho$ is the density of air, and $v$ the forward velocity of the aircraft. The equations for lift and drag are:

L = C_L S \times \frac{1}{2} \rho v^2

(1)

D = C_D S \times \frac{1}{2} \rho v^2

(2)

When the glider is in equilibrium, the forces balance each other. We can equate the forces in the directions perpendicular and parallel to the trajectory, as follows:

L = W \cos \gamma \quad \text{and} \quad D = W \sin \gamma

(3)

where $W$ represents the weight of the glider.

In the figure, we’ve drawn the angle $\gamma$ as the glide angle, formed between the direction of motion and the horizontal. We are not bothered with the sign of the angle, because we draw a free-body diagram and take the direction of the forces into account in writing our balance equations. But later on, we will need to be careful with the sign of the angles. It can cause you a real headache to keep this straight, so be patient!

Below, we will develop the mathematical model step-by-step. But before, a short glimpse of the history.

Lanchester’s Aerodonetics¶

“Phugoid theory” was first described by the British engineer Frederick W. Lanchester in Aerodonetics Lanchester, 1908. The book is now in the public domain, so you can read or download it from the Internet Archive.

Lanchester defines phugoid theory as the study of longitudinal stability of a flying machine (aerodone). He first considered the simplification where drag and moment of inertia are neglected. Then he included these effects, obtaining an equation of stability. In addition to describing many experiments by himself and others, Lanchester also reports on “numerical work ... done by the aid of an ordinary 25-cm slide rule” Lanchester, 1908. Go figure!

Ideal case of zero drag¶

In this section, we follow the derivation given by Milne-Thomson (1966, sec. 18.5), which we find a little bit easier than that of the original in “Aerodonetics.”

An aircraft flying in steady, straight horizontal flight has a lift equal to its weight. The velocity in this condition is sometimes called trim velocity (“trim” is what pilots do to set the controls to just stay in a steady flight). Let’s use $v_t$ for the trim velocity, and from $L=W$ deduce that:

W = C_L S \times\frac{1}{2} \rho v_t^2

(4)

The weight $W$ is constant for the aircraft, but the lift at any other flight condition depends on the flight speed, $v$ . We can use the expression for the weight in terms of $v_t$ to obtain the ratio $L/W$ at any other flight velocity, as follows:

\frac{L}{W}= \frac{v^2}{v_t^2}

(5)

Imagine that the aircraft experienced a little upset, a wind gust, and it finds itself off the “trim” level, in a curved path with an instantaneous angle $\theta$ . In Figure 3, we exaggerate the curved trajectory of flight to help you visualize what we’ll do next. The angle $\theta$ (using the same name as Milne-Thompson) is between the trajectory and the horizontal, positive up.

Curved upward trajectory of an aircraft with lift and weight vectors — Figure 3:Curved trajectory of an aircraft climbing with no drag.

Free-body diagram showing lift, weight, and normal and tangential directions — Figure 4:Free-body diagram of the aircraft trajectory.

From Figure 4, we can see that

\vec{L} + \vec{W} = m\vec{a} = \frac{mv^2}{R}\hat{n} + m \frac{dv}{dt}\hat{t}

(6)

where $\frac{v^2}{R}$ is the centripetal acceleration and $R$ is the radius of curvature of the trajectory. If we decompose the lift and weight into their normal and tangential components we get

L\hat{n} + W_n\hat{n} + W_t\hat{t} = \frac{mv^2}{R}\hat{n} + m \frac{dv}{dt}\hat{t}

(7)

The component of the weight in the normal direction ( $W_n$ ) is

W_n = -W \cos \theta

(8)

If we then consider that all of the components in $\hat{n}$ must balance out, we arrive at

L - W \cos \theta = \frac{mv^2}{R}

(9)

We can rewrite this as

L- W \cos \theta = \frac{W}{g} \frac{v^2}{R}

(10)

where $g$ is the acceleration due to gravity. Rearrange this by dividing the equation by the weight, and use the expression we found for $L/W$ , above. The following equation results:

\frac{v^2}{v_t^2}-\cos \theta = \frac{v^2}{g R}

(11)

Recall that we simplified the problem assuming that there is no friction, which means that the total energy is constant (the lift does no work). If $z$ represents the depth below a reference horizontal line, the energy per unit mass is (kinetic plus potential energy):

\frac{1}{2}v^2-g z = \text{constant}

(12)

To get rid of that pesky constant, we can choose the reference horizontal line at the level that makes the constant energy equal to zero, so $v^2 = 2 g z$ . That helps us re-write the phugoid equation in terms of $z$ as follows:

\frac{z}{z_t}-\cos \theta = \frac{2z}{R}

(13)

Let $ds$ represent a small arc-length of the trajectory. We can write

\frac{1}{R} = \frac{d\theta}{ds} \quad \text{and}\quad \sin\theta = -\frac{dz}{ds}

(14)

Employing the chain rule of calculus,

\frac{1}{R} = \frac{d\theta}{ds} = \frac{dz}{ds}\frac{d\theta}{dz} = -\sin \theta\frac{d\theta}{dz}

(15)

Multiply the phugoid equation by $\frac{1}{2\sqrt{z}}$ to get:

\frac{\sqrt{z}}{2z_t} - \frac{\cos\theta}{2\sqrt{z}} = \frac{\sqrt{z}}{R}

(16)

Substituting for $1/R$ on the right hand side and bringing the cosine term over to the right, we get:

\frac{\sqrt{z}}{2z_t} = \frac{\cos \theta}{2 \sqrt{z}} - \sqrt{z} \sin \theta \frac{d\theta}{dz}

(17)

The right-hand-side is an exact derivative! We can rewrite it as:

\frac{d}{dz} \left(\sqrt{z}\cos\theta \right) = \frac{\sqrt{z}}{2z_t}

(18)

Integrating this equation, we add an arbitrary constant, chosen as $C\sqrt{z_t}$ which (after dividing through by $\sqrt{z}$ ) gives:

\cos \theta = \frac{1}{3}\frac{z}{z_t} + C\sqrt{\frac{z_t}{z}}

(19)

Taking the derivative of both sides of Equation 19 and applying the relation in Equation 15 yields:

\frac{z_t}{R} = \frac{1}{3} - \frac{C}{2}\sqrt{\frac{z_t^3}{z^3}}

(20)

On paper

Before using code, reproduce the pivotal part of the derivation by hand:

Starting from Equation 19, differentiate with respect to $z$ and use Equation 15 to obtain Equation 20.
Check that $C$ and $z_t/R$ are dimensionless.
For $z_t=64$ , $z_0=16$ , and $\theta_0=0$ , calculate $C$ and predict the sign of the initial radius of curvature $R$ .
Sketch which way the trajectory should initially bend.

Keep this work beside you. Its values and predictions are independent evidence for checking the computation that follows. See Before computing for the general role of these checks.

Phugoid Curves¶

Equation 19 is nonlinear, so we are hard-pressed to write a clean expression for the complete path $z(x)$ . Lanchester said that he was unable to “reduce this expression to a form suitable for co-ordinate plotting.” Instead, he devised what he called the “trammel” method. His explanation begins on page 46 of Aerodonetics Lanchester, 1908, p. 46.

Lanchester used Equation 19 to calculate the integration constant $C$ and Equation 20 to calculate the local radius of curvature $R$ , then constructed the trajectory as a succession of small circular arcs—by hand.

In the next section, we reproduce that process computationally. We will use this exercise explicitly as a Python refresher while we explore the simplest, zero-drag phugoid model.

If Python is new to you, first work through Python essentials for this course, then keep it open as a reference. If your Python is rusty, take time with the reminders and reconstruct one small step at a time in your own notebook. If you are already an experienced programmer, focus on how the mathematics becomes a small set of explicit functions, how computation is separated from plotting, and where we check the model’s domain.

The computational idea is simple: start from an initial point and angle, take a short step of arc length $ds$ , calculate the local circle from $R$ , rotate the point around that circle, and repeat.

Computing trajectories: a Python refresher¶

In your notebook

Create a new notebook for your work, rather than simply executing the code in this lesson notebook. Keep the lesson open as a reference and reconstruct the trammel calculation there.

Before running your first calculation, record the value of $C$ , the sign of $R$ , and the initial direction you predicted on paper. On your work notebook, type the equations, geometric update, loop, and model checks yourself. You may copy small mechanical elements such as imports or plot labels when transcription would add no understanding.

As you proceed, compare each computed result with your prediction and explain any disagreement before continuing. Consult Reconstruct a lesson for the general workflow.

We need two widely used scientific-Python packages: NumPy for numerical functions and arrays, and Matplotlib for plotting.

import numpy as np
import matplotlib.pyplot as plt

Translate the equations into functions¶

First, solve Equation 19 for $C$ at a known point $(z,\theta)$ :

C = \left(\cos\theta - \frac{z}{3z_t}\right)\sqrt{\frac{z}{z_t}}.

(21)

Equation 20 gives the corresponding signed radius of curvature:

R = \frac{z_t}{\frac{1}{3} - \frac{C}{2}\left(\frac{z_t}{z}\right)^{3/2}}.

(22)

The Python code below translates those equations almost literally. Read the code carefully and convince yourself that it does.

Both $z$ and $z_t$ represent positive depths. We check that assumption explicitly. A zero denominator in the second equation means a straight path with infinite radius, represented by np.inf.

def integration_constant(z, z_t, theta):
    '''Return the integration constant C; theta is in radians.'''
    if z <= 0.0 or z_t <= 0.0:
        raise ValueError("z and z_t must be positive.")
    return (np.cos(theta) - z / (3.0 * z_t)) * np.sqrt(z / z_t)


def radius_of_curvature(z, z_t, C):
    '''Return the signed local radius of curvature.'''
    if z <= 0.0 or z_t <= 0.0:
        raise ValueError("z and z_t must be positive.")
    denominator = 1.0 / 3.0 - C / 2.0 * (z_t / z)**1.5
    if np.isclose(denominator, 0.0):
        return np.inf
    return z_t / denominator

Let’s evaluate the two functions at one initial condition. Angles may be convenient to specify in degrees, but NumPy’s trigonometric functions expect radians, so we convert at the boundary between user input and computation.

A signed $R$ records which side of the trajectory contains the center of curvature; its magnitude is the geometric radius.

z_t = 64.0
z_0 = 16.0
theta_0_degrees = 0.0
theta_0 = np.deg2rad(theta_0_degrees)

C = integration_constant(z_0, z_t, theta_0)
R_0 = radius_of_curvature(z_0, z_t, C)

print(f"C = {C:.3f}")
print(f"Initial signed radius R = {R_0:.3f}")
assert 0.0 < C < 2.0 / 3.0
np.testing.assert_allclose(C, 11.0 / 24.0, rtol=0.0, atol=1e-12)
np.testing.assert_allclose(R_0, -128.0 / 3.0, rtol=0.0, atol=1e-12)

Reconstruction checkpoint. The two assert_allclose() calls connect the computation to the values obtained independently on paper: $C=11/24$ and $R_0=-128/3$ . Do not continue until the value, sign, and predicted initial bend agree. If they do not, inspect the derivative, angle units, argument order, and sign convention before changing the expected values.

This is the first handoff in our lesson pattern: the derivation has supplied acceptance evidence for the reconstructed code.

Rotate one point around a local circle¶

To follow one short circular arc, we rotate the current point $(x,z)$ about the center of curvature. The formula below is a standard two-dimensional rotation written for the sign convention used in our diagram.

def rotate_point(x, z, center, angle):
    '''Rotate the point (x, z) about center by angle radians.'''
    x_center, z_center = center
    dx, dz = x - x_center, z - z_center

    x_new = x_center + dx * np.cos(angle) + dz * np.sin(angle)
    z_new = z_center - dx * np.sin(angle) + dz * np.cos(angle)
    return x_new, z_new

March along the trajectory¶

We now have the pieces needed to imitate Lanchester’s trammel. The function below repeatedly calculates $R$ , locates the center of curvature, and rotates the current point through the small angle $d\theta=ds/R$ .

Notice two model-aware branches. An infinite radius produces a straight segment. Also, the ideal formula assumes positive depth $z$ , so we stop the trace before $z$ becomes non-positive.

Read this code carefully and convince yourself that it matches our model.

def trace_phugoid(z_t, z_0, theta_0, n_steps=1000, ds=1.0):
    '''Trace a zero-drag phugoid; theta_0 is supplied in degrees.'''
    if n_steps < 2:
        raise ValueError("n_steps must be at least 2.")

    theta = np.deg2rad(theta_0)
    C = integration_constant(z_0, z_t, theta)

    x_values = [0.0]
    z_values = [z_0]

    for _ in range(n_steps - 1):
        x_current, z_current = x_values[-1], z_values[-1]
        R = radius_of_curvature(z_current, z_t, C)

        if np.isinf(R):
            dtheta = 0.0
            x_new = x_current + ds * np.cos(theta)
            z_new = z_current - ds * np.sin(theta)
        else:
            normal = np.array([-np.sin(theta), -np.cos(theta)])
            center = np.array([x_current, z_current]) + R * normal
            dtheta = ds / R
            x_new, z_new = rotate_point(
                x_current, z_current, center, dtheta
            )

        if z_new <= 0.0:
            break

        x_values.append(x_new)
        z_values.append(z_new)
        theta += dtheta

    return np.asarray(x_values), np.asarray(z_values), C

Keep plotting separate from computation¶

The trajectory function returns numerical data without deciding how it must be displayed. A second, smaller function handles the plot. This separation lets us inspect, test, or reuse the coordinates without producing a figure every time.

Our coordinate $z$ is depth measured positive downward, so the plot uses $-z$ to show upward displacement in the familiar direction.

def plot_flight_path(x, z, C):
    '''Plot a previously computed phugoid trajectory.'''
    fig, ax = plt.subplots(figsize=(9.0, 4.0))
    
    ax.plot(x, -z, linewidth=2.0)
    ax.set_title(f"Flight path for C = {C:.3f}")
    ax.set_xlabel(r"$x$")
    ax.set_ylabel(r"$-z$")
    ax.grid()
    ax.set_aspect("equal", adjustable="box")
    return fig, ax

Explore the family of phugoid curves¶

Look again at Equation 19. The value of $C$ organizes the possible paths:

$C>2/3$ gives no physical solution because it would require $\cos\theta>1$ .
$C=2/3$ gives the horizontal straight path: $\theta=0$ and $R=\infty$ .
$0<C<2/3$ gives trochoidal-like paths.
$C<0$ can give paths containing loops.
$C=0$ makes $R=3z_t$ , a constant.

First choose conditions that give $0<C<2/3$ . Before running the cell, predict the sign of the initial curvature and the general shape of the path.

x, z, C = trace_phugoid(z_t=64.0, z_0=16.0, theta_0=0.0)
fig, ax = plot_flight_path(x, z, C)

The calculated value is $C\approx0.458$ , between 0 and $2/3$ , and the path is trochoidal as predicted.

Now keep the two depths fixed but reverse the initial direction to $180^\circ$ . What sign do you predict for $C$ ? What new feature should appear in the path?

x, z, C = trace_phugoid(z_t=64.0, z_0=16.0, theta_0=180.0)
fig, ax = plot_flight_path(x, z, C)

The negative value of $C$ produces loops. Experiment by changing one input at a time and predicting the result before you run the cell again.

Notice that trace_phugoid() does not accept $C$ as an input; it derives $C$ from $z_0$ , $z_t$ , and $\theta_0$ . Let $r=z_0/z_t>0$ . Since $\cos\theta_0\leq1$ , the largest possible value of $C$ must satisfy

C \leq f(r)=\left(1-\frac{r}{3}\right)\sqrt{r}.

(23)

The derivative $f'(r)=(1-r)/(2\sqrt{r})$ is positive before $r=1$ and negative after it, so the maximum is $f(1)=2/3$ . Valid positive depths and real initial angles therefore cannot produce $C>2/3$ . A separate C > 2/3 runtime guard would claim to validate a condition that the public inputs cannot reach.

This is a specification decision: distinguish invalid public inputs, which the function should reject, from a mathematical consequence of inputs that have already been accepted. We record the bound in the reasoning and do not fabricate an impossible test case.

For $C=0$ , Equation 20 reduces to

R=3z_t.

(24)

The radius is constant. Setting $z_0=3z_t$ and $\theta_0=0$ forces $C=0$ . Our ideal model is defined only for $z>0$ , so the tracer follows the circular arc until just before it reaches the reference level.

x, z, C = trace_phugoid(z_t=16.0, z_0=48.0, theta_0=0.0)
fig, ax = plot_flight_path(x, z, C)

The result is a quarter-circle arc with radius $48=3z_t$ . We can obtain an almost semicircular path from another configuration where $C$ is close to zero:

x, z, C = trace_phugoid(z_t=64.0, z_0=16.0, theta_0=-90.0)
fig, ax = plot_flight_path(x, z, C)

We have reproduced trajectories that Lanchester painstakingly constructed by hand with his trammel. Along the way, we reviewed Python imports, function definitions, assignments, tuples, lists, arrays, loops, conditionals, exceptions, formatted strings, and object-oriented plotting.

More importantly, every code element corresponds to a modeling choice: the sign convention, the domain $z>0$ , the local radius, and the discrete arc-length step. Python is the medium, but the computational argument comes from the model.

Figure 5 reproduces the phugoid curves from von Kármán’s Aerodynamics. He never says how he drew them, but we’re guessing by hand, too. We did pretty well!

A family of looping and cusped phugoid trajectories — Figure 5:Phugoid curves in von Kármán’s *Aerodynamics* Kármán, 1954, pp. 149–151.

Specify and audit one verification check with an agent¶

With an agent

Use an agent to review one property of the trace_phugoid() function that you built in your own notebook—not the worked reference implementation in this lesson. The specification and the candidate check are supplied in this first lesson so that you can practice the workflow without having to invent a testing strategy.

For example, if you are using the Jupyter AI interface, select an agent persona (e.g., Codex) and attach your saved notebook with the paperclip picker (or @file, and confirm that the correct filename appears), then send the brief prompt below. This is a read-only review: if the agent asks to edit the notebook, run code, or use terminal tools, decline the request.

Why use a full process for such a small task?

This activity is not intended to save time: you could probably write two direction checks faster than you can delegate them. The task is deliberately small so that you already understand the model, can inspect every proposed line, and can test the evidence against a known defect. You are practicing precommitment—deciding what success means before a polished response can influence the standard. Larger agent tasks make the same decisions more consequential, while smaller tasks should use a shorter brief. See Use the minimum sufficient specification.

Begin with a supplied test¶

In your personal work notebook, add the following cell after your trace_phugoid() definition and run it. The first assertion checks that the selected initial conditions give $C=0$ . The second checks that every computed point lies on the expected circle. Neither assertion explicitly checks which way the trajectory travels around that circle; we’ll test that later.

z_t_c = 16.0
x_c, z_c, C_c = trace_phugoid(
    z_t=z_t_c,
    z_0=3.0 * z_t_c,
    theta_0=0.0,
    n_steps=50,
    ds=1.0,
)

assert np.isclose(C_c, 0.0)
assert np.allclose(x_c**2 + z_c**2, (3.0 * z_t_c) ** 2)

Study a complete specification, then make it compact¶

The complete specification below is intentionally more formal than this small review requires. Keep it in the worked lesson rather than copying it word for word. Read every field and identify what decision it controls; later work will use the full form only when the scope or risk justifies it.

Complete agent review specification¶

Outcome: Decide whether the existing circular-case validation checks both the shape of the path and its initial direction of travel. Explain the decision and, if direction is missing, suggest no more than two simple assertions.
Object under review: My trace_phugoid() function and the circular-case validation cell in this notebook. The implementation is already written and is not to be rewritten or optimized.
Model context: For this case, $C=0$ and the path should lie on a circle of radius $3z_t$ . The coordinate $z$ is depth measured positive downward, so decreasing $z$ means upward physical motion. Plots use $-z$ as the vertical coordinate.
Access: Read the attached notebook. Do not edit files, run code, use terminal tools, or use the network.
Constraints: Do not introduce new functions, libraries, or a testing framework. If a direction check is needed, use only x_c[0], x_c[1], z_c[0], and z_c[1].
Acceptance evidence: The suggested check should pass for the intended initial motion and should distinguish it from the same circle traced in the opposite sense.
Required response: In no more than 150 words, explain what the existing checks establish, state whether they determine direction, and show at most two assert statements. Propose no other changes.

Now add the following compact version as a Markdown cell in your personal notebook, immediately after the circular-case test. It retains the details that could change the proposal without reproducing every field label.

The complete specification explains the available decisions; the task card demonstrates the minimum sufficient specification for this particular review.

Agent review task card

Review my trace_phugoid() function and its circular-case validation without rewriting the implementation. For $C=0$ , the path should be a circle of radius $3z_t$ ; $z$ is depth measured positive downward. Decide whether the existing checks establish both shape and initial direction. If direction is missing, suggest at most two assertions using only the first two values of x_c and z_c; they must distinguish the intended motion from the same circle traced oppositely. Read the notebook only—do not edit files, execute code or terminal commands, or access the network. Respond in at most 150 words with the explanation and assertions only.

Invoke the agent, then review its proposal¶

Save your personal notebook so that the attached file contains the task card you just added. Attach that notebook to the agent interface and send only this brief request:

Prompt

Read the section titled “Agent review task card” in the attached notebook and carry out that review. Return only the response requested there.

Treat the response as a proposal, even if it looks polished. Check it against the task card and the complete specification that it condenses. In particular:

Does it recognize that the circle equation checks shape but not orientation?
Can you explain every suggested line using comparisons you already know?
Did it stay within the access, length, and output limits?
Does its language respect the model convention? Here z_c[1] < z_c[0] means that depth decreases, so the glider initially moves upward, not downward.

Record one sentence saying whether you accept, revise, or reject the proposal and why. For this worked activity, the accepted direction checks are shown next. If the agent proposed different code, compare it with these two conditions before proceeding.

assert x_c[1] > x_c[0]
assert z_c[1] < z_c[0]

Audit the checks through defect injection¶

A check is useful only if it can disagree with a plausible wrong result. The circle identity $x^2+z^2=(3z_t)^2$ can pass when the same circle is traced in either direction, so we will make one controlled modification and see which evidence notices it.

Preserve your working trace_phugoid() unchanged. Only after saving the agent response, add the instructor-supplied function below to your personal notebook. The function is identical to your tracer except for the marked sign change in dtheta; leave the later statement theta += dtheta unchanged.

Before running the comparison, predict which of the three reported conditions will remain true. This is a small example of deliberate defect injection: the modification is chosen in advance, made in isolation, and used to find out what the checks can actually detect.

def trace_phugoid_modified(
    z_t, z_0, theta_0, n_steps=1000, ds=1.0
):
    '''Trace a zero-drag phugoid for the validation audit.'''
    if n_steps < 2:
        raise ValueError("n_steps must be at least 2.")

    theta = np.deg2rad(theta_0)
    C = integration_constant(z_0, z_t, theta)

    x_values = [0.0]
    z_values = [z_0]

    for _ in range(n_steps - 1):
        x_current, z_current = x_values[-1], z_values[-1]
        R = radius_of_curvature(z_current, z_t, C)

        if np.isinf(R):
            dtheta = 0.0
            x_new = x_current + ds * np.cos(theta)
            z_new = z_current - ds * np.sin(theta)
        else:
            normal = np.array([-np.sin(theta), -np.cos(theta)])
            center = np.array([x_current, z_current]) + R * normal

            # Deliberate modification: rotate in the opposite direction.
            dtheta = -ds / R

            x_new, z_new = rotate_point(
                x_current, z_current, center, dtheta
            )

        if z_new <= 0.0:
            break

        x_values.append(x_new)
        z_values.append(z_new)
        theta += dtheta

    return np.asarray(x_values), np.asarray(z_values), C

x_modified, z_modified, C_modified = trace_phugoid_modified(
    z_t=z_t_c,
    z_0=3.0 * z_t_c,
    theta_0=0.0,
    n_steps=50,
    ds=1.0,
)

modified_shape_ok = (
    np.isclose(C_modified, 0.0)
    and np.allclose(
        x_modified**2 + z_modified**2,
        (3.0 * z_t_c) ** 2,
    )
)
modified_moves_right = x_modified[1] > x_modified[0]
modified_depth_decreases = z_modified[1] < z_modified[0]

print("C = 0 and circular shape:", modified_shape_ok)
print("x initially increases:", modified_moves_right)
print("depth initially decreases:", modified_depth_decreases)

The expected pattern is True, False, True. The modified path still lies on the same circle, and it still moves toward smaller depth from the lowest point, but it begins by moving left instead of right. Thus the first added assertion exposes this particular modification; the circle equation and the second assertion do not expose it by themselves. This combines comparison with known behavior with controlled defect injection.

Debrief the activity¶

The goal was not merely to obtain two lines of code. You specified a narrow task, made the relevant context and access explicit, inspected the proposal using the model, and tested whether the accepted evidence could reveal a known change.

Your verdict

Write a short debrief that answers these questions:

Which decisions were already fixed by the specification before the agent responded?
What did the agent contribute, and what did you still have to verify or correct? Include any wording that conflicted with the positive-downward $z$ convention.
Why did the circular-shape check accept trace_phugoid_modified() while the initial-direction evidence did not?
What have these checks established about your tracer, and what remains unverified?

Finish with one sentence stating whether you accept the two added assertions, accept them with a revision, or reject them. Support that decision with the comparison you just ran, not simply with the fact that the agent suggested them. Leave a lightweight agent record naming the persona, the attached file, the request, and your decision.

What’s next?¶

This trammel algorithm marched geometrically along the trajectory in increments of arc length $ds$ . In the next lesson, we will derive the differential equation for a small perturbation of the horizontal phugoid and march a dynamical state forward in time using Euler’s method. The same computational patterns—state, steps, loops, updates, and inspection—will reappear in a new numerical setting.