Table of Contents
- Introduction
- Radar Cross Section
- Incremental Length Diffraction Coefficients
- RCS of the F-117
- Timeline of the F-117 Development
- References
- ILDC Plots
Introduction
The development leading up to the F-117 aircraft in 1981 was heavily dependent on software simulation of its radar visibility. Different stories have been told about the methods used in the ECHO computer program by the famed Skunk Works division of the US aerospace manufacturer Lockheed. The work by Ufimtsev, published in the Soviet Union, was an important part. Using the theoretical methods available at that time, a computer program was built to calculate the radar return of the F-117 in a similar way to the initial ECHO program.
A |\ ^
/|\ | \ |
/<^>\ |\V| 2
/ / \ \ | |/ 0
/|/\ /\|\ TOP |/ | LEFT m
/ /\ | /\ \ | / |
\/ /|\ \/ |A \ |
V V \_\ _
|<- 13 m -->|
_
\_-^-_/ FRONT 4m
___/_/|\_\___ _
Lockheed F-117 Nighthawk
Radar Cross Section
To quantify the visibility of an object on radar the Radar Cross section (RCS) 𝜎 in units of m² is used. It depends on the radar frequency, polarization, location of antenna and the orientation, material and shape of the target object. But it is independent of the distance from the radar to the object (far field approximation).
Usually the send and receive antennae are assumed to be at the same location and their polarization to be the same. This is the monostatic (backscatter) case, in contrast to the general bistatic.
Due to the large range involved it is usually reported on a logarithmic scale as dBm² with 10 log₁₀ 𝜎.
incident -> target
( ) ) ) ) ( ) ( <====//
__|__ <- reflection
/radar\_____________________________
Incremental Length Diffraction Coefficients
The method of the Incremental Length Diffraction Coefficients (ILDC) described by Mitzner (1974a) in his report is a generalization of Ufimtsev’s theory. Thus we will use it to retrace the methods which could have been employed at the time of the development. The ILDC was developed at Northrop, the competitor, to the winner Lockheed with its F-117. Thus it was most certainly not used for the ECHO software. But having the most general approach allows specialization easily, while the converse does not. Also the ILDC report by Mitzner explicitly treats the handling of the apparent singularities, as we will see later. While e.g. Knott et al. (2004, p. 214) does not give their formulas.
The report was written in 1974 and was publicly released in 1976. The existence of the result, but no details where already known in 1974. (Knott 1974)
Physical Optics and Edges
The Physical Optics (PO) approximation allows the calculation of the high frequency scattering of a body. It only takes the radar reflection of flat surfaces into account.
Thus the first improvement is made by also considering edge effects. For the 2D case of an edge the Physical Theory of Diffraction (PTD) can be used to increase the accuracy of the result (Ufimtsev 1971). For the 3D case, it remains restricted to the Keller cone. This is unfortunate for our use case, since thus the edge effect will be restricted to one specific direction and be zero everywhere else. For the monostatic case we only get a return when the edge is perpendicular to the incident (and thus scattering) direction. (David C. Aronstein 1997, p. 218)
The ILDC is an generalization of the PTD coefficients to arbitrary directions and thus ideal for 3D calculations. (Mitzner 1974a)
Implementation
Common Lisp was chosen as the programming language. (it already has exp(i x) as cis function and supports complex numbers.) The code follows closely the report. The implementation details and references to the equations of the report can be found in the source code cl-rcs.tar.bz2.
No external dependencies are needed besides ASDF and the fiveam testing library for running the tests.
Some plotting functions with output to ASCII where developed and are used for the visualizations. (Getting nice axis ticks is not so easy due to a Moiré interference of the tick position and the low resolution character position.)
Data Structures
Some linear algebra functions for supporting 3D vectors were needed. The ILDC formulas further require dyads. To combine it with the PO calculation a special mesh structure is advantageous.
The dyadic diffraction coefficient is implemented as a 3x3 (complex) matrix.
While the vector is a struct with the three components x, y and z (More
on dyads in Lindell (1992, p. 17)).
Some of the values can be complex (polarization vector & diffraction
coefficient), thus some of the functions need to support this.
The mesh of flat polygons is stored as a half-edge data structure. This allows access to all the edges of each face to calculate the PO contribution, as well as the two faces for each edge for its (PTD/ILDC) effect.
Multiple polygons
The solution given in the report omit the phase factor (Mitzner 1974a, p. 81). To calculate the return of several polygons their returns have to be summed, taking into account their respective phases (eq. 7.1 Knott et al. 2004). Thus this phase factor needs to be added to the solution.
Singularities and Numerical Problems
The main problem encountered is related to numerical instability due to floating point numbers. Unfortunately this is not treated in the report. The use of double precision reduces the problem somewhat (Gordon 1988). Furthermore the calculation of angles from vectors was optimized and a pairwise sum over the individual facets is used.
Apparent Singularities
The treatment of the apparent singularities in the calculation is more complicated than apparent from the report, since comparing floating point number one needs to consider rounding errors.
For example the special case f(cos 𝜓, 𝜓) in Mitzner (eq. 3-71 1974a)
naively reads as cos 𝜓 = V but for implementation we use a floating point
equal: (float-equal (cos V-big) psi epsilon) where epsilon is e.g. 1e-3.
The epsilon is defined for each case based on inspecting the plot of the
function near the point under consideration. This is of course not perfect,
since the behavior of the general case depends on all variables and thus might
approach the singularity more rapidly than anticipated by our choice of epsilon.
For comparison: The Matlab code fun_fg.m in Ufimtsev (2014, p. 219)
uses a epsilon of 1e-5.
There is an error in the last term eq. 3.64A: the function sin[(𝛽ₛ + 𝛽ᵢ)/2] appearing therein must be replaced by sin[(𝛽ₛ + 𝛽ᵢ)/2]. (Knott et al. 2004, p. 215)
Furthermore the special case for g(cos 𝜓, 𝜓) in eq. 3-73 should not have the 1/𝛎 factor. This was numerically discovered and verified. It can also be seen by comparing the special case for the function g eq. 3-73 to same special case for f eq. 3-71 and their general cases eq. 3-66 and 3-65. For our case 𝜓 = v we see that the only difference is the extra factor of -1/sin v in f. Thus the special case for g should be like f but without the factor. (The difference in sign apparently disappears; but I was unable to derive the formulas)
The numerical evaluation for the relevant part of g can be seen in the plots below. It is clear, that the first one is wrong. In the second plot we see, that the function is still a bit discontinuous at the boundary from the special case to the general one. This is because for our choice of parameters the chosen epsilon is to large. That is, we could extend the general case a bit more.
g(cos 𝜓, 𝜓) as defined in 3-73
2.0 .|
1.5 .|
|'''''''' ''''''''
1.0 '|
0.50 -|
0.0 .| '''''''''''''''
|
+------------------------------
| | |
0.490 0.500 0.510
g(cos 𝜓, 𝜓) not divided by the wedge factor ν
1.5 .| '..
| ''..
1.5 .| '
| ...............
1.4 .|
| ..
1.4 .| '':..
+------------------------------
| | |
0.490 0.500 0.510
The sin x / x parts of the formulas (4-23) and (4-29) are implemented as sinc x to remove the singularity at x = 0.
Grazing incidence
When the incidence and scattering direction are close to grazing, the PTD and ILDC have singularities and thus can not predict the physical result. (Mitzner 1974a, p. 59; Gordon 1988, p. 2)
This is the result of not taking into account traveling waves. Besides using the modified theory it is not clear what can be done. (Ufimtsev 2014, p. 245)
The problem is exacerbated when triangulating a curved body. A sphere will thus always have triangles near grazing incidence (around the apparent edge), which lead to a wrong result. Returning zero for near grazing incidence improves the result.
Verification
The correctness of the code is check with unit tests using the fiveam library. The first verification is done using the symmetries and special cases listed in the report.
Further verification is done using reports of measured and calculated cross sections of test objects. The report of Gordon (1988) is especially useful, since it contains the calculation of only the edge contributions. Unfortunately the used reports only have the data as plots and not as table. WebPlotDigitizer Lite was used to extract data points.
The plots Figure 7 and 8 (Diamond) in the report of Mitzner (1974a, p. 89) could not be reproduced. Figure 7 has an apparently wrong RCS scale with the maximum at 5 dB; but our calculation has a maximum of 10 dB consistent with the broad side return and Figure 8. Furthermore the peaks of the two Figures do not align with our calculation. Since the trapezoid plots (Figures 9-11) match visually and our calculation also match data from other sources, there is probably some error in the parameters for the Figures 7 and 8. (see: 7)
Features & Limitations
The code supports bistatic (untested) and monostatic RCS calculation using PO, PTD (on the Keller cone) and ILDC.
Object meshes can be loaded from Wavefront .obj files. Self-Shadowing and multiple diffraction are not supported; thus limiting the objects to convex bodies. Thus it is important that the mesh is a manifold (though a opening could be used as e.g. a fully absorbing air intake).
The surfaces are assumed to be perfect electric conductors (PEC), with no support for radar absorbent coatings, further limiting its usefulness for stealth design and evaluation.
The execution speed is not optimized and uses only a single thread, despite, in principle, being easily able to exploit parallelism. In practice this is not a concern due to the processing power of contemporary computers.
Comparison with ECHO
The initial version 1 of ECHO only used PO and assumed metallic (PEC) surfaces. An analytic solution by Bill Schroeder to the PO integral for a triangle was used. (David C. Aronstein 1997, p. 19)
This can be calculated in the same way in the implemented program. Surprisingly the solution by Schroeder was only for triangles, even though this was probably just a special case of the solution for polygons. Gordon (e.g. 1975)
Edge diffraction
There was supposedly initially also some, albeit inaccurate, form of edge contribution in ECHO:
The ECHO 1 solution also included contributions from edges, because the edges form the boundaries of the surface integral. So the designers could use ECHO 1 to manage edge orientations.
(David C. Aronstein 1997, p. 19)
While later in the same text it is claimed to have only used PO (David C. Aronstein 1997, p. 218). The PO integral and its solution is along the border of the surface, i.e.: the edges of the triangle. This is probably what is meant by the above quote and not an additional effect.
This can be seen in the following RCS diagram of a quadratic plate using PO at an incidence angle of 20° to the normal. The peaks are perpendicular to the edge direction and thus a form of edge contribution.
.
. .
.::' '::.
:.' '.:
' : : + : : .
:'. .':
'::. .::'
' '
'
r: [-130.00 50.00]
The following quote in the same text can be read to support this:
Schroeder’s solution provided information on the strength of the surface returns even for non specular angles when the surface in question is not perpendicular to the incident radiation.
(David C. Aronstein 1997, p. 19)
So this is actually a comparison of PO to Geometric Optics, which only gives returns for specular reflection. Other sources note that there was initially no edge contribution and it was later added using PTD. (Westwick 2010a, p. 47).
The initial version was quite accurate since a special edge treatment reduced their radar return and thus contribution to the RCS (David C. Aronstein 1997, p. 20; Westwick 2010a, p. 48).
Later, for the development of the F-117, edge diffraction based on the PTD and empirical corrections for radar absorbing materials were added (David C. Aronstein 1997, p. 20; Westwick 2010a, p. 65). This might be the modified PTD in Knott et al. (eq. 5.44 2004) extended to scattering outside the Keller cone. The difference to the ILDC tends to disappear as more edges are added (Gordon 1988). Thus result of the ILDC will probably be similar to the version of of the PTD used in ECHO. (The PTD algorithm also available in the program uses no extension and thus restricts the edge diffraction to the Keller cone.)
Since radar absorbing materials are not taken into account, our values will be higher that those of ECHO.
RCS of the F-117
Disregarding any radar absorbent coating we calculate the monostatic RCS of the F-117 based on a simplified model from K987 (2018). The air inlet and cockpit are assumed to be fully reflective. (Westwick 2010a, p. 53; Crickmore 2014, p. 27)
RCS Plots
The RCS of the F-117 in dBm² for 2 GHz is plotted in the horizontal plane and the vertical symmetry plane. With the front of the aircraft is up and left respectively; as shown with the ASCII art in the center.
Using just the PO we get a nice butterfly RCS plot (Westwick 2010a, p. 79).
The return in the important frontal and back threat sector is greatly reduced at the expense of the sides (Knott et al. 2004, p. 277). The requirement was for a 90° wide front and back threat sector. The frequency range was 200 MHz to 10 GHz. Lockheed focused on the high frequencies. (Westwick 2010a, p. 79)
F-117 monostatic RCS 2 GHz HH using PO
.. ..
. '' : ' ' : '' .
. ' : ..' . .. .. . '.. : '
' ' ' ' ::.:':.:: ' ' ' '
. ' ' ' ' ''/^\'' ' ' ' ' .
. /\|/\ .
'.' ' /^\ ' '.
. ' .. .: :. .. ' .
' ' ... '' '.'.' '' ... ' '
' :.. ''' :'''.''': ''' ..: '
. .. ' : ' ' : ' .. .
'' :.' '.: ''
' '
F-117 monostatic RCS 2 GHz HH using PO
'
.. '
. .
. ':. .
' .:' '::':
' .' ' .
.'. :' :.
:: :' ::'
' . .::
':':'. :::
':._/=>_//::
::: :':.
::'' '::
.::' '':
':' :
.' ':.:: .
:'. :' ' .
. .
. '
' :
Using the ILDC the butterfly is less noticeable and there is more variance.
F-117 RCS 2 GHz HH using ILDC
' . . .
' . .' ' . . . '
. : '.'. ...'... .. : . '
. .:.: .::'..'' .: :'.. ..''. .
: : ' ' . '' ' . . '. . .
. :' . ' ' .
. '. ' ' ' : ' :
: . ' ' /^\ ' :
. /\|/\ .
': . /^\ .'.
. ' ' ' ' ' ' .
' ' . ' .' . ' ''' '
'''' . .. ' ' ' '
. ..' : '.''.. : . . :. .
. '.' . '::.. : .''. ..'. ' '
'.' . '''.' ''' .'
' . . : ' '' '
. .
' '
F-117 RSC 2 GHz HH using ILDC
'
'. '
. .
. '.. . .
' .'' '::''
'' ': ' . '
.'. :.. :
.': ' '::
'. ::.:. '
':':' .:.:: .
:'.:::'' _/=>_// ::.''. .
. ''''::: ' .: :: '
. :::. .'::
' '':..' '':
'' ' '. ..:
.: ::: :
::. .' ' .
. . '
. '
' :
Comparing some selected directions in dBm² with ± 2° samples reported as mean and standard deviation:
| method | front | back | side | bottom | top |
|---|---|---|---|---|---|
| PO | -43 ± 3 | -36 ± 5 | 5 ± 4 | 32 ± 16 | 26 ± 12 |
| ILDC | -17 ± 9 | -14 ± 2 | 4 ± 7 | 32 ± 16 | 26 ± 12 |
Calculation using vertical polarization leads to similar results.
As expected the ILDC values are higher, since now also the edges contribute to the radar return. The model is not fully convex. For example, when viewed from the front, the two vertical stabilizers on the back are partially obstructed. But for the calculation their full return is added, thus overestimating the RCS.
Frequency Response
F-117 frequency response frontal aspect HH in dbm²
10. '| :
| .
0.0 -| '
| ::'..
-10. -| ' ....... .
| ' ' .'.'' . .'... . . .
-20. -| ''.'.''.: ' :..:. . .
| ''' ' '' :': '
-30. .| ' . .
+-------------------------------------
| | |
0.00 5.000e+09 1.000e+10
As can be seen in the frequency response plot (Hz vs dBm²), the stealthiness is less for lower frequencies.
Expected RCS
We can check the bottom RCS easily, since it is mostly the broad side return of
flat triangle about 16 m by 13 m
(eq 5.17 Knott et al. 2004).
For 2 GHz we get
68 dBm²
while using the ILDC gives
67 dBm².
Reported values are -16 dBm² Dough (2001, p. 57) or the size of a golf ball (also an eagles eye) for a prototype (i.e. 0.042 m; RCS = 𝜋 r²), about -22 dBm² by Rich and Janos (1994, p. 36). These will probably be for the thread sector (i.e. front or back). The values, due to being logarithmic, cover a wide range. Thus we can only say that our value is not too wrong.
Calculation Speed
The fastest computer in 1969 could do about 3 MFLOPS while my student laptop can
do about 328 MFLOPS on a
single core.
The calculation was reported to take over night (Westwick 2010a, p. 48).
Decks of punched data cards were initially used (David C. Aronstein 1997, p. 18).
The calculation for a full RCS pattern in 1° steps on my laptop
takes about 142 seconds.
The left-right symmetry might have reduced the calculation by half; but having
both horizontal and vertical polarization will have increased it again.
Comparing with overnight (maybe 8 h =
28800 s),
the estimated computing power used at Lockheed was
1.6 MFLOPS.
The actual speed was probably quite a bit lower, since my program is not at all
optimized for speed and the model has more facets than the initially used models.
Timeline of the F-117 Development
- 1950s PO theory developed
- 1962 Publication of Ufimtsev (1962) in the Soviet Union
- 1966-1970 Development of first RCS computer program at Northrop
- 1971 Translation of Ufimtsev’s work (Ufimtsev 1971)
- 1970s explicit solution for the PO from a polygonal flat plate e.g. Gordon (1975) or Mitzner (1974b) as reported in Mitzner (1974a, pp. 3, 86)
- 1974 DARPA program for RCS reduction to Northrop and McDonnell Douglas
- 1975 Lockheed joins program, Development of ECHO 1 based on PO, Hopeless Diamond
- 1976 DARPA awards Phase II to Lockheed, Building of the Have Blue prototype
- 1976-1978 Added edge diffraction to ECHO 2 using the PTD and non-metallic materials
- 1981 First flight of the F-117
- 1990 First public showing of the F-117
(David C. Aronstein 1997, p. 269; Westwick 2010b, p. 71)
References
ILDC Plots
Recreation of the plots in Mitzner (Figure 7-11 1974a).
Figure 7
10. '| '''..
5.0 '| ':.
0.0 -| ':
-5.0 .| :
-10. .| :
| .
-15. '| ' ....
-20. '| ' :' :.
-25. -| . : '. ..
-30. .| . ' .'''':.
-35. .| ' ' : '. ....
| ' ' . '. .'
+--------------------------------------------------
| | | |
0.00 10.0 20.0 30.0
Figure 8
10. '| '':.
5.0 '| '.
0.0 -| :
-5.0 .| :
-10. .| :
| .
-15. '| . ....
-20. '| . .' ':
-25. -| .' :
-30. .| ' . ' .:'''.
-35. .| ' ' : ..
| '. ' ' : .' ''.
+--------------------------------------------------
| | | | |
0.00 5.00 10.0 15.0 20.0
Figure 9
20. '| .
15. '| ''
10. -|
5.0 .| ' .
0.0 .| . .
| .'
-5.0 '| ..' ''.
-10. '| .........''' '.
-15. -| ''''''' ''.
-20. .| ''.
-25. .| '..
| '..
+--------------------------------------------------
| | |
-50.0 0.00 50.0
Figure 10
20. '| .
15. '| ''
10. -|
5.0 .| . '
0.0 .| . .
| ':.
-5.0 '| :' ':..
-10. '| : ''':........
-15. -| :' ''''''''
-20. .| .:'
-25. .| .:'
| .:'
+--------------------------------------------------
| | |
-50.0 0.00 50.0
Figure 11
20. '| .
15. '| ..
10. -|
5.0 .|
0.0 .| : '
| '
-5.0 '|
-10. '| . .
-15. -|
-20. .| .'' ''
-25. .| . '''
| .. ...
+--------------------------------------------------
| | |
-50.0 0.00 50.0