Navigation Through Fog: A Mathematical Analysis of 

Folded Path Propagation

By 

David H. Tofsted 

Computational and Information Sciences Directorate 

Battlefield Environment Division

ARL-TR-2648 September 2002 

Approved for public release; distribution unlimited. 

Army Research Laboratory 
NOTICES 

Disclaimers 

The findings in this report are not to be construed as an official 

Department of the Army position, unless so designated by other 

authorized documents. 

Citation of manufacturers’ or trade names does not constitute an official 

endorsement or approval of the use thereof. 
Army Research Laboratory 

White Sands Missile Range, New Mexico 88002-5513 

ARL-TR-2648 August 2002

Navigation Through Fog: A 

Mathematical Analysis of Folded 

Path Propagation 

David H. Tofsted 

Computational and Information Sciences Directorate 

Approved for public release; distribution unlimited. 
Preface

The problem of navigation through fog was brought tragically to the fore during

the summer of 1999 when John F. Kennedy Jr.’s plane crashed into the Atlantic off

the coast of Martha’s Vinyard. The key factor here was the inability to determine

the vertical direction due to the lack of visual cues. Less catastrophic, but of more

common concern, are frequent delays caused by fog at America’s major airports.

It is estimated that one major air package carrier alone incurs yearly costs measured in millions of dollars due to delayed or rerouted flights caused by airport

fog conditions. Such delays would not be necessary if alternatives to the current

methods of handling fog conditions were available.

The current method of coping appears simply to involve reducing the rate of aircraft arrivals and departures such that collisions are avoided. However, such an

approach does not avoid the problems encountered by John Kennedy Jr., nor does

it curtail the possibility of collision with obstacles on the runway. (Incidences

of wildlife appearing on runways under low visibility conditions have also been

noted in the past.)

Though several methods for improving the ability to see through fog have been

attempted, the one of interest in this study involves the use of stereoscopic vision

devices coupled with laser illumination and realtime image deblurring software.

But a fundamental question remains as to the capabilities of this technology when

confronting real world fog environments. While models of fog scattering properties exist, the propagation models for use of these models have often relied on

Monte Carlo analysis of scattering by fog particles. This analysis is usually incomplete — it cannot adequately address the spatial and angular structure of the

radiation that reaches objects within the fog volume and the return signals from

objects in the fog. In particular, if corner reflectors were mounted at the sides of

the runways, could such reflective objects be viewed through the fog? To answer

this question, the current analysis provides the physical model of one-way propagation through the fog to a retroreflector embedded in the fog, and a return path

back to camera receivers mounted on each of the aircraft wings.

i
ii
Contents

Preface i

Executive Summary v

1 Introduction 1

2 Forward Propagation Through Fog 3

2.1 The Radiative Transfer Equation .................. 3

2.1.1 Terms within the Radiative Transfer Equation ....... 4

2.1.2 Solid Angle ......................... 4

2.2 The Small Angle Approximation .................. 5

2.3 Phase Function Representation ................... 7

2.4 Small Angle RTE . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.5 Fourier Transform Operation . . . . . . . . . . . . . . . . . . . . 12

2.6 Gaussian Beam Model . . . . . . . . . . . . . . . . . . . . . . . 15

2.6.1 Complex Beam Wave . . . . . . . . . . . . . . . . . . . 15

2.6.2 Radiance Form for a Diverging Gaussian Beam . . . . . . 16

2.7 Radiative Transfer Solution . . . . . . . . . . . . . . . . . . . . . 18

2.7.1 Implicit Solution Method . . . . . . . . . . . . . . . . . . 19

2.7.2 Fourier Solutions . . . . . . . . . . . . . . . . . . . . . . 22

2.7.3 Scattered Solution . . . . . . . . . . . . . . . . . . . . . 24

2.7.4 Error Function Application . . . . . . . . . . . . . . . . . 25

2.7.5 The Inverse Transform . . . . . . . . . . . . . . . . . . . 26

2.8 Numerical Integration . . . . . . . . . . . . . . . . . . . . . . . . 28

3 Retroreflector Interactions 33

4 Return Path Propagation 38

iii
5 Diffuse Radiation 41

5.1 Diffuse Radiation Scale Transform . . . . . . . . . . . . . . . . . 41

5.2 Dual Direction Radiance Equations . . . . . . . . . . . . . . . . . 43

6 Conclusions 46

References 47

SF 298 49

Figures

1 Angular coordinate system using two angular measures, θ and φ ... 6

2 Geometric relationship between angle θ and tangent ω ........ 7

3 Advection fog phase functions for 1.5 μm and 1.54 μm wavelengths . 9

4 Forward hemisphere portion of phase function for advection fog at

1.5 μm wavelength . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

5 Backward hemisphere portion of phase function for advection fog at

1.5 μm wavelength . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

6 Gaussian beam transverse irradiance relationship . . . . . . . . . . . 17

7 Gaussian beam transverse spread (w) and curvature (R) as functions

of range (z) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

8 Angular relations between various vectors in transform and real coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

9 Orientation of vector U relative to component vectors −κr and −η

−Zω . . 28

10 Division of the quarter plane in p and q, using regions delineated by

constant surfaces in the cylindrical coordinates ρ and θ . . . . . . . . 31

11 Set of gridded mini-reflectors as they appear on reflective material . . 33

12 Geometry of a single mini-reflector . . . . . . . . . . . . . . . . . . . 34

13 Effective geometry of a retroreflector considering the fraction of incident energy that returns in the retro-direction . . . . . . . . . . . . . 35

14 Geometry of the folded path problem where we treat the return path

as a reversed version of the forward propagation path . . . . . . . . . 38

Tables

1 Gaussian expansion coefficients for forward and backward hemisphere

phase function representations of the 1.5 μm advection fog phase

function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

iv
Executive Summary

This report discusses a mathematical analysis for assessing the capabilities of a

system to navigate through fog. This study was conducted as part of an FY99-

funded U.S. Army Research Laboratory (ARL) Director’s Research Initiative. The

main project objective was to determine the optimal means of implementing a

system designed to enhance the capability of seeing objects through low visibility

atmospheres. From a military standpoint, all weather forces will never become a

reality unless aircraft can navigate safely to landing under adverse ground visibility conditions. Often ground conditions are the limiting concern due to low-lying

haze or fog. In the majority of these cases, it is not the act of taking off that is of

concern but that of landing under reduced visibility. Additionally, from a civilian

standpoint, reduced visibility conditions can often hinder pilots from landing at

major metropolitan airports. On a yearly basis millions of dollars of business are

lost simply due to aircraft rerouting, delays in landing due to reduced visibility

conditions, or both.

A recent patent application by Wendell Watkins of ARL holds the potential of

solving this perennial problem through the use of active stereo illumination and

observation of a landing field. This system envisions use of a near-infrared laser

mounted on an aircraft using pulsed beams to illuminate reflectors outlining the

runway. A portion of the laser light is presumed to penetrate the fog and be reflected from retroreflective corner cubes located along the sides of the runway.

Cameras mounted on each wing of the aircraft are then used to acquire an image

of these returning signals. The received signal is passed through deblurring software and presented to the pilot in a headmounted display using stereo goggles.

The basic physical model of this system involves propagating a Gaussian laser

beam forward through a fog layer to a retroreflector. The reflected beam then

passes through the fog a second time to a receiver camera. We study the characteristics of the return path beam in terms of positional and directional variations.

The mathematical simulation of the propagation process utilizes the small angle

approximation when dealing with the largely forward-scattering fog aerosols studied. This small angle approximation results in a simpler version of the equation of

radiative transfer being used in the propagation analysis.

The fog aerosol used was obtained from the U.S. Army’s Electro-Optical Systems

Atmospheric Effects Library Phase Function Database (Tofsted et al., 1997).1

Analysis of data in this database provided a summed Gaussian form for the forward and backward scattering portions of the phase function. The phase function

1Tofsted, D. H., B. T. Davis, A. E. Wetmore, J. Fitzgerrel, R. C. Shirkey, and R. A. Sutherland,

EOSAEL92, Aerosol Phase Function Data Base PFNDAT, ARL-TR-273-9, U.S. Army Research

Laboratory, White Sands Missile Range, NM, 1997.

v
forms were then then coupled to a Gaussian beamwave model of the laser source

beam. The effects of the retroreflector material were modeled as a spatial clipping coupled with a nondiffraction-limited spreading (angular convolution) and

reflection of the incident beam. The return path equations were then assessed using similar equations used to assess the outbound beam. The combination of outbound propagation, retroreflector interaction, and return beam propagation then

allow for overall analysis of the propagation problem.

Approach

To model this propagation problem we used a series of related approximations

as well as processing techniques made possible by the approximations. The first

step involved invocation of the so-called small angle approximation. This step

simplifies the radiative transfer equation (RTE) and permits the scattering aerosol

to be modeled as varying only along the z axis. This approximation is also related

to the use of a sum of Gaussian terms expressing the angular structure of the

scattering phase function for the aerosol studied. Using the small angle RTE, a

fourth-order Fourier transformation in the off-axis spatial position vector [r =

(x, y)] and a two-dimensional small-angle vector [ω ] is performed that allows

a general theory (c.f. Smirnov, 1964)2 to be used in solving for the propagated

radiance as a function of position and direction at a retroreflector plane.

This solution is then used in conjunction with an equation for the form of a Gaussian beamwave in this four-dimensional system. The net results of this analysis

can then be transformed back from the Fourier domain using numerical techniques

to study typical scenarios.

Conclusion

The methods described here are also useful for addressing other propagation modeling problems and questions and should prove useful in answering questions

about the capabilities of the proposed navigation through fog device.

2Smirnov, V. I., A Course of Higher Mathematics, Volume II, Pergamon Press/Addison-Wesley,

Reading, MA, 1964.

vi
1. Introduction

The performance of any electro-optical imaging system depends on the medium

through which signals are received. As optical systems are continually being upgraded for special purposes, it is reasonable to augment system design by studying

the impact of environment on system performance. One means of studying this

problem is to simulate the system its environment such that system-controlling

parameters can be modified and tested against a variety of atmospheric conditions

one can reasonably expect to encounter. In the case studied here, we specifically

want to look at optically dense fog conditions which would hinder the normal

performance of most sensor systems. This study was conducted as part of an

FY99-funded U.S. Army Research Laboratory (ARL) Director’s Research Initiative (DRI). The DRI program provides a means of funding breakthrough technologies at a startup level within ARL.

This report discusses a portion of this overall program involving a mathematical

analysis of the capabilities of said system. The main project objective was to determine the optimal means of implementing such a system for landing aircraft (to

include helicopters) under low-visibility ground conditions. From a military standpoint, all weather forces will never become a reality unless aircraft can safely land

under adverse visibility conditions such as ground fog or heavy haze.

To facilitate this advanced capability, Wendell Watkins of ARL has submitted a

patent (Watkins et al., 2000) to allow vision through fog via the use of active stereo

illumination and observation of a landing field. This system envisions use of a

near-infrared (IR) laser mounted on an aircraft fuselage that uses pulsed beams to

illuminate reflectors outlining a runway. A portion of the laser light is presumed to

penetrate the fog and be reflected. Cameras mounted on each wing of the aircraft

are then used to acquire an image of these returning signals. The received signal is

passed through deblurring software and presented to the pilot in a head-mounted

display using stereo goggles. The type of deblurring kernel used and the effectiveness of the deblurring software will in some respects depend on the variability of

the fog aerosols encountered, but may be anticipated knowing the general characteristics of the fog and its effects on propagation of the Gaussian laser source

beam through the fog.

In this report we study the characteristics of the return path beam in terms of positional and directional variations. The mathematical simulation of the propagation

process utilizes the small angle approximation since the fog aerosols studied are

largely forward scattering. The small angle approximation means that the majority

of the energy is scattering in a small angular cone around the initial propagation

direction. Use of this approximation results in a simpler version of the equation of

1
radiative transfer being used in the propagation analysis. The fog aerosol used was

obtained from the Army’s Electro-Optical Systems Atmospheric Effects Library

(EOSAEL) (Shirkey et al., 1987) Phase Function Database (PFNDAT) (Tofsted

et al., 1997). Analysis of data in this database provided a sum of Gaussian terms

forming the forward-scattering portion of the phase function.

The specific propagation problem considered here is the degree to which the

backscattering from the fog can be mitigated by looking at the return signal from

a source laser propagating from one wing of the aircraft if the camera viewing the

return is located on the opposite wing. Such a separation will do two things: (1)

Use of stereo goggles will allow an observer to better distinguish features in the

vision field; (2) By locating the radiation laser source on the one wing and the

receiver on the opposite wing it will tend to reduce the amount of backscatter. A

third advantage results from image processing of the received signals from both

wing cameras. There will be a differential change in the amount of backscatter

that spans both received images, but this backscatter signal will have an opposite

gradient on each (the image from the left wing will see more backscatter to the

right side of the image and visa-versa for the opposite wing camera images). There

could be advantages in differential processing of the images obtained from each

wing. However, to determine the extent of these effects, one would need to model

the propagation problem. This model is outlined in this text.

2
2. Forward Propagation Through Fog

In this chapter we consider the problem of mathematically modeling the single

path propagation of a Gaussian laser source beam through a fog aerosol to a

retroreflector-type object. Following the solution of this problem, chapter 3 will

deal with the interaction of this propagated beam with a retroreflective surface

embedded in the fog. Chapter 4 will then consider the propagation of the reflected

beam back through the fog to an off-axis camera. Finally, chapter 5 considers direct beam energy scattering at large angles - the diffuse radiation field. This latter

calculation is critical in evaluating the relative complicating effects of the diffusely

scattered energy in masking objects of interest in the field of view within the fog.

2.1 The Radiative Transfer Equation

To model the propagation of laser light through a fog layer the small angle approximation can be utilized. This technique allows for the solution of the radiative

transfer equation (RTE) using Fourier methods, permitting the calculation of both

the directional and positional structure of the forward propagating radiance field

at a plane. This approximation assumes the majority of the energy involved in the

radiative transfer process is only flowing into the forward hemisphere.

The small angle approximation involves a modification to the radiative transfer

equation from a fully three-dimensional form to a single hemisphere form. But to

discuss this change we must first consider the fully three-dimensional form of the

radiative transfer equation, given by,

Ωˆ ·∇I(s, Ω)+ ˆ σ I(s, Ω) = ˆ σs



4π

I(s, Ωˆ) P(Ωˆ, Ωˆ) dΩˆ+(σ−σs)b(T(s), λ). (1)

The primary quantity generated by this equation is a radiance field, I(s, Ω) ˆ . This

field varies as a function of position, denoted by the three-dimensional vector s,

and direction, denoted by the unit vector Ωˆ. Ωˆ is composed of directional cosines

oriented along the x, y, and z coordinate axes.

Radiance is measured in units of W/(m2-sr), which specifies the rate of energy

passage through a differentially small surface into a differentially small solid angle. In this analysis the spectral dependence of the radiation is suppressed, but in

general, a series of monochromatic radiative transfer calculations must be spectrally integrated to obtain results for band-averaging sensors. For laser radiation,

the waveband over which the laser operates is normally so narrow, and the spectrally varying characteristics of aerosols are normally so slowly varying with

wavelength, that the radiative transfer of such energy can be considered relatively

monochromatic.

3
2.1.1 Terms within the Radiative Transfer Equation

The first term on the left hand side (LHS) of the radiative transfer equation (RTE)

represents the rate of change in radiance due to propagation in direction Ωˆ. This

term is thus a directional derivative, dI/ds, where ds is the magnitude of ds, a

differential change in position in the direction Ωˆ.

The second term on the LHS represents the rate of extinction of the radiance field

due to both absorption and scattering, where σ is called the extinction coefficient

and is given in units of inverse distance. Relating this term with the first term

on the LHS, the product σ ds is thus the probability of interaction between the

radiance and particles embedded in the medium within the incremental path length

ds. For fog aerosols, σ ranges between 5 km−1 and 30 km−1, implying several

collisions (virtually all involving scattering phenomena) per photon per kilometer

travelled.

Due to radiation scattering, energy lost from the radiance field due to extinction

will most likely reappear as scattered energy via the first term on the right hand

side (RHS) of the RTE. This returned energy re-enters the propagating radiance

field by scattering from different directions into the path of the original photon.

The nature of this modified directional field is quantified in terms of the scattering phase function (P(Ωˆ, Ωˆ

)), which describes the probabilistic distribution of

scattering based on both the incident (Ωˆ

) and scattered (Ωˆ) directions.

The rate of scattering is proportional to the scattering coefficient, σs, given once

again in inverse distance units. σs is related to the extinction coefficient σ through

a quantity called the single scattering albedo, , which quantifies the probability

of a scattering event (as opposed to an absorption event) given that an interaction

between a photon and a particle occurs. Hence, we normally write,

σs =  σ. (2)

The final term on the RHS of the RTE represents emission of thermal radiation,

which depends on the radiation wavelength (λ), the temperature of the atmosphere

at that point (T(s)), and the absorption coefficient, given as,

σa = (1 − ) σ. (3)

The function b(T, λ) represents the radiation emitted by a perfect blackbody (implying an assumption of thermodynamic equilibrium). The absorption coefficient

modifies the blackbody emission to account for variations in emissivity associated with the properties of the aerosols and gaseous constituents of the medium.

However, we will ignore emission effects in our analysis becauseemission influences are insignificant at near-IR wavelengths for normal ranges of atmospheric

temperatures.

2.1.2 Solid Angle

The integration contained within the first term on the RHS of the RTE is required to evaluate the energy scattered into the output direction of interest, Ωˆ.

This integration is performed with respect to a variable known as solid angle.

4
Solid angle is used to evaluate how radiant energy spreads out following passage

through a point, or conversely, how radiant energy converges toward a point as it

propagates. The meaning of solid angle, and similarly, how one measures solid

angle, can be thought of as an extension of the means for defining and measuring

two-dimensional angles in Euclidean space: A two-dimensional angle is measured

by drawing two straight lines that originate at a central point and extend to two

different points on a circle a radial distance r away. The angle subtended between

these two lines is computed by dividing the length of the arc (a) between these two

points by the radial distance (r): a/r. Solid angle is a three-dimensional analog

of this concept: Instead of extending two lines from a point to a circle, a conical

surface is extended from a point to a circular region on a sphere. We then measure

the surface area within this region on the surface of the sphere and divide by the

square of the radius of the sphere. In a more general sense, we are not restricted to

a circular region of the sphere, but rather may define some arbitrarily connected

area on the surface of a sphere that subtends the solid angle. The boundary of this

surface can be connected by an infinite set of straight lines to the center of the

sphere. In either case the solid angle (Ω) associated with this surface is evaluated

by dividing the surface area (A) of the region by the square of the radius (r2) of

the sphere: Ω = A/r2. The units of solid angle are called steradians (sr), though

from a unit analysis standpoint the length dimensions appear to cancel.

Because of the definition of solid angle as the ratio of an area on a sphere to the

radius of the sphere squared, one often encounters the concept of a “solid angle

sphere.” This term refers to the solid angle when the area A consists of the entire

surface area of the sphere of radius r. For this particular case we have A= 4πr2,

and thus

Ω= A/r2 = 4π sr, (4)

is the solid angle associated with a sphere, and can be referred to as a solid angle

sphere. At the other extreme, for a differentially small region of solid angle (dΩ),

there will be essentially only a single propagation direction (Ωˆ) associated with all

the radiation propagating away from a single point (s) that falls within that solid

angular region.

2.2 The Small Angle Approximation

The small angle approximation recasts the three-dimensional RTE into a form

that emphasizes propagation within a small solid angle cone about the forward

direction. However, though we are really interested only in propagation near the

forward direction, we shall see that we may apply the results of this analysis to

propagation in both the forward and backward hemispheres, but the results will be

most accurate in the forward direction.

The rationale behind the small angle approximation is that for many aerosols, fogs

in particular, the majority of scattering of radiation occurs in nearly the forward

direction. Coincidently, this behavior is termed forward scattering (van de Hulst,

1980b). We take the forward direction along the positive z axis.

5
In general, the propagation unit vector Ωˆ is defined in terms of angular variables

θ and φ, as illustrated in figure 1. In this system, θ is a polar/zenith angle and φ is

an azimuthal angle, such that,

Ω=Ω ˆ xˆı + Ωyˆ+ Ωz ˆ

k

= cos φ sin θ ˆı + sin φ sin θ ˆ+ cos θ ˆk. (5)

z

θ

φ

y

x

Ω

Figure 1. Angular coordinate system using two angular measures, θ and φ.

In the forward hemisphere the triplet (Ωx, Ωy, Ωz) can be expressed as,

(Ωx, Ωy, Ωz) =





ωx 

1 + ω2

x + ω2y

, ωy



1 + ω2

x + ω2y

, 1



1 + ω2

x + ω2y



 , (6)

where −∞ < ωx < ∞, −∞ < ωy < ∞ represent components of the vector,

ω = ωx ˆı + ωy ,ˆ (7)

lying in the (ωx, ωy) plane. For propagation into the backward hemisphere the

sign on Ωz is reversed.

Examining the meanings of ωx and ωy, we consider figure 2 which shows the

relationship between θ and ω = |ω | = 

ω2

x + ω2y.

An analysis of this figure’s geometry reveals that for differential dθ, triangles abc

and cde are similar, leading to the relations,

cos(θ) = 1

√1 + ω2 =

√1 + ω2 dθ

dω , sin(θ) = ω

√1 + ω2 . (8)

6
a 1 z

d e

c

b

ω dθ

dω

θ

Figure 2. Geometric relationship between angle θ and tangent ω.

Hence the integral over scattering angle in the RTE can be transformed into an

integration over ω

x and ωy in the ω forward scattering plane:

dΩ = sin θ

dθdφ = ω dω dφ

(1 + ω2)

3/2 = dω

x dωy

(1 + ω2)

3/2 . (9)

We now consider the way in which the phase function representation is also affected by the small angle approximation and the use of ω . The issue of back scattering will also be addressed, but the results of this discussion will not be applied

until chapter 5 which discusses diffuse radiation propagation.

For the moment it is instructive to mention that purely forward scattering may

be sufficient to handle our propagation considerations. Through proper filtering it

may be unnecessary to consider diffuse scattering at all. In particular, techniques

such as (1) time gating, (2) doppler filtering, and/or (3) spatial filtering may eliminate large scattering angle effects. Assuming one or more of these is used, only

a small region in the ω plane need be considered as containing relevant information. Another consideration is that the forward scattering portion of the phase

function is usually very large, and is even larger with aerosols such as fog, where

large mean droplet sizes result in substantially increased probability of forward

scattering.

2.3 Phase Function Representation

The scattering properties of the aerosols considered in addressing the problem of

navigation through fog are detailed in the Army’s Electro-Optical Systems Atmospheric Effects Library (EOSAEL) (Shirkey et al., 1987) Phase Function Database

(PFNDAT) (Tofsted et al., 1997). These scattering properties were evaluated based

on particle size distributions appropriate for fog aerosols and a Mie scattering calculation needed to determine the directional scattering properties of the aerosols.

7
These scattering probabilities are contained in a phase function. The fog aerosols

considered were developed using the Deirmendjian C.1 aerosol model (Deirmendjian, 1969), as applied in Shettle and Fenn (1979) to fog aerosols. This particle size

distribution information was then combined with complex index of refraction information for liquid water and used in A. Miller’s (1983) AGAUS Mie scattering

code. The resulting phase function information could then be analyzed and interpreted in its forward scattering form as detailed in the remainder of this section.

As seen from equation (1), the phase function, P(Ωˆ, Ωˆ

), determines the structure

of the scattered radiance field. This function represents a directional probability

distribution of scattering in any particular outbound direction (Ωˆ) given a value of

incident radiance arriving from direction Ωˆ

. The normalization condition used is

that when integrating a phase function for a fixed direction of incidence (Ωˆ

) over

all outbound directions (Ω= 4π sr) the result is unity:



Ω=4π

P(Ωˆ, Ωˆ) dΩ=1. (10)

Another property common to all phase functions is based on the reciprocity principle (van de Hulst, 1980a): The phase function is constrained such that P(Ωˆ, Ωˆ

) =

P(Ωˆ

, Ω) ˆ . In other words, all phase functions must be symmetric between input

and output directions.

For some species of scatterers, the phase function may depend on complete knowledge of both the incident and scattered propagation directions (e.g., snowflakes,

raindrops). But for hydrosols, such as fogs, the particle shapes are nearly spherical and, hence, symmetry arguments reduce the phase function dependence to

only the scattering angle θs. Sometimes this dependence is written in terms of

μs = cos θs = Ωˆ · Ωˆ

.

For the particular problem of interest, scattering through heavy haze/fog, while

models of the phase function in its complete form exist (Tofsted et al., 1998), a

more usable representation is desired. One such form is a purely Gaussian phase

function (c.f., Tam and Zardecki, 1980; Zardecki and Gerstl, 1986). In this form

the phase function is modeled as purely forward scattering, using a Gaussian function, modeled here as

Pg(∆ω) = α2 exp 

−π α2 ∆ω2



, (11)

where

θs ≈ ∆ω =

	

	

	

−→∆ω

	

	

	 = |ω − ω 

| , (12)

and α is an inverse width parameter.

A useful feature of the Gaussian phase function is that it is separable:

α2 exp 


−π α2 ∆ω2



= α2 exp 
−α2∆ω2

x



exp 
−α2∆ω2

y



. (13)

The advantage of this feature is seen more clearly when the RTE is modified for

the small angle approximation and transformed into the Fourier domain. Also note

8
that Pg exhibits the following normalization condition:

2π

∞

0

Pg(ω, α) ω dω = 1, (14)

such that the integral of Pg over the ω plane normalizes to unity.

However, it is difficult to use a purely Gaussian phase function representation

to model an actual phase function because there is usually more than one critical scale size responsible for the behavior of the phase function. Typically, phase

function profiles exhibit a highly peaked and narrow forward-scattering region,

due to diffraction, and a broader secondary peak extending to 50–60◦. These features are illustrated in figure 3, which shows plots of the 1.5 μm and 1.54 μm

wavelength phase functions for advection fog.

103

102

10

1

0.1

10–2

10–3

Phase Function

0 30 60 90 120 150 180

Scattering Angle θs

1.50 μm

1.54 μm

Figure 3. Advection fog phase functions for 1.5 μm and 1.54 μm wavelengths.

To accurately portray the phase function structure it will be necessary to formulate

the phase function in a form which is both accurate and amenable to computation.

Both goals are achievable if we can express the phase function as a sum of Gaussian terms for each hemisphere individually using,

P ̃f (∆ω) = 

M

m=1

Am Pg(∆ω, αm); (15)

P ̃b(∆ω) = 

N

n=1

Bn Pg(∆ω, βn). (16)

9
As a normalization condition we require,

1 = 

M

m=1

Am + 

N

n=1

Bn, (17)

such that the combined Gaussian components from both hemispheres integrate to

unity, which reflects our phase function normalization condition.

However, there are two complications which tend to cancel each other: (1) In

assessing the fit coefficients for the Gaussian series in both hemispheres we note

that the phase function approaches a nonzero value at θs = 90◦; (2) The phase

function integral contains the additional factor (1 + ω2)−3/2 in the integration

kernel:



4π

I(s, Ωˆ) P(Ωˆ, Ωˆ) dΩ =



∞

−∞

Pf (|ω − ω |)

(1 + ω2)

3/2 If (s, ω

) dω

+



∞

−∞

Pb (|ω + ω |)

(1 + ω2)

3/2 Ib (s, ω) dω. (18)

Here, ω  denotes integration over the backward hemisphere where the phase function peak occurs where ω  = −ω .

In this form the factors (1 + ω2)−3/2, appearing in the denominators, reflect a decreasing weight with angular separation. But they would significantly complicate

the problem when the Fourier transform is taken (section 2.5). To avoid this situation, since |ω | is small due to the small angle approximation, we may substitute

the integration with respect to ω  and ω  with integrations with respect to −→∆ω



and −→∆ω



such that



4π

I(s, Ωˆ) P(Ωˆ, Ωˆ) dΩ =



∞

−∞

Pf

	

	

	

−→∆ω



	

	

	






1 + −→∆ω

2

3/2 If (s, ω

) d

−→∆ω



+



∞

−∞

Pb

	

	

	

−→∆ω

	

	

	






1 + −→∆ω

2

3/2 Ib (s, ω) d

−→∆ω



, (19)

allowing us to define the modified phase functions,

P ̃f (∆ω

) = Pf (∆ω

)

(1+∆ω2)3/2 ; P ̃b(∆ω) = Pf (∆ω)(1+∆ω2)3/2 , (20)

and attempt to fit Gaussian series to these functions. Gaussian approximations to

the functions P ̃f and P ̃b are shown in figure 4 for the forward hemisphere and in

figure 5 for the backward hemisphere.

In the forward hemisphere approximation we used nine Gaussian components.

In the backward hemisphere, we used eight Gaussian components. Both of these

10
104

102

1

10–4

10–2

10–8

10–6

Modified Phase Function

012345

Square Root of Tangent (∆ω¢

1/2)

Forward Phase Function

Gaussian Approximation

Figure 4. Forward hemisphere portion of phase function for advection fog at 1.5 μm wavelength.

104

102

1

10–4

10–2

10–8

10–6

Modified Phase Function

0.0 1.0 2.0 3.0 4.0 5.0

Square Root of Tangent (∆ω≤1/2)

Backward Phase Function

Gaussian Approximation

Figure 5. Backward hemisphere portion of phase function for advection fog at 1.5 μm wavelength.

11
functions appropriately approach zero at infinity, but the series fit is only designed

to be valid for scattering angles ranging from 0–78◦. The modified phase function

approach is completely accurate in the forward and backscattering directions, but

exhibits increasing error for off-axis directions. Since we are primarily interested

only in propagation along the main axis, this should not be a critical problem. The

Gaussian expansion coefficients used are listed in table 1.

Table 1. Gaussian expansion coefficients for forward and backward hemisphere phase function representations of the 1.5 μm advection fog phase function.

m, n Am αm Bn βn

1 0.3423 29.38 0.00009 16.22

2 0.1144 9.911 0.0158 1.7150

3 0.2676 1.584 -0.0463 1.2810

4 0.1884 0.7401 0.0554 0.8540

5 0.0377 0.3741 0.0051 0.4076

6 0.0082 0.2079 0.0025 0.1775

7 0.0034 0.1243 0.0016 0.0800

8 0.0017 0.0667 0.00093 0.0290

9 0.00107 0.0290 — —

2.4 Small Angle RTE

The key change in the RTE occurs when the directional variables in the RTE are

replaced by functions of ω and ω 

. When these functions are substituted into the

RTE, it may be written as,

ω √1 + ω2 ·

∂I

∂r

+

1

√1 + ω2

∂I

∂z +σ I =  σ 

∞

−∞

P ̃f

	

	

	

−→∆ω



	

	

	



I (r, ω, z) d

−→∆ω



. (21)

Note that ω  = ω − −→∆ω



such that the phase function integral appears as a convolution.

This equation can be further simplified by assuming that in general ω will be small

in magnitude, allowing us to ignore the √1 + ω2 factors appearing on the LHS.

This simplification leads to the equation,

ω ·

∂I

∂r

+

∂I

∂z + σ I =  σ 

∞

−∞

P ̃f

	

	

	

−→∆ω



	

	

	



I (r, ω, z) d

−→∆ω



. (22)

2.5 Fourier Transform Operation

The key advantage of writing the RTE in the small angle form is that the angular

integration on the RHS of the modified RTE is in the form of a two-dimensional

12
convolution integral. Convolutions Fourier transform into products of the transforms of the two convolved functions. It would therefore seem expedient to Fourier

transform the RTE. The transform process is conducted in four dimensions: x, y,

ωx, and ωy. Let κ be the frequency space representation of r = xˆı + yˆ, and let η

be the frequency representation of ω .

To perform these transformations, Gaskill’s (1978) version of the Fourier transform will be used. In this version,

F(κx) = ∞

−∞

f(x) exp [−j2πκxx] dx. (23)

f(x) = ∞

−∞

F(κx) exp [j2πκxx] dκx, (24)

with j = √−1. F(κx) is the Fourier transform of the real valued function f(x).

Gaskill’s corresponding Hankel transforms are given by

G(κ)=2π

∞

0

g(r) J0 [2πκr] r dr, (25)

g(r)=2π

∞

0

G(κ) J0 [2πκr] κ dκ, (26)

where G(κ) is a radially symmetric frequency representation of the radially symmetric function g(r).

Given these forms for the Fourier and Hankel transforms, we consider some general transform rules: First, if g(r) has the Hankel transform given by G(κ), then

g(r/b) transforms as

H



g


r

b

 = |b|

2 G(bκ). (27)

This rule is useful in transforming the phase function, since

H 

exp 
−π ω2

 = exp 


−π η2



. (28)

Hence, we can write the Fourier transform of the phase function as,

Pˆf (η) = H{P ̃f (ω)} = 

M

m=1

Am exp 
−π η2/α2

m



. (29)

Several other rules relating to Fourier transforms are also of significance: If f(x)

has the Fourier transform F(κx), then,

F {f(x ± x0)} = F(κx) exp [±j2πx0κx] ; (30)

F 

f(k)(x)



= (j2πκx)

k F(κx), (31)

13
where f(k)

(x) represents the kth derivative of f(x) with respect to x. Due to the

similarity between the Fourier and inverse Fourier transforms,

F {f(x) exp [±j2πx0κx]} = F(κx ∓ κ0); (32)

F {(−j2πx) f(x)} = F(k)

(κx); (33)

A final rule applies to convolutions: If f(x) and g(x) Fourier transform to F(κx)

and G(κx), respectively, then their convolution,

f(x) ∗ g(x) = ∞

−∞

f(β) g(x − β) dβ, (34)

transforms as,

F {f(x) ∗ g(x)} = F(κx) G(κx). (35)

These rules are used in transforming the small angle form of the RTE. Using the

derivative rule, we transform the term ωx ∂I/∂x with respect to x as,

Fx



ωx

∂I

∂x

= j2πκxωx  ̃Ix, (36)

where  ̃Ix is the transform of I with respect to x. The resulting j2π factor is used

in transforming with respect to ωx,

Fωx



j2πκxωx  ̃Ix



= −κx

∂ ˆIxωx

∂ηx

, (37)

where ˆIxωx is the transform of I with respect to both x and ωx. Completing this

derivation, we find,

F4



ωx

∂I

∂x

= −κx

∂




F2

y ωy



ˆIxωx



∂ηx

= −κx

∂ ˆˆI

∂ηx

, (38)

where ˆˆI is the fourth order Fourier transform of the radiance function. A similar

procedure is followed to obtain the transform of the term ωy ∂I/∂y.

The convolution rule is used when transforming according to ω on the scattering

integral term. In addition, we use linear properties of the operators to reverse the

order of operation similar to the procedure used above where the F2

y ωy operator

was switched with the differential operator. Here, following the F2

ωx ωy operation,

the transformed phase function is factored out of the F2

x y transformation step since

it is considered a constant with respect to the x and y dimensions:

F4









∞

−∞

P ̃f (|ω − ω 

|) I (r, ω, z) dω





 = Pˆf

ˆˆI. (39)

14
The net result of the Fourier transform process is that the RTE can now be expressed as,



−κ · ∂

∂η

+

∂

∂z + σ − σPˆf



ˆˆI = 0, (40)

where the dot product on the LHS is an abbreviation for,

−κ · ∂

∂η = −κx

∂

∂ηx

− κy

∂

∂ηy

. (41)

The notation in equation (40) uses a singly carated variable to symbolize a doubly transformed function (in particular, ωx and ωy transformed), while a doubly

carated variable has been quadruply transformed (x, y, ωx, and ωy). Solving for

the radiance in Fourier space is viewed as considerably simpler than solving the

RTE in its original form. It is then a matter of recovering the radiance field via an

inverse Fourier transform process:

I(r, ω, z) = 

∞

−∞

dκ



∞

−∞

dη ˆˆI(κ, η, z) exp {j2π (κ · r + η · ω )} . (42)

2.6 Gaussian Beam Model

The method proposed for solving the navigation through fog problem uses a laser

source operating in the eye-safe near-IR spectral region, with a wavelength λ of

approximately 1.5 μm. This laser source will emit a beam focussed approximately

1 m from the laser exit aperture to enable a wide beam divergence. This condition

greatly simplifies the mathematics of the propagating laser beam as shown in the

following section.

2.6.1 Complex Beam Wave

Siegman (1986) provides a standard numerical model for describing the propagation of a TEM00 mode Gaussian beam. This form should be equivalent to the

vacuum solution of the RTE. However, due to differences in terminology it will be

necessary to translate out of the common description used in laser literature to a

form that is compatible with the radiance description commonly used in radiative

transfer.

We begin by describing the commonly used complex beamwave form. Siegman

describes his beamwave in terms of u ̃, the beamwave complex amplitude:

u ̃(x, y, z) = 

2

π

q ̃0

w0 q ̃(z) exp 

−jkz − jk r2

2 ̃q(z)



, (43)

where q ̃(z) is the complex radius of curvature for the Gaussian beamwave, w0

is the beam width at zero range, q ̃0 is the complex radius at zero range, k is the

radiation wavenumber (2π/λ), and j = √−1 and r = √x2 + y2 as before. This

15
beamwave form assumes that the minimum beam width occurs at a range of zero

and the beamwave travels with its center along the positive z axis.

In Siegman’s derivation the paraxial approximation is used. While recent observations by Zeng et al. (1999) indicate some correction to the paraxial approximation

is necessary, we can avoid these corrections by the choice of vertical polarization

of the emitted beam and the consideration that typical off-axis positions of interest

will lie primarily along the horizontal axis.

Of considerable interest is the form taken by q ̃(z), given by Siegman as a function

of w(z), the beam width at range z, R(z), the real-valued radius of curvature at z,

and wavelength λ through the expression,

1

q ̃(z) = 1R(z) − j

λ

πw2(z)

, (44)

where the radius R and width w parameters are given by,

R(z) = z +

z2

R

z

, w2

(z) = w2

0



1 + 
 z

zR

2



. (45)

These expressions, in turn, hinge on definitions of zR, termed the Rayleigh range,

and w0, the minimum beam width which is assumed to occur at range z = 0. zR

and w0 are related through the equation

zR = πw2

0

λ . (46)

Through these definitions, the equation for q ̃(z) can be written as

q ̃(z) = z +  ̃q0 = z + j zR. (47)

2.6.2 Radiance Form for a Diverging Gaussian Beam

To use functions u ̃ and q ̃(z) in a radiative transfer application we will need to deconstruct these functions to determine their relationship to the radiance function

I used by the radiative transfer method. To begin, we multiply the complex amplitude by its complex conjugate. This produces a real-valued flux density for the

propagating beamwave:

I(r, z) = 2P

πw2 exp 

−2

r2

w2



, (48)

where we use w = w(z) from equation (45). Integration of I(r, z) across the

(x, y) plane reveals that P (power) represents the energy flowing past a point z

through the (x, y) plane per unit time (given in Watts). I(r, z) thus has units of

irradiance (W/m2). A normalized form of the I function is shown in figure 6,

which is essentially equivalent to Siegman’s figure 17.2. This figure includes the

dashed outline of a “top hat” beam. The top hat beam has the property of constant

16
1.5

1.0

0.5

0.0

I(r)/I(0)

–2 –1 0 1 2

r/w

equivalent

“top hat”

beam

circular

gaussian

beam 

ATH = πw2/2

e–2r2/w2

Figure 6. Gaussian beam transverse irradiance relationship.

irradiance across a circular beam profile of radius w/√2, with irradiance equal to

the maximum of the related Gaussian beam, and exhibits the same total flux as

the Gaussian beam. Hence, the top hat beam is sometimes used to approximate

properties of the Gaussian beam yet is simpler to work with mathematically.

The I(r, z) form, however, only provides magnitude information. It does not provide the directional information needed to specify ω functional dependence. To

assess this additional dependence, the complex wavefront curvature, q ̃(z), equation must be analyzed. A key part of this analysis involved examination of Siegman’s figure 17.1 (see figure 7). This figure appears to show that a single beam

curvature, R(z), applies for every point in each plane perpendicular to the z axis,

where R(z) is given in equation (45).

More rigorously, introducing the definition of q ̃−1(z) from equation (44) into

equation (43) results in

u ̃(x, y, z) = 

2

π

q ̃0

w0 q ̃(z) exp 

−jkz − jk r2

2 R(z) − r2w2(z)



. (49)

This form reveals a separation between real and imaginary structures as a function

of off-axis distance, r. The exp(−r2/w2) term accounts for the off-axis amplitude

reduction, while exp(−jkr2/2R) accounts for the off-axis phase modulation. This

form is based on the Fresnel approximation, utilizing a wavefront curvature of R

and the assumption that r  R.

Given this constraint, it is nevertheless clear that use of the phaseterm

exp(−jkr2/2R) implies that at any given point (x, y, z) the wavefront of the prop17
w0

waist

(z=0)

w(z)

R(z)

Figure 7. Gaussian beam transverse spread (w) and curvature (R) as functions of range (z).

agating beam will have some normal vector which we can model using the triplet

{x/R(z), y/R(z), 1}. From the standpoint of adding a directional dependence to

the radiance, this normal vector leads to the addition of a two-dimensional Dirac

delta function: δ {ω − [r/R(z)]} to describe the directional properties:

δ



ω − r

R(z)



= δ



ωx − x

R(z)



δ



ωy − y

R(z)



. (50)

The vacuum form (IV ) for the radiance of a propagating Gaussian beamwave

should then be written as,

IV (r, ω, Z) = 2P

π w2(Z) exp  −2r2w2(Z)



δ



ω − r

R(Z)



, (51)

where we again mention the caveat that the beam is rapidly diverging and travelling in the positive z direction.

2.7 Radiative Transfer Solution

Equation (40) in section 2.5 provided a Fourier version of the small angle approximation of the RTE. In this section this equation is solved using a method similar

to that outlined by Smirnov (1964). Further, the solution of this process must reduce to the equation (51) vacuum Gaussian beamwave solution in the absence of

aerosols (σ = 0);

Equation (40) consists of a series of linear operators applied to the transformed

radiance function. The results of this analysis are similar to those of Zardecki

18
and Tam (1982), but utilize a different Fourier transform kernel, modified source

radiation, and a more accurate phase function description. Also here, a detailed

description of the mathematics is provided.

2.7.1 Implicit Solution Method

Smirnov’s approach involves the solution of a single linear partial differential

equation through its relationship to the solution of a set of linear differential equations. This approach begins by considering a set of variables: x0, x1, ..., xn. These

variables are related to one another within a system of differential equations:

dx0

X0

= dx1

X1

= dx2

X2

= ... = dxn

Xn

. (52)

The Xi are functions of the xi (X1 = X1 (x0, x1, ..., xn), X2 =

X2 (x0, x1, ..., xn), etc.). There are n + 1 equal ratios within this system of

relations.

Within this system there are effectively n dependent variables and one independent variable, though the form used for the relations in the system is the same for

all n + 1 variables. Nonetheless, one of these variables could be chosen as the

independent variable. For example, if the x0 variable is chosen to be the independent variable, then the other variables can be written as x1(x0), x2(x0), ..., xn(x0).

If we are able to write out explicit equations for each of the dependent variables

as functions of the single independent variable then the system of equations (52)

is solved.

Smirnov posits such a system of equations and its solution and then makes a second point: He assumes that given such a system of equations one can produce

a function, call it Ξ(x0, x1, ..., xn). This function has the property that it is not

identically zero. Yet, if any solution to the system of equations (52) is substituted

for all of the dependent variables in that system, then Ξ evaluates to a constant.

Smirnov calls such a function an integral of the system of equations.

Viewed another way, the function Ξ can be considered to define a hypersurface

within the hyperspace (x0, x1, x2, ..., xn). This hypersurface corresponds to a

solution of the system of equations (52) that satisfies the condition,

Ξ(x0, x1, x2, ..., xn) = C. (53)

In the analysis to follow, Smirnov exploits the properties of this notional function,

Ξ. However, given the nature of this analysis, he never needs to actually solve

explicitly for function Ξ in terms of the set of variables x0, ..., xn. The first step in

this exploitation process involves generating the total derivative of Ξ with respect

to the chosen independent variable x0.

dΞ

dx0

= ∂Ξ

∂x0

+

∂Ξ

∂x1

dx1

dx0

+

∂Ξ

∂x2

dx2

dx0

+ ... +

∂Ξ

∂xn

dxn

dx0

. (54)

19
Or in condensed notation,

dΞ

dx0

= n

i=0

∂Ξ

∂xi

dxi

dx0

. (55)

In this equation the terms dxi/dx0 merely specify the rate of change in one of the

variables (xi) relative to the change in the variable x0. One may imagine some

random positional displacement, −→dx, within the hyperspace. Essentially this displacement is completely arbitrary, and hence the change in Ξ may also be nonzero.

However, if the direction of displacement is determined by calculating the quantities dxi/dx0 using solutions of the system of equations for the xi in terms of x0,

then the net change in Ξ must equal zero! Multiplying both sides of equation (55)

by dx0 we generate,

dΞ

	

	

	

S = n

i=0

∂Ξ

∂xi

dxi

	

	

	

S

= 0, (56)

where the use of the subscript S implies that the dependent dxi are constrained to

vary only along the solution hypersurface.

To exploit the development in equation (56) an important auxiliary development is

now required. In this development we assign a quantity, dλ, to all the equal ratios

given in equation (52):

dλ = dx0

X0

= dx1

X1

= ... = dxi

Xi

= ... = dxn

Xn

. (57)

Since dλ is a nonzero quantity, it can be divided into equation (56), and terms can

be cancelled, to obtain,

n

i=0

∂Ξ

∂xi

dxi|

S

dλ = n

i=0

∂Ξ

∂xi

Xi = 0. (58)

We now compare the appearance of equation (58) to equation (40), the Fourier

version of the RTE. The form of equation (58) is nearly identical to the equation

to be solved except that the partial derivatives in equation (58) only appear as

operators on the integral function Ξ, whereas the equation (40) has the form,

Y0

∂xn

∂x0

+ Y1

∂xn

∂x1

+ ... + Yn−1

∂xn

∂xn−1

+ Yn = 0, (59)

where the Yi are functions of x0, x1 , ... , and xn, and the partial differentials

operate on the variable xn, which itself appears in the Yi functions. By contrast, the

Xi functions do not depend on Ξ. It is therefore necessary to relate equation (59)

to the form of equation (58) by replacing the partial derivatives of xn with partials

involving Ξ.

To perform this transformation, equation (56) is again examined under the condition of holding all but two variables constant. Let the two free variables be chosen

as xi and xn. Therefore, from equation (56),

∂Ξ

∂xn

dxn +

∂Ξ

∂xi

dxi = 0. (60)

20
or,

dxn

dxi

= − ∂Ξ/∂xi

∂Ξ/∂xn

(61)

In this case only xn and xi are allowed to vary and thus this result is a partial

derivative (∂xn/∂xi). Smirnov calls this the rule for differentiation of implicit

functions. This rule allows the partial derivatives of xn to be replaced by partials

with respect to Ξ. Equation (59) is thus transformed into,

−Y0

∂Ξ

∂x0

− Y1

∂Ξ

∂x1

− ... + Yn

∂Ξ

∂xn

= 0. (62)

We thus have a method of relating equation (59), which has the same form as the

differential equation we wish to solve, to equation (58), an equation which we

have causally connected to the system of equations (52). In particular, we have

established connections between −Yi ↔ Xi, for all i<n; and Yn ↔ Xn. Any

partial differential equation that can be written as a system of linear operators in

the form of equation (59) may be related to a system of equations in the form of

equation (52). The expectation is then that this related system of equations will be

simpler to solve than the original partial differential equation.

When we apply this series of steps to our problem the first step is to assign z as

the independent variable (x0). The remaining variables are assigned as, x1 = ηx,

x2 = ηy, and x3 = xn = ˆˆI. Next, Yi functions are assigned:

−X0 = Y0 = 1;

−X1 = Y1 = −κx;

−X2 = Y2 = −κy;

X3 = Y3 = 


σ − σ Pˆf

 ˆˆI. (63)

Introducing expressions for the Xi’s into the related system of equations generates

three independent relations, the first two of which are

−dz = dηx

κx

, −dz = dηy

κy

. (64)

Integration of these first two functions produces a relationship between the angular

frequency vector, η, the spatial frequency vector, κ, and range z:

η(z) = Cη − κz, (65)

where Cη is a constant of the integration.

This result appears to indicate that η is a completely dependent function of κ and

z, but this is not the case. Rather, η is only dependent on z and κ in terms of its

positional variability. In fact, Cη may freely vary as long as it is independent of z.

In particular,

η(0) = Cη = η0. (66)

21
Essentially the restrictions on η imply that propagation under the radiative transfer

model involves a corpuscular-type assumption — that once a photon is launched in

a given direction it does not arbitrarily begin propagating in some other direction

as it travels forward in z. Thus, diffraction effects are not treated directly. This is

a limitation of the radiative transfer method which ignores large-scale diffraction

effects. However, small-scale diffraction is modeled indirectly through use of the

scattering phase function.

A second choice for defining Cη is at path position z = Z. At this range,

η(Z) = ηZ = η0 − κZ, (67)

Cη = η0 = ηZ + κZ. (68)

These two means of expressing Cη are significant because we will deal with specific examples of propagation where we are interested in the propagated radiance

field at ranges 0 and Z, and defining Cη in terms of these ranges is important in

understanding the nature of the propagating radiance relative to these particular

ranges. Once Cη is determined, equation (65) indicates η’s range dependence at

all other ranges. In the forward propagating laser beam case, ηZ is used as the angular frequency variable since the properties of the propagated signal as it enters

the retroreflector will be of most interest.

Having thus chosen the ηZ representation for Cη, η(z) is written as,

η(z) = ηZ + κ (Z − z). (69)

2.7.2 Fourier Solutions

The third and final relationship derived from equation set (63) is given by,

dˆˆI

ˆˆI

= −σ (1 − Pˆf ) dz, (70)

which has the solution,

ˆˆI = ˆˆI0 exp





−



Z

0

σ




1 − Pˆf



dz





 . (71)

Note that, similar to the case with η, ˆˆI0 need only be a constant with respect to the

z variable. This is illustrated as we consider several solutions to this equation for

various propagation conditions.

2.7.2.1 Vacuum Solution

As observed previously, the formula given for IV in equation (51) must be the

solution of equation (71) in the case of a vacuum atmosphere (i.e., when σ = 0).

22
But when σ = 0 the argument of the exponential in equation (71) equals zero. In

this case ˆˆI = ˆˆI0 in the Fourier representation of the solution. Thus, the Fourier

transform of IV must be ˆˆI0, a constant with respect to z according to the radiative

transfer method.

Rewriting the vacuum solution at range z = Z,

IV (r, ω, Z) = 2P

π w2(Z) exp  −2r2w2(Z)



δ



ω − r

R(Z)



, (72)

this function Fourier transforms to,

ˆˆIV (κ, η, Z) = P exp 

−π2 w2 Λ2

2



, (73)

where Λ is the magnitude of the vector,

Λ =  κ + η/R. (74)

Yet both R (the radius of wavefront curvature from equation (45)) and η are functions of position, z, along the path. At first glance, then, IV appears to violate

the requirement that its transform be constant with respect to range. This difficulty arises because the actual equation for the propagating beamwave does not

truly follow radiative transfer theory’s corpuscular propagation assumption. While

this difference cannot be completely resolved, wΛ is approximately constant with

range as long as zR  z. Considered another way, as long as πw2

0  λz, such

that the characteristic width of the beam at its focal point is small compared to a

Fresnel zone length, then the beam can be considered to be nearly a point source.

This condition is possible for the laser system contemplated since the beam is

focused at a short range from the laser aperture.

Under these conditions,

w(z) = w0



1 +

z2

z2

R

≈ w0

z

zR

, (75)

where w0 varies as,

w2

0 = λ zR

π −→ w0zR

= λ

π w0

. (76)

Conversely, Λ varies approximately inversely with range,

Λ =  κ + [ηZ + κ (Z − z)]

z [1 + (z2

R/z2)] ≈ κ + ηZ

z

+ κ

Z

z − κ = (ηZ + Z κ)z = η0z . (77)

Combining the influences of w and Λ, the argument to the exponential in the

transformed vacuum radiance is given by,

−π2

2

w2 Λ2 ≈ −

λ2 η2

0

2 w2

0

, (78)

which is independent of range z. Thus, under the limitations mentioned, the vacuum solution of equation (72) is suitable as a solution to the RTE in Fourier space.

23
2.7.2.2 Unscattered Solution

Once the vacuum form of the solution is known, the unscattered solution may be

immediately written as,

IU (r, ω, Z) = IV (r, ω, Z) exp [−τ (Z)] , (79)

where

τ (Z) = 

Z

0

σ dz, (80)

is the optical depth encountered by a beam travelling through distance Z.

In Fourier space this equation is written as,

ˆˆIU = ˆˆIV exp [−τ (Z)] , ˆˆIV = ˆˆI0. (81)

This result is obtained by setting the single scattering albedo, , to zero.

2.7.3 Scattered Solution

The solution for the scattered radiation case is considerably more complex. We

begin by resolving it into a summation of scattered and unscattered components:

ˆˆI = ˆˆIU + ˆˆIS

;

= ˆˆIU + ˆˆIU [S − 1] ;

= ˆˆIU + ˆˆIU



exp









Z

0

σPˆf dz





 − 1



 , (82)

where ˆˆIS is the scattered solution and S characterizes the spatial and angular

frequency effects of the scattering.

It will be necessary to evaluate S in order to determine its influence on the radiance field arriving at the retroreflector. To simplify this task, thereby reducing the

number of free parameters involved in this study, it is assumed that the scenario

studies will be one involving an aircraft landing at an airfield enshrouded in fog.

As the aircraft approaches the field, at some point the plane will be above the fog.

Then, as the plane makes its final approach, it becomes completely immersed in

the fog. Thus the optical path may be modeled to first order as consisting of an

initial portion free of fog and a latter portion (containing the retroreflector) immersed in fog. The fraction of the modeled volume that is free of fog will range

from zero (completely free of fog) to unity (completely immersed in fog). For

simplicity, the volume is assumed to have uniform density and scattering properties throughout the fog region. Zf is designated as the distance from the laser focal

point to the forward edge of the fog enshrouded portion of the modeled volume.

The dimensionless quantity Q = (Z − Zf )/Z measures the fraction of the optical

24
path containing fog. The fog aerosol properties are assumed constant in both x

and y, as characterized by parameters σ, , and the coefficients constituting the

phase function description. The remainder of the optical path is assumed to be

filled with an optically thin aerosol whose properties represent a virtual vacuum.

Given these limiting approximations, the scattering solution reduces to an evaluation of the function S and the inverse transform of the resulting function. Because

the forward-scattering phase function was assumed to be a sum of Gaussian terms,

it is possible to exactly represent S.

2.7.4 Error Function Application

In solving for S, integrations involving Gaussian functions must be performed.

These integrals use the error function, given by

erf(x) = 2

√π

x

0

exp −t

2



dt. (83)

Defining a Gaussian function according to Gaskill’s (1978) definition,

Gaus(x) = exp 

−π x2



,

the integral of the Gaussian can be written as,

Υ (x2, x1) = x2

x1

Gaus(u) du = 1

2




erf √π x2



− erf √π x1

 . (84)

Integrating along the path for the scattering effects in S,

S = exp







 σ 

Z

Zf

Pˆf (|ηZ + κ [Z − z]|) dz







= $

M

m=1

exp 

Am  τ

Qqm

exp 
−π p2

m



1 − μ2

 Υ [Qqm + μpm, μpm]



, (85)

where τ = σQZ, qm = Zκ/αm, pm = ηZ/αm, μ = (κ · ηZ)/(κηZ) is the cosine

of the angle between the κ and ηZ vectors, and κ and ηZ are magnitudes of their

respective vectors.

Using this result, the solution in Fourier space may be written as

ˆˆI = P exp 


−AG



q2 + 2μqp + p2

 exp(−τ ) ×

$

M

m=1

exp 

Amτ

Qqm

Υ [Qqm + μpm, μpm] exp 
−π p2

m



1 − μ2



, (86)

where

AG = λ2

2w2

0

. (87)

P, AG, Q, τ , and  are constants of the laser propagation scenario, while q, p, and

μ provide the variability associated with κ, and ηZ in terms of unitless parameters.

25
2.7.5 The Inverse Transform

Because the Fourier representation of the radiance solution has factors whose

magnitudes peak at different points in frequency space, numerical methods will

be necessary to evaluate the inverse-transformed radiance field. In its full form,

the inverse transform is given by

I(r, ω, Z) = 

∞

−∞

dκ



∞

−∞

dηZ

ˆˆI(κ, ηZ, Z) exp {j2π (κ · r + ηZ · ω )} . (88)

However, processing this transform can be simplified by reducing the order of

the integral by one. This step is possible because equation (86) only involves the

quantities p, q, and μ. Here, let μ = cos(φ) define the angle, φ formed between

the vectors κ and η. Given the symmetries of the problem, it is possible to solve

for the radiance field in terms of r = |r|, ω = |ω |, and the angle formed between r

and ω (call it a), instead of r and ω in fully four-dimensional form. This reduction

acknowledges the radial symmetries that exist at off-axis positions.

To integrate out one of these four dimensions, consider figure 8. From the geometry in this figure, the dot products contained in the complex argument of the

exponential in equation (88) can be written,

κ · r = κ r cos(a + c − φ − ψ), (89)

ηZ · ω = ηZ ω cos(c − ψ). (90)

ω

y

c

a b

r

x

ψ

φ

κ

ηz

Figure 8. Angular relations between various vectors in transform and real coordinates.

The form that these dot products take suggests that if the multiple integral in

equation (88) were expressed in polar form that the ψ variability in the integral

26
can be immediately integrated out. In polar form the inverse transform becomes



∞

−∞

dκ



∞

−∞

dηZ[...] = ∞

0





∞

0







2π

0







2π

0

[...] dψ



 dφ



 η dη



 κ dκ, (91)

where the ψ variability has been carried to the innermost integration level. What

it means to integrate ψ first can be explained best with reference to figure 8. Performing the ψ integration while keeping all the other variables constant consists of

holding the angle of separation between vectors κ and ηZ constant while rotating

this fixed structure through 2π radians relative to fixed r and ω variables. This rotation covers a complete 2π radians, implying that the start and finish points of the

integration are irrelevant, a fact that is useful in evaluating the resultant integral.

To accomplish the integration, the first step is to note that the Fourier transformed

radiance field does not depend on ψ. It can therefore be moved outside the ψ

integral, leaving only the complex exponential to integrate. Second, the quantity

β is introduced, where

β = a + c − φ − ψ. (92)

Since φ is constant while ψ is varied, dβ = −dψ. Therefore, upon replacing ψ

with β, the innermost integral becomes



2π

0

exp {j 2π [κ r cos(β) + ηZ ω cos(β + φ − a)]} dβ. (93)

This integral is solved analytically using a result from Goodman (1968):



2π

0

exp {−j L cos(θ − φ)} dθ = 2π J0(L), (94)

where J0 is a zero order Bessel function of the first kind. Integral (93) can be

related to integral (94) by observing that the quantity

− 1

2π

L cos(θ − φ

)

is simply the x-component of a vector,

u = 1

2π

(−L cos [−φ

] , −L [− sin φ]) (95)

that is being rotated about the origin, with an offset angle θ

.

Equation (93) geometry is more complex, but the result can still be expressed in

terms of Goodman’s solution. The geometry of equation (94) is shown in figure 9.

The vector U in that figure is composed of the vector sum of

−κr = (κr, 0), (96)

27
y

U

κ r β

x

φ – a

ηzω

Figure 9. Orientation of vector U relative to component vectors −κr and −η

−Zω.

and −η

−Zω = (ηZω cos(φ − a), ηZω sin(φ − a)). (97)

Goodman’s integral result is independent of the value of φ since a unit circle is

being traversed. Similarly, U ignores any initial offset between the x axis and the −κr vector. The only quantity of interest is the magnitude of U , given by

U(κr, ηZω, φ − a) = 	

	

	U

	

	

	 =



κ2 r2 + 2 κrηZ ω cos(φ − a) + η2

Z ω2. (98)

Hence, the innermost integration of equation (91) evaluates to



2π

0

exp {j 2π [κ r cos(γ) + ηZ ω cos(γ + φ − a)]} dγ = 2π J0(2π U). (99)

Eliminating this inner integral significantly simplifies the calculations necessary

to evaluate the radiance field for a specific (r, ω, a) triplet. Simultaneously, integration with ψ eliminates the need to work with a complex exponential kernel

in the inverse transform. The resulting Bessel function ensures that the resulting

integration will be purely real-valued.

This result also has immediate implications for the resulting inverse transform

when both r and ω are zero. For this case the J0 term evaluates to unity and all

dependence on the angle a drops out of the resulting inverse transform.

2.8 Numerical Integration

Having removed one order of integration, we must now evaluate the triple integral

involving κ, ηZ, and φ.

28
To accomplish this task various numerical methods were considered, such as those

contained in Press et al. (1992). However, of the available sources of routines none

were found to be applicable. Hence, this section describes a new technique for

treating this task.

First, as previously stated, different terms in the integration kernel have peaks

at different locations in the integration volume. For instance, the J0 term peaks

where U(κr, ηZω, φ − a)=0. This condition occurs along the plane κr = ηZω

since there will always be a φ value for which cos(φ−a) = −1. Hence, following

integration by φ, there will be a resonance along κr = ηZω. (This resonance also

encompasses the specific condition: κr = ηZω = 0.)

On the other hand, the term

exp 

− λ2

2w2

0

(Z κ + ηZ)

2



(100)

has a resonance when Zκ = ηZ.

Comparing these two expressions, there is a double resonance for

ηZ

κ = rω

; Z = ηZ

κ −→ ω = rZ . (101)

Not surprisingly, this peak describes the direction of geometrical propagation for

unscattered radiation.

The remaining S term in equation (86) provides the scattering dependence for ˆˆI.

This function has a very complex behavior, depending on the particular combination of parameters describing the aerosol and path. Basically, all the features

describing the optical path are here: the thickness of the aerosol, the scattering

properties, the portion of the path in the fog, and the angular and positional structure.

Even so, a few generalizations are possible. For example, when μ = −1, the Υ

function tends to exhibit a broad peak. Also, when κ → 0,

Υ [Qqm + μpm, μpm]

Qqm

→ exp 
−π μ2 p2

m



. (102)

Thus, some cancellation occurs with the exp [−π p2

m (1 − μ2)] term, resulting in a

exp (−πp2

m) = exp (−πη2Z/α2m) dependence around q = 0. Unfortunately, these

general observations do not assist in the particular evaluation of the Fourier kernel. Rather, a triple integral must be evaluated numerically, including κ and ηZ

integrations out to infinity.

Press et al. (1992) discuss these types of improper integrals in their numerical

recipes book. Their method proposes dividing the integral into two parts, delimited

by a midpoint. In the lower section the integral is evaluated normally. In the upper

section a change of variables is used, wherein the inverse of the original variable is

29
introduced. Through this change, the limits of integration are modified to a finite

interval. Mathematically, this is written

∞

0

f(x) dx =



xm

0

f(x) dx +

∞

xm

f(x) dx

=



xm

0

f(x) dx +

x−1

m

0

f



y−1



y−2 dy (103)

where y = x−1. But this analysis only considered one-dimensional integration.

Since the triple integral of ˆˆI involves integration to infinity in two dimensions,

to use separate divisions of ηZ and κ space to perform this integration would

result in four different integration regions. The problem with this approach arises

in terms of the resonance regions in the integrand. Because these regions tend

to emanate from the origin, they become difficult to resolve accurately at high

argument values. At distances further from the origin the nonzero region tends

to become closer to the resonant line itself, meaning that only a small portion of

the overall integration region is actually responsible for the integral results. The

remainder of the area has an effective integrand value of zero. Using rectangular

integration regions, much larger numbers of samples are required to accurately

capture the behavior of the function at high argument values.

To overcome this difficulty, variables p and q, which already represent polar magnitudes of ηZ and κ variables, are combined into yet another polar variable, ρ. To

permit this representation it first seemed logical to modify the integration space

from the ηZ and κ variables to the dimensionless p and q variables:



∞

−∞

f1 (κ, ηZ) κ dκ ηZ dηZ = Z−2



∞

−∞

f1


 q

Z , p

p dp q dq. (104)

Introducing ωr = r/Z, the arguments of U become

U(κr, ηZω, φ − a) = U(qωr, pω, φ − a). (105)

Thus,

I = 2π

Z2

∞

0

q dq ∞

0

p dp 

2π

0

dφ ˆˆI [p, q, cos(φ)] J0 [2πU] . (106)

Introducing polar coordinates ρ and θ, such that p = ρ cos(θ), q = ρ sin(θ), we

can rewrite the integration with respect to p and q as,

∞

0

q dq ∞

0

p dp... =

∞

0

ρ3 dρ

π/2



0

1

2 sin(2θ) dθ... (107)

Having obtained this result, we can now make a few observations. First, let us consider the integration kernel on the RHS of equation (105). Due to the ρ3 weighting

30
we do not expect to find the most significant wavenumber response near the origin

nor along either the +p or +q axes. Rather, we expect the most significant region

to lie along the diagonal p = q. This consideration will to some degree dictate the

means used as we attempt to numerically evaluate the resulting integral. Second,

the use of the q variable implies a similarity of results measured at different scales.

Since off-axis variables are all expressed relative to the overall path length, we can

study the effects of scattering of fogs in the laboratory and rescale the lengths to

apply the measured results to large scale problems.

The computer procedure for evaluating the remaining (ρ, θ) representation of the

integral divides the quarter plane integration region into sections according to figure 10. The first region, labeled “1” in this figure, covers the full range of θ variation. The limit of ρ for this region is chosen such that ˆˆI has been significantly

reduced below its value at ρ = 0. However, because of the ρ3 dependence, the

remainder of the quadrant must also be evaluated using type 2 and type 3 regions

as labeled in figure 10. The latter region types are created by subdividing the θ

integral into a series of segments, ∆θ. In each of these segments successive integration regions of type 2 are evaluated, beginning at the boundary of the type

1 integration region, and continuing outward until the segment integration result

falls below some threshold. Once this threshold is crossed, the remainder of the

∆θ interval can be evaluated, from ρm (the maximum ρ value for this ∆θ segment) out to infinity, via introduction of the variable ξ = ρ−1. This final integration avoids the problems associated with integrating over both ηZ and κ because a

single diagonal region is considered.

q

1

2

3

∆ρ

∆θ

p

Figure 10. Division of the quarter plane in p and q, using regions delineated by constant surfaces in the

cylindrical coordinates ρ and θ.

The practical reason for dividing the integration task between a single type 1 region and many type 2 regions follows from the nature of the integrand being eval31
uated. Because many of the diagonals will make little or no contribution to the

overall integral, it is possible along many of these ∆θ regions to quickly move

from making type 2 calculations to making a single type 3 calculation. However,

there will probably be one or two diagonals that have significant contributions

out to large ρ values. Then, most of the computation time will be focused on the

few diagonals where more of the non-negligible integrand lies. This optimizes the

computational accuracy and speed of the technique.

32
3. Retroreflector Interactions

In the previous chapter the propagation of forward scattered radiation to a retroreflector was considered. In this chapter the interaction between the incident laser

beam and the retroreflector is considered. In the analysis of this interaction a specific retroreflective material, similar to the materials used in roadside signs, is considered. This material is composed of a series of miniature retroreflectors (minireflectors) embedded in a plastic sheathing. This retroreflective material can be

used in a standalone role or may be fashioned into a larger retroreflector that combines the microproperties of the retroreflective material with the properties of a

macro retroreflector.

At the material level the retroreflective material (retro-material) is composed of

a hexagonal grid structure of adjacent mini-reflectors similar to that shown in

figure 11.

Honeycombed

Retroreflectors

Interreflector

Gaps

Figure 11. Set of gridded mini-reflectors as they appear on reflective material.

In this material the mini-reflectors are embedded within a layer of plastic that

coats the entire surface of the material. The mini-reflectors appear to be circular

in cross-section and have some separation distance between adjacent retroreflecting regions. Thus, there is some fraction of the reflective material surface over

33
which the incident signal will not be retroreflected. Tests of retro material samples also indicate that there is a significant specularly reflected component when

the material is illuminated with a handheld laser pointer. Lastly, there is a substantial retroreflective component, but it is apparent that the combined effects of

the plastic sheathing and the inter-reflector gaps may reduce the total amount of

energy retroreflected from the material to below 50 percent.

The properties of the retroreflective regions can be understood by assuming that

each individual mini-reflecting region is a circular entrance aperture opening onto

three orthogonally oriented mirrored surfaces. The geometry of one of these reflectors is seen in figure 12. This figure shows views of the three orthogonal mirrored surfaces from both the front and the side.

A

B

O

r ́

mirror

surfaces

A

q ́ r ́

P

B

O

s ́

Figure 12. Geometry of a single mini-reflector.

In the figure, variables q, r, and s are introduced: r is the radius of an individual mini-reflector as measured in the main surface plane of the retro-material; s

represents the depth of the mini-reflector; and q measures the diagonal distance

between points P and A, representing the seam between adjacent mirrored sections

of the mini-reflector.

Variables q and s are related to r through the relations, q = r



3/2 and s =

r



1/2. These relations are found by computing the distance between points A

and B. The angle AOB has an arc measuring 120◦. The two legs, OA and OB,

in triangle AOB are both of length r

√ . Thus the distance AB can be computed as

3 r. The right triangle APB has two legs AP and BP of length q and the third leg

of length AB. The q value immediately follows, and s can be evaluated through

computations made on triangle AOP.

While complicated analyses would likely produce the same result, simple peering

through a large scale retroreflector reveals that beyond the reflection point, one

sees an image of the entrance aperture a distance 2s “into” the retroreflector.

Figure 13 illustrates an equivalent model of this situation. The fraction of light that

exits the system following passage through the entrance pupil at a given angle of

incidence is modeled as though the light were passing through a circular aperture

of radius r

, through a short tunnel of length 2s, and out a second circular aperture

34
Figure 13. Effective geometry of a retroreflector considering the fraction of incident energy that returns

in the retro-direction.

of radius r. In effect, the mirror properties of the retroreflector produce an image

of the entrance aperture a distance s beyond the deepest point (P) within the

retroreflector. There will thus also be a maximum angle for which a portion of

the light can enter the system at the edge of the entrance aperture and reach the

center of the retroreflector. All light entering at a greater angles of incidence will

not strike all three reflective surfaces and cannot retroreflect.

Mathematically, the portion of the energy which returns in the retro-direction is

given as the fractional overlap area between two circles of equal diameter. This

fraction is given by the formula,

fV (u) = 2

π




cos−1(u) − u

√

1 − u2



cyl(u/2), (108)

where fV is the vignetting factor, u = tan(θ

)/

√2, and cyl(u) is Gaskill’s cylinder

function:

cyl(u) =  1, if u < 1/2;

0, if u > 1/2. (109)

Thus fV (u) is an average function over the entire entrance aperture of each

mini-reflector. But because the retro-material contains many thousands of minireflectors, rather than compute individual lines of sight through the retro-material

geometry, it is more reasonable to apply fV (u) as a statistical quantity for any

radiation reaching the retro-material plane at z = Z. In effect, then, fV (u) acts to

reduce the reflected radiance depending on the angle of incidence of the received

radiance relative to a normal vector to the plane of the retro-material. And because

fV (u) has a multiplying effect on the radiance magnitude, upon Fourier transformation, the result is a convolution operation in angular frequency that is applied

to the return beam. This angular filtering also tends to reduce the importance of

large angle scattering because energy arriving from large angles is reduce greatly

in magnitude on retroreflection.

In addition to this magnitude modulation as a function of incidence angle, a spatial

truncation occurs due to the finite size of the retro-material. We assume the retro35
material to be circular in shape, with a diameter DR, and oriented perpendicular

to the z axis. Let the center of this reflective material be located off axis at point

rR in the z = Z plane. Due to scenario considerations we will place this off axis

position somewhere along the x axis. Mathematically, this reflector region could

be represented by another filter function,

fW (r) = cyl %

r − rR

DR

&

. (110)

However, from a practical consideration, the size of the reflector region will likely

be small. Thus, relative to the overall spread in the beam itself and the distance

off axis, we will consider the directional and positional properties of the incident

beam to be approximately constant over the surface of the retroreflector. We can

then just as easily compute the reflected radiation via use of a δ function and an

equivalent reflector area.

Finally, diffraction effects at the retroreflector will also influence the angular distribution of the retroreflected energy. While the total retroreflector area has a diameter DR, the retro-material is actually composed of thousands of mini-reflectors

with effective aperture radii of only rm, as in figure 12. We can model the diffraction spread of the return beam due to the area of these mini-reflectors as an angular

point spread function. However, crude experiments indicate the actual beamspread

may be 1.5–3 times larger than the diffraction limited results. In frequency space

this point spread function results in an additional modulation transfer function

applied to the return beam.

These influences of the retroreflective surface in spreading the retroreflected beam

are advantageous since the system envisioned would locate the observation camera at a distance offaxis from the laser provideing the source illumination. Some

beam spread thus increases the amount of energy that can be received by the offaxis camera. Insufficient spreading would result in no received energy, while too

much spread would likely result in further blurring of the received image.

Combining these considerations, it may be possible to approximately represent

the return beam as a Gaussian beam, with the return path adjusted in such a way

that the retroreflector appears on axis and the reflected “Gaussian” beam appears

to be travelling along the x = y = 0 axis. However, the receiver camera would

need to be offset from this new optical axis. To handle these considerations we

would need to assess the direction associated with the main retroreflecting beam

centroid, call it ω R. Let rX be the offset position of the receiver at the z = 0 plane

from the source laser at rL = 0. Let the main direction of incident energy arriving

at the retroreflector position be ω R. Then the direction of the main radiance reflecting back from the retroreflector along the return path would be travelling in the

direction associated with the triplet (−ω R, −1). Projecting this main propagation

centroid back to the z = 0 plane we would find an intercept at the point

rC = rR − Z ωR. (111)

36
Let us then define a new coordinate system in which we relocate rC onto the main

optical axis x = y = 0. In this new system the retroreflected energy appears to

be travelling similar to the original beam. The only changes that are needed are

to assess the characteristics of the new beam and to displace the position of the

receiver:

r

X = rX − rC. (112)

Having made this change in geometry all that remains in describing the return

radiation field is to examine the width of the Gaussian beam. We then consider

the propagation problem along the return path in chapter 4.

To describe the additional spreading of the return beam due to the retroreflector,

we begin by assuming we can characterize the incident beam’s angular structure

as a Gaussian about some mean direction ω R. We similarly assume we can characterize the nondiffraction-limited effects of the retroreflector aperture as a Gaussian

as well. These effects are then combined to generate a third Gaussian form for the

output from the retroreflective material about the mean retro-direction which we

have remapped to ω 

R = 0. Let v1 represent the angular width of the incident

radiation, and let v2 represent the spreading effect of the aperture on a monodirectional incident beam. Then we characterize the total output beam angular width

parameter, v, by

v2 = v2

1 + v22. (113)

In evaluating the reliability of this model we note that when the fog is thick, there

will be significant effects on the spreading of the beam with location and propagation direction on the return passage through the fog. It is therefore not critical

to exactly model the nature of the initial angular structure of the return beam. Of

much more critical interest is that we have captured the positional offset of the

retroreflector from the optical axis by the use of rR.

The remaining task is to evaluate v2. An approximation of this quantity can be

made based on the observation that there are approximately 60 mini-reflectors

per inch. We thus have rm ≈ 210 μm. For an aperture of diameter D the far

field diffraction pattern will have a zero at an angle of 1.22λ/D. The amplitude of

this pattern has a form called an Airy disc, or in Gaskill’s terminology, a sombrero

function. In either case, because the light passing through the mini-reflectors actually passes through two circular apertures, it is believed that the actual pattern for

the amplitude follows the square of a sombrero function. Therefore, the irradiance

pattern generated by each mini-reflector is modeled by the sombrero function to

the fourth power. This pattern may be approximated by a Gaussian function with

an angular width parameter of approximately v2 = 0.61λ/rm, where v2 appears

as a width parameter in the irradiance form: I(r

, z) ∝ exp[−2r2/(z2 v2

2)]. To

use this result we will also need to evaluate the angular properties of the incident

energy to determine a central peak in the propagation direction and to determine

the angular deviation that results in a reduction of exp(−2) from the central peak

value. The angular value of this result then becomes v1 and can be used to compute

v in equation 111.

37
4. Return Path Propagation

Figure 14 considers the overall geometry of the problem we are modeling. Though

we appear to have placed the retroreflector along the optical axis of the source,

this is not necessary. The key point is that the total path need not be considered

folded. Rather, the path can be viewed as three separate segments of one continuous path. The first portion leads from the laser source to the entrance pupil of the

retroreflector. The second stage (treated in chapter 3) passes the beam through the

retroreflector. The third segment includes the return path, which is a mirror image

of the outbound propagation path.

Retroreflector

Laser Source Fog Layer Fog Layer Wing Camera

Forward Path Return Path

zm Z – zm Z – zm zm

Figure 14. Geometry of the folded path problem where we treat the return path as a reversed version of

the forward propagation path.

In studying the return portion of the propagation path we note that the solution found for the outbound propagation path was radially symmetric. The return

beam will also share this symmetry property. In our problem, the laser source is

mounted somewhere on the fuselage or wing of an aircraft. Receiver cameras will

be mounted on each of the wings. Hence, if the laser and two cameras are mounted

along an axis of the wing, and the aircraft is flying level to the terrain generally in

the direction toward the landing field, it is likely that at least to a first approximation the cameras, laser, and retroreflector can all be considered to be lying within

a single plane. And somewhere along the line containing the laser source and the

cameras the return signal will be a maximum. This will simplify the subsequent

analysis since we will only need to look along a single line (call it the x axis) in

determining the position of the maximum return signal.

Based on the discussion in chapter 3, we will assume the receiver position has

already been adjusted to the location r

X. It should also be understood that in this

chapter we are going to be dealing with the return path. Nevertheless, to avoid

complications with notation we are not going to adopt a completely different system to denote the different axes of propagation used in analyzing the return propagation. It will have to be assumed that when using such terms as x, y, z, r, and ω .

38
in this chapter that we are nevertheless talking about a different problem than was

considered in either chapters 2 or 3. We will see that these problems are still here

to haunt us in chapter 5 where we treat radiation propagating in both hemispheres

at once, but we will leave that as a separate problem. For this chapter we simply

adopt the same terminology as was used in chapter 2, but apply it to the return

path.

We also mentioned at the close of chapter 3 that we neededto apply a factor v,

which represented the total effects of scattering on the outbound path and diffraction interaction with the retro-material. In chapter 2 we characterized the original

Gaussian beam emitted by the source laser expressed in radiance form in equation (72). For the return path we parameterize the return beam by,

IV (r, ω, Z) = 2PR

π v2 Z2 exp 

−2r2

v2 Z2



δ



ω − r

Z



, (114)

where PR is the total reflected power of the return beam. Note that if v1 was the

effective angular spread of the beam arriving at the retro-material, and RR was

the radiance arriving at the material in the peak direction of incidence, then the

arriving radiance as a function of angle would be

IR(ω, Z  ) ≈ RR exp 

−2

(ω − ω R)2

v2

1



, (115)

where we assume the incident irradiance is roughly constant over the surface of

the retroreflector. Integrating with angle and multiplying by the surface area of the

retroreflector, we have

PR ≈

π

2

v2

1 RR AE, (116)

where AE represents the effective area of the retroreflector. This factor may take

into account the foreshortening of the reflective surface if it is not perpendicular to

the incident beam and the effects of off-axis angle of arrival of the incident beam

to the retroreflective material.

Once the initial return beam is postulated, we may proceed to evaluate the net

returned beam passing through the scattering media similar to the method used in

chapter 2. The key difference here is that the fog aerosol will occupy the initial

portion of the path whereas previously it occupied the final segment of the path.

In solving for the return beam we first note that due to our change of geometry

we have the same basic problem to solve for both the vacuum and unscattered

problems. For completeness, we include these results here. First, for the vacuum

solution, given the slight changes involved with the introduction of v, we have

ˆˆIV

R(κ, η, Z) = PR exp 

−π2 v2 Z2 Λ2(Z)

2



. (117)

We then introduce a simplify Λ form,

Λ(  z) ≈ κ + η/Z = (ηZ + Zκ)

z = η0z , (118)

39
such that we can express this result as

ˆˆIV

R(κ, η, Z) = PR exp 

−π2 v2 η2

0

2



(119)

where zero range is now viewed as lying at the retroreflector plane.

The unscattered solution then follows simply as

ˆˆIU

R = ˆˆIVR exp [−τ (Z)] , τ (Z) = σ QZ. (120)

It then remains to evaluate the scattered solution. This form is slightly different

from the forward propagation problem because SR depends on aerosol that occupies the first portion of the optical path:

SR = exp







 σ 

Zf

0

Pˆf (|ηZ + κ [Z − z]|) dz







= $

M

m=1

exp Am  τ Υ [qm + μpm, (1 − Q)qm + μpm]

Qqm



. (121)

However, as expected, the resulting equation for SR generates the same results

when the aerosol occupies the full path (Q = 1).

The same techniques are used to evaluate the inverse transform of this result as

were discussed in chapter 2 for handling the forward propagation to the retroreflector plane.

40
5. Diffuse Radiation

In previous chapters we considered the forward-scattered radiation propagating

both from the laser source to the retroreflector and from the retroreflector back

to one of the wing cameras. In addition radiation may scatter from the fog volume itself and be detected by one of the wing cameras. This scattered radiation

appears as a noise source that must eventually be distinguished from the return energy propagating back from the retroreflector. Several strategies were mentioned

for removing this additional radiation. The simplest solution involves image processing to discriminate the high angular resolution return signal from the diffuse

backscattered radiance arising from scattering off the fog alone.

To determine the difficulty of the discrimination problem, the diffuse scattering radiance field must be modeled. The initial Gaussian beam forms the energy source

for this radiation field. However, unlike the forward scattering problem we need

significantly more information to perform this modeling task. To be able to model

the backscattered radiance field we will need to know details of the forward scattered field at all locations. Determining this information using the full resolution

phase function does not appear practical. Hence, we use a simplified version of

the forward scattered field. Second, we will ignore the small amount of energy

retroreflected from the corner reflector. Relative to the energy scattering within

the fog cloud this amount will be minimal.

5.1 Diffuse Radiation Scale Transform

To perform this analysis, the total diffuse radiation field will be divided into two

portions: (1) Idf flowing in the positive z direction (Ωz > 0); (2) Idb flowing in

the negative z direction (Ωz > 0). Previously we used ω to represent the x and y

components of the propagation vector in the positive z direction only. Here, we

will also use ω to denote propagation in the negative z direction but will use the

following convention:

(Ωx, Ωy, Ωz) =



 −ωx 

1 + ω2

x + ω2y

, −ωy



1 + ω2

x + ω2y

, −1



1 + ω2

x + ω2y





≈ (−ωx, −ωy, −1). (122)

Hence, Idf and Idb radiances assigned equal ω arguments represent energy propagation in opposite directions. This property will be useful when evaluating the

scattering integral in the RTE.

In addition to Idf flowing in the positive z direction, we will introduce the quantity

If (r, ω, z) representing energy that is propagating in nearly the forward direction.

41
Idf and Idb thus represent generalized diffuse fields with low angular resolution,

while the If field constitutes both unscattered energy from the laser source and

forward-scattered energy in a small solid angular region about the forward direction. To model If we use a considerably simpler RTE than was used in chapter

2.

ω ·

∂If

∂r

+

∂If

∂z + σ If = Af σIf , (123)

where the quantity Af represents the probability of scatter in the “forward” direction. This quantity can be viewed as one of the Am weighting factors used

in chapter 2. The remainder of the phase function can then be modeled using a

shorter Gaussian expansion. In particular, we will use the same expansion terms

as are contained in table 1, except that the first two terms have been removed. The

corresponding Af is thus the sum of the first two Am terms in table 1 and has the

value 0.4567. Because α1 and α2 are significantly greater than the terms, we are

justified in assuming a zero width for the forward scattering portion.

This form for the radiative transfer of the nondiffusely scattered radiation assumes

the forward scattered energy is scattered into a cone so small that even scattered

energy can be assumed to be propagating virtually along the path of the original

unscattered energy. To solve the If equation we can perform a scale transformation on the extinction coefficient (McKellar and Box, 1981). This approach involves the replacement of the original extinction coefficient by the reduced quantity:

σ' = (1 − Af) σ. (124)

Through this replacement the RTE for If is reduced to the unscattered RTE with

σ ̃ as the extinction coefficient.

ω ·

∂If

∂r

+

∂If

∂z + σ I ' f = 0, (125)

which has the solution,

If = IV exp [−τ

'(z)] , (126)

where IV is the vacuum solution introduced in equation (72), and

τ

'(z) = z

0

σ' ds. (127)

A scaling transformation also entails a change to the single scattering albedo in

which



( = (1 − Af ) 

(1 − Af)

, (128)

such that the product 

(σ' = (1 − Af )σ appears as a scattering coefficient (c.f.

Tofsted and O’Brien, 1998).

The next task is to write RTE’s for both the forward and backward oriented diffuse

radiance fields. To write these equation, however, we must recognize that scattering occurs between radiances flowing into opposite hemispheres. But since the

42
phase functions are generally written in terms of scattering angles measured by

∆ω, we must make allowances for scattering into opposite hemispheres through

the use of two different phase functions. The first,

Pdf (∆ω) = 

M

m=3

A

m α2m exp 

−α2

m ∆ω2



, (129)

represents scattering in the forward direction (within the same hemisphere into

which the scattering energy was originally travelling). The second,

Pdb(∆ω) = 

N

n=1

B

n β2n exp 

−β2

n ∆ω2



, (130)

represents the scattering of energy travelling in the direction of one hemisphere

into the opposite hemisphere. It both these cases the A

m and Bn quantities have

been adjusted from the values listed in table 1 by the following adjustments:

A

m = Am/(1 − Af ) and Bn = Bn/(1 − Af ). The reason for these adjustments is

that due to the scale transformation, the diffuse scattering terms represent 100 percent of the scattering in the scale transformed equations. But the sum of the scattering components () Am + ) Bn) now only accounts for (1 − Af ) probability

of scattering. To renormalize the phase functions we need to divide by this factor.

5.2 Dual Direction Radiance Equations

We can then write the coupled differential equations describing the propagating

forward and backward diffuse radiance fields as

ω ·

∂Idf

∂r

+

∂Idf

∂z + σ I ' df = 

( σ'



∞

−∞

Pdf (∆ω)



I

f + Idf

dω 

+ 

( σ'



∞

−∞

Pdb(∆ω) I

db dω 

, (131)

−ω ·

∂Idb

∂r − ∂Idb∂z + σ I ' db = 

( σ'



∞

−∞

Pdb(∆ω)



I

f + Idf

dω ,

+ 

( σ'



∞

−∞

Pdf (∆ω) I

db dω 

, (132)

where ∆ω = |ω − ω 

| as usual. The terms I

f , Idf , and Idb refer to radiances

propagating in directions associated with the vector ω 

. Note that the scattering of

radiance out of the opposite hemisphere is centered on the backscattering direction. Also note that the directional derivative signs in the Idb equation are switched

to represent the opposite direction of propagation for the backward propagating

43
radiance field. Evaluating both these diffuse radiance fields, and coupling the results of this analysis to the results of computations based on analysis contained

in chapters 2–4, we create a complete picture of the radiation flowing through the

fog cloud.

To solve this pair of equations, we first proceed as previously by taking the fourth

order Fourier transforms (in r and ω ) of both equations:



−κ · ∂

∂η

+

∂

∂z +  ̃σ



ˆˆIdf =  ̃σ ̃



Pˆdf *

ˆˆIdf + ˆˆIf

+

+ Pˆdb ˆˆIdb, (133)



+κ · ∂

∂η − ∂∂z +  ̃σ



ˆˆIdb =  ̃σ ̃



Pˆdb *

ˆˆIdf + ˆˆIf

+

+ Pˆdf

ˆˆIdb

. (134)

The physical geometry of the problem to be solved can thus be stated as follows:

We assume that the leading edge of the cloud begins at some distance from the

laser source. This distance was called Zf in chapters 2 and 4. Let us then assume

that the cloud ends at distance Zb. The length of the cloud is then ZC and at

distance Z embedded within the cloud is the retroreflector.

Since no forward scattering may occur before the forward edge of the cloud we

must have ˆˆIdf (Zf )=0. (135)

Likewise, there must be no backward scattering beyond the farthest edge of the

cloud: ˆˆIdb(Zb)=0. (136)

We may also attempt to simplify the system of equations by adding and subtracting

eqs. (131) and (132) from one another. Combining terms we obtain the equation

set



−κ · ∂

∂η

+

∂

∂z 

ˆˆI1 +  ̃σ ˆˆI2 =  ̃σ ̃ PˆD

*

ˆˆI2 + ˆˆIf

+

, (137)



−κ · ∂

∂η

+

∂

∂z 

ˆˆI2 +  ̃σ ˆˆI1 =  ̃σ ̃ PˆS

*

ˆˆI1 + ˆˆIf

+

, (138)

where ˆˆI1 = ˆˆIdf + ˆˆIdb, ˆˆI2 = ˆˆIdf − ˆˆIdb, (139)

PˆS = Pˆdf + Pˆdb, PˆD = Pˆdf − Pˆdb. (140)

For this system we have boundary conditions,

ˆˆI1 = −ˆˆI2, z = Zf ; ˆˆI1 = +ˆˆI2, z = Zb. (141)

Given these boundary conditions we consider the general solution to this set

of equations. The first question we need to ask is whether the implicit solution

method introduced in chapter 2 can be used in solving this problem. It would

seem that the implicit method should work individually on each pair of equations.

44
Hence, we have two analyses similar to the development leading to eqs. (63), (64),

and (70). In particular, we generate the same relationship in both eqs. (135) and

(136) to describe the behavior of η as a function of κ and z:

η(z) = η0 − κz. (142)

The remaining relationships provide differential equations for the radiances in

terms of distance:

dˆˆI1

dz = −σ ̃

*

ˆˆI2 −  ̃ PˆD


ˆˆI2 + ˆˆIf+ , (143)

dˆˆI2

dz = −σ ̃

*

ˆˆI1 −  ̃ PˆS


ˆˆI1 + ˆˆIf+ . (144)

This system of equations can be solved through use of the Newton-Raphson method.

Solution involves determining a starting value for ˆˆI1 at Zf that when propagated

through to the far end of the cloud at range Zb produces the correct boundary condition listed in equation (139). Since there is no change in either variable when σ ̃

is zero, and σ ̃ is zero prior to the beginning of the cloud, the values of both ˆˆI1 and ˆˆI2 at zero range are also their values at the forward edge of the cloud.

The Newton-Raphson method can be used by setting the function we wish to find

a zero for as ˆˆI1− ˆˆI2 evaluated at Zb as a function of ˆˆI1 at Zf . The means of evaluating the two output parameters at the far end of the cloud involves implementation

of a numerical integration procedure. Press et al. (1992) provide listings of procedures to accomplish these tasks, including routines rtsafe to do a combination

Newton-Raphson and bisection procedure, odeint to integrate a set of ordinary

differential equations, and codes rkqs and rkck to do quality-controlled fourthorder Cash-Karp Runge-Kutta integration of the differential equations.

Once the Newton-Raphson approach has been used to evaluate ˆˆI1, we can evaluate

Idb at z = 0 using the same inverse transform techniques described in sections

2.7.5 and 2.8. The only difference is in the interpretation of the outputs: While

the r system will be the same, the ω sense of the incoming radiance field will be

the opposite of the sense of the outbound radiances of the same magnitude and

direction.

We recover the original diffuse fields through the equations

ˆˆIdf =

ˆˆI1 + ˆˆI2

2 , ˆˆIdb =

ˆˆI1 − ˆˆI2

2 . (145)

With this result we have shown the means of computing all the fields and quantities of interest. We leave the analysis of these results to a future publication.

45
6. Conclusions

This document develops a theoretical basis for the problem of dual path propagation through scattering aerosol advection fog. To solve this problem our approach

has been rigorous — we have considered both propagation through a fog layer to

a retroreflector and return propagation through the fog. We have also considered

the issue of scattering within the fog layer in both forward and backward directions. A recasting of the small angle approximation has also been accomplished

as well as a complete elucidation of the methods employed by Smirnov to solve

partial differential equations, and this approach has been applied to the propagation problem at hand. In addition, we have considered means useful in reducing

the inverse transform problem from a fourth order inverse Fourier transform process to a third order calculation involving Bessel functions. This method greatly

simplifies the inverse transform process since it guarantees that one will generate

a real-valued transform function as a result. We then considered the process of

retroreflection itself and provided a model useful in approximating the effects of

retroreflection from two commercially available retroreflective materials. Lastly,

we considered the calculation of diffuse radiation scattering within the fog cloud

volume. While complicated, these procedures provide an accurate means of evaluating the effects of propagation within fog aerosols to assess a navigation through

fog system.

46
References

Deirmendjian, D., Electromagnetic Scattering on Spherical Poly-dispersions, Elsevier, New York, 1969.

Gaskill, J. D., Linear Systems, Fourier Transforms, and Optics, J. Wiley & Sons,

New York, 1978.

Goodman, J. W., Introduction to Fourier Optics, McGraw-Hill, New York, 1968.

McKellar, B. H. J., and M. A. Box, “The scaling group of the radiative transfer

equation,” J. Atmos. Sci., 38, p. 1063-1068, 1981.

Miller, A., Mie Code AGAUS 82, ASL-CR-83-0100-3, U.S. Army Atmospheric

Sciences Laboratory, White Sands Missile Range, NM, 1983.

Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical

Recipes in C, Cambridge Univ. Press, Cambridge, England, 1992.

Shettle, E. P., and R. W. Fenn, Models for the Aerosols of the Lower Atmosphere and the Effects of Humidity Variations on Their Optical Properties,

AFGL-TR-79-0214, U.S. Air Force Geophysics Laboratory, Hanscom Air Force

Base, MA, 1979.

Shirkey, R. C., L. D. Duncan, and F. E. Niles, EOSAEL 87, Executive Summary,

ASL-TR-0221-1, U.S. Army Atmospheric Sciences Laboratory, White Sands Missile Range, NM, 1987.

Siegman, A. E., Lasers, University Science Books, Mill Valley, CA 94941, 1986.

Smirnov, V. I., A Course of Higher Mathematics, Volume II, Pergamon

Press/Addison-Wesley, Reading, MA, 1964.

Tam, W. G., and A. Zardecki, “Multiple Scattering Corrections to the Beer-Lambert

Law. 1: Open Detector,” Appl. Opt., 21, 2405-2412, 1980.

Tofsted, D. H., and S. G. O’Brien, Three-Dimensional Radiative Transfer Modeling of Tropospheric Atmospheres, ARL-TR-1629, U.S. Army Research Laboratory, White Sands Missile Range, NM, 1998.

Tofsted, D. H., B. T. Davis, A. E. Wetmore, J. Fitzgerrel, R. C. Shirkey, and R. A.

Sutherland, EOSAEL92, Aerosol Phase Function Data Base PFNDAT, ARL-TR273-9, U.S. Army Research Laboratory, White Sands Missile Range, NM, 1997.

van de Hulst, H. C., Multiple Light Scattering: Tables, Formulas, and Applications, Volume 1, Academic Press, New York, NY, 1980a.

van de Hulst, H. C., Multiple Light Scattering: Tables, Formulas, and Applications, Volume 2, Academic Press, New York, NY, 1980b.

47
Watkins, W. R., D. H. Tofsted, V. G. CuQlock-Knopp, J. B. Jordan, and J. O.

Merritt, Navigation through fog using stereoscopic active imaging, Proceedings

of SPIE 4023: Enhanced and Synthetic Vision 2000, pp. 20-28, 2000.

Zeng, X., C. Liang, and Y. An, “Far-field propagation of an off-axis Gaussian

wave,” Appl. Opt., 38, p. 6253-6256, 1999.

Zardecki, A., and S. A. W. Gerstl, “Multi-Gaussian Phase Function Model for

Off-Axis Laser-Beam Scattering,” J. Opt. Soc. Am. A, 3, p. 119, 1986.

Zardecki, A., and W. G. Tam, “Multiple Scattering Corrections to the Beer-Lambert

Law. 2: Detector with a Variable Field of View,” Appl. Opt., 21, 2413-2420, 1982.

48
49

REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and 

maintaining the data needed, and completing and reviewing the collection information. Send comments regarding this burden estimate or any other aspect of this collection of information, 

including suggestions for reducing the burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson 

Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to 

comply with a collection of information if it does not display a currently valid OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

September 2002 

2. REPORT TYPE 

FINAL 

3. DATES COVERED (From - To) 

1 OCT 1999 - 30 DEC 2000 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

4. TITLE AND SUBTITLE 

Navigation Through Fog: A Mathematical Analysis of Folder Path Propagation 

5c. PROGRAM ELEMENT NUMBER 

5d. PROJECT NUMBER 

5e. TASK NUMBER 

6. AUTHOR(S) 

David H. Tofsted 

 5f. WORK UNIT NUMBER 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

U.S. Army Research Laboratory 

Information Science and Technology Directorate 

Battlefield Environment Division (ATTN: AMSRL-CI-EW) 

White Sands Missile Range, NM 88002-5501 

8. PERFORMING ORGANIZATION 

 REPORT NUMBER 

ARL-TR-2648 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR'S ACRONYM(S) 

U.S. Army Research Laboratory 

2800 Powder Mill Road 

Adelphi, MD 20783-1145

11. SPONSOR/MONITOR'S REPORT 

 NUMBER(S) 

ARL-TR-2648 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

Approved for public release; distribution unlimited. 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

Wendall Watkins of the U.S. Army Research Laboratory has a currently pending patent, which provides a technique for solving the problem of 

navigation through fog. This technique entails the use of laser illumination of retroreflective targets on airfield runways, stereo detection of the 

return signal, and deblurring processing of the received imagery. The current analysis attempts a mathematical analysis of this same problem. 

Here, the propagation equations necessary to model the propagation of a Gaussian beamwave through a forward scattering fog aerosol are 

developed. Interaction of the propagated wave with a retroreflective material is postulated, and the returned reflected energy is evaluated. Also, 

backscattered diffuse energy is analyzed. Through this analysis, it should be possible to predict the angular structure of energy returning from a 

retroreflector embedded within the fog field. It should also be possible to predict the level of diffuse radiation which must be dealt with in order to 

discriminate between the directly reflected radiation and the diffusely scattered radiation, both in the forward and backward hemispheres. 

15. SUBJECT TERMS

Aerosols, Scattering Phase Functions, Fog, Radiative Transfer, Propagation, Detection 

16. SECURITY CLASSIFICATION OF: 

19a. NAME OF RESPONSIBLE PERSON

David H. Tofsted 

a. REPORT 

UNCLASSIFIED

b. ABSTRACT 

UNCLASSIFIED

c. THIS PAGE 

UNCLASSIFIED

17. LIMITATION 

OF 

 ABSTRACT 

SAR 

18. 

NUMBER 

 OF PAGES 

 59 19b. TELEPHONE NUMBER (Include area code) 

(505) 678-3039 

Standard Form 298 (Rev. 8/98) 

Prescribed by ANSI Std. Z39.18