NBSIR 81-1653 publications

S'tt®

NBS ME

NAT'L INST. OF STAND & TECH

AlllOb TTfl'Tfl 4

MODELLING OF OIL SHALE RETORTS FOR

ELECTROMAGNETIC SENSING TECHNIQUES

H. Chew

Electromagnetic Fields Division

National Engineering Laboratory

National Bureau of Standards

Boulder, Colorado 80303

Prepared for:

U.S. Department of Energy

Laramie Energy Technology Center

Laramie, Wyoming 82071

'QO—

100

.1156

I 81-1653

1981

Lc. 2

November 1981

NBSIR 81-1653

4TXOHAL BUREAl

OV STANDARDS

LIBRARY

JAN 2 5 1982

not <<cc. - a

MODELLING OF OIL SHALE RETORTS FOR

ELECTROMAGNETIC SENSING TECHNIQUES

c.>

H. Chew

Electromagnetic Fields Division

National Engineering Laboratory

National Bureau of Standards

Boulder, Colorado 80303

November 1981

Prepared for:

U.S. Department of Energy

Laramie Energy Technology Center

Laramie, Wyoming 82071

U.S. DEPARTMENT OF COMMERCE, Malcolm Baldrige, Secretary

NATIONAL BUREAU OF STANDARDS. Ernest Ambler. Director

CONTENTS

Page

I. INTRODUCTION 1

II. FORMULATION OF THE PROBLEM 5

III. NUMERICAL RESULTS AND DISCUSSION 9

IV. CONCLUSIONS AND RECOMMENDATIONS 18

ACKNOWLEDGMENTS 20

REFERENCES 21

APPENDIX 1 22

APPENDIX 2 24

APPENDIX 3 Listing of the Program 27

LIST OF FIGURES

Figure 1. Schematic drawing of an oil shale retort 2

Figure 2. A spheroidal retort. The major axis of the prolate

spheroid is along the vertical z-axis. The transmitter

lies in the x-z plane with spherical coordinates

(r d

, 0 d’ 0) anc* the receiver coordinates are

(r , 6 0

, <j) 0

). The complex dielectric constants

inside and outside are denoted by e and e

2

, respectively 6

Figure 3.

Figure 4.

Figure 5.

Figure 6.

Figure 7.

Figure 8.

Figure 9.

Figure 10.

The geometry used in the computations: a = 45.7 m (150 ft),

b = 25.1 m (82.5 ft), r d

= 1.1 a, e d

= 90°; rQ= 1.05 a,

6 0

= 0° (10°) 360°, f = 4 MHz, k Q

a = 2WA = 3.84

|E

Q

| 2 vs scattering angle; dipole along x-axis

| E ^

| 2 vs scattering angle; dipole along x-axis

I^totl

2 vs scattering angle; dipole along x-axis

|E | 2 vs scattering angle; dipole along y-axis

<P

| E 1 2 vs scattering angle; dipole along z-axis

| E ^

| 2 vs scattering angle; dipole along z-axis

|£ | 2 vs scattering angle; dipole along z-axis

10

11

12

13

14

15

16

17

i i i

MODELLING OF OIL SHALE RETORTS

FOR ELECTROMAGNETIC SENSING TECHNIQUES

H. Chew*

We report here some work on the modelling of oil shale re- torts for electromagnetic sensing techniques. The aim is to

obtain useful information about the contents of the retort (e.g.,

rubble size, void ratio, etc.) by means of electromagnetic

probes. In this work, the retort is modelled by a spheroid with

an average dielectric constant which depends on the void ratio.

The near field due to a radiating dipole source in the vicinity of

a spheroidal retort is computed using the Extended Boundary

Condition Method due to Waterman, Barber, and Yeh. Numerical

results are given at 4 MHz for a retort with major axis 45.7 m (150 ft), minor axis 25.1 m (82.5 ft), bulk dielectric constant

8.8 + 3 . 7 j , and various void ratios. The results indicate

feasibility of determining the void ratio by remote

electromagnetic measurements. It is also believed that this work

may be of interest beyond the immediate context of oil shale

retort model ling.

Key words: oil shale retorts; remote sensing; scattering.

I. INTRODUCTION

In situ processing of oil shale offers many environmental advantages.

For example, the waste products largely remain underground and are not

released into the immediate environment. There are also many technical

problems connected with in situ processing, one of them being the gathering of

information about the contents and the state of the oil shale retorts. A

promising method for obtaining such information is electromagnetic remote

sensing. In this approach, transmitters and receivers are introduced to the

vicinity of the retort via boreholes (figure 1), and one attempts to extract

information about the contents of the retort by analyzing the received

signals. For this purpose, it is necessary to have a specific model which

relates the relevant physical quantities and allows the interpretation of the

signals.

The precise modelling of a retort of irregular shape containing rubble of

irregular size and shape is a difficult task both in principle and

numerically. To obtain tractable results, many simplifying assumptions are unavoidable. In this work, we model the retort by a spheroid embedded in an infinite medium of different electromagnetic properties (there is no

* On sabbatical leave ( 1980-1981) from Department of Physics, Clarkson College,

Potsdam, NY 13676.
Figure 1. Schematic drawing of an oil shale retort.

2
difficulty in treating media of different magnetic permeabilities, although in

the numerical work to be presented here both media are assumed to be non- magnetic), and compute the near field due to a radiating dipole source in the

vicinity of the retort using the Extended Boundary Condition Method (EBCM) of

Waterman [1], Barber and Yeh [2]. The effects of the boreholes are neglected,

and both the spheroid and the outside medium are assumed to be homogeneous and

isotropic. In practice, the shale rocks generally have some layered

structure, but the effects of anisotropy are probably not large, as laboratory

measurements at the National Bureau of Standards [3] show that the value of

the dielectric constant of shale rocks is essentially independent of

orientation. With these assumptions, the EBCM formalism is exact in the sense

that the resulting series is a solution to Maxwell's equations satisfying the

appropriate boundary conditions and therefore contains all the electromagnetic

effects (surface waves, body waves, etc.). Moreoever, the spheroid, with its

two geometric parameters (the major and minor axes) is sufficiently flexible

to simulate a variety of shapes and yet is such that the mathematical analysis

involved is manageable, if complicated. A drawback of this approach is the

complexity of the calculations, which are time consuming even when done on a fast electronic computer. The near field calculated, which includes both the

incident dipole field and the scattered field in the near zone, is a function

of the characteristics of the retort and of the surrounding medium. The

mixture of rubble and void inside the retort is described here by an average

dielectric constant which depends, among other things, on the void ratio

(defined as the ratio of the void volume to the total retort volume) in a way

to be discussed later. Thus, the dependence of the field on the average

dielectric constant may be used to extract information about the contents of

the retort. In this work, we are able to treat only the question of void

ratio and not that of rubble size. To gain some idea of this dependence

without elaborate formulas, we carried out a preliminary calculation (Appendix

1) for a spherical retort in the Rayleigh limit, and found that the scattered

field is a sensitive function of the average dielectric constant inside, being

roughly proportional to the difference between the dielectric constants inside

and outside. This strong dependence appears to persist in the much more involved spheroidal calculations as well.

The relation between the average dielectric constant and the void content

is a difficult subject, and a large number of workers [4,5,6] have examined

3
the problem of the effective dielectric constant of two-component systems.

For example, Sillars [7] obtained an expression for the dielectric constant of

a two-component system in terms of those of the constituents and the volume

ratios for the case when one component consists of spheroids of uniform size

embedded in the other medium. His result also depends on the ratio of the

spheroidal axes. Because the rubble is very unlikely to be spheroids of

uniform size, it is uncertain whether his result would be applicable to our

case, inasumch as it introduces an additional parameter (the ratio of the

axes). More complicated and frequency-dependent results are also available

[5]

, again under assumptions of doubtful applicability to oil shale retorts.

In this work, we shall use a simple empirical relation due to Lichtenecker

[6]

, wherein the logarithm of the dielectric constant is averaged in

proportion to the volume. If two media of dielectric constants e 1 and e

2 and

volumes and V 2 , respectively, form a composite medium whose average

dielectric constant is e (throughout this paper, dielectric constants refer to

dielectric constants relative to that of vacuum, except in eq (3)), then it

has been found empirically that in a large number of cases [6,8], e is given

to a good approximation by

Vi v

2 In e=-yln e 1 + —yin e 2 , (1)

where V = V x + V 2 is the total volume. In the case of an oil shale retort,

which consists of air (dielectric constant e

L = 1) and rubble (dielectric

constant e

2 = bulk value for shale rocks), the void ratio is V x /V and we may

rewrite (1) as

V - Vx

In e = — In e 2

= (1 - V

x /V) In e2

. (la)

This relation, which will be referred to as Lichtenecker 's formula, will

be used to relate the void ratio V x /V to the average dielectric constant of

the interior of the retort (rubble plus void), with the bulk dielectric

constant e

2 of the exterior medium being assumed known. In the case of

absorbing media, eq (la) is assumed to hold for both the real and imaginary

parts of the dielectric constant. This relation neglects the dependence of

the average dielectric constant on the size and shape of the rubble. This

neglect may not be very serious [6], as there is some evidence from the recent

work of Warne and Uhl [5], who concluded from some one-dimensional computer

simulation calculations that scattering effects depend largely on void

dimensions rather than rock sizes.

4
II. FORMULATION OF THE PROBLEM

In applying the Extended Boundary Condition Method [1,2], the scatterer

is replaced by a set of equivalent surface currents. The incident and

scattered fields are both expanded in series of vector spherical harmonics.

After a lengthy analysis [1,2], a transition matrix (T-matrix) is computed

which converts the known coefficients of the incident field into the

scattering coefficients. The elements of the T-matrix are surface integrals

of certain combinations of Bessel and Legendre functions which are computed

numerically. For the convenience of the reader, we summarize the key steps

invol ved.

Let the spheroid be centered at the origin with its axis of symmetry (zaxis) vertical (figure 2), and the oscillating dipole source with dipole

moment £ be located at coordinate r,. If the observer is at the

d

coordinate r = r, then for r>r

c) the incident dipole field may be expanded in

a series of vector spherical harmonics (if r<r

cj one needs only to interchange

the superscripts 1 and 3 in eqs (2) and (3)),

where k 2

t.(r) = l

D [a fi3(k,?) + b fi3(k,r)L

1 y

V V V

2 V V2

= wave number in medium 2

( 2 )

= e. (2 p+1

|

I (n-mjl! 4n (n+1 )

1 (n+m

)

I! =

rl (m=0)

i ^2 (m>0)

i - i

00

l

m=n

i

V e,o n=l m=0

a

v

jk

2

3

TT£ 2

P

•

J k

2

3 M 1 (k r ) , b = p

• ft 1

( k 2 f',). v

2 d v tt£ 2

v

2 d

(3)

The label v stands for three indices: m=0,l, • • «n; n=l,2,3««»; a=odd, even. Explicitly, we have

vx

cos md> ^m

[r . P sin m<}> n (cose) j n (k

2 r)]>

replaced by h

( 1) n

( k 2 r).

ramn

( k

2

P ) = TT VX M

z k ^nn

^ omn

a = 1,3.

(4a)

(4b)

(4c)

5
Figure 2. A spheroidal retort. The major axis of the prolate spheroid is

along the vertical z-axis. The transmitter lies in the x-z plane with spherical coordinates (r^, e^, 0) and the receiver

cordinates are (r Q

, e , <j> ). The complex dielectric constants

inside and outside are denoted by e and e 2 , respectively.

6
For a scatterer of finite conductivity, the EBCM requires as an intermediate step also the field inside the scatterer. Let this field be

£(r), with the expansions (k = wave number inside the scatterer),

t(r) - l [c

y (kf) + dv Ri (kr)], (5)

V

and

£se

(f) = l [p^ (k

2 ?) + q

v^ 3 (k

2 r)].

V V

( 6 )

Then the internal field coefficients c

y and dy may be found from the

linear equations

V { [K + J— J 3 c + [L + J—— I J dv } - - j a J

MV M

r yV J V

L yV y

p yV J J y (7a)

y { [I + / — L ] c + [J +/— K ] d } = - j b , (7b)

J yv y

p yv v yv yp yv

J

v

J y

v '

where e p

= e/e 2, y

p = y/u 2

> anc*

k 2 k 2

I = — /

n . [M3 (k r') x Mi (kr')] dS,

yV tt y

z V (8a)

k 2

2

- j = — /

n • [M3 (k r ' ) x fo (kr 1 )] dS,

yv tt y

z V (8b)

k 2

z K = — /

n - [N3 (k r ' ) x Mi (kr 1 )] dS,

y V tt (8c)

k 2 K

2 L = — /

n . [N3 (k r') x Ni (kr 1 )] dS.

y V tt

( 8d

)

The integrations are over the surface S of the spheroid (area element dS and

unit outward normal n ). After the internal field coefficients c^ and d^ are determined in terms of the known incident coefficients a^ and b^, the required

scattering coefficients p^ and q^ are given by

7
P = - j D

m y

y {[k* +

v

J ' ] c + [L 1 + / — I' ] d } (9a)

yV J V yV u

r yV V

x '

q,, - - j

/ / j 0 / 0

y {[I ' + /— L ' ] c + [J ' + / — K ' ] d } , L yV y yV V

L yV y yV V (9b)

where I' , J ' , K

1 and L' are obtained from the corresponding unprimed

yV yV yV yV

quantities by replacing the upper index 3 in the first vector spherical

harmonics in the surface integrals by the index 1.

The total field at the observer coordinate r = r is the sum of the

incident field i:.(r) and the scattered field £ (r). Because the receiver is

i sc not in general in the far zone, the total field must be calculated with the

exact M 3 and fa's rather than the asymptotic forms used in most scattering

calculations. In the near zone, the radial component of the field is in

general comparable to the e- and ^-components. While the total field can be

calculated from the series expansion in principle, it proved virtually

impossible to obtain convergence in general when the incident dipole field was included. The reason is that the terms in the dipole field expansion (see eqs

(3) and (4)) contains products of spherical Bessel and Neumann functions which

are usually such that when one factor becomes small, the other becomes large,

so that the product stays almost constant even when a large number of terms

are used. Accordingly, a separate subroutine was written to calculate the

dipole term separately from the closed-form expression [12], and the result

was added to the scattered field after the convergence of the latter was ascertai ned.

In the calculation of the scattered field, it is necessary to check

convergence with respect to three quantities: the number of m values, the

number of n values (NRANK), and the number of sections (NDPS) in the numerical

integration. Because the number of m values does not exceed n + 1, the first

check is strai ghtforward. Convergence with respect to the last two quantities

can be ascertained only by the slow process of increasing each until the

results converge at all angles. For 4 MHz and the physical parameters used

here, convergence was obtained for N^~ 12, NRANK ~ 32, and NDPS < 100. All

these will go up if either the frequency or the dielectric constants are

i ncreased.

8
III. NUMERICAL RESULTS AND DISCUSSION

In figures 4-10, we show the magnitude squares of the spherical

components of the total (dipole plus scattered) electric field and their sum

at the receiver (observer) coordinate (r

Q

, 0 Q

, 4 0

), where <j> o

= 0 . The

coordinate system is chosen so that the origin is at the center of the

spheroid and the z-axis along the symmetry axis. The dipole transmitter lies

in the x-z plane so that its spherical coordinates are (r^, e^, 0 ). This

choice of the x-z plane, which involves no loss of generality, simplifies the

incident field coefficients greately, because with 4^

= 0 , all the terms with

the factor sin m^ (see eqs (3) and (4)) are now zero while cos m 4^becomes

unity independent of m. The curves have been plotted for r Q

= 1.05 a, or 1.05

times the major axis, and r^ = 1.1 a, 0 ^

= 90°, while 0 Q runs in steps of 10°

from 0° to 360°. Results are shown for unit dipole strength and three

orthogonal dipole orientations: p along the x, y, z axes, respectively. For

all these orientations (but of course not in general), the curves are symmetric with respect to the x-axis (e Q

= 90° and 0 Q

= 270°), so only the

angular interval 90° < 0 q

< 270° is shown. Note that at 0 Q

= 90°, both the

transmitter and receiver are on the positive x-axis at a short distance 0.05 a apart, so the fields are very large in general. Because of the way the

coordinate axes are chosen, the electric field at the receiver has only 0 - and

r-components when the dipole source oscillates along the x- or z-axis, and

only a ^-component when the dipole oscillates along the y-axis. In the last

case, | E | 2 is simply equal to |E I 2 . L O L (p The other physical parameters are shown in figure 3. The axes a and b are taken to be the vertical and horizontal dimensions of an operating

Occidental Petroleum oil shale retort at Logan Wash, Colorado. In figures 4-

10, each figure shows four curves correspondi ng to void ratios 10%, 15%, 20%,

and 25%, respectively. The correspondi ng values of the average dielectric

constant e are calculated with the Lichtenecker formula from the bulk value

[3] 8.8 + 3.7 j

: 7.08 + 3.25 j (10% void), 6.35 + 3.0 j (15% void), 5.70 + 2.85 j (20% void), and 5.11 + 2.67 j (25% void).

It will be seen that the curves have generally more structure for 0 Q

< 180° (the transmitter side) than for 0 Q

> 180° (the far side). Much of this

structure arises from the complex interference between the scattered field and

the incident dipole field. Because the outside medium is highly lossy (Im

e

2 = 3.7), the dipole field is much attenuated on the far side, so the field

9
x

Figure 3. The geometry used in the computations: a = 45.7 m (150 ft), b

25.1 m (82.5 ft), r d = 1.1 a, 0 . = 90°, r = 1.05 a, 6 = 0°

(10°) 360°, f = 4 MHz, k Q

a = 2ira/A = 3.84.

10
ANGLE 0O

Figure 4. |E |2 V s scattering angle; dipole along x-axis 0

11
ANGLE G0

Figure 5. |

E

p j 2 vs scattering angle; dipole along x-axis

12
Figure 6. vs scattering angle; dipole along x-axi s .

13
ANGLE Q0

Figure 7. IE I 2 vs scattering angle; dipole along y-axis.

<l>

14
ANGLE Q0

Figure 8. |

E

Q

| 2 vs scattering angle; dipole along z-axis.

15
Figure 9. | E i 2 vs scattering angle; dipole along z-axis„

16
ANGLE 0O

Figure 10. 2 vs scattering angle; dipole along z-axis„

17
there is primarily the scattered field, which has travelled a much shorter

distance and therefore suffered less attenuation. It will also be noted that

on the far side, where the scattered field is dominant, the field intensity

increases with the void ratio. This is to be expected because a larger void

ratio means a greater difference between the dielectric constants inside and

outside, which in turn causes more scattering. This is also true for the

spherical case, as can be seen from the factor e

2 in eq (A2). For the cases shown, a 5% change in the void ratio changes the output intensity by

roughly a factor of two for most angles on the far side. There is thus hope

that the void ratio can be determined with some accuracy by measurement of the

output intensity, provided the bulk dielectric constant is known accurately.

This is also limited by the accuracy of the Lichtenecker formula, for which

more measurements with shale rocks would be desirable.

The numerical results presented here were obtained on the CDC Cyber 750

computer at NOAA, Boulder. The general description of the program and program

listing is given in appendices 2 and 3 respectively. For each dipole

orientation and fixed values of the dielectric constants at 4 MHz, it took

approximately 200 seconds to compute the squares of the components of the

fields over the whole angular range. For higher frequencies or larger

dielectric constants, both the number of terms needed (and therefore the

matrix rank in the solution of eqs (7) and (9)) and the number of integration

sections for the surfaces integrals (8) must be increased, thereby driving up

the computer time and storage capacity needed steeply.

IV. CONCLUSIONS AND RECOMMENDATIONS

The present work suggests the possibility of determining the void ratio

of an oil shale retort by remote electromagnetic measurements under suitable

conditions. When the shale rock is highly lossy, the source field is largely

damped out on the side of the retort away from the transmitter, so the field

there consists mainly of the scattered field. This field depends sensitively

on the void content of the retort and thus provides a useful tool for the

extraction of information on the contents of the retort. In practice, an average dielectric constant for the inside of the retort may be determined

using the observed signals and the void ratio obtained from the average

dielectric constant with the help of relations like Li chtenecker 1 s formula.

In the cases studied here, a small increase in the void ratio produces a

18
significant increase in the intensity of the scattered field. This

circumstance is very helpful in deducing the void ratio from the intensity

measurements. However, many assumptions and approximations are involved in

this approach, and limitations due to some of them are discussed below.

(1) Geometric approximations. After blasting, the sharp edges and corners of the retort are expected to be rounded off and the retort can probably be adequately modelled by a spheroid. Perturbation methods are available to handle small departures from spheroidal geometry, but additional

physical measurements will be needed to characterize the departures and

considerably more programming and computer time will be needed. We have not

been able to ascertain the effects of neglecting the boreholes.

(2) Limitations due to the approximate nature of the relation between

the dielectric constant and the void ratio. This is a sticky question and we have little to add to the existing literature. The Lichtenecker relation is

used here partly because of its simplicity, and it is in approximate agreement

with all the relevant data we are aware of. There is also another technical

advantage. In the Lichtenecker formula, eq (la), the void ratio V

x /V is to be

found from the measured values of e and e

2

. Because only the logarithms of

e and e

2 enter the equation, error in e and e

2 will tend to be suppressed,

giving a good determination of V^V. We believe more experiments with

accurate volume measurements on the relevant shale rocks at the operating

frequencies will be helpful.

(3) Possible presence of moisture, etc. The determination of the void

ratio from the measured signals in this model depends on the assumption that

the content of the retort consists of rubble and void. Because water has a large dielectric constant even at quite high frequencies, a small amount of

moisture may change the average dielectric constant significantly and

complicate the interpretation of the data. This is a serious potential

problem. The presence of pyrites would cause similar problems. Although

there is no difficulty in handling magnetic media in our approach, an unknown

amount of magnetic material would introduce additional uncertainty to the

analysis of the signals.

(4) Computer time. Even at 4 MHz, the amount of computer time needed to

generate the data presented here ran into hundreds of seconds for each void

ratio and dipole orientation. Mainly for this reason, we have not fully

explored the nature of the signal as a function of transmitter orientation, or

19
the effects of phase differences between the components of the dipole source,

etc. For general orientation of the dipole source, the field intensity will

not have the symmetry with respect to 90° and 270° in figures 4-10. It is

possible that certain orientations may be particularly favorable for the

determi nation of the void ratio. In the cases studied here, it appears that

the total intensity for 0 q

> 180° (figure 10) with a dipole source oscillating

vertically provides clearer separation between the curves corresponding to

different void ratios.

ACKNOWLEDGMENTS

The author would like to thank Mr. Doyle Ellerbruch and his group for

their hospitality at the NBS, Boulder, and for his kind permission for the use

of figure 1. He would also like to thank Profs. P. Barber and D.S. Wang for

making available to him important parts of the program and for invaluable

advice in the programming. The generous help of the N0AA computing staff and

Ms. J. Page's expert word processing are also gratefully acknowledged.

20
REFERENCES

[1] Waterman, P.C. Matrix methods in potential theory and

electromagnetic scattering. J. App. Phys. 50: 4550-4566; 1979. This

paper contains references to earlier work on the EBCM.

[2] Barber, P.j Yeh, C. Scattering of electromagnetic waves by

arbitrarily shaped dielectric bodies. App. Optics 14: 2864-2872;

1975.

[3] Jesch, R. Private communication.

[4] For a recent review see, for example: Lamb, W; Wood, D.M.; Ashcroft,

N.W. Long-wavelength electromagnetic propagation in heterogeneous

media. Phys. Rev. B21: 2248-2266; 1980.

[5] For a review of works more closely related to our problem see: Warne, L.K.; Ohl , J.E. Preliminary study of the remote sensing of a modified in situ oil shale retort by electromagnetic methods. Sandia

Report SAND79-1988; 1980 March. 92 p.

[6] Earlier literature, including Li chtenecker ' s work, is summarized in: von Hippel, A. Dielectrics and waves. New York: John Wiley and

Sons; 1954. 284 p.

[7] Sillars, R.W. The properties of a dielectric containing

semiconducting particles of various shapes. J. Instr. Elect. Eng.

(London) 80: 378-394; 1937.

[8] Bussey, H.; Roop, D. Private communication.

[9] Chew, H.; Kerker, M. ; Cooke, D. Electromagnetic scattering by a dielectric sphere in a diverging radiation field. Phys. Rev. 16A:

320-323; 1977.

[10] Edmonds, A.R. Angular momentum in quantum mechanics. New Jersey:

Princeton University Press; 1974. 146 p.

[11] Abramowitz, M. ; Stegun, I. A., eds. Handbook of mathematical

functions. Nat. Bur. Stand. (U.S.); 1965. 1046 p.

[12] Stratton, J.A. Electromagnetic theory. New York: McGraw-Hill Book

Co.; 1942. 615 p.

21
APPENDIX 1.

This appendix gives the derivation of the Rayleigh limit of the scattered

field for a spherical retort referred to in the Introduction.

For a spherical retort of radius a, analytical expressions can be found

for the near field when the frequency is sufficiently low so that ka = 2ira/x

(X = wavelength) is small compared to unity. They can be obtained from the

results of Chew, Kerker, and Cooke [9] generalized to absorbing media by

expanding the field coefficients in powers of ka and retaining only the lowest

order terms.

Let the retort be centered at the origin and the coordinate axes chosen

so that the dipole source (dipole moment £) lies outside the retort at

cartesian coordinates (0, 0, -d) (d>a). We expand the scattered field in a series of vector spherical harmonics in the notation of [9]:

oo m=£

£ sc

(?) - l l {— B i=l m=-£ n 2 2 w U,m) vx [h ^(k r)

am

(r) ]

+ h

(1) (k„r) 1 (r)}. (Al)

M 1 %

2 urn

1

Here, ^ (r) denotes a vector spherical harmonics as defined in Edmonds [10],

h^ v ; (x) denotes the spherical Hankel function of the first kind [11],

oj = 2irf, n 2 and k 2 being the complex index of refraction and wave number in

the (outside) medium 2. It is then shown in [9] that the expansion

coefficients are given by

B

E Ui rn) = -b

£

a d

£ (£ 1 m),

eM Uim) = -V

d

M(*i m )»

where a^ and

and the only

b„ are Mie scattering coefficients as defined in Stratton [12],

^ d d nonvanishing a

E

's and a^'s in our case are

a

d

E

(*. 1 ±l) = ±jk

2 3/„(-l)

£

(2t+ l)'

1/2

(P

x ±jP

y

)[)lh[]j(k

2 d) - (t+l)h[)J(k 2 d)],

a

d

E (2!±0) = jk 2 3(-l) £ P

z [4,e(t+l)(2t+l)] 1/2 h[1

) (k 2 d)/k 2 d.

22
a d

M Ui±l) = jk 2

3 (-l) t ["(2t+l)] 1/2

(P x±jP

y

)

Making the small ka expansion indicated above, we find for the leading

terms in the expansion

n ^

a

3 f — p ± + 2

2 i 2 m m

Kc {r) '

3

• — • 7-^7 'l [h

a ’( k 2

d > ' 2h

o V)]

(rh ( P(k 2 r))

2?

((P coS(j) - P sin(j))cosee - (P sincj) + P cos

x y x y

h

( J} (k 2 r) - h ^ i

^ ( k 2 d

)

(P C0S<(> - P sin<j>)sin0 r] + 3P

J z ^2

(rh^ ) (k 2 r))‘ 2h^ (k 2 r)cose

sinee - r )}> (A2)

where (r, 0, <f> ) denote the spherical coordinates of the receiver, the prime

indicates derivative with respect to r, e 19 e

2 denote the complex dielectric

constants inside and outside respectively, and n 2 denotes the complex

refractive index of the external medium. For an oil shale retort, where the

void ratio is typically of the order of 15 percent, e 1 does not differ very

much from e

2 and the factor e

l - e2 is a sensitive function of the void

ratio. If the medium outside is highly absorbing, then on the side of the

retort away from the source the total field is dominated by the scattered

field £ . Thus, measurment of the field there will provide valuable

sc information about the void ratio.

23
APPENDIX 2.

This appendix contains a general description of the program and

instructions for its use. The program computes the incident field coefficients, the elements of the

T-matrix, and from these the scattering coefficients. Two subroutines, GENBKR

and GENLGP, generate the necessary spherical Bessel functions of complex

arguments and Legendre functions. Subroutine GAUSS carries out the numerical

integration over the spheroid surface to give the needed matrix elements.

After the scattering coefficients are calculated, they are multiplied by the

appropriate Bessel and Legendre functions of the receiver coordinates to give

the spherical components E sc sc and E sc of the scattered field. The

dipole source has been taken, without loss of generality, to be in the x-z plane with spherical coordinates (r^, e^, ^

= 0). The scattered field is

calculated in the vertical plane <j>o = constant (taken to be zero in the

program listed in Appendix 3) for a fixed radial coordinate ro (taken to be

1.1a in Appendix 3), while the angle 6o is varied from 0° in steps of A0

(denoted by DLTANG and with the value 10° in Appendix 3) up to 360°. At each

angle subroutine DIPOLE computes the dipole field due to the same source and

adds it to the scattered field to form the spherical components of the total near field Ej'

0

*', Eg

0t

:

i E

t ot

| 2> | E t«t|

2> and

U I

tot

£ |2 = tot

1

Finally, the program computes and prints IE ..tot . _ , _tot , „ , „tot

E

2

+

tot

E t0t

|2 + 2 as functions of

the scattering angle e Q (denoted by SCANG in Appendix 3). The geometry is

shown in figure 3.

Following is a list of the input quantities needed to obtain numerical

results:

1.

The real and imaginary parts of the dielectric constants inside

and outside the retort, denoted by DCNR, DCNI; DCNR2, DCNI2,

respectively.

2. The frequency. This is expressed in terms of the vacuum wave number WN = k Q

= 2-rrf /c in the input.

3. The size parameter k Q

a = 27Tfa/c = CONK (a = major axis of the

spheroid, which is taken to be prolate in Appendix 3), and the

ratio of the spheroidal axes a/b (denoted by AOVRB).

4. The radial coordinates of the dipole transmitter (r^) and of the

receiver (e Q

).

24
a. In subroutine DIPOLE, ( =RD ) is entered as (r

d /a ) • (

k

Q

a /k

Q

)

,

or RD = (r d /a )*C0NK/WN, r d /a being taken to be 1.05 in

Appendix 3. The coordinate r Q (denoted by R0) is entered as

(r

Q /a)*a = (rQ /a ) • (k Q a/k

0

), or RO = (

r

Q /a )*C0NK/WN, with r Q /a

being set equal to 1.1 in the program. Both rQ /a and r

d /a

may be set at any value > 1. However, if r Q

< r d

, it will be

necessary to make the interchange of Bessel and Hankel

functions indicated in section III.

b. In subroutine AUDPRC, they are entered through the arguments

of the radial functions p '(eRH0P) = (r d /a)*k 2 a = (r d /a)

k Q

a /e2= (r

d /a) *C0NK*CSQRT (DCN2) , and p (eRH0) =

(r

Q /a) • k 2

a = (r Q /a) *C0NK*CSQRT(DCN2)

.

5. The angular coordinates of the transmitter e d

= THETAD, <j> d being

set equal to zero, and the azimuthal angle of the receiver <(> e PH. They are given the values THETAD = 90° and PH = 0 in

Appendix 3.

6. The cartesian components of the dipole moment P x

= PX, P^, = PY,

P

z = PZ. These are allowed to be complex so that phases between

the components may be introduced. They have to be entered in

both subroutines ADDPRC and DIPOLE.

7. The number of values of m (Nm) and n (NRANK).

8. The number of sections used in the integration over the spheroid

surface (NDPS).

9. The angular increment (denoted by DLTANG) in the scattering angle

(denoted by SCANG), and the number of angles for which the

squares of the field components are calculated. The last

quantity is denoted by NUANG and is equal to 360°/DLTANG + 1.

Most of these input quantities are entered on four data cards at the end

of the program (see listing. Appendix 3).

The first card lists Nm, NRANK, and three numbers (1,8,1) which are to be

left alone.

The second card lists CONK, AOVRB , WN, DCNR , DCNI.

The third card gives NDPS.

The fourth card lists THETAD, PH, DCNR2 , DCNI2, DLTANG, NUANG.

25
The data on these cards must be entered so that the last digit of the

first entry is on the 12th space, the second on the 24th space, the third on the 36th space, etc. These data are read by subroutine RDDATA and stored in

various common blocks for use in the other subroutines.

26
APPENDIX 3. Listing of the Program

PPPGR4M SHALE ( INPUT, OUTPUT, TAPE 5 = INPUT, TAPE 6 * OUTPUT)

r OM Pi cy A,B,nCM»S0,0S,PH,PB(40),HRK,RRK.»HFPS»BEPS,S*$QP,QSR,BKK,

1HTBKX , RTBKK , S Al, cri, SQP? ,D,CCC, RSSL SP ( 41 ),CNEUMN( 41) » PBESSLI40)

,

?SQ?»QS?,DCM?,CI»RP,HH,CCKP,CIKP,DCKK,BSLC M P,CNFUM,ACANSR

C PERFORM TUP NUMERICAL INTEGRATION AND FILL THE A AND R MATRICES.

rnMMPM OTP , PTn , r P T COMMON /MTXCOM/ NR ANK, NR ANKI , A ( 80, 80 , B ( 80, 80) ,C M XNRM( 80)

common /en^C dm/ PNVL [ G( Al) . PSLCMP ( A1 ) , CNFUM( A1 ) , BSLKPR ( 41 ), bslkpi

1(0)

COMMON /r MVC n M/ MM, K M\/,r MI ( AO ) , c MV, CM?, TWM, proom

in MM ON /ROVCOM/ DCNP,0CNI,CKPRR,CKPRI,CKR,0CKR,C0NK,A0VR8,WN,IB

Common /TUT COM / THFTA,STNTH,CnSTH,C0H(6),F°PSf6),NSFCT,NDPS(fc)

1, TMFTAO

,

0 H , KC FCt

COMMON /UVC com/ AC ANS (361 , ?? ) , AC ANSR( 361 ), DLTANG,DCNP2, DCN I2,NUANG

DIMENSION n L

p M T X ( ? 5 6 C 0 ) » r

L R T 0 T (1444) ,RH(40),WT (300),ASC(300)

F OUTVAL pnc F (A(l,l),CLRMTX(l)),(ACANS(i,i,i),CLPTOT(l))

C set program constants.

Cl = ( 0,0, 1 .0)

DTP = .017^532QP519QA3

PTO = 57. P9577Q5131

CPI = p .1^159?6535P9R

C CALL ROUTINE TP dead data and print hfadtngs FOR OUTPUT

20 CALL p D D AT A C CLFAP THF ACCUMULATING ANSWER REGISTER (USED IN ADDPPC).

on AD J=1,1444

ClPTpT(J) = 0.0

40 CONTINUE

DO 41 J = 1 , 3 61

ACANSP(J) = 0.0

41 CONTINUE

DCM? = CMPLX (nCNR?,DCNI2)

so? = C SORT ( OCN?

)

0 S 3 = 1 .D/SO?

dcn = c m plv(dcnp,dcnt)

SOP = r SORT f ncN/DCN? ) OSP = 1.0/SOP

SOP? - SQP*SOP

so = ^S OPT ( DCN

)

0* = 1,0/SO

c set multiplier p 8 Q DFPENDenjt ON IP VALUE (SVMMFTPV INDICATOR).

RP9 = 1.0

IF(IP.FQ.A) R 8 9 = ? , 0

rdyect = 1,0

C SET Ud a LOOP FOP EACH M VALUE.

DP QOO I M = 1 , mm

C SFT V nFPFNDFNT VARIABLES.

27
C M V = CMT( IM)

K M v s r mv CM? = c mv* 3

)1 ? prfidm = l.o

IF(kmv.gt.O) go in 4 4

FM * 1.0

GO Tn 6

C

44 FM z ?, o OIJANM * CMV

on 5? I FCT = 1,KMV

0 U A N M z 0'.JANM+1,0

pPnOM = OUANM*PPOOM/?,0

52 CONTINUE

60 OEM a ? .O/FM

TWM = ?.0*CMV

C TNTTJALI7F ALL MATRTX AREAS TO 7 F R P

n n R A I = i,?5frC0

CL pmty ( T ) = 0.0

RQ CONTINUE

€ S F T U D A L n 0P FAR ALL VALUES OF THETA.

c SFT IIP GFNEPAL l OOP FOP CORRFCT NUMRER OF I NTF GR A TT ON SECTIONS.

AO ROO TSFCT = 1,NSFCT

k R F C T = T S P C T NTHFTA = NOPS(ISFCT)

IF (TSFCT. FO.l) PALL GAUSS(WT # ASC*NTHETAfO.OfFPPS(ISECT))

TF(TSFCT.NF.l) CALL GAUSS ( WT . ASCjNTHETA* EPPS ( ISFCT-1

)

t EPPS ( I $ EC T )

)

C FNTFD THFTA LOOP F P o EACH SECTION,

0° 7 0 A ITMT A = 1, NTHFTA

THFTA = ASC(ITHTA)

CPSTH = COS(THFTA)

RINTH = STM (THFTA)

c gfmfratf thf lfgenorf polynomials.

CALL G F M L G P c EVALUATE kp and ITS DEPTVATTVE AS a FUNCTION nF THFTA.

3 4 R CALL GFNKP

c GFNFPATF ail NFCESSARY RFSRFL and NEUMANN FUNCTIONS AND THEIR PATIOS.

CKPRD = REAL ( SQ?*CKR)

AKPRT a AIM AG(SQ2*CKP)

CALL G E N 3 K

R

A C K P = CMPL X (AKPRR ,CKPPI

)

CTKP z 1 . 0 / C C K

P

DC*K * RQ?*DCKP

DO 340 I*1,NRANKI

RSRLSPf T ) = R S L C

M ° ( I

)

CNFUMN(I) = CNFUN(I)

3 4Q CPNTINI1F

CKPRR » PFALfSQ*CKR)

CKPPT a ATMAG (SO+CKR

)

PALI. GENRKR

I F ( T THT A ,NF. NTHFTA) GO TO 7Q

70 conttnuf

DO 3R0 K * l.NPANK

28
r»

o

r>

o

n

o

r>

R R E $ S L ( K ) = R S S L S P (K ) /PSSLSP(K+1

)

PH(K1 = (R5SL^P(K)+CT*0NFUm N(K) ) /

( R EEL* P ( K + l ) +C I*CNEUMN ( K + l ) ) pp (K)=$OP*CMPLX (RSLKPP(K) » R S L K P I (K) ) /CMPLX(RSLKPP(K.+ 1)*BSLKPI ( K 1

)

1)

3^0 C0NTTNUF

CP^SIM s wt ( ITMTa ) ** tnth

c ^ F T IJD A !_nr)D ppp EACH ppw OF THE MATRICES.

CPHW = 0.0

CROW* * CMV

n P *no T p OW = 1,K<PANK

TPnui * T d 0 W + N P A N K c p o w = c p n w +

i

. o CPnv/M » CPOWM + 1.0

r onwi =CP0W+1 .0

pp = pr ( T P OW ) HH = PM ( TPOW

)

C EFT hd A L 0 fl D Fop faCH COLUMN OF THE MATPICES.

Cmt * 0.0

C C OL m = C

M V

no aoo tool = 1 , N R A N

K

IC0L1 = TCOL+NRANK

CCOL * CCOL+1.0

ccolm = ccni.M + i,o

cc on = croi + i.o

c CALCULATF FPFCUFNTl Y USED VARIABLE COMBINATIONS.

CPIJ = CPOW+CCPL

CPEETJ = CPOW*CC0L

C M C P C 0 = CM? + OEM*CP EET J*COETH**?

DNPnn 0 = PNMLLG(TROW)*PNMLLG(ICOL ) PNPOCl = p N m l L C ( I P 0 W ) * p N M L t G ( T C 0 L + 1

)

PM PICO = P N M l l G ( I P 0 W + 1 )*PNMLLG( ICOL)

P

N

0 1 C 1 = PNMLLG-( IPnw + 1 )*PNMLLG ( ICOl +1

)

R 1 A = CPnW*C r'ETHtPNPlCl“CPnWM=4'PNP0Cl

RIB r CC0l.*C0STH*PNQ 1C1-CCOLN*PNP1CO

R k K = P R ( ICOL ) HRK=(RGELGP(IP0W+1) +CT+CNFUMN ( IPOW+1) )*CMPLY(RSLKPR(IC0L+1)»B$

ILK PT ( TC PL +1 ) ) rrk x RECLEP(IPPW + n*CM D LV(PELKPPfICnL + l)*RSLKDi(icOL+l})

mc PS = q^p^hPK

P F PS x qfp^bpk

IF(IR.EO.Q) GO TO 380

c IF TR = R (MIPPOP symmftpv R 00 y ) , I =L *0 IF IPOW AND ICOL ARE BOTH

c 000 OP pnTM FV/FN, J=K=0 IF IPOW AND ICOL APF 000, EVFN OR EVEN, ODD.

TF( ( IPOW-nCOL ) . FO. ( { IRPW+ICOL ) /? )*?) GO TO 393

C TFFT F n p MjQ (IF M«0 TMF i ano L EURMATRTCFS ARF ZERO).

380 T E (KMV. FQ. 0 ) GO TO 3 00

CALCULATE THF x , L

,

I , J ANn k,L»T,J (PRI M E) MATPTY ELEMENTS AND PLACE

THFM TN THF a AND R MATRICES RFS D ECTTVELY SO AS TO FORM A-TRANSVERSE

ANn r-tpam^vfp^f MATRICES.

CALCULATE THF TFPM epp THF CURRENT ELEMENT IN TMF I MATRIX,

29
r>

o

r>

r>

o

o

r>

o

o

B 1 R 1 A + B 1 B MTRKK = HM*8XK

RTBKK = RP*PKK

SA1 * PNPlCl*(CROW*CROVl*BKK+CCOL*CCOll*HH-CRSSIJ*(CRIJ+2.0)*CIKP)

i^nrK^^s imth

S = SAl.+ f CC* D *( 1 .C + HTRKK )-CC0L*HH-CPnw*8KK4CPSSIJ*CIKP )*R1*CCKP

A ( Tent, IP0W1 ) =A f TOOL, TR0W1 ) +B89*CMV*S

P

M $ IN*S*HBK

S PI * PNP1C1* ( CPnW*CRnwi*BKK + CCOt*CCnLl*BP-CR$S IJ* ( CP IJ+2.0) *C IK P)

l+OC^^^STNTH

S=(COKP*(1.0+BTRKK)-CCr)L*RR-CRnW*BKK+CRSSIJ*CIKP)*Pl*CCKP+SBl

R ( TCni# TPHWl ) =P f irnL , TPQWl ) RP9*rMV*SPNSIN*S*BBK

CALCULATE T HP TFRM FOP THF CURRENT ELEMENT IN THF L MATRIX.

5 = (CPKP*(S0P? + HTBKK) -CCOL*HH-CPOW*B KK + C p$S I J*C IRP)*bi*CCKP+SA1

a ( TCHL I » 15 nw ) =A ( icon. IROW ) -B89*CMV*SPMSIN*S*HEPS

c = ( cc k P* ( PQRP + PTPKK )-CC0L*RP-CP0W*BKK+CPSST J *C I K P )*R1*CCKP+$B1

B (TCHLlflROW) *B(ICOLl* IPnw) -B89*CMV*$RMSIN*S* BEPS

3^0 TF(TP.PQ.B) GO tq ADO

CALCULATE THF TFPM POP THF CIJPPFNT ELEMFNT in THF J MATPIX.

3 9? Al? = CMCPCn«PNPiri-OFM«(CPnw*CCniM«cnSTH*PNPlCO + CCOL«CROWM^COSTH 1»

[ PN

lROCl-CRnvM*CCOLM*PNROCO>

B 1 A * CCnL*CCnil*BlA

RIB = C PPW*CP CW1*R IB

n = ORM+ncKK

CCC = SOP?*HH

S=(CrKP*fRKK-CCC)+SQP?*CPnw-CCOL)*Al?«CCKP+(BlA-SOR?«BlB)«GINTH*D

A(TCOll,IPnwi) »A(TC0L1*IR0W1) +aBQ*FRM$TN*S*HFPS

CCC = RP*PQR?

S*(CCKP*(BKK-CCC)+S0R3*CR0W-CCnL)*A12*CCKP-(S0R2*BlB-BlA)*SINTH*D

P•( ICDL1 TR0W1 ) =R ( TCHL 1 . TR0W1 ) + B 8 9*S P MS I N*S * RF PF

CALCULATE the T F P n FPO THF CURRENT ELEMFNT IN THF K MATRIX.

81 = (R1A-B1P)*SINTH

S * (CCKP*< 8KK-HH) +CPnW-CCOL)*A12*CCKP+Bl*0

A ( ICUL * TP OW ) = A ( T COL • T PO W ) +B89*SR MS T N*S*HBK

p = (CfKP« (RKK-PP )+CPOW-CCnt )*A1?*CCKP + Rl*n

o ( TCOL , TPPW ) = R ( ICOL > TPOW) +r89*$PMSIN*S*RRK

AGO CONTINUE

C calculate thf normalization FACTOR (USED IN AO0PPC).

Cppow * ip nw IF(KMV.GT.O) GO to

FCTKI = 1.0

GO TO 440

4?fc IP (IROW .GF.KMV) GO TO *R0

C M YMP M ( IP PW ) = 1.0

GO TH 6 CO

430 T B F C T = JPOW-KMV + 1 TFPCT * TPOW+KMV

30
FPPnn * trfct

FOTKT = 1.0

nn A3? LFCT = TBFCT.IEFCT

FCTKT = FCTKT*FPPPD

p ppoo = p dr nn + i ,o

A3? CONTINUE

A AO CMXNPM(IPnW) = 4.0*CKP0W* (CKROW+1. 0) *FCTKI/ ( EM*(2 .0+CKR0W+1.0))

600 CONTINUF

7 00 CONTTMIJF

800 rnMTIMUF

C ppccpps C 0 N P M T F 0 MATRICES

CALL PPCSS M 900 continue

O p T 0 20

FNO

SUP POUT TK'P GAUSS ( WT» ASC»N, A A, BB ) niMFUS TPM WT( N)

,

ASC ( N

)

HOUBLF PRECISION PT, CONST

0 ata oj .CONST, T PL / 3. 1415926535897900, .14067881635700, 1. E-l?/

OATA C1,C?, C3,C4/ . 125, -.08 07291666, .246C 2864 58,-1 .824438767/

I F(N ,NF .1 ) PP TO i

A S C ( 1 ) = 0.5773502692

W T ( 1 ) = 1.0

P p T U p N

1 ON = N

NO I V 2 * N/7

N P 1 * N + l N N P 1

= N * N P 1

A P P F P T = 1. /O^OPT( (N + 0.5)**2 + CONST)

CPN1 * 0.5*(BB-AA)

0 ON? = 0.5*(BB + AA)

on ioo k = i , not v?

R = (K-.?S)*PI

B T S 0 = 1. / ( B*B)

PPPPHT = R*(l .+RISQ*(C1+RTSQ*(C2+BISG*(C3+C4*BISQ) ) )

)

X = CnS(APPP0T*PF9PPT)

113 p M 2 = 1.

DM1 = X OP 110 IN = ? * N

P = ((2*tn-1)*x*pmi-(in-h*pm2)/IN

DM? = o M 1 110 P M 1

X P

P M 1 - p M ? Al'X = 1 « / (l.-XTX)

0 F P 1 D = QM* ( DM1_X*P ) *AIJX

0 f p d d * (?.*X*0FP1P-NNP1*D)*AIJX

P ATTO x P / 0 F P 1 D XI = y-patio*(1.+PATTQ*OFP?P/(2.*DER1P))

TF( A«S( XT-x)-TOL) 111,111,11?

11? X = XT

GO Tn 113

111 A F ^ f K ) = -X

W T ( K

}

x ?.*(l.-XTX)/(0N*PNl)Tt2

31
non

oooo

o

o

n

r>

r>o

ASCfNPX-K) « - A $C ( K ) 100 UT(MPl-K) = WT(K)

I F ( MOf) ( N , 3 ) »F0*0) OP T F H4

ASCfNDXVP+l) * 0.0

N m 1 = N-l

MM? * n-2

prod = dn

DO 1?0 K * 1,NM?,2

120 P*? n n = FL0AT( NM1-K ) /FLOAT (N-K )*PR0D

WT ( NO TV 2 + 1 ) * 2./PR0D**?

1 1 A DO 130 K = 3,N

A^C(K) * CONI* ASC ( K ) 4C0N?

130 WT(K) * C0M1*WT(K)

return

FND

SUBROUTINE podata

C A ROUTINE TO READ INPUT F n P THE PROGRAM.

COMPLEX ACANSP

COMMON /MTXCOK/ NPANK,NRANKI

COMMON DTP , RTH, CPI

COMMON /CMVCOM/ NM,KMV,CNI ( AO)*CMV,CM?,TWM,PRODM

common /THTCOM/ THETA,SINTH,CnSTH,C0H(6),F p PS(6),NSECT,NOPS(6)

l, thftad,ph,ksfct

common /bdycom/ dcnp,dcni,ckprr,ckpri,ckp,dckr,conk,aovrb,wn,ib

COMMON / UVCC0M/ACANS(361> 2, 2 ) , AC AN S R ( 36 1 ) , DLT ANG * D CNR ?, D CN I 2, NU ANG

COMMON /OUTCOM/ I OUT

DIMENSION F p D F G ( 6

)

READ necessary INPUT DATA.

C ARD1 NM = NUMBER OF M VALUES, NRANK = N VALUE (MATRIX ORDER),

N^ECT = NUMBEP OF sections IN THE BODY, IB = SYMMETRY CODE IB * 8

FOR MJPROP SYMMFTPY ABOUT THETA = 90 DEGREES, IB * 9 FOR GENERAL

PEAD(B,BO) nm.npank,nsect,tb,iout

TF (FOF (5) .NE.O) GO TO 190

NPANKI = NPANK+1

W P I T f ( 6,89)

WPTTR(6,9?) NM,MRANK,NSECT, IB

CARD ? — CONK » KA OF RODY,A0VRB = A/B, PATIO OF MAJOR TO MINOR

A X T S, WN = VACUUM WAVF NUMBER USED IN ADDPPC,DCNR » REAL PART OF

DIELECTRIC CONSTANT INSIDE, DCNI = IMAGINARY PART OF SAME.

PFAD( 5, 96 ) CPNK,AOVRB,WN,DCNR,DCNI

IF (FOF(5).NF.0) GO TO 190

WRITE (6, 100)

WPITF (6,104) CONK, AOVPB,WN,DCNR, DCNI

HO 5 1=1,40

CMim * FLOAT(I-l)

5 CONTINUE

TFfNM.EO.l) CMI(l) = 1.0

CARD 3 —- NDPS = NUMBER OF INTEGRATION DIVISIONS FOP EACH SFCTICN

HE THE BODY (MUST BE A MULTIPLE OF 4).

RFADC 5, BO ) (NDPS ( I ) » I-1,NSECT ) IF (EOF (5) .NE.O) GO TO 190

WeiTE(6,l?0) (NDP3(I),I»1,NSECT)

32
o

o

a

o

o

c r A 0 D 4

THETAO = THFTA r F DIPOLE.

PM * AZIMUTHAL ANGLF OF OBSERVATION PLANE,

DCNP? + C I *DC NI 2 = COMPLEX DIFLFCTPIC CONSTANT OUTSIDE.

Dl.TANG s INCREMENT OF SCATTFR TNG ANGLF IN SCATTERING PLANE IN DEGREE

NUANG = NO. hf SFCTTHN^ IN SCATTERING PLANE WITH 360 OFGRFFS

RFAD(5,Q6) THFTAD, PH, DC NR 3, DCNI2, DLTANG, NUANG

IP (PDF (5) .NP.O) GO TO 190

WPTTF (6,117)

WPJTF (6*104 ) THPTAD,PH,DCNR?,DCNI2

C COMPUTF PND POINTS FOR THFTA.

CALL CALENP

DO 140 I * ltNSFCT

FPDEG(I) = PTD*EPPS(I)

140 CONTINUE

WRITE (6,148) (EPDFCm,I = l,NSECT)

R FTIJPN

190 WRITE (6,201)

STOP

80 FORMAT ( 8112)

88 FOPMAT( 1H144X, 60H CASES MATRIX RANK

1 BODY SHAPF)

97 FORMAT ( 1H044X,4H5 ) 96 FORMAT ( 5E 1 2 • 5, 112

)

117 PORMATf 1H025X, 79HVARI0US PARAMETERS THFTAD

1 DCNR? DCNI 2

)

100 FORM AT( 1H029X , 76HB0DY PAPAMFTFRS K ( A

)

1 WN DIE LECTR IC1 ) 104 FORMAT( 1H044X, 5P15, 3)

120 FORMAT ( 24H0 I NT E GR A TI ON S / S E C T 1 0N8 1 1 2 , / < 1H02 3X , 8 1 1 2 )

)

148 F0RMAT(24H0 END P0INTS8F12.4, /( 1H023X,8F12.4) ) 201 FORMAT( 1H0, 23H***** fnd OF DATA 4****)

END

SUB POUT INF GFNLGP

C A ROUTINE TO GENFRATF LEGENDRE POLYNOMIALS.

C THE INDEX ON THF FUNCTION IS INCREMENTED BY ONE.

common / mtxcom/ nrank,nranki

COMMON DTR,RTD,CPI

COMMON / FN C COM / PNMLLG(41)

COMMON /CMVCOM/ NM, KMV,CMI ( 40) ,CMV, CM2, TWM, PRODM

COMMON /THTCOM/ THFTA, S I NTH, C OSTH, C DH ( 6 ) , F P PS ( 6 ) , N S EC T , NO PS ( 6

)

DTWM * TWM+1.0

I F(THFT A) 16, 4, 16

4 IF(KMV-1)6, 12,6

6 DO 8 ILG - 1, NRANK

I

PNMLLG ( ILG) * 0.0

8 CONTINUE

GO TO 88

12 PNMLLG ( 1) = 0.0

PNMLLG(7) « 1.0

PLA - 1.0

GO TO 48

SECTIONS

PH

AOVRB

33
r>

r>

o

16 IF(KMV) 20,20,40

C THE SPECIAL CASE WHEN M » 0.

20 PLA = 1.0/SINTH

P L R » C0STH*PLA

PNMLLG(l) 3 PLA

PNMLLG ( 2) = PL B

T B EG = 3

GO Tn 60

C GENERAL CASE FOP M NOT EQUAL TO 0.

40 00 44 T LG = 1,K*V

PNMLLG* ILG) =0.0

44 CONTINUE

PLA = p RODM*S INTH**( KMV-l

)

pnmllG( KMV+1) = PLA

48 P L B = r)TWM*COSTH*PL A PNMLLGL KMV+a) = PLB

IBEG = KMV + 3

C 00 PEGU PSION FORMULA EOP ALL REMAINING LEGENDRE POLYNOMIALS.

60 CNMUL = IBFG+IBEG-3

C N M = 2.0

CNMM a DTWM

00 80 ILGR a IBEG.NRANKI

P LC = (CNMUL«rOSTH*PLB-CN M M«PL A ) /CNM

PNMLLG(ILGP) = P LC

PLA * PLB

P L 8 = PLC

CNMUL = CNMtJL+2.0

CNM = cnm+1,0

CM mm r CNMM+1.0

80 CONTINUE

88 RETURN

END

SUBROUTINE GENBKR

C GENERATE BEESEL FUNCTIONS FOP COMPLEX ARGUMENTS.

C THE INDEX ON THE FUNCTION IS INCRFMENTFD BY ONE.

COMPLEX CKPR,RJ(410)»A,PX,P$LCMP,CNEUM,CKP2

COMMON DTP » P TD » C P T COMMON /MJXCOM/ NR ANK, NR ANIX I COmmun /FNCCOM/ PNMLLG(41 ) ,BSLCMP(41 ),CNEUM(41 ) , BS LK P P ( 4 1 ) , RS LK P I K 41)

COMMON /BDYCHM/ dCNR,0CNI,CKPPP,CKPPI,CKR,DCKP,CONK, aovrb,wn,tr

common /THTCOM/ THETA, SINTH, C0STH,CDH(6), EPPS ( 6)> NSEC T>N DPS ( 6)

DIMENSION KM(410) ,L( 41)

CKPP a CMRLXf CKPRP,CKPRI ) PM = CAB? (CKPP ) DO 45 NRFS = 1,5

NVAL = 2*NR ANK I*N3 FS

GENERATE BESSEL FUNCTIONS.

PJ(NVAL+1) » (0.0, 0.0)

PJ(NVAL) = (1.0, 0.0)

34
oooooo

ooo

TE(RM.GT.?.0) RJ(NVAL) = ( 1. 0E-10* 0.0)

I FCRM.GT. 10.0) RJ(NVAL) * ( 1 • 0E-20> 0.0)

IE (RM.GT.25.0) RJ(NVAL) = ( 1 . 0E-30 * 0 . 0

)

TF*NVAL+2

A = C*IN< CK PP ) /CKPP

* * 0

do 10 I=?*NV&L

T J = T F — T RJ(TJ-l) = RJ ( IJ ) * FL OAT ( 2* I J— 1 ) /CKPP-RJ ( I J + l

)

IF(CABS(RJ(IJ-1)).GT.1.0E10) GO TO 8

GO TO Q

8 K = K +

1

PJ(IJ-l) = RJ (I J-1)*1.0E-10

PJ(TJ) = PJ ( TJ )*1 .0E-10

KM C I J ) * K 9 KM ( I J-l ) * K 10 CONTINUE

py = A/ p J ( 1

)

L ( 1) = 0

DO 1*5 J * ?,NRANKI

L ( J ) = L ( J-l)

I F ( KM f J ) .NF »KM(J— 1) ) L(J) = L ( J-l ) +1

1*5 continue

nn ?o I=1 .npanki

bslcmp(I) * p j

(

n*py*io.o**(-L (nno)

RSLKPR(I) = REAL(BSLCMP(I )

)

BSLKPI(I) = AIMAG(BStCMP( I)

)

20 CONTINUE

GFNFPATF NEUMANN FUNCTIONS FOP TEST.

CNEUM(l) * -CCOS (CK°P ) /CKPR

CNEUM ( 2. ) * C N F U M ( 1 ) / CKPR—

A

DO 30 I=3*nrankI

CNFUM(I) = CNFUM( 1-1 )*FL0AT(?*I-3)/CKPR-CNEUM(I-2

)

30 CONTINUE

PERFORM TWO TF STS ON RFSSFL AND NEUMANN FUNCTIONS. FIRST TEST IS

MOST ACCURATE FOR LARGE ARGUMENTS AND THE SECOND IS MOST ACCURATE

FOP SMALLER ARGUMENTS. IF EITHER TEST IS PASSFD, FUNCTIONS ARE GOOD.

FOP LARGE ARGUMENTS APS(BESSEL) SHOULD EQUAL ABS (NEUMANN)

.

C = 1.0E-05

QUART = CARS(PSLCMp(l) )/CA8S(CNEUM(1) )-1.0

OUANT > CABS( BSLCMP(NRANKI) ) /CA8S(CNEUM(NRANKI ) )-l,0

IF ( (ABS (QUART) .GT.C) .OP. (ABS(OUANT) .GT.C) ) GO TO 32

RETURN

c PERFORM WPONSKIAN Tf S T IF LARGE ARGUMENT TEST FAILS.

32 CKP? = CKPR**?

C BFSSEL TEST

OUANR T « CABS(CKR?*(PSLCMP(3)*CNFUM(1)-PSLCMP(1)*CNEUM(2))-1.0)

C N F U M A N N TEST

35
or>r>

o

r»

or»r»o

OUAMNT = CABS(CKR2*( BSLCMP(NRANKI )*CNEUN(NRANK )-BSLCMR(NRANK)*CNEU

IMfMPANK I) )-l ,0)

IF( (OUANBT.GT.C) .OP. ( CUANNT.GT.C) ) GO TO 45

RETURN

45 CONTINUE

46 THTPPT = PTO^THFTA

60 PFTIJPN

F NO

SUBPOUT INF PPCSSM

C A ROUT INF TO SOLVP TH F EQUATION T « ( A- IN VF P S F ) *B ( ALL MATRICES

C A P F TPANSPOSFD) USING GAUSS-JORDAN ELIMINATION.

COMPLEX A,B,«IJNAy,ARAT,TMAT(flO,BO)

COMMON /MTX CO m / NP, NR I, A { BO, BO) R( 80* 80 ) COMMON /OUT COM/ IOUT

EQUIVALENCE ( L » FL ) , ( K , FK ) ( P ( 1 , 1 ) , T M AT ( 1 , 1 )

)

N = ?*NP

C START REDUCTION OF THE A MATRIX.

DO 80 I : 1 »

C SEARCH FOP THF MAXIMUM ELEMENT IN THE ITH ROW OF THE A-MATRIX.

ATJM&X - A ( 1,1)

J M A X = 1 on 10 J * 2> N

TF(OABS( A(T,J )).LF.CAPS( ATJMAX) ) GO TO 10

AIJMAX = A ( I , J ) JMAX * J 10 CONTINUE

TF AIJMAX IS 7FR0 ( AS IT WILL BE FOR ANY ROW (OP COLUMN) WHERE THE

INDEX M IS , G T . THE INDEX N, I.E.* THF LEGENDRE FUNCTIONS FORCE THOSE

MATRIX Fl.FMPNTS TO ZEPOJ.THFN THE MATRIX IS SINGULAR SO SOLVE THE

REDUCED MATRIX (ORDER = 2* ( N P ANK-M ) ) . TE(CA«S (AIJMAX) .GT. 0.0) GO TO 20

JMAX = I GO TQ 75

NORMALIZE THE ITH ROW BY AIJMAX (JMAX ELFMENT OF THE ITH ROW).

20 00 30 J = 1,N

A ( I » J ) = A(T.J)/AIJMAX

NORMALIZE THF ITH ROW 0 P B.

B ( I, J ) = B( I* J ) /AIJMAX

20 CONTTNUF

USE ROW TRANSFORMATIONS TO GET ZEROS ABOVE AND BELOW THE JMAX

PLEMPNT OF THE ITH ROW OF A, APPLY SAME ROW TRANSFORMATIONS

TO TH C B MATRIX.

DO 70 K = 1,N

TF(K.EQ.I) GO TO 70

APAT « -A(K,JMAX)

DO 50 J = 1,N

IE(CABS(A (I, J ) ) .LE.O.O) GO TO 50

A (K, J ) = AR AT* A ( I, J ) +A (W , J)

50 CONTINUE

A (K, JMAX) = 0.0

00 60 J=1,N

TE(CABS(B( I,J )) .LE.O.O) GO TO 60

R(K,J) = APAT*R( I, J)+B(K, J)

36
r>

r>

o

r>

o

oo

r>

->o

r>

60 CONTINUE

70 CONTINUE

STORE ROW COUNTER (I) IN TOP ELEMENT OF JMAX COLUMN. THUS,

THE TOP POW OF A WILL CONTAIN THE LOCATION OF THE PIVOT

(UNITY) Ft EM ENT OF EACH COLUMN (AFTER REDUCTION).

75 L » I STORE THE INTEGER I TN THE TOP ROW OF A.

A ( 1, J MA X ) = FL

80 CONTINUE

THE REDUCTION OF A IS COMPLFTE. PERFORM POW INTERCHANGES

AS INDICATED IN THE FIRST ROW OF A.

DO 1?0 I = 1,N

K = I PUT THF INTEGER VALUE IN A INTO K.

<?0 FK = A ( 1 , K ) IF(K-I) QO, 120, 100

IE K ( 1 » I ) is LESS than I, THEN THAT RQW HAS ALREADY BEEN

INVOLVED IN AN INTERCHANGE, AND WE USE K(1,K) UNTIL WE GET

A VALUE OE K GREATER THAN I (CORRESPONDING TO A ROW STORED

BELOW THE TTH ROW )

.

100 Dn 171 J * 1 ,

N

A RAT = R ( I , J ) P ( T, J ) = B ( K, J ) B(K,J) = A PAT

171 CONTINUE

120 CONTINUE

C the TRANSPOSED T MATRIX IS STORED IN B. TRANSPOSE TO GET THE T

C matrix AND STORE IN A,

DO 140 I = 1 ,

N

DO I BO J * 1 ,

N

A ( I, J ) * B( J, I

)

130 CONTINUE

140 CONTINUE

C TRANSFER THE T MATRIX FROM A INTO TMAT.

DO 160 I » 1 , N

DO 150 J = 1 N

T M AT ( I , J ) = A ( I , J ) 150 CONTINUE

160 CONTTNUF

CALL ADDPRC

RETURN

END

SUBROUTINE ADDPRC

C A ROUTINE TO OBTAIN THE SCATTERED FIELD COEFFICIENTS AND CALCULATE

C THE TOTAL N FA R FIELD IN THE AZIMUTHAL PLANE PHI = CONSTANT.

COMPLEX A,TMAT,AD1(80),AD?(80),FNGANS(361,2),H1(41),H2(41),BJ1(41)

1»BJ?(41 ),CI,THC,PHC,RC,CKPp,0CN,DCN2,RH0f PHOP, HANK (41 ) , ACANS,

2ACANSP»BSSLSP,CNEUMN,PC0MP(361),S1,S2,HANKP(41) ,BSSLPP(41),

3CNF!JMP(41),AD1X(80).A017(80)»A02Y(80),SI,PX,PY,P7,FG1(80),FG2(80),

4W,ETHETA,EPHT,ER»FTH(361)»FPH(361),FPC(361)

COMMON DTR, RTD, CP I COMMON /MTXCOM/ N R ANK , NR ANK I . A ( 80 , 80 ) , TM A T ( 80 , 8 0 ) , C MX N P M ( 80

)

COMMON /FNCCOM/ PNMLIG(41),RSSLSP(41),CNEUMM(41),BSLKPR(41),

1BSLKPT ( 41

)

37
r»o

o

COMMON / ROYCOM/ DC NR, DCN I , CK PRR , CK P R I * CK R, DCKR , CONK , AO VR B , WN , IB

COMMON /CM V COM/ NM,KMV,CMI< AO) *CMV»CM2,TWM» PRQDM

COMMON /THTCQM/ THETA, SINTH,C0STH,CDH(6),EPPS(6),NSECT,NDPS{6)

IdHFTAD, PH, KSFCT

COMMQN/UVCCOM/ACANS (361,2) » AC AN S R ( 361 ) , DLT A NG, DCN P2» DCNI2»NUANG

common /OUTCOM/ tout

DIMENSION ZXOLD(361)»ZYOLD(361),ZPOLD(361)

LOGICAL TEST

DATA TFST/.TPUE,/

NP? 2 ?*NPANK

GENERATE THE INCIDENT DIPOLE COEFFICIENTS ADI AND A D 2

CT = 10.0,1.0)

DCN = C MPL X

(

DCNP, DCN I ) nCNp = CMPLX(DCMP?,ncNI2)

W * WN*CSORT{ DCN2)

SI * CI*(W**3) /(CPITDCN2+R.85AE-12)

C RHOP IS THE PADIAL COORDINATE OF THE DIPOLE TIMFS THE WAVE NUMBER*

PHOP « CONK*1.0E*CSOPT (DCN2)

CKPRR = REAL ( RHOP)

C K PR I a AIMAG

(

RHOP ) CALL GENPKR

0^ 36 T*1,NPANKI

PSSIPP(I) = BSSLSP(I)

CNFUMO(I) a CNCIJMN ( I

)

HANK p ( I ) = BSSLPP<X) + CI*CNFIJMP(I)

36 CONTINUE

PHO « C0NK*1.1*CS0RT(DCN2)

CKPRR = REAL(RHO)

CKPRI « ATMAG(RHO)

CALL GENBKQ

Dn 37 I =1,NPANKI

HANK(T) = BSSLSD(I) +CI+CNEUMN { I ) 37 CONTINUE

CN a 0.0

C GENERATE LEGENDRE FUNCTIONS FOP DIPOLF ANGULAR COORDINATE THETAD.

TWPTA * DTP+THFTAD

CALL TR TG( THFT AD, S INTM,COSTH

)

CALL GENLGP

DO 35 N«1,NRANK

NP ® N+NPANK

CN a CN+1 .0

N1 « N+l

PI - CN*C0STH*PNMLLG(N1 )-(CN+CMV)*PNMLLG(N)

P? x cmv*PNMLLG(N1)

CKDR a RHOP

BJKN) = NtNl*BSSLPP(Nl) /CKPR

BJ2CN) x BSSIPP(N) - (N/CKPR)*BSSLPP(N1)

PX x l.o

py x o.o

P 7 x o.O

ADIX(N) * C0STH+P2 *BSSLPP (NX)

A D 1 7 ( N ) x -BSSLPP ( N1 )*SINTHtP?

38
o

o

o

r>

o

r>

ADIX(NP) - P J1 (N)*PNMLLG( N 1 ) ( SINTH**2) 4 B J 2 ( N ) *C OSTH* PI

AD17(NP ) = B J 1 C N)*S INTH*C0STH*PNMUG(N1)-BJ?(N)*SINTH*P1

AD2Y(N) » -P1*BSSLPP(N1

)

AD 2Y ( N P ) * -P?*PJ?(N)

A 0 1 ( N! ) = (PX*AD1X(N) 4 PZ + AQ1Z (N) )*SI

AOl(MP) = (PXtADlX(NP) 4 PZ*AD1Z(NP) )*SI

A D 2 ( N ) * PY*A[)?V( N ) + S I A D 2 ( N P ) = py*AD?Y( NP )*ST

35 CDMTTNUF

C

T HP SCATTERED P T F

1

0 COEFFICIENTS = THE TRANSITION MATRIX TIMES THE

INCinPNT FIELD COP FF ICIFNTS

.

DO 45 I = 1,NR2

SI = 0.0

S? = 0.0

DO 40 J = 1,NP2

SI = si + TMAT( I, J )*AD1< J

)

^ 2 = S 2 + TMATd, J )*AD2< J)

40 CONTINUE

FC-l(I) = SI

FG2 ( I ) * S

2

45 CONTINUE

FVALUATF THE SCATTFPFD FTFIO at EACH SCATTERING ANGLE.

DO 170 IU = 1,NIJANG

C GENERATE THF LEGENDRE MULTIPLIERS.

THFTT « DLTANG*(IU - 1)

IF (THFTT. LE. 181 .0) GO TO 6?

PHP = ph + 1«0.0

THET = 360,0 - THFTT

THETA » 0TRTTHFT

SINTH « SIN(THFTA)

COSTH a r os (TFFTA

)

CALL GFNLGP

PHI = CYV^PHP

GO TO 41

62 IF(THFTT) 05*85*95

P 5 COSTH » 1.0

KODE = 0

91 SINTH * 0.0

T H FT A « 0.0

GO TO 119

95 TF(THFTT- 180.0) 105*101*105

101 COSTH = -1.0

KODE » 180

GO TO 01

105 THETA = DTP*THETT

SINTH = SIN (THFTA

)

COSTH a COS(THETA)

119 CALL GENLGP

PHI « CMV4PH

39
o

r>

o

or»o

ooo

61 CAM TRTG(PHT,SIMPHMCOSPHT)

FNGANS<IU»1) * 0.0

FNGANS ( HU 2 ) * 0.

0

RCOMPflU) = 0.0

CN - 0.0

DH 160 N - 1,NPANK

NP * N+NPANK

N 1 = N + l CM = CM+1.0

PI = CN*cnSTH*PNMLLG(Nl)-(CN+CMV)*PNMLLG(N)

P 2 * CMVtPNMLLG(Nl)

AA*SINPHI*P1

Pr s CHSPHT«pi

cc=stnphi*p2

DD*COSPHT *P 2 EE « PNMLLG(N1)+SINTH*C0SPHI

F F = PNMLLG(N1 ) *SI NTH*SINPHI

IF(KMV.NE.O) GO IQ A H SGN = 1.0

IF(TMFTA) 68 ,AA»A8

46 I F f K 0 D

c .F0.0) GO TQ 66

FGN * (-1.0)**N

46 PF = COSPHT+SGN

F P * SINPHI+SGN

A® CONTINUE

CALCULATE IMF THETA COMPONENT OF THE SCATTERED FIFLD.

CKPO 3 PHO

H 1 ( N ) * N*N1*HANK(N1)/CKPP

H ? ( N ) = MANK(N) - (N /CK PP ) + HANK ( N1

)

fTHFTAMDD*F01(N)-CC*FG?(N))*HANK(N1)+(PP*FG1(NP)~AA*FG?(NP))*H2(N

1)

FNGANS(IIUI) = FMGANS(TUM) + F THP T A / CM XN PM ( N ) CALCULATE THF P«I COMPONENT OP THE SCATTERED FIELD.

FPHI = -<88*FG?(N} + AA*FG1{N))*HANK(N1 ) - ( C C * FG1 { N P ) +DD* F G2 ( N P ) ) *H 2 ( N

)

FNGANMHU2) = PNGANS ( IIJ> 2) + F PH I / C MX NR M ( N ) CALCULATE the R COMPONENT OF THF SCATTERED FIELD.

pp a (EE*FG1(NP)-FF*FG2(NP))*H1(N)

PCOMP(TU) * PCOMP(TU) + FP/CMXNRM(N)

160 CONTINUE

170 CONTINUF

C ACCUMULATE THF PENULTS EOP EACH M VALUE.

DO 17? ILJP * 1,NUANG

ACANSfIUP,l) * A C A N S ( T U P , 1 ) + FNGANS ( IUP» 1)

A C A N S ( I U p * 2 ) * ACANS(TUP»2) + FNGANS ( IUP» 2

)

ACANSR(IUP) * ACANSR (IUP 1+RCOMP ( IUP

)

172 CONTINUE

40
r>

r>

o

c PRINT TMP FIELD COMPONENTS AND THEIR MODULI SQUARED.

VP TTE C 6 . 17*5 ) KMV

175 «=OPMAT( 1H1, 35X, 35H********** ACCUMULATED SUMS EOP M =,I3,11H ***** 1*****/1HO,40X,17HTOTAL NEAR FIELDS/lHO, 1X,5HANGLE,12X, 15HTHETA COM

2 P ONE NT, 1 A V, 13 HD HI COMPONENT, 16X> 11 HR -COMPONENT, 2 2 X,5HT0TAL//

)

N C 0 N V > o

M CONV 0

LCONV * 0

SCANG = 0.0

CALL OT POLE(ETH,EPH,EPC)

C ADD DIPOLE FIELDS ETH, ERH, ER C TO OBTAIN TOTAL FIELDS.

DO 1Q0 JIJP = 1 , NU A NG

T«C « A C ANS { J UP, 1 ) ETH(JUP)

PMC = AC ANS ( JUP, 2 ) + EPH(JUP)

RC = ACANSR(JMP) + E R C ( J U D ) THC * THETA COMPONENT OF TOTAL FTELD.

PMC * PHI COMPONENT OF TOTAL FIFLD.

RC = R COMPONENT OF TOTAL FIELD,

y * c AP S (THC ) **?

V = CABS(PHC)**?

R * CABS(RC)**?

FSO = X+Y+R

WRITE (6,1*1) SCANG,X,Y,P,ESO

181 FORMAT (1H ,F6.2,A(1?X,F16.6) ) C tfst fop convergence at fach angle.

T F ( TEST ) GO TO 184

I F ( A o S ( X - 7V0LD( JUP ) ) .LF. (Y + l .0E-03) ) NCONV - NCONV + 1

I F ( ARS(Y - ZYOLD ( JUP ) > .1 E. ( Y*1 .0F-03 ) ) MCONV - MCONV + 1

I F ( A BS ( R - ZPHLD( JUP) ) .LE.(R*1.0F-03) ) LCONV = LCONV 1

184 ZXOLD(JUP) = X ZYOLD(JUP) * Y

7P OLD (JUP) * R SCANG * SCANG+DLTANG

100 CONTINUE

C TFST fop COMPLETE convergence of solution.

198 IF(NCONV.EO.NIJANG. AND. MCONV. EO.NUANG. AND. LCONV. FQ.NUANG) GO TO 194

TEST = .FALSE.

C

P FTURN

194 WP T TF ( 6, 200)

200 FORMAT( 1H0, 30H*** SOLUTION HAS CONVERGED ***)

14 CONTINUE

STOP

FND

SUBROUTINE TRIG( A, SINN,CQSN)

common dtr,rtd,cpi

SINN = S IN ( DTP* A)

COSN * COS(DTP*A)

I p (A— 180.0) 5,10,15

5 IF( A+180.0) 15,10,15

10 SINN = 0.0

15 I F ( A— 90,0) 20,25,30

20 IE(A+90.0) 30,25,30

25 COSN = 0.0

30 RETURN

END

41
ooo<"»r>oor>

SUBROUTINE CALENP

C CALCULATE TM(= ANGULAR FNOPOINTS FOR EACH SECTION OF THE BODY,

COMMON DTR, RTD» CP X COMMON /ROYCOM/ DC NR, DC NT ,CKPQR,CKPR I, CKR,DCKP , CONK, AOVPB* WN, IB

COMMON /TMTCOM/ THFTA,STNTM,C0STH,CDH(6), E0PS(6).NSECT,NDPS(6)

F P R S ( 1

) *CPI/?.0

C D H C ! ) * FPPS(1)/NDPS(1)

RETURN

END

subroutine gpnkr

C CALCULATE CKP and DCKP as a function OF thfta FOP A PROLATE sphfroid.

COMMON / BDY COM / DCNP, DCNI,CKPRR,CKPPI, CKP, DCKR, CONK, AOVPBjWN, IB

COMMON / THT COM / TH E T A , S

I

NT H

,

C OS T H

,

C DH ( 6 )

,

E P P $ ( 6 ) * NS EC T

,

ND PS ( 6

)

0 R = l.0/S0RT(CCSTH**2+(A0VRB*$INTH)**2>

C K R 2 rnNK«QR

0 C K R = -CONK*CnSTH*SINTH* ( AO VR R ** 2- 1 . 0 ) * C R**3

P t Ti l R N

END

SUBROUTINE DIROLE(ETH,FPH,ERC)

A SUBROUTINE TO CALCULATE THE FIELD DUE TO A DIPOLE AT COORDINATE

(DD,THD,0.0) AND COMRONFNTS (PX,py,P7)

FTH = THFTA component OF ELECTRIC F I F ID AT ORS F RV FR COORDINATE (PO

, SCANG, PH/PHP )

,

F PH * PHI COMPONENT

E P C = R-COMPONENT

W = WAVE NUMBER IN OBSERVERS MEDIUM,

D 3 DISTANCE RFTWFFN DIPOLF AND OBSERVER,

COMPLEX i|,V,W,CI,Py,PY,P7,FTH(361),EPH(361),ERC(361),DCN2,FAC,CE>

1 A C A N S , A C A N S P COMMON DTR,RTD,CPI

COMMON / BD.Y CO M / DC NR , DC N I , CK P R P, CK P R I , CK R , DC KR , CONK , AO VP B , WN, I B COMMON /THTCOM/ THETA, SINTH, COSTH,CDH( 6) , EPPS( 6 )> NSECT, NDPS ( 6)

1 ,THETAD» PH, KS FCT

C0MM0N/UVCCOM/ACANS( 3 61* ? ) , A C AN S P ( 36 1 ) » DLT ANG , DCN R 2, DCN 1 2, NU ANG

THD 3 THFT ADADTP

CT S o o O •v H • o PX 3 1 .

0

pv S 0.0

p 7 3 0. 0

PD 3 1 . 05*C0NK / WN

p 0 s 1 . 1*C0NM / WN

DO 11 I = 1 » NUANG

SCANG * (T-D + DLTAMG

IE (SCANG. LE.1B1.0) GO TO 8

SCAN = 360 . O-S C ANG

THETA « SCAN*DTR

S I NTH * EJN (THFTA

)

COSTH = COS(THFTA)

P H P - 180.0+PH

PHO 3 PHO+OTP

GO Tn q

42
B THETA * SCANG*OTR

STNTH * SIN(THETA)

CnSTH * CHS (THETA)

P H P * PH*DTP

q pY = pp*C TNTH*CPS ( PHP)-Rn*STN ( THO)

RY = RO*SINTH*S IN( PHC)

R 7 = RO*CD$TH - P 0 *CO S ( THO

)

p = SQPT (QX**2+RY**?+P7**2)

PT=PX*C!NTH*CPS(PHO)+RY*STNTH*$IN(PHO)+R7*COSTH

P?apy*C n ETHtCPS( PHO)«-RY*CHSTH*SIN ( PHO ) -P 7 *S I NTH

P3*PY*rOMPHn)-PX*SIN(pHn)

D0N2 = CMPL X (D^NR2,DCNT?)

y » WM*C$ORT( OCN2)

U * (W**?)/D + CI*W/0**2 - 1 / 0** 3

V * 3.0*(l/0**5 - CI*W/0**4) - ( W**2 ) /0**3

FAC = 1/ (4.0*f PT*DCN2*P.P5*E-12)

r E = CE XP ( c T*W*D

)

FTH(T) = FAC* ( °Y* ( U*CnSTH*COS (PHP)+V*RX*R2)+PY* (U*COSTH*S IN ( PHO)

1+V*PY*P?)+P7*( V*RZ*P?-U*SINTH) ) *CE

FPH(I ) = E A C* ( PX* ( V* R X* P 3-U* S IN ( P HO ) ) + P Y* ( U*CO S ( PHO ) + V*R Y*R3 ) + PZ

1V*R7*R3)*CE

EPC(T) = FAC*( PX*(U*SINTH*CaS (°HO ) +V* RX* R 1 ) +PY * (U *$ INTH*$IN< PHO)

IV* PY* PI ) + P7*(U*COETH + V*PZ*Pl) )*CE

11 CONTINUE

RETURN

FNP

43
N BS-1 14A (REV. 2-80

U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO.

NBSIR 81-1653

2. Performing Organ. Report No 3. Publication Date November 1981

4. TITLE AND SUBTITLE

Modelling of Oil Shale Retorts for Electromagnetic Sensing Techniques

5. AUTHOR(S)

H. Chew *

G. PERFORMING ORGANIZATION (If joint or other than NBS, see instructions) 7. Contract/Grant No. NATIONAL bureau of standards

PR No. 20-802C10417. 000

DEPARTMENT OF COMMERCE 8. Type of Report & Period Covered

WASHINGTON, D.C. 20234 Final Report

9. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street, City, State. ZIP)

U.S. Department of Energy

Laramie Energy Technology Center

Laramie, Wyoming 82071

10.

SUPPLEMENTARY NOTES

*0n sabbatical leave (1980-1981) from Department of Physics, Clarkson College,

Potsdam, NY 13676.

1 Document describes a computer program; SF-185, FIPS Software Summary, is attached.

11.

ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant

bibliography or literature survey, mention it here)

We report here some work on the modelling of oil shale retorts for

electromagnetic sensing techniques. The aim is to obtain useful information

about the contents of the retort (e.g., rubble size, void ratio, etc.)

by means of electromagnetic probes. In this work, the retort is modelled

by a spheroid with an average dielectric constant which depends on the

void ratio. The near field due to a radiating dipole source in the

vicinity of a spheroidal retort is computed using the Extended Boundary

Condition Method due to Waterman, Barber, and Yeh. Numerical results

are given at 4 MHz for a retort with major axis 45.7 (150 ft ), minor

axis 25.1 m (82.5 ft ), bulk dielectric constant 8.8 + 3 . 7j , and various

void ratios. The results indicate feasibility of determining the void

ratio by remote electromagnetic measurements. It is also believed that

this work may be of interest beyond the immediate context of oil shale

retort modelling.

12.

KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolon s)

Oil shale retorts; remote sensing; scattering

13. AVAILABILITY 14. NO. OF

PRINTED PAGES

Unlimited

| ] For Official Distribution. Do Not Release to NTIS

T Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C.

20402.

48

15. Price

[~X] Order From National Technical Information Service (NTIS), Springfield, VA. 22161 $6.50

Yr U.S. GOVERNMENT PRINTING OFFICE: 1981-576-001/1214 Region No. 8 USCOMM-DC 6043-P80