^ecltnlcciL ^lote 92e. 336

CALCULATION OF THE ADMITTANCE OF A

PARALLEL PLATE CAPACITOR CONTAINING

A TOROID-SHAPED SAMPLE

ERIC G. JOHNSON, JR.

U. S. DEPARTMENT OF COMMERCE

NATIONAL BUREAU OF STANDARDS
THE NATIONAL BUREAU OF STANDARDS

The National Bureau of Standards is a principal focal point in the Federal Government for assuring

maximum application of the physical and engineering sciences to the advancement of technology in industry and commerce. Its responsibilities include development and maintenance of the national stand- ards of measurement, and the provisions of means for making measurements consistent with those standards; determination of physical constants and properties of materials; development of methods

for testing materials, mechanisms, and structures, and making such tests as may be necessary, particularly for government agencies; cooperation in the establishment of standard practices for incorporation in codes and specifications; advisory service to government agencies on scientific and technical problems; invention and development of devices to serve special needs of the Government; assistance

to industry, business, and consumers in the development and acceptance of commercial standards and

simplified trade practice recommendations; administration of programs in cooperation with United

States business groups and standards organizations for the development of international standards of

practice; and maintenance of a clearinghouse for the collection and dissemination of scientific, technical, and engineering information. The scope of the Bureau's activities is suggested in the following

listing of its four Institutes and their organizational units. Institute for Basic Standards. Electricity. Metrology. Heat. Radiation Physics. Mechanics. Applied Mathematics. Atomic Physics. Physical Chemistry. Laboratory Astrophysics.* Radio Standards Laboratory: Radio Standards Physics; Radio Standards Engineering.** Office of Standard Ref- erence Data. Institute for Materials Research. Analytical Chemistry. Polymers. Metallurgy. Inorganic Materials. Reactor Radiations. Cryogenics.** Office of Standard Reference Materials.

Central Radio Propagation Laboratory.** Ionosphere Research and Propagation. Troposphere

and Space Telecommunications. Radio Systems. Upper Atmosphere and Space Physics.

Institute for Applied Technology. Textiles and Apparel Technology Center. Building Research.

Industrial Equipment. Information Technology. Performance Test Development. Instrumentation.

Transport Systems. Office of Technical Services. Office of Weights and Measures. Office of Engineering Standards. Office of Industrial Services.

* NBS Group, Joint Institute for Laboratory Astrophysics at the University of Colorado.

* * Located at Boulder, Colorado.
NATIONAL BUREAU OF STANDARDS

technical ^cte 336

ISSUED April 21, 1966

CALCULATION OF THE ADMITTANCE OF A PARALLEL

PLATE CAPACITOR CONTAINING A TOROID-SHAPED SAMPLE

Eric G. Johnson, Jr. Radio Standards Laboratory

Institute for Basic Standards

National Bureau of Standards

Boulder, Colorado

NBS Technical Notes are designed to supplement the Bureau's regular publications program. They provide a means for making available scientific data that are of transient or limited interest. Technical Notes may be

listed or referred to in the open literature.

For sale by the Superintendent of Documents, U.S. Government Printing Office Washington D.C. 20402

Price: 30<f

CONTENTS

PaEl

1 . Introduction 1

2. Initialization of the Problem 2

3. Analysis 7

4. Discussion 15

5. Acknowledgements 15

Appendix A i£

Appendix B (Fortran Listing) 19

Appendix C 33

Appendix D 35

in

Calculation of the Admittance of a Parallel

Plate Capacitor Containing a Toroid- Shaped Sample

Eric G. Johnson, Jr.

The theoretical admittance of a toroid- shaped sample

in a parallel plate capacitor system is developed. This calculation represents a first and second order perturbation correction to the usual quasistatic formula. The FORTRAN

listing of the machine program is included. The analysis is

done under the condition that the sample is not in contact

with the capacitor plates and can be described by a complex

permeability, \x-^, and complex permittivity, e^.

Key Words:

Fortran

Lossy materials

Parallel plate capacitor

Permeability

Permittivity

Theoretical admittance

1. Introduction

The purpose of this calculation is to permit an accurate determination of the permittivity and permeability of a material sample for

the frequency range in which the capacitor and inductor system are

generally applicable.* This particular calculation considers

*This technical note is a companion to the NBS Technical Note #311,

"Computation of the Permeability and Permittivity of a Relatively Small

Ring Sample in a Toroidal Coil," which contains a similar analysis for

the inductor structure of a toroid- shaped sample. Both the capacitor

and inductor measurements are necessary if one is to get a complete

determination of the permeability and permittivity of a sample.
the parallel plate capacitor system in which there has been placed a

toroidal -shaped sample. Our purpose is to correct the usual quasistatic formulae (zero order solution) for characterizing a capacitor to

at least second order in perturbation theory. The zero and first order

calculation is performed under the assumption that there is no height

dependence of the electromagnetic field. Only radial dependence is

permitted. However, because we wish to be able to study lossy materials (complex |jl and e A/r ) we have chosen the geometry of the experiment to be such that the sample is not in contact with either of the

metal plates of the capacitor system. Because of this configuration, it

is necessary to make a guess of an effective permeability [x and permittivity g for the zero and first order solution to the problem. The

second order term primarily attempts to correct for the fact that there

is indeed height dependence for the electromagnetic field. This second

order is derived under a relaxation of the requirement that tangential

component of the magnetic field is continuous. It is expected that when

this second order term is important there is the corresponding circumstance that the sample tends to not have an internal electromagnetic

field, i. e. , it is as though there is a surface current acting on the

sample surface. So to aid in the convergence of the second order term

and hopefully to make it a more reliable estimate of the corrections

necessary to the zero and first order term, this calculation errs in

the direction of assuming that there is a surface current in the second

order term.

2. Initialization of the Problem

Assume an experimental system with the following boundary

conditions:

1. We have two finite radius parallel plates. Each plate

has infinite conductivity.
2. A rectangular cross- sectioned toroid is placed between the

plates.

3. Two complex frequency dependent parameters- -permeability, [i and permittivity, e, ,-- are used to characterize

the material in the toroid.

4. The excitation mechanism of the capacitor system is such

that no angular field dependence exists and is such the Hep is

the only magnetic variable. This means that the basic

electromagnetic field equations namely (1) will be reduced

to a single equation, (2) below. Equation (1) assumes the

exp (ia)t) time dependence. Thus,

V X E = - i cu M-H, and

V XH =i u3sE

i =JT"

(1)

(natural units are used, namely, c = 1. Hence all units are

powers of cm. The uu is cm" ; \i, e are dimensionless).

5. The coordinates of interest are shown in Figures 1 and 2.

SIDE VIEW

r =

i r

1 2

c

o o A

^M' M

C

Guard'

Ring

(The system is symmetric around the r = position.)

Figure 1

3
TOP VIEW OF SYSTEM

GUARD RING

Figure 2
The r, and r^ are the radial coordinates of the toroid sample A.

1 2

The r is the radius at which the system has its voltage drop prescribed.

The z and z are the z axis coordinates of the toroid sample A.

The z , z and z , z are the z axis coordinates of two toroid-shaped

material B. They have no loss in their e and u . The relative

z - z z - z

2 1 4 3

size of and is about 2. 5%. The structure is z - z z - z

3 2 3 2

used to provide a precise measure of the distances z , - z~ and

z - z . We characterize the permeability and permittivity of the

space r > r > -0 and z . < z < z by:

|jl = 1 and 6=1, for region C;

u = u and e= e for region B;

and u = u._, and 6 = e ,,. for region A. M M»

The regions A, B, and C are indicated in Figure 1. A functional form

for u is u = P(r) [P(z) (u M

- Uj) + u

} - l] + 1 .

The same form for e is possible. Here we have defined

P(r)= 8( r - r ) - 8(r - r ) and

P(z)= 6(z - z

3

) - 6(z - z

2

),

where is the step function.
3. Analysis

Upon reducing the (1) with assumptions 1 to 5 and

setting f s rH„ the following differential equation results:

— °r re

1 °r u z e12z

(2)

The boundary conditions are g f = at z = z , z -and — -* as r -* 0.

The admittance, Y,of this system at r = r is defined as the ratio of the

peak current to peak driving voltage. Thus,

Y = 2 rriujf(r

3 )/ J.*

-3 f dz

r r r = r.

2TriU)r.

Z4' z

l

[f(r

3 )/fer

f)(-

3 )] <3)

The position r is chosen so that the influence of the toroid

sample on the field structure has no resulting z dependence. That is,

r is sufficiently far from the sample so that the only type of field that

exists at r is pure propagating waves with no z dependence

Jr > r + 4(z - z )l. By fixing the voltage drop point at r we solve

I -J U ^t J. I '-J a basic problem for making the admittance of a parallel plate capacitor

unique. Thus the voltage drop at r is given by the line integral

1U) J er

fdz

r = r

3*

This means we have chosen a path at r = r that is only along the z direction. With this specification and noting that the admittance contribution outside r = r is the same with or without the sample, we can

subtract the two admittance measurements and equate the result to the

corresponding subtraction of the two theoretical admittances.
In order to eliminate mathematical consideration of fringing

fields, we assume that the capacitor is infinite. We then arbitrarily

choose r as the point where we define the interior of the capacitor

system. We must experimentally fix r and still be able to neglect the

effects of fringing fields. Practically, this is accomplished by having

a capacitor of radius r and by using a guard ring. We now solve the

equation for f. A first order guess to f is to assume that the sample

completely fills both regions A and B with the weighted average

parameters e, p., as is defined below. We can solve this f exactly.

We next assume that f may be written as f = f Q+ e(r, z) g (r, z)„where

g is assumed to be a perturbation connection to the f solution. By

allowing this transformation and making g and its derivative continuous

we are granting the existence of surface current on the material sample.

We thus have the following equation for g

^^nJ^—^J-H-^f*-*)*.*-^^-^,^ (4)

where f obeys the equation

n -4. »n

r o + °* ^ "** ~ °- (5) "• ^n 4. ^n ^o +

""^ -?e, - O

Here, £ = PCnJC^-l) f | ^= PCn)(Aa

-i) + |,

where ^ = €|/(| + ( «, /^ -Qfy^.))

and

The solution to (5) is

"fo = n. L A T, ^ < «o n. Ja! ) +• £ K", C * *o £s? O J for each region.

Thus,

-fa - a. T, (cort) n<rv,

7
We require f and — 9 f to be continuous across the r, and r^ boundary,

o e r o 12

This implies with definitions x = iW r and x = iw r /|j. e

that T,C^J) = AT,(\J + !!(,«,)

T (X,)= /^[ AT.G,) -8K„(x,)l > and that

We solve for A, B as follows:

The A, B is given by

where

The Wronskian is

V/ = x.c*/tc.c*) -t,(%) W D c*/ ^cov^+T.ooMO (9) - -L for all z including z -» 0.

2

s

Thus we can write D - "*" anc^ ^k. ~ ~~—~ • (10)

x

X

Next we reduce the inhomogeneous term as follows:

ip r (11)

8
Equation (4) can be solved in terms of the expansion

3 v\=o

where

*„•=- nir

V*.

and W« =|f a- ~C

The equations to fix R (r) are now given as

_ V. Here, A ,, D , B , D are defined:

nn n n n

)

"en, € ° n.-i?. c n, where

Here 5 =0 n ^ and

no »

5 =1 n = . no

We further compute A /to get

(13)

(14)

K** ^jf*' C^S

M' €

-j [^^:<^^ -^)

+ ac^sCwh^vj^,)

a.

•- *w

- V? * »*>' ^v» + VP H

. (15)

when n ^ n'. When n

/

- n , we have

r*w»i

ST c 1 - L) 1 +-

(16)

+ *i *, -i
Thus we now define: "\^ B

** n - «* ^ = 4- \f X^-

X A

X

r >0-«x 4,n , x M

=• + f5^ (17)

i 7 n ^ / '

l "vi and the Green s function equations

tn^i$i ~ \^n)lGn ^,0 - -^-SV, (18)

^n^*)!^- _l;T_ -AJOO]- -S<n-*i . (19)

The solution to these equations is given as

The normalization condition on T given by

n

or since the term in the bracket is the Wronskian, W, and its solution is

M ~^~ W or ^ ~ ~ , we have the requirement that C = - 1.

To find W, we fix T (r) and C (r ) as follows:

n n

Cn <*) - Cj T,CAn n) * c * K,cvn n), n, < * < n t_

T^n) " S^x.CV) + Sn

W K,(Xh a),

~ (20)

T„<0 = S*T,(V) + S

1; K,CXh n), wx ^ 'n. .

To check if rW : 1, we consider the case r -* 0. We have

v_ h

- - .- J

'

hence rW = -1 as desired.

^\.«To -C^ )-T, -^XJCT,^,,)} c

b J-o vv h

. x, ^L^ C,,j^ - J-, (21)

10
The conditions for connection of the three regions are that T

n' C , and their derivatives must be continuous across the r , r boundaries. As a consequence, we have the following equations:

(22)

An

Upon solving these equations, we find that:

u

^

J r i See (10) for this reduction.

11
c

/ K u.«,; w,(A^)\/c-x,a,n ,)+^K,a^,)

>c 1 \ ^n

V.

Oa

= K.T, 4T4 K, - J-Hr.I.'-X.lc'l - "^ J) and finally,

I" »

An

We next construct the integral equation for R . That is

n

(23)

Now we assume R has only outgoing waves; therefore, we have the

n

following integral equation:

I? GO - S^Gjn^JKa)^- 6*1 + E^Gv^/Vn.) + -pj; CjVl ) * n ^

' ~ (24)

The solution of this integral equation by standard iteration is for n f 0:

^(25)

and for n = 0:

The final form for f,at r = r is

3»

f '-^f?D^ ^[bH,(^^) +I/IUA,)] ^. ^

12
The admittance is given with r = r as follows:

Y

2_-rr tu3

V*.

£%,<, MW^

) Vl (;«Bj) 1

- l*° k.CJ-^.) vtd c ^ ^ (27)

b a The basic computational problem is to determine B and B . We make

the expansion of the f and R function explicit so that we can determine

B ;

(28)

(29)

f?o Cn,) = UjC.Cn^J £„ S^-TjiO K £ ( n ) <*

•d + ^ l.<0 ^.T.C^)

(30)

+ ^2 ^^ Oj) ^T Cn)R^CK3 <$:

Here x = X r, n ^ 0, and X = ill) /

Thus we define

2 2

^o

= Ui ^iCi^itj ) ^

1-1

' + «** H

1

' J and note that B

b

i

The FT

1

', H

(2 ' are defined as follows:

s now

(31)

13
- - $ A

We define

^ (TO f J* ^Ca) T,C^ ain.,

The integrals are:

p/Vn.) = [ ^(X,(Xn)Xe Cx)x-X,T CXJT| Cx))]^ /?„,

Px\0 -" [-tu(x,CxJ K.C*)1 + >h T ( xn) tC,(x))]72n ,

We have defined:

(32)

(33)

(34)

X„Ciu) = n. [ C„<V) l-t^f*) +-T/IO H^f-l],

» f^J ^.a _

" (36)

hence, l

K

te = Z hn H^. (37) "O o Now, in addition we have

^'J = B HTM + X^T.fN) + X5

L T C^)

(3g) ^ = *r i»?< K, +. c h^>.

(39)

14
4. Discussion

The analysis has gone as far as possible without numerical

integration of (36). We summarize the basic necessary results in

Appendix A. In Appendix B we have a FORTRAN 63 listing of the

machine program. In Appendix C is a sample calculation. (Here \x = UM +IUM* and e^ - EM + IEM* in the computer program. UM, EM are

the real parts and UM*, EM* are the imaginary parts . They must be

negative to describe a lossy system. Appendix D gives a flow chart of

the input data necessary to run the machine program.

The quasistatic formula used is (r£ - r^)

qs »

4

_ ^

2

The computation has explicitly separated out the highest order

(2) correction, H , for the determination of the admittance of the sample.

This correction can be used to get an estimate of the accuracy of the

admittance calculation. This correction is calculated by Simpson

(2) intervals to an accuracy of about 3%. Thus if this H term is less than

10%, and subject to discussion of the introduction we can conclude that

the admittance is known to at least 0. 5%. In addition, one has an option

of varying r by 10%. Under this option one can get an estimate of how

sensitive the experimental system is to the inexact knowledge of r . As a consequence we have a simple measure of the degree the radiation

part of the electromagnetic field is present.

The program prints out the results in mho's for the admittance.

5. Acknowledgements

I wish to thank Dr. S. Jarvis for many helpful discussions,

A.Rasmussen for providing me data to use in testing the machine program,

G. Hardy for help in programming, and Drs. H. Bussey, D. Kerns , and

J. Wait for useful editorial comments.

15
APPENDIX A

SUMMARY OF FORMULAE USED IN MACHINE PROGRAM

(

A)=*j*i<*^ ^(^)Vf i ** ( iw^^))

TCxO X (xjA A x , c \) + g k.otj

aK x ( **ix ^ *.<*!> \f r,oo

6 Te C*0 -x.cx,; A /S

(Al)

(A2)

P."- 4 x,<0 + 8 k.co V = cco\Twlz

n. T?,<n. < nz

(A3)

(A4)

^«r-*,

tf* ^ D

wh

E *~~ \>(* -|) Fn'+ (ll-O Jn. I Vo = -L_

f; - £

»V

-<j--r,

>„*- = If^-w1

" j \ = + V^

l*»\

n, < n. < tlv

OG K'

_ A.>-«^C$,-*,)H

(A5)

(A6)

(A 7-)

(A 8)

(A9)

(A 10)

(All)

v,

16
(!

CL. ^ ^^

IC C\V) MVO V M\0 \

A T.(\av) -x.a^OA 4^ K c^x)

(A12)

A.

c-t v^ocx^,) -r.a.njA ^ 1(^ n)

(A13)

Ah

C

C

s

- \*Y (A 14)

- \*H

(A 16)

6^=8+ H

C0

+

^

H

^v

K° - Bjw t -D*TD Cno +Oo

b To rnO

H

t v ^ g€ c ^^ , h

^Ql) j + sV hJ%j - H'.%3)

><a; * A., *( *

H>,

5" '-"T^tO R„ro Jo.

(A17)

(A 18)

(A 19)

(A 2 0)

(A21)

(A22)

17
-v s pxYo

pr^)-[-n.CT,CxJ K.Olft 4\ £<*J ^<7) )]/^ (A28)

(A23)

(A24)

(A25)

(A2 6)

(A27)

PjVn )=["*-( K^J T/xjA +\ KP Cxh ; X(( x j )J /^ (A29)

RjViO = [-^(k.c* )Kw.i- x

H K c*,)k,(*-y)y ?n (Aso)

2Tv, A^ — A x^ = ^n. * = X n_

18
APPENDIX B

FORTRAN LISTING

PROGRAM CAPACTOR CAP

COMMON R,Z»F,Ul, SK,ALM,BB,BA, AL ,AA ,AB, AD ,ADB,BKO, BK1 ,BIO »BI 1 , CAP

1CA,CB,SA,SB,X1,X2» Z21 .Z31 ,Z41 *UM.E1 ,EM,RATIO,Yl »Z32 , IB ,ONA , CAP

2SNl»SN2,DU»DE,AKtCFT»SUM»F2tICYLE» TEMP,ST,C,B, SUB ,11 CAP

3 ,STEP2

DIMENSION R(4) ,Z(4)

,

SK ( 100 , 2 ) »ALM( 100, 2 ) » AL ( 100 ) AA ( 100 ) » CAP

1AB( 100 ) ,AD( 100) ,ADB( 100) , BK0(5 0O) ,BK1 (500) ,BI 0(500) ,BI 1(500)

,

CAP

2CA( 100,2) ,CB( 100,2) ,SA ( 100 ,2 ) ,SB ( 100,2 ) ,X1 ( 2 ) , X2 ( 2 ) RAT 10(2) ,Y2(3)CAP

3,SN1( 100,2 ),SN2( 100,2 ) ,CFT ( 6 ) ,SUM(7) ,SUB( 2 ) ,ST (6) AP

4 ,C(3,2)

TYPE COMPLEX Ul, ALM,BB,BA, AL » AA, AB , AD, ADB »BK0» BK1 , B I , BI 1 , CA

,

CAP

1CB,SA,SB,X1,X2, X, UM,E1,EM, RAT 10 ,Y 1 , AREA , Y2

,

CAP

2Y3, ONA,DU,DE, SUM, TEMP, ST ,SUB CAP

3,C,CFT ,E2,U2

DATA ( CON=2.0944E-10) ,(PI=1.668E-02 ) C THIS PROGRAM COMPUTES THE ADMITTANCE IN MHOS OF A TOROID SAMPLE CAP

C IN PARALLEL PLATE CAPACITOR. THE INPUT DATA IS THE THREE RAD1(CM),CAP

C THE FOUR Z AXIS VALUES. THE FREQUENCY (HZ), AND THE TWO SETS OF CAP

C VALUES FOR THE PERMITTIVITY AND PERMIBILITY. THE OUTPUT DATA CAP

C IS THE ADMITTANCES— SECOND ORDER CORRECT ION, QUAS ISTAT IC

,

CAP

C AND TOTAL. CAP

SENSE LIGHT 1 ONA=(1.»0.

)

2 READ 1,(R( I ) ,1=1,3) ,Z, IA,IB,E1,U1 CAP

1 FORMAT(3E18.5/4E18.5/ 214/ 2C ( E18.5 ,E18.5 )) CAP

C R IS RADIUS VALUES,Z IS HEIGHT VALUES IA=-1 IS ADMITTANCE TEST

C IA=0 SKIPS TEST OF SENSITIVITY IB=1 SKIPS HIGHER MODE =0 DOES IT IT=-1 CAP

101 READ 102,E,IC,IS,ITES CAP

C E IS FREQUENCY IN HERTZ IC= NUMBER OF DIFFERENT EM AND UM

C IS IS PAGE NUMBER ITES=-1 TO EXIT AFTER NEXT BATCH OF DATA

C ITES=0 NEW SPACIAL DIMENSION =+1 NEW FREQUENCY

F= CON*E $ F2=F**2 CAP

102 F0RMAT(E18.5,3I4) CAP

6 PRINT 7,IS,E,(R( I ) ,

I =1 ,3 ) »Z,E1 »U1 CAP

7 FORMAT* 51H1ADMITTANCE FOR CAPACITOR WITH TOROID GEOMETRY PAGE I4/CAP

118H THE FREQUENCY IS E18.4, 4HHZ. / CAP

227H THE RADIAL DIMENSIONS ARE 3E18.4 / CAP

327H THE Z AXIS DIMENSIONS ARE 4E18.4 /14H SPACER El IS CAP

4C(E18.4,E18.4) ,7H Ul IS C ( E18 .4 , E18 .4 )

/

CAP

591H THE LENGTH MEASUREMENTS ARE IN CM. ADMITTANCES ARE IN MHOS. TCAP

6HE REST ARE DIMENSIONLESS. //) CAP

PRINT 701 CAP

701 FORMAT(115H MU EPSILON QUACAP

1SI STATIC ADM CORRECT ADM DBL SUM ADM / CAP

2125H UM +1 UM* EM +1 EM* Yl +ICAP

3 Y2 YC1 +1 YC2 DY1 +1 DY2 CNT >CAP

ITT=1 CAP

9 IS=IS+1 CAP

IF(IT)71,72,13 CAP

71 Z41= Z(4)-Z(l) $ B= 1./Z41 CAP

STEP2=(R(2)-R(1) )/32.

FIX= PI /Z41 CAP

212 Z21= Z(2)-Z(l)$ Z32 = Z(3)-Z(2 ) CAP

19
Z31 = Z(3)-Z(l) $ AK=(3.14159265)/Z41 CAP

72 AREA= FIX *( R < 2 ) **2-R ( 1

) **2 > *< 0. ,0.5 ) *F CAP

12 DO 13 11 = 1, IC CAP

211 F0RMAT(2C(E18.5»E18.5)

)

CAP

225 READ 211, EM,UM CAP

Z5=Z32/Z41

E2=E1/(0NA+Z5*(E1/EM-0NA)

)

U2=U1+Z5*(UM-U1

)

Y1=AREA*(E2-0NA)

15 CALL ADMIT(E2,U2)

216 X=R(3) *F*(0.,1.) $ CALL BESSL2 ( X ,500 ) $ F3=FIX*R(3)

DO 215 J=l,2 CAP

Y2( J)=F3*( (RAT 10 ( J ) *BK 1 ( 500 ) +BI 1 (500) ) / ( RAT I 0( J ) * ( -BKO( 5 00 ) > + CAP

1BI0( 500) )-BI 1(500)/BI0( 500)

)

IF( IB)215, 215,213

215 CONTINUE

IF( IA) 217,218,218 CAP

217 X= 0.9*X S F3=0.9*F3 SCALL BESSL2 ( X , 500 ) $ ITT=ITT+1

Y2(3 )=F3*( (RATIO(2)*BK.1(500)+BI 1 (500) )/( RATIO ( 2)*<-BK0( 5 00) )+ CAP

1BI 0(500) )-BI 1 (500) /BIO ( 500)

)

PER= CABS( ( Y2(2)-Y2(3) ) /Y2(2) ) CAP

219 FORMAT160H NEXT ADMITTANCE FRACTIONAL CHANGE WITH R(3)= 0.9*R(3) CAP

IIS E12.4) CAP

PRINT 219, PER CAP

GO TO 218

213 Y2(2)=Y2(1)

218 Y3 = Y2(2)-Y2('l ) CAP

302 PRINT 8,UM,EM,Y1,Y2(2 ) ,Y3,ICYLE CAP

8 F0RMAT(C(E13.4,E12.4) ,4C( E12.4.E12.4) ,14) CAP

ITT=ITT+1 CAP

IF( ITT-46) 13,6,6 CAP

13 CONTINUE CAP

IF( ITES)240.2,103 CAP

103 IT=0 CAP

GO TO 101 CAP

240 CALL EXIT CAP

END CAP

SUBROUTINE SAVE(J) SAVE

COMMON R,Z,F,U1, SK ,ALM , BB , BA , AL , AA , AB, AD ,ADB, BKO , BKl , BI , BI 1 , CAP

1CA,CB,5A,SB,X1,X2, Z21,Z31,Z41,UM,E1,EM,RATI0,Y1,Z32,IB,0NA, CAP

2SN1,SN2,DU,DE,AK,CFT,SUM,F2»ICYLE, TEMP,ST,C,B, SUB ,IT CAP

DIMENSION R(4),Z(4), SK ( 100 , 2 ) , ALM { 100 , 2 ) , AL ( 100 ) , AA ( 100 ) » CAP

1AB( 100 ) ,AD( 100) ,ADB( 100) , BKO (500) ,BKl ( 500 ) ,BI < 500) ,BI 1(500) , CAP

2CA( 100,2) ,CB( 100,2 ) ,SA( 100,2) ,SB( 100,2 ) »XK2) ,X2<2) ,RATI0(2) ,Y2( 3 )CAP

3,SN1 ( 100,2 ) ,SN2( 100,2 ) ,CFT (6) ,SUM(7 ) ,SUB(2 ) ,ST(6) AP

4 ,C(3,2)

TYPE COMPLEX Ul, ALM,BB,BA, AL , AA , AB , AD, ADB ,BK0 , BKl , B I , BI 1 ,CA , CAP

1CB,SA,SB,X1,X2, X, UM,E1,EM, RAT 10 , Yl , AREA , Y2 , CAP

2Y3, 0NA,DU,DE, SUM, TEMP, ST ,SUB CAP

3.,C,CFT

IF( IT)50,101,101

50 J1=0

I T = l B1=B*1. 414213562 $B2=S0RTF(b) $B3=B2*1 .414213562 $B4=B/2.0

20
X6=Z32*B $SN2(1»1)=X6 $SN2 < 1 ,2 ) =B2* ( Z32)

2 GO TO 9 SAVE

101 IF( J-Jl) 10,10,3

3 SK(J,1)=S<(J-1,1)+AK SAVE

SK(J,2)=SK(J,1)**2 SAVE

X4=SK( J,1)*Z21 SAVE

X5=SK( J,1)*Z31 SAVE

SN1(J»2)=(SINF(X5)-SINF(X4> )/SK(J»l) $SN2(J,2)= SN1(J,2)*B3 SAVE

SN1I J,l )=( (SINF( 2.*X5)-SINF(2.*X4) ) /( SK(J,1)))*B4 SAVE

SN2(J,1>= (SN1(J,1)+X6 SAVE

SN1(J»2)= SNKJ,2)*B1 SAVE

9 J1=J1+1 SAVE

10 RETURN SAVE

END SAVE

SUBROUTINE I NTGR ( J

)

INTGR

COMMON R,Z,F,U1, SK, ALM , BB»BA, AL » AA , AB, AD, ADB, BKO » BK 1 , BI »B I 1 , CAP

1CA,CB,SA,SB,X1,X2, Z21 ,Z31 ,Z41 ,UM,E1 , EM, RAT 10, Yl ,Z32 , IB, ON A, CAP

2SN1,SN2,DU,DE,AK,CFT,SUM,F2,ICYLE, TEMP,ST,C,B, SUB CAP

3,IT,STEP2»T,HC,AF,P

DIMENSION R(4),Z(4), SK ( 100 , 2 ) ,ALM ( 100 , 2 ) , AL ( 100 ) , AA ( 100 ) , CAP

1AB( 100 ) ,AD( 100) ,ADB( 100) , BKO (500) ,BK1(500) ,BI 0(500) tBI 1(500)

,

CAP

2CA( 100,2) ,CS( 100,2) ,SA( 100,2) ,SB(100,2),X1(2),X2(2),RATIO(2),Y2(3)CAP

3,SN1 ( 100,2 ),SN2( 100,2) ,CFT(6) ,SUM(7) ,SUB ( 2 ) *ST ( 6

)

4,H(8,2),T(6,2),HC(5), P ( 4 ) ,C( 3,2 ) ,AF( 33

)

TYPE COMPLEX Ul, ALM,BB,BA, AL ,AA, AB ,AD, ADB,BK0,BK1 , BIO.BI 1 ,CA, CAP

1CB,SA,SB,X1,X2,T0,P1,P2,P3,P4,X,X3,UM,E1,EM, RAT 10 ,Y1 , AREA , Y2 , CAP

2Y3, 0NA»DU,DE,CFT,SUM,H,T,CR,HC,TEMP.ST ,P,C,AF,SUB CAP

TYPE COMPLEX HF1 ,HF2 ,HF3 ,HF4,HF5 , HFB3 , HFB4 ,HFB5 ,HFB6

,

INTGR

1 RF,TF,CF,RFB,TFB,CFB,X0,X5,RT,RTB,DM1,DM2 ,RA1 ,RA2,RC,AL4 INTGR

TF(K,K1)= SA(K,2)*BI1(K1)+SB(K,2)*BK1(K1

)

INTGR

CF(K,K1)= CA(K.2)*BI1(K1)+CB(K,2)*BK1(K1

)

INTGR

P1(K1,K2)= BI1(K1)*BI0(K2)*X0- X*BI0 < Kl ) *BI 1 ( K2

)

INTGR

P2(K1,K2)=-(BI1(K1)*BK0(K2)*X0+X*BI0(K1)*BK1(K2) ) P3(K1,K2)= BKKK1)*BI0(K2>*X0+X*BK0(K1 )*BI1 (K2) INTGR

P4(K1,K2)=-(BK1(K1)*BKO(K2>*XO-X*BKO(K1)*BK1(K2) > HF5( I )=(SA( 100,2 )*P( I )+SB( 100,2 )*P( 1+1 )

)

HF3(K) =(SA(K,2)*H(5,1 ) +SB < K ,2 ) *H ( 6 , 1 ) )*ADB(K) INTGR

HF41K) =(CA(K,2)*H(5»1 >+CB ( K ,2 ) *H ( 6 » 1 ) ) *ADB ( K

)

INTGR

IF( J-l ) 1,1 .3 84 SUM(4)=(AF(1)+AF(33) )/(2.,0.) INTGR

DO 85 L=5,30,4 INTGR

85 SUM(4)=SUM(4)+AF(L) INTGR

SUM(5)= 2.*(SUM(4)+2.*SUM(3)

)

INTGR

SUM(6)= SUM(4)+SUM(3)+2.*SUM(2

)

INTGR

SUM(1)= (SUM(6)-SUM(5) ) / ( SUM ( 6 ) *15 .

)

INTGR

TOO = CABS(SUMd)) $ J3=J-1 INTGR

IF( TOO -0.02) 87,87,86

86 PRINT 90,TOO,J3 INTGR

90 FORMAT!* THE INTEGRATION ERROR IS* E10.2 *THE MODE NUMBER IS* 14) INTGR

87 ST(4)= (STEP*0. 666666667 )*SUM(6) INTGR

RETURN INTGR

1 18=1 S GO TO 2 3 18 = 2 2 JB1=J+100 INTGR

21
JB2=J+300 INTGR

R(4)=R(1) INTGR

JR=JB1 $ JS=200

16=1 1 17=2 INTGR

GO TO 10 INTGR

20 R(4)=R(2) INTGR

H(l,l)= H(3,l) $ H(7,l )=H(4,1) INTGR

T(l,l)= T(3,l) INTGR

16=2 INTGR

HC(2 )=AL(J)*T(1,1)-H(1,1)

JR=JB2 $JS=400 INTGR

GO TO 10 INTGR

30 H(2,1)=H(3,1) $H(8,1 )=H{4,1) INTGR

T(2,l)=T(3f 1 ) INTGR

GO TO (90.91 ) ,18

91 C(2,1)=C(3,1) INTGR

HC(3)= AB( J)*C(2.1)+H(8»1) INTGR

17=1 INTGR

R(4)=R<2) INTGR

RA1 = (C(3»1)*(H(3,1)+HC(2) >+T < 3 » 1 ) * ( HC ( 3

) -H( 4, 1 ) ) ) *R ( 4

)

JR= J+400 $JS=500 INTGR

AF(1)= RA1*HC(5)

STEP=0.15/ CABS( ALM(J,2)

)

INTGR

IF(STEP-STEP2)32,32,31

31 STEP=STEP2

32 DO 80 L=2.33 INTGR

R(4)=R(4)- STEP INTGR

16=3 INTGR

GO TO 10 INTGR

40 X=ALM( 1»2)*R(4)

CALL BESSL2(X,401) $ T ( 6, 1 ) =TF ( 1 .40 1

)

RA1 = (C(3.1)*(H(3,1 )+HC(2) )+T ( 3 . 1 ) * (HC ( 3 >~H( 4. 1 ) ) ) *R ( 4

)

80 AF(l.)=RAl*T(6»l)

DO 81 L=2.3 INTGR

81 SUM(L)=(0..0.) INTGR

DO 82 L=2.33,2 INTGR

82 5UM(2)=SUM(2)+AF(L) INTGR

DO 83 L=3.33.4 INTGR

83 SUM(3)=SUM(3)+AF(L) INTGR

GO TO 84 INTGR

10 X0= ALM(100.2)*R(4)

X = ALM( J,2)*R(4) INTGR

GO TO ( 11 .12) ,17

11 CALL BESSL2(X0,500)

CALL BESSL2U .JR ) INTGR

12 P(l)= PKJR.Js ) P(2)= P2( JR,Js ) P(3)= P3(JR,JS ) P(4)= P4(JR.JS ) H(5,1)=HF5(1) INTGR

H(6,l )=HF5(3) INTGR

H(3,1)=HF3(J) $ T(3,1)=TF( J,JR) INTGR

GO TO (851,861) , 18

861 H(4,1)=HF4(J) $ C ( 3 , 1 ) =CF ( J, JR

)

INTGR

22
851 GO TO (20.30,401,16 INTGR

90 HC(3)=AB(1)*T(2,1)

ST(5)=HC(3)+HC(2)+H(2,1)

HC(5)=T(2,1)

RETURN

END INTGR

SUBROUTINE COEF(J) COEF

COMMON R»Z,F,U1, SK ,ALM, BB»BA, AL ,AA ,AB, AD, ADB, BK0»BK1 ,BI ,BI 1

,

CAP

1CA,CB,SA,SB,X1,X2, Z21 ,Z31,Z41 ,UM ,E1 .EM, RAT 10, Yl ,Z32 , IB.ONA, CAP

2SN1,SN2,DU,DE,AK,CFT,SUM,F2,ICYLE, TEMP,ST,C,B, SUB CAP

DIMENSION R(4),Z(4), SK( 100 ,2 ) ,ALM< 100. 2 ) » AL ( 100 ) , AA ( 100 )

,

CAP

1AB( 100) ,AD( 100) ,ADB( 100) ,BK0(500) ,BK1 (500) ,BI 0(500) ,BI1(500)

,

CAP

2CA( 100,2),CB(100,2),SA(100,2),SB(100,2).X1(2),X2(2) ,RATI0(2) ,Y2(3)CAP

3.SN1 ( 100,2 ) ,SN2( 100,2 ) ,CFT(6) ,SUM(7) ,SUB(2) ,ST (6) AP

4 ,C(3,2)

TYPE COMPLEX Ul, ALM,BB,BA, AL .AA.AB , AD, ADB,BK0,BK1 ,BI ,BI 1 ,CA

,

CAP

1CB,SA,SB,X1,X2, X, UM,E1,EM, RAT 10 , Y 1 AREA , Y2

,

CAP

2Y3, ONA,DU,DE, SUM, TEMP, ST ,SUB CAP

3.CCFT

Jl= J COEF

JB1=J+100 COEF

J2 =J+200 COEF

JB2=J+300 COEF

CALL BESSL2(X1( 1) ,J1) COEF

CALL BESSL2(X1(2) ,JB1) COEF

CALL BESSL2(X2( 1) ,J2 ) COEF

CALL BESSL2(X2(2) ,JB2 ) COEF

TEMP= ALM( J,l )/ALM(J,2

)

COEF

IF( J-100)5,7,7

7 TEMP=EM*TEMP

CB(100,1)=<1.,0. ) GO TO 6

5 CA( J,2)=(BK0( JB2)*BK1

(

J2J-BK1 ( JB2 ) *TEMP*BK0 ( J2 ) )*X2<2) COEF

CB( J,2)=(BI0( JB2)*BK1( J2)+BI1( JB2 ) *TEMP*BK0 < J2) )*X2(2

)

COEF

CB(J,1)= ( BI0(J1)*(CA( J,2)*BI1 ( JB1 )+CB( J,2)*BK1( JB1 )

)

COEF

1 -BI 1( Jl )*<CA( J»2)*BI0( JB1 )-CB< J,2)*BK0< JB1 ) )/TEMP)*Xl( 1) COEF

6 SA( J.2)=(BK0( JB1)*BI1 ( Jl >+BKl ( JB1 ) *TEMP*BI 0( Jl ) )*X1 ( 2 ) /CB ( J, 1

)

COEF

SB( J»2) = (BI0( JB1)*BI1 ( Jll-BIK JB1)*TEMP*BI0( Jl ) )*X1(2 )/CB( J,l ) COEF

IF(J-l) 1,1,2 COEF

1 ST( 1 )=(SA( J,2)*BI 1( JB2 )+SB( J, 2 )*BK1( JB2 ) ) *X2(1) COEF

ST(2)=(SA( J,2)*BI0( JB2)-SB( J,2)*BK0( JB2 ) )*X2(1 ) /TEMP COEF

SA(J,1)= BK0( J2)*ST(1 )+BKK J2)*ST( 2) COEF

SB(J»1)= BI0( J2)*ST(1)-BI 1(J2)*ST(2) COEF

3 RETURN COEF

2 IF( J-100)3,l,3

END COEF

SUBROUTINE ADMIT(E2,U2)

COMMON R,Z,F,U1, SK , ALM, BB , BA , AL , AA , AB, AD , ADB , BKO , BK1 , B I , B I 1

,

CAP

1CA,CB.SA,SB,X1,X2, Z21 ,Z31 ,Z41 ,UM , El , EM, RAT 10, Yl ,Z32 , IB,ONA, CAP

2SN1,SN2,DU,DE,AK,CFT,SUM,F2,ICYLE, TEMP, ST, CB, SUB ,IT CAP

DIMENSION R(4) ,Z(4) , SK ( 100 ,2 ) ,ALM( 100,2 ) , AL ( 100 ) , AA( 100 )

,

CAP

1AB( 100 ) ,AD( 100) ,ADB( 100) , BKO (500) ,BKl (500) ,BI0(500).BI1(5 00). CAP

2CA( 100,2 ) ,CB( 100,2) ,SA( 100,2) ,SB( 100,2) ,X1 (2) ,X2( 2) ,RATIO(2) ,Y2(3)CAP

3.SNK 100,2) ,SN2( 100,2 ) ,CFT(6 ) ,SUM( 7 ) ,SUB(2 ) ,ST(6) AP

23

M
4 »C(3*2)

TYPE COMPLEX Ul, ALM,BB»BA» AL AA , AB » AD, ADB ,BK0 »BK1 . BIO »BI 1 »CA , CAP

1CB»SA.SB,X1.X2* X. UM»E1,EM» RAT 10 tYl » ARE_A»Y2, CAP

2Y3. ONA»DU,DE, SUM, TEMP, ST ,SUB CAP

3.CCFT ,DE1 ,DE2,AAA,DE3,AAB

4,DE4,DE5,DE6»E2»U2,E3,AAC,AAD »AL4,Al_5 ,AL8

TYPE REAL MSI ,MS2 ,M IC1 ,MI C2

IF( IT)58,63,63

58 Z6=SQRTF(Z41)

63 DU=UM*EM-U1*E1

DE=E1*U1-1. $DE6=U2*E2 $DE5=E2/E1-1.$DE3=CSQRT ( DE6

)

DE1 = E2*(1./EM-1./E1 ) $DE2 = -F2* ( DE1*U2-UM+U1 ) $DE4= (U1-DE6/E1 )*F2

DE4=DE4*Z6 $DE5=DE5*Z6

C FO CASE HAS 1=100

C BAR STUFF IN EQUATIONS IS LABLED BY TWO ADMIT

SK(1,1)=0. ADMIT

IF DIVIDE CHECK 101,101

101 CALL SAVE(l)

AD(1)= DU*SN2(1,1) +DE ADMIT

ALM(1,1)= (0.,1.)*F $ALM( 100,1)=ALM( 1,1) ADMIT

E3=EM $EM=E2

ALM( 100,2 )=ALM( 100, 1)*DE3

11 AL4=ALM(100,2)**2

ALM(1»2)= CSQRT( F2* ( 1 .+AD( 1 ) )

) * ( 0. , 1 . ) AL8=ALM(1,2)**2-AL4

IF(CABS(AL8/AL4)-1.0E-6)600»600,604

600 PRINT 602

602 FORMAT!* CONDITIONS OF SYSTEM ARE SUCH THAT THE ACCURACY IS */

1* TOTALLY UNRELIABLE *)

CALL EXIT

604 ADB( 1) =DE2*SN2< 1,2) +DE4

ADBU)=ADB< D/AL8

DO 150, J=l, 100,99

Xl(l)= ALM( J,1)*R( 1

)

ADMIT

Xl(2)= ALM( J,2)*R( 1) ADMIT

X2(l)= ALM( J,1)*R(2) ADMIT

X2(2)= ALM( J,2)*R(2) ADMIT

150 CALL COEF (J) ADMIT

EM=E3

AAC= Xlll) *BI0(100)

AAD = -X2(1) *(SA( 100,1 ) *BI ( 3 00 ) -SB( 100,1 )*BK0(300) ) AAA=DE1*AAC $AAB=DE1*AAD

AL(1 >=AAA*SN2< 1,2) +AAC*DE5

AB(1)=AAB*SN2(1,2) +AAD*DE5

IF DIVIDE CHECK 40,200

200 BA= SA(1,1) -1. ADMIT

IF(CABS(BA)-l.E-7)20 1,20 1,22

220 PRINT 42

42 FORMAT(*NORMALIZATION ERROR HAS OCCURRED*)

201 BB= SB(100,1)

ICYLE =0

IA=-1

IF EXPONENT FAULT 499,499

499 CALL INTGR(l)

2U
IF EXPONENT FAULT 500,501

500 RATIO( 1)=RATI0(2 )=(0.,0. ) S GO TO 40

501 BB=BB+ST(5)

SUB(2) =BB

RATI0( 1)=BB/SA(100,1)

IF( IB) 1,1,2

1 DO 100,1=2,99 $ 1=1

CALL SAVE( I

)

ADMIT

ICYLE =ICYLE+1

AD( I )=DU*SN2( 1,1 )+DE

AL( I )=AAA*SN2( 1,2) $AB( I )=AAB*SN2( 1,2)

ALM(I,1)= ( SMI, 2) -F2) ADMIT

ALM(I,2)= ( ALM(I.l) -F2*AD(I)) ADMIT

AL5=ALM( I,2)-AL4 $ALM( I ,2 ) =CSQRT ( ALM { 1,2))

IF <CABS( AL5/AL41-1.0E-6) 600,600,601

601 ADB( I

) =DE2*SN2( I ,2)/AL5

ALM(I,1)= C5QRT( ALM( I ,1) ) $AA ( I

) =SN1 ( I , 2 ) *DU

Xl( 1 )=ALM( I ,1)*R( 1) ADMIT

X1(2)=ALM( I,2)*R( 1) ADMIT

X2( 1 )=ALM( I,1)*R(2) ADMIT

X2(2)=ALM( I ,2)*R(2) ADMIT

CALL COEF( I) ADMIT

IF DIVIDE CHECK 40,300

300 IF EXPONENT FAULT 301,301

301 CALL INTGR(I)

IF EXPONENT FAULT 40,400

400 ST(6) = ST(4)*AA(I) *F2 ADMIT

MS2= CABS(SUB<2) )*(0.001)

MIC2=CABS(ST(6)

)

ADMIT

SUB(2)=SUB(2)+ST(6)

20 IF( MS2-MIC2) 100,100,99

99 IF( IA) 102,103,2

102 14=1+1 $IA = S GO TO 100

103 IF( 1-14)2,2,97

97 14=1+1

100 CONTINUE ADMIT

2 RATIO( 2)=SUB(2)/SA( 100,1

)

RETURN ADMIT

40 PRINT 41

41 FORMAT(*NEXT LINE HAS HAD AN OVERFLOW OR DIVIDE CHECK*)

RETURN

END ADMIT

SUBROUTINE BESSL2(X7,I) BESSL2

COMMON R,Z,F,U1, SK , ALM, BB , BA, AL , AA , AB, AD , ADB. BKO , BK1 ,BI ,B I 1 , CAP

1CA,CB,SA,SB,X1,X2, Z21 ,Z31 ,Z41 ,UM ,E1 , EM, RAT 10, Yl ,Z32 , I B.ONA, CAP

2SN1, SN2,DU,DE,AK,CFT, SUM, F2, ICYLE, TEMP,ST,C,B, SUB CAP

DIMENSION R(4),Z(4), SK ( 100 , 2 ) ,ALM( 100 , 2 ) , AL ( 100 ) , AA ( 100 ) , CAP

1AB( 100 ), AD (100) , ADB (100) , BKO (500) , BK1 ( 500 ) ,B I ( 500) ,BI1 (500) , CAP

2CA(100»2),CB(100,2),SA(100,2),SB(100,2),X1(2),X2(2) ,RATIO(2) ,Y2(3)CAP

3.SNK 100,2 ) ,SN2( 100,2 ) ,CFT(6) , SUM ( 7 ) ,SUB ( 2 ) , ST ( 6

)

AP

4 ,C(3,2)

TYPE COMPLEX Ul, ALM,BB,BA, AL , AA , AB , AD , ADB , BKO , BK1 , B I , BI 1 ,CA , CAP

1CB,SA,SB,X1,X2, X, UM,E1,EM, RAT I , Yl , ARE A , Y2 , CAP

2Y3, ONA,DU,DE, SUM, TEMP, ST ,SUB CAP

25

1^
100

200

222

21

4

3,C»CFT

TYPE DOUBLE X4,X5 ,H01 A ,H0 1B,H02A,H02B ,H1 1A ,H11B H12A,H12B

TYPE COMPLEX X3,PH,X7

DATA (PH=(0. »1. ) ) IF DIVIDE CHECK 2.1

X3=X7*PH

X4 = X3 $ X5=X7

CALL ALLHU4,X5»H01A,H01B,H02A,H02B,H11A,H11B,H12A,H12B)

BIl(l)=lO. »-0.5>*( (HllA+H12A)+PH*(HilB+H12B) ) BK1II ) = (-1.5 70 79633,0. ) * ( HI 1A+PH*H11B ) IF( 1-401) 100,200,100

BIO( I )=0.5*( (H01A+H02A)+PH*(H01B+H02B ) ) BK0< I

) =(0. ,1.57079633 )*( HOI A+PH*H0 IB)

IF DIVIDE CHECK 22,3

J=l $ GO TO 21

J = 2 PRINT 4, I,

J

F0RMAT(* NEXT, LINE HAS BESSL2 DIVIDE CHECK 1= *,2I4)

RETURN

END

SUBROUTINE ALLH ( ZA ,ZB,H01 A ,H01B ,H02A ,H02B,H1 1A , HI IB ,H12A ,H12B ) THE REMAINING SUBROUTINES ARE A PACKAGE FOR HANKEL FUNCTION

GENERATION . THIS PACKAGE WAS DEVELOPED BY GARNEY HARDY

TYPE DOUBLE ZA ,ZB , HOlB ,H0 1A ,H02 A,H02B ,H1 1 A ,H1 1 B ,H12A

2,H12B,X,Y»ABSZ

X=ZA

Y = ZB

ABSZ = DSQRT <X*X + Y*Y)

IFIABSZ-7.7D) 1 ,1,3

1 CALL HSMALL(X,Y,H01A,H01B,H02A,H02B,H11A,H11B,H12A,H12B)

RETURN

3 CALL HLARGE(X,Y,H01A,H01B,H02A,H02B,H11A,H11B,H12A,H12B)

RETURN

END

SUBROUTINE CHISI CHI ) DIMENSION CHI (20)

TYPE DOUBLE CHI

CHI ( 1

) =1.57D

CHI (2)=2.D

CHI (3) =2.36D

CHI (4) =2.67D

CHI ( 5)=2.95D

CHI (6) =3. 2D

CHI (7)=3.44D

CHI (8)=3.66D

CHI (9)=3.87D

CHI'( 10)=4.06D

CHI ( 11 )=4.26D

CHI ( 12 ) =4.43D

CHI ( 13 )=4.61D

CHI ( 14)=4.77D

CHI ( 15 1=4. 94D

CHI ( 16)=5.09D

CHI ( 17) =5.25D

BESSL2

BESSL2
CHI (18)=5.39D

CHI ( 19)=5.54D

CHI (20)=5.68D

RETURN

END

SUBROUTINE AS(A0,A1)

DIMENSION AO(20) ,A1(20)

TYPE DOUBLE A0,A1

TYPE DOUBLE AN

A0(1 )=-.125D

A1(1)=.375D

DO 1 N=2,2U

AN=N

AO(N)=-( ( (2.D*AN-1.D)**2)*A0(N-1) )/(AN*8.D)

1 A1(N)=( (4.D0-( (2.D0*AN-1.D0)**2) )*A1(N-1) )/( AN*8.D0)

RETURN

END

SUBROUTINE J SMALL ( X,Y AJ0A,AJ0B,AJ1A,AJ1B)

TYPE DOUBLE X,Y» AJOA,AJOB ,AJ1A»AJ1B*ZA,ZB, AM ,ZMA»

1ZMB .SIGN . SUMA » SUMB , FACTM » TERMA , T ERMB

ZA = X ZB=Y

ZA=ZA/2.D

ZB=ZB/2.D

CALL ZMPY(ZA,ZB»ZA,ZB.ZA,ZB)

AM=0.D

FACTM=1.D

ZMA=1.D

ZMB=O.D

SUMA=1.D

SUMB=O.D

SIGN=1.D

1 CALL ZMPY(ZA,ZB,ZMA,ZMB»ZMA,ZMB)

SIGN=-SIGN

AM=AM+1.D

FACTM=FACTM*AM

TERMA=ZMA/ ( FACTM*FACTM

)

TERMB=ZMB/(FACTM*FACTM)

SUMA=SUMA+SIGN*TERMA

SUMB=SUMB+SIGN*TERMB

2 IF(DABS( TERMA )+DABS ( TERMB ) -l.D-13 ) 3 » 1 »1

3 AJOA=SUMA

AJOB=SUMB

AM=1.D

FACTM=1.D

ZMA=1.D

ZMB=O.D

SUMA=1.D

SUMB=O.D

SIGN=1.D

5 CALL ZMPY(ZA,ZB.ZMA,ZMB»ZMA*ZMB)

SIGN=-SIGN

FACTM=FACTM*AM

AM=AM+1.D

27
TERMA=ZMA/ ( FACTM*FACTM*AM

)

TERMB=ZMB/ ( FACTM*FACTM*AM

)

SUMA=SUMA+SIGN*TERMA

SUMB=SUMB+SIGN*TERMB

4 IF(DABS(TERMA)+DABS(TERMB)-1.D-13)7,5,5

7 AJ1A=SUMA/2.D

AJ1B=SUMB/2.D

CALL ZMPY(X»Y,AJ1A,AJ1B»AJ1A»AJ1B)

RETURN

END

SUBROUTINE HSMALL ( X , Y » H01A.H01B ,H02A,H02B ,H11 A ,H11B ,H12A ,H12B

)

TYPE DOUBLE X ,Y , HOIA ,HO IB ,H02A,H02B ,H11A,H11B ,H12A,

1H12B,ZA,ZB,AJ0A»AJ0B,AJ1A,AJ1B,PSI»SIX,SUMA,

2SUMB»S» FACTS »TERMA,TERMB»ALOG,BLOG.FUNNYA»

3FUNNYB»ZMA,ZMB,SIGN

ZA = X ZB = Y CALL JSMALL(ZA,ZB»AJ0A,AJ0B»AJ1A.AJ1B)

PSI=~. 5772 1566490 1532 8606065 12

D

S I X=. 6 3661 9772367 58 1343075 53 50D

ZA=ZA/2.D

ZB=ZB/2.D

CALL ZMPY(ZA,ZB»ZA,ZB»ZA,ZB)

ZMA=1.D

ZMB=O.D

SUMA=PSI

SUMB=O.D

S = O.D

FACTS=1.D

SIGN=1.D

1 CALL ZMPY(ZA,ZB»ZMA,ZMB,ZMA,ZMB)

SIGN=-SIGN

S=S+1.D

FACTS=FACTS*S*S

PSI=PSI+1.D/S

TERMA=PSI*ZMA/FACTS

TERMB=PSI*ZMB/FACTS

SUMA=SUMA+SIGN*TERMA

SUMB=SUMB+SIGN*TERMB

IF(DABS(TERMA)+DABS(TERMB)-1.D-13) 3.1,1

3 H01A=AJ0A+SUMB*SIX

H01B=AJ0B-SUMA*SIX

H02A=AJ0A-SUMB*SIX

H02B=AJ0B+SUMA*SIX

PS I =-.5772 156649 01 532 8606065 12

D

ZMA=1.D

ZMB=O.D

SUMA=PSI+PSI+1.D

PSI=SUMA

SUMB=O.D

5 = 0.

D

FACTS=1.D

SIGN=1.D

5 CALL ZMPY(ZA,ZB,ZMA,ZMB»ZMA,ZMB)

28
SIGN=-SIGN

S=S+1.D

FACTS=FACTS*S*(S+1.D)

PSI=PSI+1.D/S+1«D/(S+1.D)

TERMA=PSI*ZMA/FACTS

TERMB=PSI*ZMB/FACTS

SUMA=SUMA+SIGN*TERMA

SUMB=SUMB+SIGN*TERMB

IF(DABS<TERMA)+DABS(TERMB)-1.D-13)7,5.5

' CALL ZMPY(ZA.ZB»SUMA»SUMB,SUMA,SUMB)

SUMA=SUMA+1.D

CALL ZMPY(X,Y,AJ1A,AJ1B,AJ1A,AJ1B)

H11A=AJ1A+SIX*SUMB

H11B=AJ1B-SIX*SUMA

H12A=AJ1A-SIX*SUMB

H12B=AJ1B+SIX*SUMA

ZA=X/2.D

ZB=Y/2.D

ALOG=DSQRT ( ZA*ZA+ZB*ZB

)

ALOG=DLOG(ALOG)

BL0G=DATAN2(Y.X)

FUNNYA=-SIX*BLOG

FUNNYB=SIX*ALOG

CALL ZDIV(AJ1A,AJ1B,X,Y.AJ1A,AJ1B)

CALL ZMPY(FUNNYA,FUNNYB,AjOA,AJOB,AJOA,AJOB)

CALL ZMPY(FUNNYA»FUNNYB»AJ1A»AJ1B.AJ1A,AJ1B)

H01A=H01A+AJ0A

H01B=H01B+AJ0B

H02A=H02A-AJ0A

H02B=H02B-AJ0B

CALL ZDIV(H11A,H11B»X.Y»H11A»H11B)

CALL ZDIV(H12A,H12B,X,Y»H12A,H12B)

H11A=H11A+AJ1A

H11B=H11B+AJ1B

H12A=H12A-AJ1A

H12B=H12B-AJ1B

RETURN

END

SUBROUTINE HLARGE ( X ,Y ,H01A .H01B H02A,H02B,H11A,H11B,H12A ,H12B

)

TYPE DOUBLE X , Y » H01A*H01B ,H02A »H02B .HI 1A.H1 1B.H12A.

1H12B»ZA,ZB»RIPI »PI , ABSZ .ARGZ » ABZM » FACTO , FACT 1 ,SUM01 A

,

2SUM02A»SUM01B,SUM02B.SUM11A»SUM11B.SUM12A»SUM12B»

3CXA,CXB,CXNEGA,CXNEGB,ZSA,ZSB,EPS»EPA»PA,PB,QA,QB»

4ZHALA,ZHALB,EY

DIMENSION CHI (20) A0(20).A1(20)

TYPE DOUBLE A0»A1,CHI

ZA = X ZB = Y IF(SENSE LIGHT 1) 1,3

L RIPI=. 5641895835477562 869480 795D

PI =3. 14159265358 9 7932 3 846 2643D

CALL CHISJCHI ) CALL AS(AO.Al)

ABSZ = DSQR T ( ZA*ZA+ZB*ZB

)

29
ARGZ=DATAN2(ZB,ZA)

IF ( DABS (ARGZ) -PI /2.D) 29*29*31

31 IF(ARGZ-PI/2.D)33. 29*35

35 ZA=-ZA

ARGZ=PI-ARGZ

M=l

GO TO 37

33 ZA=-ZA

ZB=-ZB

ARGZ=ARGZ+PI

M = 2 GO TO 37

29 M = 3 37 ABZM=1.D

FACT0=2.D*DEXP(PI/(8.D*ABSZ)

)

FACT 1 = 2. D*DEXP< (3.D*PI )/(8.D*ABSZ)

)

SUM01A=1.D

SUM01B=0.D

SUM02A=1.D

SUM02B=0.D

SUM11A=1.D

SUM11B=0.D

SUM12A=1.D

SUM12B=0.D

CXA=1.D

CXB=O.D

CXNEGA=1.D

CXNEGB=O.D

ZSA=1.D

13 ZSB=O.D

EPS=1.D

DO 5 K=l*20

EPA=EPS

ABZM=ABZM*ABSZ

EPS = FACTO*CHI ( K. ) *DABS ( AO < K. ) ) /ABZM

IF(EPS-EPA)19,7,7

19 CALL ZMPY(ZSA,ZSB.ZA*ZB*ZSA»ZSB)

CALL ZMPY(CXA,CXB,O.D,l.D,CXA,CXB)

CALL ZMPY(CXNEGA,CXNEGB»0.D»-1.D.CXNEGA,CXNEGB)

PA=CXA*AO(K)

15 PB=CXB*AO(K)

CALL ZDIV(PA,PB*ZSA,ZSB»QA.QB)

SUM01A=SUM01A+QA

SUM01B=SUM01B+QB

PA=CXNEGA*AO<<

)

17 PB=CXNEGB*AO(K)

CALL ZDIV(PA,PB*ZSA,ZSB.QA»QB)

SUM02A=SUM02A+QA

SUM02B=SUM02B+QB

18 IF(EPS-1.D-10)7,7,5

5 CONTINUE

7 CALL ZMPY(RIPI ,-RIPl ,SUMO 1A.SUM01B »H01 A »H01B

)

CALL ZMPY(RIPI.RIPI,SUM02A»SUM02B»H02A*H02B)

CXA=1.D

30
CXB=O.D

CXNEGA=1.D

CXNEGB=O.D

ZSA=1.D

ZSB=O.D

ABZM=1.D

EPS=1.D

DO 9 K=l»20

EPA=EPS

ABZM=ABZM*ABSZ

EPS=FACT1*CHI (K)*DABS( A1(K) ) /ABZM

IF(EPS-EPA)25,11,11

25 CALL ZMPY(CXA,CXB,O.D»l.D,CXA,CXB)

CALL ZMPY(CXNEGA»CXNEGB.O.D»-l.D»CXNEGA,CXNEGB)

CALL ZMPY(ZSA,ZSB,ZA,ZB,ZSA,ZSB)

PA=CXA*A1 IK.) 21 PB=CXB*A1(K)

CALL ZDIV(PA,PB.ZSA,ZSB»QA»QB)

SUM11A=SUM11A+QA

SUM11B=SUM11B+QB

PA=CXNEGA*A1 (K)

23 PB=CXNEGB*A1 (K)

CALL ZDIV(PA,PB,ZSA,ZSB,QA,QB)

SUM12A=SUM12A+QA

SUM12B=SUM12B+QB

24 IF(EPS-1.D-10)11,H,9

9 CONTINUE

11 CALL ZMPY(-RIPI -RIPI ,SUM11A,SUM1 IB ,H1 1A.H118

)

CALL ZMPY(-RIPI ,RIPI ,SUM12A .SUM12B .H12A >H12B

)

ZHALA = DSQRT( ABSZ)*DCOS( .5D*ARGZ ) ZHALB=DSQRT( ABSZ )*DSIN( .5D*ARGZ)

EY=DEXP(ZB)

PA=DCOS(ZA)/EY

PB=DSIN <ZA)/EY

CALL ZDIV(PA,PB,ZHALA,ZHALB,QA,QB)

CALL ZMPY(H01A,H01B,QA,QB,H01A,H01B)

CALL ZMPY(H11A.H11B,QA,QB»H11A»H11B)

CALL ZMPY(PA,PB,ZHALA,ZHALB,QA,QB)

CALL ZDIV(H02A,H02B,QA,QB,H0 2A,H02B)

CALL ZDIV(H12A,H12B,QA,QB,H12A,H12B)

GO TO (39»41,43) »M 39 H02A=2.D*H01A+H02A

H02B=-2.D*H01B-H02B

H12A=-2.D*H11A-H12A

H12B=2.D*H11B+H12B

H01A=-H01A

H11B=-H11B

RETURN

41 PA=H01A

PB=H01B

QA=H11A

QB=H11B

H01A=2.D*H01A+H0 2A

H01B=2.D*H01B+H02B

31
H11A=-2.D*H11A-H12A

H11B=-2.D*H11B-H12B

H02A=-PA

H02B=-PB

H12A=QA

H12B=QB

43 RETURN

END

32
t>

•*

(VI (VI >- • X a rl a i < 1-1 * o X O 3 O

u> 1pH

_i co m r- o m -* •

>-

H (VI a oo >* 1 © o o

in I in o (VI in o «-4 • o ro

>o (VI 00 o • o • o«/> X >- H on a o uj «t •-• • -1 •» to oz i- © 1 o <-> o •-• Ul 1 © cvi en X r-« o o Z x o O O UJ o ro 1 X <_> x0 z O O t-t r-l • o o in q o <•>

t—i >*f-

© >-

u

H in rvi o ui

<;

-i • • o <r •* (VIO + < oo o o o o X s OH- O 1

1—1 Q

u Otfl

r -t

<!

O Ul Z o —

z • <r a roro H

•—

< voao

Oh

U UJ (VI ao • i (J >- «rt W O (VI 1— m m

<c

J o o i— t-i a, o o < * ao 1 • 1- o

vo o in GO

O <

-*

m o CO 1 I- I H4 !-• xfj (VI (VI CO X CO 00 r4 • • »-t — <*> z 3 < m 3 ~>* O ro • ^>-»ro<

UJ UJ >-

X o ox * < o < om a. o

10 • o o o >- o ui » a ONOU II o »— o o o z ~ o

Ul o o o < roo> x 1 • h- * — o o >0 © O 1— Z X X • UJ MO |h O UJ -4

CD • -•in

x -I

X I M w-t in a l-t >-• 1— OX • • < t/i Min

•i >o *i ,o

X 2 o O O • • Ul O x o X Ul + o u (S O

t-

©

© z <o o O Z <X©

x o Ol-l xrvi

\— o * Ul •

*-l • o Ul

X r-. 2 -H ©X Ul Ul o <

UJ -1< X tr x © z^ o < < • CO O O »- Nl- »-l o -• 00 CO z t— u Z Z UJ o o < O O X «© a. •—1 »-H UJ x o < VI CO X « ii.au O C/> Z Z 3 x •

»H Ul UJ 10 3 uiin

X XX < <_> i O V i-t »—i UJ 3 -• Z IlUQQWS x «tro z HI 1— o

uj ui _i to

x 1— o

z-g

«t i-1

—1 1— *-t t-t X Ul O X

O < Ul o <t z in »-X < X Ul «s O

h-U.

X

M Ul -I r- »H (J t— • 2$ UJ UJ 4 UJ X (VI X X X X X Ul « r- 1- t— «/» r- 3 Z

33

(iK

APPENDIX d

FLOW CHART OF MACHINE PROGRAM

BEGIN

READ j CARDS

© r , , r

2

, r

3 ® Z,,Z 2 ,Z S ,Z 4 IA , IB © € l»Mi

READ | CARD

© FREQ. , IC, IS, ITES

€ M « MM

PRINT OUT RESULTS

RADIATION

TEST

<D

<D

CALCULATE

IB =

CALCULATE H

(2)

-I

EXIT 35 BPO nMII


U.S. DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20230

POSTAGE AND FEES PAID

U.S. DEPARTMENT OF COMMERC

OFFICIAL BUSINESS