Clojure from the ground up: modeling

Until this point in the book, we’ve dealt primarily in specific details: what an expression is, how math works, which functions apply to different data structures, and where code lives. But programming, like speaking a language, painting landscapes, or designing turbines, is about more than the nuts and bolts of the trade. It’s knowing how to combine those parts into a cohesive whole–and this is a skill which is difficult to describe formally. In this part of the book, I’d like to work with you on an integrative tour of one particular problem: modeling a rocket in flight.

We’re going to reinforce our concrete knowledge of the standard library by using maps, sequences, and math functions together. At the same time, we’re going to practice how to represent a complex system; decomposing a problem into smaller parts, naming functions and variables, and writing tests.

So you want to go to space

First, we need a representation of a craft. The obvious properties for a rocket are its dry mass (how much it weighs without fuel), fuel mass, position, velocity, and time. We’ll create a new file in our scratch project–src/scratch/rocket.clj–to talk about spacecraft.

(defn mass"The total mass of a craft."[craft](+ (:dry-masscraft)(:fuel-masscraft)))

What about the position and velocity? We could represent them in Cartesian coordinates–x, y, and z–or we could choose spherical coordinates: a radius from the planet and angle from the pole and 0 degrees longitude. I’ve got a hunch that spherical coordinates will be easier for position, but accelerating the craft will be simplest in in x, y, and z terms. The center of the planet is a natural choice for the coordinate system’s origin (0, 0, 0). We’ll choose z along the north pole, and x and y in the plane of the equator.

Let’s define a space center where we launch from–let’s say it’s initially on the equator at y=0. To figure out the x coordinate, we’ll need to know how far the space center is from the center of the earth. The earth’s equatorial radius is ~6378 kilometers.

(def earth-equatorial-radius"Radius of the earth, in meters"6378137)

How fast is the surface moving? Well the earth’s day is 86,400 seconds long,

(def earth-day"Length of an earth day, in seconds."86400)

which means a given point on the equator has to go 2 * pi * equatorial radius meters in earth-day seconds:

(def earth-equatorial-speed"How fast points on the equator move, relative to the center of the earth, in meters/sec."(/ (* 2Math/PIearth-equatorial-radius)earth-day))

So our space center is on the equator (z=0), at y=0 by choice, which means x is the equatorial radius. Since the earth is spinning, the space center is moving at earth-equatorial-speed in the y direction–and not changing at all in x or z.

(def initial-space-center"The initial position and velocity of the launch facility"{:time0:position{:xearth-equatorial-radius:y0:z0}:velocity{:x0:yearth-equatorial-speed:z0}})

:position and :velocity are both vectors, in the sense that they describe a position, or a direction, in terms of x, y, and z components. This is a different kind of vector than a Clojure vector, like [1 2 3]. We’re actually representing these logical vectors as Clojure maps, with :x, :y, and :z keys, corresponding to the distance along the x, y, and z directions, from the center of the earth. Throughout this chapter, I’ll mainly use the term coordinates to talk about these structures, to avoid confusion with Clojure vectors.

Now let’s create a function which positions our craft on the launchpad at time 0. We’ll just merge the spacecraft’s with the initial space center, overwriting the craft’s time and space coordinates.

(defn prepare"Prepares a craft for launch from an equatorial space center."[craft](merge craftinitial-space-center))

Forces

Gravity continually pulls the spacecraft towards the center of the Earth, accelerating it by 9.8 meters/second every second. To figure out what direction is towards the Earth, we’ll need the angles of a spherical coordinate system. We’ll use the 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(:xc)2)(Math/pow(:yc)2)(Math/pow(:zc)2))))(defn cartesian->spherical"Converts a map of Cartesian coordinates :x, :y, and :z to spherical coordinates :r, :theta, and :phi."[c](let [r(magnitudec)]{:rr:theta(Math/acos(/ (:zc)r)):phi(Math/atan(/ (:yc)(:xc)))}))(defn spherical->cartesian"Converts spherical to Cartesian coordinates."[c]{:x(* (:rc)(Math/sin(:thetac))(Math/cos(:phic))):y(* (:rc)(Math/sin(:thetac))(Math/sin(:phic))):z(* (:rc)(Math/cos(:phic)))})

With those angles in mind, computing the gravitational acceleration is easy. We just take the spherical coordinates of the spacecraft, and replace the radius with the total force due to gravity. Then we can transform that spherical force back into Cartesian coordinates.

(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(masscraft))](-> craft; Now we'll take the craft's position:position; in spherical coordinates,cartesian->spherical; replace the radius with the gravitational force...(assoc :rtotal-force); and transform back to Cartesian-landspherical->cartesian)))

Rockets produce thrust by consuming fuel. Let’s say the fuel consumption is always the maximum, until we run out:

(defn fuel-rate"How fast is fuel, in kilograms/second, consumed by the craft?"[craft](if (pos? (:fuel-masscraft))(:max-fuel-ratecraft)0))

Now that we know how much fuel is being consumed, we can compute the force the rocket engine develops. That force is simply the mass consumed per second times the exhaust velocity–which is the specific impulse :isp. We’ll ignore atmospheric effects.

(defn thrust"How much force, in newtons, does the craft's rocket engines exert?"[craft](* (fuel-ratecraft)(:ispcraft)))

Cool. What about the direction of thrust? Just for grins, let’s keep the rocket pointing entirely along the x axis.

(defn engine-force"The force vector, each component in Newtons, due to the rocket engine."[craft](let [t(thrustcraft)]{:xt:y0:z0}))

The total force on a craft is just the sum of gravity and thrust. To sum these maps together, we’ll need a way to sum the x, y, and z components independently. Clojure’s merge-with function combines common fields in maps using a function, so this is surprisingly straightforward.

The acceleration of a craft, by Newton’s second law, is force divided by mass. This one’s a little trickier; given {:x 1 :y 2 :z 4} we want to apply a function–say, multiplication by a factor, to each number. Since maps are sequences of key/value pairs…

user=>(seq {:x1:y2:z3})([:z3][:y2][:x1])

… and we can build up new maps out of key/value pairs using into…

user=>(into {}[[:x4][:y5]]){:x4, :y5}

… we can write a function map-values which works like map, but affects the values of a map data structure.

And that allows us to define a scale function which scales a set of coordinates by some factor:

(defn scale"Multiplies a map of x, y, and z coordinates by the given factor."[factorcoordinates](map-values(partial * factor)coordinates))

What’s that partial thing? It’s a function which takes a function, and some arguments, and returns a new function. What does the new function do? It calls the original function, with the arguments passed to partial, followed by any arguments passed to the new function. In short, (partial * factor) returns a function that takes any number, and multiplies it by factor.

So to divide each component of the force vector by the mass of the craft:

Note that (/ m) returns 1/m. Our scale function can do double-duty as both multiplication and division.

With the acceleration and fuel consumption all figured out, we’re ready to apply those changes over time. We’ll write a function which takes the rocket at a particular time, and returns a version of it dt seconds later.

OK. Let’s save the rocket.clj file, load that code into the REPL, and fire it up.

user=>(use'scratch.rocket:reload)nil

use is like a shorthand for (:require ... :refer :all). We’re passing :reload because we want the REPL to re-read the file. Notice that in ns declarations, the namespace name scratch.rocket is unquoted–but when we call use or require at the repl, we quote the namespace name.

Great; there it is on the launchpad. Wow, even “standing still”, it’s moving at 463 meters/sec because of the earth’s rotation! That means you and I are flying through space at almost half a kilometer every second! Let’s step forward one second in time.

In evaluating this expression, Clojure reached a point where it could not continue, and aborted execution. We call this error an exception, and the process of aborting throwing the exception. Clojure backs up to the function which called the function that threw, then the function which called that function, and so on, all the way to the top-level expression. The REPL finally intercepts the exception, prints an error to the console, and stashes the exception object in a special variable *e.

In this case, we know that the exception in question was a NullPointerException, which occurs when a function received nil unexpectedly. This one came from clojure.lang.Numbers.ops, which suggests some sort of math was involved. Let’s use pst to find out where it came from.

This is called a stack trace: the stack is the context of the program at each function call. It traces the path the computer took in evaluating the expression, from the bottom to the top. At the bottom is the REPL, and Clojure compiler. Our code begins at user/eval1478–that’s the compiler’s name for the expression we just typed. That function called scratch.rocket/step, which in turn called Numbers.add, and that called Numbers.ops. Let’s start by looking at the last function we wrote before calling into Clojure’s standard library: the step function, in rocket.clj, on line 125.

123(assoc craft124; Time advances by dt seconds125:t(+ dt(:tcraft))

Ah; we named the time field :time earlier, not :t. Let’s replace :t with :time, save the file, and reload.

Look at that! Our position is unchanged (because our velocity was zero), but our velocity has shifted. We’re now moving… wait, -9.8 meters per second south? That can’t be right. Gravity points down, not sideways. Something must be wrong with our spherical coordinate system. Let’s write a test in test/scratch/rocket_test.clj to explore.

Oh, wait, that looks odd. {:x 1 :y 1 :z 0} is on the equator: phi–the angle from the pole–should be pi/2 or ~1.57, and theta–the angle around the equator–should be pi/4 or 0.78. Those coordinates are reversed! Double-checking our formulas with Wolfram MathWorld shows that we mixed up phi and theta! Let’s redefine cartesian->polar correctly.

This time, our velocity is increasing in the +x direction, at half a meter per second. We have liftoff!

Flight

We have a function that can move the rocket forward by one small step of time, but we’d like to understand the rocket’s trajectory as a whole; to see all positions it will take. We’ll use iterate to construct a lazy, infinite sequence of rocket states, each one constructed by stepping forward from the last.

Notice that each map is like a frame of a movie, playing at one frame per second. We can make the simulation more or less accurate by raising or lowering the framerate–adjusting the parameter fed to trajectory. For now, though, we’ll stick with one-second intervals.

How high above the surface is the rocket?

(defn altitude"The height above the surface of the equator, in meters."[craft](-> craft:positioncartesian->spherical:r(- earth-equatorial-radius)))

Now we can explore the rocket’s path as a series of altitudes over time:

The million dollar question, though, is whether the rocket breaks orbit.

(defn above-ground?"Is the craft at or above the surface?"[craft](<= 0(altitudecraft)))(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(* 1003600)]; 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 (flighttrajectory))))(defn apoapsis"The highest altitude achieved during a trajectory."[trajectory](apply max (map altitudetrajectory)))(defn apoapsis-time"The time of apoapsis"[trajectory](:time(apply max-key altitude(flighttrajectory))))

If the rocket goes below ground, we know it crashed. If the rocket stays in orbit, the trajectory will never end. That makes it a bit tricky to tell whether the rocket is in a stable orbit or not, because we can’t ask about every element, or the last element, of an infinite sequence: it’ll take infinite time to evaluate. Instead, we’ll assume that the rocket should crash within the first, say, 100 hours; if it makes it past that point, we’ll assume it made orbit successfully. With these functions in hand, we’ll write a test in test/scratch/rocket_test.clj to see whether or not the launch is successful:

When we exhaust the fuel reserves of the primary stage, we’ll de-couple the main booster from the Centaur. In terms of our simulation, the Atlas V will be replaced by its next stage, the Centaur. We’ll write a function stage which separates the vehicles when ready:

We’re using cond to handle three distinct cases: where there’s fuel remaining in the craft, where there is no stage to separate, and when we’re ready for stage separation. Separation is easy: we simply return the next stage of the current craft, with the current craft’s time, position, and velocity merged in.

Finally, we’ll have to update our step function to take into account the possibility of stage separation.

Still crashed–but we increased our apoapsis from 750 kilometers to 4,598 kilometers. That’s plenty high, but we’re still not making orbit. Why? Because we’re going straight up, then straight back down. To orbit, we need to move sideways, around the earth.

Orbital insertion

Our spacecraft is shooting upwards, but in order to remain in orbit around the earth, it has to execute a second burn: an orbital injection maneuver. That injection maneuver is also called a circularization burn because it turns the orbit from an ascending parabola into a circle (or something roughly like it). We don’t need to be precise about circularization–any trajectory that doesn’t hit the planet will suffice. All we have to do is burn towards the horizon, once we get high enough.

To do that, we’ll need to enhance the rocket’s control software. We assumed that the rocket would always thrust in the +x direction; but now we’ll need to thrust in multiple directions. We’ll break up the engine force into two parts: thrust (how hard the rocket motor pushes) and orientation (which determines the direction the rocket is pointing.)

(defn unit-vector"Scales coordinates to magnitude 1."[coordinates](scale(/ (magnitudecoordinates))coordinates))(defn engine-force"The force vector, each component in Newtons, due to the rocket engine."[craft](scale(thrustcraft)(unit-vector(orientationcraft))))

We’re taking the orientation of the craft–some coordinates–and scaling it to be of length one with unit-vector. Then we’re scaling the orientation vector by the thrust, returning a thrust vector.

As we go back and redefine parts of the program, you might see an error like

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)

This is a stack trace from the Clojure compiler. It indicates that in scratch/rocket.clj, on line 69, column 11, we used the symbol unit-vector–but it didn’t have a meaning at that point in the program. Perhaps unit-vector is defined below that line. There are two ways to solve this.

Organize your functions so that the simple ones come first, and those that depend on them come later. Read this way, namespaces tell a story, progressing from smaller to bigger, more complex problems.

Sometimes, ordering functions this way is impossible, or would put related ideas too far apart. In this case you can (declare unit-vector) near the top of the namespace. This tells Clojure that unit-vector isn’t defined yet, but it’ll come later.

Now that we’ve broken up engine-force into thrust and orientation, we have to control the thrust properly for our two burns. We’ll start by defining the times for the initial ascent and circularization burn, expressed as vectors of start and end times, in seconds.

(def ascent"The start and end times for the ascent burn."[03000])(def circularization"The start and end times for the circularization burn."[40001000])

Now we’ll change the thrust by adjusting the rate of fuel consumption. Instead of burning at maximum until running out of fuel, we’ll execute two distinct burns.

We’re using cond to express four distinct possibilities: that we’ve run out of fuel, executing either of the two burns, or resting with the engines shut down. Because the comparison function <= takes any number of arguments and asserts that they occur in order, expressing intervals like “the time is between the first and last times in the ascent” is easy.

Finally, we need to determine the direction to burn in. This one’s gonna require some tricky linear algebra. You don’t need to worry about the specifics here–the goal is to find out what direction the rocket should burn to fly towards the horizon, in a circle around the planet. We’re doing that by taking the rocket’s velocity vector, and flattening out the velocity towards or away from the planet. All that’s left is the direction the rocket is flying around the earth.

With the mathematical underpinnings ready, we’ll define the orientation control software:

(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)(:timecraft)(last ascent))(:positioncraft); 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)(:timecraft)(last circularization))(rejection(:velocitycraft)(:positioncraft)); Otherwise, just point straight ahead.:else(:velocitycraft)))

For the ascent burn, we’ll push straight away from the planet. For circularization, we use the rejection function to find the part of the velocity around the planet, and thrust in that direction. By default, we’ll leave the rocket pointing in the direction of travel.

With these changes made, the rocket should execute two burns. Let’s re-run the tests and see.

Review

(ns scratch.rocket);; Linear algebra for {:x 1 :y 2 :z 3} coordinate vectors.(defn map-values"Applies f to every value in the map m."[fm](into {}(map (fn [pair][(key pair)(f(val pair))])m)))(defn magnitude"What's the radius of a given set of cartesian coordinates?"[c]; By the Pythagorean theorem...(Math/sqrt(+ (Math/pow(:xc)2)(Math/pow(:yc)2)(Math/pow(:zc)2))))(defn scale"Multiplies a map of x, y, and z coordinates by the given factor."[factorcoordinates](map-values(partial * factor)coordinates))(defn unit-vector"Scales coordinates to magnitude 1."[coordinates](scale(/ (magnitudecoordinates))coordinates))(defn dot-product"Finds the inner product of two x, y, z coordinate maps. See http://en.wikipedia.org/wiki/Dot_product"[c1c2](+ (* (:xc1)(:xc2))(* (:yc1)(:yc2))(* (:zc1)(:zc2))))(defn projection"The component of coordinate map a in the direction of coordinate map b. See http://en.wikipedia.org/wiki/Vector_projection."[ab](let [b(unit-vectorb)](scale(dot-productab)b)))(defn rejection"The component of coordinate map a *not* in the direction of coordinate map b."[ab](let [a'(projectionab)]{:x(- (:xa)(:xa')):y(- (:ya)(:ya')):z(- (:za)(:za'))}));; Coordinate conversion(defn cartesian->spherical"Converts a map of Cartesian coordinates :x, :y, and :z to spherical coordinates :r, :theta, and :phi."[c](let [r(magnitudec)]{:rr:phi(Math/acos(/ (:zc)r)):theta(Math/atan(/ (:yc)(:xc)))}))(defn spherical->cartesian"Converts spherical to Cartesian coordinates."[c]{:x(* (:rc)(Math/cos(:thetac))(Math/sin(:phic))):y(* (:rc)(Math/sin(:thetac))(Math/sin(:phic))):z(* (:rc)(Math/cos(:phic)))});; The earth(def earth-equatorial-radius"Radius of the earth, in meters"6378137)(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."(/ (* 2Math/PIearth-equatorial-radius)earth-day))(def g"Acceleration of gravity in meters/s^2"-9.8)(def initial-space-center"The initial position and velocity of the launch facility"{:time0:position{:xearth-equatorial-radius:y0:z0}:velocity{:x0:yearth-equatorial-speed:z0}});; Craft(defn centaur"The upper rocket stage. http://en.wikipedia.org/wiki/Centaur_(rocket_stage) http://www.astronautix.com/stages/cenaurde.htm"[]{:dry-mass2361:fuel-mass13897:isp4354:max-fuel-rate(/ 13897470)})(defn atlas-v"The full launch vehicle. http://en.wikipedia.org/wiki/Atlas_V"[next-stage]{:dry-mass50050:fuel-mass284450:isp3050:max-fuel-rate(/ 284450253):next-stagenext-stage});; Flight control(def ascent"The start and end times for the ascent burn."[0300])(def circularization"The start and end times for the circularization burn."[4001000])(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)(:timecraft)(last ascent))(:positioncraft); 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)(:timecraft)(last circularization))(rejection(:velocitycraft)(:positioncraft)); Otherwise, just point straight ahead.:else(:velocitycraft)))(defn fuel-rate"How fast is fuel, in kilograms/second, consumed by the craft?"[craft](cond; Out of fuel(<= (:fuel-masscraft)0)0; Ascent burn(<= (first ascent)(:timecraft)(last ascent))(:max-fuel-ratecraft); Circularization burn(<= (first circularization)(:timecraft)(last circularization))(:max-fuel-ratecraft); Shut down engines otherwise:else0))(defn stage"When fuel reserves are exhausted, separate stages. Otherwise, return craft unchanged."[craft](cond; Still fuel left(pos? (:fuel-masscraft))craft; No remaining stages(nil? (:next-stagecraft))craft; Stage!:else(merge (:next-stagecraft)(select-keys craft[:time:position:velocity]))));; Dynamics(defn thrust"How much force, in newtons, does the craft's rocket engines exert?"[craft](* (fuel-ratecraft)(:ispcraft)))(defn mass"The total mass of a craft."[craft](+ (:dry-masscraft)(:fuel-masscraft)))(defn gravity-force"The force vector, each component in Newtons, due to gravity."[craft]; Since force is mass times acceleration...(let [total-force(* g(masscraft))](-> craft; Now we'll take the craft's position:position; in spherical coordinates,cartesian->spherical; replace the radius with the gravitational force...(assoc :rtotal-force); and transform back to Cartesian-landspherical->cartesian)))(declare altitude)(defn engine-force"The force vector, each component in Newtons, due to the rocket engine."[craft]; Debugging; useful for working through trajectories in detail.; (println craft); (println "t " (:time craft) "alt" (altitude craft) "thrust" (thrust craft)); (println "fuel" (:fuel-mass craft)); (println "vel " (:velocity craft)); (println "ori " (unit-vector (orientation craft)))(scale(thrustcraft)(unit-vector(orientationcraft))))(defn total-force"Total force on a craft."[craft](merge-with + (engine-forcecraft)(gravity-forcecraft)))(defn acceleration"Total acceleration of a craft."[craft](let [m(masscraft)](scale(/ m)(total-forcecraft))))(defn step[craftdt](let [craft(stagecraft)](assoc craft; Time advances by dt seconds:time(+ dt(:timecraft)); We burn some fuel:fuel-mass(- (:fuel-masscraft)(* dt(fuel-ratecraft))); Our position changes based on our velocity:position(merge-with + (:positioncraft)(scaledt(:velocitycraft))); And our velocity changes based on our acceleration:velocity(merge-with + (:velocitycraft)(scaledt(accelerationcraft))))));; Launch and flight(defn prepare"Prepares a craft for launch from an equatorial space center."[craft](merge craftinitial-space-center))(defn trajectory[dtcraft]"Returns all future states of the craft, at dt-second intervals."(iterate #(step%1)craft));; Analyzing trajectories(defn altitude"The height above the surface of the equator, in meters."[craft](-> craft:positioncartesian->spherical:r(- earth-equatorial-radius)))(defn above-ground?"Is the craft at or above the surface?"[craft](<= 0(altitudecraft)))(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 10 hours are up?"[trajectory](let [time-limit(* 103600)]; 10 hours(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 (flighttrajectory))))(defn apoapsis"The highest altitude achieved during a trajectory."[trajectory](apply max (map altitude(flighttrajectory))))(defn apoapsis-time"The time of apoapsis"[trajectory](:time(apply max-key altitude(flighttrajectory))))

As written here, our first non-trivial program tells a story–though a different one than the process of exploration and refinement that brought the rocket to orbit. It builds from small, abstract ideas: linear algebra and coordinates; physical constants describing the universe for the simulation; and the basic outline of the spacecraft. Then we define the software controlling the rocket; the times for the burns, how much to thrust, in what direction, and when to separate stages. Using those control functions, we build a physics engine including gravity and thrust forces, and use Newton’s second law to build a basic Euler Method solver. Finally, we analyze the trajectories the solver produces to answer key questions: how high, how long, and did it explode?

We used Clojure’s immutable data structures–mostly maps–to represent the state of the universe, and defined pure functions to interpret those states and construct new ones. Using iterate, we projected a single state forward into an infinite timeline of the future–evaluated as demanded by the analysis functions. Though we pay a performance penalty, immutable data structures, pure functions, and lazy evaluation make simulating complex systems easier to reason about.

Had we written this simulation in a different language, different techniques might have come into play. In Java, C++, or Ruby, we would have defined a hierarchy of datatypes called classes, each one representing a small piece of state. We might define a Craft type, and created subtypes Atlas and Centaur. We’d create a Coordinate type, subdivided into Cartesian and Spherical, and so on. The types add complexity and rigidity, but also prevent mis-spellings, and can prevent us from interpreting, say, coordinates as craft or vice-versa.

To move the system forward in a language emphasizing mutable data structures, we would have updated the time and coordinates of a single craft in-place. This introduces additional complexity, because many of the changes we made depended on the current values of the craft. To ensure the correct ordering of calculations, we’d scatter temporary variables and explicit copies throughout the code, ensuring that functions didn’t see inconsistent pictures of the craft state. The mutable approach would likely be faster, but would still demand some copying of data, and sacrifice clarity.

More imperative languages place less emphasis on laziness, and make it harder to express ideas like map and take. We might have simulated the trajectory for some fixed time, constructing a history of all the intermediate results we needed, then analyzed it by moving explicitly from slot to slot in that history, checking if the craft had crashed, and so on.

Across all these languages, though, some ideas remain the same. We solve big problems by breaking them up into smaller ones. We use data structures to represent the state of the system, and functions to alter that state. Comments and docstrings clarify the story of the code, making it readable to others. Tests ensure the software is correct, and allow us to work piecewise towards a solution.

Exercises

We know the spacecraft reached orbit, but we have no idea what that orbit looks like. Since the trajectory is infinite in length, we can’t ask about the entire history using max–but we know that all orbits have a high and low point. At the highest point, the difference between successive altitudes changes from increasing to decreasing, and at the lowest point, the difference between successive altitudes changes from decreasing to increasing. Using this technique, refine the apoapsis function to find the highest point using that inflection in altitudes–and write a corresponding periapsis function that finds the lowest point in the orbit. Display both periapsis and apoapsis in the test suite.

We assumed the force of gravity resulted in a constant 9.8 meter/second/second acceleration towards the earth, but in the real world, gravity falls off with the inverse square law. Using the mass of the earth, mass of the spacecraft, and Newton’s constant, refine the gravitational force used in this simulation to take Newton’s law into account. How does this affect the apoapsis?

We ignored the atmosphere, which exerts drag on the craft as it moves through the air. Write a basic air-density function which falls off with altitude. Make some educated guesses as to how much drag a real rocket experiences, and assume that the drag force is proportional to the square of the rocket’s velocity. Can your rocket still reach orbit?

Notice that the periapsis and apoapsis of the rocket are different. By executing the circularization burn carefully, can you make them the same–achieving a perfectly circular orbit? One way to do this is to pick an orbital altitude and velocity of a known satellite–say, the International Space Station–and write the control software to match that velocity at that altitude.

This is a nice intro to functional techniques, but there are hundreds of these around the web. What’s really missing is more and better examples of dealing with resource management and side effects in FP.

Just a question, do you have anything against records and protocols? As far I understand this is be suitable for model data, as presented here. It would be safer than using only maps. Maybe, you consider this approach is too OO? Or you just want to keep it simple?

not sure if this would be possible using records (fresh Clojure noob here). So probably you prefer maps for flexibility? Is there a way to combine this with records or protocols to make it safer, or would you stick with this in a serious / bigger app?

Protocols are strictly for polymorphism, and since we don’t have different types of records here, there’s no need for a protocol. We’ll be addressing polymorphism in a later chapter, though. :)

Records have two advantages over maps: first, lookup for the basis fields in a record is an order of magnitude faster than a hashmap, and second, records can specify implementations of polymorphic functions (namely, via protocols and interfaces). They don’t provide type safety for their fields (you can happily assign a Cat to a Dog field) and don’t constrain their fields (you can assoc arbitrary keys into a record, just like a map), so they don’t really provide the sort of type safety you’d be looking for in a typed Object or algebraic datatype.

Because Records print differently and have a few subtle API differences, and because I’m trying to give learners a chance to settle in to the core data abstractions before moving up to polymorphism and types, this chapter sticks to maps. :)

ascent and circularization have typos in the narrative which prevent the circularization stage from ever being triggered - 3000 instead of 300 and 4000 instead of 400

I would also recommend promising the orientation function later while you’re talking about unit-vector and engine-force under Orbital insertion. While I was working on this I thought you had just forgotten to include it, and tried to work it out myself.

Thanks again!

Post a Comment

Please avoid writing anything here unless you are a computer:
Captcha
This is also a trap:
Comment

Name

Email

Http

Body

Supports github-flavored markdown for
[links](http://foo.com/), *emphasis*, _underline_, `code`, and >
blockquotes. Use ```clj on its own line to start a Clojure code block,
and ``` to end the block.