# Computing an Energy-Maneuverability (E-M) Diagram for an Aircraft

Computing an Energy-Maneuverability (E-M) Diagram for an Aircraft

The Energy-Maneuverability (E-M) Diagram is a helpful way to show the horizontal turning performance of an aircraft at various speeds.

While I have found many discussions of these diagrams online, I have not found anything explaining how to compute a theoretical diagram for a specific aircraft. Here is my attempt to create one for a WWII fighter aircraft - the FM2 Wildcat:

QUESTION

Am doing this right?

To simplify discussion, I have listed the Variables in Appendix 1 and the Equations in Appendix 2.

MAXIMUM ANGLE OF ATTACK (AOA) LINE

The left blue line (currently labelled AoAmax and partly overwritten by the green line) shows the turn rates at maximum Angle of Attack (AOA) at various speeds.

Simply put, the Angle of Attack (AOA) produces Lift. The equation for Lift max is:

  Lift = cL*Q*S [eq. 2.0]


where cL is the coefficient of Lift, Q is the Dynamic Pressure and S is the Wing Area. For a given speed, the maximum cL produces the maximum Lift. The cL is a function of the Angle of Attack. The AOA is "the angle at which the chord of an aircraft's wing meets the relative wind. The chord is a straight line from the leading edge to the trailing edge". The maximum AOA is generally around 16 degrees. The cL is generally around 10% of the AOA.

At take-off speed, Lift is just greater than Weight, so the aircraft can fly. But the aircraft is unable to turn because the aircraft needs all of the lift to fly. As aircraft speed increases, lift at maximum AOA increases, which means that there is more lift available for turning.

To create this line, I started by computing the take-off speed (Vto). At this speed, Lift max = W (Weight) and the turn rate is zero. The equation for Lift max is:

  Lift max = cLmax*Q*S [eq. 2.1]


The equation for Q (Dynamic Pressure) is:

  Q = 1/2*ρ*V^2 [eq. 4.0]


We can combine the above into

  W = Lift max = (1/2*ρ*V^2)*(cLmax*S)


Solving for V gives:

  Vto = SQRT(W/(1/2*ρ*cLmax*S)) [5.1]
Vto = SQRT(7048/(1/2*0.00237*1.6*260)) = 119.57 ft/sec = 81.53 mph


You can then compute the AOA max speed for each different G-Force (GF). At each different GF:

  Lift = GF*W


Since everything else remains constant, you can modify the Vto equation to:

  Vgf = SQRT((GF*W)/(.5*ρ*cLmax*S))


Since all of the variables remain unchanged, except for GF, you can simplify this to:

  Vgf = Vto*SQRT(GF)


Thus, where GFmax = 7:

  Vgf = 119.57*SQRT(7) = 316.36 ft/sec = 215.69 mph.


You might notice that this is the same equation used to compute speed for accelerated stalls. Thus, for each GF, the Vmin value on the E-M Diagram equals the minimum speed for an accelerated stall.

Since the E-M diagram computes horizontal turn rate and since the maximum horizontal turn is in level flight, you can use the standard equations to compute the turn rate at each Vgf:

  BankAngle = DEGREES(ACOS(1/GF)) [eq. 6.0]


Thus where GF = 7:

  BankAngle = DEGREES(ACOS(1/7)) = 81.79 deg
TurnRate = DEGREES(TAN(RADIANS(81.79)))*32.174/316.36 = 40.39 deg/sec


This means that the shape of the Vmin curve will be about the same for all aircraft. For a given aircraft, the size of the curve will vary depending on S (wing area) and standard W (weight). For a given flight, the size will vary depending on actual W and p (air density).

The curve will extend from Vto (where TurnRate is zero) to the speed and turn rate at maximum GF.

G-FORCES (GF) LINES

The red lines show the G-Forces at various speeds. You can use the equations above to plot these lines.

In level flight, each G-Force (GF) has a corresponding BankAngle, which remains constant for the entire speed range. You can compute this BankAngle using the BankAngle equation (above):

  BankAngle = DEGREES(ACOS(1/GF)) [eq. 6.1]


You can then use TurnRate equation (above) to compute the TurnRate at various speeds:

  TurnRate = DEGREES(TAN(RADIANS(BankAngle)))*g/V [eq. 7.0]


Ideally, the top GF Line will show the upper boundaries of what the aircraft and pilot can tolerate. For modern aircraft, this line is generally at 9 g's, which is what the pilot can tolerate on a sustained basis. Those aircraft can generally operate near this upper limit in level flight.

I am using a lower GF Line for the FM2 because the aircraft simply does not have the power to reach this limit in a sustained turn and rarely (if ever) in an instantaneous turn. (Anecdotally, a test pilot who dove the plane straight down was never able to pull more than 7.2 g's.)

VNE LINE

Normally, as noted below, there should a vertical line at Vne (the Never-Exceed Speed). However, the Wildcat was one of the extremely rare aircraft for which there was no real Vne. You could apparently dive the aircraft all the way to terminal speed (around 550 mph).

SUSTAINED TURN (PS) CURVES

The green line (labeled PS = 0) shows the maximum possible sustained TurnRates at various speeds in level flight.

The first three lines: the Max AOA Line, the Max GF Line and the VNE Line should set the boundaries within which the aircraft can operate. The Sustained Turn Lines show the TurnRates available on a sustained basis within those boundaries.

The maximum TurnRate is determined by a combination of Excess AOA and Excess Net Thrust above what is required for straight and level flight. Where the Line matches the Max AOA Line, the TurnRate is limited by the AOA Available. There is still Excess Net Thrust Available, which can be used to accelerate. At the peak of the Line, where the TurnRate is greatest, the aircraft is using all of the available AOA and Net Thrust. At higher speeds, the TurnRate is limited by Net Thrust Available. The reason is that, as speed increases, Thrust decreases (on a prop aircraft) and Dp increases to the point that the excess of Thrust over Dp is not enough to offset the Di at the maximum BankAngle for that speed. This means that BankAngle and TurnRate must decrease so that Thrust = Total Drag (Dt). At the top speed all available Thrust is required to offset Dt in level flight. There is no excess Thrust to handle the increase in Di associated with a turn.

Since the max TurnRate is computed assuming we are using all excess Thrust to increase TurnRate, there is no excess Thrust remaining to increase speed. Thus, to increase speed, you must decrease TurnRate.

The computation of the maximum TurnRate for different speeds involves several equations and takes several steps (or, as below, you could combine them all into a single equation).

(1) Compute the maximum Di available using:

  Di= Thrust max - Dp
where:
Thrust max = 550*BHPmax*η/V [eq. 1.0]
Dp = cD0*Q*S [eq. 3.1]


(2) Compute the maximum cDi using:

  cDi = Di/(Q*S) [eq. 3.2.1]


(3) Compute the related cL using:

  cL = SQRT(cDi*PI()*e*AR) [eq. 2.1]


If the cL is greater than cLmax, then cL = cL max and there is no reduction in max TurnRate.

(4) Compute the maximum BankAngle at that speed using:

  BankAngle=DEGREES(ACOS(W/Lift)) [eq. 6.1]
where:
Lift = cL*Q*S [eq. 2.0]


(5) Compute the TurnRate for that speed using:

  TurnRate = DEGREES(TAN(RADIANS(BankAngle)))*g/V [7.0]


You can also compute the PS curves for situations where you climb or descend. You can descend to use gravity to supplement Thrust and increase TurnRate. The resulting curve line will be higher than the PS0 curve. Conversely, climbing will decrease TurnRate. The resulting curve line will be lower than the PS0 curve. However, to keep things simple, I am just computing the PS0 curve.

Several posts below discuss how to compute the PS curve in greater detail. As those discussions indicate, those computations are made using various assumptions, some of which are necessary to simplify the computation. This is an area where one would expect an E-M diagram based on computations to vary from an E-M diagram based on actual flight data.

CONTOUR MAPPING (added Jan 17, 2023)

Several of the responses below refer to contour graphing, which is a new concept to me. Although, I am familiar with contour maps, but never considered how those contour lines were created. For others who are in the same boat, here is a simple explanation of what that involves.

In the question above, I am plotting lines created using specific formulae. In contrast, contour mapping involves creating arrays of data and then calling a contour-plotting formula which will parse the data, perform some interpolations and plot the data.

For example, the MatLab program creates several different "arrays" of data indexed by reference to increments of Speed and GLoad. For example, the PS value for Speed=300fps and GLoad=5 is stored in a PS array at location PS[50Speed_Increment,5GLoad_Increment]. The program then calls a built-in contour graphing function in which you specify the values that you want plotted (e.g. PS=0). The program will then parse the data to estimate where PS would equal 0. It correlates GLoad to TurnRate and plots the Speed and TurnRate values on the chart.

My research indicates that Excel contour charts are unable to draw the kind of contour charts shown below. One major limitation is that Excel apparently cannot create a single chart that uses more than one data source, e.g. GForce, cLmax and PS.

APPENDIX 1. VARIABLES

AR = Aspect Ratio of Wing = 5.55
cL = coefficient of Lift; cLmax = 1.6
cD = coefficient of D = sum of cD0 and cDi
cD0 = coefficient of Dp = 0.021
cDi = coefficient of Di (equations below)
D = Drag
Dp = Parasitic Drag (equation below)
Di = Induced Drag (equation below)
e = Wing Efficiency = 0.75
BHP = Brake Horsepower; BHPmax = 1300
g = Gravity = 32.174 ft/sec^2
GF = G-Force
η = Propeller Efficiency = 0.75
ρ = Air Density at Sea Level = 0.002377 slug/ft^3
Q = Dynamic Pressure (equation below)
S = Wing Area = 260 ft^2
V = Velocity or Speed (ft/sec)
Vgf = Speed for Given G-Force
Vmin = Minimum Speed for Given Turn Rate
Vto = Takeoff Speed (equation below)
W = Weight = 7048 lbs


APPENDIX 2. EQUATIONS

[1.0] Thrust = 550*BHP*η/V
Thrust max = 550*BHPmax*η/V
[2.0] Lift = cL*Q*S = (1/2*ρ*V^2)*(cL*S)
Lift max = cLmax*Q*S = (1/2*ρ*V^2)*(cL*S)
[2.1] cL = SQRT(cDi*PI()*e*AR) (from equation 3.3.1)
[3.0] Drag (D) = Dp+Di = cDt*Q*S
[3.1] Parasitic Drag (Dp) = cD0*Q*S
[3.2] Induced Drag (Di) = cDi*Q*S
[3.2.1] cDi = Di/(Q*S)
[3.3.1] cDi = cL^2/(PI()*e*AR)
[4.0] Dynamic Pressure (Q) = 1/2*ρ*V^2
[5.0] V = SQRT(W/(.5*ρ*cL*S))
Vto = SQRT(W/(.5*ρ*cLmax*S))
[5.2] Vgf = SQRT((GF*W)/(.5*ρ*cLmax*S))
[6.0] BankAngle = DEGREES(ACOS(1/GF))
[6.1] BankAngle = DEGREES(ACOS(W/Lift))

• One aspect of an EM chart is that you don't just want to plot Ps=0, but also other contours of Ps. This can be done algebraically (like you have done) for simple drag polars and engine models, but it is a real pain. Instead, it is normal to use a contour plot to establish these lines implicitly. I.e. calculate Ps for every point on a grid of turn rate and velocity, then contour that data. Excel's contour plots are terrible, you'll want to use Matlab or Python. Jan 9 at 16:39
• The interesting aspect here is you get an airframe EM curve readily, but it is power (thrust) constrained. The EM diagrams have lines where they say "you can do this but you are gaining or losing energy (speed or altitude)" based on available thrust and drag. (Diving adds thrust). Sustained turning would be no net energy loss in the manuver. Naturally, more drag from turning requires more thrust. No wonder best rate of turn is around Vy. Jan 10 at 7:55
• @Rob, Yes I can add some of those lines. I was just concerned that the diagram would get a bit messy, especially since the results seemed to simply parallel the Ps = 0 line. To compute those numbers for a particular vertical speed (VrtSpd), I determined the angle required to generate the required climb or descent (Angle=DEGREES(TAN(VrtSpd/Spd)) and then determined the portion of the Gravity force to add or subtract from Di (Force=W*SIN(RADIANS(Angle)). Jan 10 at 10:05
• @PhilCrowther well the key would be to have the lines for both planes. I recently saw an EM for a Sabre and Mig 15. The Sabre had automatic slats (like Me109), which helped its low speed ROT. The Mig climbed very well. Power/Weight will factor in. Jan 10 at 13:12
• @PhilCrowther I would argue that without different lines of Ps, it isn't really an EM chart. The lines of Ps are what let you compare two aircraft at a particular point and determine who has the advantage. The greater point is the power of using contour plots instead of algebraically solving for the curves. This is a must-do when the drag polar or engine model are more complex. It is also a lot easier. Jan 10 at 21:37

I call those diagrams turn rate diagrams, because that is plotted over airspeed. Never mind.

Am doing this right?

I think so, your result looks fine. Now let's start with the explanation how to compute a theoretical diagram for a specific aircraft which you are missing on the interwebs. Please allow me to use the standard symbols in my equations which differ from what you used.

First you draw lines of equal turn radii which radiate away from the origin. Next, you draw lines of equal load factor which are hyperbolic. This is what you get:

The equations for those lines are $$\text{turn rate}=\frac{\text{speed}}{\text{radius}}$$ and $$\text{turn rate}=\frac{g\cdot\sqrt{n_z^2-cos\gamma}}{\text{radius}}$$ $$n_z$$ is the load factor and $$\gamma$$ is the flight path angle, positive values of which denote a dive. The plot above is for horizontal flight, so $$cos\gamma$$ = 1. This is to account for the vertical component of lift.

You will notice that I used odd load factors below $$n_z$$ = 2. They are chosen for specific bank angles, namely $$\Phi$$ = 30° and 45°.

There are two ways to continue, one for instationary turns and one for stationary turns. Stationary means you can fly this without losing altitude, whereas instationary turns trade height for additional energy to maximize the turn rate. Of course, you could as well calculate a climbing turn where less energy is available for turning.

Instationary turns are easy: You plot two segments, one for minimum speed and one for maximum load factor. The turn rate at minimum speed is determined by the wing loading and the maximum lift coefficient $$c_{L_{max}}$$ and is $$\text{turn rate}=\sqrt{\frac{(n_z^2-cos\gamma)\cdot g\cdot c_{L_{max}}\cdot S\cdot\rho}{2\cdot m\cdot n_z}}$$ Expressed over speed, this becomes $$\text{turn rate}=\frac{\text{speed}}{R_{min}}=\frac{g\cdot\sqrt{\left(\frac{\rho\cdot\text{speed}^2\cdot S\cdot c_{L_{max}}}{2\cdot m\cdot g}\right)^2-1}}{\text{speed}}$$ with $$m$$ the airplane mass, $$S$$ the wing area and $$\rho$$ air density. When you vary $$n_z$$ (upper equation) or speed (lower equation), you get the line for maximum lift. If the term in the square root is negative, speed is lower than minimum speed at 1 g and no solution is possible.

The next limit is the maximum load factor the structure can withstand and applies above the minimum speed. The result should look like this:

where I added coloured lines for different engines with different thrust over speed behavior. A piston engine is represented by the blue line. Above the coloured lines the airplane cannot sustain the turn rate with engine power alone while below it the airplane needs to throttle the engine in order to not climb away. In other words: The coloured lines are the limits for stationary turns.

You will note that this is not for the Wildcat but some generic plane. Also, $$c_{L_{max}}$$ is usually not constant but goes up with load factor because of the increasing Reynolds number. To do this correctly for a specific airplane will need to include all speed-dependent adjustments.

If your airplane is limited in angle of attack (like the Phantom II was), replace $$c_{L_{max}}$$ by $$c_{L\alpha}\cdot\alpha_{max}$$.

• This is an excellent explanation and expansion of my question. To assist readers, I have tried to change my variables to match yours. I did not change GF to n because I thought there might be some confusion with the Greek "n" used to represent prop efficiency. Jan 10 at 21:33
• @PhilCrowther The Greek $\eta$ looks distinct from an n, having that longer right leg. Jan 10 at 22:40

Generally, these diagrams are produced from extensive flight test and data collection. You cannot get sufficient accuracy from the simplified equations you are using. Computational Fluid Dynamics applications software use linear approximations of Navier Stokes equations and mathematical models of the complete airframe and wing to perform an iterative (millisecond by millesecond) computation, which still needs to be validated by flight test data.

That is not to say that the approach you are using cannot generate a qualitatively accurate picture. I think it does. I would just be hesitant to put any actual values on the graph that imply accuracy that it does not have.

1. There should be a vertical line at Vne.
2. The left boundary is not a minimum speed line. It is a maximum Angle of Attack (AOA) line. A 1G line would reflect level flight. The 1 G line, where it intersects with the far left (Max AOA) line would designate 1 G stall speed, or minimum Takeoff speed. The Max AOA line would continue down and left of that point to the zero G line.
3. It is important to understand that Specific power is about sustainability, not instantaneously capability. The zero Ps line should not be aligned at the left (AOA) boundary at very high Gs. Can FM-2 Wildcat really sustain (continuous) 3 Gs at Maximum AOA (Clmax)? The zero Ps line should intersect the far left edge of the envelope at the point you would end up at if you select full power and pull on the pole to just below stall AOA.
4. In your text, you say 

The green line shows the Specific Power (PS) curve for full power in level flight (PS0). This shows the maximum possible TurnRate at various speeds

This is not correct. First, it shows the zero Specific excess power line at each airspeed, specifying the maximum sustainable G-load at that airspeed. It does not show level flight. That is simply the 1G line. Secondly, it is not the maximum possible turn rate at that airspeed. It is the maximum sustainable turn rate at that airspeed. I emphasize this because this is the essence of the entire concept of Energy maneuverability. It is about sustainability, not about instantaneous capability.

• "The left boundary is not a minimum speed line. It is a maximum Angle of Attack (AOA) line" only when your airplane is AoA limited, like the Phantom II. Normally it is indeed a minimum speed line at the respective load factor. Jan 9 at 15:36
• @Peter, All aircraft are AOA limited. It's called Cl max or stall Angle of Attack. Minimum speed at at some specific load factor is just another way of looking at it, except that using that perspective does not take into consideration gross weight. (The less you weigh, the higher the allowable load factor for the same stress on the airframe structure). But Max AOA is the same, regardless of gross weight, bank angle, or airspeed. Jan 9 at 22:08
• @Phil, I found your statement about Vne surprising. Do you know what it's terminal velocity was when pointed straight down with full throttle? Do you know- was this design constraint you mention accidental or deliberate? Jan 9 at 22:21
• "But Max AOA is the same, regardless of gross weight, bank angle, or airspeed." Right only when you are AoA limited, but wrong for most airplanes. In all airplanes the maximum lift coefficient and maximum AoA go up with speed due to Reynolds number effects. This becomes even more pronounced for high pitch rates. Jan 10 at 5:32
• @PhilCrowther This answer shows a large collection of windtunnel data for maximum lift coefficient over Re. Above 2 million you can see a fairly steady increase over log(Re) which should also be valid for the Wildcat wing. Jan 10 at 12:39

Here is a sample EM diagram for an F-16.

Note that we do not expect all of the Ps lines to intersect the CLmax limit.

It looks to me like you may have a problem with your induced drag term not being sensitive to increased CL required for the turns.

• CDi as a function of CL is OK -- but the CL needs to be CL=nW/(qS) taking into account n - the load factor during the turn. Jan 11 at 22:43
• In climbing & turning flight, L=nWcos(theta). You appear to be using a sin instead of cos. When computing Ps, we usually ignore the cos(theta) in the formula for L (and thereby CL). With it, the calculation becomes iterative -- the climb angle depends on the drag, drag depends on the lift, lift depends on the climb angle. While this seems contradictory when you derive Ps from forces, it makes more sense when you work from an energy perspective. A pilot can choose to put the excess power into climb or acceleration (or some of both). The angle is not determined. Jan 12 at 6:15
• I agree that the 1g line should correspond to zero turn rate. Not sure what is going on here. This figure was taken from a F-16 SAC chart -- prepared by GD and delivered to the USAF. It is certainly a reputable source. Jan 12 at 6:26
• I checked several other SAC charts -- for the F16 and other aircraft. This one was the only one with an EM diagram. The n=1 contour is absolutely wrong. Good eye. All I can figure is some idiot manager thought the figure looked 'bad' without a n=1 line and told some poor draftsperson to insert one. Somehow it made it out the door. Jan 12 at 7:21
• When generating performance charts (flight envelopes, skymaps, doghouse charts, etc.) (i.e. when computing Ps) we assume theta=0. Ps is not simply rate of climb. It is the rate of change in energy height (ze=h+0.5*V^2/g). You can put Ps into altitude or velocity (or both). At Ps=0, you can trade altitude for speed. At Ps<0, you can dive or slow down. Jan 12 at 16:39

I wrote a simple turn code in Matlab.

Here is the EM / Doghouse plot at sea level.

It looks like it agrees with your chart. The reason we don't see the Ps curves 'close' inside the CLmax line is insufficient thrust. If you increase the engine power, you should see the center of the 'eye' move to the right.

Moving to a non-constant propeller efficiency model would also be a nice update.

I've calculated Ps to the left of the stall line -- but that should be taken with a heap of salt. Not only will the wing not provide enough lift to turn at those points, there would be tremendous drag that is not accounted for in the simple parabolic polar.

I think you will see how using a contour algorithm to draw the lines greatly simplifies the code. All of the equations are in their canonical 'forward' mode. In addition to eliminating all the algebra required to put things into other forms, this approach also enables more complex models that can not be algebraically re-arranged -- such as a complex drag polar or thrust model.

Yes there are a lot of lines -- the chart can be improved by labeling the contours less frequently and less conspicuously. In addition, use of different line weights, styles, and shades of grey / colors can help a lot.

In general, I find that complex plots like this are still usable/readable as long as the various lines are crossing one another. If they run parallel, then the different trends become very confusing to discern.

clear all
format compact

Sref = 260;
W = 7000;

nmax = 7;

% Maximum speed to be considered.
Vmax = 500; % ft/s

etap = 0.75;
Pmax = 1300 * 550; % lbf ft / s

CLmax = 1.6;

% Simple parabolic drag polar.
AR = 5.5;
CD0 = 0.021;
e = 0.75;

% Only perform calculations at sea level.
rho = 0.0023769;

g = 32.1740; % ft/s^2

% Resolution to perform calculations, increase for smoother lines.
nn = 50;
nV = 51;

ns = linspace( 1, nmax, nn );
Vs = linspace( 1, Vmax, nV );

for in = 1:nn
n = ns( in );

for iV = 1:nV
V = Vs( iV );

q = 0.5 * rho * V^2;

L = n * W;

CL = L / ( q * Sref );

CD = CD0 + ( CL^2 ) / (pi * e * AR );

D = CD * q * Sref;

T = Pmax * etap / V;

Ps = V * ( T - D ) / W;

omega = g * sqrt( n^2 - 1.0 ) / V;

nmat( in, iV ) = n;
Vmat( in, iV ) = V;
Psmat( in, iV ) = Ps;
CLmat( in, iV ) = CL;
omegamat( in, iV ) = omega * 180.0 / pi;
end
end

figure( 1 )

[c,h] = contour( Vmat, omegamat, nmat, 1:1:nmax, 'k' );

axis( [0 Vmax 0 50] )

clabel( c, h );
hold on

[c,h] = contour( Vmat, omegamat, CLmat, CLmax * [1 1], 'k' );
clabel( c, h );

[c,h] = contour( Vmat, omegamat, Psmat, 0:50:1000, 'k' );
clabel( c, h );
hold off

xlabel( 'V (ft/s)' )
ylabel( '\omega (deg/s)' )

• Since contour mapping is new to me, I have added a section to my question where I try to explain what it means. Let me know if my explanation is correct. Jan 17 at 10:37
• @PhilCrowther I think your understanding of contour plotting is pretty much correct. I've created a simple contour plot in Excel, but it is very awkward and I haven't found much capability to adjust things. There probably is a way, but I would recommend using a tool with better plotting capability (Matlab, Octave, Python/Matplotlib, Gnuplot, TecPlot, etc.) Jan 18 at 0:07

Here is a simple contour plot example done in Excel... It is a form of a 3D surface plot -- where your viewpoint is from directly above.

The red column ends up being the X-data. The blue row ends up being the Y-data.

These can probably be swapped/transposed somehow.

The green data is a 'matrix'. Each cell contains a computation corresponding to the X,Y of the matching row,column.

In this case, the equation is:

=(F$$1-0.5)^2+($$A6-1)^2-(0.5)^2
`

I.e. f(x,y)=0 will be the equation of a circle centered at 0.5,1 with radius 0.5.

The fact that Excel draws the grid makes it particularly easy to understand how contouring works. Look at an individual square with a contour line drawn through it like this one...

Each of the corners has a value associated with it. This is the value from the green data in the matrix.

In this case, we are interested in drawing a contour at zero -- but we could draw any other value. We start by searching for any edge that has endpoint values that differ in sign -- i.e. one is greater than zero and the other is less than zero.

Whenever we find such an edge, we perform interpolation to create a point along that edge that estimates where the zero value would be. Then for each square, we connect those interpolated points together.

We keep doing this for the entire matrix and the result is a contour plot.

A Mathematician might call this an implicit function.