Modeling the launch vehicle
For this exercise, we’ll build a simulation of an Atlas V launch vehicle in a new file, src/scratch/rocket.clj. We’ll stick to SI units and keep the model simple: constant thrust, no atmospheric effects, and spherical coordinates for position.
The rocket’s key parameters are its dry mass (50,050 kg), fuel mass (about 284,450 kg, the conversion from 627,105 lbs), thrust (4,152 kN), burn time (253 seconds), and specific impulse (3.05 km/s, which we treat as exhaust velocity). Here’s the physical setup:
(defn atlas-v
[]
{:dry-mass 50050
:fuel-mass 284450
:time 0
:isp 3050
:max-fuel-rate (/ 284450 253)
:max-thrust 4.152e6})
The dry mass is the non-fuel mass, so the total vehicle mass is the sum of the two:
(defn mass
"The total mass of a craft."
[craft]
(+ (:dry-mass craft) (:fuel-mass craft)))
We’ll use a geocentric Cartesian coordinate system: x and y lie in the equatorial plane, z points toward the north pole, and the origin is the planet’s center. That choice makes forces and accelerations straightforward to compute, but positions and velocities are easier to express in spherical coordinates. We’ll convert back and forth as needed.
Our launch site sits on the equator at y = 0. The earth’s equatorial radius is about 6,378 km, so we can set up the initial position. The rotate speed at the equator comes from the earth’s rotation rate:
(def earth-day
"Length of an earth day, in seconds."
86400)
(def earth-equatorial-speed
"How fast points on the equator move, relative to the center of the earth,
in meters/sec."
(/ (* 2 Math/PI earth-equatorial-radius)
earth-day))
At the equator, that speed is entirely in the y direction. The initial conditions for the rocket at time zero are assembled by merging those values with the craft’s map:
(def initial-space-center
"The initial position and velocity of the launch facility"
{:time 0
:position {:x earth-equatorial-radius
:y 0
:z 0}
:velocity {:x 0
:y earth-equatorial-speed
:z 0}})
(defn prepare
"Prepares a craft for launch from an equatorial space center."
[craft]
(merge craft initial-space-center))
Note that we’re using maps with :x, :y, and :z keys to represent vectors—these are logical coordinates, not Clojure vectors.
Forces and motion
Gravity pulls the craft toward the planet’s center at 9.8 m/s². To compute that direction, we need the spherical coordinate angles of the rocket’s position, which we derive using trigonometric functions from java.lang.Math:
(defn magnitude
"What's the radius of a given set of cartesian coordinates?"
[c]
; By the Pythagorean theorem...
(Math/sqrt (+ (Math/pow (:x c) 2)
(Math/pow (:y c) 2)
(Math/pow (:z c) 2))))
(defn cartesian->spherical
"Converts a map of Cartesian coordinates :x, :y, and :z to spherical coordinates :r, :theta, and :phi."
[c]
(let [r (magnitude c)]
{:r r
:theta (Math/acos (/ (:z c) r))
:phi (Math/atan2 (:y c) (:x c))}))
(defn spherical->cartesian
"Converts spherical to Cartesian coordinates."
[c]
{:x (* (:r c) (Math/sin (:theta c)) (Math/cos (:phi c)))
:y (* (:r c) (Math/sin (:theta c)) (Math/sin (:phi c)))
:z (* (:r c) (Math/cos (:phi c)))})
Then gravity is simply the spherical vector with the radius replaced by the force magnitude, converted back to Cartesian:
(def g "Acceleration of gravity in meters/s^2" -9.8)
(defn gravity-force
"The force vector, each component in Newtons, due to gravity."
[craft]
; Since force is mass times acceleration...
(let [total-force (* g (mass craft))]
(-> craft
; Now we'll take the craft's position
:position
; in spherical coordinates,
cartesian->spherical
; replace the radius with the gravitational force...
(assoc :r total-force)
; and transform back to Cartesian-land
spherical->cartesian)))
Thrust burns fuel at a constant maximum rate until it’s exhausted:
(defn fuel-rate
"How fast is fuel, in kilograms/second, consumed by the craft?"
[craft]
(if (pos? (:fuel-mass craft))
(:max-fuel-rate craft)
0))
The thrust force is the fuel burn rate times the exhaust velocity (:isp). For fun, we’ll point the engine entirely along the x axis:
(defn thrust
"How much force, in newtons, does the craft's rocket engines exert?"
[craft]
(* (fuel-rate craft) (:isp craft)))
(defn engine-force
"The force vector, each component in Newtons, due to the rocket engine."
[craft]
(let [t (thrust craft)]
{:x t
:y 0
:z 0}))
Total force is gravity plus thrust. We can combine the two maps with merge-with to add corresponding components:
(defn total-force
"Total force on a craft."
[craft]
(merge-with + (engine-force craft)
(gravity-force craft)))
Acceleration comes from Newton’s second law: force divided by mass. Since our coordinates are maps, we need a utility to apply a function to each value. Map entries are key/value pairs, and into can rebuild a map from those pairs:
user=> (seq {:x 1 :y 2 :z 3})
([:z 3] [:y 2] [:x 1])
user=> (into {} [[:x 4] [:y 5]])
{:x 4, :y 5}
That gives us a map-values function that is like map for map values. With it, a scale function multiplies coordinates by a factor. The partial operator creates the multiplier function:
(defn map-values
"Applies f to every value in the map m."
[f m]
(into {}
(map (fn [pair]
[(key pair) (f (val pair))])
m)))
(defn scale
"Multiplies a map of x, y, and z coordinates by the given factor."
[factor coordinates]
(map-values (partial * factor) coordinates))
Dividing by mass is just scaling by the reciprocal:
(defn acceleration
"Total acceleration of a craft."
[craft]
(let [m (mass craft)]
(scale (/ m) (total-force craft))))
Now we can step the rocket forward in time. A step function takes the current state and returns a new state dt seconds later:
(defn step
[craft dt]
(assoc craft
; Time advances by dt seconds
:t (+ dt (:t craft))
; We burn some fuel
:fuel-mass (- (:fuel-mass craft) (* dt (fuel-rate craft)))
; Our position changes based on our velocity
:position (merge-with + (:position craft)
(scale dt (:velocity craft)))
; And our velocity changes based on our acceleration
:velocity (merge-with + (:velocity craft)
(scale dt (acceleration craft)))))
Debugging the first launch
Load the code into the REPL and launch. The use form with :reload forces a fresh read of the file:
user=> (use 'scratch.rocket :reload)
nil
user=> (atlas-v)
{:dry-mass 50050, :fuel-mass 284450, :time 0, :isp 3050, :max-fuel-rate 284450/253, :max-thrust 4152000.0}
Preparing the rocket on the pad shows an interesting detail: even “at rest,” it’s moving at 463 m/s eastward from the earth’s rotation. Stepping forward one second, however, throws a NullPointerException:
user=> (-> (atlas-v) prepare pprint)
{:velocity {:x 0, :y 463.8312116386399, :z 0},
:position {:x 6378137, :y 0, :z 0},
:dry-mass 50050,
:fuel-mass 284450,
:time 0,
:isp 3050,
:max-fuel-rate 284450/253,
:max-thrust 4152000.0}
user=> (-> (atlas-v) prepare (step 1) pprint)
NullPointerException clojure.lang.Numbers.ops (Numbers.java:942)
A stack trace (pst) shows the failure originates in our step function, which references a missing :t field—we named it :time in the craft map:
user=> (pst *e)
NullPointerException
clojure.lang.Numbers.ops (Numbers.java:942)
clojure.lang.Numbers.add (Numbers.java:126)
scratch.rocket/step (rocket.clj:125)
user/eval1478 (NO_SOURCE_FILE:1)
clojure.lang.Compiler.eval (Compiler.java:6619)
clojure.lang.Compiler.eval (Compiler.java:6582)
clojure.core/eval (core.clj:2852)
clojure.main/repl/read-eval-print--6588/fn--6591 (main.clj:259)
clojure.main/repl/read-eval-print--6588 (main.clj:259)
clojure.main/repl/fn--6597 (main.clj:277)
clojure.main/repl (main.clj:277)
clojure.tools.nrepl.middleware.interruptible-eval/evaluate/fn--589 (interruptible_eval.clj:56)
123 (assoc craft
124 ; Time advances by dt seconds
125 :t (+ dt (:t craft))
After fixing the field name and reloading, the velocity shifts by -9.8 m/s — but south, not down. That points to a bug in the spherical-to-Cartesian conversion. A unit test on a point in the equatorial plane reveals swapped angles:
(ns scratch.rocket-test
(:require [clojure.test :refer :all]
[scratch.rocket :refer :all]))
(deftest spherical-coordinate-test
(let [pos {:x 1 :y 2 :z 3}]
(testing "roundtrip"
(is (= pos (-> pos cartesian->spherical spherical->cartesian))))))
aphyr@waterhouse:~/scratch$ lein test
lein test scratch.core-test
lein test scratch.rocket-test
lein test :only scratch.rocket-test/spherical-coordinate-test
FAIL in (spherical-coordinate-test) (rocket_test.clj:8)
roundtrip
expected: (= pos (-> pos cartesian->spherical spherical->cartesian))
actual: (not (= {:z 3, :y 2, :x 1} {:x 1.0, :y 1.9999999999999996, :z 1.6733200530681513}))
Ran 2 tests containing 4 assertions.
1 failures, 0 errors.
Tests failed.
(deftest spherical-coordinate-test
(testing "spherical->cartesian"
(is (= (spherical->cartesian {:r 2
:phi 0
:theta 0})
{:x 0.0 :y 0.0 :z 2.0})))
(testing "roundtrip"
(let [pos {:x 1.0 :y 2.0 :z 3.0}]
(is (= pos (-> pos cartesian->spherical spherical->cartesian))))))
user=> (cartesian->spherical {:x 0.00001 :y 0.00001 :z 2.0})
{:r 2.00000000005, :theta 7.071068104411588E-6, :phi 0.7853981633974483}
user=> (cartesian->spherical {:x 1 :y 2 :z 3})
{:r 3.7416573867739413, :theta 0.6405223126794245, :phi 1.1071487177940904}
user=> (spherical->cartesian (cartesian->spherical {:x 1 :y 2 :z 3}))
{:x 1.0, :y 1.9999999999999996, :z 1.6733200530681513}
user=> (cartesian->spherical {:x 1 :y 2 :z 0})
{:r 2.23606797749979, :theta 1.5707963267948966, :phi 1.1071487177940904}
user=> (cartesian->spherical {:x 1 :y 1 :z 0})
{:r 1.4142135623730951, :theta 1.5707963267948966, :phi 0.7853981633974483}
Cross-checking with reference formulas confirms phi (angle from the pole) and theta (equatorial angle) were reversed. The corrected conversion resolves it:
(defn cartesian->spherical
"Converts a map of Cartesian coordinates :x, :y, and :z to spherical
coordinates :r, :theta, and :phi."
[c]
(let [r (Math/sqrt (+ (Math/pow (:x c) 2)
(Math/pow (:y c) 2)
(Math/pow (:z c) 2)))]
{:r r
:phi (Math/acos (/ (:z c) r))
:theta (Math/atan2 (:y c) (:x c))}))
aphyr@waterhouse:~/scratch$ lein test
lein test scratch.core-test
lein test scratch.rocket-test
Ran 2 tests containing 5 assertions.
0 failures, 0 errors.
With that fixed, stepping forward now accelerates the rocket in the +x direction. Liftoff is achieved.
Simulating the full flight
Rather than stepping one second at a time, we can build an infinite lazy sequence of states with iterate, where each frame is one step past the previous:
(defn trajectory
[dt craft]
"Returns all future states of the craft, at dt-second intervals."
(iterate #(step % 1) craft))
user=> (->> (atlas-v) prepare (trajectory 1) (take 3) pprint)
({:velocity {:x 0, :y 463.8312116386399, :z 0},
:position {:x 6378137, :y 0, :z 0},
:dry-mass 50050,
:fuel-mass 284450,
:time 0,
:isp 3050,
:max-fuel-rate 284450/253,
:max-thrust 4152000.0}
{:velocity
{:x 0.45154055666826204,
:y 463.8312116386399,
:z -6.000769315822031E-16},
:position {:x 6378137, :y 463.8312116386399, :z 0},
:dry-mass 50050,
:fuel-mass 71681400/253,
:time 1,
:isp 3050,
:max-fuel-rate 284450/253,
:max-thrust 4152000.0}
{:velocity
{:x 0.9376544222659078,
:y 463.83049896253056,
:z -1.200153863164406E-15},
:position
{:x 6378137.451540557,
:y 927.6624232772798,
:z -6.000769315822031E-16},
:dry-mass 50050,
:fuel-mass 71396950/253,
:time 2,
:isp 3050,
:max-fuel-rate 284450/253,
:max-thrust 4152000.0})
Each element is a snapshot at one-second intervals. The altitude above the surface is the spherical radius minus the earth’s radius:
(defn altitude
"The height above the surface of the equator, in meters."
[craft]
(-> craft
:position
cartesian->spherical
:r
(- earth-equatorial-radius)))
We can plot the altitude at successive times, but the meaningful question is whether the rocket breaks orbit or crashes. Because the trajectory is infinite, we can’t check the entire sequence. We’ll instead assume the craft should crash within the first 100 hours; if it survives past that, we call it a successful orbit. That logic goes into a test in test/scratch/rocket_test.clj:
user=> (->> (atlas-v) prepare (trajectory 1) (map altitude) (take 10) pprint)
(0.0
0.016865378245711327
0.519002066925168
1.540983198210597
3.117615718394518
5.283942770212889
8.075246102176607
11.52704851794988
15.675116359256208
20.555462017655373)
(defn above-ground?
"Is the craft at or above the surface?"
[craft]
(<= 0 (altitude craft)))
(defn flight
"The above-ground portion of a trajectory."
[trajectory]
(take-while above-ground? trajectory))
(defn crashed?
"Does this trajectory crash into the surface before 100 hours are up?"
[trajectory]
(let [time-limit (* 100 3600)] ; 1 hour
(not (every? above-ground?
(take-while #(<= (:time %) time-limit) trajectory)))))
(defn crash-time
"Given a trajectory, returns the time the rocket impacted the ground."
[trajectory]
(:time (last (flight trajectory))))
(defn apoapsis
"The highest altitude achieved during a trajectory."
[trajectory]
(apply max (map altitude trajectory)))
(defn apoapsis-time
"The time of apoapsis"
[trajectory]
(:time (apply max-key altitude (flight trajectory))))
(deftest makes-orbit
(let [trajectory (->> (atlas-v)
prepare
(trajectory 1))]
(when (crashed? trajectory)
(println "Crashed at" (crash-time trajectory) "seconds")
(println "Maximum altitude" (apoapsis trajectory)
"meters at" (apoapsis-time trajectory) "seconds"))
; Assert that the rocket eventually made it to orbit.
(is (not (crashed? trajectory)))))
aphyr@waterhouse:~/scratch$ lein test scratch.rocket-test
lein test scratch.rocket-test
Crashed at 982 seconds
Maximum altitude 753838.039645385 meters at 532 seconds
lein test :only scratch.rocket-test/makes-orbit
FAIL in (makes-orbit) (rocket_test.clj:26)
expected: (not (crashed? trajectory))
actual: (not (not true))
Ran 2 tests containing 3 assertions.
1 failures, 0 errors.
Tests failed.
The test finds the rocket reaches 750 km altitude but crashes back to earth 982 seconds after launch. More thrust will be required.
Adding a Second Stage
The Atlas V needs help to reach orbit. Real launch vehicles carry a smaller, more efficient upper stage—the Centaur—that sits atop the main booster and takes over once the first stage burns out.
(defn centaur
"The upper rocket stage.
http://en.wikipedia.org/wiki/Centaur_(rocket_stage)
http://www.astronautix.com/stages/cenaurde.htm"
[]
{:dry-mass 2361
:fuel-mass 13897
:isp 4354
:max-fuel-rate (/ 13897 470)})
We can model this by making our atlas-v function accept the next stage as an argument.
(defn atlas-v
"The full launch vehicle. http://en.wikipedia.org/wiki/Atlas_V"
[next-stage]
{:dry-mass 50050
:fuel-mass 284450
:isp 3050
:max-fuel-rate (/ 284450 253)
:next-stage next-stage})
Our launch test then constructs the vehicle with both stages in place:
(let [trajectory (->> (atlas-v (centaur))
prepare
(trajectory 1))]
When the primary stage runs dry, the vehicle needs to shed that dead weight. We introduce a stage function that checks the vehicle's state and separates the booster, letting the Centaur continue the flight:
(defn stage
"When fuel reserves are exhausted, separate stages. Otherwise, return craft
unchanged."
[craft]
(cond
; Still fuel left
(pos? (:fuel-mass craft))
craft
; No remaining stages
(nil? (:next-stage craft))
craft
; Stage!
:else
(merge (:next-stage craft)
(select-keys craft [:time :position :velocity]))))
The cond in stage covers three situations: fuel remaining, no second stage to separate, and the moment of separation itself. In that final case we return the Centaur's state but carry over the current time, position, and velocity. We also need to write the separation check into our physics loop:
(defn step
[craft dt]
(let [craft (stage craft)]
(assoc craft
; Time advances by dt seconds
:time (+ dt (:time craft))
; We burn some fuel
:fuel-mass (- (:fuel-mass craft) (* dt (fuel-rate craft)))
; Our position changes based on our velocity
:position (merge-with + (:position craft)
(scale dt (:velocity craft)))
; And our velocity changes based on our acceleration
:velocity (merge-with + (:velocity craft)
(scale dt (acceleration craft))))))
Re-running the launch shows the extra stage buys us a much higher trajectory—our apoapsis has grown from 750 kilometers to 4,598 kilometers. Still, the rocket crashes down. The problem is no longer altitude; it's the flight path.
Horizontal Velocity
The rocket flies straight up and comes straight back down. An orbit requires sideways motion around the planet, produced by a second engine burn once the vehicle has climbed high enough. We don't need a perfect circularization maneuver—any trajectory that misses the ground counts as an orbit. The key is burning toward the horizon.
That requires finer control over the engine. Instead of a single thrust vector that always pushes in the +x direction, we break engine output into thrust (magnitude) and orientation (direction):
(defn unit-vector
"Scales coordinates to magnitude 1."
[coordinates]
(scale (/ (magnitude coordinates)) coordinates))
(defn engine-force
"The force vector, each component in Newtons, due to the rocket engine."
[craft]
(scale (thrust craft) (unit-vector (orientation craft))))
The code normalizes the orientation vector to length one with unit-vector, then scales it by the thrust magnitude to produce a thrust vector. If you reorganize namespaces as you go, you might see an error like this one when a function is used before it's defined:
Exception in thread "main" java.lang.RuntimeException: Unable to resolve symbol: unit-vector in this context, compiling:(scratch/rocket.clj:69:11)
at clojure.lang.Compiler.analyze(Compiler.java:6380)
at clojure.lang.Compiler.analyze(Compiler.java:6322)
Clojure resolves symbols in order as the file loads. When a needed dependency appears later in a file, you have two options: order simple functions first so that each name has a meaning before it is referenced, or place a (declare unit-vector) form near the top to notify the compiler that the definition comes later in the file.
With the engine force separated out, we define the burn schedule as start/end time pairs in seconds:
(def ascent
"The start and end times for the ascent burn."
[0 3000])
(def circularization
"The start and end times for the circularization burn."
[4000 1000])
The burn logic itself changes from "full throttle until empty" to two explicit maneuvers:
(defn fuel-rate
"How fast is fuel, in kilograms/second, consumed by the craft?"
[craft]
(cond
; Out of fuel
(<= (:fuel-mass craft) 0)
0
; Ascent burn
(<= (first ascent) (:time craft) (last ascent))
(:max-fuel-rate craft)
; Circularization burn
(<= (first circularization) (:time craft) (last circularization))
(:max-fuel-rate craft)
; Shut down engines otherwise
:else 0))
This cond tests for three states: propellant exhausted, burning in the ascent or circularization window, and engines off. Because Clojure's <= accepts any number of arguments, expressing time intervals such as "between the first and last ascent times" is straightforward.
The final piece is choosing which direction to burn. That takes some vector math: we need a heading tangential to the planet, not toward or away from it. Taking the rocket's velocity and removing the component pointing toward or away from Earth leaves the component that travels around the globe:
(defn dot-product
"Finds the inner product of two x, y, z coordinate maps.
See http://en.wikipedia.org/wiki/Dot_product."
[c1 c2]
(+ (* (:x c1) (:x c2))
(* (:y c1) (:y c2))
(* (:z c1) (:z c2))))
(defn projection
"The component of coordinate map a in the direction of coordinate map b.
See http://en.wikipedia.org/wiki/Vector_projection."
[a b]
(let [b (unit-vector b)]
(scale (dot-product a b) b)))
(defn rejection
"The component of coordinate map a *not* in the direction of coordinate map
b."
[a b]
(let [a' (projection a b)]
{:x (- (:x a) (:x a'))
:y (- (:y a) (:y a'))
:z (- (:z a) (:z a'))}))
The orientation logic then mirrors the burn schedule:
(defn orientation
"What direction is the craft pointing?"
[craft]
(cond
; Initially, point along the *position* vector of the craft--that is
; to say, straight up, away from the earth.
(<= (first ascent) (:time craft) (last ascent))
(:position craft)
; During the circularization burn, we want to burn *sideways*, in the
; direction of the orbit. We'll find the component of our velocity
; which is aligned with our position vector (that is to say, the vertical
; velocity), and subtract the vertical component. All that's left is the
; *horizontal* part of our velocity.
(<= (first circularization) (:time craft) (last circularization))
(rejection (:velocity craft) (:position craft))
; Otherwise, just point straight ahead.
:else (:velocity craft)))
During ascent we push directly away from the planet's center. For circularization we use the rejection function to isolate the horizontal component of velocity and burn along it. As a default, the vehicle points in its current direction of travel.
With these changes, the test launches shows two distinct burns and finally a successful orbit.
aphyr@waterhouse:~/scratch$ lein test scratch.rocket-test
lein test scratch.rocket-test
Ran 2 tests containing 3 assertions.
0 failures, 0 errors.
A Structured Simulation
Looking back at the complete program, the pieces build from abstract foundations to a working flight. It opens with linear algebra utilities for coordinates, constants representing the simulated universe, and the initial spacecraft definition. Tooling for the rocket comes next: burn timing, throttle logic, orientation, and staging. From those pieces flows the physics engine applying gravity and thrust, integrated with the Euler method, and finally the analysis functions that answer whether the vehicle reached orbit, how long the flight took, and whether the rocket exploded.
The design leans on immutable maps representing the state of the vehicle, loaded with pure functions that interpret state and produce new state. The iterate function expands a single moment into a lazy history of every future frame, which the analysis functions consume. That trade-off—performance for clarity—makes the cause and effect within the simulation much more visible.
The approach exposes choices unique to Clojure's style. In a language with object-oriented idioms, that state would probably be scattered among classes: Craft with Atlas and Centaur subclasses, a Coordinate type split into Cartesian and Spherical variants. The added structure prevents mixing types by accident, but introduces rigidity for little benefit in a simulation like this. Those languages encourage mutating one state object in place; synchronization concerns get handled by scattering variable copies through the code, with a price in clarity and a modest gain in memory use. Likewise, avoiding lazy sequences would force us to precompute a state history and manually step through indices to find crash events and maxima.
For all the differences, the fundamentals carry over between approaches. Break a large problem into parts, store each piece of the world in a structure, and write functions that describe how to progress from one state to the next. Clear docstrings, narrative structure, and unit tests keep the story of a program legible as it grows.
Exercises
- The plotted trajectory has both the maximum and minimum altitude documented, but not their orbital locations. Because the orbit history is infinite, you cannot apply
maxto it directly. Instead, find the apoapsis by locating the frame in which the altitude delta stops increasing, then write a correspondingperiapsisfunction. Report both values in the test output. - Force of gravity in the current simulation is constant acceleration toward the Earth. Real gravity falls off with the square of distance. Update the gravitational component of the force model using Earth's mass, the craft's mass, and the gravitational constant. Compare the resulting apoapsis to the original.
- Air drag is not included. Add an altitude-dependent density function to the simulation, assume the response grows with the square of the rocket's velocity, and make a basic guess at the numbers. Does the flight still reach orbit?
- The periapsis and apoapsis of the final orbit are not equal. Tune the circularization burn condition—try values near the International Space Station's altitude and orbiting speed—so that the two altitudes are close enough that the eccentricity dips below 0.2.
