Appendix G

Using the Sommerfeld–Watson

Transformation

G.1 Introduction
The complex spectral number summation technique already has been
described briefly in Section 3.7. Here, for illustrative purposes, we outline the
technique for computing the scattering from a perfectly reflecting sphere and
from a transparent refracting sphere. The technique offers efficient convergence
provided that (1) a complete set of poles can be isolated, (2) their residues can
be calculated accurately, and (3) the set of poles offers a converging solution.
By use of contour integration, we replace the sum for the scattering along the
real axis in real-integer spectral number space with another sum in complex
spectral number space (see Fig. 3-9). From Eq. (3.12-3), the radial component
of the scattered wave in series form is
∞
E
Er( j ) = − 2o
x
∑ Sl( j ) (2l + 1)il −1ξl+ ( x ) Pl1(cosθ ) (G-1)
l =1

Here the superscript “ ( j ) ” denotes the j th degree scattered wave and Sl( j ) is
its scattering coefficient. Eq. (3.9-3) gives (
Sl( 0 ) = − 1 + Wl− / Wl ) / 2 ,
2 + 3
Sl(1)= n( −2i / Wl ) / 2 , Sl( 2 )
= 2 nWl / (Wl ) , etc. For the transverse component,
appropriately simplified per discussion after Eq. (3.12-4), a similar expression
is obtained:
∞
Eo 2l + 1 ′ ′
Eθ( j ) = ∑ Sl( j ) l(l + 1) il −1ξl+ ( x )Pl1 (cosθ )sin θ (G-2)
x l =1

587
588 Appendix G

Let us consider just the radial component in Eq. (G-1).

We define the function

F[l ] = Sl( j )i l (2l + 1)ξl+ ( x ) Pl1 (cos(θ + π )) (G-3)

Consider the contour integral I + in the complex plane along a closed path
enclosing the positive real axis and lying at an infinitesimal distance ε above
and below the real axis [see Fig. 3-9(a)]:

1 ∞− iε F[l ′ − 1 / 2] 0 + iε F[l ′ − 1 / 2] 
I+ = ∫0 −iε dl ′ + ∫∞+ iε dl ′  (G-4)
2  cos(πl ′ ) cos(πl ′ ) 
We can use the theory of residues in contour integration of analytic functions to
evaluate this integral by summing the residues at the simple poles of the
integrand, which on the real axis are located only at the half-integer points in l' .
Upon summing the residues from the complex contour integrations around each
of the poles in Eq. (G-4), one obtains
∞ ∞
I + = −i ∑ ( −1)κ F[ k ] = ∑ (−1)l −1iF[l]
κ =0 l =1
∞
= ∑ Sl( j )il +1(2l + 1)ξl+ ( x )(−1)l −1 Pl1(cos(θ + π )) (G-5)
l =1
∞
x2
=− ∑ Sl( j )i l −1 (2l + 1)ξl+ ( x ) Pl1 (cos θ ) = Er( s )
l =1
Eo

To obtain the result in Eq. (G-5), note that Pl1 (cos θ ) = ( −1)l −1 Pl1 (cos (θ + π )) ,
and also that F[0] = 0 because P01 (cos θ ) ≡ 0 . A similar relation holds for Eθ( s ) .
In this manner, we convert the sum in Eq. (G-1) into the contour integral given
in Eq. (G-4).
Next, we need to extend the contour integration in Eq. (G-4) to include the
entire real axis so that we can use the vanishing property of the integrand when
it is evaluated on a semi-circular arc bounding the upper complex plane whose
radius approaches infinity. Therefore, we will need to invoke the symmetry
properties of the integrand in this contour integral. From the defining
differential equation for the spherical Hankel functions of order l − 1 / 2 [see
Eq. (3.6-2)], it follows that these functions must be either symmetric or anti-
symmetric in l . It is easily shown from the defining equation for the Bessel
function that

ξl±−1/ 2 ( x ) = ( −1)l ξ−±l −1/ 2 ( x ) (G-6)

Using the Sommerfeld–Watson Transformation 589

It follows from Eq. (3.5-11) that all of the Wronskian scattering terms,
Wl±−1/ 2 and Wl −1/ 2 , each of which involves products of a pair of spherical
Hankel functions, are symmetric in l . Therefore, it follows from Eq. (3.5-11)
that all scattering coefficients of spectral number l − 1 / 2 , Sl(−1
0) (1)
/ 2 , Sl −1/ 2 ,
Sl(−21) / 2 , ⋅⋅⋅ , are symmetric in l . Next, from the defining differential equation for
the Legendre functions, we have

Pl1−1/ 2 (cos θ ) = P−1l −1/ 2 (cos θ ) (G-7)

Also, P01 (cos θ ) = P−11 (cos θ ) ≡ 0 . Noting that 2(l − 1 / 2) + 1 = 2l , it follows that

( −2l )i − l −1/ 2 = ( −1)l +1 (2l )i l −1/ 2 (G-8)

Assembling all of the parts, we see that

F[ −l − 1 / 2] F[l − 1 / 2]
=− , F[0] = F[ −1] = 0 (G-9)
cos( −πl ) cos(πl )

Now, the contour integral I − enclosing the negative real axis is to be taken in
the same counter-clockwise sense that also applied to I + in Eq. (G-4).
Therefore, we obtain for I −

1 −∞+ iε F[l ′ − 1 / 2] 0 − iε F[l ′ − 1 / 2] 

I− = ∫0 + iε dl ′ + ∫−∞−iε dl ′  (G-10)
2  cos(πl ′ ) cos(πl ′ ) 

Upon making a change of integration variable to v = −l ′ and applying the anti-

symmetry property in Eq. (G-9), the resulting contour integrals in Eq. (G-10)
become identical with those in Eq. (G-4). Hence, I − = I + . It follows that we
may extend our contour integration in Eq. (G-4) to enclose the entire real axis.
We obtain

Eo  ∞− iε F[l − 1 / 2] −∞+ iε F[l − 1 / 2] 

Er( s ) = ∫−∞−iε dl + ∫∞+ iε dl  (G-11)
4 x 2  cos(πl ) cos(πl ) 

The same symmetry holds for Eθ( s ) .

It is convenient to change from the spectral number variable l in
Eq. (G-11) to the argument of the Airy functions y because we will be using
the asymptotic forms for the Hankel functions in terms of Airy functions and
we will be concerned with the zeros of certain combinations of these functions.
The defining relationships between y and l are given by Eqs. (3.8-2) and
(3.8-3) (here l = ν − 1 / 2 ), which can be greatly simplified in our case where
590 Appendix G

xo >> 1 to the near-linear relationships given in Eq. (3.13-15). In this case,

Eq. (G-11) becomes

Eo K x o  ∞− iε F[l( y) − 1 / 2] −∞+ iε F[l( y) − 1 / 2] 

Er( s ) = ∫

4 x 2  −∞− iε cos(πl( y))
dy + ∫∞+ iε cos(πl( y))
[ ]
dy + O xo−1  

l = xo + yK x o 

(G-12)

and similarly for Eθ( s ) .

Next, we deform this closed contour encompassing the entire real axis for
y into one that excludes the real axis as interior points but encompasses the
upper complex plane [see Fig. 3-9(b)]. For Hankel functions of the first kind, it
can be shown that the integral along the outer boundary of this deformed
contour, whose radius approaches infinity, is zero. Thus, the integral around the
closed path of the deformed contour may be evaluated in terms of the residues
at the poles of the integrand located anywhere in the complex plane except on
the real axis. The potential advantage of this approach may be a much more
rapidly converging series in the complex plane than the original series along the
positive real axis.
To obtain expressions for Er( s ) and Eθ( s ) , we must find the poles of
F[l( y) − 1 / 2] / cos(πl( y)) other than those on the real axis and sum up the
residues at these poles. Let yκ , κ = 1, 2,L define the location of those poles of
F[l( y) − 1 / 2] / cos(πl( y)) in the complex plane away from the real axis. From
Eqs. (3.5-11) and (G-3), we see that these poles occur at the zeros of Wl −1/ 2 .
Moreover, Eq. (3.5-11) also shows that when F[l − 1 / 2] represents the
scattering coefficients for the j th -degree scattered wave, it has poles of order
j + 1 at these zeros. Thus, F[l − 1 / 2] for an externally reflected wave has only
simple poles; F[l − 1 / 2] for the primary wave, which is refracted twice as it
passes through the sphere without any internal reflections, has second-order
poles; and so on.
Let us write
F[l( y) − 1 / 2] U[ y]
= (G-13)
cos(πl( y)) V [ y]

where U[ y] has no poles for Im[ y] > 0 and is not zero at the zeros of V [ y],
which occur at the points y = yκ , κ = 1, 2,L. If F[l − 1 / 2] has poles of order m
at these points, then V [ yκ ] = V ′[ yκ ] = L = V ( m −1) [ yκ ] = 0 , but V ( m ) [ yκ ] ≠ 0 .
From the theory of analytic functions, it follows that the residue for U[ y] V [ y]
at a pole of order m is given by
Using the Sommerfeld–Watson Transformation 591

U[ y ] mU ( m −1)
∫ V [ y]
dy = 2πi
V (m)
(G-14)
y = yκ y = yκ

G.2 Application to a Perfectly Reflecting Sphere

Let us first consider the case of the externally reflected wave from a large,
perfectly reflecting sphere, which has only simple poles. From residue theory, it
follows from Eqs. (G-12) and (G-14) that the radial component of the scattered
field is given by
∞
iπEo U[ yκ ]
Er( s ) = 2 K xo ∑ (G-15)
2x κ =1
V ′[ yκ ]

Here the sequence yκ , κ = 1, 2,L yκ , κ = 1, 2,L defines the locations in the

complex plane where V [ y] is zero.
Recalling the asymptotic forms for the scattering coefficients, we have
from Eqs. (3.9-3) and (3.17-1)
−
1  Wl −1/ 2  ψ ′ (x ) Ai ′[ y]
− Sl(−01) / 2 =  1 +  → l +−1/ 2 o → (G-16)
2  Wl −1/ 2  n →∞ ξl′−1/ 2 ( xo ) Ai ′[ y] − i Bi ′[ y]

where y and l are given by Eq. (3.13-15) with ν replaced by l + 1 / 2 , i.e.,

y = K x−1o (l + 1 / 2 − xo ) + L and l + 1 / 2 = xo + K x o y + L . It follows that U and
V are given by

U ( y) = −2i µ lξµ+ ( x ) Pµ1 (cos(θ + π )) Ai ′[ y] / cos(πl ) 


V ( y) = Ai ′[ y] − i Bi ′[ y]  (G-17)

µ = l −1/ 2 
For xo >> 1, the asymptotic forms for the Hankel functions in terms of Airy
functions [see Eq. (3.8-1)] apply. Also, for an observer some distance from the
reflecting sphere, x >> xo >> 1, the asymptotic forms for the Airy functions in
terms of complex exponential functions apply [see Eq. (3.8-10)]. In this case,
1/ 2
 x  π  1
i µ ξµ+ ( x ) →   exp i Dl + lθ l −   , µ = l − (G-18)
x > x o >>1  Dl   2   2

where Dl and θ l , now complex, are still defined in Eq. (3.10-3) (see Fig. 3-14).
Similarly, the product Pµ1 [cos(θ + π )] / cos(lπ ) has the asymptotic form
592 Appendix G

Pµ1 [cos(θ + π )] 2 µ  1/ 2 cos[l(θ + π ) + π / 4]

→  , µ = l −1/ 2 (G-19)
cos(lπ ) µ >>1  π sin θ  cos(lπ )

The zeros of [ Ai ′[ yˆ ] − i Bi ′[ yˆ ]] only lie in the first quadrant of the complex

plane where Re[ y] > 0 and Im[ y] > 0 . I t f o l l o w s t h a t
Re[l ] = xo + Re[ y]K x o >> 1 and Im[l ] = Im[ y]K x o >> 1 . Therefore,
exp(iπl ) → 0 , and we can simplify Eq. (G-19) to

Pµ1 [cos(θ + π )] 2 µ  1/ 2  π  1
→ exp −i lθ +   , µ = l − (G-20)
cos(lπ )  π sin θ    4  2

It follows that
1/ 2
 x 3 sin 3 θ l 
U ( y ) → 4  exp[i( Dl + l(θ l − θ ))] Ai ′[ y] (G-21)
| y |>> ~ 2 / K  2πi cos θ l sin θ 
x >> x o >>1

When this asymptotic form for U ( y) plus V ( y) is used in Eq. (G-11), it is

easily shown that the result is virtually identical to the scattering integral for
Er( So ) using the stationary phase approach as given in Eq. (3.12-5) (with the
phasor exp(iΦ + ) deleted). Thus, we could have deformed the scattering
integral in Eq. (3.12-5) directly into one that spanned the upper complex plane
to arrive at Eq. (G-15) without enduring the foregoing discussion leading up to
Eq. (G-15).
To evaluate Eq. (G-15), we need the zero crossings of ( Ai ′[ yˆ ] − i Bi ′[ yˆ ]) ,
which are exhibited in Fig. G-1. They lie only along the straight line defined by
y = β exp(iπ / 3) in the positive half of the complex plane, and their values,
given by βκ , κ = 1, 2 ⋅⋅⋅, are shown in this figure. Because of the high
sensitivity of l and µ for large xo to changes in y ,
| ∂l ∂y | = | ∂µ ∂y | =˙ K x o ≈ 475 , high precision is required in determining the
values of yκ in order to obtain accurate values for the phasor part of U[ yκ ] .
Asymptotic formulas for obtaining the roots of Ai ′[ y] − i Bi ′[ y] = 0 are found in
[1]. The values − βκ , κ = 1, 2 ⋅⋅⋅ are also the zeros of Ai ′[ y] along the negative
real axis [see Eq. (G-25)]. From the defining differential equation for the Airy
function, z ′′ = xz , it follows that to calculate the residue at the pole we have
V ′( y) = y( Ai[ y] − i Bi[ y]) . We have, therefore, all of the parts required to
calculate the electric field of the reflected wave at the low Earth orbiting (LEO)
satellite as a function of θ using the theory of residues and Eq. (G-15).
Using the Sommerfeld–Watson Transformation 593

In order for the summed series in Eq. (G-15) to converge practicably,

U ( yκ ) must approach zero with increasing κ in an efficient way. Inspection of
Eq. (G-21) shows that this convergence question hinges on the behavior of the
imaginary part of the complex phase term Φ l = Dl + l(θ l − θ ) . By expanding
Φ l in powers of (l − xo ) , we obtain

( K x y )2
o
Φ l = Φ x o + (θ o − θ ) K x o y+ (G-22)
2 Do

A sign change in d (Im[Φ l ]) / d (Im[ y]) occurs at

Im[ y] = Dx o (θ − θ o ) / K x o ≈ 200(θ − θ o ) , where (θ − θ o ) is expressed in
milliradians. Therefore, when the angular position of the LEO is above the
geometric shadow boundary, i.e., (θ − θ o ) > 0 , ever so slightly, one obtains a
very slowly converging series, which is impracticable. On the other hand, for
(θ − θ o ) ≤ 0 , the series converges rapidly. The latter is, of course, the shadow
region where no stationary phase points exist. Thus, two methods for summing
the spectral series, the stationary phase technique and the contour integration
technique in the complex spectral number plane, complement each other to
some extent.

G.3 Application to a Refracting Sphere

The extension to a refracting sphere with a finite index of refraction is
straightforward. Here one replaces the derivatives of the Airy functions shown
in Eqs. (G-17) and (G-21) with the corresponding Wronskian forms from
Eq. (3.5-11). For example, for the external reflected component, j = 0 , and one
obtains
1/ 2 
 x 3 sin 3 θ l 
U ( y) =  
 2πi cos θ l sin θ 
[( )](
exp i Dl + l(θ l − θ ) Wl− + Wl ) 

 (G-23a)
1 
V ( y ) = Wl , l + = x o + K x o y + L 
2
For the refracted wave passing through the sphere without internal reflections,
j = 1 , and U ( y) and V ( y) become
594 Appendix G

1/ 2 
 x 3 sin 3 θ l 
U ( y ) = 4 
 2πi cos θ l sin θ 
[(
exp i Dl + l(θ l − θ ) 
 )]
 (G-23b)
2 1 
V ( y) = (Wl ) , l + = xo + K x o y + L 
2

which has poles of order two at the zeros of Wl .

In the case of either j = 0 or j = 1 , we need to isolate the zeros of Wl , as
defined in Eq. (3.5-11). Using the Airy function asymptotic forms for the
Hankel functions, it follows from Eqs. (3.8-1)–(3.8-4), that the Wronskian
scattering forms in Eq. (3.5-11) are given by

Wl± = πn1/ 2 (n(Ai ′[ y] m i Bi ′[ y])(Ai[ yˆ ] m i Bi[ yˆ ]) 


− n −1 (Ai[ y] m i Bi[ y])(Ai ′[ yˆ ] m i Bi ′[ yˆ ]) + O[ N , ν −5 / 3 ] 
) 
1/ 2
Wl = πn (n(Ai ′[ y] − i Bi ′[ y])(Ai[ yˆ ] + i Bi[ yˆ ]) 

− n −1 (Ai[ y] − i Bi[ y])(Ai ′[ yˆ ] + i Bi ′[ yˆ ]) + O[ N , ν −5 / 3 ] 
)

ν
y = ν 2 / 3ζ   =˙ K x−o1 (ν − xo ) 
x
 o 

ν 
ˆy = ν 2 / 3ζ 
 nx  =˙ n
 o
−1/ 3 −1
( [
y − Nxo Knx o 1 + O (nxo ) −1
]) 


(G-24)

The zero points of Wl also are shown in Figure G-1 for refractivity values
ranging continuously from Nxo = 83π down to Nxo = π / 10 . As N increases,
( )
we note from Eq. (G-24) that Wl + Wl− / Wl → 2 Ai ′[ y] / ( Ai ′[ y] − i Bi ′[ y]) ;
therefore, the zero points of Wl should approach the zeros on the line
Im[ y] = tan[π / 3]Re[ y], the zeros for the perfectly reflecting case.
However, as N grows large, the form for Wl given in Eq. (G-24) loses
numerical precision. Although the magnitude of y may be small,
yˆ =˙ n −1/ 3 y − Nxo Knx
−1
o
can be large if Nxo >> K x o . Because y is complex,
Ai[ yˆ ] and Bi[ yˆ ] can grow very large. Equation (G-24), which involves
differences between these terms, can lose numerical precision. This can be
alleviated to a certain extent by noting several identities involving the Airy
functions in the complex plane [1]. These are given by
Using the Sommerfeld–Watson Transformation 595

k=
10
10

8 0]
]=
y
i '[ k
iB

Im [ y ]
6 ]- p
83 10
i' [
yk = p/
=
[A Nx o xo
4 N
k=2 W l ( yk) = 0
2 xo = 6378k
k=1

2 4 6
Re[ y ]

Fig. G-1. Location of zeros for Ai' [ y ] - iBi' [ y ] and also for Wl .

Ai[ y] ± i Bi[ y] = 2( p ± )1/ 2 Ai yp m

[ ] 
Ai ′[ y] ± i Bi ′[ y] = 2( p m )1/ 2 Ai ′[ yp m ] 

Ai[ y] + p + Ai[ yp + ] + p − Ai[ yp − ] = 0  (G-25)

Bi[ y] + p + Bi[ yp + ] + p − Bi[ yp − ] = 0 

p ± = exp[ ±2πi / 3] 

and
1/ 4
[ [ ] ]
W Ai yp + , Ai[ yp − ] = i(2π ) −1, W Ai[ y], Ai yp ± = ±( p m ) [ [ ]] (2π ) −1 
 (G-26)
W [ Ai[ y], Bi[ y]] = π −1 

where W [ X , Y ] is the Wronskian of X and Y . Thus, Eq. (G-24) can be

rewritten in the form

Wl± = 4πn1/ 2 (n Ai[ yp ± ]Ai ′[ yp ˆ ± ])

ˆ ± ] − n −1 Ai ′[ yp ± ]Ai[ yp 
 (G-27)
Wl = 4πn1/ 2 (np − Ai[ yp + ]Ai ′[ yp ˆ − ]) 
ˆ − ] − n −1 p + Ai ′[ yp + ]Ai[ yp

When N is small, even though Nxo >> K x o , we can factor out the n −1 term in
Eq. (G-27) and expand the remaining n 2 term in powers of N. The problem
term is the first involving the zeroth power of N. We can write this term in the
form
596 Appendix G

ˆ − ]
d  Ai[ yp
p − Ai[ yp + ]Ai ′[ yp
ˆ − ] − p + Ai ′[ yp + ]Ai[ yp
ˆ − ] = Ai[ yp + ]2   (G-28)
dy  Ai[ yp + ]

Also, in this case when Im[ yp ˆ − ] ≠ 0 , one can use the exponential asymptotic
forms for Ai[ yp ˆ − ] and Bi[ yp
ˆ − ] given in Eqs. (3.8-7) and (3.8-8).
For decreasing N , the zeros of Wl for a fixed value of κ drift
logarithmically downward and to the right in Fig. G-1. For very small N,
Nxo >> K x o , the locations of the zeros are given by the asymptotic condition

(
πNK x2o y( Ai[ y]2 + Bi[ y]2 ) − Ai ′[ y]2 − Bi ′[ y]2 = i ) (G-29)

Figure G-1 shows that, for a fixed value of N, the zeros of Wl migrate
upward from the real axis for increasing integer values of the index κ . As was
the case for the perfectly reflecting sphere, here it also can be readily shown
that U ( y) → 0 exponentially for increasing Im[ yκ ] when (θ − θ o ) ≤ 0 . In this
case the individual contribution of the zeros to the scattering coefficients
diminishes exponentially as their distance above the real axis increases.

G.4 Aggregate Scattering

For the aggregate scattering coefficient, which is given by Eq. (3.5-15b),
we see that its poles are located where Wl + Wl+ = 0 . From Eqs. (G-24) and
(G-26), it follows that

2πn1/ 2 (n Ai[ yˆ ](Ai ′[ y] − i Bi ′[ y]) − n −1 Ai ′[ yˆ ](Ai[ y] − i Bi[ y]))


or

Wl + Wl+ =˙ 4πn1/ 2 (( p + )1/ 2 n Ai[ yˆ ]Ai ′[ yp + ] − ( p − )1/ 2 n −1 Ai ′[ yˆ ]Ai[ yp + ]) =

  d  Ai[ yp + ] ( 
 4π ( np + −1/ 2
) Ai[ yˆ ] Ai[ ˆ
y ]   + 2 N + N 2 ) p + Ai ′[ yp + ]
  dy  Ai[ yˆ ]  
(G-30)

Unfortunately, the zeros of Wl + Wl+ are located in an infinite string along and
slightly above the negative real y-axis. Figure G-2 shows the first several zero
points for two fixed values of N . Here, contrary to the case shown in Fig. G-1,
Im[ yκ ] either decreases or does not increase sharply with increasing κ . Also,
the numerator term U ( y) for Sl does not converge to zero for large values of
κ . Therefore, numerical convergence problems result in attempting to evaluate
the aggregate scattering field using the poles of Sl in complex spectral number
Using the Sommerfeld–Watson Transformation 597

1.0

1
2 k= 0.8
3 k=
k=
Nxo = 83 p 0.6

Im [ y ]
+
W (yk ) + W ( yk ) = 0
l l
0.4
Nxo = 830p

3 0.2
2
k= k= 1
k=
0.0
-10 -8 -6 -4 -2 0
Re[ y ]

Fig. G-2. Zeros of Wl + Wl + in the complex plane, from

Eq. (G-30)

space. One can sum up the contributions from the individual j th -degree
scattered waves up to some specified truncation degree using the zeros of Wl
and the residues from the j th -order pole approach outlined earlier. This would
require successive numerical differentiation of the appropriate expressions for
U ( y) and V ( y) per the requirements given in Eq. (G-14) for a jth-order pole.
A numerical integration approach aided by the stationary phase technique,
which has been used extensively in this monograph, seems simpler for smaller
magnitudes of N, roughly N < ~ 2 K x−o2 . But it requires high precision and also
dense sampling to avoid aliasing, and when multiple rays are involved, it also
requires care in the vector addition of the field contributions from these
separate stationary phase neighborhoods in spectral number. For larger
magnitudes of N, the complex spectral number summation technique is
preferable in the region where (θ − θ o ) ≤ 0 , even with the aforementioned
numerical precision problems and the convergence problem for the aggregate
scattering coefficient. Also, near shadow boundaries, numerical integration of
the scattering integrals when N is large becomes a struggle because of the high
accelerations in phase in the scattering coefficients and the large number of
stationary points in spectral number.
The bottom line is that there seems to be no free lunch in wave theory for
calculating accurately in all regimes the complete field vector of an
electromagnetic wave that has passed through a transparent refracting sphere.
Parabolic equation methods such as the multiple phase screen approach to
propagate the wave through a refracting sphere is perhaps the most promising,
but multiple internal reflections from the sphere ( j = 2, 3L ) would require
598 Appendix G

careful treatment. For limb sounding from a LEO, a thin-screen/scalar

diffraction model offers an attractive alternative, provided in-screen caustics are
avoided.

Reference
[1] M. Abramowitz and I. Stegun, eds., Handbook of Mathematical
Functions With Formulas, Graphs, and Mathematical Tables, National
Bureau of Standards Series 55, Washington, DC, 1964.