izunfl
inf...

 

 

3. an ,
xx. . . 3 5 a ..
Quadrsxufl . ;§QWMQ
. _ . . mm

4 v > 1“ 4-

.3.“

5‘! uhpliz
.u 1.2.5,
24%;.»2.
a 1...: 3 .
:Lflfiufiufla

Anna». .1... .

.
’..J.v I!
.3.) my...“

 

 

95v
3... .
34:..qu
l

a.”

L:

. .I

54“}

\Vkrltarx»

I: 1.4 2
8...!!- v.
:t

. 9
at...» a
.91....

51.3.3.4.
1\yvflar.uu..v.rfirfl

a... . fldufibuufiu
ammrflgflcazrfifl
vi 1))... c... I.

3 $3931.. V
03. uteruuuwufimgfim
.41 wriniwnm .fi):

{mm ; [.13

9.... . i. .. ‘ .
41...»: 5 w... . 3.9.91
4113...?! ~ .5. (all!
.1; :

 

 

 

a. 93.... 34% .5?

 

 

T4815

3.00}

This is to certify that the
thesis entitled

TRANSIENT REFLECTION OF PLANE WAVES FROM A
LORENTZ-MEDIUM HALF-SPACE

presented by
STEVEN MICHAEL COSSMANN

has been accepted towards fulfillment
of the requirements for the

MS. degree in Electrical Equineerigq

(ZS/2%” /-——/

 

 

Major Professor’s Signature

0 5’ WL 2006

 

Date

MSU is an Affirmative Action/Equal Opportunity Institution

 

LIBRARY
Michigan State
University

 

 

-.-.---.---.-.-.-.-,_.- -.- -

PLACE IN RETURN BOX to remove this checkout from your record.
To AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE

DATE DUE

DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2/05 p:/ClRC/DaleDue.indd-p.1

 

TRANSIENT REFLECTION OF PLANE WAVES FROM A
LORENTZ-MEDIUM HALF-SPACE

By

Steven Michael Cossmann

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE
Electrical and Computer Engineering

2006

ABSTRACT

TRANSIENT REFLECTION OF PLANE WAVES FROM A
LORENTZ-MEDIUM HALF-SPACE
By

Steven Michael Cossmann

For materials that exhibit resonances, or have parameters that vary rapidly with
frequency, the Lorentz model is becoming increasingly more popular, especially when
dealing with signals with content at optical frequencies. The propagation and reflec-
tion of transient pulses by these media is of particular interest.

It has been shown previously that the transient plane-wave field reflected from a
Lorentz-medium half-space can be represented as an infinite sum of fractional-order
Bessel functions. In order to gain more physical insight into the physical behavior
of the problem, in this thesis a closed form solution is formulated which has no
infinite sums. The solution is obtained by taking an inverse Laplace transform of
the frequency-domain reflection coefficient. This is accomplished by rearranging the
frequency-domain reflection coefficient into separate terms which have inverse Laplace
transforms which are found in standard tables.

The results obtained using the closed-form solutions are verified through com-
parison with an inverse fast Fourier transform of the frequency-domain reflection

coefficient.

For My Parents

iii

ACKNOWLEDGMENTS

There are many people who I would like to thank for all their help and support
throughout my six years at Michigan State. First, I would like to thank my par-
ents, Greg and Marianne. You paid my tuition throughout my undergrad career and
supported me every step of the way.

I would like to thank Leo Kempel, for making sure I had work and a paycheck
and helping to get me interested in electromagnetics, which may not have happened
without your guidance. Thanks for convincing me to stay on for my masters degree,
it has been worth the time and effort.

I would also definitely like to thank Ed Rothwell for supporting me through my
masters degree and helping guide me throughout. The entire basis for this thesis
was your idea and our weekly meetings have helped me to no end. It was a great
advantage for me to have someone willing to give me advice whenever I needed it and
who was always willing to answer questions.

And thanks Shanker B, since almost all of my graduate EM classes were taught
by you and I feel like I’ve definitely learned something during them.

I have to thank all the guys in the office, Andrew, Brad, Dan, Gary, Greg and
Jeong. It has been great working at a place with people who are easy to get along
with and make work actually fun. Our weekly trips to Oodles (for those that went)
kept me happy, fed and fat for 2 years.

And I would like to thank my girlfriend Alexis, for sticking with me through these

years when I have had no money to take you anywhere nice.

iv

TABLE OF CONTENTS

LIST OF FIGURES ................................ vii
CHAPTER 1

Introduction and Background ........................... 1
CHAPTER 2

Frequency Domain Reflection Coefficients .................... 3

2.1 The Frequency Domain Wave Equation ................. 3

2.2 Reflection from a Half-Space ....................... 5

2.2.1 TE Polarization .......................... 6

2.2.2 TM Polarization ......................... 8
CHAPTER 3

TE Plane Wave Reflection for the Optical Case ................. 14

3.1 Laplace Domain Representation ..................... 15

3.2 Time-Domain Reflection Coefficient ................... 17

3.2.1 Inversion of the 91(3) and {12(8) Terms ............. 18

3.2.2 Final Expression for P(t) ..................... 25

3.3 Numerical Results ............................. 26
CHAPTER 4

TM Plane Wave Reflection for the Optical Case ................. 31

4.1 Laplace-Domain Representation ..................... 31

4.2 Time-Domain Reflection Coefficient ................... 32

4.2.1 Inversion of the 971(3) Terms ................... 35

4.2.2 Final expression for G (t) ..................... 38

4.2.3 Inversion of the C (3) Term .................... 39

4.3 Numerical Results ............................. 41
CHAPTER 5

TE Plane Wave Reflection for the General Case ................. 51

5.1 Laplace Domain Representation ..................... 52

5.2 Time-Domain Reflection Coefficient ................... 54

5.2.1 Inversion of the C(s) Term .................... 56

5.2.2 Inversion of the C (9) Term .................... 58

5.3 Numerical Results ............................. 59

CHAPTER 6

TM Plane Wave Reflection for the General Case ................. 65
6.1 Laplace Domain Representation ..................... 65
6.2 Time-Domain Reflection Coefficient ................... 67

6.2.1 Inversion of the C(s) Term .................... 71
6.2.2 Inversion of the C (5) Term .................... 73
6.3 Numerical Results ............................. 75

CHAPTER 7

Conclusions ..................................... 85
7.1 Suggestions for Future Work ....................... 85

APPENDIX A

Fortran Code for Optical TE Case ........................ 87

APPENDIX B

Fortran Code for Optical TM Case ........................ 103

APPENDIX C

Fortran Code for General TE Case ........................ 122

APPENDIX D

Fortran Code for General TM Case ........................ 141

BIBLIOGRAPHY ................................. 164

vi

Figure 2.1

Figure 2.2

Figure 2.3

Figure 3.1

Figure 3.2

Figure 3.3

Figure 4.1

Figure 4.2

Figure 4.3

Figure 4.4

Figure 4.5
Figure 4.6
Figure 4.7

LIST OF FIGURES

Geometry for an incident wave at a discontinuity between two ma-
terial regions

ooooooooooooooooooooooooooooo

Fields for a TE—polarized plane wave incident on a material half-
space.

ooooooooooooooooooooooooooooooooo

Fields for a TM-polarized plane wave incident on a material half-
space.

ooooooooooooooooooooooooooooooooo

Time-domain reflection coefficient with incidence angle 6 = 30°
and material parameters em = 4.0 x 101634, b2 = 20.0 x 10323'2,

6 z 0.28 x 10163—1. This choice of parameters corresponds to case
1, Eq. (3.56). .............................

Time-domain reflection coefficient with incidence angle 6 = 30°
and material parameters em = 2.0 x 10153‘1, b2 = 20.0 x 10293’2,
(5 = 0.28 x 10163-1. This choice of parameters corresponds to case
2, Eq. (3.53). .............................

Time-domain reflection coefficient with at an angle of 6 = 30° with
parameter choices of em = 2.0 x 10153—1, b2 = 20.0 x 10323—2,
6 = 0.28 x 10163’1. This choice of parameters corresponds to case
3, Eq. (3.58). .............................

Time-domain reflection coefficient with incidence angle 6 = 30°
and material parameters em = 4.0 x 10163—1, b2 = 20.0 x 10323-2,
6 = 0.28 x 10163'1. This choice of parameters corresponds to case
1, Eq. (4.33). .............................

Time-domain reflection coefficient with incidence angle 6 = 30°
and material parameters em = 2.0 x 10153’1, b2 = 20.0 x 10293-2,

6 = 0.28 x 10163-1. This choice of parameters corresponds to case
2, Eq. (4.31). .............................

Time-domain reflection coefficient with at an angle of 6 = 30° with
parameter choices of 0’0 = 2,0 x 10153—1, b2 = 20.0 x 10323-2,

6 = 0.28 x 10163—1. This choice of parameters corresponds to case
3, Eq. (4.34). .............................

Time-domain reflection coefficient with incidence angle 6 = 50°
and material parameters em = 4.0 x 10163-1, b2 = 20.0 x 10323’2,
6 = 0.28 x 10163—1. This choice of parameters corresponds to case
1, Eq. (4.33) with non-causal term from Eq. (4.41) .........

Non-causal C(t) —- 6(t) for case 1 results ...............
Non-causal [C(t) — 6(t)] * C(t) for case 1 results.
Close—up of the non-causal [C(t) — 6 (15)] * C(t) for case 1 results. .

vii

11

12

13

28

29

30

44

45

46

47
48
49
50

Figure 5.1

Figure 5.2

Figure 5.3

Figure 6.1

Figure 6.2

Figure 6.3

Figure 6.4

Figure 6.5
Figure 6.6

Figure 6.7

Time-domain reflection coefficient with incidence angle 6 = 30°
and material parameters em = 4.0 x 10113-1, b = 6.24 x 10113—1,
(5 = 2.5 x 10108—1, 600 = 2, M = 1. This choice of parameters cor-
responds to case 1, Eq. (5.34) .....................

Time-domain reflection coefficient with incidence angle 6 := 30°
and material parameters wo = 4.0 x 10103-1, b = 6.24 x 10108—1,
6 = 2.5 x 10113-1, 600 = 2, pr = 1. This choice of parameters cor-
responds to case 2, Eq. (5.32) .....................

Time-domain reflection coefficient with at an angle of 6 = 30°
with parameter choices of wo = 4.0 x 1093"1, b = 6.24 x 10113'1,
6 = 2.5 x 10103-1, 600 = 2, p1. = 1. This choice of parameters cor-
responds to case 3, Eq. (5.35) .....................

Time-domain reflection coefficient with incidence angle 6 = 30°
and material parameters wo = 4.0 x 10113-1, b = 6.24 x 101134,
6 = 2.5 x 10103-1, coo = 2, Mr 2 1. This choice of parameters cor-
responds to case 1, Eq. (6.27) .....................

Time-domain reflection coefficient with incidence angle 6 = 30°
and material parameters wo = 4.0 x 10108—1, b = 6.24 x 10103—1,
(5 = 2.5 x 10113‘1, 600 = 2, ,ur = 1. This choice of parameters cor-
responds to case 2, Eq. (6.26) .....................

Time-domain reflection coefficient with at an angle of 6 = 30°
with parameter choices of tag = 4.0 x 1098—1, b = 6.24 x 10113_1,
5 = 2.5 x 10103‘1, 600 = 2, or = 1. This choice of parameters cor-
responds to case 3, Eq. (6.28) .....................

Time-domain reflection coefficient with at an angle of 6 = 60°
with parameter choices of Lao = 4.0 x 1011.3‘1, b = 6.24 x 10113-1,
(5 = 2.5 x 10103—1, 600 = 2, pr = 1. This choice of parameters cor-
responds to case 1, Eq. (6.27) with a non-causal term from Eq.
(6.34). .................................

Non-causal C (t) term for case 1 results. ..............

Non-causal C (t) term convolved with G (t) minus the delta function
term for case 1 results .........................

Non-causal C(t) term convolved with the delta function term for
case 1 results ..............................

viii

62

63

64

78

79

80

83

84

CHAPTER 1

INTRODUCTION AND BACKGROUND

Simple electromagnetic materials are oftenmodeled as having constitutive parame-
ters that are constant with frequency. More complicated materials with frequency-
dependent parameters are often represented using the Debye model, which is valid
for many liquids at moderate frequency. For materials that exhibit resonances, or
have parameters that vary rapidly with frequency, the Lorentz model is becoming in-
creasingly more popular, especially when dealing with signals with content at optical
frequencies [1].

The propagation and reflection of short duration transient pulses in these disper-
sive materials has generated much interest, and the effects of dispersion on the shape
of the propagating wave, particularly in the Sommerfeld and Brillouin precursors, has
been extensively studied [2H5]. Recently, the use of short pulses to probe materials
has prompted the investigation of the reflection of transient waves from material half
spaces of various types [6]—[9].

The transient field reflection of a plane wave from a Lorentz-medium half-space
has been calculated previously for both TE polarization [10] and TM polarization
[11] by expanding the frequency-domain reflection coefficient using a simple expan-
sion method and taking the inverse Laplace transform of each term in the infinite
series. This produces a result which is given as an infinite sum of Bessel functions
of fractional-order. While this method produces a convenient solution for numerical
techniques, it provides little insight into the behavior of the reflected field.

This thesis proposes a method for formulating the transient reflection coefficient
by rearranging the frequency-domain reflection coefficient into a form which enables

the inverse Laplace transform to be taken without making an expansion as done

previously. This leads to a solution containing a finite number of convolutions of
Bessel functions and exponentials, which gives a clearer picture of the behavior of
the reflection coefficient. Material parameter and incidence angle choices determine
whether the convolution terms contain ordinary or modified Bessel functions. In
order to obtain a complete solution to the problem, different sets of parameters are
chosen which are representative of all the possible combinations of standard Bessel
functions and modified Bessel functions. The various combinations of standard Bessel
functions and modified Bessel functions allow for the general behavior of the reflected
wave to be predicted from the material choices. A more oscillatory reflection would
be expected from a solution containing only standard Bessel functions, while a less
oscillatory reflection would be expected from a solution containing only modified
Bessel functions.

With certain material parameters, there are some angles of incidence in which the
final transient expression for a TM-polarized wave contains a non-causal term. For a
TE-polarized wave this is not the case as long as the material is not highly magnetic.
Explanation for these non—causal terms are examined in more detail as they arise in
formulation.

Once the result is formulated, the transient reflection coefficient is computed nu-

merically and then compared to the inverse fast Fourier transform (FF T )

CHAPTER 2

FREQUENCY DOMAIN REFLECTION COEFFICIENTS

In order to derive the time-domain reflection coefficient for a Lorentz-medium half-
space, the frequency domain reflection coeflicient must first be formulated. The most
general way to accomplish this is to enforce tangential field continuity across the
boundary of the half-space. In this chapter the frequency domain reflection coefficient
of an obliquely incident plane wave from a material half-space is derived. This solution

is valid for any homogeneous, isotropic material.

2.1 The Hequency Domain Wave Equation

Maxwell’s equations for a source-free region of linear, isotropic, homogeneous material
for a time-harmonic field are given in point form in terms of E(r,w) and H(r,w) as

[12]

V x E = -—jw,uH, (2.1)
V x H = jch, (2.2)
V - E = 0, (2.3)
V - H = 0, (2.4)

where c is the frequency dependant complex permittivity. Taking the curl of (2.1)

produces

VXVXEz—jwp.(VxH)

= w2ucE, (2.5)

where (2.2) has be substituted for V x H. Utilizing both (2.3) and the vector identity

V x V x E = V (V - E) — V2E, allows (2.5) to be rewritten as
v2}; + 1813 = 0. (2.6)

where the wave number in the medium is given as k = , Me. This is the homogeneous
vector Helmholtz equation for the electric field. By a similar process, the homogeneous

vector Helmholtz equation for the magnetic field can be written as
V2H + k2H = o. (2.7)

In cartesian coordinates, each component of E and H must satisfy the scalar
Helmholtz equation

V211; + W) = 0. (2.8)

The expression for the electric field of a plane wave can be written as
E = EOe—ik", (2.9)
where k is the wave vector defined as
k = xkx + yky + 21.3,, (2.10)
with |k| = k and r is the position vector defined as

r=x.r+yy+zz. (2.11)

The magnetic field for a uniform plane wave is related to the electric field by

k E
H = x . (2.12)
w

 

2.2 Reflection from a Half-Space

The problem of interest is a plane wave in free space (region 1) incident on a material
half-space. A portion of the wave is reflected and a portion is transmitted into the
material half-space (region 2). The problem geometry can be seen in Figure 2.1.
Examining the geometry, it can be seen that since the wave is propagating in the x-z
plane, the field is invariant in the y-direction which means kg = 0. It can also be seen

from examining the geometry that. the wavenumber for the incident wave is given as

kw = k0 sin 6,, (2.13a)

kzo = 1:0 cos 6,, (2.13b)

and the wavenumber for the transmitted wave is given as

kg, = k sin 6t, (2.14a)

k. = kcos6t. (2.14b)

Any uniform plane wave incident on a planar surface may be decomposed into two
orthogonal polarization components, one perpendicular and one parallel to the plane
of incidence. The perpendicular polarization is also known as TE polarization, and the
parallel polarization is also known as TM polarization. Each of these polarizations
can be solved for separately and their solutions can be added together to form a

complete solution for the total field.

2.2. 1 TE Polarization

The fields for a TE—polarized plane wave incident on a material half-space can be seen

in Figure 2.2. The fields for the incident plane wave in free space can be given as

E: = yEae—jkma: sin 6i+z c056,), (2,15)
. Ei . ‘.
1 = 7’30 (—x cos 6,- + 23in 6,) €_JA0(I ””1 01'” COS 6’)- (2-16)

were 710 = \/[1.0/€0 and k0 = w,/,u()c(). The total field in region 1 can be expressed as

the incident plane wave plus a reflected plane wave. The reflected field is given as

E1 : yEse—jk0(.r sin 61-—:: cos6r), (2.17)
E" . .

H: = —Q (x cos 6,. + 25m.) e—J‘0(1‘S”‘0T—ZCOS°’I). (2.18)
710

In region 2, the area to the right of the interface, the transmitted field can be repre-

sented as a plane wave given as

 

Et = yEée—jk(.rslngt+zCOSOt), (219)
Et ., ,.

Ht ___ 0 (—)°(COS 6t + 28111 6t) e-jk(.1‘bm 9t+ZCOSOt), (220)
77

were 17 = (/ a/c. In order to have phase continuity across the interface for all values

of 1:, the exponential terms must be equal at the interface, 2 = 0. This implies that

kg sin 6., = k0 sin 6r = ksin 6t. (2.21)

It can be easily deduced from this equation that

 

sin 9t kn /60/-10
sin 6,- k 611 ( 22)

The boundary conditions at the interface requires the tangential electric fields to be

continuous across the interface. This requires

mm

 

zx(1+E1)

 

2:0

Using (2.23) and the fact that the phase is continuous across the interface, and using

the expressions for the electric fields it can be seen that
%+%=%. (we

Dividing this by the incident field amplitude leads to

H Et

0 _ O

1+E‘E’

0 0
I+F_L =71, (2.25)

where Pi is the reflection coefficient and T .L is the transmission coefficient.
The boundary conditions also require that the magnetic field be continuous across
the interface. This requires

=2xHj (2%)

 

 

2x( 1+Hfl) .
z=0

~
Ar

Using (2.26) and the fact that the phase is continuous across the interface, and using

the expressions for the magnetic field it can be seen that

Er Et
cos 6,- + —° cos 6,» = ————Q cos 6:. (2.27)
710 77

i
.&

770

Dividing by the incident wave amplitude again yields

 

1 E’" Et
——cos6,~+ 2.0 cos6r=——iicos6t,
720 one Eon
1 F T
—— cos6, + i cos6r = ——i cos6t.
720 710

Combining (2.25) and (2.28) allows both F_L and TI to be solved for as

 

T _ 2Zi
i _ Z i + 20’

1“ = M
i zL + 20’

with the definitions
770
Z =
0 cos 6.,"
7) kn

 

 

 

Z 1 = .
COS 6‘ \/k2 — kg sin2 6,-

In the above expression, using (2.22) cos 6t is given as

 

k2 -— kg sin2 6,-
k

 

COS 6t =

2.2.2 TM Polarization

(2.28)

(2.29)

(2.30)

(2.31)

(2.32)

(2.33)

For a TM-polarized plane wave, the method to formulate the reflection and transmis-

sion coefficients is similar to the method for the TE case. The fields for a T M-polarized

plane wave incident on a material half-space can be seen in Figure 2.3. The fields for

the incident plane wave in free-space can be given as

f, = E5 (2 cos 6,- — 23m 0,) e‘fkonsingz'Hwa),

(2.34)

”_ _yE_0€-jk0(:r sin62- +z c056 ,)

The reflected and transmitted fields are given as

EH — —E0 (xcos6r — zsin6

(i = _y§_(_)_e—jk0(.rsi116r—zcos6r),
Tl0

and

Elt|_ — E6 (xcos 6t — zsin 6t) e—jk($°i°°t+z°°°0t),

t
Hltl ___ yfle—jk(2:sin6t+zcos6t),
77

)e-jk0(;r $11107~—Z cos 67-),

(2.35)

(2.36)

(2.37)

(2.33)

(2.39)

respectively. J list. as in the TE case, enforcing phase continuity across the boundary

 

 

 

leads to
61 = 61‘,
sin 6t __ k0 _ €0,140
sin 6,- - k _ ecu.
Using
X( II+EII) ‘2XEIli ,
2:0 =0
yields

E6 cos 6,- + E6 cos 6,. = E8 cos 6t.

Dividing by the incident wave amplitude leads to

E6 t
0036, +— 0003 6r— — —9 cos 6t,
E0 E0

cos 6,- + F“ cos 6T = T“ cos 6t.

(2.40)

(2.41)

(2.42)

(2.43)

Enforcing tangential field continuity for the magnetic field across the interface requires

 

2:0

 

ix(i+HU

’

2:0

which leads to
no no 7) '

Dividing by the incident. wave amplitude leads to

1 E51_E51
770 E6770 E577,
i+fl=fl,
720 770 77

Combining (2.43) and (2.46) allows F“ and T“ to be solved for as

 

 

with the definitions

20 = 770 C08 61',

 

Z“ = ncos6t = %\/k2 — kgsin2 6,.

10

an)

(2.45)

aw)

(2.47)

(2.48)

(2.49)

(2.50)

Region 1
(#0 7 60)

Region 2

 

(Me)

U

Figure 2.1. Geometry for an incident wave at a discontinuity between two material

regions

11

kt

 

H

 

 

Figure 2.2. Fields for a TE—polarized plane wave incident on a material half-space.

12

z'

2'

ilfL'

'
"
O

’01‘9
‘9 |
‘ 0" .I t

 

N"

Figure 2.3. Fields for a TM-polarized plane wave incident on a material half-space.

13

CHAPTER 3

TE PLANE WAVE REFLECTION FOR THE OPTICAL CASE

The transient reflection from a Lorentz-medium half-space for a plane wave can be
found analytically using an inverse Laplace transform technique. The simplest case
to examine is that of a TE—polarized plane wave. A steady-state TE—polarized plane
wave of frequency w is obliquely incident on an interface separating free space (region
1) from a homogeneous Lorentz medium (region 2). The angle of incidence 6 is
measured from the normal to the interface. This geometry can be seen in Figure 2.2.
As shown in Chapter 2, the reflection coefficient for a plane wave from a material

half-space can be expressed as

(3.1)

For a TE-polarized plane wave, the impedance of the incident wave is ZO = 710/ cos 6

and the wave impedance of the transmitted wave is given by

(3.2)

where

720 = Vito/60a (3.3a)

 

77 = 9/6, (3-3b)
k; = \/k2 -— kg sin2 6, (3.3e)
k0 = Law, (3.3d)

k = w\/p—.c, (3.3e)

c = 606,.(w). (3.3f)

14

The relative permittivity of a single-resonance Lorentz medium is given by [1]

2
w (e —c )
(7(a)) = 630 — 2 0— 2:w6 fwg, (3.4)

 

where 220 is the resonance frermency and 6 is the damping coefficient. 63 is the static
permittivity which is the value of 6,- when w = 0, and 600 is the optical permittivity
which is the value of 6,- as w —> 00. In optical problems, the case when 600 = 1 and
the relative permeability n,- = 1 is of the most interest. In this case, (3.4) can be

rewritten as [4]

 

()2

3 3.5

where b is the plasma frequency of the medium.

3.1 Laplace Domain Representation

A frequency domain quantity can be generalized to the Laplace domain using 8 = jw.
Using this, (3.5) can be represented in the Laplace domain as
b2

6. 3 21+ . 3.6
r() 52+263+wg ( )

 

The wave number in the Lorentz medium can then be calculated as

13(3) = —js\/p_.c
—js,/p.ooc \/1+

 

 

 

 

32 + 263 + (4)8
+ 268 + “’0 + b2
= — 's,/ . c 3.7
J WNW/82 302+265+w ( )
It is then convenient to make the substitution
32 + 263 + tug = (s — 31)(s — 32), (3.8)

15

where
512 = —6 :l: [/62 -— (428, (3.9)

and also the substitution
2 , ,2 2 _
s +263+w0+b — (S—S3)(S—S4), (3.10)

where
33,4 = —6i (M? —w3—b2. (3.11)

The Laplace domain reflection coefficient can then be written as
__ k30 " kz
_ A‘20 + 19::
cos6 — \/(k/k0)2 — sin2 6
cos6 + \/(k/k0)2 — sin2 6
cos6 — (/[(s — s3><s — saws — 31)(s — 82)] — sin2 6
C059 + \/I(8 - 53)(s - SUI/[(5 - 81)(8 ~ 82)] - Sin2 9
cos6\/s — sh/s — 32 —- \/(s — 33)(s — 34) — sin2 6(3 — sl)(s — 32)

= . (3.12)
cos 6\/s — sh/s — 82 + (/(s —- S3)(S —— S4) — sin2 6(3 — sl)(s — 32)

F(s)

 

 

 

 

 

 

 

 

 

 

 

Expanding (s — sl)(s — 32) and (s -— 33)(s -— 54) using (3.8) and (3.10) respectively

allows the reflection coefficient to be written as

 

[(9) = cos6\/s — 51\/s — 32 — [32(1— sin2 6) + 265(1— sin2 6)+

 

(“’8 + b2 — tug sin2 6)]1/2] / [cos6\[s — 81\/8 — 82 +[s2(1— sin2 6)+

268(1 — Sin2 6) + (W8 + b2 _ “’8 $1112 9)]1/2]

16

 

 

\/s—sl\/s—52-— fi+268+wg+b2/cosz6

 

 

 

 

— \/8—81\/S-82+\/;2+258+w8+b2/COS29. (3.13)
It is then convenient to make another substitution with the definition
32 + 268 + 918 + b2/ cos2 6 = (s — 35)(s —— 86), (3.14)
where
855 = —6 i \/62 — (.08 — b2/cos2 6. (3.15)

Substituting (3.14) into (3.13) allows the reflection coefficient in the Laplace domain

to be rewritten as

 

 

_ ¢s_.s,,/S_.2_¢s_35,/s_s,
\/s—31\/s—32+\/s—35\/s—56'

 

F(s) (3.16)

 

 

3.2 Time—Domain Reflection Coefficient

lsing the frequency—domain reflection coefficient found in Section 3.1, the time-
domain reflection coefficient can be found using an inverse Laplace transform [13].
In order to perform the inverse Laplace transform, the frequency-domain reflection
coefficient may be rearranged into a more manageable form. Rationalizing the de-

nominator in (3.16) leads to

 

 

_(Js—81¢§—82 — «o—ost/s—M _ fl
”3) ‘ (8-81)(8—S2) — (s —85)(8—86) “ D‘ ‘3'”)

Examining the denominator term in more detail it can be seen that

D = (8 — 51)(8 - 82) — (8 — 85)(3 7* 86)
2

b
= (32+263+w8) _(S2+268+w3+C—0:‘2_6)

17

= —b2/ cos2 6

= -B2. (3.18)

The numerator can then be expanded to

 

 

N = (s— 31)(s—32) —2\/s —S1\/8—82\/S— 35\/5_ 56+(5‘35)(3—36)1
= [(.s—31)(3—32) — \/8—81\/8—32\/3‘35\/3—36l +
[(8-85)(8—86)-\/S-81\/S-S2\/S-85\/8-88],

= 91(8) + 92(9)- (3-19)

 

 

This last separation is taken in order to make the inverse Laplace transform easier to
perform.

The invertible form of the reflection coefficient is then given by
4321‘s) = 91(5) + 92(3), (320)

where each term, 91(3) and 92(8), can then be inverted separately.
3.2.1 Inversion of the 91(3) and 92(3) Terms

W hen inverting the term 91(8), it is convenient to make the substitutions

 

812 = —d :1: A1, (3.21)
S536 = —d :i: A5, (3.22)
where
A1 = 62 — Lug, (3.23a)
A5 = (/62 -— tag - 32, (3.2310)

and B2 is as defined in (3.18). In order to perform the inverse Laplace transform,

g1(s) may be rearranged into

91(8) = (8 - 51)(S - 82)(8 - 85)(8 - 86)><
1 l 1
(8 - 85)(8— 86) _ \/3 —81\/8——82\/S - Ssx/S— 86 '

 

 

(3.24)

 

Using partial fraction expansion, the first term in the brackets can be rewritten in

the form

1 K K
= 1 + 2 , (3.25)
(8—85)(5-86) 8-85 8-86

 

 

where K1 and K 2 are given by

1 1
K = = — .2 -
l 1
K = = ___. 3.26b
2 86 — 35 2A5 ( )

 

 

Since tag and 82 are both positive numbers, Re{s5,6} < 0 always. This means the

standard inverse Laplace transform [14]

1
3+5

 

+—> e_(3tu(t), (3.27)

can be used, where u(t.) is a unit step function which is zero for t < O and has unit

amplitude for t 2 0. Using this identity, (3.25) can be transformed into

 

1 1 s t t]
__ 5 _ 36‘ t
(8-35)(8-86) “‘3 2A5 [6 e u”
___ 3:5: [ekst _ e-Ast] u(t)
2A5 '
{at
= Tsinh()\5t)u(t). (328)
5

19

The second term in the bracket in (3.24) can be transformed into a convolution

between two terms if the following inverse Laplace transform is used [14]:

 

1 -1 1
mm ‘3 10 [2“ a”) “W (3‘29)

where 122(1) is the modified Bessel function of the first kind of order 71. Using this,

the second term can be transformed into

1
Fs—slvs-sws—sm—ss

 

 

t—* 5‘” [{10()\1t)11(t)} * {10(/\5t)U(t)}l- (330)

Using (3.8), (3.14), (3.28) and (3.30), and using the differentiation theorem, it is

possible to write gl(t) as

d2 d 2 d2 d 2 2
2‘: 26—
g1() (Cl—{2+ dt+w0)<d——t2+26dt+w0+8 x

6“” [3%— sinh(/\5 t)1 1(t—) {10(A1t) u (t)} * {IO(/\5t)u(t)}] . (3.31)

5

The next step is to carry out the derivatives. It is important to note at this point
that when taking a derivative in time of two signals convolved with respect to time,
the derivative only needs to be taken over one of the signals in the convolution. This

can be expressed mathematically as

d df dg
"CE” *9): dt *9— - f* I? (3.32)
Defining
E105) = 6"” isinhfls 011 (t )- {10(/\1l5)vl(t )} * {Io(/\5t)U(t)} . (3-33)

20

makes the solution easier to compute. The first derivative of §1(t) computed using

the product rule is given by

d

32W) = —6§1(t) + m). (3.34)

where §’1(t) is the exponential term multiplied by the derivative of the term in the

brackets in (3.31). Using the bessel function identity [15]

16(2) = 11a), (3.35)
and the identity

d

— t: t .

dtU() 6( )1 (336)

allows g’1(t) to be written as

EN) = {at [00811053140 — A5{10(/\11t)11(1f)} * {119501100} - 10(A1t)11(t)] - (3-37)

In the convolution term, the derivative was taken on the second modified Bessel
function containing the argument A515. The modified Bessel function of order zero
came from the fact that {10(A1t)u(t)} * {10(A5t)6(t)} = I0()\1t)u(t). The second

derivative is given by

Ema) = —6 (Edi-951(0) + 3:261“)

= -5 (—6§1 (t) + 3(a)) — 6330) + a’flt)

= 62M) — 26-91(0’ + 'g-a’m, (3.38)

21

where y’l’ (t) is defined as the exponential term multiplied by the derivative of the term

in the brackets in (3.37). Using the Bessel function identity [15]

14.0) = 112-1(1) + 12.1w, (3.39)

NIH

and again takin the derivative on the second modified Bessel function in the convo-

lution term allows 57(1‘) to be solved for as

2
iii“) = 8—“ [)5 Sth‘stNU) — % {10(A1f)u(t)} * {10(A5t)11(t)} -

2
% {10(/\1t)11(t)} * {12(A5t)U(t)} —— A111(/\1t)u(t) . (3.40)

Combining (3.33), (3.34) and (3.38), it can be seen that

d2 d _ _.
(2172 + ”a + w?) + B?) m) = 5221(1)- 26930 + atm—
262m> + 26130) + (wé + 82mm

= g’,’(t) — (62 — .33 — B2)§1(t). (3.41)
Using the definition for /\5 given in (3.23), (3.41) can be rewritten as
(l2 d 2 2 _ _[I 2_
afi+2éfi+w0+3 91(t)=gl(t)_)‘591(t)

2
= 5‘” [52‘5- {10()\1t)11(t)} * UoOsUUU) -12()\5t)u(t)} - A111(A1t)11(t)]

IlUSUUU)

= 8—6t [A5{10()\1t)u(t)} * { t

} — A111(/\1t)u(t)]

2 31(1), (3.42)

22

which takes advantage of the identity [15]

2n.
711201“)=In—1($)—Irz+1(I)- (3.43)
It is then possible to find
d — _ _,
5’11“) = —6h1(t) + h1(t), (3.44)

where 5,10) is the exponential term multiplied by the derivative of the term in the
brackets in (3.42). Performing the derivative on the first modified Bessel function in

. —I .
the convolution allows h1(t) to be written as

-I

hut) = 6“” [A1A5{11(A1t)u(t)} * {M

' 2
t } + ASS-[100x50 — 12(/\5t)]U(t)—

2
%[10(A1t) + 12(A1f)]11(t):l . (3'45)

The second derivative is then found to be

d2—

d _ d —/
37,711“) = ——6 [d—th1(t)] + a—ihdt)

= _5 [—551(t>+ Him] — 533(1) + Wt)

6251(1) — 265302) + E’l’u), (3.46)

where 5,1,“) is the exponential term multiplied by the derivative of the term in the
brackets in (3.45). When taking the derivative of the convolution term, the derivative

is taken on the first modified Bessel function. This leads to

—” —6t ”\fAS

mm = e —— (110(41):) + 120101140} *{

11050110)
. ——————}+

t

23

41’
713110113 + 13(31t)lu(t) +

 

*3
7111050 —13(/\5t)lu(t) -

Combining (3.42).(3.44) and (3.46) it can be seen that

d2 d _ _ __ _
(E? + 25a; + .33) M) = 62h1<t>— 2614(2) + hi’m—
25251u) + 2551(1) + .1331“)

= 751’“) _ (52 — 423) Wt). (348)

Using the definition for A1 given in (3.23), (3.48) can be rewritten as

d2 d 2 — —II 2-
a? + 26a +w0 111(13): (11(4) " Alhlm

which leads to

_ . . . . 9—)?
21(1):. 5‘ [—Aagxlalt).11(151)+1§12(.\,t)+Agx2(.\51)]+ 5, 16013.50)

 

where

 

In(.r) = I u(:c). (3.51)
The same approach can be taken to find

_ - . . - 22—12
22(1):. 4 [—Afxéwit)*Ilust)+4112<11t>+4312<45t)]— 5, 13013.52)

 

24

3.2.2 Final Expression for F(t)

Substituting gl(t) from (3.50) and 92(t) from (3.52) into (3.20) yields the final ex-

pression for the reflection coefficient

2 _ - . . .
I‘(t) = we 5t [A§/\§11(,\11)411(A5t) — A§I2(.\1t) — 2312mm] . (3.53)

The solution (3.53) is generally valid for all cases, but when A1 or A5 is imaginary,
the modified Bessel functions can be replaced with ordinary Bessel functions according
to three possible cases. Case 1 occurs when tag > 62. In this case both /\1 and A5 are

purely imaginary and defined as

 

42,: X1: tug — 62, —j/\5 = is = 313+ 13? — 62. (354)
Then, using the property [15]

I,,(j;r) =jan(17), (3.55)

allows the reflection coefficient to be rewritten as

2 _ —2—2~ — . — -3~ — —3~ —
W") = 13—28 M [A1A5J1()‘1t)* J1()‘5t) — 4132910 - A5J2(/\5t)] , (3-56)

where, similar to (3.51),

 

32(2) = u($). (3.57)

Case 2 occurs when (.08 + B2 < 62. In this case, both A1 and A5 are purely real and
the same expression from (3.53) can be used. Case 3 occurs when —B2 < tag —62 < 0.

In this case, A1 is purely real and A5 is purely imaginary. X5 is defined as it was in

25

(3.54) and the expression for the reflection coefficient can be written as

2 _ _2. . — . —3~ —
F(t) 2 —E6 6t [/\¥2\511(/\1t) * J1()\5t) + A¥12()\1t) + A5J2(A5t)] . (3.58)

3.3 Numerical Results

In order to validate the expressions derived in the previous section, the time—domain
reflection coefficient is evaluated numerically in Fortran, and then compared to the
inverse FF T of the frequency domain reflection coefficient given in (3.16). The For-
tran code is included in Appendix A. The inverse FF T was done using WaveCalc.
The frequency-domain data was zero-padded up to the maximum limit allowed by
WaveCalc, 32,768, before the inverse FFT was taken. Since the signal had already
decayed to zero by this point, no windowing was necessary. A set of parameters cor-
responding to each of the three possible cases is used. The number of points and step
size for the results varied between the different cases.

The first set of parameters is the same as those chosen by Brillouin [16]: wo =
4.0 x 1016 3’1, b2 = 20.0 x 1032 3'2, 6 = 0.28 x 1016 3‘1. This choice of parameters
corresponds to case 1. When computing the numerical results, for the frequency-
domain data 16,384 frequency points were calculated with a step size of 8,000 GHz.
In the time-domain, 4,096 points were calculated with a step size of 1x 10‘9 ns. Using
9 = 30°, (3.56) has been plotted in Figure 3.1 and compared to the inverse FFT. The
results show excellent agreement. Since this function includes only standard Bessel
functions, which are highly oscillatory, and no modified Bessel functions, which are
not oscillatory, the waveform is highly oscillatory and only lightly damped.

The next choice of parameters is: rug 2 2.0 x 1015 5’1, b2 = 20.0 x 1029 3‘2, 6 =
0.28 x 1016 3’1, which corresponds to case 2, (3.53). When computing the numerical
results, for the frequency-domain data 4,096 frequency points were calculated with a

step size of 400 GHz. In the time-domain, 4,096 points were calculated with a step size

26

of 1 x 10‘8 ms. The results are shown in Figure 3.2. Again, the closed-form expression
and the inverse FFT compare well. For this choice of parameters 62 > c128 + 82, and
the resulting waveform is overdamped, showing no oscillatory behavior and only a
single negative peak. Since the expression for case 2 only involves modified Bessel
functions, which do not have the oscillatory behavior of ordinary Bessel functions, this
observed behavior is easily predicted from the mathematical form of the expression.

The final choice of parameters, which corresponds to case 3, is: mo 2 2.0 x
1015 3’1, b2 = 20.0 x 1032 3‘2, 6 = 0.28 x 1016 3’1, which corresponds to case 3,
(3.58). When computing the numerical results, for the frequency-domain data 16,384
frequency points were calculated with a step size of 8,000 GHz. In the time-domain,
4,096 points were calculated with a step size of 1 x 10‘9 us. The analytic expression
again matches the inverse F FT, as seen in Figure 3.3. As expected, since 6 > tag, but
62 < c128 + B2, there is more damping and less oscillation than with case 1, but more
oscillation than with case 2. Here the expression for the reflection coefficient has a

combination of ordinary and modified Bessel functions.

27

x 1016 Case 1

 

 

 

 

 

 

 

0.5 -
E
.92
s 0-
3
0
C.
.9
3 -o.5
c
Q)
a:
_1 _
Closed Form Expression
-~X- -- Inverse FFT
_1.5 1 1 1 4 1
0 0.2 0.4 0.6 0.8 1

Time (fs)

Figure 3.1. Time-domain reflection coefficient with incidence angle 0 = 30° and mate-
rial parameters wo = 4.0 x 10168—1, b2 = 20.0 x 1032s”2, 6 = 0.28 x 10163-1. This
choice of parameters corresponds to case 1, Eq. (3.56).

28

x1013 Case 2

Reflection Coefficient

 

 

 

 

 

 

 

 

-9 ~ Closed Form Expression
Wx Inverse FFT
_10 1 1 1 1 1
0 1 2 3 4 5
Time (fs)

Figure 3.2. Time-domain reflection coefficient with incidence angle 0 = 30° and mate-
rial parameters wo = 2.0 x 10153-1, b2 = 20.0 x 10293—2, (5 = 0.28 x 10163’1. This
choice of parameters corresponds to case 2, Eq. (3.53).

29

x1015 Case 3

 

 

 

 

 

 

 

5 _
O .-
E
.92
s —5
a)
O
O
C
.9
8 ~10 —
:
OJ
Ct
_15 >-
Closed Form Expression
~ - 'X‘ ‘ - Inverse FFT
_20 l 1 l l J
O 0.2 0.4 0.6 0.8 1
Time (fs)

Figure 3.3. Time-domain reflection coefficient with at an angle of 6 = 30° with param-
eter choices of wo = 2.0 x 10153—1, b2 = 20.0 x 10325-2, 6 = 0.28 x 10163—1. This
choice of parameters corresponds to case 3, Eq. (3.58).

30

CHAPTER 4

TM PLANE WAVE REFLECTION FOR THE OPTICAL CASE

The reflection of a TM plane wave from a Lorentz medium half-space can be found
with a similar method as for the TE case. The problem is defined as in Chapter
3 with a plane wave from free space obliquely incident on a homogeneous Lorentz
medium at an angle 0 measured normal to the interface. This geometry can be seen
in Figure 2.3. As shown in Chapter 2, the reflection coefficient for a plane wave from
a material half-space can be expressed as

_ ZHW) — ZO

F||(LU) — m.

(4.1)

For a TM-polarized plane wave, the impedance of the incident wave is Z0 2 170 c056

and the wave impedance of the transmitted wave is given by

Z“(..;) = '” (4.2)

where the terms 132,77, and k are defined in (3.3). The relative permittivity is the

same as given in (3.5).

4.1 Laplace-Domain Representation
For a TM—polarized plane wave, the Laplace domain reflection coeflicient can be
written as

= k2 — 6T'kzO
19: + 6rig/2'0

\/(k/k0)2 — sin2 9 — er cosB
\/(k/k0)2 — sin2 6 + 6r cos 0

I‘(s)

 

 

 

31

e,- — sin2 6 — er cos6

\/ er — sin2 6 + 6,- cos 6

 

(4.3)

Combining the relative permittivity 61- into one fraction, it can be substituted into

(4.3) to give

 

[(9) = [\/(52 + 263 + 068 + b2)/(s2 + 268 + (.63) — sin2 6 — [(s2 + 263 + 1123 + b2)/

 

(52 + 265 + w8)] cos 6] / [\/(52 + 26.9 + 028 + b2)/(s2 + 26.9 + (12(2)) — sin2 6+
[(32 + 265 + 112(2) + ()2)/(32 + 263 + 1.03)] cos 6]

12
= [fi—les—sg [32+263+w3+b2— (82+26s+w3)sin26] / —

 

 

(s — 83)(8 — 34)cos6]/[\/s — ran/s — .52 [82 + 263 + 411(2) + b2—

]1/2

(32 + 268 + L118)sin2 6 + (s — 53)(s — 34) cos 6]

 

 

\/s—.91\/s—32\/32+26s+wg+b2/c0326—(3—33)(s—34) ( )
= , 4.4
\/s —— sh/s —— 32\/s2 + 263 +468 + b2/cos26 + (s — 33)(s - S4)

 

 

 

where 512 are defined as in (3.9) and 33,4 are defined as in (3.11). Using the substi-
tution defined in (3.14), allows (4.4) to be rewritten into the final frequency-domain

reflection coefficient

 

NS) __ \/3 — 81\/S - 32\/S - 85\/8 - 86 - (8 - S3)(8 - 84)

_ \/3’31\/3‘32\[9—35\/3—36+(3—33)(3-34). (4'5)

 

4.2 Time-Domain Reflection Coefficient

Using the frequency-domain reflection coefficient found in Section 4.1, the time-
domain reflection coefficient can be found using an inverse Laplace transform. In order
to perform the inverse Laplace transform, the frequency-domain reflection coefficient

may be rearranged into a more manageable form. Rationalizing the denominator in

32

(4.5) leads to

 

 

 

S. : lx/S- 51\/3—32\/3 -85\/S—86 - (3-83)(S--S4)l2 _ E
F“) (s—sixs-s2><s—ss><s—s6) —<s—s3>2<s—s4>2 ‘ 0' (4'6)

Noting that

 

 

 

b2
_ _ = _ _ —, 4.7
(S 83)(8 S4) (8 81)(8 82) + COSQQ ( )
allows the denominator to be expanded into
b2 2 2
D = (s — si><s — 32) (s — sue — 321+ .2 — [<3 — 81)(s - s21+b ]
cos 6
1
b (cos26 2) (s 81)(8 32) b
= 62(tan2 6 —1)(s2 + 263 + (.68) — b4
2 2 2 2 ’92
2b (tan 6—1)(S +26S+UJ0 — m)
= 62(tan2 6 —1)(s — 3A)(s - 83), (4.8)
where
5 52 2 b2
= _ :1: _ ___——
SA’B (110 + tan26 — 1
= —6 :1: AA. (4.9)

Using this, (4.6) can be rewritten as

 

 

2 , 2 _ [\/s—sl\F—32\/s—35\/s—36—(.9—33)(.s—34)]2
1) (tan 6 — 1)F(s) — (s _ 5A)(3 _ SB) . (4.10)

Using (4.7), the numerator can be expanded into

N = (s — me - .2). — s5><s — s6)—

33

 

2(3 - 83)(8 - S4)\/S - S1x/8 - S2x/8 - Ssx/S — 86 + (S - 33)2(8 - 84)2
= (S - 81)(8 — 82)(8 - 85)(8 - S6)-

2(8 ~ 81)(8 - 82)\/8 - 51\/3 - .9ng - Ssx/S - 36-

%W%mthfi-mkflwh—m%-mh

232(3 — 31)(s — 32) +194. (4.11)

 

 

At this point, factoring out (3 — 31)(s — 52) allows the numerator to be rewritten in
a form that includes both 91(5) and 92(5) which are defined and inverted previously

in Section 3.2.

An:e—sne—saks—axs—ay+e—sae—sa—

 

 

 

 

 

2 (8 - 85)(S - 86)
2\/s — sh/s — 32\/s — s5\/s — 56 — 2b \/3 _ sh/s _ 82¢8 _ s5¢s _ 36+
b2 2
e—sne—sa+2bl
=<4—axs—a>bma+gxa—2¥ae)+wmer+%fl. (4m)

In the above equation, the 912(3) terms are defined as

 

91(8) = (8 — 81)(8 - 82) - W9 — 31\/3 — 82¢: 85\/8 — 86, (4-13a)
92(8) = (S - 85)(S - 86) - \fS - 81\/S — Szx/E- 85\/§ — 86, (4-13b)

(8 — 85)(8 - 36)

 

 

 

 

 

93(8) = \/s —31\/s—32\/s—s5\/s—35’ (4'13C)
1
94(5) = (S _ S1)(8 _ 82). (4.13d)
By defining
C(S) = (S_SI)(S—S2) (414)

34

the reflection coefficient can be rewritten substituting (4.8), (4.12) and (4.14) into

(4.6) as

(22(tan2 6 —1)F(.s) = C(s) [91(3) + 92(8) — 252g3(8) + 5494(3) + 262]

= C(s)G(s). ' (4.15)

Each term in this expression can be inverted separately to give a final transient

reflection coefficient.
4.2.1 Inversion of the 922(3) Terms

Inversion of the 91(5) and 92(3) terms are solved for previously in Section 3.2.1 and

are given as

- . . . A2 —— 2?
91(0 = 64’ [-AfA§11(/\1t)*11(65t)+ 611201101” 6312951)] — #50),
(4.16)
. A . .. A2 _ A2
92(t) = 67‘” [—616511()~1t)*11(/\5t) + 613120115) + 6512950] + 43-35(0-
(4.17)

In these equations, /\1 and A5 are defined in (3.23) and I72(:r) is defined in (3.51).
The identity from (3.29), redefined below, can be used to perform the inversion

on 93(5) to find 93(t).

 

1 1 1
.__4 "20”")t _ _ . ,
m s + a 6 I0 [2(1) 0):] u(t). (418)

Using this, it is possible to write

b2

cos2 6

 

d2 d
93(t) =( + 26— +423 +

a .11 ).—6t [{IO(A1t)u(t)} * {Ioemuum

35

d2 +26d + 2+——b2 — (t) (419)
= — -— (4} ‘ , . -
dt2 dt 0 cos2 6 93
Taking the first derivative of §3(t) using the product rule yields

$530) = mm) + an), (4.20)

where §g(t) is the exponential term multiplied by the derivative of the term in the

brackets in (4.19). Taking the derivative on the second term in the convolution allows

630) to be written as

630) = e_‘St[/\5{10(61i)11(t)}*{11()\51‘/)U(t)}+10(/\1t)11(t)l- (4-21)

From this, the second derivative can be written as

d d I
d_t2_§3(t) = —(5 [azjfifl] + $93“)

= -5 [—6§3(t) + 1130)] - 6.630) + fig“)

= 5253(1) — 2352(1) + 320:), (4.22)

where §g(t) is the exponential term multiplied by the derivative of the term in the
brackets in (4.21). Using the identity given in (3.39), it can then be shown that
I! 6t )‘5

63(0) = 8’ *2-{10(/\1t)11(()} * {110(651?)+12(65()lu(t)}+ 6111910110) + 5(t) ,
(4.23)

where the derivative was again taken on the second term in the convolution. Com—

bining (4.19), (4.20) and (4.22) allows g3(t) to be written as

930) = 5223(1) - 23—22(1) + 533) — 262530) + 263%) + (4% + 821533)

36

= §§’(t) — (.52 — 413— 82mm.

(4%)

Using the definition for A5 given in (3.23), the final expression for g3(t) is found to

be

93m = 33%) — 4353M
' 2

—6t 55

2

b

 

= e—(St :—/\g {10(A1t)u(t)} * 11(A5t) + 3111(A1t)u(t)] + 6(t).

This took advantage of the Bessel identity defined in (3.43).

The final g-term to invert, 94(3), is given by (4.12) as

1
(8 - 81)(8 - 82)

 

Using partial fraction expansion, this term can be rewritten as

 

 

(s)— 1
9“ ‘<s—s1><s—s2)
_ K1 K2
_3—31 3—32’
where
1
K'=——
1 2A1’
I
A2=—Eq=—K1

 

{1(1()\1t)11(t)}* {112050 —10(/\5t)l11(t)} + 6111910140] + 5(t)

(4%)

(4%)

(4m)

(4.28a)

(4mm

(4%)

Using the identity as defined in (3.27), g4(t) is found to be

gm) <—> 94m = —1— [es1t — esz’] u<t>

2A1
: _e—_6t [CAlt — e—Alt] u(t)
2A1 ’
e—6t .
= —/\— s1nh(A1t)u(t). (4.30)
1

Using (4.16), (4.17), (4.25) and (4.30), C(t) from (4.15) can be written as

0(1) = 2.2-5’ — A’f’xgilolt) *11(A5t) + 111312611) + x§i2(.\51)+

4

522% {10(A1t)u(t)}*11(A5t)—62A111(A1t)u(t)+£\—lsinh(A1t)u(t) . (4.31)

4.2.2 Final expression for C(t)

Just as in the TB case in Chapter 3, there are cases where A1 or A5 can be found
to be imaginary and the modified Bessel functions in (4.31) can be replaced with
ordinary Bessel functions. These cases occur for the same material properties as in
the TE case. Case 1 occurs when (12(2) > 62. In this case both A1 and A5 are purely
imaginary and defined as in (3.54). Then, using the property defined in (3.55) to

replace modified Bessel functions with standard Bessel functions and noting that
sinh(j:r) =jsin(.r.), (4.32)

it can be seen that

0(1) = 2.3—5’ — X55166) * 31(151) + $32511) + £32650—
_ _ . _ _ _ b4 _
62Ag{J0(A1t)u(t)}*J1(A5t)+bQA1J1(A1t)u(t)+fi—sin(A1t)u(t) , (4.33)
1

38

where .ln(;1:)is defined in (3.57).

Case 2 occurs when 1123+ 82 < 62. In this case, both A1 and A5 are purely real and
the same expression from (4.31) can be used. Case 3 occurs when —82 < 10(2) —62 < 0.
In this case, A1 is purely real and A5 is purely imaginary. A5 is defined as it was in

(3.54) and the expression for C(t) can be written as

em = 26‘” 452mm * 31652:) + Aii2(A1t) + 4232650—

4

521?, {10(A1t)u(t)} * 31(A5t) - 52A111(A1t)u(t) + $1- sinh(A1t)u(t) . (4.34)

4.2.3 Inversion of the C(s) Term

The C(s) term is defined in (4.14) as

 

0(3) = (s — s1><s — s2)

(8 - SA)(8 - 3B). (435)

Through algebraic manipulation and partial fraction expansion, it can be rewritten

into an invertible form as

_ (8-81)(8-82)- (S-SA)(8-SB)
0(3) ‘ H (s — we — s3)

 

 

 

 

 

2 2 2 2 b2
(s +26s+w0)—(s +26s+w0—m)
(S-SA)(S -SB)
b2
(tan2 6 —— 1)(s — 3A)(s —- 83)
b2 1 1
:1 — . 4.36
+2AA(tan26—1)[s—3A s—sB ( )

The inversion of C (3) depends on the values of s A and 83. If 0J8 > b2/(tan2 6 — 1),

C (s) can be inverted using the identity (3.27) into

b2
2AA(tan2 6 — 1)

 

C(s) 4—4 C(t) = 6(t) + [.234t — eth]u(t)

39

026—6t
AA(tan2 6 — 1)

 

= 6(t) + sinh(AAt)u(t). (4.37)

If 1128 > 62 + b2/(tan2 6 — 1), A A is purely imaginary and defined as

 

. — b2

In this case, using (4.32) C(t) can be simplified into

b2e—6t
AA(tan2 6 — 1)

sin(AAt)u(t). (4.39)

 

0(1) = 4(1) +

If 0J3 > 62 + 122/(tan2 6 — 1), AA is real and less than 6. In this case, (4.37) can be
used directly. If Lug < bz/(tan26 — 1), A A is real and greater than 6. In this case

Re{sA} > O, which requires the following inverse Laplace transform identity to be

used [14]:
1

+——+-—efi‘u —- . .
S_fi t( 0 @4m

 

Using this identity and the identity defined previously in (3.27), C (t) can be written

as

 

 

528%” A t A 1
0(1) = 5(1) + 2221(141126— 1) [—e A u(—t) — e- A 11(1)]
b28—6t
==5e) e-AAVL (441)

— 2AA(tan26 — I)

The crucial issue to address with this decomposition is the fact that this term is
non-causal (ie, has non-zero value for t < 0). Further examination of the terms in
question show that this case can only occur at angles greater than 45°and between
some upper angle limit determined by the material properties. It appears that this
non-causality is due to a problem with the model for 6r. A true plane wave is a

non-physical construct that is convenient in many cases, but does not exist. When

40

an infinite plane wave is obliquely incident on an infinite half-space, it can be seen
that there is no point in time when the wave is not intersecting the half-space. This
has the potential to cause problems with causality such as in this case. Ultimately,
the non-causal term exists in this case solely due to this fact.

Using the final expressions for G (t) and C (t), the final expression for the reflection

coefficient can be written as

 

r(1) — 0(1) 4 0(1). (4.42)

4.3 Numerical Results

In order to validate the expressions derived in the previous section, the time-domain
reflection coefficient is evaluated numerically in Fortran, and then compared to the
inverse FF T of the frequency domain reflection coefficient given in (4.5). The Fortran
code is included in Appendix B. The inverse FFT was done using WaveCalc. The
frequency-domain data was zero-padded up to the maximum limit allowed by Wave-
Calc, 32,768, before the inverse F FT was taken. Since the signal had already decayed
to zero by this point, no windowing was necessary. A set of parameters corresponding
to each of the three possible cases for C(t) is used. The angle is then adjusted for
the case 1 parameters to produce results that contain a non-causal contribution. The
number of points and step size for the results varied between the different cases.
The first set of parameters is the same as those chosen by Brillouin: 020 = 4.0 x
1016 3‘1, b2 = 20.0 x 1032 3‘2, 6 = 0.28 x 1016 3‘1. This choice of parameters
corresponds to case 1. When computing the numerical results, for the frequency-
domain data 16,384 frequency points were calculated with a step size of 8,000 GHz.
In the time-domain, 4,096 points were calculated with a step size of 1 x 10‘9 ns. Using
6 = 30°, (4.33) has been plotted in Figure 4.1 and compared to the inverse FFT. The

results show excellent agreement. Since this function includes only standard Bessel

41

functions, which are highly oscillatory, and no modified Bessel functions, which are
not oscillatory, the waveform is highly oscillatory and only lightly damped.

The next choice of parameters is: (.00 = 2.0 x 1015 3'1, 52 = 20.0 x 1029 3’2, 6 =
0.28 x 1016 3‘1, which corresponds to case 2, (4.31). When computing the numerical
results, for the frequency-domain data 32,768 frequency points were calculated with a
step size of 400 GHz. In the time-domain, 4,096 points were calculated with a step size
of 1 x 10’9 ns. The results are shown in Figure 4.2. Again, the closed-form expression
and the inverse FFT compare well. For this choice of parameters 62 > L08 + 82, and
the resulting waveform is overdamped, showing no oscillatory behavior and only a
single negative peak. Since the expression for case 2 only involves modified Bessel
functions, which do not have the oscillatory behavior of ordinary Bessel functions, this
observed behavior is easily predicted from the mathematical form of the expression.

The next choice of parameters, which corresponds to case 3, is: 100 = 2.0 x
1015 3‘1, 1)2 = 20.0 x 1032 3’2, 6 = 0.28 x 1016 3‘1, which corresponds to case 3,
(4.34). When computing the numerical results, for the frequency-domain data 32,768
frequency points were calculated with a step size of 6,000 GHz. In the time-domain,
16,384 points were calculated with a step size of 1 x 10—10 ns. The analytic expression
again matches the inverse FFT, as seen in Figure 4.3. As expected, since 6 > (.00, but
62 < (.08 + B2, there is more damping and less oscillation than with case 1, but more
oscillation than with case 2. Here the expression for the reflection coeflicient has a
combination of ordinary and modified Bessel functions.

To examine the non-causal result which is possible for the C (1) term, Brillouin‘s
choice of parameters are used again, but at an angle of 50°. This causes the C(t)
term to be non-causal as shown in (4.41). When computing the numerical results, for
the frequency-domain data 16,384 frequency points were calculated with a step size
of 8,000 GHz. In the time-domain, 4,096 points were calculated with a step size of

1 x 10’9 ms with a starting time value of —4.096 x 10‘6 ns. As seen in Figure 4.4, the

42

expression again matches well with the inverse FFT for t > 0. For t < 0, non-causal
term in the analytic expression is extremely small. Unfortunately, the inverse FFT
does not have the necessary resolution to capture the small non-causal portion of the
reflection coefficient, so it cannot be compared to the analytic expression.

In order to have a better insight into the non-causal term, the delta function
term was subtracted from the non-causal C-term and plotted in Figure 4.5. It can
immediately be seen that this term is extremely small when compared to the final
transient field. In order to examine the non-causality even further, this term is
convolved with G (t), as it is in the final expression for the reflection coefficient, and
plotted in Figure 4.6. Just as in Figure 4.4, the non-causal term is not even visible
in this plot. Taking a closer View of this plot in Figure 4.7 it can be seen that the
non-causal portion of the signal is around 3-orders of magnitude less than the causal

portion.

43

x 1015 Case 1

Reflection Coefficient
O

 

 

 

 

 

 

 

-2
-4 r.
-6 - .
Closed Form Expressuon
---><- -- Inverse FFT
_8 4 1 m 1 J
0 0.2 0.4 0.6 0.8 1

Time (is)

Figure 4.1. Time-domain reflection coefficient with incidence angle 6 = 30° and mate-
rial parameters wo = 4.0 x 10163—1, b2 = 20.0 x 10323’2, 6 = 0.28 x 10163—1. This
choice of parameters corresponds to case 1, Eq. (4.33).

44

x1013 Case 2

Reflection Coefficient

 

 

 

 

 

 

 

 

-4.5 - Closed Form Expression
- - -><- ~ Inverse FFT
_5 a 1 1 1 1
O 1 2 3 4 5

Time (fs)

Figure 4.2. Time-domain reflection coefficient with incidence angle 6 = 30° and mate-
rial parameters 110 = 2.0 x 10155-1, 52 = 20.0 x 10293-2, 6 = 0.28 x 10161-1. This
choice of parameters corresponds to case 2, Eq. (4.31).

45

x1015 Case 3

Reflection Coefficient

_L
o
r

 

I

-12

 

Closed Form Expression
- - -><~ -- inverse FFT

_14 l l l l
0 0.2 0.4 0.6 0.8 1

Time (is)

 

 

 

 

Figure 4.3. Time-domain reflection coefficient with at an angle of 6 = 30° with param-
eter choices of 100 = 2.0 x 10153’1, b2 = 20.0 x 10323-2, 6 = 0.28 x 10163-1. This
choice of parameters corresponds to case 3, Eq. (4.34).

46

x 1016 Case 1 - Non-Causal

I

0.5

Reflection Coefficient

 

 

 

 

Closed Form Expression
- --x1~- Inverse FFT

 

 

 

 

l J

-0.2 O 0.2 0.4 0.6 0.8 1
Time (is)

Figure 4.4. Time—domain reflection coefficient with incidence angle 6 = 50° and mate-
rial parameters 100 = 4.0 x 10163—1, ()2 = 20.0 x 10323—2, 6 = 0.28 x 10163-1. This
choice of parameters corresponds to case 1, Eq. (4.33) with non-causal term from Eq.
(4.41).

47

Reflection Coefficient

x 10’18 Non-Causal C-term

 

 

 

 

 

 

-0.5 0 0.5
Time (is)

Figure 4.5. Non—causal C(t) — 6(t) for case 1 results.

48

Reflection Coefficient

x 1015 Non-Causal C term convolved with 6 term

 

 

 

 

 

 

 

 

1 L L

-0.5 0 0.5
Time (is)

Figure 4.6. Non-causal [C(t) — 6(t)] =1: C(t) for case 1 results.

49

x 10‘? Non-Causal C term convolved with G term

Reflection Coefficient

 

 

_5 1 1 1 L
-0.05 -0.04 -0.03 -0.02 -0.01 O 0.01
Time (fs)

Figure 4.7. Close-up of the non-causal [C(t) — 6(1)] * C(t) for case 1 results.

50

CHAPTER 5

TE PLANE WAVE REFLECTION FOR THE GENERAL CASE

In Chapter 3, the transient reflection of a TB. plane-wave from the Lorentz-medium
half-space is examined for a Lorentz-medium simplified to the optical case. In this
chapter, the reflection coefficient will be examined for the most general case. As in
the Chapter 3, a steady-state TE—polarized plane wave of frequency (.0 is obliquely
incident on an interface separating free space (region 1) from a homogeneous Lorentz
medium (region 2). The angle of incidence 6 is measured from the normal to the
interface. This geometry can be seen in Figure 2.2. As shown in Chapter 2, the

reflection coefficient for a plane wave from a material half—space can be expressed as
(5.1)

For a TE—polarized plane wave, the impedance of the incident wave is Z0 = 710/ cos 6
and the wave impedance of the transmitted wave is given by

k(w)n(w)

21(9)): kz(LU) 1

(5.2)

where the terms k2, r), and k are defined in (3.3). The relative permittivity is the
same as given in (3.4). For the general case, the simplifications made in Chapter 3
cannot be made, but the relative permittivity can be rewritten as

b2
(08 — 102 + 23106,

 

where b is the plasma frequency of the medium.

51

5.1 Laplace Domain Representation

Equation (5.3) can be represented in the Laplace domain as

b2
82 + 268 +103.

 

67(8) = 600 +

Using this, the wave number in the Lorentz medium can then be calculated as

11s) = 43052

 

b2
32 + 263 + L123

2 | l 2' 2
= — 's,/1€ C . 5.5
3 (00°[/ s2+263+,.8 ( )

 

= —js,/,u€0 (.30 +

 

 

For a TE—polarized plane wave, the Laplace domain reflection coefficient can then be
written as
= #1630 _ k:
kirk/:0 + k:
,11,~ cos6 — \/(k/k0)2 — sin.2 6
— 112 cos6 + \/(k/k0)2 — sin2 6

_ 4. case — (Meme - s3><s — 141/[<3 — sl><s - s21] — sin2 0

F(s)

 

 

 

 

 

 

 

, (5.6)
112 cos6 + (/11,-coo(s — 53)(3 — S4)[(8 — 31)(s — 52)] -— sin2 6
where 31,2 are the same as defined in (3.9) and 3314 are defined as
(,2
834 = —6 :f: (52 - (41(2) — —, (5.7)
500
which is derived from the fact that
2 2 b2
3 +26s+c00+— = (3—33)(s—34). (5.8)

600

52

 

Multiplying the numerator and denominator in (5.6) by \/(s — 31)(s — 5'2) and ex-

panding (s — sl)(s — 32) and (3 — 33)(s — 84) leads to

 

 

_,..1051¢(. — sl><s — .2) — was — s3><s —- s4) — sin26<s — s1><s — 32)

1(5)
,ur cos6/(s — 81)(8 — 52) + (Al/reeds — 33)(s — 34) — sin2 6(3 — 31)(s — 32)

 

 

 

 

= [111‘ cos6\/(s — 31)(s — 52) — (32012-600 — sin2 6) + 265012600 — sin2 6)+

 

1 2
108011-600 — sin2 6) + 11,152) / ]/[/1,- cos 6\/(s — 31)(s — 82)+

 

 

 

 

 

 

1 2
(32012100 — sin2 6) + 268(nrcoo — sin2 6) + 108012600 — sin2 6) + 14202) /
(5.9)
For convenience, define
K2 = 11.2600 — sin2 6. (5.10)
This allows (5.9) to be written as
11.2 cos 6\/(s — 31)(s — 5'2) — K\/s2 + 268 + 1.03 + p.7b2/K2
s = . (5.11)
#7: cos6\/(s — 31)(s —— 82) + K\/32 + 263 + (08 + urb2/K2
It is then convenient to make the substitution
2 - 2 b2
3 + 263 + (4)0 + (”I-{3 = (s — 35)(s — 35), (5.12)
where
2 b2
35,6 = —6 :1: 62 - 0J0 - [141-F. (5.13)

Substituting (5.12) into (5.11) allows the reflection coefficient in the Laplace domain

to be rewritten as

 

 

F( [l~,~COSg\/S—Sl\/5—SQ—Kx/S—55\/S_36
s = .
11,. cos 6\/s — sM/s — 32 + K\/S - S5\/8 — 86

 

(5.14)

 

 

53

5.2 Time-Domain Reflection Coefficient

Using the frequency—domain reflection coefficient found in Section 5.1, the time-
domain reflection coefficient can be found using an inverse Laplace transform. In
order to perform the inverse Laplace transform, the frequency-domain reflection co-
efficient. may be rearranged into a more manageable form. Before this can be done
however, it must be noted that as s —+ 00, the reflection coefficient does not tend

to zero. Instead it can be seen that

11., cos 6 — K

as s—>oo Fs —.
()—) p,.cos6+K

(5.15)

This term results in a delta function contribution in the final time-domain reflection
coefficient. In order to perform the inverse Laplace transform, this term, defined as

Foo, needs to be subtracted from the total reflection coefficient. This leads to

[(9): F(s) — FOO
_11.4,-cos6\/s — .91\/s — 82 — K(/s — s5\/s — 35 hr cos6 — K
— 11,- cos 6\/s — sl\/s — 32 + K\/s — 35¢s — 35 ,u,- cos6 + K
2Kur cos6[\/s — sh/s — 32 — \/s — 35\fs — 35]

= . 5.16
(112- cos6 + K) (pr cos 6\/s — sh/s — 52 + K\/s — 55\/s — 36) ( )

 

 

 

 

 

 

 

 

 

 

Rationalizing the denominator leads to

 

 

 

F(.s) = 2K,u,- cos6 (\[S —— 81\/S — 32 — J3 — s5\/s — 85) (11.1- cos6\/s — 31\/s — 32—

 

K\/s — s5\/s — 5'5) ] / [(112.cos6 + K) [11.3 cos2 6(3 — 31)(s — 32)—
4'20 - s5)<s — 411]]

N
=F1

5, (5.17)

54

where

2K/1rcos6
F =———-. '.1
l prcos6+K (a 8)

Expanding the terms in the denominator and examining them in further detail, it can
be seen that
”1.1)2 )

D = 11? cos2 6(32 + 26.9 + 1.03) — K2(s2 + 263 + 103 + 72—-

: 32(113cos2 6 — K2) + 263(113cos2 6 — K2) + 108(113cos2 6 — K2) — #252

2
2 2 '2 2 2 #75?
= 1.cos 6—11 5 +26s+" —
(l7 )( “/0 )12cosz6—K2)

 

= (#311,052 6 — K2) (82 + 268 + ’72) , (5'19)

with the definition

 

1 = w _ . .
A, 0 )rgcosz6—K2
At this point, it is convenient to define
32+26s+72 = (s—SC)(s—SD), (5.21)
where
5C,D = —6i V62 - ’72
= —6 :1: A0. (5.22)

Substituting (5.21) into (5.19) allows the denominator to be written in its final form
as

D = (11? cos2 6 — K2)(s — SC)(s — SD). (5.23)

Expanding the numerator and examining it in further detail leads to

 

 

N = )1,- cos6 [(s — 81)(S — 32) — \/s — shfls — 82\/S — s5\/s — 36] +
K [(8 — «85)(8 — 86) - \/8 — 31\/3 — 82x/8 - 85\/8 — 36]

= 11,. cos 691(3) + Kgg(s). (5.24)

 

Applying the following definition

2K112-cos6
(11,. cos6 + K)2()u.r cos6 — K)’

 

F = (5.25)

allows the Laplace domain reflection coefficient to be rewritten by substituting (5.23)

and (5.24) into (5.17) as

~ _ #rC05691(3)+K92(3)

 

 

F .9 — F
( ) (s — 80)(8 — 81)
= FC(s)G(s), (5.26)
where it is defined that
0(1) — 1 (5 27a)
(8—80)(8-SD)’ '
C(s) = 11,. cos 6g1(3) + Kgg(s). (5.27b)

5.2.1 Inversion of the C(s) Term

Both terms 91(3) and 92(3) in 0(3) were inverted previously in Section 3.2.1. The

solutions are given as

2 2
41—45
2

 

91(() = 8’6] [‘6f45i10‘1t) * i1(/\5t) + 61122910 + 6312950] - (W),

56

2 2
)‘1_)‘5

 

92(1) = e-6t [5951511) 111(151) + 11112011) + 2212551)] + 6(1).
(5.29)
In these equations, A1 and A5 are defined from
812 = —(5 :1: A1, (5.30)
55,6 = —6 :1: A5, (5.31)

and 122(2) is defined in (3.51). The final expression for C(t) can then be given as

0(1) = (11,. cos6 + roe—5* [411311510 111051) + A¥i2(,\11)+ 2312051)] —

“Tb?

2K2

 

(112- cos6 — K)6(t). (5.32)

Just as in the optical TE case in Chapter 3, there are cases where A1 or A5 can be
found to be imaginary and the modified Bessel functions in (5.32) can be replaced
with ordinary Bessel functions. Case 1 occurs when 103 > 62. In this case both A1

and A5 are purely imaginary and defined as

_ _ ..b2
—jA1=A1=(/1.03—62, —jA5=A5=\/w3+u7;—2-—62. (5.33)

Then, using the identity given in (3.55) to replace modified Bessel functions with

 

ordinary Bessel functions and the identity given in (4.32), (5.32) can be expressed as

0(1) = ()1. cos a + 1101-8 [511331510 13,1151) + 11.12511) + 1232651)
2

£5728 cos6 — K)6(t), (5.34)

b—J

where Jn(;r) was defined in (3.57).

57

Case 2 occurs when wg+prb2/K2 < 62. In this case, both A1 and A5 are purely real
and the same expression from (5.32) can be used. Case 3 occurs when —p.rb2 / K 2 <
Lug — 62 < O. In this case, A1 is purely real and A5 is purely imaginary. A5 is defined

as it was in (5.33) and the expression for C(t) can be written as

C(t) = ()1... c056 + roe-5t [A¥A§I1(A1t)*31(A5t)+ A?I2(A1t) + $332655] —

[tr ()2

WULTCOSQ — K)6(t). (5.35)

5.2.2 Inversion of the C(s) Term

Using partial fraction expansion, C (s) can be rearranged into

1
(8 - Sc)(8 - 3D)

1 1
2A0 s — SC 3 — 81)

As in previous cases, the values of SC and s D affect the inversion of C (5) If 72 > 0,

C(8) =

 

 

 

then Re{sc,D} < 0 and the inversion identity (3.27) can be used. Using this identity,
C(t) is found to be

—6t
= :3— Act _ —Act
C(t) 2A0 (e e )u(t)
e—ét .
= —— smh(ACt)u(t). (5.37)
*0

If 72 < (52, AC is purely real and (5.37) can be used directly. If 72 > 62, AC is purely

imaginary and can be defined as

—J')\C = XC = 72 - 52. (5.38)

58

In this case, using the identity given in (4.32), (5.37) can be rewritten into

e—(ft _
C(t) = T sin(ACt)u(t). (5.39)

If '72 < 0, AC is again purely real, but in this case, Re{sc} > O, which requires the

identity given in (4.40) to be used. In this case, C (t) is given as

—6t
. :_C____ ACt _ —ACt
C(t) 2A0 (e u( t)+e u(t))
6—6t

= —— ”‘0'”. 5.40
2M6 ( )

As in Chapter 4 for the TM optical case, there is a problem with causality in this
term, however an interesting result can be observed by examining the 72 term. In

order for ")2 < 0 to occur, the following relationship must hold:

2
2 llrb
w < , .
0 #2 cos2 (9 — K 2

 

(5.41)

Examining the denominator on the right side of the equation in closer detail reveals
that as long as #r < 6.30, the denominator is negative at any angle. Since this is not
the optical case, it holds that 600 > 1. The Lorentz model is generally used for a
dielectric material, so it is expected that pr < 600 will hold in general. This implies

that the causality issue is not a concern in this case.

5.3 Numerical Results

In order to validate the expressions derived in the previous section, the time-domain
reflection coefficient is evaluated numerically in Fortran, and then compared to the
inverse F FT of the frequency-domain reflection coefficient given in (5.16). The For-

tran code is included in Appendix C. The inverse F FT was done using WaveCalc.

59

The frequency-domain data was zero-padded up to the maximum limit allowed by
WaveCalc, 32,768, before the inverse FFT was taken. Since the signal had already
decayed to zero by this point, no windowing was necessary. A set of parameters
corresponding to each of the three possible cases is used.

The first set of parameters chosen are: wO = 4.0 x 1011 3‘1, b = 6.24 x 1011 5‘1,
6 = 2.5 x 1010 3‘1, 600 = 2, M = 1. This choice of parameters corresponds to case
1. When computing the numerical results, for the frequency-domain data 4,096 fre-
quency points were calculated with a step size of 400 MHz. In the time-domain, 4,096
points were calculated with a step size of 1 x 10’4 ns. Using 6 = 30°, (5.34) has been
plotted in Figure 5.1 and compared to the inverse FFT. The results show excellent
agreement. Since this function includes only standard Bessel functions, which are
highly oscillatory, and no modified Bessel functions, which are not oscillatory, the
waveform is highly oscillatory and only lightly damped.

The next choice of parameters is: wo = 4.0 x 1010 3'1, b = 6.24 x 1010 3'1,
6 = 2.5 x 1011 3‘1, 600 = 2, Hr = 1, which corresponds to case 2, (5.32). When
computing the numerical results, for the frequency-domain data 32,768 frequency
points were calculated with a step size of 4 MHz. In the time-domain, 8,192 points
were calculated with a step size of 2 x 10‘4 ns. The results are shown in Figure
5.2. Again, the closed-form expression and the inverse FFT compare well. For this
choice of parameters tug + prb‘? / K 2 < (52, and the resulting waveform is overdamped,
showing no oscillatory behavior and only a single negative peak. Since the expression
for case 2 only involves modified Bessel functions, which do not have the oscillatory
behavior of ordinary Bessel functions, this observed behavior is easily predicted from
the mathematical form of the expression.

The final choice of parameters, which corresponds to case 3, is: too = 4.0 x 109 3‘1,
b = 6.24 x1011 8'1, 6 = 2.5 X 1010 3‘1, 600 = 2, Hr = 1, which corresponds to case 3,

(5.35). When computing the numerical results, for the frequency-domain data 8,192

60

 

 

 

frequency points were calculated with a step size of 200 MHz. In the time—domain,
4,096 points were calculated with a step size of 1 x 10‘4 ns. The analytic expression
again matches the inverse F F T, as seen in Figure 5.3. As expected, since 62 > tag, but
62 < 0.23 + urbz/K2, there is more damping and less oscillation than with case 1, but
more oscillation than with case 2. Here the expression for the reflection coefficient

has a combination of ordinary and modified Bessel functions.

61

X 1011 Case 1

 

 

 

 

 

 

 

1 r
0.5 ~
‘5
.92
£3 0 _ 5...; .......................
d)
o
0
C
.9
E -0.5
d)
o:
-1
Closed Form Expression
---X- -- Inverse FFT
_15 1 1 1 4 I
0 0.02 0.04 0.06 0.08 0.1

Time (ns)

Figure 5.1. Time-domain reflection coefficient with incidence angle 0 = 30° and mate-
rial parameters wo = 4.0 x 10113‘1, b = 6.24 x 10113-1, 6 = 2.5 x 10103—1, 600 = 2,
m = 1. This choice of parameters corresponds to case 1, Eq. (5.34).

62

x 103 Case 2

‘(
X

Reflection Coefficient

 

 

 

Closed Form Expression
---><' -- Inverse FFT

_15 1 1 4L 1 1

0.2 0.4 0.6 0.8 1
Time (ns)

 

 

 

 

Figure 5.2. Time-domain reflection coeflicient with incidence angle 0 = 30° and mate-
rial parameters wO = 4.0 x 10103-1, b = 6.24 x 10103—1, 6 = 2.5 x 10113-1, 600 = 2,
11? = 1. This choice of parameters corresponds to case 2, Eq. (5.32).

63

x 101° Case 3

 

 

 

 

 

 

 

5 ..
O __ A A
E
.9
g -5
a)
O
0
C
.9
3 -10 ~
C
Q)
Q:
-15 _
Closed Form Expression
~~ ~><- - ‘ Inverse FFT
_20 1 1 1 1 1
0 0.02 0.04 0.06 0.08 0.1

Time (ns)

Figure 5.3. Time-domain reflection coeflicient with at an angle of 9 = 30° with param-
eter choices of (.20 = 4.0 x 1093-1, b = 6.24 x 10113—1, 6 = 2.5 x 10103—1, 600 = 2,
11.7. = 1. This choice of parameters corresponds to case 3, Eq. (5.35).

64

CHAPTER 6

TM PLANE WAVE REFLECTION FOR THE GENERAL CASE

In Chapter 4, the transient reflection of a TM plane-wave from a Lorentz-medium
half-space is examined for a Lorentz-medium simplified to the optical case. In this
chapter, the reflection coefficient will be examined for the most general case. As in
the Chapter 4, a steady-state TM-polarized plane wave of frequency w is obliquely
incident on an interface separating free space (region 1) from a homogeneous Lorentz
medium (region 2). The angle of incidence 0 is measured from the normal to the
interface. This geometry can be seen in Figure 2.3. As shown in Chapter 2, the
reflection coefficient for a plane wave from a material half-space can be expressed as

_ anW) — 20

(6.1)

For a TM—polarized plane wave, the impedance of the incident wave is no c089 and

the wave impedance of the transmitted wave is given by

new)

2“”) : k<w>

(6.2)

where the terms 1%,”, and k are defined in (3.3). The relative permittivity is the

same as given in (5.3).

6.1 Laplace Domain Representation

The Laplace domain reflection coefficient for a TM-polarized plane wave can be writ-

ten as

65

 

\/(k/k0)2 — sin2 6 -— €7- cos6

 

 

\/(k/k0)2 — sin2 6 + 6p cos 6

_ \/ urer — sin2 6 — 6,- cos6
\/ 111.6,. — sin2 6 + 6r cos 6

= [ill-Mods — 83)(8 - S4)/l(8 - 81)(S - 82)] - 31112 6ill/2"

 

coo cos6(8 — 33)(s — S4)/[(8 — 31)(3 ._ 32)]]/
[inroads — s3)(s — S4)/[(S - 81)(s _ 32)] _ sin? all/2+

€30 COS 9(8 - S3)(8 - S4)/[(8 - 81)(8 - 82“], (5-3)

where 312 are the same as defined in (3.9) and 33.4 are the same as defined in (5.7).

 

lV‘Iultiplying the numerator and denominator by \/(s — sl)(s — 82) and expanding

(s — 31)(s — 32) and (s — 53)(s — 34) leads to

 

re) = [v - slvs — same — s3>(s — s.) _ sin2 9(3 _ 3,)(3 _ 3251/2-

 

600(8 - 33)(s — 54) cos 6] / [$9 — 31\/8 — 32l/lr€oo(3 — 33)(3 - 34)—

sin2 6(s — 31)(s — 52))”2 + 600(8 —- 33)(s — .94) cos 6]

 

= [Vs — sn/s' —- 32 [(preoo — sin2 6)s2 + 263(preoo — sin2 6)+

2 - 2 2 ”2
w0(/1.,~€oo — Sin 6) + ,11.,.b ] — (.00 cos 6(3 — S3)(S — 34) /

 

[\/s —— 31\/s - .92 [wrap — sin2 6)s2 + 268(ureoo — sin2 6)+

1/2
wgmreoo — sin2 6) + #752] + 600 cos 6(3 — 33)(s — 34)]. (6.4)

66

Substituting K2 as defined in (5.10) allows (6. 4) to be written as

 

 

K\/s — sl\/s — 82\/32 +268 ~l-wf)Z +urb2/K2 — coo cos6(s — 33)(s - 34)

 

 

NS) =

 

K\/s — sh/s — 52\/32 + 263 + (11(2) + per/K2 + £00 cos 6(3 — 33)(s — 34)

(6.5)

Substituting (5.12) into (6.5) allows the final form of the reflection coefficient to be

rewritten as

 

F(.s)- K\/s — ssh/s — 52¢s — 55\/s — .95 — coocos6(s — 33)(3— 34)

A \/s — sh/s — 82\/8 — s5\/s — 86 + 60C cos6(s — S3)(8 — 84) (6'6)

 

6.2 Time-Domain Reflection Coefficient

Using the frequency-domain reflection coefficient found in Section 6.1, the time-
dornain reflection coefficient can be found using an inverse Laplace transform. In
order to perform the inverse Laplace transform, the frequency—domain reflection co-
eflicient may be rearranged into a more manageable form. Before this can be done
however, it must be noted that as s ——> 00, the reflection coefficient does not tend

to zero. Instead it can be seen that

K— eoocos6
K+eoocos6’

 

as s ——> oo F(s) ——>

(6.7)

which produces a delta function in the final time—domain reflection coefficient. In
order to perform the inverse Laplace transform, this term, defined as FOO, needs to

be subtracted from the total reflection coefficient. This leads to

F(s) = F(s) — Foo

_ Kfi— sh/s —— sub/s — sh/s —- 52 — coo cos6(s — 3:3)(3 — 34) K — coo cos6
_ K\/s — 31\/3 — 52¢s — sh/s — 82 + 600 cos6(s — 53)(s — 54) K + £00 cos6

 

 

 

 

67

 

 

2Keoo cos6 [\/s — sl\/s — .9ng — S5\[S — 36 — (s — 53)(s — 34)]
(K + 600 cos 6) [K\/s — sh/s — 52\/s — 55¢}; — 35 + 600 cos6(s — 33)(s — 84)].

 

 

 

 

(6.8)
Rationalizing the denominator allows the denominator to be written as
D = (K + 600 cos 6) [K2(s -— 31)(s — 32)(s — s5)(s — 36)—
630 cos? 6(3 — .93)2(s — 34)2] , (6.9)
and the numerator to be written as
N = 2Kcoc. cos6[I\"(s — 31)(s — 32)(s — s5)(s — 35)—
(K + (3C, cos 6)(s — 33)(s — .34)\/s - sn/s — 52¢s -— 35\/s —— 35+
(00 cos 6(s — 33)2(s — 54)?) . (6.10)
Using the definition
b2( 1. — 26 .cos2 6)
2 l 7‘ so
= 6.11
7 K2 — ego cos? 6 ( )
the denominator term can be expanded and written as
D = (K + 600 cos 6)(K2 -— (30 cos 6) [s4 + 4633 + (21128 + 462 + 72)s2+
(14 cos2 6
4612+262s+ 14+ '22-— . 6.12
( WO 7 )5 “0 “/07 K2 -— ego cos26 ( )

The next step is to factor the denominator term. This is accomplished by dividing
the quartic function in the denominator by s2 + 263 + 111(2) + 72/ 2 + x using a long

division technique. The remainder term is set equal to zero and solved for x. This

68

allows the denominator to be written as

,2
D = (K+€oocos6)2(K—coocos6) (32+263+w5 + %— —X) x

2
(32 + 263 + to?) + 7? + x) , (6.13)

where X is defined as

 

b20
= . . . 6.14
X 2(K2 -— 65C cos2 6) ‘ ( )
and a is defined as
02 = )1? + sin2(26). (6.15)

It is then convenient to redefine the denominator as
D = (K + 600 cos 6)2(K - coo cos 6)(s — 8E)(S — 3F)(s — 30)(s — 3H), (6.16)

where

2
(s—sE)(s—3F)=S2+263+w8+12-—X, (6.1721)
,72
(s—sG)(.s—3H) =32+263+w3+3+x (6.17b)

The roots can then be found to be

72
.SE,F=—6:t 62—(wgi-3—x)

 

 

= —5iAE, (6.188.)
72
SG,H="6i 62— (w8+-é—+x)
= —-5 :1: A0. (6.1813)

69

The numerator also needs to be rearranged in order for the inverse Laplace trans-

form to be taken. Noting that

b2
(8-83)(8-84) = (8-81)(8-52)+—1 (6-19)

600

it is possible to rewrite the numerator term as

N = 21(600 cos 6(3 — 31)(s — 52) [60C cos 6[(s — 31)(s — 52)—

 

 

 

x/S - Six/8 - 82x/8 — 85\/8 - Sol + KKS — 85)(8 - 86)—

 

 

 

 

 

 

b2
\/3 — ssh/.9 — 82\/S —— s5\/s — .95] — -€—(K + (.00 c0s6)(s — s5)(s - 36)x
00
1 b4cos6 1 2
+ +2b cos6
\/'3’51\/3_52\/5_55\/3—50l foo (8—81)(8—82)

= 2K600 cos 6(3 —— sl)(s —- 32) [600 cos 691(3) + Kgg(s)—

b2 b4 9
—<K + co. sesame) +

6m 6C30

 

94(3) + 21)2 cos 6]

= 2K600 cos 6(3 — 81)(8 — 32)G(s). (6.20)

The gn(s) terms are the same as in the previous cases. The reflection coefficient in

the frequency domain can then be written as

2KeoO cos 6(3 — 31)(s — 32)G(s)

 

 

f . : . 6.21
(5) (K + 600 cos 6)2(K — coo cos6)(s — 3E)(s — 3F)(s — SG)(S — 5H) ( )
Making the convenient substitutions
2Keoo cos 6
F = , 6.22
2 (K + 600 cos 6)2(K - €00 cos 6) ( a)
C(s) = (3 — 31x3 _ 32) (6.22b)

 

(8 - 5E)(3 - 8F)(8 - SG)(3 — SH)’

70

allows (6.21) to be rewritten as
fig) = FQC(8)G(3). (6.23)

The Laplace transform of both the C (3) and the G (3) term can be taken separately

and the time-domain reflection coefficient is simply a convolution of those terms.
6.2.1 Inversion of the C(s) Term

As shown in (6.20), C(s) is given as

b2 b4 :6
0(8) = (00 cos 691(3) + K92(3) —— —(K + 600 cos 6)93(s) + C08

600 CCX)

 

94(3) + 2b2 cos6 .

(6.24)

All of the terms in 0(3) are solved in previous sections. The inversion of the 91(3)
and 92(3) terms is done in Section 3.2.1 and the inversion of the 93(3) and 94(3) is

done in Section 4.2. These terms are given as

A A A .. A2 _ /\2
91(1) = 6‘” {-AfA511(A1t)*11(A5t) + Ailgolt) + £12056] — #2—566).
(6.25a)
6t 2 2* * 3* 3* Af _ ”)2
92(0 = (f— [*A1A511(/\1t) *110‘50 + Ail2()\1t)+ (\512050) + 74350).
(6.25b)
93(t) = e76t [—Ag{ I()(A1t)u(t) }*11(A5t) + A1I1(A1t)u(t)] + (5(1‘), (6.25c)
e—(St
94(t) 2 TI sinh(A1t)u(t). (6.25d)

In these equations, A1 and A5 are defined as in (5.30) and In(:r) is defined as in (3.51).

The final expression for C(t) can then be given as

. . b2 .
C(t) = e_6t{(1{ + 500 COS 9) [~A¥A§11(A1t)*11()\5t) + C—Ag {I0(A1t)u(t)} *11(A5t)—
oo

71

 

()2 3,. 3,. ()4 C086
—A111(A1t)l_l(t) + A'112(/\1t) + A512(x\5t):| + sinh(A1t)u(t)} +
€x EmAl
A _2K3
b2(K -- 60c cos 6) (Mi—2T) 6(t). (6.26)
00

Just as in the optical TM case in Chapter 4, there are cases where A1 or A5 can be
found to be imaginary and the modified Bessel functions in (6.26) can be replaced
with ordinary Bessel functions. Case 1 occurs when (12(2) > (52. In this case both A1
and A5 are purely imaginary and defined as in (5.33). Then, using the properties

(3.55) and (4.32), it can be seen that (6.26) can be expressed as

 

__. _ - _ fl- _ ._
C(l‘.) = e—6t {(K' + 500 C08 6) [~A¥A§J1(A1t) * 31(A5t) - :Ag {J0(A1t)u(t)} * J1(A5f)+
(X)

112— _ _,._ _._ b4 (9 _
—A1J1(A1t)u(t) + 1112011) + 1212951)] + 00.3 sin(A1t)u(t)} +

£00 €oo 1

, 1'6. _ 2A’2

02([1 - 600 COS 6) (ll—201:7?) 6(t), (6.27)

where 3,,(1) is defined in (3.57).

Case 2 occurs when w8+urb2/K2 < 62. In this case, both A1 and A5 are purely real
and the same expression from (6.26) can be used. Case 3 occurs when —/1,-b2 / K 2 <
613 — 62 < 0. In this case, A1 is purely real and A5 is purely imaginary. A5 is defined

as in (5.30) and the expression for C(t) can be written as

b2

C(t) = 6"” {(K + .00 cos6) (13211011) 1.11651) — 6.0.3; {10(A1t)u(t)} * 31651)—

 

2 A _ A _ 4
-I—)—A1I1(A1t)u(t) + A¥I2(A1t) + A§.12(A5t) + b C086 sinh(A1t)u(t) +
Cm 500 1
. — 2K2
b2(K — coo cos 6) (EL—623121?) 6(1‘). (6.28)
m

72

6.2.2 Inversion of the C(s) Term

As described in (6.22), C(s) is given in the Laplace domain as

 

 

 

 

(8 -81)(8 - 82)

(7(3) = . (6.29)

(S - SE)(8 - 8F)(8 - SG)(S - 8H)

Using partial fraction expansion, this term can be rewritten as

C(s) : a — 11,. + 26x cos2 6 1 g 1
4U/\E S — 3E S — 3F
a + 11,. —— 26QC cos2 6 1 1
40AG 3 — 3G 3 — 3H

= C1(S) + c2(3). (6.30)

An inverse Laplace transform of each term in (6.30) can be taken separately. When
inverting 01(3), if (.128 + 72/2 -— X > O the standard transform used in (3.27) can be

applied directly resulting in the expression

 

 

_ . 2.. ,, ‘26
C1(S) 9—) C1(t) : (7 HT :UA;OC C013 e—6t [eAEt _ e—AEt] 11(t)
— .. 2 . :26
z 0 “r + 600“” -5t 3111116539119). (6.31)
QUAE

If 1113 + 72 / 2 — X > 62, A E is purely imaginary and defined as

 

_ ,2
—jAE=AE=\/w3+%—X—62. (6.32)

This allows (6.31) to be rewritten as

 

a — 11.7- + 2600 cos2 6 —<St —
e

c1(t) = _ ’sin(AEt)u(t). (6.33)
20AE

73

If (11(2) + 72/2 -— X < 0 however, BASE} > 0. In this case the identity given in (4.40)

must be used. This produces the result

 

a — 11,- + 260C. cos2 6 —6t
— e
40A};
_ . -2
_U [1.7. + 2600 (10b 66—6te—AEltl,
40AE

c1(t) -—— [eAEtu(-t) + e_’\Etu(t)]

 

(6.34)

which is non-causal. As argued in Chapter 4, this non-causality is due to the infinite
nature of the problem and is not a physical phenomenon.

Taking the Laplace inversion of c2(3) can be done in much the same manner as
for c1(3). If 1.113 + 72 / 2 + X > 0 the standard transform used in (3.27) can be applied

directly resulting in the expression

0 +61: — 2600 cos2 6 —6t
6
40AG

. —2 ,,.>26
: “+62 16°C “0‘ 6"” 51611000119). (6.35)
0 G

 

62(3) 1—» 62(1) = [6*0t — e-AG‘] 11(1)

 

If 0J8 + 72 / 2 + X > 62, AG is purely imaginary and defined as

 

_ 2
—jAG=AG=\/wg+%+x—62. (6.36)

This allows (6.35) to be rewritten as

__ — ‘2 —
a +11r 2cm cos 6 “(St Sill()\0t)u(t)' (6.37)
QUAG

 

C20) =

If 111(2) + 72/2 + X < 0 however, Re{3G} > O. In this case the identity given in (4.40)

must again be used. This produces the result

a + 11-,- — 2600 cos2 6 —6t
— e
40AG

 

62(1) = [Jam—1:) + e—AGtu(t)]

74

_ 2
= _0' + [17‘ 2600 C08 66—6te—AGItl,
4aAG

 

(6.38)

which is again non-causal. The same argument from earlier can be used to explain
this non—causality.
Analyzing these two terms, it can be seen that in general, C (t) is causal whenever

X < —|w3 + 72/2) and non-causal when X > —|w(2) + 72/2l.

6.3 Numerical Results

In order to validate the expressions derived in the previous section, the time-domain
reflection coefficient is evaluated numerically in Fortran, and then compared to the in-
verse FFT of the frequency domain reflection coefficient given in (6.21). The Fortran
code is included in Appendix D. The inverse FFT was done using WaveCalc. The
frequency-domain data was zero-padded up to the maximum limit allowed by Wave-
Calc, 32,768, before the inverse FFT was taken. Since the signal had already decayed
to zero by this point, no windowing was necessary. A set of parameters corresponding
to each of the three possible cases for G (t) is used. The angle is then adjusted for
the case 1 parameters to produce results that contain a non-causal contribution. The
number of points and step size for the results varied between the different cases.
The first set of parameters chosen are: 1110 = 4.0 x 1011 3’1, b = 6.24 x 1011 3’1,
6 = 2.5 x 1010 3’1, 600 = 2, M = 1. This choice of parameters corresponds to case
1. When computing the numerical results, for the frequency-domain data 8,192 fre-
quency points were calculated with a step size of 400 MHz. In the time-domain, 8,192
points were calculated with a step size of 1 x 10‘4 ns. Using 6 = 30°, (6.27) has been
plotted in Figure 6.1 and compared to the inverse FFT. The results show excellent
agreement. Since this function includes only standard Bessel functions, which are
highly oscillatory, and no modified Bessel functions, which are not oscillatory, the

waveform is highly oscillatory and only lightly damped.

The next choice of parameters is: L120 = 4.0 x 1010 3‘1, b = 6.24 x 1010 3‘1,
6 = 2.5 x 1011 3‘1, 600 = 2, 11,. = 1, which corresponds to case 2, (6.26). When
computing the numerical results, for the frequency-domain data 32,768 frequency
points were calculated with a step size of 4 MHz. In the time-domain, 8,192 points
were calculated with a step size of 1 x 10“; ns. The results are shown in Figure
6.2. Again, the closed-form expression and the inverse FFT compare well. For this
choice of parameters 111(2) + urb2 / K 2 < 62, and the resulting waveform is overdamped,
showing no oscillatory behavior and only a single negative peak. Since the expression
for case 2 only involves modified Bessel functions, which do not have the oscillatory
behavior of ordinary Bessel functions, this observed behavior is easily predicted from
the mathematical form of the expression.

The next choice of parameters, which corresponds to case 3, is: 1120 = 4.0 x 109 3’1,
b = 6.24 x1011 3‘1, 6 = 2.5 ><1010 {1,600 = 2,11T = 1, which corresponds to case 3,
(6.28). When computing the numerical results, for the frequency-domain data 32,768
frequency points were calculated with a step size of 400 MHz. In the time—domain,
8,192 points were calculated with a step size of 5 x 10"5 ns. The analytic expression
again matches the inverse FFT, as seen in Figure 6.3. As expected, since 6 > (120, but
62 < (.125 + nrb2/ K 2, there is more damping and less oscillation than with case 1, but
more oscillation than with case 2. Here the expression for the reflection coefficient
has a combination of ordinary and modified Bessel functions.

To examine the non-causal result which is possible for the C (t) term, the material
parameters corresponding to case 1 are used again, but at an angle of 60°. This
causes the 61(t) term in C (t) to be non-causal as shown in (6.34). When computing
the numerical results, for the frequency—domain data 32,768 frequency points were
calculated with a step size of 400 MHz. In the time-domain, 8,192 points were cal-
culated with a step size of 5 x 10’5 us with a starting time value of —1.024 x 10"1

ns. As seen in Figure 6.4, the expression again matches well with the inverse FFT

76

for t > 0. For t < O, non—causal term in the analytic expression is extremely small.
Unfortunately, the inverse FFT does not have the necessary resolution to capture the
small non-causal portion of the reflection coefficient, so it cannot be compared to the
analytic expression.

In order to gain better insight into the non-causal term, the non-causal C (t) term
is plotted in Figure 6.5. The non-causal part is much more significant in this C(t)
than in the TM optical case. To investigate further, C(t) was convolved with C(t)
minus the delta function term. This can be seen in Figure 6.6. Again, the non-causal
portion has a significant contribution to this term. However, when C (t) is convolved
with the delta function term, it can be seen in Figure 6.7 that the non—causal portion
has the opposite sign of the term in Figure 6.6, and when those terms are added, the

non-causal portion becomes extremely small in the final transient reflected field.

77

x 10‘0 Case 1

Reflection Coefficient
O

 

 

 

 

 

 

 

-2
-4
-6 . '
Closed Form Expressuon
~--><- ~ Inverse FFT
_8 1 1 1 m 1
O 0.02 0.04 0.06 0.08 0.1

Time (ns)

Figure 6.1. Time-domain reflection coefficient with incidence angle 6 = 30° and mate-
rial parameters 1.10 = 4.0 x 10113-1, b = 6.24 x 10113-1, 15 = 2.5 x 10103—4, 600 = 2,
11¢ = 1. This choice of parameters corresponds to case 1, Eq. (6.27).

78

x 103 Case 2

Reflection Coefficient

 

 

 

 

 

 

 

 

_7 .1
Closed Form Expression
. - -><- - - Inverse FFT
_8 ‘ 1 1 1 1 1
0 0.1 0.2 0.3 0.4 0.5

Time (ns)

Figure 6.2. Time-domain reflection coefficient with incidence angle 6 = 30° and mate-
rial parameters (.110 = 4.0 x 10103-1, b = 6.24 x 10103—1, 6 = 2.5 x 10113—1, 600 = 2,
11¢ = 1. This choice of parameters corresponds to case 2, Eq. (6.26).

79

x1010 Case 3

Reflection Coefficient

-10

I

 

-12

I

 

Closed Form Expression
- ~ . X- > . Inverse FFT

 

 

 

 

J

0 0.02 0.04 0.06 0.08 0.1
Time (ns)

-14

Figure 6.3. Time—domain reflection coefficient with at an angle of 6 = 30° with param—
eter choices of wo =2 4.0 x 1093-1, b = 6.24 x 10113‘1, 6 = 2.5 x 10103—1, 600 = 2,
11,. = 1. This choice of parameters corresponds to case 3, Eq. (6.28).

80

x 1010 Case 1 - Non-Causal

Reflection Coefficient
O

 

 

 

 

 

 

 

-2 - r
-4 _
-6 .. .
Closed Form Express10n
- - .x. ~ Inverse FFT
_8 1 1 1 L 1 1
-0.02 0 0.02 0.04 0.06 0.08 0.1

Time (ns)

Figure 6.4. Time—domain reflection coefficient with at an angle of 6 = 600 with param-
eter choices of we = 4.0 x 10113-1, b = 6.24 x 10113-1, 6 = 2.5 x 10103—1, 600 = 2,
11, = 1. This choice of parameters corresponds to case 1, Eq. (6.27) with a non-causal
term from Eq. (6.34).

81

Reflection Coefficient

“3 Non-Causal C-term

 

 

 

 

 

 

 

 

 

 

 

 

 

l L l J

 

-0.15 —0.1 —0.05 0 0.05 0.1 0.15 0.2
Time (ns)

Figure 6.5. Non-causal C(t) term for case 1 results.

82

x 1010 Non-Causal C-term convolved with G-term

.0
U1
1

 

Reflection Coefficient
O

 

 

 

-1.5r

 

1 l l l

-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2
Time (ns)

Figure 6.6. Non-causal C (t) term convolved with G (t) minus the delta function term
for case 1 results.

83

x 109 Non-Causal C-term Multiplied by Delta Function Term
3 -

 

 

 

 

 

 

 

 

 

Reflection Coefficient
O

 

 

 

 

 

 

_ 1 1 1 1

3 L l l l
-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2
Time (ns)

Figure 6.7. Non-causal C (t) term convolved with the delta function term for case 1
results.

84

CHAPTER 7

CONCLUSIONS

In this thesis an analytic solution for the transient reflection coefficient of plane waves
from a Lorentz-medium half-space is formulated analytically. Using an inverse Laplace
transform technique, the solution was found to be a combination of Bessel functions
and convolutions of Bessel functions with a decaying exponential factor proportional
to 6, the damping coefficient of the material.

Unlike previous work [10],[11] which represent the transient field as an infinite
sum of fractional—order Bessel functions, the formulation in this thesis has no infinite
sums and presents much more physical insight into the behavior of the transient
field. Depending on field polarization, material properties, and angle of incidence,
the expression for the transient field changes and gives a clearer view of the time-
domain behavior. It is also clear that for certain angles there is a problem with
causality that occurrs in the TM-polarized case. This non-causal function is due to
the infinite nature of the problem and is not a physical phenomenon. It is, however,

an interesting discovery.

7.1 Suggestions for Future Work

The Lorentz model is a single resonance model that is closely related to the permit-
tivity of a cold plasma. A very similar reflection coefficient can be developed using
plasma equations. Using this same method, developing the transient field reflection
coefficient from a cold plasma half-space is a problem that is of some interest.
'IYansmission into dispersive material is also a topic of interest. It would be
interesting to formulate the transient transmitted field using the same basic techniques

outlined in this thesis for the reflected field.

85

APPENDICES

86

1XP19EHVIJIXLIX

FORTRAN CODE FOR OPTICAL TE CASE

program specia1_TE

implicit none

integer nsize

parameter<nsize = 10000)

real*8 BESSI,om0,delta,b,theta,rlam1,rlam5,dt,ct,t,pi
real*8 j1,j2,j1p,j2p,il,i2,df,f,om,gammat(nsize)

rea1*8 j11h(nsize),j12h(nsize),j21h(nsize),j22h(nsize)
rea1*8 111h(nsize),112h(nsize),i21h(nsize),i22h(nsize)
real*8 j1j1(nsize),ilil(nsize),ilj1(nsize),expon(nsize)
integer i,nt,nf

complex*16 zl,z2,s,zj,gamma,zzl,222,sl,s2,35,56

omO = 4.d16

b = dsqrt(20.d32)
delta = 0.28d16

theta = 0.d0

pi = 4.d0*atan(1.d0)

ct = cos(theta*pi/180.d0)

zj=dcmp1x(0.d0,1.d0)

87

df = 8d3

nf = ((4095*2+1)*2+1)
nt = 4096
dt = 1.d-18

221 = cdsqrt(dcmplx(delta*delta-om0*om0,0.dO))

81 = -delta+221

s2 = -delta-221

222 = cdsqrt(dcmplx(delta*delta-omO*omO-b*b/(ct*ct),0.d0))
$5 = -de1ta+222

$6 = -de1ta-222

open (10,file=’gamma.dat’,status=’unknown’)
do i=1,nf

f = df*i

om = 2.d0*pi*f*1.d9

s = zj*om

21 cdsqrt(s-sl)*cdsqrt(s-s2)

22

cdsqrt(s-s5)*cdsqrt(s-SS)

gamma = (21-22)/(21+z2)

write (10,*) f,rea1(gamma),aimag(gamma)
end do

close (10)

do i=1,nt
t = (i-1)*dt+.00001*dt

expon(i) = exp(-delta*t)

88

end do

write(*,*) ’om‘2 = ’,om0*om0
write(*,*) ’om‘2+B‘2 = ’,om0*om0+b*b/(ct*ct)
write(*,*) ’delta“2 = ’,delta*delta

CASE 1

om0“2 > delta“2

 

if((om0*om0).gt.(delta*de1ta)) then
write(*,*) "CASE 1"
rlaml = dsqrt(omO*omO-delta*delta)
r1am5 = dsqrt(omO*omO+(b*b/(ct*ct))-de1ta*de1ta)
do i=1,nt
t = (i-1)*dt+.00001*dt
call jbess(rlam1*t,2,j2,j2p)
j21h(i) = j2/(rlam1*t)
call jbessCrlam5*t,2,j2,j2p)
j22h(i) = j2/(rlam5*t)
call jbess(r1am1*t,1,j1,j1p)
j11h(i) = j1/(r1am1*t)
call jbess(rlam5*t,1,j1,j1p)
j12h(i) = jl/(rlam5*t)
end do

call convfjllh,j12h,nt,dt,j1j1)

open(10,file=’gammat.dat’,status=’unknown’)

do i=1,nt

89

gammat(i) = expon(i)*(2.d0*(rlam1**3)*j21h(i)+
2 2.d0*(r1am5**3)*j22h(i)-2.d0*r1am1*rlam1*r1am5*rlam5*
3 3131(1))
gammat(i) = -gammat(i)*ct*ct/(b*b)
write(10,*) ((i-1)*dt+.00001*dt)*1.d9,gammat(i)
end do
close(10)
CASE 2
om0‘2 + B“2 < delta“2

elseif((omO*omO+(b*b/(ct*ct))).lt.de1ta*de1ta) then

write(*,*) "CASE 2"
rlaml = dsqrt(delta*de1ta-om0*om0)
r1am5 = dsqrt(-om0*omO-(b*b/(ct*ct))+de1ta*delta)
do i=1,nt

t = (i-1)*dt+.00001*dt

i2 = BESSI(2,rlam1*t)

i21h(i) = i2/(rlam1*t)

i2 = BESSI(2,rlam5*t)

i22h(i) = i2/(rlam5*t)

ii = BESSI(1,rlam1*t)

111h(i) = i1/(r1am1*t)

i1 = BESSI(1,r1am5*t)

112h(i) = il/(rlam5*t)

end do

90

call conv(i11h,il2h,nt,dt,ilil)

open(10,file=’gammat.dat’,status=’unknown’)
do i=1,nt
gammat(i) = expon(i)*(2.d0*(rlam1**3)*i21h(i)+
2 2.dO*(r1am5**3)*i22h(i)-2.dO*rlam1*rlam1*rlam5*rlam5*
3 1111(1))
gammat(i) = -gammat(i)*ct*ct/(b*b)
write(10,*) ((i-l)*dt+.00001*dt)*1.d9,gammat(i)
end do
close(10)
else
CASE 3
om0“2 + 8‘2 > delta‘2 & om0“2 < delta‘2
write(*,*) "CASE 3"
rlaml = dsqrt(delta*de1ta-om0*om0)
rlamS = dsqrt(omO*omO+(b*b/(ct*ct))-delta*delta)
do i=1,nt
t = (i-1)*dt+.00001*dt
12 = BESSI(2,r1am1*t)
i21h(i) = i2/(rlam1*t)
call jbess(r1am5*t,2,j2,j2p)
j22h(i) = j2/(r1am5*t)
11 = BESSI(1,r1am1*t)

ilthi) = il/(rlam1*t)

91

call jbess(r1am5*t,1,j1,j1p)
j12h(i) = jl/(rlam5*t)
end do

call conv(111h,j12h,nt,dt,iljl)

open(10,file=’gammat.dat’,status=’unknown’)
do i=1,nt
gammat(i) = expon(i)*(2.d0*(rlam1**3)*i21h(i)+
2 2.dO*(rlam5**3)*j22h(i)+2.d0*r1am1*r1am1*rlam5*r1am5*
3 ilj1(i))
gammat(i) = -gammat(i)*ct*ct/(b*b)
write(10,*) ((i-l)*dt+.00001*dt)*1.d9,gammat(i)
end do
close(10)

endif

end

subroutine conv(f,g,n,de1,c)

performs convolution of waveforms f and g using

linear interpolation

parameter (nsi2e=10000)
parameter (C13=0.33333333333333333333d0,

2 c16=0.16666666666666666666d0)

92

implicit rea1*8 (a-h,o-2)
real*8 f(nsize),gfnsize),c(nsize)
c(1)=0.d0
do k=2,n
sum = O.d0
do l=1,k-1

sum = sum + c13*f(1+1)*g(k-1)

2 + c13*f(l)*g(k-l+1)

3 + c16*f(l)*g(k-l)

4 + c16*f(1+1)*g(k-1+1)
end do

C(k) = del*sum
end do
return

end

FUNCTION BESSI(N,X)

This subroutine calculates the first kind modified Bessel function
of integer order N, for any REAL X. We use here the classical
recursion formula, when X > N. For X < N, the Miller’s algorithm
is used to avoid overflows.

REFERENCE:

C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS,

MATHEMATICAL TABLES, VOL.5, 1962.

93

PARAMETER (IACC = 40,BIGNO = 1.010, BIGNI = 1.D-10)
REAL *8 X,BESSI,BESSIO,BESSI1,TOX,BIM,BI,BIP
IF (N.EQ.O) THEN

BESSI = BESSIO(X)

RETURN

ENDIF

IF (N.EQ.1) THEN

BESSI = BESSIlCX)

RETURN

ENDIF

IF(X.EQ.O.DO) THEN

BESSI=O.DO

RETURN

ENDIF

TOX = 2.D0/X

BIP

O.DO

BI = 1.DO

BESSI = O.DO

M = 2*((N+INT(SQRT(FLOAT(IACC*N)))))

D0 12 J = M,1,-1

BIM = BIP+DFLDAT(J)*TOX*BI
BIP = BI
BI = BIM

IF (ABS(BI).GT.BIGNO) THEN

BI BI*BIGNI

BIP BIP*BIGNI

94

BESSI = BESSI*BIGNI
ENDIF
IF (J.EQ.N) BESSI = BIP
12 CONTINUE
BESSI = BESSI*BESSIO(X)/BI
RETURN

END

 

! Auxiliary Bessel functions for N=O, N=1
FUNCTION BESSIO(X)
REAL *8 X,BESSIO,Y,P1,P2,P3,P4,P5,P6,P7,
* Qi,Q2,Qs,Q4,Qs,Qs,Q7,Qa,Qe,Ax,Bx
DATA P1,P2,P3,P4,P5,P6,P7/1 D0,3.5156229DO,3.0899424D0,1.2067429D0
* ,0.2659732D0,0 360768D-1,0.45813D-2/
DATA Ql,Q2,QB,Q4,QS,06,Q7,Q8,09/0.39894228D0,0.1328592D-1,
* 0.225319D-2,-0.157565D-2,0.916281D-2,-0.2057706D-1,
* O.2635537D-1,-0.1647633D-1,0.392377D-2/
IF(ABS(X).LT.3 7500) THEN
Y=(X/3.75D0)**2
BESSIO=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
ELSE
Ax=ABs(x)
Y=3.75D0/AX
BX=EXP(AX)/SQRT(AX)
AX=Ql+Y*(Q2+Y*(QB+Y*(Q4+Y*(QS+Y*(QS+Y*(Q?+Y*(08+Y*09)))))))
BESSIO=AX*BX

ENDIF

95

RETURN

END

 

FUNCTION BESSI1(X)

REAL *8 X,BESSIl,Y,P1,P2,P3,P4,P5,P6,P7,

* Ql,Q2,03,Q4,QS,QS,Q7,Q8,09,AX,BX

DATA P1,P2,P3,P4,P5,P6,P7/O.5D0,0.87890594D0,0.51498869D0,
* 0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/

DATA Q1,Q2,QB,Q4,QS,QG,Q7.08,09/O.39894228DO,-0.3988024D-1,
* -o 362018D-2,0.163801D-2,-0 1031555D-1,0.2282967D-1,
* -o.2895312D-1,o.1787654D-1,-0.420059D-2/
IF(ABS(X).LT.3.75DO) THEN

Y=(X/3.75D0)**2
BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))

ELSE

Ax=ABS(x)

Y=3.75D0/Ax

BX=EXP<Ax)/SQRT(AX)
AX=Ql+Y*(Q2+Y*(QB+Y*(Q4+Y*(QS+Y*(06+Y*(Q7+Y*(08+Y*09)))))))
BESSIl=AX*BX

ENDIF

RETURN

END

subroutine jbess (x,n,bj,bjp)

96

c calculates jn and jn’ for n positive or negative

c
implicit real*8 (a-h,o-2)
c
if (n .eq. 0) then
bj = bessj(0,x)
bjp = -bessj(1,x)
return
endif
c
if (n .ge. 0) then
m = n
else
m = -n
endif
c

if (x .eq. 0.d0) then

bj = 0.d0

bjpl = bessj(m+1,x)

bjml = besstm-1,x)

bjp = (bjml-bjp1)/2.d0
else

bj = bessj(m,x)

bj1 = bessj(m+1,x)

bjp = -bj1 + n*bj/x
endif

97

if (n .1t. 0) then

bj = bj*(-1)**n
bjp = bjp*<-1)**n
endif
return
end

FUNCTION BESSJO(X)

implicit rea1*8 (a-h,o-z)

REAL*8 Y,P1,P2,P3,P4,P5,Ql,QZ,Q3,Q4,QS,R1,R2,R3,R4,R5,R6,
* Sl,S2,S3,S4,SS,SG

DATA P1,P2,P3,P4,P5/1.D0,-.1098628627D-2,.2734510407D-4,

* -.2073370639D-5,.2093887211D-6/, 01,02,QB,Q4,Q5/-.1562499995D-
*1,
* .14304887650-3,-.6911147651D-5,.7621095161D-6,-.934945152D-7/

DATA R1,R2,R3,R4,R5,R6/57568490574 DO,-13362590354.D0,651619640.7D
*0,
* -11214424.18DO,77392.33017D0,-184.9052456D0/,
* Sl,S2,SS,S4,SS,SS/57568490411.D0,1029532985.D0,
* 9494680.718D0,59272.64853DO,267.8532712D0,1.D0/
IF(ABS(X).LT.8.)THEN
Y=X**2

BESSJO=(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))

98

* /(Sl+Y*(S2+Y*(83+Y*(S4+Y*(SS+Y*86)))))
ELSE
Ax=ABS(x)
z=8./Ax
Y=Z**2
xx=Ax-.785398164
BESSJO=SQRT(.636619772/AX)*(COS(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y
* *P5))))-Z*SIN(XX)*(QI+Y*(02+Y*(QB+Y*(Q4+Y*QS)))))
ENDIF
RETURN

END

FUNCTION BESSJ1CX)

implicit real*8 (a-h,o-z)

REAL*8 Y,F1,P2,P3,P4,P5,Q1,Q2,QS,Q4,Q5,R1,R2,R3,R4,R5,R6,
* Sl,SQ,S3,S4,SS,86

DATA R1,R2,R3,R4,R5,R6/72362614232.DO,-7895059235.D0,242396853.1D0

1k

9

* -2972611.439DO,15704.48260D0,-30.1603660600/,
* Sl,S2,SB,S4,SS,36/144725228442.D0,2300535178.D0,
* 18583304.74D0,99447.43394D0,376.9991397D0,1.DO/

DATA P1,P2,P3,P4,P5/1.D0,.183105D-2,-.3516396496D-4,.2457520174D-5

:1:

2

* -.240337019D-6/, Ql,QZ,Q3,Q4,QS/.O4687499995DO,-.2002690873D-3

*

99

* .8449199096D-5,-.88228987D-6,.105787412D-6/
IF(ABS(X).LT.8.)THEN
Y=X**2
BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))
* /(SI+Y*(82+Y*(S3+Y*(S4+Y*(SS+Y*86)))))
ELSE
AX=ABS(X)
Z=8./AX
Y=Z**2
XX=AXf2.356194491
BESSJ1=SQRT(.636619772/AX)*(COS(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y
* *P5))))-Z*SIN(XX)*(01+Y*(02+Y*(QB+Y*(Q4+Y*QS)))))
* *SIGN(1.,X)
ENDIF
RETURN

END

FUNCTION BESSJ(N,X)
implicit real*8 (a-h,o-2)
PARAMETER (IACC=40,BIGNO=1.le,BIGNI=1.d-lO)

if (x .ge. O.d0) then

ff = 1.dO
else

ff = (-1.d0)**n
end if

100

11

xx = abs(x)
if (n .eq. 0) then
bessj = ff*bessj0(xx)
return
endif
if (n .eq. 1) then
bessj = ff*bessj1(xx)
return
endif
TOX=2./Xx
IF(Xx.GT.FLOAT(N))THEN

BJM=BESSJO(XX)

BJ=BESSJ1<XX)

DO 11 J=1,N-1
BJP=J*TOX*BJ-BJM
BJM=BJ
BJ=BJP

CONTINUE

BESSJ=ff*BJ

ELSE

M=2*((N+INT(SQRT(FLOAT(IACC*N))))/2)

BESSJ=O.

JSUM=0

SUM=O.

BJP=O.

BJ=1.

DO 12 J=M,1,-1

101

BJM=J*TOX*BJ-BJP
BJP=BJ
BJ=BJM
IF(ABS(BJ).GT.BIGNO)THEN
BJ=BJ*BIGNI A
BJP=BJP*BIGNI
BESSJ=BESSJ*BIGNI
SUM=SUM*BIGNI
ENDIF
IF(JSUM.NE.O)SUM=SUM+BJ
JSUM=1-JSUM
IF(J.EQ.N)BESSJ=BJP
12 CONTINUE
SUM=2.*SUM-BJ
BESSJ=ff*BESSJ/SUM
ENDIF
RETURN

END

102

1AFHPEEJEHD(]3

FYDIEFILAJV CXDI)EJFYDIL(DIVFI(1AJQTERA CHASHE

program special_TM

implicit none

integer nsize

parameter(nsize = 10000)

real*8

rea1*8

real*8

real*8

real*8

rea1*8

real*8

real*8

BESSI,om0,de1ta,b,theta,dt,ct,tt,t,pi,df,f,om,YY
rlaml,rlam5,r1ame,expon(nsize),gams,temp,term1,term2
101(nsi2e),ill(nsize),115(nsize),i21(nsize),i25(nsize)
illh(nsize),115h(nsi2e),i21h(nsize),i25h(nsize),iOi2(nsi2e)
i1i1(nsi2e),i2i2(nsize),i011(nsize),iOiO(nsize),iOS(nsize)
c3ft(nsize),c3f_t(nsize),f_t(nsize),temp2(nsize)
c3(nsize),ft(nsize),gammat(nsize)

c33hift(nsize),ftshiftCnsize),f_tshift(nsize)

integer i,nt,nf,c3case3,shift

complex*16 21,22,s,zj,gamma,zz1,222,sl,s2,sS,36,zze,se,sf,num

omO

b

delta

theta

4.d16
dsqrt(20.d32)
0.28d16

30.d0

pi = 4.d0*atan(1.d0)

103

ct cos(theta*pi/180.d0)

tt tan(theta*pi/180.d0)
2j=dcmplx(0.d0,1.d0)
gams = b*b/(tt*tt-1)

c3case3 = 0

shift = 4

df = 8d3

nf = 4096*2*2*2
nt = 4096

dt = 1.d-18

221 = cdsqrt(dcmplx(delta*delta-om0*om0,0.d0))
31 = -delta+221
s2 = -delta-zz1
222 = cdsqrt(dcmplx(delta*delta-omO*omO-b*b/(ct*ct),O.d0))
s5 = -delta+222
$6 = -delta-222
zze = cdsqrt(dcmplx(delta*delta-om0*om0+gams,O.d0))
se = -de1ta+22e

sf = -delta-zze

open (10,file=’gamma.dat’,status=’unknown’)
do i=1,nf
f = df*i

om = 2.d0*pi*f*1.d9

104

s = 2j*om

21

cdsqrt(s-sl)*cdsqrt(s-s2)*cdsqrt(s-85)*cdsqrt(s-sG)

22 (s-sl)*(s-s2)+b*b

gamma = (21-22)/(21+z2)

write (10,*) f,rea1(gamma),aimag(gamma)
end do
close (10)
do i=1,nt

t = (i-1)*dt+.00001*dt

expon(i) = exp(-delta*t)

end do

if((om0*om0).gt.(delta*delta+gams)) then
write(*,*) "C term 1"
rlame = dsqrt(omO*om0-delta*de1ta-gams)
do i=1,nt
t = (i-1)*dt+.00001*dt
c3(i) = expon(i)*sin(rlame*t)/rlame
end do
elseif((om0*om0).gt.(gams)) then
write(*,*) "C term 2"
rlame = dsqrt(delta*de1ta-om0*om0+gams)
do i=1,nt

t = (i-1)*dt+.00001*dt

105

c3(i) = expon(i)*sinh(rlame*t)/r1ame
end do
else
write(*,*) "C term 3"
rlame = dsqrt(delta*delta-om0*om0+gams)
c3case3 = 1
rlame = dsqrt(delta*delta-om0*om0+gams)
do i=1,nt/shift
t = (i-l-nt/Shift)*dt+.00001*dt
c3(i) = -exp((-delta+r1ame)*t)/(2.d0*r1ame)
end do
do i=nt/shift+1,nt
t = (i-l-nt/Shift)*dt+.00001*dt
c3(i) = -exp(-(de1ta+r1ame)*t)/(2.d0*rlame)
end do

endif

 

 

Case 1: om0‘2>delta“2

 

 

if((om0*om0).gt.(delta*delta)) then
write(*,*) "CASE 1"
rlaml = dsqrt(omO*om0-de1ta*delta)
rlam5 = dsqrt(om0*om0+(b*b/(ct*ct))-de1ta*delta)
do i=1,nt

t = (i-1)*dt+.00001*dt

106

call jbess(rlam1*t,0,101(i),temp)
call jbess(rlam1*t,1,ill(i),temp)
call jbess(rlam5*t,1,115(i),temp)
call jbess(r1am1*t,2,i21(i),temp)

call jbess(rlam5*t,2,i25(i),temp)

i11h(i) = 111(i)/t

i15h(i) = 115(i)/t

121h(i) = i21(i)/t

i25h(i) = i25(i)/t
end do

call convCillh,115h,nt,dt,ilil)
call conv(i01,i15h,nt,dt,iOi1)
do i=1,nt
t = (i-1)*dt+.00001*dt
ft(i) = 2.d0*expon(i)*(-rlam1*r1am5*ili1(i)+r1am1*rlam1*
2 i21h(i)+r1am5*r1am5*i25h(i)-b*b*rlam5*iOil(i)+

3 b*b*r1am1*il1(i)+b*b*b*b*sin(r1am1*t)/(2.d0*r1am1))

 

elseif((om0*om0).1t.(delta*delta-(b*b/(ct*ct)))) then
write(*,*) "CASE 2"
rlaml = dsqrt(delta*de1ta-om0*om0)
r1am5 = dsqrt(-om0*omO-(b*b/(ct*ct))+de1ta*delta)
do i=1,nt

t = (i-1)*dt+.00001*dt

107

101(1) = BESSI(O,rlam1*t)
i11(i) = BESSI(1,r1am1*t)
115(1) = BESSI(1,r1am5*t)
i21(i) = BESSI(2,r1am1*t)
125(1) = BESSI(2,r1am5*t)
111h(i) = 111(i)/t
115h(i) = i15(i)/t
i21h(i) = i21(i)/t
i25h(i) = i25(i)/t
end do

call conv(111h,115h,nt,dt,1111)
call conv(i01,115h,nt,dt,iOil)
do i=1,nt
t = (i-1)*dt+.00001*dt
ft(i) = 2.d0*expon(i)*(Erlam1*rlam5*ilil(i)+rlam1*rlam1*
2 i21h(i)+r1am5*r1am5*125h(i)+b*b*r1am5*iOil(i)-

3 b*b*r1am1*ill(i)+b*b*b*b*sinh(r1am1*t)/(2.d0*r1am1))

__—-_-—--—-————u.--—_—————.——a—c—u—o—_-—--—~—___————————————-——-————————-—-—-——

 

else
write(*,*) "CASE 3"
r1am1 = dsqrt(delta*delta-om0*om0)
rlam5 = dsqrt(omO*om0+(b*b/(ct*ct))-de1ta*delta)
do i=1,nt

t = (i-1)*dt+.00001*dt

108

101(1)

BESSI(O,r1am1*t)

111(1) BESSI(1,rlam1*t)

call jbess(r1am5*t,1,i15(i),temp)
i21(i) = BESSI(2,r1am1*t)

call jbess(rlam5*t,2,i25(i),temp)

111h(i) = il1(i)/t

i15h(i) = 115(i)/t

i21h(i) = i21(i)/t

i25h(i) = i25(i)/t
end do

call conv(illh,115h,nt,dt,ili1)
call conv(i01,ilSh,nt,dt,iOil)
do i=1,nt

t = (i-1)*dt+.00001*dt

ft(i) = 2.d0*expon(i)*(r1am1*rlam5*ilil(i)+r1am1*rlam1*

2 i21h(i)+r1am5*r1am5*i25h(i)-b*b*rlam5*iOil(i)-

3 b*b*r1am1*ill(i)+b*b*b*b*sinh(r1am1*t)/(2.d0*r1am1))
end do

endif

 

term2 = b*b*(tt*tt-1)

if (theta.eq.45.d0) then
call conv(ft,ft,nt,dt,f_t)
do i=1,nt

gammat(i) = -f_t(i)/(b*b*b*b)

109

end do
elseif (c3case3.eq.0) then
call conv(c3,ft,nt,dt,c3ft)
do i=1,nt
gammat(i) = (ft(i)+gams*c3ft(i))/term2
temp2(i) = gammat(i)
end do
else
do i=1,(nt/shift)
ftshift(i) = 0
end do
do i=(nt/shift+1),nt
ftshift(i) = ft(i-nt/shift)
end do
call conv(c3,ftshift,nt,dt,c3ft)
do i=1,(nt/shift)
f_tshift(i) = 0
end do
do i=(nt/shift+1),nt
f_tshift(i) = ftshiftCi-nt/shift)
end do
do i=1,nt
gammat(i) = (f_tshift(i)+gams*c3ft(i))/term2
end do
endif
open(10,file=’gammat.dat’,status=’unknown’)

do i=1,nt

110

if (c3case3.eq.0) then
write(10,*) ((i-l)*dt+.00001*dt)*1.d9,gammat(i)
else
write(10,*) ((1-1-2*nt/(Shift))*dt+.00001*dt)*1.d9,
2 gammat(i)
endif
end do
close(10)

end

subroutine conv(f,g,n,de1,c)

performs convolution of waveforms f and g using

linear interpolation

parameter (nsize=10000)
parameter (C13=0.33333333333333333333d0,
2 C16=O.16666666666666666666d0)
implicit rea1*8 (a-h,o-z)
rea1*8 f(nsize),g(nsize),c(nsize)
c(1)=0.d0
do k=2,n
sum = 0.d0
do 1=1,k-1
sum = sum + c13*f(1+1)*g(k—1)

2 + c13*f(l)*g(k-1+1)

111

3 + c16*f(1)*g(k-l)
4 + c16*f(l+1)*g(k-l+1)
end do
C(k) = del*sum
end do
return

end

FUNCTION BESSI(N,X)

This subroutine calculates the first kind modified Bessel function
of integer order N, for any REAL x. We use here the classical
recursion formula, when X > N. For X < N, the Miller’s algorithm
is used to avoid overflows.

REFERENCE:

C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS,

MATHEMATICAL TABLES, VDL.5, 1962.

PARAMETER (IACC = 40,BIGNO = 1.D10, BIGNI = 1.D-10)
REAL *8 X,BESSI,BESSIO,BESSI1,TOX,BIM,BI,BIP

IF (N.EQ.0) THEN

BESSI = BESSIOCX)

RETURN

ENDIF

IF (N.EQ.1) THEN

112

BESSI = BESSIlCX)
RETURN

ENDIF
IF(X.EQ.O.DO) THEN
BESSI=O.DO

RETURN

ENDIF

TOX = 2.DO/X

BIP

II
0

.00

BI = 1.D0

BESSI 0.00
M = 2*((n+INT(SQRT(FLOAT(IACC*N)))))

D0 12 J = M,1,-1

BIM = BIP+DFLDAT(J)*TDX*BI
BIP = BI
BI = BIM

IF (ABS(BI).GT.BIGNO) THEN

BI

BI*BIGNI

BIP

BIP*BIGNI
BESSI = BESSI*BIGNI
ENDIF
IF (J.EQ.N) BESSI = BIP
12 CONTINUE
BESSI = BESSI*BESSIO(X)/BI
RETURN

END

! Auxiliary Bessel functions for N=0, N=1
FUNCTION BESSIO(X)
REAL *8 X,BESSIO,Y,P1,P2,P3,P4,P5,P6,P7,
* Q1,02,Q3,Q4,QS,QG,Q7,Q8,QQ,AX,BX
DATA P1,P2,P3,P4,P5,P6,P7/1.DO,3.515622QDO,3.0899424D0,1.2067429D0
* ,0 2659732D0,0 360768D-I,O 45813D-2/
DATA Ql,02,Q3,Q4,QS,QS,Q7,QB,QQ/O.39894228D0,0.1328592D-1,
* 0.225319D-2,-O.1575650-2,0.916281D-2,-0.20577060-1,
* O.2635537D-1,-O.1647633D-1,0.392377D-2/
IF(ABS(X).LT.3.75DO) THEN
Y=(X/3.75D0)**2
BESSIO=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
ELSE
AX=ABS(X)
Y=3.75DO/AX
BX=EXPCAX)/SQRT(AX)
AX=Ql+Y*(Q2+Y*(Q3+Y*(Q4+Y*(QS+Y*(Q6+Y*(Q7+Y*(Q8+Y*QQ)))))))
BESSIO=AX*BX
ENDIF
RETURN

END

FUNCTION BESSIiCX)

REAL *8 X,BESSIl,Y,P1,P2,P3,P4,P5,P6,P7,

* Ql,Q2,Q3,Q4,QS,QG,Q7,Q8,Q9,AX,BX

DATA P1,P2,P3,P4,P5,P6,P7/O.5D0,0.87890594D0,0.51498869D0,

* O.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/

114

DATA Q1,Q2,QS,Q4,QS,QS,Q7,Q8,QQ/O.39894228DO,-0.3988024D-1,
* -O.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1,

* -O.2895312D-1,0.1787654D-1,-O.420059D-2/
IF(ABS(X).LT.3.75DO) THEN

Y=(X/3.75DO)**2
BESSII=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))

ELSE

AX=ABS(X)

Y=3.75DO/AX

BX=EXP<AX)/SQRT(AX)
AX=QI+Y*(Q2+Y*(QB+Y*(Q4+Y*(Q5+Y*(QG+Y*(Q7+Y*(Q8+Y*QQ)))))))
BESSI1=AX*BX

ENDIF

RETURN

END

subroutine jbess (x,n,bj,bjp)

calculates jn and jn’ for n positive or negative

implicit rea1*8 (a-h,o-z)

if (n .eq. 0) then

bj bessj(0,x)

bjp -bessj(1,x)

115

return

endif

if (n .ge. 0) then

endif

if (x .eq. 0 d0) then

bj = 0.d0

bjpl = bessj(m+1,x)

bjml = bessj(m-1,x)

bjp = (bjml-bjp1)/2.d0
else

bj = bessj(m,x)

bj1 = bessj(m+1,x)

bjp = -bj1 + n*bj/x
endif

if (n .1t. 0) then

bj = bj*(-1)**n
bjp = bjp*<-1>**n
endif
return
end

116

FUNCTION BESSJO(X)

implicit rea1*8 (a-h,o-z)

REAL*8 Y,P1,P2,P3,P4,P5,Ql,Q2,03,Q4,QS,R1,R2,R3,R4,R5,R6,
* 31,82,33,S4,S5,36

DATA P1,P2,P3,P4,P5/1.D0,-.10986286270-2,.2734510407D-4,

* -.20733706390-5,.2093887211D-6/, Ql,QZ,QB,Q4,QS/-.1562499995D-
*1,
* .14304887650-3,- 6911147651D—5, 7621095161D-6,- 934945152D-7/

DATA R1,R2,R3,R4,R5,R6/57568490574 DO,-13362590354.DO,651619640 7D

*0,

* -11214424.18D0,77392.33017DO,-184.9052456D0/,

* Sl,82,S3,S4,SS,SG/57568490411.DO,1029532985.D0,

* 9494680.718DO,59272.64853DO,267.8532712D0,1.DO/

IF(ABS(X).LT.8.)THEN

Y=X**2
BESSJO=<R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))

* /(S1+Y*(82+Y*(SB+Y*(S4+Y*(SS+Y*36)))))

ELSE
AX=ABS(X)
Z=8./AX
Y=Z**2
XX=AX-.785398164
BESSJO=SQRT(.636619772/AX)*(COS(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y

* *P5))))-Z*SIN(XX)*(01+Y*(Q2+Y*(QB+Y*(Q4+Y*05)))))

117

ENDIF
RETURN

END

FUNCTION BESSJICX)

implicit real*8 (a-h,o-z)

REAL*8 Y,P1,P2,P3,P4,P5,Ql,Q2,QB,Q4,QS,R1,R2,R3,R4,R5,R6,
* Sl,SZ,S3,S4,SS,36

DATA R1,R2,R3,R4,R5,R6/72362614232.DO,-7895059235.D0,242396853.1DO

*

’

* -2972611.439DO,15704.48260D0,-30.16036606D0/,
* Sl,32,S3,S4,SS,SG/144725228442.DO,2300535178.D0,
* 18583304.74DO,99447.43394D0,376.9991397D0,1.DO/

DATA P1,P2,P3,P4,P5/1.DO,.183105D-2,-.35163964960-4,.2457520174D-5

*

’

* -.240337019D-6/, Q1,QQ,QB,Q4,QS/.04687499995D0,-.2002690873D-3

*

* .8449199096D-5,-.88228987D-6,.105787412D-6/
IFCABS(X).LT.8.)THEN
Y=X**2

BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))

* /(Si+Y*(32+Y*(SB+Y*(S4+Y*(SS+Y*SG)))))
ELSE
AX=ABS(X)
Z=8./AX
118

Y=Z**2
XX=AX-2.356194491
BESSJ1=SQRT(.636619772/AX)*(COS(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y
* *P5))))-Z*SIN(XX)*(Ql+Y*(Q2+Y*(QB+Y*(Q4+Y*QS)))))
* *SIGN(1.,X)
ENDIF
RETURN

END

FUNCTION BESSJ(N,X)
implicit real*8 (a-h,o-z)
PARAMETER (IACC=40,BIGNO=1.d10,BIGNI=1.d-10)

if (x .ge. 0 d0) then

ff = 1.d0
else

ff = (-1.d0)**n
end if

xx = abs(x)
if (n .eq. 0) then
bessj = ff*bessj0(xx)
return
endif
if (n .eq. 1) then
bessj = ff*bessj1(xx)

return

119

endif
TOX=2./Xx
IF(Xx.GT.FLOAT(N))THEN

BJM=BESSJO<XX)

BJ=BESSJ1(XX)

DO 11 J=1,N-1
BJP=J*TOX*BJ-BJM
BJM=BJ
BJ=BJP

11 CONTINUE

BESSJ=ff*BJ

ELSE

M=2*((N+INT(SQRT(FLOAT(IACC*N))))/2)

BESSJ=O.

JSUM=O

SUM=O.

BJP=O.

BJ=1.

DO 12 J=M,1,-1
BJM=J*TOX*BJ-BJP
BJP=BJ
BJ=BJM
IF(ABS(BJ).GT.BIGNO)THEN

BJ=BJ*BIGNI
BJP=BJP*BIGNI
BESSJ=BESSJ*BIGNI

SUM=SUM*BIGNI

120

ENDIF
IF(JSUM.NE.O)SUM=SUM+BJ
JSUM=1-JSUM
IF(J.EQ.N)BESSJ=BJP
12 CONTINUE
SUM=2.*SUM-BJ
BESSJ=ff*BESSJ/SUM
ENDIF
RETURN

END

121

APPENDIX C

FYDIKTILAJV CXDI)E3FT)FL(3EHVEHRJ\L.TUE(CJKSIE

program genera1_TE

implicit none

integer nsize

parameter(nsize = 50000)

real*8 BESSI,om0,de1ta,b,theta,dt,ct,tt,t,pi,df,f,om,YY
rea1*8 rmur,epsinf,st,rkap,f1,f2

rea1*8 rlaml,rlam5,rlama,expon(nsize),gams,temp

real*8 i11h(nsize),115h(nsize),i15(nsize),i11(nsize)
real*8 i21h(nsize),i25h(nsize),i1i1(nsize)

rea1*8 c3ft(nsize),c35hift(nsize),ftshift(nsize)

real*8 c3(nsize),ft(nsize),gammat(nsize)

integer i,nt,nf,c3case3,shift

complex*16 zI,z2,s,zj,gamma,zzl,zz2,sl,s2,35,86

rmur = 1.d0

epsinf = 2.d0

om0 4.d11

b 6.24d11

122

delta

theta

pi

ct

St

1:13

4

C

S

t

2.5d10

60.d0

.d0*atan(1.d0)
os(theta*pi/180.d0)
in(theta*pi/180.d0)

an(theta*pi/180.d0)

zj=dcmplx(0.d0,1.d0)

c3case3 = 0

shift
rkap
gams

f1 =

f2

df

nf =

nt=

dt

2

f

0

4

4

2

=8
dsqrt(rmur*epsinf—st*st)
om0*om0-(rmur*b*b)/(rmur*rmur*ct*ct-rkap*rkap)
.d0*rkap*rmur*ct/(rmur*ct+rkap)
1/(rmur*rmur*ct*ct-rkap*rkap)
.4d0
095*2
096*2

.d-13

write(*,*) omO*omO

write(*,*) delta*delta-rmur*b*b/(rkap*rkap)

221 =

81 =

s2 =

222 =

SS =

$6 =

cdsqrt(dcmplx(delta*delta-om0*om0,0.dO))

-delta+zzl

-de1ta-zzl
cdsqrt(dcmplx(delta*delta-om0*omO-rmur*b*b/(rkap*rkap),O.d0))
-de1ta+zz2

-de1ta-zz2

123

open (10,file=’gamma.dat’,status=’unknown’)
do i=1,nf

f = df*i

om = 2.d0*pi*f*1.d9

s = zj*om

21 rmur*ct*cdsqrt(s-sl)*cdsqrt(s-s2)

22 rkap*cdsqrt(s-s5)*cdsqrt(s-36)
gamma = (zl-z2)/(21+z2)-(rmur*ct-rkap)/(rmur*ct+rkap)
write (10,*) f,real(gamma),aimag(gamma)

end do

close (10)

do i=1,nt
t = (i-l)*dt+.00001*dt
expon(i) = exp(-de1ta*t)

end do

if(gams.gt.(delta*delta)) then
write(*,*) "C term 1"
rlama = dsqrt(gams-delta*delta)

do i=1,nt

124

t = (i-1)*dt+.00001*dt
c3(i) = expon(i)*sin(rlama*t)/rlama
end do
elseif(0.gt.gams) then
write(*,*) "C term 3"
rlama = dsqrt(delta*delta-gams)
c3case3 = 1
do i=1,nt/shift
t = (i-l-nt/Shift)*dt+.00001*dt
c3(i) = -exp((-de1ta+rlama)*t)/(2.d0*rlama)
end do
do i=nt/shift+1,nt
t = (i-1-nt/shift)*dt+.00001*dt
c3(i) = ~exp(-(delta+rlama)*t)/(2.d0*r1ama)

end do

else
write(*,*) "C term 2"
rlama = dsqrt(delta*delta-gams)
do i=1,nt
t = (i-1)*dt+.00001*dt
c3(i) = expon(i)*sinh(rlama*t)/rlama

end do

Calculate G Term

125

Case 1: om0“2>delta“2

 

if((om0*om0).gt.(delta*de1ta)) then

write(*,*) "CASE 1"

r1am1 = dsqrt(omO*omO-delta*delta)
rlam5 = dsqrt(om0*om0+(rmur*b*b/(rkap*rkap))-delta*de1ta)
do i=1,nt

t = (i-1)*dt+.00001*dt

call jbess(r1am1*t,1,111h(i),temp)
call jbess(rlam5*t,1,115(i),temp)
call jbess(rlam1*t,2,i21h(i),temp)

call jbess(r1am5*t,2,i25h(i),temp)

i11h(i) = i11h(i)/t

115h(i) = i15(i)/t

i21h(i) = i21h(i)/t

i25h(i) = i25h(i)/t
end do

call conv(i11h,i15h,nt,dt,i111)
do i=1,nt
ft(i) = expon(i)*(rmur*ct+rkap)*(-r1am1*r1am5*ilil(i)+

2 rlam1*rlam1*i21h(i)+rlam5*rlam5*i25h(i))

 

 

elseif((om0*om0).1t.(delta*delta-(rmur*b*b/(rkap*rkap)))) then

126

write(*,*) "CASE 2"

rlam1 = dsqrt(delta*delta-om0*om0)
rlam5 = dsqrt(de1ta*de1ta-om0*omO-rmur*b*b/(rkap*rkap))
do i=1,nt

t = (i-1)*dt+.00001*dt

111(1) = BESSI(1,rlam1*t)
115(1) = BESSI(1,r1am5*t)
i21h(i) = BESSI(2,rlam1*t)
i25h(i) = BESSI(2,rlam5*t)
111h(i) = 111(i)/t

115h(i) = i15(i)/t

i21h(i) = i21h(i)/t
125h(i) = i25h(i)/t

end do

call conv(111h,115h,nt,dt,1111)
do i=1,nt
ft(i) = expon(i)*(rmur*ct+rkap)*(-r1am1*rlam5*ilii(i)+

2 r1am1*rlam1*i21h(i)+r1am5*rlam5*i25h(i))

write(*,*) "CASE 3"
rlaml = dsqrt(delta*delta-om0*om0)
rlam5 = dsqrt(om0*om0+(rmur*b*b/(rkap*rkap))-delta*delta)

do i=1,nt

127

t = (i-1)*dt+.00001*dt

i11h(i) = BESSIC1,rlam1*t)

call jbess(rlam5*t,1,115(i),temp)
121h(i) = BESSI(2,r1am1*t)

call jbess(r1am5*t,2,i25h(i),temp)

111h(i) = 111h(i)/t

115h(1) = 115(i)/t

121h(i) = 121h(i)/t

125h(1) = 125h(i)/t
end do

call conv(i11h,115h,nt,dt,1111)
do i=1,nt
ft(i) = expon(i)*(rmur*ct+rkap)*(rlam1*rlam5*1111(i)+
2 rlam1*rlam1*121h(1)+rlam5*rlam5*125h(1))
end do

endif

 

YY = (rmur*b*b)*(rkap-rmur*ct)/(2.d0*rkap*rkap)
if (ct.eq.rkap/rmur) then
write(*,*) ’Special Denom’
do i=1,nt
gammat(i) = -f1*ft(i)/(rmur*b*b)
end do
elseif (c3case3.eq.0) then

call conv(c3,ft,nt,dt,c3ft)

128

do i=1,nt
gammat(i) = f2*(c3ft(i)+YY*c3(i))

end do

write(*,*) ’non-causal’
do i=1,(nt/shift)
ftshift(1) = 0
end do
do i=(nt/shift+1),nt
ftshift(i) = ft(i-nt/shift)
end do
call conv(c3,ftshift,nt,dt,c3ft)
do i=1,(nt/Sh1ft+1)
C38h1ft(i) = 0
end do
do i=(nt/shift+2),nt
c3shift(1) = c3(i-nt/shift-1)
end do
do i=1,nt
gammat(i) = f2*(c3ft(i)+YY*c33hift(i))

and do

endif

open(10,file=’gammat.dat’,status=’unknown’)

do i=1,nt

if (c3case3.eq.0) then

write(10,*) ((1-1)*dt+.00001*dt)*1.d9,gammat(i)

129

else
write(10,*) ((i-1-2*nt/(shift))*dt+.00001*dt)*1.d9,gammat(i)
endif
end do

close(10)

end

subroutine conv(f,g,n,de1,c)

performs convolution of waveforms f and g using

linear interpolation

parameter (nsize=10000)
parameter (€13=0.33333333333333333333d0,

2 c16=0.16666666666666666666d0)
implicit real*8 (a-h,o-z)
real*8 f(nsize),g(nsize),c(nsize)
c(1)=0.d0
do k=2,n

sum = 0.d0

do 1=1,k-1

sum = sum + c13*f(1+1)*g(k-1)

2 + c13*f(l)*g(k-1+1)
3 + c16*f(1)*g(k-1)
4 + c16*f(1+1)*g(k-1+1)

130

end do

C(k) = del*sum
end do
return

end

FUNCTION BESSI(N,X)

This subroutine calculates the first kind modified Bessel function
of integer order N, for any REAL X. We use here the classical
recursion formula, when X > N. For X < N, the Miller’s algorithm
is used to avoid overflows.

REFERENCE:

C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS,

MATHEMATICAL TABLES, VOL.5, 1962.

PARAMETER (IACC = 40,BIGNO = 1.D10, BIGNI = 1.D-10)
REAL *8 X,BESSI,BESSIO,BESSI1,TOX,BIM,BI,BIP

IF (N.EQ.0) THEN

BESSI = BESSIO(X)

RETURN

ENDIF

IF (N.EQ.1) THEN

BESSI = BESSI1(X)

RETURN

131

ENDIF
IF(X.EQ.0.DO) THEN
BESSI=0.D0

RETURN

ENDIF

TOX = 2.D0/X

BIP = 0.DO

BI = 1.D0

BESSI 0.00
M = 2*((N+INT(SQRT(FLOAT(IACC*N)))))

DO 12 J = M,1,-1

BIM = BIP+DFLOAT(J)*TOX*BI
BIP = BI
BI = BIM

IF (ABS(BI).GT.BIGNO) THEN

BI

BI*BIGNI

BIP BIP*BIGNI
BESSI = BESSI*BIGNI
ENDIF
IF (J.EQ.N) BESSI = BIP

12 CONTINUE
BESSI = BESSI*BESSIOTX)/BI
RETURN

END

! Auxiliary Bessel functions for N=O, N=1

FUNCTION BESSIO(X)

132

REAL *8 X,BESSIO,Y,P1,P2,P3,P4,P5,P6,P7,

* Ql,Q2,Q3,Q4,QS,06,Q7,Q8,QQ,AX,BX

DATA P1,P2,P3,P4,P5,P6,P7/1 D0,3.5156229D0,3 0899424DO,1.2067429DO
* ,0.2659732D0,0.360768D-1,0.45813D-2/

DATA Ql,Q2,Q3,Q4,QS,Q6,Q7,QB,QQ/O.39894228D0,0.1328592D-1,
* 0.225319D-2,-0.157565D-2,0.916281D-2,-O.2057706D-1,

* O.2635537D-1,-0.1647633D-1,0.392377D-2/
IF(ABS(X).LT.3.75DO) THEN

Y=(X/3.75DO)**2
BESSIO=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))

ELSE

AX=ABS(X)

Y=3.75DO/AX

BX=EXP(AX)/SQRT(AX)
AX=Ql+Y*(Q2+Y*(QB+Y*(Q4+Y*(QS+Y*(06+Y*(Q7+Y*(08+Y*QQ)))))))
BESSIO=AX*BX

ENDIF

RETURN

END

FUNCTION BESSIICX)

REAL *8 X,BESSI1,Y,P1,P2,P3,P4,P5,P6,P7,
* Ql,Q2,QS,Q4,QS,Q6,Q7,Q8,09,AX,BX

DATA P1,P2,P3,P4,P5,P6,P7/0.5D0,0.87890594D0,0.51498869D0,
* 0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/

DATA Ql,Q2,03,Q4,Q5,06,Q7,Q8,09/0.39894228D0,-0.3988024D-1,

* -0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1,

133

* —0 2895312D-1,0 1787654D-1,-O.420059D-2/
IF(ABS(X).LT.3.75DO) THEN

Y=(X/3.75D0)**2
BESSII=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
ELSE

AX=ABS(X)

Y=3.75DO/AX

BX=EXP(AX)/SQRT(AX)
AX=Ql+Y*(Q2+Y*(Q3+Y*(Q4+Y*(05+Y*(Q6+Y*(Q7+Y*(08+Y*Q9)))))))
BESSIl=AX*BX

ENDIF

RETURN

END

subroutine jbess (x,n,bj,bjp)

calculates jn and jn’ for n positive or negative

implicit real*8 (a-h,o-z)

if (n .eq. 0) then

bj = bessj(0,x)
bjp = -bessj(1,x)
return

endif

134

if (n .ge. 0) then

endif

if (x .eq. 0 d0) then

bj = 0.d0

bjp1 = bessj(m+1,x)

bjml = bessj(m-1,x)

bjp = (bjml-bjp1)/2 d0
else

bj = bessj(m,x)

bj1 = bessj(m+1,x)

bjp = -bj1 + n*bj/x
endif

if (n .1t. 0) then

bj = bj*(-1)**n
bjp = bjp*<-1>**n
endif
return
end

135

FUNCTION BESSJO(X)

implicit rea1*8 (a-h,o-z)

REAL*8 Y,P1,P2,P3,P4,P5,Ql,Q2,QB,Q4,QS,R1,R2,R3,R4,R5,R6,
* Sl,S2,SB,S4,SS,S6

DATA P1,P2,P3,P4,P5/1.D0,-.1098628627D-2,.2734510407D-4,

* -.2073370639D-5,.2093887211D-6/. Ql,Q2,Q3,Q4,Q5/-.15624999950-
*1,
* .1430488765D-3,-.69111476510‘5,.7621095161D-6,-.9349451520-7/

DATA R1,R2,R3,R4,R5,R6/57568490574.DO,-13362590354.DO,651619640.7D
*0,
* -11214424.18D0,77392.33017D0,-184.9052456D0/.
* S1,S2,S3,S4,S5,S6/57568490411.DO,1029532985.D0,
* 9494680.718D0,59272.64853DO,267.8532712D0,1.DO/
IF(ABS(X).LT.8.)THEN
Y=X**2
BESSJO=<R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))
* /(Sl+Y*(S2+Y*(SB+Y*(S4+Y*(SS+Y*SG)))))
ELSE
AX=ABS(X)
Z=8./AX
Y=Z**2
xx=Ax-.785398164
BESSJO=SQRT(.636619772/AX)*(COS(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y
* *P5))))-Z*SIN(XX)*(Ql+Y*(QQ+Y*(Q3+Y*(Q4+Y*QS)))))
ENDIF

RETURN

136

END

FUNCTION BESSJICX)

implicit rea1*8 (a-h,o-z)

REAL*8 Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,QS,R1,R2,R3,R4,R5,R6,

* SI,S2,S3,S4,SS,S6

DATA R1,R2,R3,R4,R5,R6/72362614232.DO,-7895059235.D0,242396853.1DO

*

2

* -2972611.439D0,15704.4826ODO,-30.16036606D0/,
* Sl,S2,S3,34,35,36/144725228442.D0,2300535178.DO,
* 18583304.?4D0,99447.43394D0,376.9991397D0,1.DO/

DATA P1,P2,P3,P4,P5/1.D0,.183105D-2,-.3516396496D-4,.2457520174D-5

*

9

* -.240337019D-6/, Ql,Q2,QB,Q4,QS/.04687499995D0,-.2002690873D-3

*

* .8449199096D-5,-.88228987D-6,.105787412D-6/
IF(ABS(X).LT.8.)THEN

Y=X**2

BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))
* /(Sl+Y*(82+Y*(33+Y*(S4+Y*(SS+Y*SG)))))
ELSE

AX=ABS(X)

Z=8./AX

Y=Z**2

XX=AX-2.356194491

137

BESSJ1=SQRT(.636619772/AX)*(COS(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y
* *PS))))-Z*SIN(XX)*(QI+Y*(Q2+Y*(QS+Y*(Q4+Y*QS)))))
* *SIGN(1.,X)
ENDIF
RETURN

END

FUNCTION BESSJ(N,X)
implicit real*8 (a-h,o-z)
PARAMETER (IACC=40,BIGNO=1.d10,BIGNI=1.d-10)

if (x .ge. 0.d0) then

ff = 1.d0
else

ff = (-1.d0)**n
end if

xx = abs(x)
if (n .eq. 0) then
bessj = ff*bessj0(xx)
return
endif
if (n .eq. 1) then
bessj = ff*bessj1(xx)
return
endif

TOX=2./Xx

138

IF(Xx.GT.FLOAT(N))THEN
BJM=BESSJO(XX)
BJ=BESSJ1(XX)

DO 11 J=1,N-1
BJP=J*TOX*BJ-BJM
BJM=BJ
BJ=BJP

CONTINUE

BESSJ=ff*BJ

ELSE
M=2*((N+INT(SQRT(FLOAT(IACC*N))))/2)
BESSJ=0.

JSUM=0

SUM=0.

BJP=0.

BJ=1.

DO 12 J=M,1,-1
BJM=J*TOX*BJ-BJP
BJP=BJ
BJ=BJM
IFCABS(BJ).GT.BIGNO)THEN

BJ=BJ*BIGNI

BJP=BJP*BIGNI

BESSJ=BESSJ*BIGNI

SUM=SUM*BIGNI
ENDIF

IF(JSUM.NE.O)SUM=SUM+BJ

139

JSUM=1-JSUM
IF(J . EQ . N)BESSJ=BJP
12 CONTINUE
SUM=2 . *SUM-BJ
BESSJ=ff*BESSJ/SUM
ENDIF
RETURN

END

140

IAPUPEHNIDLXLI)

PWDIUTILAJV CKDIJE)FT)FL(3EHHEH&AHLTERATCUXSE)

program genera1_TM

implicit none

integer nsize

parameter(nsize = 50000)

real*8 BESSI,om0,delta,b,theta,dt,ct,st,tt,t,pi,df,f,om,YY
real*8 rmur,epsinf,rkap,f1,sigma,chi

real*8 rlaml,rlam5,rlamc,rlame,expon(nsize),gams,temp
real*8 111h(nsize),i15h(nsize),115(nsize),111(nsize)
real*8 i21h(nsize),i25h(nsize),ili1(nsize),101(nsize)
real*8 c3ft(nsize),c38hift(nsize),ftshift(nsize)

real*8 c3(nsize),ft(nsize),gammat(nsize),iOii(nsize)
real*8 c1(nsize),c2(nsize),ta(nsize),scale

integer i,nt,nf,c1nc,c2nc,shift,c3case3

complex*16 zl,z2,s,zj,gamma,zzl,zz2,sl,s2,sS,s6

rmur = 1.d0
epsinf = 2.d0

omO = 4.d11

141

b = 6.24d11
delta = 2.5d10
theta = 60.d0

pi = 4.d0*atan(1.d0)

ct = cos(theta*pi/18O d0)
st = sin(theta*pi/180.d0)
tt = tan(theta*pi/180.d0)
zj = dcmplx(0.d0,1.d0)
c1nc = O

c2nc = 0

c3case3 = 0

scale = 1.d-2O

rkap dsqrt(rmur*epsinf-st*st)

gams b*b*(rmur-2.d0*epsinf*ct*ct)/(rkap*rkap-epsinf*epsinf*ct*ct)
sigma = dsqrt(rmur*rmur-sin(2.dO*theta*pi/180.dO)*

2 sin(2.dO*theta*pi/180.d0))

chi = (b*b*sigma)/(2.d0*(rkap*rkap-epsinf*epsinf*ct*ct))

shift = 8

YY b*b*(rkap-epsinf*ct)*(rmur*epsinf-2.d0*rkap*rkap)/

2 (2.d0*rkap*rkap*epsinf)

f1 (2.d0*rkap*epsinf*ct)/((rkap+epsinf*ct)*(rkap+epsinf*ct)*

2 (rkap-epsinf*ct))

df 0.4d0

nf 4095*8

142

4096*2

nt

dt 1d-13
221 = cdsqrt(dcmplx(delta*delta-om0*om0,0.dO))

sl = -de1ta+zzl

s2 = -de1ta-zzl

222 = cdsqrt(dcmplx(delta*de1ta-om0*omO-rmur*b*b/(rkap*rkap),0.d0))
$5 = -delta+zz2

$6 = -delta-zz2

 

open (10,file=’gamma.dat’,status=’unknown’)
do i=1,nf

f = df*i

om = 2.d0*pi*f*1.d9

s = zj*om

21 cdsqrt(s-sl)*cdsqrt(s-s2)*cdsqrt(s-sS)*cdsqrt(s-s6)

z2 epsinf*ct*((s-sl)*(s-s2)+b*b/epsinf)
gamma = (rkap*zl-z2)/(rkap*zl+z2)-(rkap-epsinf*ct)/(rkap+epsinf*ct)
write (10,*) f,rea1(gamma),aimag(gamma)

end do

close (10)

do i=1,nt

t = (1-1)*dt+.00001*dt

143

expon(i) = exp(-delta*t)

end do

if(om0*om0+(gams/2 dO)-chi.gt.(delta*delta)) then
write(*,*) "Cl term 1“
rlamc = dsqrt(omO*om0+(gams/2.dO)-chi-delta*delta)
do i=1,nt
t = (i-1)*dt+.00001*dt
c1(i) = (sigma-rmur+2.d0*epsinf*ct*ct)*expon(i)*
2 sin(r1amc*t)/(2.d0*sigma*rlamc)
end do
elseif(omO*om0+(gams/2.dO)-chi.1t.0) then
write(*,*) "C1 term 3"
rlamc = dsqrt(-om0*omO-(gams/2.d0)+chi+de1ta*delta)
c1nc = 1
c3case3 = 1
do i=1,nt/Shift
t = (i-1-nt/shift)*dt+.00001*dt
c1(i) = -(sigma-rmur+2.d0*epsinf*ct*ct)*
2 exp((-delta+rlamc)*t)/(4.d0*sigma*rlamc)
end do
do i=nt/shift+1,nt
t = (i-1-nt/shift)*dt+.00001*dt

c1(i) = -(sigma-rmur+2.d0*epsinf*ct*ct)*

144

2 exp(-(delta+rlamc)*t)/(4.d0*sigma*rlamc)

end do

else
write(*,*) "Cl term 2"
rlamc = dsqrt(-om0*omO-(gams/2.d0)+chi+delta*delta)
do i=1,nt
t = (i-l)*dt+.00001*dt
c1(i) = (sigma-rmur+2.d0*epsinf*ct*ct)*expon(i)*
2 sinh(rlamc*t)/(2.d0*sigma*rlamc)

end do

 

if(om0*om0+(gams/2.d0)+chi.gt.(delta*de1ta)) then
write(*,*) "C2 term 1"
rlame = dsqrt(omO*omO+(gams/2.dO)+chi-delta*delta)
do i=1,nt
t = (i-1)*dt+.00001*dt
c2(i) = (sigma+rmur-2.d0*epsinf*ct*ct)*expon(i)*
2 sin(rlame*t)/(2.dO*sigma*rlame)
end do
elseif(omO*om0+(gams/2.d0)+chi.1t.0) then
write(*,*) "C2 term 3"

rlame = dsqrt(-om0*omO-(gams/2.d0)-chi+de1ta*delta)

145

c2nc = 1
c3case3 = 1
do i=1,nt/shift
t = (i-l-nt/shift)*dt+.00001*dt
c2(i) = -(sigma+rmur-2.d0*epsinf*ct*ct)*
2 exp((-delta+rlame)*t)/(4 d0*sigma*r1ame)
end do
do i=nt/shift+1,nt
t = (1-1-nt/shift)*dt+.00001*dt
c2(i) = -(sigma+rmur-2.d0*epsinf*ct*ct)*
2 exp(-(delta+rlame)*t)/(4.d0*sigma*r1ame)

end do

else
write(*,*) "C2 term 2"
rlame = dsqrt(-om0*omO-(gams/2.dO)-ch1+delta*delta)
do i=1,nt
t = (i-1)*dt+.00001*dt
c2(i) = (sigma+rmur-2.d0*epsinf*ct*ct)*expon(i)*
2 sinh(rlame*t)/(2 d0*sigma*rlame)

end do

if (c3case3.ne.0) then

146

if (c1nc.eq.0) then
write(*,*) ’shift c1’
ta = Cl;
do i=1,(nt/shift)
c1(i) = 0
end do
do i=(nt/shift+1),nt
c1(i) = ta(i-nt/shift)
end do
endif
if (c2nc.eq.0) then
write(*,*) ’shift c2’
ta = c2;
do i=1,(nt/shift)
c2(i) = 0
end do
do i=(nt/shift+1),nt
c2(i) = ta(i-nt/shift)
end do
endif

endif

do i=1,nt

c3(i) = c1(i)+c2(i)

end do

147

 

Case 1: om0‘2>delta“2

if((om0*om0).gt.(delta*delta)) then

write(*,*) "CASE 1"

r1am1 = dsqrt(omO*omO-delta*delta)
rlam5 = dsqrt(om0*om0+(rmur*b*b/(rkap*rkap))-delta*delta)
do i=1,nt

t = (i‘l)*dt+.00001*dt

call jbess(r1am1*t,0,101(i),temp)
call jbess(r1am1*t,1,111(i),temp)
call jbess(r1am5*t,1,115h(i),temp)
call jbess(rlam1*t,2,i21h(i),temp)

call jbess(r1am5*t,2,i25h(i),temp)

111h(1) = 111(i)/t

115h(1) = 115h(i)/t

i21h(i) = i21h(i)/t

125h(i) = 125h(1)/t
end do

call conv(iOl,ilSh,nt,dt,iOil)
call conv(111h,115h,nt,dt,1111)
do i=1,nt
t = (i-1)*dt+.00001*dt
ft(i) = expon(i)*((epsinf*ct+rkap)*(-r1am1*rlam5*1111(i)-

2 (b*b/epsinf)*rlam5*1011(i)+(b*b/epsinf)*rlam1*111(i)+

148

3 r1am1*rlam1*i21h(i)+r1am5*r1am5*i25h(i))+

4 (b*b*b*b*ct/(epsinf*rlam1))*sin(rlam1*t))

 

elseif((om0*om0).1t.(delta*de1ta-(rmur*b*b/(rkap*rkap)))) then
write(*,*) "CASE 2"
rlaml = dsqrt(delta*delta-om0*om0)
rlam5 = dsqrt(delta*de1ta-om0*omO-rmur*b*b/(rkap*rkap))
do i=1,nt
t = (i-1)*dt+.00001*dt

101(1)

BESSI(O,r1am1*t)

111(1)

BESSI(1,r1am1*t)
115h(i) = BESSI(1,r1am5*t)
i21h(i) = BESSI(2,r1am1*t)
i25h(i) = BESSI(2,rlam5*t)
111h(i) = 111(1)/t
115h(i) = 115h(i)/t
i21h(i) = i21h(i)/t
i25h(i) = i25h(1)/t

end do

call conv(iOl,115h,nt,dt,1011)

call conv(illh,115h,nt,dt,1111)

do i=1,nt
t = (i-1)*dt+.00001*dt

ft(i) = expon(i)*((epsinf*ct+rkap)*(-r1am1*rlam5*ilil(i)+

149

2 (b*b/epsinf)*r1am5*i0i1(i)-(b*b/epsinf)*rlam1*111(i)+
3 rlam1*r1am1*i21h(i)+rlam5*rlam5*i25h(i))+

4 (b*b*b*b*ct/(epsinf*rlam1))*sinh(rlam1*t))

else
write(*,*) "CASE 3"
r1am1 = dsqrt(delta*delta-om0*om0)
rlam5 = dsqrt(om0*om0+(rmur*b*b/(rkap*rkap))-delta*delta)
do i=1,nt
t = (i-1)*dt+.00001*dt

101(1)

BESSI(0,rlam1*t)

111(1) BESSI(1,r1am1*t)
call jbess(rlam5*t,1,115h(i),temp)
121h(i) = BESSI(2,r1am1*t)

call jbess(r1am5*t,2,i25h(i),temp)

111h(i) = 111(1)/t

115h(1) = 115h(1)/t

i21h(i) = i21h(i)/t

i25h(i) = 125h(1)/t
end do

call conv(101,115h,nt,dt,iOi1)
call conv(111h,115h,nt,dt,1111)
do i=1,nt

t = (i-1)*dt+.00001*dt

150

ft(i) = expon(i)*((epsinf*ct+rkap)*(rlam1*r1am5*ilil(i)-

2 (b*b/epsinf)*rlam5*1011(i)-(b*b/epsinf)*rlam1*i11(i)+
3 rlam1*r1am1*i21h(i)+rlam5*r1am5*i25h(1))+
4 (b*b*b*b*ct/(epsinf*r1am1))*sinh(r1am1*t))
end do
endif

if (c3case3.eq.0) then
call conv(c3,ft,nt,dt,c3ft)
do i=1,nt
gammat(i) = f1*(c3ft(1)+YY*c3(i))
end do
else
write(*,*) ’non-causal’
do i=1,(nt/shift)
ftshift(i) = 0
end do
do i=(nt/shift+1),nt
ftshift(i) = ft(i-nt/shift)
end do
call conv(c3,ftshift,nt,dt,c3ft)
do i=1,(nt/shift+1)
c3shift(1) = 0
end do

do i=(nt/shift+2),nt

151

c38hift(i) = c3(i-nt/shift-1)
end do
do i=1,nt

gammat(i) = f1*(c3ft(i)+YY*c3shift(1))
end do

endif

open(10,file=’gammat.dat’,status=’unknown’)
do i=1,nt
if (c3case3.eq.0) then
write(10,*) ((i-1)*dt+.00001*dt)*1.d9,gammat(i)
else
write(10,*) ((1-1-2*nt/(shift))*dt+.00001*dt)*1.d9,gammat(i)
endif
end do

close(10)

end

subroutine conv(f,g,n,de1,c)

performs convolution of waveforms f and g using

linear interpolation

parameter (nsize=10000)

parameter (c13=0.33333333333333333333d0,

152

2 c16=0.16666666666666666666d0)
implicit real*8 (a-h,o-z)
real*8 f(nsize),g(nsize),c(nsize)
c(1)=0.d0
do k=2,n

sum = 0.dO

do l=1,k-1

sum = sum + c13*f(l+1)*g(k-l)

2 + c13*f(l)*g(k-1+1)

3 + c16*f(1)*g(k-1)

4 + c16*f(1+1)*g(k-1+1)
end do

c(k) = del*sum
end do
return

end

FUNCTION BESSI(N,X)

This subroutine calculates the first kind modified Bessel function
of integer order N, for any REAL X. We use here the classical
recursion formula, when X > N. For X < N, the Miller’s algorithm
is used to avoid overflows.

REFERENCE:

C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS,

153

MATHEMATICAL TABLES, VOL.5, 1962.

PARAMETER (IACC = 40,BIGNO = 1.D10, BIGNI = 1.D-10)
REAL *8 X,BESSI,BESSIO,BESSI1,TOX,BIM,BI,BIP
IF (N.EQ.0) THEN

BESSI = BESSIO(X)

RETURN

ENDIF

IF (N.EQ.1) THEN

BESSI = BESSI1(X)

RETURN

ENDIF

IF(X.EQ.0.D0) THEN

BESSI=0.D0

RETURN

ENDIF

TOX = 2.D0/X

BIP = 0.D0

BI = 1.D0

BESSI 0.D0

M = 2*((N+INT(SQRT(FLOAT(IACC*N)))))

DO 12 J = M,1,-1

BIM = BIP+DFLOAT(J)*TOX*BI
BIP = BI
BI = BIM

IF (ABS(BI).GT.BIGNO) THEN

BI = BI*BIGNI

154

BIP = BIP*BIGNI
BESSI = BESSI*BIGNI
ENDIF
IF (J.EQ.N) BESSI = BIP
12 CONTINUE
BESSI = BESSI*BESSIO(X)/BI
RETURN

END

! Auxiliary Bessel functions for N=O, N=1
FUNCTION BESSIO(X)
REAL *8 X,BESSIO,Y,P1,P2,P3,P4,P5,P6,P7,
* 01,Q2,03,Q4,QS,06,Q7,Q8,09,AX,BX
DATA P1,P2,P3,P4,P5,P6,P7/1.D0,3.5156229D0,3.0899424DO,1.2067429DO
* ,0.2659732D0,0.360768D-1,0.45813D-2/
DATA Ql,Q2,03,Q4,QS,QG,Q7,08,QQ/O.39894228D0,0.1328592D-1,
* 0.225319D-2,-O.157565D-2,0.916281D-2,-0.2057706D-1,
* 0.2635537D-1,-0.16476330-1,0.392377D-2/
IF(ABS(X).LT.3.75DO) THEN
Y=(X/3.75D0)**2
BESSIO=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
ELSE
AX=ABS(X)
Y=3.75DO/AX
BX=EXP(AX)/SQRT(AX)
AX=01+Y*(02+Y*(Q3+Y*(Q4+Y*(QS+Y*(06+Y*(Q7+Y*(Q8+Y*QQ)))))))

BESSIO=AX*BX

155

ENDIF

RETURN

END

FUNCTION BESSI1(X)

REAL *8 X,BESSI1,Y,P1,P2,P3,P4,P5,P6,P7,
* Ql,Q2,03,Q4,QS,06,Q7,08,09,AX,BX

DATA P1,P2,P3,P4,P5,P6,P7/0.5D0,0.87890594D0,0.51498869D0,
* 0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/
DATA Ql,Q2,QB,Q4,05,06,07,Q8,09/0.39894228D0,-0.3988024D-1,
* -0.362018D-2,0.1638010-2,-0.1031555D-1,0.2282967D-1,
* -0.2895312D-1,0.1787654D-1,-0.420059D-2/
IF(ABS(X).LT.3.75DO) THEN

Y=(X/3.75D0)**2
BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))

ELSE

AX=ABS(X)

Y=3.75D0/AX

BX=EXP(AX)/SQRT(AX)
AX=Ql+Y*(QZ+Y*(QB+Y*(Q4+Y*(QS+Y*(QG+Y*(Q7+Y*(08+Y*09)))))))
BESSI1=AX*BX

ENDIF

RETURN

END

156

subroutine jbess (x,n,bj,bjp)

c calculates jn and jn’ for n positive or negative

c
implicit real*8 (a-h,o-z)
c
if (n .eq. 0) then
bj = bessj(0,x)
bjp = -bessj(1,x)
return
endif
c
if (n .ge. 0) then
m = n
else
m = -n
endif
c

if (x .eq. 0.d0) then

bj = O.dO
bjp1 = bessj(m+1,x)
bjml = bessj(m-1,x)
bjp = (bjml-bjp1)/2.d0
else
bj = bessj(m,x)
bj1 = bessj(m+1,x)
bjp = -bj1 + n*bj/x

157

endif

if (n .1t. 0) then

bj = bj*(-1)**n
bjp = bjp*<-1>**n
endif
return
end

FUNCTION BESSJO(X)

implicit real*8 (a-h,o-z)

REAL*8 Y,P1,P2,P3,P4,P5,Ql,Q2,QB,Q4,QS,R1,R2,R3,R4,R5,R6,
* S1,S2,SB,S4,SS,S6

DATA P1,P2,P3,P4,P5/1.D0,-.10986286270-2,.2734510407D-4,

* -.2073370639D-5,.2093887211D-6/, QI,Q2,QB,Q4,QS/-.1562499995D-
*1,
* .14304887650-3,-.6911147651D-5,.7621095161D-6,-.934945152D-7/

DATA R1,R2,R3,R4,R5,R6/57568490574.DO,-13362590354.D0,651619640.7D
*0,
* -11214424.18D0,77392.33017D0,-184.9052456D0/,
* Sl,32,S3,S4,S5,86/57568490411.D0,1029532985.D0,
* 9494680.718D0,59272.64853D0,267.8532712D0,1.DO/
IF(ABS(X).LT.8.)THEN

Y=X**2

158

BESSJO=(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))
* /(Sl+Y*(S2+Y*(SB+Y*(S4+Y*(SS+Y*SG)))))
ELSE
AX=ABS(X)
Z=8./Ax
Y=Z**2
XX=AX-.785398164
BESSJO=SQRT(.636619772/AX)*(COS(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y
* *Ps))))-Z*SIN(XX)*(01+Y*(Q2+Y*(QB+Y*(Q4+Y*QS)))))
ENDIF
RETURN

END

FUNCTION BESSJ1(X)

implicit real*8 (a-h,o-z)

REAL*8 Y,P1,P2,P3,P4,P5,Ql,Q2,QB,Q4,QS,R1,R2,R3,R4,R5,R6,

* S1,S2,83,S4,S5,36

DATA R1,R2,R3,R4,R5,R6/72362614232.D0,-7895059235.D0,242396853.1DO

a:

* -2972611.439D0,15704.48260D0,-30.16036606D0/,
* Sl,S2,S3,S4,SS,S6/144725228442.D0,2300535178.D0,
* 18583304.74D0,99447.43394D0,376.9991397D0,1.D0/

DATA P1,P2,P3,P4,P5/1.D0,.183105D-2,-.3516396496D-4,.2457520174D-5

*

3

* -.240337019D-6/, Ql,Q2,Q3,Q4,QS/.04687499995D0,-.2002690873D-3

* .8449199096D-5,- 88228987D-6,.105787412D-6/
IF(ABS(X).LT.8.)THEN
Y=X**2
BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))
* /(SI+Y*(82+Y*(SB+Y*(S4+Y*(SS+Y*SG)))))
ELSE
AX=ABS(X)
Z=8./AX
Y=Z**2
XX=AX-2.356194491
BESSJ1=SQRT(.636619772/AX)*(COS(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y
* *P5))))-Z*SIN(XX)*(01+Y*(Q2+Y*(QB+Y*(Q4+Y*QS)))))
* *SIGN(1.,X)
ENDIF
RETURN

END

FUNCTION BESSJ(N,X)
implicit real*8 (a-h,o-z)
PARAMETER (IACC=40,BIGNO=1.le,BIGNI=1.d-lO)

if (x .ge. O.d0) then

ff = 1.d0
else
If = (-1.d0)**n

160

11

end if
xx = abs(x)
if (n .eq. 0) then
bessj = ff*bessj0(xx)
return
endif
if (n .eq. 1) then
bessj = ff*bessj1(xx)
return
endif
TOX=2./Xx
IF(XX.GT.FLOAT(N))THEN

BJM=BESSJO(Xx)

BJ=BESSJ1(Xx)

DO 11 J=1,N-1
BJP=J*TOX*BJ-BJM
BJM=BJ
BJ=BJP

CONTINUE

BESSJ=ff*BJ

ELSE

M=2*((N+INT(SQRT(FLOATCIACC*N))))/2)

BESSJ=0.

JSUM=0

SUM=0.

BJP=0.

BJ=1.

161

DO 12 J=M,1,-1
BJM=J*TOX*BJ-BJP
BJP=BJ
BJ=BJM
IF(ABS(BJ).GT.BIGNO)THEN

BJ=BJ*BIGNI
BJP=BJP*BIGNI
BESSJ=BESSJ*BIGNI
SUM=SUM*BIGNI
ENDIF
IF(JSUM.NE.0)SUM=SUM+BJ
JSUM=1-JSUM
IF(J.EQ.N)BESSJ=BJP
12 CONTINUE
SUM=2.*SUM-BJ
BESSJ=ff*BESSJ/SUM
ENDIF
RETURN

END

162

BIBLIOGRAPHY

163

BIBLIOGRAPHY

[1] RM. Joseph, S.C. Hagness, and A. Taflove, “Direct time integration of Maxwell’s
equations in linear dispersive media with absorption for scattering and propaga-
tion of femtosecond electromagnetic pulses,” Opt. Lett. 16, 1412-1414 (1991).

[2] KB. Oughstun and CC. Sherman, “Propagation of electromagnetic pulses in a
linear dispersive medium (the Lorentz medium),” J. Opt. Soc. Am. B 5, 817-849
(1988).

[3] KB. Oughstun and CC. Sherman, “Uniform asymptotic description of electro-
magnetic pulse propagation in a linear dispersive medium with absorption (the
Lorentz medium),” J. Opt. Soc. Am. A 6, 1394-1420 (1989).

[4] EL. Mokole and SN. Samaddar, “Transmission and reflection of normally inci-

dent, pulsed electromagnetic plane waves upon a Lorentz half—space,” J. Opt.
Soc. Am. B 16, 812-831 (1999).

[5] J .G. Blaschak and J. Franzen, “Precursor propagation in dispersive media from
short-rise-time pulses at oblique incidence,” J. Opt. Soc. Am. B 12, 1501-1512
(1995)

[6] OJ Stenholm, E.J. Rothwell, D.P. Nyquist, L.C. Kempel, and LL. Frasch, “E-
Pulse diagnostics of simpled layered materials,” IEEE Trans. Antennas Propag.
51, 3221-3227 (2003).

[7] H.-Y. Pao, S.L. Dvorak, and DC. Dudley, “An accurate and efficient analysis for
transient plane waves obliquely incident on a conductive half space (TM case),”
IEEE Trans. Antennas Propag. 44, 925-932 (1996).

[8] J .C. Oh, E. Rothwell, B.T. Perry, and M.J. Havrilla, “Natural Resonance Repre-
sentation of the 'ITansient Field Reflected by a Conductor-Backed Layer of Debye
Material, J. Electromagnetic Waves and Applications 18, 571-589 (2004).

[9] J .W. Suk and E.J. Rothwell, “Transient Analysis of TM-Plane Wave Reflection
from a Layered Medium,” J. Electromagnetic Waves and Applications 16, 1195-
1208 (2002).

[10] KG. Gray, “The reflected impulse response of a Lorentz medium,” Proc. IEEE
68, 408-409 (1980).

164

[11] RV. Stanic, D.R. Milanovic, and J.M. Cvetic, “Pulse reflection from a lossy
Lorentz medium half-space (TM polarization),” J. Phys. D: Appl. Phys. 24, 1245-
1249 (1991).

[12] E.J. Rothwell and M.J. Cloud, Electromagnetics (CRC Press, Boca Raton, 2001).

[13] VV.R. LePage, Complex Variables and the Laplace "Transform for Engineers
(Dover Press, New York, 1980).

[14] GA. Campbell and RM. Foster, Fourier Integrals for Practical Applications,
2nd ed. (D. Van Nostrand Co., New York, 1951).

[15] M. Abramowitz and 1.3. Stegun, Handbook of Mathematical Functions (Dover,
New York, 1965).

[16] L. Brillouin, Wave Propagation and Group Velocity (Academic Press, New York,
1960).

165

MICHIGAN STATE U

[ll/[[N/Mllllllfl]]l]]]]]]]]][]l

3 1293 02845 716