AiiiQM saaaib

NBSIR 86-3362

SCAT: A Vector Program to

Solve A Transient MFIE

PUBLICATIONS

Egon Marx

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratories

Center for Manufacturing Engineering

Gaithersburg, MD 20899

April 1986

U.S. DEPARTMENT OF COMMERCE

QC

100

.056

86-3362

1986

C. 2

BUREAU OF STANDARDS

tres KSEAP.CH K^FOiy.IAT!0:« CENTE?.

NBSIR 86-3362

9 « « SCAT: A VECTOR PROGRAM TO

SOLVE A TRANSIENT MFIE

Egon Marx

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratories

Center for Manufacturing Engineering

Gaithersburg, MD 20899

April 1986

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

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director

Abstract

The FORTRAN program SCAT is used to solve the magnetic -field

integral equation (MFIE) to determine the -fields scattered by a

perfectly conducting sphere. The incident field is a plane-wave

pulse, and a steppi ng-i n-ti me procedure is used to determirie the

surface current density induced on the sphere. The program does

not take advantage of the special symmetry of the scatterer

because it is intended to serve as a verified starting point for

more general programs. The output is compared to that of the

program PERF, which computes the same fields via a Fourier

transform of the monochromatic fields obtained from the Mie

formulas. The contributions of the self-patch and neighboring

patches to the singular integral are optionally computed by using

their expansions in the linear size of the patches. The selfpatch term is important for the solution of other integral

equations that may be of the first kind. For the MFIE, these

corrections are small but not negligible. The program takes

advantage of the vector programming features of the CYBER 205.

Key words: computer programs; el ectromagnet i c pulse scattering;

magnetic field integral equation; Mie scattering;

perfectly conducting sphere; stepping-in-time

procedure; transient el ectromagnet i c fields; vector

pr ogrammi ng

.

1
sf Qsq&sq&s

Page

1. Introduction ........ ..... 1

2. Program PERF ....... ....... 5

3. Program SCAT 8

3. 1 The Main Loop 9

3.2 The Interpolation Formulas 9

3.3 Subroutine RUNT 12

3.4 Subroutine INFL 13

3.5 Subroutine FARFL 14

4. Self-Patch Corrections .... 15

5. Neighboring—Patch Corrections 17

6. Vector Programming Considerations 23

7. Ancillary Programs 25

7. 1 Program PERFPLT 25

7.2 Program CURCHK 26

7.3 Program SCATP 26

7.4 Subroutine VINUFT 27

8. Program Inputs and Outputs 28

9. Conclusions 38

References ..... 39

Append!

x

41
i • iDtCQdyct i_on

There are many procedures related to the metrology mission

of the National Bureau of Standards that require the analysis of

waves scattered by an obstacle. Comparison of the scattered

waves with those computed for known obstacles can serve to

determine parameters related to the scatterer. Codes used in

inverse scattering problems also may require the solution of the

direct problem. If the scatterer is a homogeneous body, the

scattered fields can be computed from tangential fields defined

on the surface of the scatterer. In turn, these fields can be

computed by solving a singular integral equation.

The el ectromagneti c fields scattered by a perfect conductor

* can be obtained form the surface current density J CID . This

s

tangential field can be determined by solving the magnetic field

integral equation (MFIE)

,

-+ i n

J (x,t) = (2/M,-)n X B (x,t)

^ n X 4) dS^R

2tt

1 S 1 (1)

-in where B is the magnetic field of the incident pulse, Mq is the

permeability of free space, n is the outside normal to the

Numbers in square brackets indicate the literature references

at the end of his report. There are further references in these

publ i cat i one.

1
scatterer, t is the retarded time

T = t - R/c, ( 2 )

and

R = X - X' , R = IRI . (3)

This is a singular integral equation because the integrand

diverges Mhere R =0, and it is an integral equation o+ the

second kind because the unknown -field also appears outside the

integral, on the le-ft side of the equation.

An example of a singular integral equation of the first kind

is the electric field integral equation (EFIE)

,

^in - 1

" 0 = nxE (x,t) + rr-f 4n X dS'

S

(x ' , T

)

+ (x ' , T ) R'^

®

Mq ,T)'

R aV (4)

1 n where E is the electric field of the incident pulse, is the

permittivity of free space, and is the surface charge density,

related to by conservation of charge.

We have shown C2, 3D how the scattered and transmitted

fields of a pulse scattered off a homogeneous dielectric body can

be determined from a single tangential vector field that obeys a singular integral equation of the first kind. We have further

2
generalized C41 this -formulation to the scattering by a

dispersive body, such as an imper-fect conductor. In some papers

we also derive equations -for the scattering of other types of

fields by homogeneous bodies, as well as the correspondi ng

equations for the scattering of monochromat i c waves.

The solution of the integral equations by means of the

stepping-in-time procedure takes advantage of causality and the

finite speed of propagation of the wave. The new values of the

unknown field are expressed in terms of an integral over the

retarded values of the field, which have been calculated in

previous steps.

The integration can be performed numerically by dividing the

surface of the scatterer into patches and adding the products of

the integrand at the center of each patch times the area of the

patch. This prescription does not apply to the contribution of

the self -patch because this is the patch where the integrand

becomes infinite as the field point x and the source point x'

coincide. For an integral equation of the second kind, the

contribution of the self—patch can be neglected in the

computation of a new value of the field in comparison to the

large number of contr i but i ons from the other patches. This is

not so for an integral equation of the first kind, because the

self-patch contribution to the integral is the only term that

depends on the new value of the field.

We have expanded the contribution of the self-patch to first

order in the linear size of the patch for an arbitrary surface

[51, and in this form it may be included in the integral. We
have similarly expanded the contribution o-f a neighboring patch,

because the integrand cannot be well approximated by its value at

2 the center o-f such a patch when 1/R or 1/R varies significantly

over the patch C63.

Me have developed the program SCAT to implement the

computation of the surface current density and the far fields

scattered by a perfectly conducting sphere by solving the MFIE.

In this program we provide the option of computing the correction

to the surface fields due to the self—patch contribution and the

more accurate computation of the nei ghboring-patch contributions.

We found that terms that contain the time derivatives of the

fields, terms that are negligible in most computations, cannot be

ignored when the fields themselves are small C73. The expansions

-4 of the surface current density J and its derivatives that we use

s in SCAT are explained in C73; they differ from those we used in

£51 and £61 in that the dependence of the retarded time r on the

curvilinear coordinates u and v is taken into account when a partial derivative with respect to u or v is taken.

We intend SCAT to serve as a validated starting point for

other programs that will determine the fields scattered by more

general types of scatterers. For this reason we do not take

advantage of special symmetry properties of the sphere. To

validate the program, we compare the far fields to those computed

by means of another program, PERF, which uses the Mie formulas

£71 for monochromati c fields and a Fourier transform. This

program is discussed in section 2.

The general structure of SCAT is presented in section 3.

Section 4 is devoted to a discussion of the issues related to the

4
computation of the contri buti ons of the patches neighboring the

self-patch, and in section 5 we examine the contribution of the

self-patch itself. In section 6 we point out where we use vector

programming as appropriate for the CYBER 205. We describe a

number of ancillary programs in section 7, and we show some

examples of results in section 8. Listings of computer programs

are given in the appendix.

We define the incident plane-wave pulse so that it reaches

the sphere at t = 0. We choose a double-exponential smoothed at

the beginning by a power factor for the profile of the plane

wave. The incident magnetic field is

= O, <5>

where 0 is the unit step function, a is the rise constant of the

pulse, ^ the decay constant, and

f=z+a— ct, (6)

where a is the radius of the sphere.

2. RC99C§[D EERE

The Mie formulas give the fields scattered by a homogeneous

sphere from an incident plane monochromatic wave. For a perfect

conductor and a wave incident along the z—axis, the scattered far

field components reduce to
E exp ( 1 kr

)

ikr

ac< 51 2

2n

n (n

n=l

+ 1 + 1 )

P

^ ^

^ (cos©

)

nsin©

+ b n

dP ^ (cos©)

n d© }-1 ) n+1

(7)

= exp ( i kr

)

ikr

'H

, Y 2n + 1 / cos«> > — —

-<a

Z. n <n + 1 ) In

n=l ^

dP ( 1 ) n (cos©

)

d©

-f- P^^^ (COS0)\

b U-l)

n sin© I

n+1

»

( 8 )

where k is the wavenumber, r is the distance from the center of

the sphere to point where the field is computed, © is the polar

or scattering angle, 0 is the azimuthal or polarization angle,

P^^^ is the associated Legendre function, and a and b are

coefficients given by

a (p) = - j (/!)) /h (/>: , (9) n -*n n

b (p) = - Cpj (p) D D'

,

n n n

( 10 )

6
where j is a spherical Bessel function, h^^^ is a spherical

n n Hankel function of the first kind, and

p = ka = 2na/X. <11)

The primes in eq. (11) indicate derivatives with respect to p.

The Fourier transform of the incident wave (5) is

= expCikCz + a) 1 (^/c) C

ik)“^ - (a - ik)“^l, (12) w

where k = oj/c. The coefficient of exp(ikz) multiplies the

amplitude of the far fields given in eqs. (7) and (8). The

steady state component of the incident pulse does not contribute

to the scattered fields.

In the program PERF, we compute first the required Bessel

and Hankel functions by means of the subroutines BESRJ and BESRY

l 83, then the scattering coefficients a and b in the subroutine

n n SCACO, and the monochromati c fields in the subroutine FLD. The

inverse Fourier transform is carried out in subroutine INUFT C9]

or VINUFT (its vector version); we do not use a fast Fourier

transform because we determined that it was better to cover low

and high frequencies at different intervals than to use an

equispaced grid. The frequencies are selected in subroutine

OMEGA C91 and the intervals increase in a geometric progression,

which gives a thorough coverage of low frequencies and a sparse

but extended coverage of high frequencies. The input to the

program requires the specification of a, ff, and the radius of the

sphere, the number of scattering angles and the scattering

angles. In addition. we have to specify the lower and upper

7
^^requency limits at which the fields are computed and the number

of frequency and time points. The output consists of the

intensity of the scattered fields for both polarizations of the

scattered fields; they are saved in a file and they are plotted

using the subroutines DRAW3 C103 or DRAW4 till.

3. SCAT

In this program we compute the surface current density

induced on the sphere by solving the MFIE (1) in the time domain

by a steppi ng—i n-ti me procedure. The scattered fields are then

obtained by integration, and we compute the intensities in the

far zone for comparison with the output of PERF. The magnetic

field in the radiation zone is

(X, t) 4nr r X U dS

aj (x',T)

s

at^

(13)

The surface current density is tangential to the surface and

has only two independent components. To decrease the amount of

memory required to run the program, we decided to store only the

components along the meridians and parallels on the sphere. To

sum the contributions from different patches, we first compute

the components of the vectors in a Cartesian coordinate system.

Once the two tangential components of the surface current density

are determined, they are saved on disk for future use.

8
3^- i It!® L9°E

Most o-f the execution time needed -for the program is spent

computing the surface current density by solving the MFIE. The

components of this tangential vector are computed at the center

of each patch starting at t = O and after each time increment At.

If we ignore the self-patch contribution to the integral in

eq. (1), the contributions from the other patches can be computed

because the fields at their centers have previously been

calculated at the retarded times. We have to choose a time

interval that is smaller than the distance between any pair of

patches (divided by the speed of light c) to obtain a retardation

that is sufficiently large. It is unlikely that the current

density will have been calculated at precisely the retarded time,

and an interpol ation is needed to obtain the necessary values.

The i nterpol at i on is discussed in the next subsection.

After the values of the surface current density are computed

in this way, we can perform the corrections for the self—patch,

for the neighboring patches, or for both; these corrections are

described in detail in separate sections.

3.2 The Interpol at i_gn Formulas

The interpolation method C123 turned out to affect

significantly the outcome of the calculations. The program SCAT

allows for two alternatives. In the first one, the fields at the

retarded times are computed from a polynomial of degree 4 that

passes through five points chosen so that they fall on both sides

9
o-f the retarded time. The relation has the -form

y(t) = a^t"^ + + a^t + ( 14 )

and we rescale the independent variable so that it takes integer

values from -2 to 2. In terms of the correspond! ng values of the

dependent variable, the coefficients are

^4 “ - iy(l) + iy< 0 ) - iy(-D + 54>'<-2>. (15)

a^ = ^y(2) -

12' + zy(-D - ^^(-2), 12 ’ ( 16)

a„ = - 24 J (2) + |y(l) - zy<0) + |y(-l) - irv<-2). 24' (17)

a^ = - Y^<2) + |y(l) - §y<-l) + (18)

a^.j = y (O) . (19)

In the second method, a polynomial of degree 2 is fitted to the

same five points by a 1 east-squares approximation. The

coefficients a. and a_ are zero, and the other coefficients are

4 3

a„ = iy(2) - Tr/<1> " kv <0) - + ^y(~2), 14’ 14

( 20 )

10
= iy<2) ( 1 ) - ( 21 )

a

O 35

12

y (-1) ( 22 )

The polynomial interpolation gives more accurate results, but

leads more easily to numerical instabilities E7, 13]. The time

derivatives of the fields are calculated from

yMt) = 4a^t" + 3a^t^ + 2a^t + a^. (23)

The values of the fields and their derivatives are set to

zero for times before the time the incident pulse reaches the

center of the correspond! ng patch on the sphere. The

interpolated values do not necessarily vanish at these times and

the error increases when they are used.

The interpolation formulas are implemented in a form

suitable for the execution of vector instructions in subroutines

INTRPO and INTRPl, and INTRP2 and INTRP3 are scalar versions

needed only for the computation of the derivatives of the fields

at the self—patch. The choice of the five points for the

interpolation also has an effect on the accuracy of the results,

which is improved if the retarded time is close to the center of

the interval. We choose the points in this way except for

patches that are so close to the self—patch that the retardation

is too small for all the needed values to be available. The

values of the components of the surface current density and of

their time derivatives at the self-patch itself are needed to

1

1
carry out the sel-f- and neighboring-patch corrections; we use the

uncorrected values at the self-patch as a first approximation

where necessary.

3.3 Subroutine RUNT

This subroutine is executed once at the beginning of the

program to perform the division of the sphere into patches and to

compute coefficients that will be needed during the execution of

the program.

We do not take advantage of the symmetry of the sphere to

simplify the integral equations, but we have to take into account

the character i st i cs of the sphere to choose a good way to divide

it into patches. We naturally use spherical coordinates, the

polar angle 0 and the azimuthal angle which form a principal

orthogonal curvilinear coordinate system C51.

We divide the surface of the sphere into patches by

meridians and parallels. We first divide the surface into n^

strips of equal increments of 0. Each strip is then divided into

a number of patches that is roughly proportional to the length of

the strip, that is, to sin0. We leave two undivided caps at the

poles, where spherical coordinates become undefined; we use

Cartesian coordinates for these caps, which requires special

formulas for the computations. We divide the strips next to the

caps into a given number n^ of patches, usually 5. Then the

other strips are divided into approx i matel

y

n = (n_ - n^)sin0 + n_ (24) 2 •:> 3

12
patches o+ equal size by lines a-f constant

This more or less asymmetrical division of the surface into

patches leads to computed nonzero values for some intermediate

variables that should vanish by symmetry. Although this

discrepancy does not seem to make much difference in the final

results, we get intermediate values that are closer to zero by

restricting the number of patches in each strip to a power of 2.

For strips near the poles, n is not large and the values of

the increment in = 2n/n, may not be small compared to 1.

The side of the patch, asinOA^, is small compared to a. To obtain patches that are approx i matel y square in form, we

would have to choose n = 2n . A choice that works essentially ^ 1 just as well is n^ = n^. The consequent reduction in memory

requirements and computing time is a made possible by the

periodicity of the fields in 0, which greatly improves the

accuracy of the integration over that variable C73.

We compute and store the components of the radius vector, of

the tangential vectors, the surface elements, and the

coefficients required for self— and nei ghbori ng—patch

corrections. These computations involve the evaluation of

tr i gonometr i c functions and logarithms, the performance of

numerical i ntegrati ons, and other time-consuming operations, but

they are done only once in each run.

3.4 Subroutine INFL

i n The incident field B has to be computed at the center of

13
each patch at each time step; this is done in subroutine INFL.

For the plane-Mave pulse, eq. (5) allows us to calculate the

magnetic -field -for arbitrary t. I-f the pulse is spatially

limited, the incident -fields have to be computed -from their

initial values by solving Maxwell's equations. This can be done

by numerical integration over the information-col lecting sphere

C 14] ,

3.5 Subroutine FARFL

The far field is calculated in subroutine FARFL by doing the

integration in eq. (13) for any number of scattering directions

given by the corresponding values of Q and ^ for the direction of

the scattered pulse. The intensities of the fields are plotted

and they are also saved on disk for comparison with the results

of other cal cul ati ons.

The time scale is referred back to the center of the sphere,

as is appropriate to compare the fields with the output from

PERF. The time the scattered pulse starts is a function of the

scattering angle,

t.^ = (a/c)Cl - 2sin<e/2)l, (25)

and the duration of the pulse is comparable to the length of the

pulse and to the size of the sphere.

14
4. Sel-f -Patch Corrections

The surface integral in eq. (9) is singular because the

integrand diverges when R = O, which happens on the self-patch.

The magnitude of the first term in the integrand is proportional

to R ^ and the integral of that term is well defined. The

-2

magnitude of the second term is proportional to R , and the

integral is defined in a limiting sense analogous to the

principal value of a one-dimensional integral. In an expansion

in the size of the patch C53, the leading term in the integral is

not absolutely integrable, but it vanishes by symmetry when the

field is computed at the center of the patch. The next term is

linear in the size of the patch.

The self-patch contribution to the integral is usually

neglected, although there is no more reason to neglect the

contribution of this patch than that of any other patch.

Actually, the self—patch contribution is linear in the size of

the patch, while those of most other terms are quadratic. Also

the seif—patch contribution is of fundamental importance in an

integral equation of the first kind, where the unknown field does

not appear outside the integral. We thus include in the program

the option to compute the self—patch correction.

The expansion of the integral of the first term is of order

2 in the size of the patch, but we cannot neglect this term when

is comparable to L73. In this case, because cAt is less

than R in the denominators, the integral of the first term is

larger than the integral of the second term. This happens, for

instance, when the pulse first reaches a patch. To make sure

15
that we do not neglect an important contribution, we keep this

term in all self-patch corrections. The importance of this term

complicates the implementation of a similar solution for an

integral equation of the first kind.

We include terms with spatial derivatives of the surface

current density, which are often ignored or assumed to be

negligible by requiring that the fields vary slowly over the

scatterer H153. We also have to take into account the spatial

variation of the tangential unit vectors.

The resulting expression for the self patch correction is

-+ AJ

a , 4n a^C^ - sin^ec^j (-J qS + J^i)

Ms - at ^ , 0^0, n

,

(26)

A 0

-+J ±

s

/aj

1 so*.-

. “-si 4 ax

" ay -•

}

a(A0)^ a so

i +

^

.A, 16c 1[at ax at

] ay atJ

j 0 = 0, n.

where the coefficients are given in C51 or C71, 0 and ^ are thi

16
tangential unit vectors along the coordinate curves, and A9 is

the angular width of a strip. To compute the corrections, we

need to know the fields and their time derivatives at the selfpatch, where the time retardation vanishes. We determine their

values from the uncorrected results. If the integral of the term

proportional to SJ^/dt is neglected, the correction is

proportional to and it can be incorporated in the coefficient

of the term on the left side of the equation.

Derivatives with respect to <6 are calculated by

approximating them by a quotient of differences of the respective

quantities for the patches on both sides. Derivatives with

respect to 0 are calculated in a similar manner. If the values

of the functions at the same ^ on the strips above and below are

not available, they have to be obtained by interpolation. There

IS an additional problem when one of the neighboring strips is a

polar cap because the 0- and fli -components of a tangential vector

are not well defined there? we define the 0-direction tangent to

the meridian that passes through the center of the self-patch,

and the correspondi ng tfk-direction from = n x 0. At the poles,

the derivatives with respect to x and y are obtained from the

values on both sides on the neighboring strip.

Since there is only one self-patch contribution to the

integral, the effect is hardly noticeable.

5. Neighboring-Patch Cgrrectigns

The integral over patches other than the self—patch is first

computed by multiplying the value of the integrand at the center

17
o-f the patch by the sur-face o-f the patch. This is a reasonable

approach it the integrand varies slowly over the patch. This is

2 not so for factors like 1/R or 1/R where R is small, that is,

for patches near the self-patch. Neighboring patches are those

that are closer than a small multiple (set in subroutine RUNT) of

the patch size to the center of the self-patch. We expand the

contributi ons from the neighboring patches in the parameters that

are of the order of magnitude of the size of the patch, which

here includes the displacements between the centers of the

patches C63.

As we explained in the previous section, the contribution of

the time-derivative term is not negligible when the current

density at the previous time is zero or small, and we also have

to expand this term L71.

We have to compute n x I, where

-4I R X, aj (x',T)

1 s

’ J (x" s

(27)

and is a neighboring patch. Since n is the unit normal at the

center of the self-patch S^, it is not per pendi cul ar to the

tangent! al vectors Q' and . If we expand n in terms of the

vectors related to the patch the expansion is great 1

y

si mpl i f i ed but the resulting contribution to J i s not

tangential. Also the change in n between neighboring patches is

not small when these patches are on a strip near a pole of the

sphere. Consequently, we expand I itself instead of n x I . We

18
abtai

n

I {[

(-1 + CQte^jA0^)C(l ,0,3) -cot9^ C (2,0,3)

a O

- 3cot0^jA0^C ( 1 , 2, 5) :ot0^jC(2,2,5)

f(-l + 2cot0lA0^ )C(0, 1 ,3) - -cot0;C(l, 1,3) + cos0:lA<i^C(l,O.

L O O • a O O O

•

- 3cot0^A9^jC<O,3,5) + |^ot0j^C(l,3,5) + j^|^(2,0,3)

- A0^C(l,O,3)]n^ + [|^sc0^;.C(O,2,3) - Afli^C (O, 1 , 3) j

io ' ) x

j^A0^,.C(l,O,3) - ^C(2,0,3)j9' + j^A0^.C (O, 1 , 3)

{[A^ - C ( 1 , O, 3) O -CSC0CC (1

a O

sin0^jA^^jC(O, 1 3) -C(0,2,3)

a

X

19

C-i4
1 + cQte^Ae^)C(i,o, :) - -cote;.c(2,o,2) d

- 2cote;;Ae^C(l,2,4) 0 0 ^cote^c<2,2,4)]e'

+ [(-1 + 2cote^Ae^j)C(0, 1,2) - |cote^c<i, 1,2) + cose^A*^jC( 1,0,2)

- 2cote-:Aa^C<0,3,4) + ^CDtO^C<l,3,4) + [1^(2, 0,2)

- Ae^C(l,0,2)]n^ + [|^sc9^C<0,2,2) Ad^C(0, 1,2) X

1c

- dJ

s at"

^|^A0^C<l,O,2) - ^C(2,0,2)je" + [Ae^C<0,l,2)

3 3

ic(l,l,2)]S'} X i

5 t? + {[^*^C(1,0,2) - icsce^C(l,l,2)]S'

CSC0;

-f aj

[sine^A(*^C (0,1,2) - ^C<0,2,2)]J"y x - ~ <28)

where 6' and are the tangential unit vectors at x", the center

20
of S^, jo' is the radial unit vector in the xy-plane, and the

current density and its derivatives are evaluated at at the

correspond! nq retards time The derivatives with respect to 0

or ^ take into account the dependence of the retarded time on 0

and The coefficients are given in C63 and C71.

Some of the coefficients become undefined when one of the

sides of the neighboring patch has the same coordinate as the

center of the sef-patch C73. We then have to evaluate these

coefficients from alternative expressions. Also the azimuthal

angles of the first and the last patches on a strip differ by an

amount close to 2nj we have to subtract 2n from such a difference

to get the proper value of for the calculations.

If the neighboring patch is the spherical cap at 0 = O,

n, we use Cartesian coordinates and obtain

I it (0,0,3) - iC' <1,0,3) - jCM0,l,3)

+ <0,0,3) - C' <2,0,3) - CM0,2,3)jy x

+ []o^jCMl,0,3) - iCM2,0,3) - jC'<l,l,3)j x 3J^/dx'

+ p^CM0,l,3) - iCMl,l,3) - jCM0,2,3)j x dJ^/dy'

21
+ ^

-(fl^CM0,0,2) - iCMl,0,2) - jCM0,l,2)

^ “ CM2, 0,2) -CM0,2,2)j^ x 3J^/3t'

+ -cLOf2.CMl,0,2) - iCM2,0,2) - jCMl, 1,2)1 x 3J /3t'

J 3x' s

+ - f2^CMO,l,2) - iCM 1,1,2) - j’CM0,2,2)l x /3t', (29) clO ’ -* J 3 y s

where

-f *’' = asine , jsin*^), (30)

and the coefficients are defined in C7J. The computation of

these coefficient involves a numerical integration C7D, for which

we use the subprogram QIDB in CMLIB C161. The integrand becomes

undefined for certain values of the variable of integration, and

an alternative expression has to be used where this happens.

The corrected nei ghboring-patch contr i but i ons replace those

obtained in the simple way. We have compared the two computed

values of the nei ghbor i ng-patch contributions and we have found

that sometimes they differ little while at other times the

difference can be large. Since a patch can have about ten

neighboring patches, the overall effect of this correction is

more significant than that of the self-patch.

22
6 „ Vector PrQacamming Considerations

The largest savings in execution time can be achieved by

vectorizing the computations in the main loop.

The values of the components of that may be needed for

the integration are kept in memory in two-dimensional arrays,

where one index represents the patch number and the other the

time. After these fields are no longer needed, they are saved on

a disk file.

An important observation that facilitates the vectori zat i on

of the i nterpol at i on in the time variable is that the location of

the retarded fields in the two-dimensional array for a fixed

field point is the same for the different components of the

fields and is shifted by one column as the time is incremented

for the interpol ation. Thus the array IND is computed once and

used in a QSVGATHR function repeatedly. This function is used to

collect the values of the components of at the different

retarded time in a single array. We originally computed the

values stored in the array IND and the related array FT once for

each field point and saved them, to be used at each time. This

procedure requires square arrays for IND and FT, increasing the

memory requi rements. Although the CYBER 205 has virtual memory,

paging generates much activity and slows significantly the

execution of the program. The arrays that contain the values of

the surface current density components are accessed in an

irregular manner and they have to be in memory at the time of

execution. The arrays IND and FT may be paged in gradually as

the field point changes, but we have chosen to recompute their
values -for each time point with only a small increase in the

running time. We save square bit arrays BT3 and BT4 to control

where the interpolation has to be shi-fted from the center of the

i nterval

.

The vector i zati on of this main loop produces a reduction of

the execution time by a factor of 7 if the sphere is divided into

70 patches, and 20 for 550 patches. These factors can probably

be reduced somewhat by improving the scalar version of the

program, where we use the three Cartesian components of the

fields instead of the two tangential components we use in the

vector program. The vectori zati on has to be done explicitly to

be able to take advantage of the Q8VGATHR function.

The subroutine RUNT is called only once at the beginning of

the program, and we made no great effort to vectorize this part

of the code. The subroutine FARFL can be called repeatedly to

find the scattered fields in different directions, and its

structure is not very different from the main loop since we are

also doing surface integrals over retarded fields.

The computation of the self—patch correction takes only an

additional one percent for execution. The additional execution

time for the nei ghbori ng-patch corrections depends on the number

of neighboring patches, which in turn depends on the criteria

used to define them and the patch distribution on the sphere. It

varies between ten and 40 percent of the basic execution time for

the loop. It may be possible to decrease this time by

vectorizing the computation of the nei ghbor i ng-patch

contributions.

24
7. Anci.1 jLary PCQQCSfDS

We have developed two programs that we use to examine the

results of the main programs PERF and SCAT. They are PERFPLT,

and CURCHK. These programs require very little execution time,

but they run on the CYBER 205 because that is where the output

-files are stored.

We have written a program SCATP that runs on the front-end

computer to determine the required dimensions of the different

arrays in SCAT and the number of large pages needed to accomodate

these arrays.

We have also written vector versions of subroutines NUFT and

INUFT, which we use to compute the direct and inverse Fourier

transform of functions that are not given at equispaced points.

7.1 Program PERFPLT

Program PERFPLT uses the output files PERFSV of PERF and

SCFLxxx of SCAT, where xxx stands for the qualifier that

identifies the run, and allows us to plot the far—field

intensities on a single plot for each outgoing direction of

propagation. We thus can compare the plotted output of PERF with

up to 4 plots produced by SCAT. We can choose printer plots by

using DRAW3 ClOl or plots that can be produced on a Versatec

pr 1 nter /pi otter by using DRAW4 Cl ID adapted to FORTRAN 77.

Actually, the information for the plots is added to a file stored

on the front-end computer. Then, when several plots have

accumulated, this information is transfered to tape which is

25
taken to our local PE 3230 minicomputer -for plotting.

Me can select the initial and -Final times on the plots to

show details o-F particular regions of the pulses.

7.2 Prograrn CURCHK

Program CURCHK selects the values of the components of

saved in a file SCSVxxx for up to 10 patches for the whole time

range or a given interval. These values can be printed or they

can be plotted on the printer or on a plotter and they provide a useful check on the solution of the integral equation.

7.3 Program SCATP

We run this program on the CYBER 180/855, the front end, and

then use the editor to pass the information to the parameter

statements in SCAT and CURCHK. The class and priority of a job

on the CYBER 205 depends on the memory and time requi rements, and

it is better to keep them as small as possible.

Program SCATP does the part of the calculations that will be

made in SCAT to find the number of patches for the given n^, n^,

and n ^5 estimates the number of time steps that have to be kept

in memory, and the number and distribution of neighboring patches

if this correction will be made.

A discrepancy arose in these calculations due to the

differences between the CYBER 170/855 and the CYBER 205. The

former returned the result cos(n/6) = .4999... while the latter

26
returned cos(n/6) = .5j rounding led to a different number of

patches and an insufficient size for the arrays; this is why the

size is incremented by 2 in SCATP.

7.4 Subroutine VINUFT

This subroutine is used by program PERF to compute the

inverse Fourier transform of the scattered fields. It is related

to subroutines NUFT and INUFT [91.

The subroutines NUFT can be used to compute the Fourier

transform of a function f (t)

,

a 5 t 5 b when the values of f (t)

are known at a set of values t^ that are not equi spaced. The

function is approximated by straight line segments joining

adjacent points, the Fourier transform of this approx i mat i on is

determined anal yti cal 1 y , and the result is evaluated numerically

at a given set of values w.

.

1 Subroutine NUFT can be used twice to obtain the inverse

Fourier transform in the analogous approximation, but we can do

this more efficiently by using subroutine INUFT if we know that

f(t) IS real, that is, when

f (-to) = f *(w:> . (31

)

Subroutines

NUFT and INUFT,

versions do not

function used in

9NUFT and VINUFT are the vectorized versions of

and they are listed in the appendix. The vector

have the small angle approximation for the sine

the scalar versions.

27
9" EHQQCSfl! iQByfes and Out guts

The inputs required to run PERF and SCAT are contained in

the files PERFIN and SClNxx (filenames on the front end can have

7 characters only, and we add another digit to indicate the type

of correction selected). The procedure used to submit a run to

the computer causes a file SCINxxx to be saved, so that the input

file SCINxk may be modified.

The information that is required for a given run relates

both to the physical problem and to the numerical computation.

The constants a and ff determine the shape of the pulse, and

we also need the radius a of the sphere. Actually, the value of

a simply fixes the scale of the problem, and we could set a = I everywhere in the program. The outputs for a = 1, a, B, and t

are the same as those for a, a/a, fi/a, and at.

We have to specify the angles 0 and ^ of the direction of

propagation of the scattered fields, and we may give times t^ and

t„ to obtain plots of some part of the scattered fields. Several

sets of directions and times may be specified, and all the

cai cui at i ons use the same previously determined surface currents.

If tj^ or t^ is not given, we make an estimate of the extreme

values based on the time of arrival of the pulse and the time

covered by the calculation of J^. Further sets of fields may be

computed in separate runs without recomputing the surface current

density, which is saved in SCSVxxx.

The parameters needed for the computation of the fields in

SCAT are the time increment At and the patch distribution on the

sphere. The distance light travels in the time At has to be

28
smaller than the separation betMeen the centers of any two

patches; the program SCAT makes this At equal to O.S times the

smallest separation between the centers if none is given or if

the given At is larger than this quantity. We can also give the

number of time increments in the calculation, which determines

trie time interval for which the surface current density is

determined; if this number is not given, the program selects a

number that corresponds to propagation over a distance 6a. The

division of the sphere into patches is controlled by the number

of strips, n^, the approximate number of patches on the

equatorial strip, n^, and the number of patches on the strip next

to a polar cap, n^, as discussed in section 3. To do the

calculation is PERF we also have to give the first frequency

increment, the highest frequency, and the number of frequency

values where the monochromat i c fields are to be calculated.

The agreement between the curves computed via the Fourier

transform and by solving the MFIE is shown in figure 1. The

output of SCAT in this figure corresponds to the largest run we

made. We needed 2S large pages, which is essentially the size of

the computer memory, the main loop executed in 10754 s with no

self- or nei ghbori ng—patch corrections, and the total CPU—time

for the job was 10812 s. The sphere was divided into 58 strips,

and this was also the maximum specified for the number of patches

in a strip (in powers of 2, the actual number goes from 8 near

the poles to 64 at the equator), for a total of 2898 patches. We

-10 selected At = 0.3 x 10 s for the time interval and computed

the current density for 900 time steps. There is some evidence

29
8C0TTERE0 INTENS11V. HORIZONTAL POLARIZATION SCATTERED INTENSITY, VERTICAL POLARIZATION

SCATTERED INTENSITY, HORIZONTAL POLARIZATION SCATTERED INTENSITY, HORIZONTAL POLARIZATION

Intensity (in

i a -Function of

Figure 1.

fields a

and (c) at 90 degrees with both pol arizati ons,

backward direction for an incident pulse with a =

f ar

(b)

the

* m . The solid curves are computed via a Fourier transform

from liie scattering formulas and the dashed curves with the

markers were computed via the MFIE. There is good agreement.

arbitrary units; of the scattered

ct in (a) the forward direction,

and_

j

^ = '^ a\ and

30
of oscillations at the later times, especially in backward

scat teri ng

.

Many factors affect the accuracy of the computation of the

fields from the integral equation, and these factors interact in

complex ways. The two effects of the i nterpol at i on formula used

for the time variable are illustrated in figure 2 for runs in

which the sphere was divided into 19 strips with 8 to 32 patches

in a strip (370 patches total) or 5 to 20 patches in a strip (258

patches total), and with At = 2.75 x 10 s, which is a

relatively large interval. Near the first peak, the curves that

correspond to the fourth—degree polynomial i nterpol at i on are much

closer to the correct solution, while near the end of the graph

these curves show the onset of oscillations. These oscillations

tend to become larger as time progresses, limiting the usefulness

of this interpolation scheme. We also see in figure 2 that the

particular way in which the strips are divided into patches,

whether we use powers of 2 or not, has little effect on the final

r esul ts.

In figure 3 we show a details of the curves about the first

peak to illustrate the effect of the self— and nei ghbor i ng-patch

corrections. The sphere was divided into 1986 patches, and the

time increment was 0.5 x 10 s. We see that the effect of

including the contribution of the self-patch is slight, while

that of correcting the contributions of the neighboring patches

is more significant.

The periodicity of the surface fields as functions of the

azimuthal angle ^ has the effect of reducing the number of

patches that are required in a strip L71. The trapezoidal rule
INTENSITY

(xl0-M

SCPITTERED INTENSITY, HORIZONTRL POLRRIZRTION

Figurs 2. Illustration of the effects of the subdivision of the

strips and the interpolation scheme. The curves marked with O

and 2 correspond to patch numbers that are powers of 2 and

curves O and 1 are computed wit!i 1 east—squares interpolation.

The curves were cut off at ct = 2.8 m.

gives a much more accurate integral for periodic than for

aroitrary functions when the integration is over one period.

if the sphere is divided into n strips, it is better to

divide the equatorial strip into n^ patches instead of the 2n^

required to make the patches more or less square. The smallness

of the effect of this or even a greater reduction is shown in

figure 4. The increase in the size of the patch has an adverse

effect in the computation of the self- and nei ghbor i ng-pat ch

32
SCATTERED INTENSITY. HORIZONTAL POLARIZATION

LO

Figure 3. Effect of seif-patch and neighboring-patch

corrections. We show a detail of the curves near the first

maximum. Curves 2 and 3, computed with nei ghbor i ng-patch

corrections, agree better with the solid curve than 0 and 1.

The inclusion of self—patch corrections in the computation of 1 and 3 does not lead to a noticeable improvement over O and 2,

respect i vei y

.

corrections, which is also illustrated in figure 4.

We show the components of the current density on the top cap

in figure 5a; the x—component remains essentially equal to zero,

as expected by symmetry. The computed surface current density

remains exactly equal to zero until the incident wave reaches the

patch; after that time, the current density is determined by the

integral equation. On the other hand, the actual fields cannot

33
INTENSITY

(xl0-M

SCATTERED INTENSITY. HORIZONTAL POLARIZATION

Figure 4. 1 1 i ustrat i on o-f the small effect of a variation in the

number of patches per strip. The sphere is divided into 29

strips, and the number of patches in the equatorial strip is 15

for curve O, 29 for curve 1, and 45 for curve 2. Curve 3

includes corrections with the same patch distribution used

for curve O, showing a deter i orat i on of the results.

reach a patch until they go around the outside of the perfectly

conducting sphere. This means that the solution of the integral

equation at the top cap of the sphere has to remain zero between

ct = 2a and ct = (1 + n)a, an effect that is illustrated in

figures 5b to 5f . The di screpanci es shown for the different

calculations are an indication of the accuracy that was attained

in the solution of the integral equation.

The fields computed for the program-sel ected value of At are

34
Figure 5. Components of the current density (in arbitrary units)

at the top pole of the sphere; the x-component remains zero.

We show the complete graph (a) in the run with 2898 patches

with 1 east-squar es i nterpol at i on , and a detail of the region

between ct = 2 m and ct = 1 + n/2 m for (b) the same run, (c)

one with 1586 patches, (d) same with corrections, (e) one with

370 patches, and (f) 370 patches with polynomial i nter pol at i on

.

35
8CRTTERED INTENSITY, HORIEONTPIL POLRRIZPITION

=»

Figure 6 . Fields computed with prog^0m—sel ected (minimum) values

of At (approx i mate! y 3.5 x 10 s) and number o-f time

increments with 370 patches. The computation without

corrections shown in curve O is compared to the same with

polynomial interpolation corrections in curve 1 and with a time increment o-f 1 . x 10 s in curve 2 .

compared with those computed with a smaller time increment in

-figure 6 . The curves show that it is important to use the more

accurate polynomial interpolation Tor a small number o-f patches.

The scattered -fields -for a sharper incident pulse, which has

constants a = 6 m^, ^ = 9m^, are shown in -figure 7. The fit

at the peaks is not very good and there are oscillations at later

times. Smaller patches are required for a sharper pulse.

36
SCATTERED INTENSITY, HORIZONTAL POLARIZATION SCATTERED INTENSITY. HORIZONTAL POLARIZATION

SCATTERED INTENSITY. VERTICAL POLARIZATION SCATTERED INTENSITY. HORIZONTAL POLARIZATION

Figure 7. Intensity of the scattered

f orward di rect i on

(b) and <c ) at

polarizations, and (d) the bac

i nci dent pulse with a = 6

X m and ^ = correspond to runs wi th 1586 patches.

oscillations after the initial peak.

far fields in (a) the

90 degrees with both

kwar^ direction for an

9 m . The graphs, which

show an increase in the

37
9. Concl usi ans

The agreement between the results of programs PERF and SCAT

shows that both programs are essentially correct. Thus, program

SCAT can be used as the basis for more general programs that

compute scattering of pulses by arbitrary perfectly conducting

surfaces. The parts of the program that would have to be

modified are mainly those related to the division of the surface

into patches and the computation of spatial derivatives. The

scattering of spatially localized pulses can also be included in

the program by making the appropriate changes in the computation

of the incident fields.

The importance of the time-derivative term in the

computation of the contribution of the self -patch causes

difficulties in the use of a similar procedure to solve integral

equations of the first kind, such as the EFIE or the integral

equation that describes the scattering by a dielectric.

The correction for the nei ghbori ng—patch contributions has a greater effect than inclusion of the self-patch, probably because

there are a number of neighboring patches for each patch. Gn the

other hand, the uncorrected integrals also gave satisfactory

solutions with a sufficiently large number of patches.

The vector i zati on of the program reduces the running time by

a factor that depends on the number of patches on the sphere.

This factor is approx

i

matel y equal to 20 for 550 patches. We

also found that straight polynomial interpolation gave better

accuracy than a least-squares polynomial interpolation at the

cost of increased instability.

38
References

C13 Marx, E. Integral Equations for Transient Ei ectromagneti

c

Fields, National Bureau of Standards Technical Note 1157,

February 1982.

C2D Marx, E. Single Integral Equation for Wave Scattering,

Journal of Mathematical Physics 23, 1057 (1982).

C31 Marx, E. Integral Equation for Scattering by a Dielectric,

IEEE Transactions on Antennas and Propagation AP—32, 166

(1984)

.

i;41 Marx, E. Transient Fields in Dispersive Media, Journal of

Mathematical Physics 24, 2602 (1983).

C53 Marx, E. Self—Patch Integrals in Transient El ectromagnet i

c

Scattering, IEEE Transactions on Antennas and Propagation

0Ez32, 763 (1985).

L6J Marx, E. Nei ghbori ng—Patch Integrals in Transient

El ectromagneti c Scattering, IEEE Transactions on

Antennas and Propagation AP-32, 768 (1985).

L71 Marx, E. El ectromagnet i c Pulse Scattered by a Sphere, IEEE

Transactions on Antennas and Propagation, submitted for

publ i cation.

183 Sookne, D. J. Bessel Functions of Real Argument and Integer

Order, Journal of Research of the National Bureau of

Standards, 77B, 125 (1973).

C93 Brandstein, A. G. and Marx, E. Numerical Fourier

Transform, Harry Diamond Laboratories Technical Report 1748

September 1976.

1103 Marx, E. Printer Version of Plots Made by an Incremental

Plotter, Harry Diamond Laboratories Technical Manual 75—33,

December 1975.

t 1 1 3 Noon, T. V. User’s Manual for Modular Anal ysi s-Package

Libraries ANAPAC and TRANL, Harry Diamond Laboratories

Technical Report 1782, November 1976.

C123 Hildebrand, F. B. Introduction to Numerical Analysis, p.

290 (McGraw—Hi 1 1 , New York, 1956).

11133 Mieras, H. Time-Domain Scattering from Dielectric Solids,

Sperry research Center, Sudbury, MA, June 1980.

C143 Marx, E. Free—Space Propagation of Light Pulses, National

Bureau of Standards IR 84-2835, May 1984.

39
II15II Mi eras, H. and Bennett, C. L. Space—Time Integral Approach

to Dielectric Targets, IEEE Transactions on Antennas and

Propagation AP-30, 2 (1982).

tl6D Boisvert, R. F. , Howe, S. E. , and Kahaner, D- K. Guide to

Available Mathematical So-ftware, National Bureau of

Standards IR 84—2824, January 1984.

40
AgEendiM

In this appendix we provide listings -for the programs SCAT,

SCATP, PERF, PERFPLT, and CURCHK together with most of the

subprograms that they require. We also provide listings for the

Fourier transform subroutines VNUFT and VINUFT.

The plotting subroutines DRAW3 and DRAW4 are documented

elsewhere, as are the integration subroutine QIDB, the interval

generating subroutine OMEGA, and the subroutines to compute

Bessel functions BESRJ and BESRY.

41
PROGRAM SCAT

rC VECTORIZED VERSION OF PROGRAM THAT COMPUTES SCATTERED FIELDS.

C SELF-PATCH AND NEIGHBORING-PATCH CONTRIBUTIONS CAN BE INCLUDED.

C MO CHARGE DENSITY IS COMPUTED

C

C XRl, XR2, XR3; ARRAYS WITH THE COMPONENTS OF THE POSITION VECTOR

C ON THE SPHERE

C SRF: (SOLID ANGLE SUBTENDED BY EACH PATCH)/(TWO PI)

C TF: TIME OF FIRST ARRIVAL OF THE WAVE AT EACH PATCH

C CO: COEFFICIENT FOR CONTRIBUTION TO SELF-PATCH

C DO. SAME AS CO.. FOR THE TIME-DERIVATIVE TERM

C TH2> TH3: COMPONENTS OF THE UNIT VECTOR ALONG MERIDIANS

C PHI, PH2: COMPONENTS OF UNIT VECTOR ALONG PARALLELS

C FUl,. FU2, FD.1, FD2, FCl, FC2: INFORMATION FOR INTERPOLATION

C THE! A: POLAR ANGLE OF PATCH

C PHIJ; AZIMUTHAL ANGLE OF PATCH

C DPHII; INCREMENT OF AZIMUTHAL ANGLE

C RB.1, RB2, R.B3, RB4; ARRAYS THAT CONTAIN BIT ARRAYS

C COJ, COJT, COJU, COJTU, COJV, COJTV; ARRAYS WITH COEFFICIENTS

C NEEDED TO COMPUTE NEIGHBORING-PATCH CORRECTIONS

C SINTl, SINT2, SI NTS: CARTESIAN COMPONENTS OF INTEGRALS

C DRl .. DR2, DR3: COMPONENTS OF R = X - X

' C DR, C.0RM2, DRM3: WORK ARRAYS

C CJTHD, CJPHD: COMPONENTS OF THE RETARDED SURFACE CURRENT DENSITY

C IN THE THETA- AND PHI -DIRECT IONS

C DJTHD, DJPHD: SAME FOR TIME DERIVATIVE OF SURFACE CURRENT DENSITY

C UTH.. UPH, AO, Al, A2, A3, A4, YM2, YMl, YO, Yl, Y2‘ WORK ARRAYS

C CJTHA, CJPHA; COMPONENTS OF SURFACE CURRENT DENSITY AT TIME T

C CJTH, CJPH: TWO-DIMENSIONAL ARRAYS WITH COMPONENTS OF ALREADY CALCULATED

C CURRENT DENSITY

C FT INFORMATION FOR TIME INTERPOLATION

C ML, NR, NUI, NU2, NDl, ND2; INFORMATION FOR SPATIAL INTERPOLATION

C NET, NPQINT; INFORMATION FOR NEIGHBOR ING-PATCH CORRECTION

C IND: INDEX ARRAY FOR VECTOR PROCESSING

C IJ: IJ(I) I, FOR VECTOR PROCESSING

C BT.l BIT ARRAY SELECTS ALL PATCHES EXCEPT THE SELF PATCH

C BT2: BIT ARRAY TO INDICATE WHETHER THE PULSE HAS REACHED A PATCH

C BT3, BT4: BIT ARRAYS TO INDICATE NECESSITY OF SHIFTING INTERPOLATION

C

PARAMETER (LL==LLG, LC=LCG, LL12=LL12G, LL20^-LL20G.

1

LB ^:LL**2/64+l, LB 1 -LL/64+1 ) BIT BTl (LL), BT2(LL), BT3(LL, LL ) , BT4(LL, LL

>

CHARACTER TITLE^l•80, SUBT^^SO, XLAB*30, YLAB*30, QAL.-»4

COMMON /RCOM/ XR 1 ( LL

)

, XR2 ( LL ) , XR3 ( LL ) , SRF ( LL ) , TF ( LL )

,

1 CO(LL), DO(LL), THl (LL), TH2(LL), TH3(LL),

2 PHI (LL), PH2(LL), FUl (LL), FU2(LL), FDl (LL), FD2(LL), FCl (2), FC2(2).

3 THETA (LL), PHII (LL), DPHII (LL), RBI (LBl ), RB2(LB1 ), RB3(LD). RB4(LD).

4 COJ (6, LL12) , C0JT(6, LL12), COJU (2, LL12) , C0JTU(2, LL12)

-

5 C0JV(2, LL12), COJTV (2, LL 1 2 ) , SI NTl ( LL ) , SINT2(LL), SINT3(LL).

6 DR 1 ( LL ) , DR2 < LL ) , DR3

(

LL ) . DR ( LL ) , CDRM2 ( LL ) , DRM3 ( LL )

,

7 CsJTHD(LL), CJPHD (LL), DJTHD(LL), DJPHD (L.L),

8 UTH ( L L ) , UPH ( LL ) , AO ( LL > , A 1 < LL ) , A2 ( LL ) , A3 ( LL ) , A4 ( L.L ) , 9 YM2 ( LL ) , YMl ( LL

)

, YO ( LL

)

, Yl ( LL ) , Y2 ( LL )

,

42
A CJTHA<LL), C-JPHA(LL), CJTH<LL, LC ) , CJPH(LL, LC)> FT(LL ),

B NL(LL), NR(LL), NUl (LL>, NU2<LL)> NDl (LL), ND2(LL),

C NEI <LL, LL20), NPOINKLL, LL20), IND(LL). IJ<LL)

COMMON /ICOM/ IX, IS, ISl, Nl, N2, NI, ICORR, IPLOT, INTER, IDEB,

1 IPATCH, N3, NU

COMMON /LABEL/ TITLE, SUBT, XLAB, YLAB

COMMON /PARAM/ ALPHA, BETA, DT. RAD, RADMl , CDTMl , THSC, SR i , SR2, SR3

COMMON /CONST/ AMUO, TWOMUO, TWOEPO, C, CMl, PI, TOPI, PITO, OS, OTW, OTF

EQUIVALENCE (RB1,BT1), <RB2, BT2), (RB3, BT3), <RB4, BT4)

DATA IZO/0/, BT1/LL*E ' 1

V

WRITE (XLAB, 7)

WRITE! YLAB, S)

LL2=2+^LL

OS--=l. /7.

OTW--^l. / 1 2.

DTF:=1, /24.

C==3. E8

CM 1=1. /C

PIT0=2. -itATANU. ) PI-^2. -«PITO

TOP 1-^4, ftp I TO

CDNV==PI/180.

AMU0-'--PI-ft4. E--7

TW0MU0=2. /AMUO

TW0EP0=2. *CM 1^*2/ AMUO

PRINT MINI)"

READ ICORR, ISKIP, IPLOT, INTER, IDEB, IPATCH. QAL

PRINT -ft, '• ICORR==^ ICORR, S ISKIP= S ISKIP, ^ IPLOT^^: IPLOT,

1 L. INTERS L. INTER, ^ IDEB= ' , IDEB, IPATCH= '

, I PATCH,. '

, QAL= ^ QAL

READ IN VALUES OF CONSTANTS FOR THE PULSE

READ 4, DT, RAD, ALPHA, BETA

READ 5, NI

READ 5, N1,N2, N3

IF<N1, GT. 100) STOP 2

IF(W3. EQ. 0) N3=^5

R ADM 1 = 1. /RAD

IF( ISKIP. EQ. 1 ) 1C0RR=0

COMPUTE COMPONENTS OF VECTORS ON

SELF-PATCH AND NEIGHBORING-PATCH

CALL RUNT

THEN

I XRl ( I

)

TF(I) CO(I)

TH2 ( I ) TH3 ( I

)

PHKI) PH2(I)

SPHERE AND INFORMATION NEEDED FOR

CORRECTIONS

IF< IDEB. EQ. 1 ) WRITE (7, ft)

1 'SRF(I)

2 ^THl(I)

WRITE <7, ft) '

XR2( I

)

Don

)

FUl ( I

XR3( I

)

1

' FU2 ( I ) FD 1 ( I ) FD2 ( I )

' WRITE(7, ft)' NL(I) NR(T) NUl < I

)

NU2(I) NDl ( I ) ND2(I)'

DO 90 1 = 1, IS

WRITE (7, 2) I, XRl ( I ), XR2( I ), XR3( I ), SRF < I ) , TF( I )

,

1 C0( I >, D0( I ) , THl ( I ), TH2< I ),

2 T H3 ( I ), PH 1 ( I ) , PH2 ( I ) , FU 1 < I ) , FU2 < I ) , FD 1 ( I ) .• FD2 ( I ) WRITE(7, 1 ) NL( I ), NR( I ), NUl ( I ), NU2( I ), NDl ( I ), N.D2( I )

43
90

9f-;

C

100

C

C

CONTINUE

IF( ICGRR. GE, 2) THEN

DC) 95 IS

WRITE (7, 1) I, <NEI < I, J), NEI (I, ,l ) ) CONTINUE

END IF

END IF

NO CALCULATIONS IF NI <' 0

IF(NI.LT. 0) STOP

DO 100 I==-l, IS

I J( I )=I

CONTINUE

DTr'11--=:l. /DT

CDTM^=CMH«-DTM1

COMPUTE THE NUMBER OF TIME STEPS TO A DIAMETER TO DIMENSION ARRAY

NU-2, *RAD^<-CMHtDTMl+6.

COMPUTE THE NUMBER OF TIME STEPS IF NOT GIVEN

IF<NI, EQ. 0) NI=^6. *RAD/(C*DT)

PRINT 'NI^MvII, ^ N1=SN1>', N2='.N2>'> N3=UN3,

1 U DT==',DT, S RAD=-'. RAD. ^ ALPHA- ALPHA, '

, DETA=',BETA

PRINT NU. ' ROWS NEEDED OUT OF '.LC

IF(NU. GT. LC) STOP 3

NUM2=^NU-2

IF< ICORR. EQ. 1 > THEN

PRINT Jt. 'SELF-PATCH CORRECTIONS ONLY'

ELSE IF( ICORR. EQ. 2) THEN

PRINT 'NEIGHBORING-PATCH CORRECTIONS ONLY'

ELSE IF( ICORR. EQ. 3) THEN

PRINT 'SELF- AND NEIGHBOR ING-PATCH CORRECTIONS'

ELSE

PRINT •«.. 'MO PATCH CORRECTIONS'

END IF

ENCODE INFORMATION FOR TITLES OF GRAPHS

IF( ICORR. EQ. 1 ) THEN

IF( INTER. EQ 0.^ THEN

WRITEv'SUBT. 11 ) ALPHA. BETA. DT. NI. Nl. N2. N3

ELSE

WRITE(SUj3T, 12) ALPHA. BETA. DT. NI . Nl , N2. N3

END IF

ELSE IF ( ICORR. EG. 2) THEN

IF ( INTER. EQ. 0) THEN

WRI TEISUBT. 13) ALPHA. BETA. DT. NI. Nl. N2. N3

ELSE

WRITETSUBT. 14) ALPHA. BETA. DT. NI. Nl. N2. N3

END IF

ELSE IFdCORR. EQ. 3) THEN

IF ( INTER. EQ. 0) THEN

WRrTE(SUBT. 15) ALPHA. BETA. DT. NI. Nl. N2. N3

ELSE

WRITE(SUBT. 16) ALPHA. BETA. DT. NI, Nl. N2. N3

END IF

ELSE

IF ( INTER. EQ. 0) THEN

WRITE (SUBT. 9) ALPHA. BETA. DT. NI, Nl. N2. N3

44
ELSE

WRITE(SUBT, 10) ALPHA, BETA, DT, NI, Nl, M2, N3

END IF

END IF

CALL Q5RQUEST( 'LFN= 'SCFL ' //QAL, 'MXR=^ S 2400. 'LEN=^50,

1 'RT=^ 'W' ) OPEN (9, FILE=='SCFL V/QAL, FORM== ^UNFORMATTED ' )

T--=--DT NF2^S-ttIS

LENGTH=^N I * I S-«-2/ 512+10

IF< ISKIP. EG. 1 ) THEN

C COMPUTE FAR FIELDS FROM SAVED FILE SCSVLQAL

CALL Q5ATTACH( 'LFN=^ •'SCSVV/QAD

OPEN (2, FILE='SCSVV/QAL, STATUS="OLD^ ACCESS:= 'DIRECT

1 FORM== 'UNFORMATTED', RECL:=NF2)

T--NI-ttDT

QO TCI 400

END IF

CALL Q5RQUEST( 'LFN=--=', 'SCSV'//QAL, ' MXR=-' ' , NF2, 'LEN= ', LENGTH,

1 'SFO=', 'D')

OPEN (2, FILE='SCSV'//QAL, FORM-= 'UNFORMATTED',

1 RECL=-NF2, ACCESS=--= 'DIRECT ' ) TL=SECOND(

)

NUM=--a

NREC^=M

I XH----0

/-• C INITIALIZE ARRAYS

r CJTHA( 1; IS)=0.

CJPHA< 1; IS)^^--0.

CJTHIX 1 ; IS)-=^0,

CJPHD( 1; IS)=-0.

DJTHD( 1; IS)===0.

DJPHDd; IS):=0,

UTH(1; IS)^0.

UPH< .1; IS)=^0.

DO 120 J=^1,LC

CJTH( 1, vJi IS)-0.

CJPPK 1. Ji IS)=0.

120 CONTINUE

r

C MAIN LOOP

C

DO 300 IZ--2,NI

NUM=NUM+1

PRINT 17, ' OF ', NI

C SHIFT THE NUMBERS IN THE ARRAYS

DO 170 1=2, NU

CJTH< ] , I~li IS)=CJTH( 1, I; IS)

C JPH < 1

,

I

-- 1 , IS) =C JPH < 1 , I ; IS)

170 CONTINUE

C COMPUTE DENSITIES AT EACH PATCH ON THE SPHERE

DO 200 IX=1, IS

45
o

o

o

o

C RETRIEVE THE RADIUS VECTOR (NORMAL.)

X=XR1 ( IX

)

Y=-vXR2(;(X ) Z~XR3( IX

)

C FIND THE INCIDENT FIELDS

CALL INFL(5INT10, SINT20, SINT30,

1 X, Y, Z, T. «-210)

DRl ( li IS)r-^X-XRl < 1; IS)

DR2(li IS)=--Y-XR2( 1, IS)

DR3(1, IS):--=Z-XR3(1, IS)

CDRM2U; IS)'^DR1 ( .1.. IS ) h^*2+DR2 (1 ; IS ) **2+DR3 < 1 ; IS)**2

DR( i; IS)=VSQRT(CDRM2< 1, IS); DR( li IS)

)

DRM3(1; IS)=CDTM1*DR< 1, IS)

IND( 1, IS)-VINT(DRM3( 1; IS); IND( 1; IS) ) FT (1 ; IS) =VA I NT ( DRM3 (1 ; IS); IS

)

-DRM3 (

1

; IS)

C SHIFT THE FIVE POINTS FOR INTERPOLATION IF NECESSARY

WHERE<BT3( 1 , IX, IS)

)

INDU, IS) = IND( 1.. IS)+1

FT( 1; IS)=FT( 1; IS) + l. END WHERE

WHERE(BT4( 1, IX; IS) ) IND( 1; IS)=IND< 1, IS)+2

FT( 1 , IS)=FT< 1

IS) +2.

END WHERE

INDd, IS)-IJ(1; IS)-IND(1, IS)*LL

,tm ( I X ) =B '0-' C DETERMINE WHETHER THE PULSE HAS REACHED THE PATCH

BT2(1; IS)==T-CM1*DR (1 , I S ) . LE. TF (1 , I S ) WHERE <BT1 ( ,t .. IS) ) DRM3< 1.. IS)^1. /DR ( 1; IS)

CDRM2(1; IS)=DRM3(1; IS)*^f2

DRM3(1; IS)=CDRM2(1; IS ) *DRM3 ( 1 , IS)

CDRM2(1; IS)-CDRM2(1; TS);«-CDTM1

C DETERMINE RETARDED INTERPOLATED FIELDS

IF (INTER. EQ. O) THEN

CALL INTRPOI IND, FT, CJTH, CJTHD, DJTHD, AO, Al, A2,

t YM2, YMl, YO, Yl, Y2)

CALL INTRPOdND, FT, CJPH, CJPHD, DJPHD, AO, Al, A2,

1 YM2, YMl, YO, Yl, Y2)

ELSE

CALL INTRPl ( IND, FT, CJTH, CJTHD, DJTHD, AO, Al, A2, A3, A4,

1 YM2, YMl, YO, Yl. Y2)

CALL INTRPl ( IND, FT, CJPH, CJPHD, DJPHD. AO. Al, A2, A3, A4.

.1 YM2, YMl, YO, Yl , Y2)

END IF

5L-T FIELDS TO ZERO IF PULSE HAS NOT REACHED PATCH

( INTERPOLATION MAY MAKE FIELDS NONZERO BEFORE THEY ARE)

WHERE (BT2d, IS) ) CJTHDd; IS)=0.

CJPHD (1; IS)=0.

DJTHD d; IS)=0.

DJPHDd; IS)=0.

END WHERE

DETERMINE THE CARTESIAN COMPONENTS OF THE CONTRIBUTION OF THE

PATCH TO THE INTE(3RAL

46
o

r;

UTH( 1, IS)=CDRM2( 1, I S ) «-DJTHD <1 , IS)+DRM3(1, IS)*CJTHD( li IB)

UPHd. IS):=CDRM2( 1.. IS)#DJPHD( li IS)+DRM3( li IS)-«-CJPHD( 1; IS)

SINTl ( 1; IS)=-SRF(1; I S ) « ( DR2 ( 1 , IS)4«-UTH< 1; IS)^^TH3( li IS)

1 -DR3(1; IS)*<UTH( 1 , IB)*TH2( li IS)+UPH( li IS)*PH2< 1, IS) )

)

SINT2( li IS) = -SRF( li IS)-?^«:DR3( li IS)-«-(UTH( li IS)*TH1 < li IS)

1 +UPH(li IS)*PH1 ( li IS) )-DRl( li IS)*UTH( li IS)*TH3( 1> IS)

)

SINT3(li IS)=-SRF(li IS)*(DR1(1, IS)*(UTH(li IS)^^TH2(li IS)

1 -+-UPH(li IS)*PH2(li IS))-DR2(li IS)*(UTH(li IS)-»THl(li IS)

2 +UPH( 1, IS)*PH1< li IS) )

)

SNT1=Q8SBUM(SINT1 < li IS). BTK li IS) )+SINT10

SNT2=-Q8SSUM<BINT2< li IS). BTl < 1, IS) )+SINT20

SNT3=Q3SSUM(SINT3( li IS).BTld. IS))+SINT30

BTl (IX )=B d ' C

C FILL IN APPROXIMATE VALUES FOR CURRENTS AND DERIVATIVES

IF( IX. EQ. 1 ) THEN

CJTHD( 1 )=SNT2

CJPHD< 1 )=-SNTl

CJTH( 1, NU)=SNT2

CJPH( 1. MU)^-SNT1

ELSE

CJTHD( rX)=-PHl ( IX)«-SNT1-PH2( IX)*SNT2

C JPHD a X ) =^TH1 (IX) *SNT1 +TH2 (IX) ^<SNT2+TH3 (IX) *SMT3

CJTH( IX, NU)=CJTHD( IX

)

CJPH( IX. NU)=CJPHD( IX

)

END IF

IF( INTER. E(3. 0) THEN

CALL I NTRP2 ( 2. . C JTH. I X . NU-2. DJTHDX

)

CALL INTRP2(2. . CJPH. IX. NU-2, DJPHDX

)

ELSE

CALL INTRP3(2. , CJTH, I X, NU-2. DJTHDX

)

CALL INTRP3(2. , CJPH. IX, NU-2. DJPHDX)

END IF

DJTHD( IX JTHDX

DJPHD( IX )=DJPHDX

COMPUTE CORRECTIONS FOR NEIGHBORING PATCHES

C IF ICORR = 0 OR 1. NEKIX.l) = 1 AND LOOP IS SKIPPED

DC) .190 IXX==2. NEI ( IX. 1 ) IXP^^NEI ( IX, IXX ) IF(T-CMt*DR( IXP). LE. TF( IXP) ) GO TO 190

NPNTn^NPOINT( IX. IXX ) CGJU^CQJ( 1. NPNT)

C0J2=C0J(2. NPNT)

CDJ3-“=C0J(3, NPNT)

CGJ4::=CGJ(4, NPNT)

COJ5=dlOJ( 5.. NPNT)

C(3J6^--C0J(6, NPNT)

COJTi-COJTv 1, NPNT)

CGJT2:^CQJT(2. NPNT)

COJT3=--C0 JT ( 3. NPNT ) CCJ JT4 =^CO JT ( 4. NPNT )

47
COJTS^^CQJKS. NPNT)

CDJT6==C0JT<6, NPNT)

COJUi=^COJU( 1, NPNT)

C0JU2=--C0JU(2. NPNT)

CDJTUl“-“COJTU( 1; NPNT)

C0JTU2-=CGJTU(2, NPNT)

COJVl=COJV( 1, NPNT)

C0JV2=CDJV(2. NPNT)

CDJTVl=--COJTV( 1. NPNT)

CGJTV2=C0JTV(2; NPNT)

IF( IXP. EQ. .1. OR. IXP. EQ. IS) THEN

COMPUTE NEIGHBORING-PATCH CORRECTIONS FROM CAPS, IF NECESSARY

IFdXP.EQ. 1) THEN

ISP •=2

ELSE

ISPr=IS.l

END IF

DETERMINE THE INCREMENTS WITH RESPECT TO X CvJlD2=^THl ( ISP)<t-CJTHD< ISP)+PH1 ( ISP ) ^^C JPHD ( I SP ) CJ2D2=--TH2( ISP )*CJTHD( ISP)+PH2< ISP ) *C JPHD ( I SP ) CJ3D2-d'H3( ISP)-K-CJTHD( ISP)

N1=^--NL( IXP ) N2==NR( IXP)

CJlDL^dHl (N1 )*CJTMD<N1 )+PHl (N1 )-»-CJPHD(Nl ) CJ2DL=--=TH2(N1 )*CJTHD(N1 )^-PH2(Nl ) «-CJPHD(Nl ) C-J3DL^=^TH3(N1 )-k-CJTHD(N1 ) CJIDR-THI (N2)«CJTHD(N2)+PH1 (N2)*CJPHD(N2)

C J2DR==^TH2 ( N2 ) ^«-C JTHD ( N2 ) +PH2 ( N2 ) -«-C JPHD ( N2 ) C J3DR:--=TH3 ( N2 ) -k-C JTHD ( N2 ) IFdXP.EQ, 1) THEN

DXCJl=^CJlD2-(CJlDLi^FCl ( 1 )+CJlDR*FC2( 1 ) ) DXCJ2=--CJ2D2-(CJ2DL*FC1 < 1 )+CJ2DR»FC2( 1 ) ) DXCJ3=CJ3D2-(CJ3DL^<•FC1 ( 1 ) +C J3DR*FC2 ( 1 ) ) ELSE

DXCJ1=^CJ1D2-(CJ1DL*FC1 ( 2 ) +C J1 DR^tFC2 ( 2 ) ) DXCv.i2=^CJ2D2~<CJ2DL-«-FCl ( 2 ) +C J2DR#FC2 ( 2 ) ) DXCJ3=CJ3D2-(C J3DLitFCl ( 2 ) +C J3DR-«-FC2 ( 2 ) ) END IF

DETERMINE THE INCREMENTS WITH RESPECT TO X N1--=NU1 ( ISP)

N2-^=NU2 d SP ) CJ1PU=(TH1 (N1 )*CJTHD(N1 )+PHl ( N1 ) *C JPHD ( N1 ) )»FU1 ( ISP)

,1. +(TH1 (N2)*C JTHD<N2)+PH1 (N2)#CJPHD(N2) )<^FU2dSP ) CJ2PU-^(TH2(N1 )*CJTHD(N1 )-*-PH2<Nl )-k-CJPHD(N1 ) )*FU1 ( ISP )

I + d'H2(N2)*C JTHD < N2 ) -<-PH2 ( N2 ) «-C JPHD (N2) )*FU2( ISP)

CJ3PU=TH3(N1 )#CJTHD(N1 )<^FU1 ( ISP )

1 +TH3<N2)^<-CJTHD(N2)*FU2dSP)

N1=NDI

(

ISP ) N2==ND2 (ISP)

CJ1PD==(TH1 <N1 )^^CJTHD(N1 )+PHl (N1 )*CJPHD(N1 ) )*FD1 ( ISP )

1 +<TH1 (N2)*CJTHD(N2)+PH1 ( N2 ) *C JPHD ( N2 ) )-«-FD2( ISP ) CJ2PD=(TH2(N1 > «C JTHD < N1 )+PH2(Nl )*C JPHD (N1 ) )#FD1 (ISP)

1 •+ ( TH2 < N2 ) *C JTHD ( N2 ) +PH2 < N2 ) *C JPHD ( N2 ) ) *FD2 (ISP)

CJ3PD-TH3(N1 )#CJTHD(N1 )*FD1 ( ISP

)

40
1 +TH3 ( N2 ) JTHD ( N2 ) ^^FD2 (ISP)

DYCJl=-=CJiPU-CJlPD

DYCJ2==CJ2PU--CJ2PD

DYCJ3:--^CJ3PU-CJ3PD

DETERMINE INCREMENTS OF TIME DERIVATIVES WITH RESPECT TO X

DJ 1 D2=^TH 1 ( I SP ) *DJTHD (ISP) +PH 1 ( I SP ) *DJPHD (ISP)

D J2D2-^TH2 (ISP) ^tDJTHD (ISP) +PH2 (ISP) -m-D JPHD (ISP)

D J3D2^^=TH3 (ISP) -»-D JTHD (IBP)

N1=--NL( IXP)

N2=^NR( IXP)

DJIDL^THI (N1 )-J^DJTHD(Nl )+PHl (N1 )*DJPHD(N1 ) DJ2DL:--TH2(N1 ) *DJTHD ( N1 ) +PH2 ( N1 )*DJPHD(N1 ) DJ3DL=-^TH3 ( N 1 ) *DJTHD ( N 1 ) DJ 1 DR^ TH 1 < W2 ) *D JTHD ( N2 ) +PH 1 ( N2 ) *DJPHD ( N2

)

DJ2DR=---TH2 ( N2 ) #DJTHD ( N2

)

+PH2 ( N2 ) *DJPHD ( N2

)

DJ3DR^-=TH3 ( N2 ) *DJTHD ( N2

)

IF( IXP. EQ. 1

)

THEN

DXDJl=--DJlD2-(DJlDL-i.^FCl ( 1 ) +DJ1DR*FC2 ( 1 ) ) DXDJ2=DJ2D2-(DJ2DL#FC1 ( 1 ) +DJ2DR*FC2 ( 1 )

)

DXDJ3=^DJ3D2“(DJ3DL#FC1 ( 1 ) +DJ3DR*FC2 ( 1 ) ) ELSE

DX DJ 1 DJ 1 D2 - ( DJ 1 DL-t^FC 1(2) +DJ 1 DR*FC2 ( 2 ) ) DXDJ2-DJ2D2-\DJ2DL*FC1 ( 2 ) +DJ2DR-»^FC2 ( 2 ) ) DXDJ3=DJ3D2- ( DJ3DL*FC 1(2) +DJ3DR*FC2 ( 2 )

)

END IF

C DETERMINE INCREMENTS OF TIME DERIVATIVES WITH RESPECT TO Y

N1===NU1 (ISP)

N2=--NU2 (ISP)

DJ1PU=(TH1 (N1 )^^DJTHD(N1 )+PHl (N1 )*DJPHD(N1 ) )-J^FUl ( ISP)

1 +(TH1 (N2)^^DJTHD(N2)+PH1 (N2)*DJPHD(N2) )-frFU2( ISP ) DJ2PU=(TH2(N1 )«DJTHD(N1 )+PH2(Nl ) J^DJPHD(N1 ) )*FU1 (ISP)

1 +(TH2(N2)^^DJTHD(N2)+PH2(N2)*DJPHD(N2) )*FU2( ISP

)

DJ3PU==TH3(N1 )*DJTHD(N1 )^^FU1 ( ISP)

1 +7 H3 ( N2 ) *DJTHD ( N2 ) *FU2 (ISP)

N 1 ===ND 1 ( I SP ) N2=^=N02 ( I SP ) DJ1PD=(TH1 (N1 )-i^DJTHD(Nl )+PHl (N1 )*DJPHD(N1 ) )*FD1 ( ISP )

;l +(TH1 (N2)*DJTHD(N2)+PH1 ( N2

)

*DJPHD ( N2 ) )-t^FD2( ISP ) DJ2PD=(TH2(N1 )^^DJTHD(N1 )+PH2(Nl ) »DJPHD(N1 ) )^«FD1 (ISP)

1 +(TH2(N2)*DJTHD(N2)+PH2(N2)*DJPHD(N2) )*FD2( ISP

)

DJ3PD==TH3(N1 ) *DJTHD ( N1 ) *FD1 ( ISP ) 1 +TH3 ( N2

)

^^D JTHD ( N2

)

*FD2 (ISP)

DYDJ 1 =--D J 1 PU~DJ 1 PD

DY DJ2=^ D J2P U- DJ2P D

C DETERMINE CORRECTION TO INTEGRAL FROM CAPS

DYDJ3:-“-DJ3PU-DJ3PD

SNT1P=^C0J3*CJPHD( I XP ) -C0JU2^^DXC J3-C0JV2-K-DYC J3

;l +C0JT3#DJPHD ( I XP ) ~CQJTU2-frDXDJ3-C0JTV2*DYD J3

SNT2P=-*-C0J3^fCJTHD( I XP ) +COJU1 *DXCJ3+C0JV 1 *DYC J3

1 -C0JT3*DJTHD( IXP)-(-CajTUl*DXDJ3 + C0JTVl*DYDJ3

SNT3P==--C0JH^CJPHD( IXP )+C0J2*CJTHD ( I XP ) -COJU U^DXC JPH

1 +C0JU2#DXCJTH-C0JVl•^tDYCJPH+C0JV2^^•DYCJTH

2 -COJTl-«-DJPHD( IXP )+COJT2»tDJTHD ( I XP ) -COJTU 1 itDXD JPH

3 +C0JTU2-frDXDJTH-C0JTVm-DYDJPH+C0JTV2^H)YDJTH

49
oo

SNT1=:SNT1-SINT1 < IXP)+SNT1P

SNT2=SNT2-SINT2 < I XP > +SNT2P

SNT3=SNT3“S I NT3 ( I XP ) +SNT3P

EL.SE

COMPUTE NEIGHBORING-PATCH CORRECTIONS FROM REGULAR PATCHES

DETERMINE INCREMENT OF FIELDS AND TIME DERIVATIVES WITH RESPECT TO THETA

CALL AVCUR( IXP, NUl. NU2. FUli FU2> CJTHUO, CJPHUO>

1 CJTHD, CJPHD)

CALL AVCURdXP, NDl, ND2, FDl, FD2, CJTHDO, CJPHDO,

1 CJTHDiCJPHD)

CALL AVCUR( IXP> NUl, NU2, FUl, FU2, DJTHUO, DJPHUO,

;l DJTHD. DJPHD)

CALL AVCURdXP, NDl, ND2. FDl, FD2, DJTHDO. DJPHDO.

.1 DJTHD,DJPHD)

DTHCJTH=-CJTHDO-CJTHUO

DTHCJPH=CJPHDO-CJPHUO

DTHDJTH==DJTHDO-DJTHUO

DTHDJPH^DJPHDO-DJPHUO

C DETERMINE INCREMENT OF FIELDS AND TIME DERIVATIVES WITH RESPECT TO PHI

IXl=-NRdXP)

IX2^WLdXP)

CJTHXl=CJTHDdXl ) CvJPHXi---CJPHDdXi ) CJTHX2=CJTHDdX2)

C JPHX2==C JPHD d X2

)

DPHC JTH=C JTHX 1 -C JTHX2

DPHC JPH==-C JPHX .1 -C JPHX2

DJTHXl=DJTHDdXl ) DJPHXl=DJPHDdXl ) DJTHX2=DJTHDdX2)

DJPHX2=DJPHDdX2>

DPHDJTH=DJTHX 1 -DJTHX2

DPHD JPH=--^DJPHX 1 -D JPHX2

CJTHDP=CJTHDdXP > CJPHDP=-CJPHDdXP)

DJTHDP---:D JTHD< IXP ) DJPHDP = DvJPHD< IXP)

C DETERMINE CONTRIBUTION TO INTEGRAL FROM NEIGHBORING PATCH

CJITH'^C0Jl*CJTHDP+C0J4<f-CJPHDP

DJITH:^C0JTH«-DJTHDP+C0JT4*DJPHDP

CJIPH^=C0J2*CJTHDP+C0U5*CJPHDP

DJIPH=COJT2^^DJTHDP+COJT5*DJPHDP

CJIN0=C0J3*CJTHDP+C0J6-»^CJPHDP+C0JUH^DTHCJTH

:l +C0JU2<^DTHCJPHh-C0JV1-«-DPHCJTH<-C0JV2*DPHCJPH

DJIN0=^^C0JT3#DJTHDP+C0JT6-t^DJPHDP+C0JTUl-«-DTHDJTH

.1. H-C0JTU2-S-DTHDJPH+C0JTV1-K-DPHDJTH+C0JTV2K-DPHDJPH

CITH==-CJITH+DJITH

CIPH^^CJIPH+DJIPH

CINO==CJINO+DJINO

THlP=-dHl dXP ) TH2P= TH2 (IXP)

TH3P:----TH3dXP ) PHlP=--PH.t < IXP ) PH2P=^PH2( IXP)

50
o

o

o

o

oo

190

C

XP=RADMi-«XRl < IXP)

YP=RADM1*XR2( IXP)

ZP=RADm-»XR3( IXP)

SNT 1 P=-TH 1 P-ft-C I TH-PH 1 P #C I PH-XP »C I NO

SNT2P=-TH2P*C I TH-PH2P*C I PH-YP#C I NO

SNT3P=-TH3P*C I TH-ZP*C I NO

SNT1=SNT1~SINT1 < IXP)+SNT1P

SNT2=SNT2-S I NT2 (IXP) +SNT2P

SNT3=SNT3-S I NT3 (IXP) +SNT3P

END IF

CONTINUE

STORE CORRECTED VALUES OF INTEGRAL

IF( IX. EQ. 1 ) THEN

CJTHA( 1 )=SNT2

CJPHA( 1 )=-SNTl

ELSE

C JTHA (IX) =-PH 1 ( I X > -tf-SNT 1 -PH2 (IX) »SNT2

C JPHA ( I X ) =TH 1 ( I X ) *SNT 1 +TH2 (IX) -K-SNT2+TH3 (IX) *SNT3

END IF

CORRECT FOR SELF-PATCH CONTRIBUTIONS IF REQUIRED

IF( ICORR. EQ. 1. OR. ICORR. EQ. 3) THEN

DO BOTTOM CAP (THETA = PI)

IF( IX. EQ. 1 ) THEN

C J3D2=TH3 ( 2 ) *C JTHD ( 2

)

N1=NL( 1

)

N2=NR ( 1 CJ3DL=TH3(N1 )*CJTHD(N1

)

C J3DR=TH3 ( N2 ) *CJTHD ( N2

)

DXCJ3=CJ3D2-(CJ3DL*FC1 ( 1 )+CJ3DR*FC2( 1 )

)

N1=NU1 ( 1

)

N2==NU2 ( 1 CJ3PU=TH3(N1 )*CJTHD(N1 )*FU1 ( 1 )+TH3(N2)*CJTHD(N2)*FU2( 1

)

N1=ND1 ( 1

)

N2=ND2 ( 1 CJ3PD=TH3(N1 )*CJTHD(N1 )*FD1 ( 1 ) +TH3 ( N2 ) *C JTHD ( N2 ) »FD2 ( 1

)

DYCJ3=CJ3PU-CJ3PD

DJ3D2=TH3 ( 2 ) *DJTHD ( 2

)

N1=NL( 1

)

N2=NR ( 1 DJ3DL=TH3 ( N 1 ) *DJTHD ( N 1

)

DJ3DR=TH3 ( N2 ) *DJTHD ( N2

)

DXDJ3=DJ3D2-(DJ3DL*FC1 ( 1 ) +DJ3DR*FC2 ( 1 )

)

N1=NU1 ( 1

)

N2=NU2( 1 DJ3PU=TH3 ( N 1 ) *DJTHD ( N 1 ) «FU 1(1) +TH3 ( N2 ) *DJTHD ( N2 ) »FU2 ( 1

)

N1=ND1 ( 1

)

N2=ND2( 1 DJ3PD=TH3 ( N 1 ) *DJTHD ( N 1 ) *FD 1(1) +TH3 ( N2 ) *DJTHD ( N2 ) *FD2 ( 1

)

DYDJ3=DJ3PU-DJ3PD

PC=CO( 1

)

51
PD^-=DO( 1 ) PDD=DO(IS)

CJTHA( 1 )=CJTHAa > ^«-PC+DJTHD ( 1 ) 1 25*DXCJ3-PDD^^DXDJ3

CJPHA( 1 )=CJPHA< 1 )<^PC+DJPHD( 1 )^tPD-. 1 25^^DYC J3-PDD^f•DYDJ3

ELSE IF (IX. NE, IS) THEN

C

C DO F^EGULAR PATCF^ES

C

C JTHA < I X > =C JTHA ( I X )

( 1 . -CO (IX)) -DJTHD (IX) «DO (IX)

CJPHA( IX )=CJPHA( IX) «•( i. +CO(IX) ) +DJPHD ( I X ) *DO ( I X ) EI.SE

P

C DO TOP CAP IF RE(3UIRED

C

CJ3D2=TH3( ISl )*CJTHD( ISl ) N1=NL( IS)

N2=^NR (IS)

C J3DL=^TH3 ( N1 ) »C JTHD ( N 1 ) C J3DR=TH3 ( M2 ) *C JTHD ( N2

)

DXCJ3=^-CJ3D2-(CJ3DL»-FC1 (2)+CJ3DR*FC2(2) ) N1==NU1 ( IS)

N2=--NU2 (IS)

CJ3PU=^TH3(N1 )-«-CJTHD(Nl ) »-FUl ( IS ) +TH3 ( N2 ) *C JTHD ( N2 ) *FU2 ( IS)

N1=^ND1 ( IS)

N2=ND2( IS)

CJ3PD:“-=TH3(N1 )^^CJTHD(N1 )*FD1 ( I S ) +TH3 ( N2 ) *C JTHD ( N2 ) *FD2 ( IS)

DYCJ3=^CvF3PU-CJ3PD

DJ3D2'--^TH3( ISl )*DJTFHD( ISl ) N1=NL( IS)

N2==NR (IS)

DJ3DL:^TH3(N1 )-«-DJTHD(Nl ) D J3DR---TH3 ( N2 ) JTHD ( N2 ) DXDJ3-- DJ3D2- ( DJ3DL*FC 1(2) +DJ3DR-»^FC2 ( 2 ) ) N1==NU1 ( IS)

N2^==NU2 (IS)

DJ3PU^TH3(N1 )^^DJTHD(N1 )-»-FU,t ( IS ) + TH3 ( N2 ) *DJTHD ( N2 ) -B-FU2 ( IS)

Nl^-NDl (IS)

N2^=ND2( IS)

DJ3PD:--=TH3(N1 )^^-DJTHD(Nl )-«-FDl ( IS ) +TH3 ( N2 ) *DJTHD ( N2 ) *FD2 ( IS)

DYDJ3==DJ3PU-DJ3PD

CJTHA( IS)=CJTHA( I S ) -tt-PC+D JTHD ( IS)^«-PD+. 1 25*DXC J3+PDDi^DXDJ3

CJPHA( IS)=CJPHA( IS)*PC4-DJPHD( IS)^^PD+. 1 25*DYC J3+PDD^^DYDJ3

END IF

END IF

200 CONTINUE

210 IXH=IX-1

C STORE NEW VALUES OF FIELDS

CJTH( 1, NUi IS)=CJTHA( 1; IS)

CJPHC 1, NU; IS)=CJPHA( 1; IS)

IF(NUM. EQ. NU) THEN

C WRITE RESULTS OUT TO DISK IF ARRAY IS ALL NEW

DO 2S0 1=1, NU

WRIT£(2, REC=NREC ) CJTH(l,IiIS)

WRITE(2, REC=NREC + 1 ) CJPHd, li IS)

52
NREC^NREC+2

280 CONTINUE

NUM=^0

END IF

C INCREMENT TIME

T=T+DT

300 CONTINUE

TL=SECOND( )-TL MAIN

PRINT 'EXECUTION TIME FOR THE LOOP IS ',TL, ' SECONDS'

C WRITE REMAINING RESULTS OUT TO DISK

DO 310 I=NU-NUM+1, NU

WR J.TE<2, REC-=NREC) CJTH<1,I;IS)

WRITE (2, REC=NREC+1 ) CJPH( 1, li IS)

NREC=NREC+2

310 CONTINUE

C MAKE DISK FILE PERMANENT

CLOSE < 2)

CALL Q5PURGE( 'LFN=- 'SC3V'//QAL, 'STATUS= ', ISTAT

)

IFdSTAT. GT. 0. AND. ISTAT. NE. 1402) STOP 4

CALL Q5DEFINE( 'LFN=', 'SCSV'//QAD

OPEN (2, FILE='SC3V'//QAL, STATUS=- 'OLD ', ACCESS^ 'DIRECT '

,

1

FORM^ ' UNFORMATTED ' , RECL=NF2

)

400 READ(*, 4, END^-500) THSC> PHSC

ITH=THSC

lPHI:--:-PHSC

WRITE(TITLE, 6) ITH, IPHI

THSC=THSC*CONV

PHSC=PHSC^«-CONV

C COMPUTE COMPONENTS OF SCALED VECTOR TO FIELD POINT

SR 1 :=^S I N ( THSC ) -J^COS ( PHSC ) *CDTM 1 SR2=^S I N ( THSC ) -J^S I N ( PHSC ) *CDTM 1 SR3==C0S(THSC )*CDTM1

C COMPUTE FAR FIELDS

CALL FARFL

GO TO 400

500 CL0SE(9)

C SAVE FILE WITH COMPUTED FIELD INTENSITIES

CALL G5PURGE< 'LFN=', 'SCFL'//QAL, 'STATUS= ISTAT ) IFdSTAT. GT. 0. AND. ISTAT. NE. 1402) STOP 5

CALL Q5DEFINE< 'LFN=- ', 'SCFL'//QAD

STOP

1 FORMAT < 20 ( IX, 15) ) 2 FORMATS IX. 14, 3X, 10(2X, IPEIO. 3)

)

4 FORMAT (El 2. 5)

5 FORMAT (314)

6 FORMAT ( 'PULSE SCATTERED IN THE DIRECTION THETA=

1 ', PHI= ', 15)

7 FORMAT ( 'TIME i«-C IN METERS')

3 FORMAT (' INTENSITY ')

9 r np|viAT( 'LST. SQ, ALPHA=^ '

. F4. 1 , '

, BETA= ', F4. 1 , '

,

i WT--=',I3, '. N1=',I2, ', N2=^',I2, N3=',I2, ', 10 EORMAT( 'INTERP. , ALPHA= ' , F4. 1 , '

, PETA~ '

, F4. 1 , '

,

1 NT=^', 13, ', Nl=- ', 12, ', N2-^-' .. 12, '. N3= '

. 12,

11 FGRMAK LST. SQ,. ALPHA^^ '

, F4. 1 , '

, BETA^-- ' , F4. 1 , '

,

', 14,

DT= ', 1PE7 1,

NO CORR '

)

DT= ', 1PE7. 1,

NO CORR '

)

DT-', 1PE7. 1.
12

13

14

15

16

C

C

c

c

c

c

1 13, ^ Nl-»', 12, S N2==S 12. ^ N3= ^ 12, ' FORMAK 'INTERP. , AL.PHA= F4. 1, ^ BETA=^',F4. 1, '

1 NT==',I3, Nl^', 12, N2-^',I2, N3=',I2, ' FORMAK'LST. SQ, ALPHATM F4. 1

,

BETA=^F4. 1,'

1 NT==',.T3, Nl==^12, N2TM',12, N3=:-SI2, '' FORMAT ( 'INTERP. , ALPHATM', F4. 1, ', BETA=^',F4. 1, '

1 NT=-^', 13, ', Ni=',I2, ', N2=-',I2, ', N3=--=',I2, ' rORMAT< 'LST. SQ, ALPHA= '

, F4. 1 , '

, BETA'^ ' , F4. 1 , '

1 ', N1=',I2, ', N2=^',I2, ', N3=',.T2, ' FORMAT ( 'INTERP. , ALPHA===SF4. 1, ', DETA= ',F4. 1, '

i NT=--', 13, ', Nl=', 12, ', N2=^', 12, ', N3=', .12, ' f:md

SLF PCH';

DT*', 1PE7. 1,

SLF PCH')

DT-', 1PE7. 1,

NEI PCH')

DT«', 1PE7. 1,

ME I PCH')

DT*', 1PE7. 1,

SLF S< N')

DT=', 1PE7. 1,

SLF ?/ N')

SUBROUTINE AVCUR ( IN, NXl , NX2, FXl , FX2, CJTHfU C-..fFHP,

1 CJTHX.CJPHX)

TO COMPUTE CURRENTS AT DESIRED PHI BY AVERAGING

PARAMETER ( LL=LLG, LC^^LCG ) DIMENSION NXl (LL), NX2(LL), FXl (LL), FX2(LL),

1 CJTHX(LL), CJPHX(LL)

M1==NX1 ( IN)

N2:=NX2( IN)

F1=FX1 (IN)

F2-=^F X2 (IN)

IF(N1.EQ. N2) THEN

IF(N1.EQ. 1) THEN

FOR BOTTOM CAP

CJTHP =-F1 «-CJTHX (N1 > ~F2^5-C JPHX ( N 1 ) ELSE

FOR TOP CAP

CJTHP=--F1*CJTHX (N1 )+F2-K-CJPHX(Nl ) END IF

CJPHP==-F2#CJTHX <N1 ) +F 1 *C JPHX ( N 1 ) ELSE

FDR OTHER PATCHES

CJTHP =^-F1 »CJTHX (N1 ) +F2*C JTHX ( N2 ) CUPHP=Fl-«-C JPHX (N1 )+F2-K-CJPHX (N2>

END IF

RETURN

END

SUBROUIINE FARFL

C THIS SUBROUTINE COMPUTES THE RADIATION FIELDS AT TIMES IN T1

C IN THE DIRECTION GIVEN BY SRI, SR2, SR3.

C TMP: TIME ARRAY

C AMP: AMPLITUDE ARRAY

r PARAMETER (LL=LLO, LC=LCG, LL12=^LL12G, LL20=LL20G,

1 LB==LLit*2/64+l, LB l=LL/64-«-l ) D I MENS I ON TMP ( 300 ) , AMP ( 300

)

COMMON /RCOM/ XR 1 < LL ) , XR2 ( LL > , XR3 < LL ) , SRF

(

LL ) , TF < LL ) .

1 CO(LL), DO(LL), THl (LL), TH2(LL), TH3(LL),

2 PHI (LL), PH2(LL), FUl (LL), FU2(LL), FDl (LL), FD2(LL). FCl (2), FC2(2),

54
3 THETA (LL), PHI I (LL) , DPHI I< LL ) , RB 1 < LB 1 ) , RB2<LB1 ) > RB3(LB) , RB4(LB ),

4 C0J<6, LL12), C0JT(6, LL12), C0vJU(2, LL12), C0JTU(2, LL12),

5 C0JV(2, LL12), C0JTV(2. LL12), SINTKLL), SINT2(LL), SINT3(LL)>

6 DR 1 ( LL) . DR2 < LL) , DR3 ( LL) , DR (LL) . CDRM2 < LL) , DRM3 ( LL) , 7 C JTHD ( LL) , C JPHD ( LL) > DJTHD ( LL) , DJPHD ( LL)

8 UTH(LL), UPHELD, AO < LL) > A1 ( LL) , A2 < LL) , A3 ( LL) , A4 ( LL)

,

9 YM2<LL), YMl (LL), YO<LL), Y1 (LL ), Y2<LL),

A C JTHA ( LL ) , C JPHA ( LL ) , C JTH ( LL, LC ) , CJPH ( LL, LC ) , FT ( LL )

,

B NL(LL), NR(LL), NUl (LL), NU2(LL), NDl (LL), ND2(LL),

C NEI (LL, LL20), NPGINT(LL, LL20) , IND(LL), IJ(LL)

COMMON /ICOM/ IX, IS, ISl, Nl, N2, NI, ICORR, IPLOT, INTER, IDEB,

1 IPATCH, N3, NU

CHARACTER TITLE*80, SUBT*80, XLAB*30, YLAB*30

COMMON /LABEL/ TITLE, SUBT, XLAB, YLAB

COMMON /PARAM/ ALPHA, BETA, DT, RAD, RADMl, CDTMl, THSC, SRI, SR2, SR3

COMMON /CONST/ AMUO, TWOMUO, TWOEPO, C, CMl, PI, TOPI, PITO, OS, OTW, OTF

PRINT *, 'START FAR-FIELD COMPUTATION'

CDTM2=CDTMl*i^-2

RCDTMi=--RAD^J-CDTMl

FAC=^(1. E-7*T0PI*RADM1 )**2

NU=1

NI2=NI*2

READ(-», 1 ) Tl, T2

COMPUTE FIRST TIME OF ARRIVAL IF Tl = 0

IF(T1. EQ. 0. ) T1=RAD*CM1*( 1. -2. *SIN(. 5*THSC ) ) -DT

COMPUTE LAST TIME IF T2 = O

IF(T2. EQ. 0. ) T2=-(NI -5)*DT-RAD#CM1

T1DT--^T1/DT

READ(*,2) N

IF(N. GT. 300) STOP 6

AMP( i; 300)=— 1. DT1 = (T2-T1 ) / ( (M-1 )-«-DT)

NMAX---T 1 DT+RCDTM 1 +4.

NMIN=T lDT-RCDTMi-3.

NTOT=^NMAX-NMIN+i

IF(NTOT, GT. LC ) THEN

PRINT 'NTOT= ', NTOT, ' GREATER THAN LC~ '

. LC

STOP 7

END IF

IF(NMIN. LE. 0) THEN

DO 130 Jl=^l, -NMIN+1

CJTH( 1, Jli IS)=0.

CJPH( 1, Jl; IS)^0.

.130 CONTINUE

WREC-0

ELSE

NREC=(NMIN-1 )*2

Jl=^l

END IF

DO 150 J=J1,NT0T

NREC=NREC+2

IF(NREC. GT. NI2) THEN

PRINT «•, 'TIME LIMIT ON SPHERE EXCEEDED DURING INTEGRATION'

PRINT 'T2= ',T2, ', T== ',T1DT*DT, '. K=- ',K, ' OUT OF ',N

55
STOP S

END IF

READ(2, REC=NREC-1 ) CJTHd, J; IS)

READ(2, REC-NREC) CJPHd, J; IS)

150 CONTINUE

DO 200 K=1,N

NX==TlDT+RCDTm4-4.

NN~TlDT-RCDTMl-3.

IF(NX~NN+1. GT. LC) THEN

PRINT 'NEW RANGE SNX-NN+l, ' GREATER THAN LC=',LC

PRINT *, 'NN=',NN, S NX^=^SNX, K-',K

STOP 9

END IF

IF<NX-NMIN4-1. GT. LX ) THEN

NA^NN-NMIN

DO 160 L-1, NMAX-NN+1

MA=NA+1

CJTH( 1, L; IS)=CJTH< 1, NA; IS)

CJPHd. L> IS)=CJPHd, NA; IS)

160 CONTINUE

NNIN==NN

ELSE

L^NMAX-NNINf-2

END IF

IF(NX. GT. NMAX ) THEN

DO 170 J==^L, LC

NREC=NREC+2

IF(NREC. GT. NI2) THEN

PRINT «-> 'TIME LIMIT ON SPHERE EXCEEDED DURING INTEGRATION'

PRINT 'T2= ',T2, T= ',T1DT»-DT, K= '.K. ' OUT OF ',N

STOP 10

END IF

READ (2, REC=NREC-1 ) CJTH< 1, J; IS)

READ(2, REC=NREC) CJPHd,J;IS)

170 CONTINUE

MMAX=^-TviMIN4-L.C~l

END IF

NMINl^NMIN-1

TMP (K)--T1DT^^DT»^C

DR d; IS)=-T1DT + XR1 ( 1; I S ) »^SR 1 -*-XR2 < li IS ) -»-SR2+XR3 ( 1 ; IS)»SR3

I ND ( 1 I S

)

I NT ( DR d i IS); I ND d ; IS))

FTd; IS)=DR d> IS) -VAINT(DR < 1; IS); IS)

INDd; IS) = I Jd, IS) + dNDd; IS) -NMINl )*LL

IF (INTER. EQ. 0) THEN

CALL. INTRPO( IND, FT, C JTH, CJTHD, DJTHD, AO, Al, A2.

1 YM2, YMl, YO, Yl, Y2)

CALL INTRPO( IND, FT, CJPH, CJPHD, DJPHD. AO, Al, A2,

1 YM2, YMl, YO, Yl, Y2)

ELSE

CALL INTRPl ( IND, FT, CJTH, CJTHD, DJTHD, AO, Al, A2, A3, A4,

1 YM2, YMl, YO, Yl, Y2)

CALL INTRPl ( IND, FT, CJPH, CJPHD, DJPHD, AO, Al, A2, A3, A4,

YM2, YMl, YO, Yl, Y2)

END IF

56

1
o

o

o

C1=Q8SSUM( <DJTHDU.; IS)*TH1 (1; IS)+DJPHD(li IS)^^PH1 (1; IS) )

1 *SRF< li IB)

)

C2=^QSSSUM( <DJTHD( 1; IS)^^TH2(1; IS)+DJPHD(1; IS)»-PH2(li IS) ) 1 ^tSRF< 1; IS) ) C3-^QSSSUM(DJTHD( 1; IS)*TH3U; IS)*SRF< 1; IS) ) AMP ( K ) = ( ( C 1 *^»•2+C2fr*2+C3^^*2 ) i^CDTM2- ( C 1 *SR 1 +C2-JtSR2

1 +C3^^SR3)**2)*FAC

T1DT=T1DT+DT1

200 CONTINUE

WRITE (9) IMP

WRITE < 9) AMP

IFdPLOT. EQ. 1. OR. IPLOT. EQ 2)

1 CALL DRAWl (1, 16, 12, 56, 00, N, 0. , 0. , XLAB, YLAB. TITLE, SUBT, TMP, AMP

)

IF( IPLOT. EQ. 0. OR. IPLOT. EQ. 2)

1 CALL DRAW2< 1, 16, 12, 56, SO, N, 0. , 0. , XLAB, YLAB, TITLE, SUBT, TMP, AMP)

RETURN

1 FORMAT < El 2. 5)

2 FORMAT (21 4)

END

C

SUBROUTINE INFL(SINT1, SINT2, SINT3, X, Y, 2, T, *)

C

C THIS SUBROUTINE PROVIDES THE INCIDENT FIELDS

r COMMON /PARAM/ ALPHA, BETA, DT, RAD, RADMl, CDTMl, THSC, SRI, SR2, SR3

COMMON /CONST./ AMUO, TWOMUO, TWOEPO, C, CMl, PI, TOPI, PI TO, OS, OTW, OTF

2CT=^2-C^fT+RAD

IF( 2CT. GE. 0. ) RETURN 1 AZ:=ALPHA^<-2CT

BZ=BETA*ZCT

IF(BZ. GE, -25. ) THEN

Bl-^-BZ*(EXP(AZ)-EXP(BZ) ) E2=^-B1#C

ELSE

BL“^0.

£2^0 . END IF

SINT1=TW0MU0KD1

SINT2==0.

SI NT 3=0.

RETURN

END

C

SUBROUTINE INTRP0< IND, FT, ARR, CA, DA, AO, Al, A2,

i YM2, YMl, YO, Yl, Y2)

LEAST SQUARES INTERPOLATION

PARAMETER ( LL=LLG, LC=LCG

)

COMMON /ICOM/ IX, IS, ISl, Nl, N2, NI, ICORR, IPLOT, INTER, IDEB,

1 IPATCH, N3, NU

COMMON /CONST/ AMUO, TWOMUO, TWOEPO, C, CMl, PI, TOPI, PITO, OS, OTW. OTF

DIMENSION ARR(LL, LC ) , CA(LL), DA(LL), YM2(LL), YMl <LL) , YO(LL),

1 Y 1 ( LL ) , Y2 ( LL ) , AO ( LL ) , A 1 ( LL ) , A2 ( LL ) , I ND ( LL ) , FT ( LL

)

57
Y2( 1; IS)=Q8VGATHR(ARR(1, NU+2i 1 ), IND< 1; IS); Y2< 1; IS)

)

Y1 ( 1; IS)==QSVGATHR(ARR( 1, NU+l; 1 ), IND< 1; IS); Y1 (1; IS) ) Y0( 1; IS)=QSVGATHR(ARR( 1, NU; 1 ), INDd; IS); Y0< 1; IS) ) YMl ( 1; IS)=-~Q8VGATHR(ARR< 1, NU-1; 1 ), INDd; IS); YMl d; IS) ) YM2d; IS)-Q8VGATHR<ARRd, NU-2; 1 )/ INDd; IS); YM2d; IS) ) AOd; IS) = (YM2d; IS)+YM1 d; IS)+YOd; IS)+Y1 d; IS)+Y2d; IS) )*. 2

Aid; IS)-(Y2d; IS)+. 5^^(Y1 d; IS) -YMl d; IS) )-YM2d; IS) )*. 2

A2d; IS) = (Y2d; IS)-. 5*(Y1 d; IS)+YM1 d; IS) )-YOd; IS)+YM2d; IS) )*OS

CAd; IS)=AOd; IS)+FTd; IS)^^A1<1; IS)

1 +<FTd; IS)*-K^2-2. )*A2d; IS)

DAd; IS)*A1 <1; IS)+FTd; IS) *2. *A2d; IS)

RETURN

END

C

SUBROUTINE INTRPl ( IND, FT, ARR, CA, DA, AO, Al, A2, A3, A4,

1 YM2, YMl. YO, Yl, Y2)

C

C POLYNOMIAL INTERPOLATION

C

PARAMETER < LL’^LLG, LC=LCG > COMMON /ICOM/ IX, IS, ISl, Nl, N2, NI.. ICORR, IPLOT, INTER, IDEB,

1 IPATCH, N3, NU

COMMON /CONST/ AMUO, TWOMUO, TWOEPO, C, CMl, PI, TOPI, PI TO, OS, OTW, OTF

DIMENSION ARR(LL, LC), CA<LL), DA(LL>, YM2(LL), YMl (LL), YO(LL),

1 Yl (LL), Y2<LL), AO<LL), AKLL), A2(LL), A3(LL). A4(LL),

2 IND(LL), FT(L.L)

Y2d; IS)=Q8VGATHR(ARRd, NU+2; 1 ), INDd; IS); Y2d; IS) ) Yl d; IS)==Q8VGATHR(ARRd, NU+l; 1 ), INDd; IS); Yl d; IS) ) YO( 1; IS)=QSVGATHR(ARRd, NU; 1 ). INDd; IS); YOd; IS) ) YMl d; IS)==Q8VGATHR(ARR(1, NU-1; 1 ), INDd; IS); YMl d; IS) ) YM2d; IS)==Q8VGATHR(ARRd, NU-2; 1 ), INDd; IS); YM2d; IS) ) C AOd; IS)=YO( 1; IS)

Aid; IS)=OTWK-( YM2d; IS)-Y2d; IS) -8 * ( YMl d ; IS ) -Yl (1 ; I S ) ) ) A2d; IS)=0TF*(-YM2d; IS)~Y2d; IS) + 16. « ( YMl d ; I S ) +Y 1 (1 ; IS) )

1 -30. -^fYOdi IS) ) A3( 1, IS)=0TW-K-(-YM2d; IS)+Y2( 1; IS) +2. *( YMl ( 1; IS)-Y1 d, IS) ) ) A4( 1 ; IS)=OTF*( YM2( 1; IS)+Y2d; IS) -4. *(YM1 d; IS)+Y1 ( 1; IS) )

i +6. *Y0( 1; IS)

)

CAd, IS)=YOd; IS)+FTd; IS)*(Ald; IS)+FTd; IS)

1 » ( A2 d ; I S ) +FT d ; I

S

) « ( A3 d ; IS) +FT (1 ; IS) -«-A4 d ; IS))))

DAd; IS)==A1 d; IS)+FTd; IS)*(2. *A2d; IS)+FT(1; IS)

1 ^( 3 . *A3d; IS)+FTd; IS) *4. #A4d; IS) ) ) RETURN

END

SUBROUTINE INTRP2(TT, ARR, I, NTP, DA)

C

C LEAST-SQUARES INTERPOLATION, SCALAR VERSION, DERIVATIVE ONLY

r PARAMETER (LL ^^LLG, LC=LCG)

DIMENSION ARR(LL, LC)

COMMON /CONST/ AMUO, TWOMUO, TWOEPO, C, CMl, PI, TOPI, PI TO, OS. OTW. OTF

YM2=ARRd, NTP-2)

YM1=--^ARR (I, NTP-1 ) YO^ARRd, NTP)

58
o

o

o

Y1=ARR( I> NTP+i

)

Y2=ARR( I. NTP+2)

C A0== ( YM2+YM 1 +YO+Y 1 +Y2 ) *. 2

A1 = (Y2+. 5->^(Yl-YMl )-YM2)*. 2

A2=<Y2-. 5*<Y14-YM1 )-Y0+YM2)*0S

C CA=A0+TT*Al + (TT*s^2-2. )*A2

DA=AH-TH^2. *A2

RETURN

END

C

SUBROUTINE INTRP3(TT, ARR, I, NTP, DA)

POLYNOMIAL INTERPOLATION, SCALAR VERSION, DERIVATIVE ONLY

PARAMETER ( LL==LLG, LC=LCG

)

DIMENSION ARR(LL, LC)

COMMON /CONST/ AMUO. TWOMUO, TWOEPO, C, CMl , P I , TOP I , P I TO, OS, OTW, DTP

YM2=ARR( I, NTP-2)

YM1=ARR( I, NTP-1

)

YO=ARR< I, NTP)

Y.l=ARR( I, NTP+1 ) Y2=ARR( I, NTP+2)

Al=OTW*( YM2-Y2-8. -fr<YMl-Yl )

)

A2=0TF*(-YM2-Y2+16. YMl+Yl ) -30. *Y0)

A3=0TW^^(-YM2+Y2+2. *(YM1-Y1 ) ) A4=0TF-k-( YM2+Y2-4. *(YMl+Yl)+6. -«-Y0)

C CA=YO+FT*< A1+TT*(A2+TT*(A3+TT#A4) )

)

DA=Al-i-TT<^(2. ->tA2-»-TT^t<3. *A3+TT^^4. *A4 ) ) RETURN

END

SUBROUTINE RUNT

C

C THIS SUBROUTINE COMPUTES COMPONENTS OF THE RADIUS VECTOR

C AND FACTORS NEEDED FOR CORRECTIONS

C

C N1 -NUMBER OF STRIPS ON THE SPHERE (VALUES OF THETA)

C N2:=MAXIMUM NUMBER OF PATCHES IN A STRIP (VALUES OF PHI ON EQUATOR)

C N3=MINIMUM NUMBER OF PATCHES IN A STRIP (NEXT TO A POLAR CAP)

C NPHI IS THE ARRAY WHERE THE NUMBER OF PATCHES IN A STRIP IS STORED

C IS=TOTAL NUMBER OF ANGLE PAIRS

C

PARAMETER (LL=LLG, LC=LCG, LL12=LL120, LL20=LL20G,

1

LB=^LL**2/64+l, LBl=LL/64+l ) DIMENSION NPHI ( 100)

BIT BT3(LL, LL), BT4(LL, LL)

EXTERNAL F003, F103. F013. F203, FU3, F023

EXTERNAL F002, F102, F012, F202, FI 12, F022

COMMON /RCOM/ XR 1 ( LL ) , XR2 ( LL ) , XR3 ( LL ) , SRF ( LL ) , TF ( LL )

,

1 CO(LL), DO(LL). THl (LL), 7H2(LL), TH3(LL),

2 PHI (LL), PH2(LL), FUl (LL), FU2(LL), FDl (LL), FD2(LL), FCl (2), FC2(2),

3 THETA (LL), PHII (LL), DPHII (LL), RBI (LBl ), RB2(LB1 ), RB3(LB), RB4(LB),

4 C0J(6, LL12), C0JT(6, LL12), C0JU(2, LL12), C0JTU(2, LL12),

5 C0JV(2, LL12), C0JTV(2, LL12 ) , SINTl ( LL ) , SINT2(LL), SINT3(LL),

6 DR 1 ( LL ) , DR2 ( LL ) , DR3 ( LL ) , DR ( LL ) , CDRM2 ( LL ) , DRM3 ( LL )

,

59
7 C JTHD ( LL) , C JPHD

LL) , DJTHD < LL) , DJPHD < LL)

.

3 UTH<LL), UPH(LL), AO ( LL ) , A1 < LL) , A2 ( LL) , A3 ( LL) , A4 < LL )

.

9 YM2<LL), YMl <LL)/ YO ( LL ) , Yl( LL ) > Y2 ( LL )

,

A C JTHA < LL ) , C JPHA ( LL ) > C JTH < LL, LC ) , CJPH < LL, LC ) , FT ( LL )

,

B NL(LL), NR(LL), NUl (LL), NU2(LL), NDl (LL), ND2(LL),

C NEI (LL, LL20), NPOINKLL, LL20), IND(LL), IJ(LL)

COMMON /ICOM/ IX, IS, ISl, Nl, N2, NI, ICORR, IPLOT, INTER, IDEB,

1 IPATCH, N3, NU

COMMON /PARAM/ ALPHA, BETA, DT, RAD, RADMl , CDTMl , THSC, SR 1 , SR2, SR3

COMMON /CONST/ AMUO, TWOMUO, TWOEPO, C, CMl, PI, TOPI, PI TO, OS, OTW, OTF

COMMON /PARFUN/ PHI , RHOQ, RHOQI , TRHOQ, TRHOQI , RHOSQ, RHOSQI , 03

EQUIVALENCE (RB3, BT3), (RB4, BT4)

DATA PC, PC2, PC3, ZR0/4*0. /, EPS/1. E-10/

C COMPUTE DELTA-THETA

03=1. /3.

PIM1=1. /PI

TOP I Ml = 1. /TOPI

F0PIM1 = . 54LrOPIMl

RCDTM1=RAD-»^CDTM1

TPCTM1=RCDTM1*T0PIM1

FPCTMl=RCDTMm-FOPIMl

DTH=PI/(N1-1

)

SDTHM1 = 1. /SIN(DTH)

DTH2=DTH*. 5

DTHM1 = 1. /DTH

DTHM1H=F0PIM1/DTH

DTHM1C=DTHM1H*RCDTM1

SDTH2=SIN(DTH2)

SDT2M1 = 1. /SDTH2

C COMPUTE NUMBER OF PATCHES IN STRIP

TH=PI

IS=2

NPHI (

1

)=1

NPHI (Nl )=1

N1H=(N1+1 ) /2

DPH=1.

IF( IPATCH. EQ. 0) THEN

C NUMBER OF PATCHES IN A STRIP HAS TO BE A POWER OF 2

NPHX=1

90 NPHX=NPHX*2

1F(NPHX. LT. N3) GO TO 90

DO 100 1=2. NIH

TH=TH-DTH

STH=SIN( TH)

NPH=(N2-N3)^^STH+N3+. 5

IF(NPH. GT. NPHX) NPHX=NPHX<^2

NPHI ( I )=NPHX

IS=IS+NPHX

DPH=AMIN1 (DPH, STH*SIN(PI/NPHX)

)

100 CONTINUE

ELSE

DO 105 1=2, NIH

TH=TH-DTH

STH=STN(TH)

60
o

o

o

on

105

C

1 10

c

c

c

120

NPH=(N2-N3)*BTH+N3+. 5

NPHI ( I )=NPH

IS=IS+NPH

DPH=AMIN1 <DPH, STH#SIN(PI/NPH)

)

CONTINUE

END IF

TO ASSURE SYMMETRY

IF(M0D(N1, 2). EQ. 1 ) N1H=N1H-1

DO 1 10 J=I, Nl-1

NPH=NPHI <N1H)

NPHI < J)=NPH

IS-=IS+NPH

N1H=N1H-1

CONTINUE

CHOOSE AN ADEQUATE (MINIMUM) TIME INTERVAL

DT1 = 1. 6*RAD*CM1*AMIN1 (SDTH2, DPH)

IF(DT. EQ. 0. ) THEN

USE DTI IF NO TIME INCREMENT IS GIVEN

DT^^DTl

ELSE

IF(DT1. LT. DT) THEN

USE DTI IF GIVEN TIME INCREMENT IS TOO LARGE

PRINT 3, DT, DTI

NI=NI*DT/DT1

DT=DT1

END IF

END IF

PRINT THE COMPUTED VALUES

FILE 'OUTPT' IS USED BECAUSE THE OUTPUT FILE TENDS TO FILL UP

IF( IDEB. EQ. 1 ) THEN

LENGTH- I S* 1 33/ 1 024+ 1

0

CALL Q5RQUEST( 'LFN=^ 'OUTPT', 'MXR=', 133, 'LEN=^ LENGTH

)

OPEN <7, FILE= 'OUTPT '

)

WRITE(7, '( IHl )

'

)

WRITE(7,-}f) 'NUMBER OF PATCHES IN STRIPS'

DO 120 1=^1, Nl, 10

WRITE (7. 1 ) NPHI ( I ) , NPHI < I + l ), NPHI (1+2) , NPHI ( 1+3) , NPHI ( 1+4)

,

1 NPHI (1+5), NPHI (1+6), NPHI (1+7), NPHI ( I +-8) , NPHI (1+9)

CONTINUE

WRITE (7, '(//)')

END IF

PRINT 'NUMBER OF PATCHES IS ',IS, ', NUMBER OF ROWS IS ',LL

IF( IS. GT. LL) STOP 11

IS1=^IS-NPH

rn=p I Bl=2. *(RAD«SIN(DTH*. 25) )^t*2

IN=1

COMPUTE COMPONENT OF RADIUS VECTOR, SURFACE ELEMENT, TANGENTIAL UNIT

VECTORS, ETC.

BOTTOM CAP, THETA == PI

XRl ( 1 )=0.

XR2( 1 )=0.

XR3( 1 )=-RAD

SRF( 1 )==S1

61
TF( 1 )=0.

THETA ( 1 ) =^P I PHI I ( 1 )~0.

DPMI I ( .1 ) -TOPI

THl ( 1 )"1.

TH2( 1 )"0.

TH3( 1 ):=0.

PHI < .1 )=0.

PH2(1 )=1.

NP2=.l

NP3-^NPHI (2)

I MB 2=1

INB3=2

DPHI2--=T0PI

DPHI3=TDPI/NP3

IF( ICORR. NE. 0) THEN

C COMPUTE COEFFICIENTS FDR SELF-PATCH CORRECTIONS FOR CAPS

C0< 1 ) = 1, +. 25^5-DTH

D0( 1 )=. 062!5*DTH-»-*2*RCDTMl

D0( IS)-=. 03125^<-DTH^^RCDTM1

CALL UD(PI. DPHI3, NP3, 1. 2, NL, NR, FCl, FC2)

CALL UD(PITO, DPMI 3, NP3, 1, 2 . NUl, NU2, FUJ , FU2)

CALL UD<PIT0«3. , DPHI3, NP3, 1 , 2, ND 1 , ND2, FD 1 , FD2

)

END IF

C REGULAR PATCHES

DO 200 1=2, Ml -1

TH=TH-DTH

STH=SIN(TH)

CTH=COS(TH)

'RSTH=RAD^P3TH

RCTH=RAD-»tCTH

T F I =CM ( R AD+RCTH ). NP1=NP2

NP2=NP3

NP3=NPHI ( I-f 1 ) INB1=INB2

INB2=INB3

INB3=INB2+NP2

DPHI 1=DPHI2

DPHI2=DPHI3

DPHI3=T0PI/NP3

PHI=0.

C

C SOLID ANGLE FOR THE PATCH

C

SOLG=P I M 1 ^^DPH

1

2#STH*SDTH2*R AD^^^^2

C

C PATCH COEFFICIENT

C

IF (ICORR. NF. 0) THEN

C COMPUTE COEFFICIENTS FOR SELF-PATCH CORRECTIONS

STH,DPH=STH-«-DPH 1 2

STHDPH1=1. /STHDPH

DGAMA2=DTH**?+STHDPH**2

62
DGAMMA=:SGRT ( DGAMA2

)

AL 1 =ALOG < ( DGAMMA+DTH ) / ( DGAMMA-DT H )

)

AL2=AL0G< (DGAMMA+STHDPH) / (DGAMMA-STHDPH)

)

PC2=F0P I M 1 -»-STHDPH*AL 1 PC 3=F0P I M i }«-DTH-»-AL2

PC-^PC2-PC3

ATDG=, 5*DGAMA2*ATAN<STHDPH*DTHMi

)

TAB=. 5*DTH*STHDPH

BNPH2=. 5»fPIT0*STHDPH^^^^-2

PC7=FPCTM 1 < TAB+ATDG-3NPH2 ) PCS=FPCTM 1 -» < TAB-ATDG^SNPH2 ) PD=PC7-PCS

END IF

DO 190 J=i,NP2

I 1 N+

1

CPH=--COB(PHI ) SPH^SIN<PHI

)

XRU IN)=RSTH-^CPH

XR2( IN.P^RSTHr^SPH

XR3( IN)^=RCTH

SRF( IN)=SOLG

THETA ( IN)=TH

PHI I ( IN)=PHI

DPHI I ( IM) -=^DPHI2

THl ( IN):=CTH->^^CPH

TH2< IN)=CTH^^-BPH

TH3-; IN)==-STH

PHI ( IN)=-^-SPH

PH2 (IN) ==CPH

TF( IN)=TFI

IF( ICGRR. NE. 0) THEN

CO(IN)=^PC

D0( IN)=PD

IF(J. EG. 1) THEN

NL( IN)=IN+NP2-1

ELSE

NL ( I N ) := I N1

END IF

IF(J. EG. NP2) THEN

NR ( IN)=IN-NP2+1

ELSE

NR ( IN)^IN+1

END IF

CALL UD(PHI, DPHIl, NP 1 , IN, INBl, NDl,. ND2, FDl, FD2)

CALL UD(PHI, DPHI3, NP3, IN, INB3, NUl, NU2, FUl, FU2)

END IF

PHI=PHI+DPHI2

190 CONTINUE

200 CONTINUE

C TOP CAP, THETA = 0

XRl < IS)=^0.

XR2( IS)=0.

XR3( IS)=RAD

SP.F( IS):^S1

63
THETA( IS):^0.

PHI I

(

IS) =0.

DPHU ( IS)=TOP:t

THl ( IS)"1.

TH2( IS)=--^0.

TH3( IS)==--0,

IF< ICORR. NE, 0) THEN

CALL UD (P.T, DPHI2, NP2, IS, IND2, NL, NR, PHI, PH2)

FCl (2)=PH1 ( IS)

FC2(2)=PH2( IS)

CALL UD(PI TO, DPHI2, NP2, IS, INB2, NUl, NU2, FUl, FU2

)

CALL UD(PITG*3, DPHI2, NP2, IS, INB2, NDl, ND2, FDl, FD2

)

END IF

PHI ( IS)=-0.

PH2( ,[S)^='-1.

TF( IS)^::-2. i^CM,1,*RAD

C CONFUTE QUANTITIES RELATED TO TWO PATCHES

DLl = DTH-«-RAD

NMAX^=-0

I N== 1 NPNT^-“-0

NP:-..=0 DO 240 1=^1, IS

X-^^^XRl ( I ) Y=--XR2 ( I ) Z=--^KR3( I ) IH~THETA( I

)

PHI=PHI Id)

STi-L^SINdH)

RSTH=RAD-J>STH

NM=1

IF(NP. EQ. 0) THEN

NPH^^^NPHI (IN)

DPH^TDPI/NPH

DL2--=DPH*RSTH

rMAX=SQRT( 2 )-i^AMAXl (DLl, DL2)

IFd. E(T 1. OR. I. EQ. IS) RNAX=,1. 5->^DTH<^RAD

NP--^NPH

I N--' I N+

1

END IF

NP--=NP-t

DRl (1; 1S) = X-XRI ( 1.; IS)

DR2( ii IS) =-Y-XR2( li IS)

DR3di IS)^Z-XR3( 1; IS)

DR( 1; IS)“V3QRT(DR1 ( li I S

) -i^-t'-;2+DR2 ( li I S

) -«-«-2-^DR3 ( 1; IS)

1 DR d ; I S ) ) DRN3( 1; JS)=CDTNH^DR di IS)

IND( li IS)=VINT (DRN3( 1; IS)i IND( 1, IS)

)

C DETERMINE WHERE SHIFTS IN ITERPOLATION ARE NEEDED

BT3( 1, I; IS)=^INDdi IS). EQ. 2

R T 4 d , I , I S ) I ND d i IS). EQ. 1 IFdCORR, LT. 2) GO TO 230

C COMPUTE COEFFICIENTS FOR NEIGHBORING-PATCH INTEGRALS

DO 220 J=l, IS

64
TF( J. £G. I. DR. DR( J». GT. RMAX) GO TO 220

IF(NM GT. LL20 OR . NPNT. GT. LL 1 2 ) THFN

PRINT NPNT-- ' , NPNT, ' , J-= '

, J

STOP 12

END IF

NE I < I , NM

)

NPNT=NPNT-H

NPOINT( I, Nh)==NPNT

IF(J. EQ. 1. OR. J. EQ. IS) THEN

NEIGHBORING PATCH IS A SPHERICAL CAP

RH0Q=^STH^SDT2M1

TRH0Q^2 *RHDQ

RHnQH=STH>5F0PIMl*SDTHMl

RH0QCT=RH0GH^^RCDTM1

R HOG! 1 = 1. /RHOQ

TRH0QI=2. frRHOOI

RH0SQ=1. f-RH0Q^t*2

RH0SQI = 1. +RH0G.I-S-S-2

RHCTm=STH«-RCDTNi^^TOPIMl

RHCTM1H=. 5*RHCTN1

USE QIDB FROM CMLIB TO DO NUMERICAL INTEGRATIONS

CP003 INCLUDES AN ADDITIONAL FACTOR RHQO

CALL 0tDB<F003, ZRO, TOPI, EPS, CP003, ER, KF, I FLAG)

IF( IFLAG. GT. 3) THEN

PRINT 'PROBLEM IN CP003= ', CP003, : ERROR= ',ER,

1 NO. OF COMPUTATIONS^ KF, ', FL.AG=', IFLAG

STOP

END IF

CALL QIDB (F103, ZRO, TOPI, EPS, CP 103, ER, KF, IFLAG)

IF (IFLAG, GT. 3) THEN

PRINT 'PROBLEM IN CP 1 03= '

, CP 1 03, ' ERROR=',ER,

1

', NO OF COMPUTATIONS^', KF, ', FLAG= ', IFLAG

STOP

END IF

CALL Q1DE(F013, ZRO.. TOPI, EPS, CP013. ER, KF, IFLAG;

IF( IFLAG GT. 3) THEN

PRINT 'PROBLEM IN CP013= ', CP013, ', ERROR= ',ER,

1 ', NO. OF COMPUTATIONS^ ', KF, FL AG= ', IFLAG

STOP

END IF

CP 203 INCLUDES AN ADDITIONAL FACTOR 1/RHOO

CALL Q1DB(F203. ZRO, TOPI, EPS, CP203, ER, KF, IFLAG

IF (IFLAG. GT. 3) THEN

PRINT 'PROBLEM IN CP203= '

, CP203, ', ERROR=',ER,

1 ', NO OF COMPUTATIONS-^', KF, ', FL.AG= ', IFLAG

STOP

END IF

CPI 13 INCLUDES AN ADDITIONAL FACTOR 1/RHOO

CALL QIDB (FI 13, ZRO, TOPI, EPS, CPI 13, ER, KF, IFLAG)

IF( IFL AG. GT. 3) THEN

PRINT 'PROBLEM IN CP 11 3= ', CP 1 1 3, ', ERROR=',ER.

60
1 NO OF COMPUTATIONS=--=S KF, S FLAG--= ^ IFLAG

STOP

END IF

CP023 INCLUDF:,S an additional FACTOR 1/RHOO

CALL QIDB (F023, ZRO, TOPI, EPS, CP023, ER, KF. IFLAG)

IFdFLAG. GT. 3) THEN

PRINT 'PROBLEN IN CP023« '

, CP023, 'i ERROR=',ER,

.1 L. NO OF CONPUTATIONS--=- KF, FLAG=

I

FLAG

STOP

END IF

CALL GtDB<F002, ZRO, TOPI, EPS, CP002, ER, KF , IFLAG)

IFdFLAG. QT. 3) THEN

PRINT «-, 'PROBLEM IN CP002=-' , CP002, '

; ERROR^',ER,

I NO OF COMPUTATIONS= ', KF, FLAG-', IFLAG

Slop

END IF

CP.1.02 INCLUDES AN ADDITIONAL FACTOR 1/RHOO

CALL Q.1DB-:F1G2, ZRO, TOPI, EPS, CP102, ER, KF, IFLAG)

IF( IFLAG. GT. 3) THEN

PRINT 'PROBLEM IN CP .1 02= ', CP .1 02, '

; ERROR- ', F:R ,

i ', NO OF COMPUTATIONS- ', KF, ', FLAG-', IFLAG

STOP

END IF

CP012 INCLUDES AN ADDITIONAL FACTOR 1/RHOO

CALL Q1DB(F012, ZRO, TOPI, EPS, CP012, ER, KF, IFLAG)

IF< IFLAG GT. 3) THEN

PRINT 'PROBLEM IN CPO 12= ', CPOl 2, '

, ERROR- ',ER,

', NO OF COMPUTATIONS- KF, FLAG-', IFLAG

STOP

END IF

CP202 INCLUDES AN ADDITIONAL FACTOR l/RH00^<-«2

CALL G;lDB(F202, ZRO, TOPI, EPS, CP202, ER, KF. IFLAG)

IFdFLAG. GT. 3) THEN

PRINT -ti-, 'PROBLEM IN CP202- ' , CP202, '

; ERROR-', ER,

1 NO OF COMPUTATIONS- ', KF, ', FLAG=',IFIAC

STOP

END IF

CPI 12 INCLUDES AN ADDITIONAL FACTOR 1/RH004»h^2

CALL G1DB(F112, ZRO, TOPI, EPS, CPJ 12, ER, KF, IFLAG)

IFdFLAG. GT. 3) THEN

PRINT 'PROBLEM IN CP 1 1 2- ' , CP d 2, '

, ERROR-', ER,

.. NO OF COMPUTATIONS- ', KF, ', FLAG- ', I FLAG

STOP

END IF

CP022 INCLUDES AN ADDITIONAL FACTOR I/RHOOm-4^2

CALL GIDB (F022, ZRO, TOPI, EPS, CP022. ER, KF, IFLAG)

IFdFLAG GT. 3) THEN

PRINT 'PROBLEM IN CP022- ' , CP022, '

, ERROR- ',ER,

i ', NO OF COMPUTATIONS- ', KF, ', FLAG-'. IFLAG

STOP

END IF

CPH=COS<PHI ) SPH-SIN(PHI ) COJ ( 1 , NPNT ) -TOP I M 1

«• ( CPH-K-CP003- CP 1 03 )

hh
COJ ( 2. NPNT ) =TOP I Ml * < SPH-J^CP003-CP0 1 3 ) COJT( 1, NPNT )=RHCTMH^<CPH*CP002-CP 102)

COJT ( 2. NPNT ) =RHCTMi * ( SPH*CP002-CP01 2

)

IF( J. EQ. 1 ) THEN

COJ ( 3 . NPN r ) =FOP I M 1 -K-S 7 H* ( CP003-CP203--CP023

)

C0JT(3; NPNT)=RHCTM1H-J^STH*<CP002-CP202-CP022>

ELSE

COJ ( 3, NPNT ) =-FOP I Ml-K-STH* ( CP003-CP2C3-CP023 ) COJT ( 3, NPNT ) = -RHCTM 1 H^^STH-*^ ( CP002-CP202-CP022

)

END IF

COJU( 1, NPNT) = RH0QH-k-<CPH-ij-CP103-CP203)

C0JU(2, NPNT)=RH0QH^^(SPH*CP103-CP1 13)

COJVd. NPNT)=RH0QH*<CPH*CP013-CP1 13)

COJJ < 2, NPNT ) =^RHOQH-{?- < SPH-»^CPO 1 3-CP023 ) COJTU< 1, NPNT)=RH0QCT^^(CPH#CF102-CP202)

COJTU ( 2, NPNT ) =^^RHOQCT-f^ ( SPH-«-CP 102--CP 112)

COJTV( 1. NPNT)=^RH0QCT-»-(CPH#CP012--CPl 12)

COJT J (2. NPNT ) =RHOQCT* < SPH-«-CP01 2-CP022 ) ELSE

C OTHER NEIGHBORING PATCHES

THP=--THETA( J)

PHIF--PHI I ( J)

GTHP^^^SINCI HP ) CTHP*---COS ( THP ) CSCTHP=^1, /STHP

C0TTHP=CTHP*C3CTHP

BPHP---SIN(PHIP ) CPi !P-=COS(PHIP ) STCP==STHP-«-CPHP

STSP- ST HP ^SP HP

CTCP=CTHPii-CPHP

CTSP=--CTHPi'cSPHP

DELTH=THF-TH

DELPHI=PHIP-PHI

IF(I. EQ. 1. OR. I. EQ. 13) THEN

DEEP H 1=0.

ELSE.

IF(DEi-PHI. GT PI ) DEl..PHI=DELPHl-TOPI

1F(DE!.PHT LT. -PI) DELPh I =DELPH I +TOP I END IF

STD£LP=STHP-)tDELPHI

THTP=DEI..TH+. 5*DTH

THTP2=THTP-«->^2

TH TP3=THTP2-«-THTP

THTP4=THTP3-«-THTP

THTM=DEI.TH-. 5*DTH

THTM2-=THTM«-»-2

THTM3=THTM2*THTM

THTM4=THTM3*THTM

DPHI=DPHI I ( J)

PHTP=( DELPHI +. S^^DPHI ) i^STHP

PHTP2=PHTP^(•«2

PHTP3=PHTP2-R-PHTP

PHTP4-PHrP3*PHTP

67
PHTM=M DELPHI-, S^^DPHI > *STHP

PHTM2-PHThH<-*2

PHTh3=PHTM2*PHTM

PHTM4-PHTM3*PHTM

GAM 1 T~SGRT ( THTM2+PHTM2

)

GAM1TI = 1. /QAMIT

GAMlTI2--=GAMlT.T*->^2

GAM 1 T 1 3-GAM 1 T 1 2^«•GAM 1 T

I

GAM2T^^SQRT ( THTP2t-PHTM2

)

GAM2Ti ==l. /GAM2T

GaN2T T 2=GAM2T I >t^2

GAM2T .[ 3---GAM2 V 1 2^^GAM2T I GAM3T ^^SGRT ( THTP24-R|f[ P2 ) GAM3T:i^ i. /GAM3T

QAM3T I 2=^-GAM3T .( **2

GAM3T ;i 3^GAM3T 1 2-t^GAM3T I GAM4T-:SGRT ( THTM2+PH f P2 ) GAM4T1--1. /GAM4T

GAM4 f T2^GAM4TI-»-h2

GAM4T :i 3 =GAM4T 1 2-R-GAM4 T I THMAX2==AMAX1 (THTM2. THTP2)*1. E-J2

PHMAX2--^AMAX1 (PH7M2, PHTP2.>-k-1. E-12

IF(THTM2. LT. PHMAX2. AND. PHTP. LT. O ) THEN

C103=--AL0G(PHTM-« (GAM2T+PHTM)/< ( GAM3T-^PHTP ) »-PHTP) ) ELSE IF(THTP2. LT. PHMAX2. AND. PHTP. LT. 0. ) THEN

C103--= AL0G( (GAM4T+PHTP)i<-PHTP/(PHTMi^(GANl T-t PHTM/ ) )

El... BP

C 1 03^AL0G ( ( GAM4T + PHTP ) * ( GAM2T * PHTM ) / ( v GArl3T-4 PHTP ) fr

t (GAhiT^-PHTM) ) ) END IF

IF(PHTM2. LT THMAX2. AND. THTP. LT. O. ) THEN

C013=-AL0G(THTM-k-(GAM4T<-THTM)/ < (GAN3T+THTP )^<THTP > > ELSE IF (PHTP2. LT. THMAX2. AND. THTP. LT 0. ) THEM

C013'^ALGG( (GAM2T4-THTP ) b-THTP / ( THTM* ( GAM 1 T + THTM ) > ) ELSE

C013=^AL0G( <GAM2T^-THTP)*(GAM4T +-THTM)/( <GAN3TfTHTP) >*

•i (GANi T-t-THTM) ,> ) END IF

IF(PHTM,2. LT. THMAX2) THEN

C203-= -PHTM*ABS ( ALOG ( THTP /THTM ) ) 2 5-sPHTP*AL0G< (GAM3T+THTP )*(GAM4T-THTM) /

3 ( ( GAN3T-Tt-n P ) * ( GAM4T +THTM ) ) ) ELSE IF(PHTP2. LT. THMAX2) THEN

C203=-. 5*PHTM-«-ALDG( ( GAN2T -t-^THTP )»>( GAN 1 T -THTM ) /

1 ( ( GAN27 -THTP ) < GAM 1 T+THT M ) ) ) 2 ^-PHTP^^ABS< ALOG (THTP /THTM) ) ELSE

C203=-^-. D-«-PHTM-frALOG< < GAM2T + THTP ) * ( GAM 1 T-TH I M ) /

1 ( (GAN2T-THTP )*(GAM1T+THTM) )

)

2 +. 5^PHTP-»^AL0G( <GAM3T-(-THTP)-»KGAM4T-THTM)/

3 ( ( GAN3T -THTP ) * ( GAM4T+THTM ) )

)

END IF

C

1

1 3=-GAMl T+GAM2T-GAM3T+GAM4T

IF(THTM2. I T. PHMAX2) THEN

68
C023= - T HTM*ABS ( ALGO ( PHTP/PHTM >

)

2 4 . 5*THTP*ALGG ( ( GAM3T4-PHTP > * ( GAM2T-PHTM ) /

3 ( ( GAM3T- PHTP ) * < GAM2T4 PHTM ) )

)

ELSE IF(THTP2. LT. PHMAX2) THEN

C023^-. 5*THTN*ALOG< ( GAM4T +PHTP ) ( GAMIT-PHTM ) /

1 ( (GAN4T-PHTP )*(GAM1T+PHTM) )

)

2 4- THTP*ABS(ALOG(PHTP/PHTN) ) ELSE

C023^~ . S-t^THTN^ALOG ( ( GAM4T4-PHTP ) ( GAMIT-PHTM ) /

1 ( (GAi'-14T-PHTP ) *(GANlT4-PHTIi) ) ) 2 4 . 54^THTP#AL OG ( ( GAM3T4 PHTP ) # < GAM2T-PHTM ) /

3 ( ( GAM3T- PHTP ) ( GAM2 T + PHTM > ) ) END ir

C125=03-« ( (PHTP4^(GAM3TI GAM4TI ) -PhTM-s- ( GAM2T I -GAM i T I ) >

1 4-C103)

C035-^03k->. (THTM4J-(GAM4TI -GANITI ) -TH f P-!?- ( GAM3T I -GAM2T I ) 4-2. *COi3)

C225:---^Q3-s-'. fHTP*(PHTP*GAM3T I -PHTPHJ-GAN2T I )

1 --(THTMiv(FHTP*GAI'14 ri-PHTri*GAMlT 1 ) ) ) C 1 35 = THTM2-«' ( 03-»-THTM2 r- ( GAMl T I 3-GAM4T 1 3 ) -GAN 1 T I

1 4-GAN4TI )-*-PHTP4*03<^(GAN4TI3-GAN3T.T3)

2 (-THTP2 « ( 03*THTP2* ( GAN3T I 3-GAN2TI 3 ) -GAN37 I 3 4-GAN2TI )4PHTM4-fr03«-<GAM2TI3-GANlTI3)

DTHT= fHTP -THIN

DPHT:-^PHTP-PHTN

PHT 1 =^ATAN2 ( PHTM, THIN ) PHT'2:^-ATAN2 ( PHTN, THTP > PH f3=ATAN2 ( PHTP THTP ) PHT4=^ATAN2(PHTP, THTN

)

PHT.12==PHT1-PHT2

PHT23==PHT2-PHT3

PHT34=PHT3-PHT4

PHT41=PHT4-PHT1

IF(PHTM LT. 0. . AND. PHTP. GT. 0. . AND. THTP. LT. O. ) THEN

PHT23 =PHT234-TGP I PHT4 i:=PHT41-T0PI

END IF

ALG .1 2=-ALGG ( GAM 1 T>GAN2T I ) ALG23=^AL0G ( GAN2TR-GAN3T I ) ALG34=ALGG(GAN3T*GAN4TI

)

ALG41=AL0G(GAN4T-»GAN1TI ) C102=“THTPsPHT234-PHTP4S<^LG34-THTM*PHT4 1+PHTN*ALG12

CO 1 2=^-PHTP )*PHT34-TH1 P*ALG23-- PHTM*PHT 1 2-THTN#ALG4

1

C202-=. 5 -k-(DTHT-«-DPHT4PHTN2*PHT12-THTP2k^PHT 23

1 4PHTP2+:PHT34-THTM24iPHf41 ) C.l 12- 5t^(PHTM2-t^ALG12-THTP2-KALG23

1 4 PHT P2-K ALG34“THTN24i ALG4 1 ) €022=--. 5 : DTHT*DPHT-PHTN2*PHT 1 2 < THTP2+^PHT23

1 -PHTP2-i<PHT344-THTM2*PHT41 ) C 1 24 = . ( PHTN3* < GAM2T 1 2-GAN 1 T 1 2 ) -THTP->s-PHT23

1 HHTP24« (PHTN*GAN2^I2-PHTP^^GAM3T 12)

2 4 PHTP3 K- ( GAM4T 1 2-GAN3T 1 2 ) -THTM*PHT4 1 3 h-THTN2« (PHTP*GAM4TI2-PHTM-«-GAMl r 121 ) C034--. PHTM*PHT12-PHTN2*(THTP*GAN2TI2-THTNkGAN1TI2)

69
1 >THTP-«-vPHTM2#0AM2TI2"PHTP2*GAM3TI2-2, -»ALG23)

2 - PHTP*PHT34-PHTP2* ( THTM^GAMATK’-THTP^GAMST 1 2 ) 3 MHTI''|>^(PHTP2*GAN4TI2- PHTM2*GAM1TI2~2. •»-AL.G4l ) ) C224^. 25* ( ~PHTM2*PHT 1 2 + PHTM3* ( THTP*GAM2T 1 2 --THTM*GAM 1 T 1 2

)

1 -THTP2*PHT23+THTP3* < PHTM*GAH2T 1 2-PHTP*GAM3T 1 2

)

2 -PHTP2*PHT34+PHTP3*( rHTM*GAM^TI2-7HTP*GAM3TI2)

3 --THTM2*PHT4 1 +THT M3* ( PHTP *GAM4T .1 2--PHTM*GAM I T 1 2 ) ) Cl 34--=. 25*(PHTM4*(GAM2TI2-GAM1T 12)

.1 i^THTP2*(PHTM2*GAM2TI2-"PHTP2*GAM3TT2-2. *ALG23 ) 2 +PHTP4* ( GAM4T I 2-GAlM;3T 1 2 ) 3 *THTM2*(PHTP2*GAM4TI2-PHTM2*GAMlT12 -2. *ALG'-'V1 ) ) DTHCOT=DELTH*COTTHP

CljETH-ClOS^i < -1. +DTHC0T) + (-C203<-l 5 «»C225 ) *COTTHP

1 -3. *C125*DTHCOT

C0£PH==C013*( -1. +2. *DTHC0T>+CTHP«-DELPH1*C103

1 *C0TTHP*(“2. *C113vl 5*C135)-3. *DTHC0T*CG35

COEIMO=. 5*C203--DE.LTH*C 1 03

C 0El R . 5 » C 0

2

3 * C SC T HP - DEL P H I *C 0 1

3

CTHTH=DELTH*C

1

03-C203

CTHPH==dELTH*CO 1 3 -C .1 1 3

CPffrH-3TDELP*C 1 03-C 1 1

3

C P HP H^C> T DEEP *C 0 1 3 -C 02 3

COJ 1: JTHTH. 2; JTHPH 3 J1 HM. 4-JPHTH, 5 JPHPH 6 JPHN

C O J < 1 . NP MT ) = “TOP I M 1 *C 1 PIP

H

COJ ( 2, NPNT ) =TOP IMl * < COENO+COEF?*S T HP + CTHTH ) COJ< 3 > NPNT)=--T0P .TM1 ( -CDEPH-t-CPHT H*COT T HP ) C0J(4, NPNT >=T0PIMl-*^r(-C0EN0-C0FR>J-5THP- CPHPH)

CGJ( 5. iMPNT)=TOPIMl*CPHTH

CO J <6.. MPNT >=^T0PIMl*(C0ElH-t-C0ER*CTHP4-CPHPH*CnTTHn

COJU' .1, MPNT)=-DTHM1H*C1 HPH

CGJU<2, NPNT )=DTHM1H*CTHTH

STDMIH^EDPIMI / (STHPi^DPtHI ) CGJV( 1.. NPNT)==-STDM1H*CPHPH

CO JV

(

2> NPNT ) =STDM1 H*CPHTH

C0ETH=C102*( -1 +DTHC0T) + ( -C202«-C22^ ?*CCn I HP

.1 --2, *C124*DTHC0T

C0EPH==C012*( -1 -t-2 *DrHCOT )+CTHP*DELPHl + C 102

] -i-CGl T HP* < -2. *C1 l2*C134)-2. « DTHCDl hC034

CQt-.NO=-^. 5*C202--DELTH*C102

COER=. 5*C022*CSCTHP -DELPHI *C012

CTHTH=DELTH*C 102-C202

CTHPH=DELTH*C012-C1 12

CPHTH=^STDELP*C102-C1 12

CPHPH=S rDELP*CO 1 2--C022

COJTd, NPNT)=-TPCTM1*CTHPH

COJT ( 2, NPNT ) =TPCTMi * ( COENO+COER *STHP-t CTH fH ) COJT < 3, NPNT ) =TPCTM1 ft < -COEPH-t-CPHTH*COT IMP ) CUJT (4.. NPNT ) =TPCTM1 * ( -COENO-COER*STHP-CPHPH ) COJT( 5. NPNT)=TPCTM1*CPHTH

C0JT(6, NPNT )=^^TPCTMl ft(COETH+COERftCTHP>CPHPHftCOJTHP > COJTU( 1 , NPNT)=-DTHM1C*CTHPH

CGJTU(2/ NPNT)=-DTHIilC*CTHTH

STDM1C=STDM1H*CRDTM1

COJTV( I, NPNT)=-3TDM1C*CPHPH

70
C0JTV(2, NPNT)=STDM1C#CPHTH

END IF

2S0 CONTINUE

230 NEI(I,. 1)=NN

I F ( NM . QT . NMA X ) NNA X = NM

240 CONTINUE

PRINT ^NUMBER OF COLUMNS NEEDED IS UNMAX-1, ^ OUT OF 19,

1

' TOTAL NUMBER OF NEIGHBORING PATCHES IS U NPNT,

2 U. W I TH A hi A X I MUM OF •- , LL 1

2

RETURN

1 FORMAT^ 10(2X, I 10) )

3

FORMAT GIVEN DT-=UE9. 2, ' GREATER THAN "MINIMUM" DT=^UE9. 2)

END

SUBROUTINE UB(PHI, DPHI, NP, IN, INB, NXl, NX2, FX1 , FX2)

C

C THIS SUBROUTINE IS USED TO COMPUTE QUANTITIES NEEDED FOR PHI -INTERPOLATION

C PARAMETER (LL==L LG, LC-LCG)

D I MENS I ON NX 1 ( LL ) , NX2 < LL ) , F X 1 < LL )

,

FX2 < LL

)

FI -PH I /DPHI

Ml -FI

Fi-Fl-REAL(N.l ) NXl ( IN) = IMB-H41

IFUMl EQ NP- i ) THEN

NX2‘: IN ^ -INB

ELSE

NX2( IN)=INB+N1+1

END IE

IFUMi . EQ. N2) THEN

C FOFv' POLAR CAPS

FXl ( IN) -COS (PHI

)

FX2( IN )- SIN (PHI ) ELSE

FXl ( IN)-1. -FI

FX2( IN) -FI

END IF

RETT'RN

END

r FUNCTION FOOCKX)

common /PARFUN/ PHI , RHOQ, RHOQI, TRHOQ. TRHOQI , RHOSQ, RHOSQI , 03

C03X-C0S( X

)

EKS-SQRT ; RHOSG-TRHOQ-frCOSX ) F003- 1 . / <; EKS-<5 ( EKB vRHOQ-COSX ) ) RETURN

END

r FUNCTION F103(X)

COMMON /PARFUN/ PHI, RHDQ, RHOQI, TRHOQ, TRHOQI, RHOSQ, RHOSQI, 03

CQGX=CGS(X)

RCQGX-RHOQB-COSX

TCGSX-2. -«-RCOSX

EKS-SQRT ( RHOSG-TCOSX

)

71
ECOSX=^EKS#CC)GX

IF<COSX, GT. . 5) THEN

F"103=-^C0S(X+PHI )-5^( ALOG( < RHOQ+RCOSX ) / < EKS-^RCOSX-1 . ) > 4-

1 2. *COSX/ (EKSXEKS+RHOQ-COSX) )-

2 (1 -TCOSX)/(EKS*( 1. -ECOSX-RCOSX) )

)

ELSE IF(CGSX. GT, 5) THEN

Fi03=^C0S<X4PHI )*(ALOG' ( RHGQ+RCOSX ) / < EKS+RCOSX“l . ) ) +

1 2. -frCOSX/ (EKS^^(EKS4 RH0Q~C0SX )

> -- 2 (1. ^-ECOSX--RCOSX)/(EKS^<•( 1. -C0SX**2> ) ) ELSE

Fi03=-C0S( X4PHI )*< ALOG( ( EKS-RCOSX + 1 . ) / <RH0Q-RC05X ) ) +

1 2 )^COSX/ (EK5*(EKS+RH0Q-CC)3X ) )-

2 (1. -TCOSX ) / (EKS^^( i -ECDSX-RCOSX ) ) ) END IF

RETURN

END

F'JNCTKDN F0i3(X)

COMMON /PARFUN/ PHI, RHOQ, RHOQI, TRHOQ, TRHOQI, RHOSG. RHDSGI, 03

COSX^COSC X .) RC0SX=:RH0Q-»^C0'3X

TCDSX=^2. -»-RC.OSX

£KS==SQRT(RHOSG-TCOBX ) ECOSX--EKS*CGSX

IF(COSX. GT. f>) THEN

FO 13=^3 IN (X+-PHI (ALOG< ( RHOQ-^-RCOSX ) / ( EKS+RCOSX - 1 ) ) +

.1. 2. -iiCOSX / ( EKS^^- < EKS+RHOQ-COSX ) ) - 2 (J, "TCOSX)/(EKS-(^( 1. -ECOSX-RCOSX ) ) ) ELSF IFCCOSX. GT - 5) THEN

F0i3=^SIN( X^•PHI ':AL0G< (RHOQ+RCOSX)/ (EKS4RC03X-I. ) )4

1 2, -«^COSX/vEKS«-(EKS+RHOQ~COSX) )-

2 (.1. 4EC0SX-RC0BX)/(EKS*v1. -C0SX*-«2)))

ELBE:

F013=SIN< X+PHI )*( AL0G( ( EKS-RCOSX + 1 . ) / ( RHOQ-RCOSX ) )

1 2, -fs-COSX/ (EKS*<EKS+RHOQ-COSX ) )-

2 (1 -TCOGX ) / ( EKS# ( 1

. -ECDSX-RCOSX ) )

)

END IF

RETURN

END

FUNCTION F203(X)

COMMON /PARFUN/ PHI. RHOQ. RHOQI, TRHOQ. TRHOQI. RHOSQ, RHOSQI. 03

COGX=COS( X

)

RCOSX=RHOQ*CC)SX

TC0SX=2. ^RCOGX

EKS^--SQRT ( RH0GQ-TC03X ) ECOSX=EKS*CGSX

IF(COSX. GT. . 5) THEN

F203-C0S(X+PHI )*-«2-«-(RH0QI*EKS-1. +

1 3. -«-COSX*ALOG( (RHOQ-t-RCOSX)/(EKS+RCOSX-l. ) ) 4- 2 (4. •s-C0SX**2-l . )/ (EKS-k-(EKG4-RH0Q-C0SX) )-

3 2, -«-COSX*( 1. -TC0SX)/<EKS-J^( 1. -ECOGX-RCOGX ) ) ) ELBE IF(COSX. GT. 5) THEN

F203--'-^C0B(X4 PHI ) *^<2* ( RHOQI *EKS- 1 . 4-

72
1 3. -itCOSX*ALOG( (RHOQ-»-RCOSX)/(EKS+RCOSX-l. )) + 2 <4. *C0SX**2~1. )/<EKS*<EKG+RHOQ--COSX) )-

3 2, 5tC0SX* ( 1 . +ECOSX-RCOSX ) / < EKS* < 1 . “C0SX^t*2 ) ) ) ELSE

F203=C0S( X+PHI )*-«2#(RH0QHtEKS-l. +

1 3. H-CDSX*ALOG( (EKS-RC03X + 1. ) / ( RHOQ-RCOSX ) )

2 (4. »^C0SX**2--1. ) / <EKS-?<-<EKS+RHOQ“COSX ) )-

3 2 *COBX*< 1. -TCOSX)/(EKS<^( 1. -ECOSX-RCOSX ) ) ) END IF

RETURN

END

FUNCTION F;ll3(X)

COMMON /PARFUN/ PHI, RHOQ, RHOQI> TRHOQ, TRHOQI, RHOSQ> RHOSQI, 03

CDSX--COS( X ) RCOSX=RHOQ>tCOSX

TCG5X=2. *RCOSX

EKS==SGRT>:RH0SQ -TCOSX ) EC05X=EKS^C0GX

IF^COSX. GT. , 5) THEN

FI 13=C03(X-^PHI )*5IN<X+PHI )*(RH0QI-frEK3“l. +

1 3 -»^COSX*ALOG< <RH0Q-i-RC0SX)/(EKS4-RC0SX-l. )) + 2 (4. ^COSX-s-K-2-1. ) / (EKS*<EKG4-RH0Q“C0SX ) )~

3 2. -«C0SX*(1. -TCOSX)/<EKS*a. “EC0SX-RCQ3X ) ) ) else IF(C0SX. GT. 5) THEN

FI 13=C0S< X+PHl )*SIN< X+PHI ) * ( RHOQI *EKS“ 1 +

1 3. «COSX-«-ALOG( (RHOQ<-RCGSX ) / (EKS+RCOSX-1. )) +

2 (4. i^COSX^-J^O-l ) / (EKS*<EKS+RHOQ-COSX ) )-

3 2. s-COSX f^ ( 1 . +ECOSX-RCOSX ) / < EKG* ( 1 . --C0SX^£”»-2 ) ) ) ELSE

FI 13=C0S( X+PHI )*SIN(X+PHI )#<RH0QI*EKS-1. +

1 3. ^^-COSX-«-ALOG( \EKS-RC0SX + 1. ) / ( RHOQ-RCOGX ) )

2 <4. -t^C0SX«-fr2-i > / (EKS*<EKS+RHOQ--CQSX ) )-

3 2 *COSXh^ ( 1 . -TCOSX ) / < EKS* ( 1 . -ECOSX-RCOSX ) ) ) END IF

RETURN

END

FUNCriON F023<X)

COMMON /PARFUN/ PHI, RHOQ, RHOQI, TRHOQ, TRHOQI, RHOSQ, RHOSQI, 03

COSX=COS( X

)

RC0SX=RH0Q*C03X

TC0SX=2. *RCOGX

EKS=SQRT ( RHOSG-TCOSX

)

EC0SX=EK3*C0GX

IPUCDSX. GT. . 5) THEN

F023=SIN( X+PHI RHOQI #EKB-1. +

1 3. *C03X*AL0G( (RH0Q+RC0SX)/(EK3+RC03X-1. ))+

2 (4, ^^C0SX t^H^2- 1 ) / (EKSXEKSh-RHOQ-COSX ) )--

3 2 «-COSX-)t < 1 . -TCOSX ) / ( EKS* ( 1 . -ECOSX-RCOGX ) ) ) ELSE IF(CGSX. GT. ~. 5) THEN

F023=BIN( X+PHI ) -«-« 2 a- < RHOQ I -«<-EKS- 1 . +

1 3. -«-CGSX-«-ALOG< (RH0Q'^RC0SX)/(EK5+RC0SX-1. )) + 2 (4 ^^C0SXi^«2-- l. ) / (EK5-«-<EKG-i-RH0Q“C0SX ) )-

73
3 2, -«-COSX»^( 1. +EC0SX-RC0SX)/<EKS<^(1. --C0SX*«2)))

ELSE

F023»SIN( X4PHI )*«2s-(RH0GI*EKS“l. +

1 3. •i^COSX'Ji-ALOGv' (EKS-RCOSX + 1. ) / < RHOQ-RCOSX ) ) + 2 <4. •»-C0SX**2-l. )/<EKS»<EKS4-RH0Q-CDSX) )-

3 2. «-CDSX*(l. -TCSX)/<EKG^^•( 1. -ECOSX-RCGSX ) ) ) END IF

RETURN

END

FUNCTION F002(X)

COMMON /F^ARFUN/ PHI, RHOQ, RHOQI, TRHOQ, TRHOQI, RHOSQ. RHOSQI, 03

SINX=SIM(X)

COGX-^COS ( X ) EKS2=RHDSGI~TRH0QI-«-C0SX

I F ( AB5 ( S I NX ) . GE. 1 . E-4 ) THEN

F02=ATAN(SINX/ (RHOQ-COSX) )/SINX

ELSE

DEN=1. /(RHOG-COSX)

F02=DEN-03*S I NX *-fr2*DEN »^-fr3 END IF

F002=. 5« ALOG ( EKS2 ) +C0SX*F02

RETURN

END

FUNCTION F102(X)

COMMON /PARFUN/ PHI , RHOQ, RHOQI , TRHOQ, TRF^OQI , RH05G, RHOSQI , 03

SINX==--^SIN( X ) COSX=-COS( X ) EKS2=^RH0SQI--TRH0QI#C0SX

IF(ABS(SINX). GE. l.E~4> THEN

F02=ATAN(SINX/ (RHOQ--COSX > ) /SI MX

ELBE

DEM=1. / <RF^0Q•~C0SX ) F02=DEN-03->^S I NX*«2*DEN-«-^^3

END IF

F 1 02=^C0S ( X +PH I ) * (RHOQ I +COSX »-AL0G ( EKS2 )

1

- ( 1 . -2. *CQSX-fr-R-2)*F02)

RETURN

END

FUNCTION F012(X)

COMMON /PARFUN/ PHI, RHOQ, RHOQI, TRHOQ, TRHOQI, RHOSQ. RHOSQI, 03

SINX==SIN(X ) COSX=COS(X)

C0SX2==C03X«-«-2

EKS2=RH0SQI-TRH0QI«-C0SX

IF(ABS(3INX). GE. I. E-4> THEN

F02=ATAN(SINX/ (RHOQ-COSX) )/SINX

ELSE

DEN=1, /(RHOQ-COSX)

F02=DEN-03*S I NX*« 2^«-DEN^<-« 3

END IF

FO 1 2=S IN ( X +PHI ) »^ ( RHOQI +COSX»ALOG ( EKS2

)

74
1 -(1. -2. t^C0SX**2)*F02)

RETURN

END

FUNCTION F202<X)

COMMON /PARFUN/ PHI, RHOQ, RHOQI, TRHOQ, TRHOQI, RHOSQ, RHOSQI, 03

SINX:^SIN( X ) COSX^=COS< X ) C0SX2=C0SXi^i^2

EKS2=Rh'0SQ I -TRHOQI *COSX

IF(ABS(SINX). GE 1. E--4) THEN

F02=ATAN(5INX/ (RHOQ-COSX) )/SINX

ELSE

DEN==1. /(RHOQ-COBX)

t~02=DEN'~03'J^SI NX‘fr*2*DEN'i^’'ti‘3

END IF

F202=C0S( X+PHI) K-*2*(RH0QH<-( . 5*RH0QI+2. *C0SX ) 1 +(2. -!^C0SX2- 5)^^AL0G(EK52)-C0SX^^<3. -4, -t^C0SX2 ) *F02 ) RETURN

END

FUNCTION F112(X)

COMMON /PARFUN/ PHI, RHOQ, RHOQI, TRHOQ, TRHOQI, RHOBQ, RHOSQI, 03

SINX=SIN(X)

COSX=COS( X

)

C03X2=C0SXfi^^2

EKS2^--=RH0SQI-TRH0QI*C03X

IF(ABG<SINX ). GE. 1. E-4) THEN

F02=ATAN(SINX/ (RHOQ-COSX ) ) /SINX

ELSE

DEN=1. / (RHOQ -CGSX ) F02=DEN- 03 K-S I NX 2*DEN**3

END IF

F112=C0S(X+-PHI )*SIN(X-('PHI )*(RHOQI*(. 5*RH0QI+2 ^tC0SX.)

1 +(2. *C0SX2-. 5) »<-AL0G<EKS2)-C0SX-ii-(3. -4. *C0SX2)*F02)

RETURN

END

FUNCTION F022(X)

COMMON /PARFUN/ PHI, RHOQ, RHOQI. TRHOQ, TRHOQI, RHOSQ, RHOSQI, 03

BINX=SIN( X

)

C0SX=C0S(X)

C0SX2=C0SX-»-«-2

EKS2=RH0SQ I -TRHOQ I -k-COSX

IF* ABS(SINX ) GE. l.E-4) THEM

F02=ATAN(SINX/ <RH0Q-C0SX ) ) /SINX

ELSE

DEN=1 /(RHOQ-COSX)

F02=DEN- 03-«-S I NX**2^H)EN> fr3 END IF

F022=SIN( X+PHI ) -K-*2-«-< RHOQI tt(

.

5*RH0QI+2. *C0SX )

1 +(2. -i^C0SX2-. 5)-«-AL0G<EKS2)-C0SX«(3. -4. *C0SX2)*F02)

RETURN

END

7D
PROGRAM SC ATP

C

C TO COMPUTE THE SIZES OF ARRAYS FOR PROGRAM SCAT AND RELATED PROGRAMS

C

CHARACTER QAL*4, SUF*3, JCAT«-2, PR I *2

COMMON /PARAM/ RAD> ALPHA, BETA, DT, DTMl . THDIR

COMMON /ICOM/ IX, IS, ISl. NTOT, Nl, N2, N3, NI, ISKIP, IPATCH, ICORR

COMMON /CONST/ C, CMl , P I , TOP I DATA C/3. E8/

CM1 = 1. /C

PI =4. #ATAN( 1. ) TOP I •=2. *PT

CONV,'=PI/lSO

OPEM( 1, FILE^-^ ^SCTQAL STATUS^ 'OLD '

) READ( 1, #) QAL.. ISKIP, NT, IPATCH, PRI, ICORR, NTL, NLPl

PRINT 'INPUT FILE; SCIN'//QAL, ' SKIP= ', ISKIP, ', PATCH= ', IPATCH.

1 CORR= ICORR

PEN < 2, F I LE-=^ ' SC I N ' / /QAL, STATUS^ ' OLD ' ) READ (2,4 .) DT, RAD, ALPHA, BETA

READ (2, 5) NI

READ (2, 5) N1,N2, N3

IF(N,t. GT. 100) STOP 1 IF(N3. EQ. 0) N3=5

C FIND NUMBER OF PATCHES AND NUMBER OF NEIGHBORING PATCHES

CALL RUNTl (LL12, LL20)

NU-2. «RAD-^CM1 «-DTMl +6.

I F ( N I . EG, 0 ) N 1=6. »R AD/ ( C*DT ) IF(NI.LT. 0) NU=1

PRINT 'DT=',DT, ', RAD='.RAD, '. ALPHA= ', ALPHA, '

. BETA=',PE7A,

1 NI=',NI,', N1=',N1,', M2=',N2, ', N3= ' , N3

PRINT -i^, NU, ' ROWS NEEDED'

C FIND DIMENSION OF ARRAY NEEDED FOR FAR-FIELD COMPUTATIONS

400 READ(2, 4,. END=500) THDIR, PHDIR

ITH=THDIR

IPHI=PHDIR

THDIR=THDIRi'rCOMV

P HD I R = P HD I R ON

V

CALL FARFL

I F ( NTOT . GT , M I'M ) NTM=NTGT

GO TO 400

500 LL=I3+2

C ADD 2 BECAUSE COMPUTATIONS ON 855 AND 205 SOMETIMES DISAGREE

LC=MAX0(NU, NTM)4-2

C PRELIMINARY CALCULATIONS ONLY

IF(NI LT. 0) LC=1

C FILE WITH DIRECTIVES FOR THE FULL SCREEN EDITOR

OPEN (5, FILE='SCTFSE', STATUS^ ' NEW '

)

LB=LL#*2/64+l

LBl=LL/64+l

C COMPUTE NUMBER OF LARGE PAGES

NLP= < LL* ( 2#LC+2*LL20+55 ) +20*LL 1 2+2* ( LB+LB 1 ) +65539 ) Z65536

PRINT *. 'LL=',LL, ', LC=',LC, LL12= ', LL12, ', LL20=',LL20

PRINT *, 'COMPUTED NLP=',NLP

IF(LL. GT. 99999. OR. LC. GT. 9999. OR. LL12. GT. 999999. OR LL20 GT 99^)

76
1 STOP 10

WRITE(5> 6 ) LL, LC, LL12.. LL20

C IF NUMBER OF LARGE PAGES IS GIVEN

IF(NLP1. NE. O) WLP=NLP1

IF (NLP. GT. 28. OR. NTL. GT. 599940) STOP

C ASSIGN TIME LIMITS (IF NOT GIVEN) AND PRIORITIES

IF (NLP. LE. 5) THEN

IF (NTL. EQ. O) NTL= 500

IFdMTL. LE. 900) THEN

JCAT= 'P2'

ELSE IF ( NTL. LE. 3600) THEN

JC AT= ^P3'

ELSE

JCAT='P5'

END IF

ELSE IF(NLP. LE. 10) THEN

IF(NTL. EQ. O) NTL=900

IFdMTL. LE. 3600) THEN

v.lCAT= 'P3'

El.SE

JC AT= 'P5

END IF

ELSE IF (NLP. LE. 20) THEN

IF (NTL. EG. 0) NTL=900

IFdMTL. LE. 900) THEN

JC AT= 'P4

ELSE

JCAT= 'P5 •' END IF

ELSE

IFdMTL. EG. 0) NTL=5000

JCAT= 'P5

END IF

IF (PR I. NE. ^XX '

) THEN

IF (PRI. LT. JCAT) THEN

PRINT "ASSIGNED PRIORITY TOO LOW'

PRI NT « . 'PRI --= APR!. ", JCAT^^ ', JCAT

STOP

ELSE

JCAT^PRI

END IF

END IF

0PEW(4, FILE= 'SCTRES ", STATUS^ " NEW '

)

WRITE(4, -K-) 'RESOURCE, TL=", NTL, ",LP=",NLP, " .. JCAT= "

, JC AT

,

1 ", NT--- ", NT, ".

' PRINT «, "PRIORITY ", JCAT, ", TIME LIMIT "

, NTL , ", ",NT. ' STOP

4 FORMAT (E 12 5)

5 FORMAT (314)

6 FORMATk "REPLACE/LLG/ ", 15, "./ALL QUIET",/, 'REPLACE/LCG/

1 "/ALL QUIET ',/, 'REPLACE/LL12G/ ", 16, "/ALL QUIET",/,

2 'REPL.ACE/LL20G/ ", 13, "/ALL QUIET ",/, "QUI T '

) END

r

TAPE "

14 ,

77
o

o

o

o SUBROUTINE RUNT I (LL12> LL20)

THIS SUBROUTINE DETERMINES THE NUMBER OF PATCHES AND THE DISTRIBUTION

NEIGHBORING PATCHES IF ICORR ~ 2 OR 3

DIMENSION NPHI(IOO)

DIMENSION XRK5000), XR2<5000)> XR3<5000), THETA < 5000) , PHI I < 5000)

COMMON /PARAM/ RAD, ALPHA, BETA, DTi DTMl . THDI

R

COMMON /ICOM/ IX, IS, ISl. NTOT, Nl, N2, N3, NI, ISKIP, IPATCH, ICORR

COMMON /CONST/ C, CMl , P I , TOP I

DT!-i=^P I / ( N 1 - 1 ) BDTH2^SIN(DTHfr. 5)

TH=PI

IS=2

NP)iJ ( 1 NPHI (Nl )=1

N1H=(N1 +1 >/2

DPH=-1.

IF ( IPATCH. EG. 0) THEN

NPHX=-:.t 90 NPHX=NPHX->i-2

IFINPHX, LT. N3) GO TO 90

DO 100 1=2, NIH

TH=TH-DTH

5TH=SIN(TH)

NPI-N(N2-N3)-i-STH+N3+. 5

IF(NPH. GT. NPHX) NPHX= NPHX-k-2

NPHI ( I )=NPHX

IS=IS4NPHX

DPH=AMIN1 (DPH, STH*SI N ( P I /NPHX )

)

100 CONTINUE

ELBE

DO 105 1=2, NIH

TH=1H~DTH

STH=3IN(TH)

NP H= ( N2 -N3 ) •i^-STH+ N3+ . 5

NPHI ( I .)=NPH

IS=I'3+NPH

DPH=AMIN1 (DPH, STH-»^-SIN ( P I /NPH ) ) 105 CONTINUE

END IF

I F ( MOD ( N 1 , 2 ) . EQ. 1 ) N 1 H=N 1 H1

DO 1 1 0 J= I , N I -

1

NPH=NPHI (NIH)

NPHI ( J)=NPH

I3=IS+NPH

N1H=N1H-1

110 CONTINUE

DT1=1 . 6*RADKCM1*AMIN1 (SDTH2. DPH)

IF(DT. EQ. 0. > THEN

DT=DT1

ELSE

IF(DT1. LT. DT) THEN

78
PRINT 3, DT.DTl

NI=NI«-DT/DT1

DT=DT1

END IF

END IF

DTM1=1. /DT

PRINT 'NUMBER OF PATCHES IS IS

C RETURN IF NO NEIGHBORING-PATCH CORRECTIONS ARE GOING TO BE MADE

IF( ICORR. LT. 2. OR. ISKIP. EQ, 1 ) THEN

LL12=1

L!.20=l

RETURN

END IF

C CHECK IF ARRAYS ARE BIG ENOUGH TO CALCULATE THE NUMBER OF NEIGHBORING

C PATCHES

IF( IS. GT. 5000) STOP 5

I3i==IS-NPH

TH^^PI

IN=1

XRK 1 )~ 0 . XR2(1 )=0.

XR3< 1 )=-~RAD

THETA ( 1 )=PI

PHI I ( 1 )=0.

NP2=.l

NP3=NPHI <2>

INB2--=1

INB3=2

DPHI2=T0PI

DPHI3=T0PI /NP3

DO 200 I=2,N1-1

TH=TH-DTH

3TH=SIN(TH)

CTH=COS(TH)

RSTH==RAD-«-STH

RCTH==RAD#CTH

NPJ.=NP2

MP2=NP3

NP3=NPHI ( I+l

)

INB1=TNB2

INB2=INB3

INB3=INB2+-NP2

DPHI 1=DPHI2

DPHI2=DPHI3

DPH 1 3==T0'" I /NP 3

PHI=0

DO 190 J=1,NP2

IN=IN+l

CPH=CGS<PHI

)

SPH=SIN(PHI

XRl < IN)=RSTH-«CPH

XR2( IN)=RSTHfrSPH

XR3( IN)=RCTH

THETA < IN)=TH

79
PHI I

(

IN)^PHI

PHI=PHI+DPHI2

IPO CONTINUE

200 CONTINUE

XRl < IS)=0.

XR2( IS)=-"0.

XR3(IS)«RAD

THETA < IS) =0,

PHII < IS)=0.

DLl-~-DTH^fRAD

MMAX==0

IN=-^1 NPNT^--0

NP=0

DO 240 1 = 1, IG

X=XR1 ( I

)

Y=XR2( I

)

Z=XR3( I

)

TH=THETA( I

)

PHI=PHI Id)

NM= 1 :tF (NP.EQ. O) THEN

NPH=NPHI < IN)

DPH=TOPI/NPM

DP HP R = DP N*l. 1 DL2=DPH*BI N ( TH > ->^R AD

RHAX=BQRTf2. ) frAMAXl (DLl, DL2 ) IFd, EQ 1. OR I EQ. IS) RMAX=1. 5*DTH^^RAD

NP=NPH

IN^=IN+1

END IF

NP=NP-1

DO 220 J=1,IS

DR=SQRT< ( X-XRl < J) ) *#2-K Y-XR2 ( J ) ) «-«-2+ ( 2-XR3 ( J ) )*«2)

IF< J. EQ. I. OR. DR. GT. RMAX) GO TO 220

NN=NM+.l

NPNT=NPM1 +1

220 CONTINUE

IF(NM. GT. NMAX ) NMAX=NM

240 CONTINUE

PRINT *, 'NUNBER OF COLUMNS NEEDED IS ',NMAX-1.

1

' TOTAL NUMBER OF NEIGHBORING PATCHES IS ',NPNT

LI..12=NPNT

I..L20=NMAX

RETURN

3 FORMAK' GIVEN DT= ',E9. 2, ' GREATER THAN "MINIMUM" DT=',E9 2)

SUBROUTINE FARFL.

C

COMMON /PARAM/ RAD, ALPHA, BETA, DT, DTM) . THDIR

COMMON /ICOM/ IX, IS, ISl. NTOT, Nl, N2, N3. NI, ISKIP, IPATCH, ICORR

COMMON /CONST/ C, CM 1 , P I , TOP I READ < 2,1) T1,T2

SO
C COMPUTE FIRST TIME. OF ARRIVAL IF T1 = 0

IF(T1. EQ O. ) T1=RAD*CM1*( 1. -2. K-SIN(. S«^THDIR) ) ~DT

C COMPUTE LAST TIME IF T2 = 0

IF<T2. EQ. 0. ) T2=(NI-5>*DT-RAD^<CM1

READ (2, 2) N

MMAX=(TH-RAD*CM1 )/DT+4.

NMIN=(T1--RAD*CM1 >/DT-3.

NT0T=NMAX-NMIN<-1

PRINT 'NTOT=MMTOT

RETURN

1 FORMAT (£12. 5)

2 FORMAT (214)

END

81
F'ROGRAM PERF"

C COMPUTES MIE SCArTERING FOR A PERFECT CONDUCTOR

PARAMETER <Ny~a 192, NTM=3000, NT4==1200)

CHARACTER XLAB«-:30, YLAB^<30. PTLABUtSO, PTLAB2-K-80. SBLAB*0O

REAL PT I( NTM ) , T I < NTM ) , AMP 1 < NTM ) , AMP2 ( NTM ) . WW ( NW )

,

1

SI < NW ) , S2 < NW ) , S3 ( NW ) , S4 ( NUI

)

DIMENSION PI (5, NT4 ) , P(NT4), TAU(D, NT4 ) , X(N14)

COMPLEX FT 1 < MW ) , FT2 < NW ) , FAC , EYE

CGMMON/PERAR/ ANGL E ( 5 ) , AG ,1 R ( 5, NW ; , AG 1 U 5, NW )

,

1 AG2R<5, NW); AG21 (3, MW)

.

2 BES ( NT 4 ) ; EWEU ( N 1 4 ) .. R A ( N T 4 ) , £ I A ( NT4 ) , RB ( NT4 ) , E I B ( NT4 )

.. 3 P r [ , n ; AMP 1 , AMP2, WW; FT 1 .. FT2; S 1 , S2; B3> S4 , P I , P , TAU. X

COMMON/UAR/ IMANGL, RHO, EM

DATA XLAB/ 'C-kTIME; (METERS) V> YLAB/ ' IMTENSl TY V

DATA PTLABl/ 'SCATTERED INTENSITY, HORIZONTAL POLARIZATION'/

DATA PTLAB2/ 'SCATTERED INTENSITY, VERTICAL POLARIZATION '/

DATA C/3. Ee/,EYE/(0. , 1. )/

CMl^l. /C

READ H-, I PLOT

OPEN ( 7, FILE^ 'PERFSV '. FORM=- 'UNFORMATTED ' ) R.EAD(^«; 1 ) CTMAX, CTMIN, Rl, ALPHA, BETA

READ ( *, 2 ) NT, NSTAR, NANGL

C,DT=(CTMAX-CTMIN) / (NT-1 ) IF (NANGL. GT. 5) STOP 1 1F(MSTAR. GT. MW) STOP 8

IF (NT. GT. NTM) STOP 3

READ *, (ANGLE( I ) , 1^1, NANGL)

READ(-a , 1) DW, WMAX

NRE(.2=1

CTIM=CTMIN

DO 90 1 = 1, NT

PTI ( I )=CTIM

TI ( I )=PTI ( I ) kCMI

CTIM=CTIM+CDT

90 CONTINUE

CALL ONEGA ( WW, WMAX , NGTAR , DW

)

B2A2C2= ( DE;:TA k «-2-ALPHA*«2) -KC* fr2 BAC = - ( BETA-ALPHA ) ^^C-«-2. ABC = - ( AL PHA+BETA )

-« C

BC2=BETA-s^C-Jt*2

ABC2=ALPHA«BC2

MMAX=WMAX*Rl/C+7.

IF(MMAX. GT. NT4-2) THEN

PRINT 'NMAX ',MMAX

STOP 20

END IF

I x -=o TL=SECOND(

)

DO 100 IW=2, NSTAR

RH0=R1 )^WW( IW) / C

EM=RH0+7.

M=EM+. 5

IF (M. GT. NT4-2) STOP 21

IX=IX+1

tl2
CALL BESRJ(RHO, . 5, M+2, 1, BES, NCALC)

IF(MCALC. NE. Mh-2) THEN

PRINT 'ERROR IN COMPUTATION OF BESSEL FUNCTIONS ' STOP 10

END IF

CALL BESRY(RHO, . 5, l'H-2; 1, ENEU. NCALC)

IF(NCALC. NE. h^2) THEN

PRINT K-, 'ERROR IN COMPUTATION OF NEUMANN FUNCTIONS ' STOP 1

1

END IF

C COMPUTATION OF SCATTERING COEFFICIENTS

CALL SCACO

C COMPUTATION OF MONOCHROMATIC FIELDS

CALL FLDdW. M)

100 CONTINUE

TL = SECOND < )-TL

PRINT 'TIME TO COMPUTE MONOCHROMAT I C FIELDS IS ',TL, ' DO 210 IAMG=1,NANGL

FTl ( 1 )=--0. FT2( 1 )=-0,

DO 200 I-- 2.. NSTAR

W=WN ( I

)

WRC=W-frRi/C

FAC=CMPLX(B2A2C2.. EAC-i^W ) /CMPLX ( ABC2- W t^^-2, ABC*W;*^^2

FAC=FAC-«-CMPLX (COS(WRC ) . SIN(WRC) ) FAC = BC2-»-FAC/Ui

FTl ( I )=EYE*CMPLX ( AG1R< lAlMG, I ) , AGl I < lANG, I ) )rFAC

FT2( I )=^EVE^^-CMPLX < AG2R ( lANG, I ). AG2I < lANG, I ) )*FAC

200 CONTINUE

TL^SECGND(

)

CALL VINUFKTI, AMP 1 , NT, NSTAR, WW, 1, FTl, SI, S2, S3, S4, -1 ) TL=SECOND( )-TL

PRINT K-, 'TIME TO COMPUTE A FOURIER TRANSFORM IS ',TL, ' WRITE (7) AMPl

C COMPUTE OTHER POLARIZATION UNLESS THETA =0, ISO

IF(AMGLE(IAMG). NF. 0. . AND. ANGLEUANG). ME. 180. ) THEN

TL=-S£COMf>( ) CALL VINUFTCTI, AMP2, NT, NSTAR, WW, i, FT2.- SI, 52, S3, S4, -1

TL=-S£COWD( )-TL

PRINT TIME 10 COMPUTE A FOURIER TRANSFORM IS '

, TL

WRITE (7) AMP

2

END IF

210 CONTINUE

REWIND 7

DO 300 IANG-1,NAMGL

WRITE(SBLAB, 3) ALPHA, BETA, Rl> ANGLE! lANG)

IF(ANGLE( lANG). NE. 0. . AND. ANGLE! lANG). NE ISO. ) THEN

READ! 7) AMP

2

READ! 7) AMPl

ELSE

READ! 7) AMPl

END IF

DO 220 I^-1,NT

AMPl ( I )=AMP1 ! I )^<•fr2

, NCALC

, NCALC

SECONDS

SECONDS •'

)

. 'SECONDS '

83
220 CONTINUE

IF< IPLOT. EQ. 1. OR. IPLOT. EQ. 2) THEN

CALL DRAW4< 1, 2, 16, 12, 44, 72, XLAB, YLAB, PTLABl, GBLAB ) CALL DRAW4A( 1, 1, NT, ' ', 0, 0, PTI, AMPl, 0. . 0. ) CALL DRAW4B< 1, 0, 0, 0, NT, PTI, AliPl, 2. , 2. ) END IF

IF (IPLOT. EG. 0. OR. IPLOT. EQ. 2) THEN

CALL DRAW3( 1, 2, 16, 12, 44, 72, XLAB, YLAB, PTLABl, SBLAB)

CALL DRAW3A( 1, 1, NT, ' ^ 0, 0, PTI , AMP 1 . 0. , 0. ) CALL DRAW3B(1, 0, 0, 0, NT, PTI, AMPl, 2. , 2. ) END IF

IF < ANGLE (lANG) . NE. 0. . AND, ANGLE(IANG). NE. ISO. ) THEN

DO 260 1=1, NT

AMP 2 ( I ) =AMP2 ( I ) •fi-s-2 260 CONTINUE

IF ( IPLOT. EG. 1. OR. IPLOT EG. 2) THEN

CALL DRAW4( 1, 2, 16, 12, 44, 72, XLAB, YLAB, PTLAB2, SBLAB)

CALL DRAW4A(1, 1, NT, ' ', 0, 0, PTI, AMP2, 0. , 0. ) CALL DRAN4D< 1, 0, 0, 0. NT, PTI, AMP2, 2 ,2. ) END IF

IF( IPLOT. EQ. 0. OR IPLOT. EQ. 2) THEN

CALL DRAW3( 1, 2, 16, 12, 44, 72, XLAB, YLAB, P7LAB2, SBLAB)

CALL DRAW3A(1, 1,NT, ' ', 0, O, PT

I

, AMP2, 0. ,0. ) CALL DRAW3B( 1, 0, 0, 0, NT, PTI, AMP2, 2. , 2. ) END IF

END IF

300 CONTINUE

CLOSE (7)

CALL G5PURG£< O..FN= O^ERFSO', ' STATUS= ', I STAT ) IFdSTAT. GT. 0. AND. ISTAT NE. 1402) STOP 4

CALL Q5DEFIME( 'LFN= 'PERFSU ')

STOP

1 FORMAT (E 12. 5)

2 FORMAT (14)

3 FORMAT ( 'ALPHA = ',F0. 2.. BETA = ',F 5. 2,

1 RADIUS U F4. 2, ', ANGLE ',r4.0)

END

SUBROUTINE SCACO

PARAMETER i NW 1 '"'2, T!T .QO-.-w, m. 4= 1200)

p ' I ‘ r.;! ; . , T 1 ( NT M > , AMP 1 ( I'n ^ ^ AM: '2 ( im i rl ; , NW ( NO ^

j. b 1 ( NW ) , S2 ( )

• ••"'3 ^ N', .: ) , S4 I'lN n T M'"t'.!C .[ OtM F * ' NT4 ^ , P • i -i i 4 ) , I AU ( 5, NT4 ) , X ( NT 4 ) COMPLEX F 'l 1 (NN; , FT2(NW) , FAC, EYE

COMMON/PERAR /ANGLE ( 5 ) , AGIR ( 5, NW )

,

AGl I ( 5, NW )

,

1 AG2R ( 5, NW) , AG21 ( 5, NW)

,

2 BES(NT4), ENEU(NT4), RA(NT4), EIA(NT4), RD(NT4). EIB(NT4),

3 PTI, TI, AMPl, AMP2, WW, FTl, FT2, SI , S2, S3, S4, P I , P. TAU, X COMMON/ VAR/ NANGL, RHO, EM

IERR=0

M=EM

C CALCULATION FOR PERFECTLY CONDUCTING PPUP-nr r~ DO 82 I = 1 , M

AY:= I
ATER=-(AY+.t. )*e£S( I + l )-RHO*BES( 1+2)

BTER=(AY+1. )*ENEU< I + l )-RHO*ENEU< 1+2)

A2=ATER**2

FRAC=1. /(A2+BTER**2)

RB< I )=-A2*FRAC

EIB ( I )=ATER*BTER*FRAC

ATER-BES( I+l

)

BTER=ENEU< I+l

)

A2=ATER*#2

FRAC=1. /(A2+BTER**2)

RA(I )=-A2*FRAC

EIA( I )=ATER>BTER*FRAC

82 CONTINUE

RETURN

END

C

SUBROUTINE FLDdW, LMAX)

PARAMETER (NW=8192, NTM=3000, NT4=1200)

REAL PTI (NTM), TI (NTM), AMPl (NTM), AMP2(NTM), WW(NW),

1

SI (NW), S2<NW), S3(NW), S4(NW)

DIMENSION PI (5; NT4), P<NT4), TAU(5, NT4), X<NT4)

COMPLEX FTl (NW), FT2(NW), FAC, EYE

COMMON/PERAR/ANGLE(5), AG1R<5, NW), AGII (5, NW),

1 AG2R(5, NW), AG2I (5, NW),

2 BES(NT4), ENEU(NT4), RA(NT4) , EIA(NT4) , RB(NT4), EIB(NT4),

3 PTI, TI, AMPl, AMP2, WW, FTl, FT2, SI, S2, S3, S4, PI, P, TAU, X

COMMON/VAR/ NANGL, RHO, EM

DO 3 J=l, NANGL

THETA=ANGLE( J)* 1. 745329251994E-2

P(l)=COS (THETA)

PI ( J, 1 ) = 1. TAU( J, 1 )=P( 1

)

P(2) = l. 5*P( 1 )**2-. 5

PI < J, 2)=3. >P< 1

)

TAU( J, 2)=4. *P(2)-1.

DO 2 L=3, LMAX

EL=L

P(L> = (2. -«-EL-l. )/EL-«-P( 1 )*P(L-1 )-(EL-l. ) /EL*P ( L-2 ) PI ( J, L)=P( 1 )*PI ( J, L-1 )+EL*P(L-l

)

2 TAU( J, L)=EL*(EL+1. ) *P < L ) -P < 1 ) *P I < J, L

)

3 CONTINUE

DO 4 L=l, LMAX

EL=L

4 X (L) = (2. *EL+1. )/(EL*(EL+l. ))

M=EM+. 0001

DO 6 J= 1 , NANGL

AGIR ( J, IW)=0.

AGII ( J, IW)=0.

AG2R( J, IW)=0.

AG2I ( J, IW)=0.

DO 6 L=l,

M

AG1R< J, IW)=AG1R( J, IW)+X(L)*(RA<L)«-PI < J, L)+RB(L) JtTAU< J, L) ) AGIKJ, IW)=AG1I(J, IW)+X(L)*(EIA(L)*PI ( J, L)+EIB(L)«TAU( J, L) ) AG2R(J, IW)=AG2R(J, I W ) +X ( L ) * ( RB ( L ) *P I ( J, L ) +RA ( L ) *TAU ( J, L )

)

85
'w.'

AG2KJ, IW)=AG2I( J, IW)+X <L)-k-<EIB<L)*PI < J, L)+EIA<L)*TAU( J, L) ) 6 CONTINUE t .

' )?';

I RETURN ; »

;

>••<;•> '*7 i END "' :-': V.. -0 ..: nfit. f»A-«( r )flfl

' ^’5.'' H'-vV fj#r c‘i 9* •• T ,. . -

,

j I

^ ':,>TJ. ,A«f*uo

I/jiifcT :''V-"'

’ .; , 5AP'^*>£A-»( I )AW

m%7::S v: ..fa '.'aUUfTMOS

r.

.8

»-iRUr3U^

0l^^

'6! DU (r09^»fTM.aQAC«^^T|1.v5j?lf^«WK>f3?3HARA^ M‘ i ^5 HTn>4fWA <<;’ f Hi I sv* ( HTH ) r TS JA3R •

'

<' WM) ^e . i wM*|a , (^Mi)SB , (uin > 1 8 I'

.

^Tii^.*uAT .

.; f tm;.:<^ v4 ATH .e > iR ,«oi8i/i.’Mia'

, /'• .#Ci. ^<a; -' ' 'K.MMO;)

''•"/

-f ’ t,C' .<• NT, i*’i i

‘ 1 3 t wk) »u I XSOA i

-V < , (^TM>9« tfATVl)AI3 .tATVDAR . t ^ V». "J>13 . ( ^

X UAT . 1 aTJ i*^3 X8 vS;e . 18 *st3 , it^. .w« .s:~-ma . 'ma . . ; . : t a c

. m3 , OMR . JO,' ^ NRAVM^ommoo

* fi'-i'

3jjf

I

• . I “C C 00

, a-af*5^*cj:v«€C , i *

,

< u .> a jom/v*at3ht

91 . 'K. ;*i'-:i ''•^ 5 (m3ht> eo>.<«>2

>

:; : A i,^Kr-4 <# V J -S'-' i:* < I .U ) I

A

. - ^Wu' ,

,

.' ,^ ,. ...

,

•p ';jo i

'

'

v'*\;;r‘NSS

, Ji'' :'u; '^ ' '/rF^

( I >R»< I ilr)UAT

e

^r;.:' V- Vtf. ,

- jj.' J'

V '\: ,

• ...'R.Jl

.

• '-*- (5 .1. t ,.x-f s

) •-:% J, 'UAT

,. , 'rslj,.v.^ ' 00

1^ -i3 , <<;•„ 'Jv„ t .13) ' 1 -J X )R« J=i ' i- j1* ';' '-J)S

i: ! •.). .‘iq

: J U >9W J)S ' . ' ' • '-I.U.'I'C'' &

J H‘ um!TV!Od r

V,’ ; . .

,

I '1 <.-•., "J A 90

K . . NM '

.

fiiirrm -**"^3

f.<i' ^ ;

: /• < ^ J3 > « Jl> \ < J * ..’.il '• r. ) • IJ 1 X >

'

. s i * s Lv fa»M > - A,-; 1 n ‘

. rts-m

iiPI'-'S' »• “'-' <w7(.s'^' •

.

, .#’2" •'• *,S

^ ,^^^,.| X,) UAT.’f^f >ti

}

a>H 'J • ‘V‘

’ ’ - •r.i.fU 3> ' - J

• * : 3 PC<

.•'•'

.1 J .1,. >'UA i -r {

J

) AR+ < I ' ' =*'

' '*“ «

0*iWl .DRS^BA

0(MI .UMfiOA

H.i«j it> on

Mm .u M^<M1

C'A-^ n.r.LT

86

ca

...OS.

'

V..'> V ,.: f . i

'

' ii-i w L'. . ^ .'i' K
PROGRAM PERFPLT (INPUT, OUTPUT)

C PLOTS PREVIOUSLY CALCULATED SCATTERED FIELDS

CHARACTER XLAD*30, YLAB«30, PTLAB1*80, PTLAB2*80, SBLAB^^SO

DIMENSION PTI (3000), AMPOl (3000), AMP02(3000),

1 T1 (300), T2(300), T3(300), T4(300), AMPl (300), AMP2(300),

2 AMP3 ( 300 ) , AMP4 ( 300 ) , ANGLE ( 5

)

LOGICAL Two, NEWT 1, NEWT2

CHARACTER ^<-4 (5L 1 , QLP.. QL3, (5L4, SL ( 4 ) CHARACTER^s-40 L(6)

DATA XLAD/ 'C«TIME (MEIERS) '/, YLAB/ ' INTENSITY '/

DATA P7LAB1/ 'SCATTERED INTENSITY, HORIZONTAL POLARIZATION'/

DATA PTLAD2/ 'SCATTERED INTENSITY, VERTICAL POLARIZATION '/

DATA C/3. E8/, NEWTl, NEWT2/2*. FALSE. /

READ «, GLl, QL2, QL3, QL4, CTl, CT2, IPLOT, INVER. (L( I ), 1 = 1, 6)

PRINT •«, 'COMPARISON OP PLOTS PRODUCED BY PERF WITH THOSE',

i

' PRODUCED BY SCAT, IF ANY'

PRINT 'QUALIFIERS '

, QL 1 , QL2, QL3, QL4

DO 100 1 = 1,6

PRINT 4

PRINT *, L(I)

100 CONTINUE

PRINT 4

IF (CTl. NE. 0. AND. CT2. NE. O. . AND. CTl. GT. CT2) THEN

PRINT 'INITIAL TIME ',CT1. ' GREATER THAN FINAL TIME ',CT2

STOP 10

END IF

IF (CTl. NE. 0. ) NEWTl-=, TRUE.

IF(CT2. NE. 0. ) NEWT2=. TRUE.

CALL Q5ATTACH( 'LFN=', 'PERFSV')

0PEN(7, FIlE= 'PERFSV', FORM= 'UNFORMATTED ' , STATUG='OLD' ) READ 1, CTMAX, CTMIN. Rl, ALPHA, BETA

READ 2, NT, NSTAR.. NANGL

C D T= ( C TMA X -C TM I W ) / ( NT - 1

)

IF (NANGL. GT'. 5) STOP 1 IF (NT. GT. 3000) STOP 3

READ ( ANGLE ( I ) , 1=1, NANGL

)

READ 1, DW, WMAX

NREC=1

NQ=0

iF(QLl. NE, 'MO' ) THEN

NQ=NQ4-1

SL(NQ)=QL1

CALL Q5ATTACH( •'LFN= ', 'SCFL'//QLD

C)PEN(e, FILE= 'SCFL '//QLl, STATUS= ' OLD '

,

1 ACCESS= 'DIRECT', FORM= 'UNFORMATTED ' , RECL=2400)

END IF

IF( QL2. NE. 'NO ' ) THEN

NQ=N(3M

SL ( NQ ) =QL2

CALL Q5ATTACH( 'LFN= ', 'SCFL'//QL2)

OPEN (9, FILE= 'SCFL '//QL2. STATUS= ' OLD '

,

1 ACCESS= 'DIRECT ', FORM= ' UNFORMATTED ' , RECL=2400)

END IF

IF(QL3. ME. 'MO ' ) THEN

87
I NQ“NG+ 1 SL ( NQ ) =QL3

CALL Q5ATTACH< 'i 'SCrL'//QL3)

OPEN< 10, FILE='SCFL V/QL3, GTATUG='LD^

1 ACCESS^ 'DIRECT', F0RM= 'UNFORMATTED '

, RECL=2400)

END IF

IF(QL4. NE. 'NO' > THEN

NO-NGtl

SL ( NQ ) =QL4

CALL Q5ATTACH( 'LFN= 'SCFL'//QL4)

OPEN( 11, FILE= 'GCFL '//QL4, STATUG= 'OLD

1 ACCESS^-- 'DIRECT FGRM== 'UNFORMATTED '

, RECL=2400)

END IF

DD 500 IANG=1,NANGL

WRrTE(SBLAB, 3) ALPHA, BETA, Rl, ANGLE! lANG), ( GL < I

)

, 1 = 1. NQ)

IF(ANGLE( lANG). NE. 0. . AND. ANGLE(IANG). NE. 180. ) THEN

TWO=. TRUE.

IF( INVER EQ. 1 ) THEN

READ (7, END=600 ) AMP02

READ(7) AMPOl

ELSE

READ ( 7, END=600 ) AMPOl

READ(7) AMP 02

END IF

EL.SE

TNO= F-ALBE.

READ ( 7 , END=600 ) AMPO

1

ENDIF

CTIM=CTMIN

DC 210 1=1 , NT

AMPO id) =AMP0 1 d ) -K- 5^2 PTX ( I )=CTIM

CTI M=CTIM-iCDT

210 CONTINUE

I F ( NEWT 1 ) THEN

DC) 222 101 = 1. NT

I

F

i P T 1 d 0 1 ) GE . C T 1 ) GO TO 223

222 CONTINUE

PRINT 'INITIAL TIME OCTl, ' TOO LARGE'

S T OP 1 0

ELBE

101 = 1 END IF

223 IF(NEWT2) THEN

DO 224 I02~-NT.. 1, "1

I F ( PT

I

< 1 02 ) . LE. CT2 ) GO TO 225

224 CONTINUE

PRINT *, 'FINAL TIME ',CT2, ' TOO GMALL

STOP 10

ELBE

I02=NT

END IF

225 IF<I01.GT. 1) THEN

DO 226 13=101, 102

88
AMPOl ( 13-101 -t-1 >=AMP01 ( 13)

PTI ( I3-I01+1)=PTI ( 13)

226 CONTINUE

END IF

NTP=I02-I01+

1

IF<QL1.NE. 'MG') THEN

READ (8, END-600) T1

READ(8) AMPl

DO 240 NT 1=1..300

IF< ANPl (MTl ). LT. 0. ) GO TO 245

240 CONTINUE

245 NT1=NT1-1

IF(NEWTl) THEN

DO 252 ,T1 = 1.NT1

IF(T1 ( 1 1 ). GE. CTl ) GO TO 253

252 CONTINUE

ELBE

11 = 1 END IF' 253 IF(NEWT2) THEN

DO 254 I2=1MT1,. 1 . -1

IF(T1 ( 12) . LE. CT2) GO TO 255

254 CONTINUE

ELSE

I2=NT1

END IF

255 IFdl. NE 1 ) THEN

DO 256 13=11, 12

AMP 1 ( 1 3-1 1 +-1 ) =AMP 1(13)

T1 < I3-I l-l )=T1 ( 13)

256 CONTINUE

END IF

NT1=12-I 1+1

END IF

IF(QL2. NE. 'NO ' ) THEN

READ(9, END=600) T2

READ(9) AMP

2

DO 260 NT2= 1.-300

IF( AMP2(NT2). LT. 0. ) GO TO 265

260 CONTINUE

265 NT2=NT2-1

IF(NEWTl) THEN

DO 272 Ii=l,NT2

IF(T2d 1 ) GE CTl ) GO TO 273

272 CONTINUE

ELBE

11 = 1 END IF

273 IF(MEWT2) THEN

DO 274 I2=NT2> 1,-1

IF(T2(I2). LE. CT2) GO TO 275

274 CONTINUE

ELSE

I2=NT2

89
E£ND IF

275 TF< II. NE. 1 ) THEN

DO 276 1

3---= I I, 12

AMP 2 ( 1 3- 1 :l. 4 1 .) ::*-AMP2 0. 3 .> T2( I3“1 1 + 1 ):^T2( 13)

276 CONTINUE

END IF

NT2=I2-I 1+1

END IF

IF(QL3. ME. 'NOM THEN

READ< 10. END==600) T3

READ (10) AMP

3

DO 200 NT3:=-1,300

IF( AMP3(NT3) . LT. 0. ) GO TO 285

230 CONTINUE

285 NT3=NT3~1

IF'INEWTl) THEN

DO 292 I 1=^1, NT

3

IF(T3< II ), GE. CTl ) GO TO 293

292 CONTINUE

ELSE

11 = 1 END IF

293 IF(NEWT2) THEN

DO 294 I2=NT3, 1,-1

IF(T3( 12). LE. CT2) GO TO 295

294 CONTINUE

ELBE

I2=NT3

END IF

295 IF( I 1 . NE. 1 ) THEN

DO 296 13=11, 12

AMP3( I3-I 1+1 )=AMP3( 13)

T3( I3-1 1 + 1 )-=T3( 13)

296 CONTINUE

END IF

NT3=I2-I 1 1 END IF' IF(0L4. NE 'NO' ) 1 HEM

R EAD ( 1 1 , END=600 ) T 4

READ<il) AMP4

DO 300 NT 4 =1,300

IF(AMP4(NT4) . LT. 0. ) GO TO 305

300 CONTINUE

305 NT4=NT4~1

IF(NEWTl) THEN

DO 312 I1=1,NT4

IF(T4(1 1 ) . GE. CTl ) GO TO 313

312 CONTINUE

ELBE

I 1 = 1 END IF

313 IF(NEWT2) THEN

DO 314 I2=NT4, 1, ~1

90
314

315

316

320

326

340

34 5

IF(T4( 12). LE. CT2) GO TO 315

CONTINUE

ELSE

I2=NT4

END IF

IFdl.NE. 1) THEN

DO 316 I3=-1 1 , 12

AMP4( 13- 1 1 + 1 )=AMP4( 13)

T4 (1 3-1 1 + 1 )=T4( 13)

CONTINUE

END IF

NT4--=I2-1 1 + 1 END IF

NREO-NREC+2

NMl^^NT 1/20+4

NM2-NM1 +1

NM3^NN2+1

NM4=NM3+1

IFdFLOT. EQ. 1. OR I PLOT. EQ. 2) THEN

CALL DRAW4( 1, 1 , 16, 12, 44, 80, XLAB, YLAD, PTLABl, SBLAB

)

CALL DRAW4Ad, 1.. NTP, ' ' , 0, 0, PTI, AMPOl.- 0. , 0. ) IF(QL1.NE. 'NO') CALL DRAW4 A < 1 , 1 , NTl , ' 0 ' , NMl , 1 , T 1 , AMP 1 , 0. , O. ) IF(QL2. NE. 3M0' ) CALL DRAW4Ad, 1, NT2, ' 1 ', NM2, 2, T2, AMP2, O. , 0. ) IF(QL3. NE. 'NO') CALL DR AW4A d , 1 , NT3, ' 2 '

, NM3, 3, T3, AMP3.. 0. , 0. ) IF(QL4. NE. 'NO') CALL DR AW4 A d , 1 , NT4 . ' 3 '

, NM4, 4, T4, AMP4, 0. , 0. ) CALL DRAW4Bd, O, 0, 0, 3000, PII, AMPOl, 2. , 2. ) END IF

IF( IPLOT. EQ. 0. OR. IPLOT. EQ. 2) THEN

CALL DRAW3( 1 . 1, 16, 12, 44, 80, XLAB, YLAB, PTLABl, SBLAB)

CALL DRAW3Ad, 1, NTP, ' ', 0, 0, PTI , AMPOl , 0. , 0. ) IF(QL1 ME. 'NO' ) CALL DR AW3A d , 1 , NT 1 ' 0 '

, NMl , 1 , T 1 , AMP 1 , 0. , 0. ) IF(QL2 NE. 'NO') CALL DR AW3A d , 1 , NT2, ' 1

'

, MM2, 2, T2, AMP2, 0. , 0. ) IF(QL3. NE. 'NO') CALL DRAW3A d , 1 , NT3, ' 2 '

, NM3, 3, T3, AMP3, 0. , 0. ) IF(QL4. NE. 'NO') CALL DRAW3A d , 1 , NT4, ' 3 '

, NM4, 4.. T4, AMP4, O. , 0. ) CALL DRAW3B d, 0, 0, 0, 3000, PTI, AMPOl, 2. , 2. ) END IF

IF < TWO) THEM

DO 320 I==1,MT

AMP 02 < I )-AMP02( I )k-k-2

CONTINUE

IFdOl.NE. 1) THEN

DO 326 I

'3--^^ 1 01.. 102

AMP02( 1 3-101 + 1 )=^AMP02( 13)

CONTINUE

END IF

IF(QL1. NE 'NOd THEN

READ (8, END-600) T1

READ(8) AMPl

DO 340 NT 1 = 1, 300

IF( AMPKNTl ) . LT. 0. ) GO TO 345

CONTINUE

|\IT1=NT1-1

IFHMEWTl) IHEN

DO 352 11=1, NTl

91
352

353

354

355

356

360

365

372

373

374

375

376

IF(T1 < 1 1 ). GE. CTl ) GO TO 353

CONTINUE

ELSE

11 = 1 END IF

IF(NEWT2) THEN

DO 354 I2=NT1, 1, -1

IF(T1 ( 12). LE. CT2) GO TO 355

CONTINUE

ELSE

I2=NT1

END IF

IF< II, NE. 1 ) THEN

DO 356 13=11, 12

AMPK I3-I1 + 1 )=AMP1 ( 13)

T1 ( I3-1 14-1 )=T1 ( 13)

CONTINUE

END IF

NTl = I2-Il4-i

END IF

IF(QL2. NE. 'NOM THEN

READ(9, END=600) T2

READ<9) ANP2

DO 360 NT2= 1,300

IF(AMP2(NT2) . LT. 0. ) GO TO 365

CONTINUE

NT2=NT2-1

IF ( NEWT 1) THEN

DO 372 11=1, NT2

IF(T2< II ). GE. CTl ) GO TO 373

CONTINUE

ELSE

11 = 1 END IF

IF(NEWT2) THEN

DO 374 I2=NT2, 1,-1

IF(T2< 12), LE. CT2) GO TO 375

CONTINUE

ELSE

I2=NT2

END IF

IF( II. NE. 1 ) THEN

DO 376 13=11, 12

AMP2< I3-I14-1 )=AMP2( 13)

72(13-11+1 )=T2( 13)

CONTINUE

END IF

NT2=I2-I 1+1

END IF

IF(QL3. NE. 'NOM THEN

READ (10, END=600 ) T3

READ(IO) AMP3

DO 380 NT3= 1,300

IF(AMP3(NT3). LT. 0. ) GO TO 385

92
380

385

392

393

394

395

396

400

405

412

413

414

415

416

CONTINUE

NT3=NT3-1

IF ( NEWT 1) THEN

DO 392 11=1, NT3

IF(T3( II ). GE. CTl ) GO TO 393

CONTINUE

ELSE

11 = 1 END IF

IF(NEWT2) THEN

DO 394 I2=NT3, 1, -1

IF(T3( 12). LE. CT2) GO TO 395

CONTINUE

ELSE

I2=NT3

END IF

IFdl.ME. 1) THEN

DO 396 13=11, 12

AMP3( I3-I1+1 )=AMP3< 13)

T3( I3--1 1 + 1 )=T3( 13)

CONTINUE

END IF

NT3=I2-I1+1

END IF

IF(QL4. NE. 'HO'

)

THEN

READ( 1 1, END=600) T4

READ(ll) AMP4

DO 400 NT4= 1,300

IF(AMP4(NT4). LT. 0. ) GO TO 405

CONTINUE

NT4=NT4~1

IF ( NEWT 1) THEN

DO 412 11=1, NT4

IF(T4< II ). GE. CTl ) GO TO 413

CONTINUE

ELSE

11 = 1 END IF

IF(NEWT2) THEN

DO 414 I2=NT4, 1, -1

IF(T4< 12). LE. CT2) GO TO 415

CONTINUE

ELSE

I2=NT4

END IF

IF( II. NE. 1 ) THEN

DO 416 13=11, 12

AMP4< I3-I1+1 )=AMP4( 13)

T4( I3-I1+1 >=T4( 13)

CONTINUE

END IF

NT4=I2-I 1+1

END IF

NREC=NREC+2

93
CJ

PJ

NMl«NTl/20+-4

NM2«NM1+1

NM3aNM2+l

NM4«NM3+1

IF< IPLOT. EQ. 1. OR. IPLOT. EQ. 2) THEN

CALL DRAW4( 1, 1, 16, 12, 44, 80, XLAB, YLAB, PTLAB2, SBLAB)

CALL DRAW4A( 1, 1, NTP, ^ 0, 0, PTI, AMP02. 0. , 0. ) IF(QL1. NE. 'NO' ) CALL DRAW4A< 1, 1, NTl, ' 0

'

, NMl , 1 , T1 , AMP 1 , 0. , 0. ) IF(QL2. NE. 'NO' ) CALL DRAW4A( 1, 1, NT2, ' 1

'

, NM2, 2, T2, AMP2, 0. , 0. ) IF<QL3. NE 'NO') CALL DRAW4A < 1 , 1 , NT3, ' 2

'

, NM3, 3, T3, AMP3, 0. , 0. ) IF(QL4. NE 'NO') CALL DRAW4A (

1

, 1 , NT4, ' 3 ' , NM4, 4, T4, AMP4, 0 , 0. ) CALL DRAW4B (1 , 0, 0, 0, 3000, PTI , AMP02, 2. , 2. ) END IF

IF < IPLOT. EQ. 0 OR. IPLOT. EQ. 2) THEN

CALL DRAW3< 1, 1, 16, 12, 44, 80, XLAB, YLAB, PTLAB2, SBLAB)

CALL DRAW3A (1 , 1 , NTP, ' '

, 0, 0, PTI . AMP02, O. , O. ) IF(QL1. NE. 'NO ' ) CALL DRAW3A( 1, 1, NTl, ' 0 ' . NMl , 1 . T1 , AMP 1 . 0. , 0. ) IF(QL2. NE. 'NO') CALL DRAW3A ( 1 , 1 , NT2, ' 1

'

, NM2, 2, T2. AMP2, 0. , 0. ) IF(QL3. NE. 'NO') CALL DRAW3A ( 1 , 1 , NT3. ' 2 '

, NM3, 3, T3, AMP3, 0 , 0. ) IF(QL4. NE. 'NO ' ) CALL DRAW3A ( 1 , 1 , NT4, ' 3 '

,

NM4, 4, T4, AMP4, 0. , 0 ) CALL DRAW3B( 1, 0, 0, 0, 3000, PTI, AMP02, 2. . 2. ) END IF

END IF

500 CONTINUE

600 STOP

1 FORMAT (E12. 5)

FORMAT ( 14)

FORMAT (

' ALPHA= '

, F5. 2, '

, BETA= '

, F5. 2,

1 ', RADIUS= ', F4. 2, ', ANGLE= ', F4. O:

2 QUAL. ; ',A4;', ',A4;', ',A4:', '

,

A4

)

4 FORMAT <//)

END

94
o

a

a

a

o

o

PROGRAM CURCHK

PARAMETER <LL=LLG)

C LCG, LL12G, AND LL20G NOT NEEDED

DIMENSION CJTH(900, 1 0 ) , C JPH < 900, 10),T<900)

DIMENSION NPCH( 10), CJA(LL)

CHARACTER QAL*4, XLAB*30, YLAB*30, TITLE*80, SUBT*80

COMMON /PASS/ Nl, N2, N3, IS, IPATCH, RAD, C, DTI

DATA XLAB/'C * TIME V, YLAB/ 'SURFACE DENSITY'/,

1 TITLE/ 'SURFACE CURRENT DENSITY'/

DATA C/3. E8/

READ INFORMATION FOR CURRENT CHECK

READ *, QAL, IPLDT, IPRINT, NN

IF<NN. GT. 10) STOP 1 READ »•, (NPCHd ) , 1 = 1, NN)

READ *, CTl, CT2, IPATCH

READ IN VALUES OF CONSTANTS FOR THE PULSE

READ 4, DT, RAD, ALPHA, BETA

READ 5, NI

READ 5, N1,N2, N3

CALL RUNT2

IF<DT. EQ. 0. ) DT=DT1

IF(NI.EQ. 0) NI=6. *RAD/ (C-t^-DT)

CDT=C*DT

NI l=CTl/CDT-f-l,

NI2=CT2/CDT+1.

CALL Q5ATTACH( 'LFN=', 'SCSV'//QAD

0PEN(2, FILE= 'SCSV '//QAL, STATUS= 'OLD ', ACCESS= 'DIRECT ',

1 FDRM= 'UNFORMATTED RECL=ISi{-S)

IF(NI2. EQ. 1 ) NI2=NI

IF(NI1. GE. N12. OR. NI2. GT. 900) THEN

PRINT «, 'NI1=', Nil, ', NI2= ', NI2, CT1=',CT1,', CT2=',CT2

STOP 2

E.ND IF

NI I=NI2-NI 1+1

NK=NI 1/12+2

NK1=NK+1

TT=(NI 1-1 )»CDT

DO 200 I=NI 1, NI2

READ<2, REC = I-k2-1 ) CJA(1;IS)

DO 1 10 N=l, NN

C JTH ( I , N ) =C JA ( NPCH ( N )

)

no CONTINUE

READ<2, REC=I*2) CJAddS)

DO 120 N= 1 , NN

C JPH d , N ) =C JA ( NPCH ( N )

)

120 CONTINUE

T( I )=TT

TT=TT+CDT

200 CONTINUE

IF( IPRINT. NE. 0) THEN

95
o

o

PRINT ^ INPUT FILE SCINSQAL

PRINT RAD=',RAD/'/ AL.PHA= ALPHA, '

• BETA=',BETA,

1 ^ MI=MMI, S N1=',N1,S N^=^N2, M3~',N3

DO 220 N= i , NN

NU=NPCH<N)

DO 210 I=NI1, NI2

PRINT 2, NU, I, CJTH( I, N), CJPH( I, N)

210 CONTINUE

PRINT 3

220 CONTINUE

END IF

IF< IPRINT, NE. 1 ) THEN

DO 250 N=l, NN

WRITE<SUBT, 1 ) QAL, NPCH(N)

IFdPLOT. NE. 1 ) THEN

CALL DRAW3( 1, 7, 8, 15, 23, 33, XLAB, YLAB, TITLE, SUBT)

CALL DRAW3A< 1, 1, NI I, ' 1

'

, NK, 0, T ( NI 1 ) , CJTH(NI 1 , N), 0. , 0 ) CALL DRAW3A(1, 1,NII, ' 2 '

, NKl , 0, T ( NI 1 ) , C JPH< NI 1 , N ) , 0 ,0. ) CALL DRAW3B( 1, 0, 0, 0, NI, CJA, T, 2. , 2. ) END IF

IFdPLOT. NE. 0) THEN

CALL DRAW4( 1, 7, 8, 15, 23, 33, XLAB, YLAB, TITLE, BUBT)

CALL DRAW4A( 1, 1, Nil, ' 1

'

, NK, 0, T < NI 1 ) , CJTH(NH, N), 0. , 0. ) CALL DRAW4A(1, 1, Nil, ' 2

, NKl , 0, T < NI

1

) , CJPH(NI1, N), 0. , 0 ) CALL DRAW4B( 1, 0, 0, 0, NI, CJA, T, 2. , 2. > END IF

250 CONTINUE

END IF

STOP

1 FORMAT ( 'QUALIFIER ',A4, ', PATCH NUMBER ',14)

2 FORMAT( IX, 218, 2(3X, 1PE16. 9)

)

3 FORMAT (///)

4 FORMAT (E12. 5)

5 F0RMAT(3I4)

END

SUBROUTINE RUNT2

C

C THIS SUBROUTINE DETERMINES THE NUMBER OF PATCHES AND

DT FROM THE PATCH DISTRIBUTION ON THE SPHERE

COMMON /PASS/ NI, N2. N3, IS, IPATCH, RAD, C, DTI

C

PI =4. *ATAN( 1. ) DTH=PI/ (Nl~l ) SDTH2=SIN(DTH*. 5)

TH=PI

N1PN(N1+1 )/2

DPH=1.

IS=2

IF ( IPATCH. EO. 0) THEN

NPHX=1

90 NPHX=NPHXit2

IF(NPHX. LT. N3) GO TO 90

DO 100 1=2, NIH

96
100

105

9

TH=TH-DTH

STH=SIN(TH)

NPH=<N2-N3) i^STH+N3+. 5

IF(NPH. GT. NPHX) NPHX=NPHX*2

IB=IS+NPHX*2

DPH=AMIN1 (DPH, STH*BIN(PI/NPHX )

)

CONTINUE

IF(M0D<N1, 2) . EQ. 1 ) IS=IB~NPHX

ELSE

DO 105 1=2, NIH

TH=TH-DTH

STH=SIN(TH)

NPH=(N2-N3)-*«-STH+N3+. 5

IS=IS+NPH»2

DPH=AMINi <DPH, STH-«SIN < P I /NPH ) ) CONTINUE

IF(M0D(N1, 2). EQ. 1 ) IS=IS-NPH

END IF

DT1 = 1. 6»^RAD/C#AMIN1 (SDTH2, DFl!)

RETURN

END

97
SUBROUTINE VNUFKX, Y, N. NU> 0M> lOM, FT, DX, CX, SX, JFLAG

)

C THIS SUBROUTINE CALCULATES THE FOURIER TRANSFORM OF THE

C FUNCTION GIVEN BY STRAIGHT LINES JOINING THE POINTS GIVEN BY THE

C ARRAYS X,Y. THE NUMBER OF IMPUT POINTS IS N, THE FREQUENCY ARRAY OM IS

C PROVIDED AND HAS DIMENSION NU. DX, CX, SX ARE REAL ARRAYS OF DIMENSION N

C THAT ARE USED FOR SCRATCH. THE FOURIER TRANSFORM IS GIVEN BY THE

C COMPLEX ARRAY FT.

C IF IOM=l, INPUT AND OUTPUT OM ARE CIRCULAR FREQUENCIES

C IF I0M=2, INPUT ARRAY IS CIRCULAR FREQUENCY, OUTPUT IS FREQUENCY

C IF I0M=3, INPUT ARRAY IS FREQUENCY, OUTPUT IS CIRCULAR FREQUENCY

C IF I0M=4.. INPUT AND OUTPUT OM ARE FREQUENCIES

C JFLAG^+l IF THE FOURIER TRANSFORM HAS A FACTOR EXP(-HWT)

C JFLAG=-1 IF THE FOURIER TRANSFORM HAS A FACTOR EXP(-IWT)

DIMENSION X(N), Y(N), DX(N.>, CX(N), SX(N), OM(NU)

COMPLEX FT(NU)

TOP I =8. *ATAN< 1. ) TPIN=1. /TOPI

IF( lOM. LE. 2) GO TO 25

OM (1 ; N ) -OM ( 1 i N ) *TOP

I

25 NP--N-1

DX( .t.; NP)=QSVDELT<X(li N)i DX(1; NP) )

NO=----l IF(OM( 1 )

.

EQ. 0. ) THEN

NO-2

S-QSSSUM < Q8VADJM < Y < 1 ; N ) ; NP ) *DX ( 1 ; NP )

)

FT( 1 )-CMPLX<S, O. ) END IF

DX( 1; NP )=Q8VDELT( Y< 1; N); NP)/DX( 1; NP

)

DO 70 I-NO, NU

W=OM( I

)

WIN-1. /W

CX( 1; N)=VCOS(W*X< liN);CX(liN))

3X(1; N)=--VSIN(W*X( 1; N)i SXd; N) ) S2R==Q8SSUM(Q8VDELT(CX( li N)i NP)*DX< 1; NP ) ) S2I -QSSSUM ( Q8VDELT ( SX (1 ; N ) ; NP ) *DX ( 1 ; NP )

)

S1R=Y( 1 )«CX( 1 )~Y(N)*CX(N)

SI I-Y( 1 )itSX( i ) -Y(N)-»J-SX (N)

FT( I )-CMPLX ( ( -SI I +S2R-«-WIN)*WIN, < S 1 R + S2I *WI N ) -«-W I N ) 70 CONTINUE

IF( lOM. EQ. 1. OR. lOM. EQ. 3) GO TO 85

0M< 1; N)-=OM< l5 N) ftTPlN

85 IF( JFLAG. EQ. 1 ) RETURN

FT ( 1 ; NU ) -VCDNJG < FT < 1 i NU > i FT < 1 ; NU )

)

RETURN

END

C

SUBROUTINE VINUFT( X, Y, N, NU, OM, lOM, FT, SCRR, SCRI, CX, SX, JFLAG)

C THIS SUBROUTINE CALCULATES THE INVERSE FOURIER TRANSFORM OF THE FUNCTION

C GIVEN BY STRAIGHT LINES JOINING THE POINTS GIVEN BY THE COMPLEX ARRAY FT

C AS A FUNCTION OF THE REAL ARRAY OM. THE NUMBER OF INPUT POINTS IS NU

C THE INVERSE F. T. IS ASSUMED TO BE A REAL FUNCTION AND GOES IN THE

C ARRAY Y, AND THE REAL ARRAY X IS GIVEN AND HAS DIMENSION N. C SCRR, SCRI ARE REAL ARRAYS OF DIMENSION NU AND IS USED FOR SCRATCH

C IF IOM-1, INPUT AND OUTPUT OM ARE CIRCULAR FREQUENCIES

90
C IF I0M=^2, INPUT ARRAY IS CIRCULAR FREQUENCY, OUTPUT IS FREQUENCY

C IF I0M=3, INPUT ARRAY IS FREQUENCY, OUTPUT IS CIRCULAR FREQUENCY

C IF I0M=4, INPUT AND OUTPUT ON ARE FREQUENCIES

C vJFLAG=h-1 if THE INVERSE FOURIER TRANSFORM HAS A FACTOR EXP(+IWT>

C JFLAG=^-1 IF THE INVERSE FOURIER TRANSFORM HAS A FACTOR EXP(-IWT)

DIMENSION X<N), Y<N), OM(NU), SCRR(NU), SCRI (NU), CX ( NU ) , SX < NU)

COMPLEX FT(NU)

TOPICS. *ATAN< 1. ) TPIN=1. /TOPI

PIN==2. «-TPIN

IF(vlFLAG. EQ. 1 ) GO TO 15

FT (1 ; NU ) =VCONJG ( FT (1 ; NU ) ; FT ( 1 j NU )

)

15 IF<JDM. LE. 2) GO TO 25

OM ( 1 ; NU ) =OM < 1 ; NU ) *TOP

I

25 MP=NU-1

SCRI ( 1.; NP) = 1. /QaVDELT(OM( li NU); NP ) SCRR ( 1 ; NP ) =Q8VDELT < VREAL < FT < 1 , NU ) ; NU ) ; NP ) *SCR I ( 1 ; NP

)

SCRI ( 1; NP)=Q8VDELT<VAIMAG(FT( 1; NU), NU); NP)*SCRI ( 1 , NP

)

DO 70 1=1,

N

W=X ( I

)

IF(ABS(W). LT. 1. E-19) THEN

Y (I ) =P I N*Q8SSUM < Q8VADJM ( VREAL ( FT < 1 ; NU ) ; NU > ; NP )

1 Q8VDELT<0M( 1; NU); NP)

)

ELSE

WIN=1. /W

CX ( 1; NU)=VCGS( W*OM( 1; NU) ; CX < 1 ; NU)

)

SX ( 1 ; NU ) =VSI N ( W*OM ( 1 ; NU ) ; SX ( 1 ; NU )

)

S2R=QaSSUM(Q8VDELT(CX ( 1; NU); NP)*SCRR(1; NP )

1 Q8VDELT ( SX ( 1 ; NU ) ; NP ) <^SCR I ( 1 ; NP ) ) S1I=REAL(FT< 1 ) )*SX( 1 )+AIMAG(FT(l ) >*CXn

)

1 ~REAL(FT<NU) )-;^SX (NU)-AIMAG(FT(NU) )*CX (NU)

Y( I )=PIN*(-S1 1-«-S2R*WIN)*WIN

END IF

70 CONTINUE

IFdOM. EQ, 1. OR. lOM. EQ. 3) GO TO 85

0M( 1; NU)=OM< 1; NU)*TPIN

35 IF( JFLAG. EQ. 1 ) RETURN

FT ( 1 ; NU ) =VCONJG ( FT ( 1 ; NU ) ; FT ( 1 ; NU )

)

RETURN

END

99
NBS*114A iREv. 2 >eo) U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO.

NBSIR 86-3362

2. Performing Organ. Report No. 3. Publication Date APRIL 1986

4. TITLE AND SUBTITLE

SCAT: A Vector Program to Solve a Transient MFIE

5. AUTHOR(S)

Egon Marx

6. PERFORMING ORGANIZATION (If joint or other than NBS, see instructions) 7. Contract/Grant No. NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE I. Type of Report & Period Covered

WASHINGTON, D.C. 20234

9.

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

10.

SUPPLEMENTARY NOTES

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

bi bliography or literature survey, mention it here) The FORTRAN program SCAT is used to solve the magnetic field integral equation (MFIE)

to determine the fields scattered by a perfectly conducting sphere. The incident

field is a plane-wave pulse, and a stepping-in-time procedure is used to determine

the surface current density induced on the sphere. The program does not take

advantage of the special symmetry of the scatterer because it is Intended to serve

as a verified starting point for more general programs. The output is compared to

that of the program PERF, which computes the same fields via a Fourier transform of

the monochromatic fields obtained from the Mie formulas. The contributions of the

self-patch and neighboring patches to the singular integral are optionally computed

by using their expansions in the linear size of the patches. The self-patch term

is important for the solution of other integral equations that may be of the first

kind. For the MFIE, these corrections are small but not negligible. The program

takes advantage of the vector programming features of the CYBER 205.

12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolon s) computer programs; electromagnetic pulse scattering; magnetic field integral

equation; Mie scattering; perfectly conducting sphere; stepping-in-time procedure;

transient electromagnetic fields; vector programming

13. AVAILABILITY 14. NO. OF

PRINTED PACES

[_^ Uni imited

1 1 For Official Distribution. Do Not Release to NTIS rn Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C.

20402.

104

IS. Price

|j^ Order From National Technical Information Service (NTIS), Springfield, VA. 22161 $16.95