NBS

PUBLICATIONS

NATL INST OF STANO 4 TECH

A11107 310537

NBSIR 84-2835

Free-Space Propagation of Light

Pulses

U S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratories

Center for Manufacturing Engineering

Washington, DC 20234

May 1984

U S. DEPARTMENT OF COMMERCE

ATIONAL BUREAU OF STANDARDS

100

. U56

34-2335

NATIONAL BUREAU OF STANDARDS

LIBRARY C l r C

NBSIR 84-2835

FREE-SPACE PROPAGATION OF LIGHT

PULSES

Egon Marx

U S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratories

Center for Manufacturing Engineering

Washington, DC 20234

May 1984

U.S. DEPARTMENT OF COMMERCE. Malcolm Baldriga. Stcrttmry

NATIONAL BUREAU OF STANDARDS. ErmM Amtator. Dktetor

Abstract

A transient electromagnetic field in free space is

completely specified when the initial values of the electric and

magnetic fields are given. Green's function for the scalar wave

equation can then be used to find the field at later times. A

group of computer programs that implement these equations and

process the output are presented in this report.

Key words: computer programs; conservation of energy and

momentum; el ectromagneti c field propagation;

finite light pulses; graphic displays of energy

density; transient electromagnetic fields.

1
Tabl^e of Contents

1. Introduction 1

2. Propagation of Localized Pulses . 2

3. Computer Programs ...... 6

3. 1 Program BALL 6

3.2 Subroutine POINT 8

3.3 Subroutine ICS ........... 8

3.4 Program BALPLT 10

3.5 Program BALDSK 11

3.6 Program BALCHK ........... 12

4. Sample Outputs 13

References 19

Appendix 20
1 . In troducti on

In order to compute transient electromagnetic fields

scattered by conductors or dielectrics, as required for the Wave

Optics program at the National Bureau of Standards, we have to

determine the incident fields specified by their initial values.

These incident fields are free fields that obey Maxwell's

equation without sources in all of space, disregarding the

scatt erer

The programs described in this report allow us to compute

the free fields at arbitrary points and times from given initial

values of the fields.

A simple example is the plane wave: the fields are constant

A

over planes perpendi cul ar to a direction of propagation n. The

4

initial electric field, is a vector field perpendi cul ar to

n that is a function of a single variable? the free electric

field at other times is then given by

E(x,t> = E

q

(£ - ct) (1)

where c is the speed of light in vacuum and

4 A A -4 £ = x * n, n Eq =0. (2)

The magnetic field of this plane wave is then given by

cB(x,t) =nxE(x,t). (3)

For fields other than plane waves, the incident field can be

4 4

computed by integrations from the initial values of E and B. The

1
solution of the free-field Maxwell equations is shown in section

2, computer programs developed to do the numerical integrations

are described in section 3, and a few sample outputs are shown in

section 4. The listings of the computer programs are given in an

appendix

2. tLEQe^gat i^on gf Localized Pulses

We have shown C 1 3 how the retarded Green function for the

scalar wave equation,

s‘

0> <*,*) =

* !t - I 'V ^' l/C)

. (4)

can be used to express the incident electromagnetic fields in

terms of their initial values given at time t = 0 for all of

space. The symbol x stands for the four-vector with components

•4 x , m = 0, 1,2,3, that is, x = (ct,x), and 5(t) is the Dirac delta

M

function. The resulting integrals are

-* f r 2 -* (0)

e

(

x > =- | y

dV' I (i/c )e

q

<x ' )as^ <x,x' ) /at'

- B

q

(x') X ?'G‘

lJ '(x,x')]

t

, =0 , (5)

Numbers in square brackets indicate the literature references at

the end of this report.
B(x> = - (1/c >| v

dV

f

dV'[l0

(^') ~ , _iy; . ,

v B

r

( x , x )

-+ -+ (O) 1 + b

q

(x / ) as

R

i*,x' )/at']

t

, =0 , <6)

where v' is the gradient operator with respect to the variable

-4x'. We substitute the Green -function (4) into eqs. (5) and (6)

to obtain

-» r ri - -» s

f E (x) = J v

dV' [^E0(x'> —(t - R/c)

4nR

+

“f -4 B ^ ( x ' ) x

u

(S ' ( t - R/c) . 5 (t - R/c ) \ £1

( 4nRC 4nR2 )

R

J (7)

B (x ) ( S ' (t - R/c)

\ 4nRc

+

S (t - R/c)

7

4ttR

R

- B.

t - R/c)

ix >

4ttR

( 8 )

-4 -4 *4 -4 *A* -4 where R = x — x' , R = IRI, and R = R/R. We integrate over R and

are left with integrals over a sphere of radius R = ct centered

”4 at the field point x, known as the information-collecting sphere.

-4 We can use the solid angle from x and express the fields as
ct E(x) = j- 5 (9)

4rr MK

- +

“4 ~T “+ "4 1 2cB0 ” R-V , cB

q

) x r| r

-r If f r -4 -+ -+ -4 -+ -+T cB(x) = 4— dO cB - R-v'cB, - (2Ert - R-y'Ert ) x R _ . . (10) ^nj LO 0 O O JR = ct

The fields and B^ are arbitrary solenoidal fields, that

is, they satisfy the constraint equations

EQ

(x) = 0, ( 11 )

v B

Q

(x) = 0. (12)

A simple way to find a solenoidal field is to take the curl of

another vector field.

We would like the field to be spatially localized, but

otherwise to look somewhat like a plane wave propagating in the

2 -direction, We restrict the initial fields to the region z < 0

and we assume that they decay essentially like a gaussian in the

x — and y-di recti ons. We choose

E„(x) = A expj^ — k (x^ + y

J- )j(yi — xj)<t>'(z), (13)

4
cBq<x) = A exp [- k (x

-^ + y2 )j^(xi + yj)$ y (z>

- 2k [l - k ix

2

+ y

2

) J

4> <z > y.

( 14 )

where $ is a function that vanishes for z £ 0 and has a

continuous first derivative. These expressions can be used to

compute the dyadics that are the gradients of and B^.

The initial values of the Roynting vector and energy density

are proportional to

Eq x cBq = A^exp^ - 2k (x*" + - K (x“ + y^)

j <x i + yj )

2 2 2 ~\ + (x + y )$>' k>. (15)

->'? r 2 2

E^T + c^B^ = A^exp — 2k (x + y ) K

2( * 2 2

+ y

)<*'

+ 4^1 - k (

2 + y2

)j

2

#

2

^>. ( 16)

This energy density is invariant under rotations about the

z—axis, and the Roynting vector has a component in the positive

z—direction and one radial component (perpendicular to the

z—axis). For a light pulse, $ has the form

= t

(

z)sin (kz)

,

$ (z ) (17)
where k is the wave number o-f the light produced by a laser.

3 . Computer Programs

A computer program named BALL allows us to calculate the

fields by performing the integrations in eqs. (9) and (10). This

program calls a subroutine ICS where the unit sphere is covered

by the appropriate number of patches required to carry out the

integrations, and a subroutine POINT that is used to step through

the field points. The output of BALL is stored in a file, that

is then used as input to separate programs that present the

results in graphic form. These programs are BALPLT, which

produces a plot of the energy density as a function of the field

points, and 8ALDSK, which prepares a binary interpolated file

that can be shown on the monitor of our group's image processing

facility or photographed by a special computer-dr i ven camera. An

additional program, BALCHK, can be used to perform a check on the

accuracy of the numerical calculations by means of the

conservation of energy and momentum. These programs are listed in

the appendix.

3. 1 Progr am BALL

We use this program to do the integrations required to find

-4 the fields at a time t and at a point x for a given set of

initial fields. In the program shown in the appendix the pulse

described in eqs. (13) and (14) has a shape given by

<|>(z) = [exp (az) — exp (£z ) Isin (kz ) ® (—z ) , (18)

6
i -f a positive k is specified, or

$ ( z ) = z

J~Cexp(az) — exp (£z ) 38 (—z ) (19)

if k vanishes. The unit step function 6 makes the pulse vanish

for positive z at t = 0. For the pulsed lasers used at present,

the underlying fields oscillate maybe one hundred times for the

duration of the pulse, but this number is being reduced in

ongoing experimental research programs. Since Maxwell's

equations in free space are invariant under changes of scale in

space and time, the same calculations can be applied to other

el ectromagnet i c fields that may be true pulses, without an

underlying sinusoidal field.

The parameters of the pulse, the region over which the field

points range, and other parameters of the calculation are given

in an input file BALIN.XXX, where XXX is the qualifier that

4

identifies a particular run. The initial values of t and x as

well as the increment of any one or two of these variables are

specified in this file. The output file BALOUT. XXX contains the

values of the components of E and cB.

The integration is performed by summing the contributions to

the field components of the integrands calculated at the center

of each patch multiplied by the area of the patch.

Different versions of the program have been tried. Some o-f the variations change the way the patches on the sphere are

determined (a simpler subroutine can be used at the expense o-f computing time or accuracy) , start from other initial fields, use

double precision variables to do the computations (improves t

v u*

7
accuracy of the results for times significantly larger than the

duration of the initial pulse), or save only the energy density

instead of the six field components (saves disk space). The

program contains a provision that allows us to restart the

calculation without much loss of time if there is an

i nterrupti on

3.2 Subroutine F'O I.NJ

This subroutine is used to read the information on the

region of space and time where the fields are to be calculated,

-4 ana to generate successive values of x and t.

In its present form, the program requires either that one

variable change over a given range, where the given number of

points is to be distributed uniformly, or that two variables

change, with points distributed over a rectangular grid uniformly

in each direction. A simple variation allows the program to

compute the fields at points located diagonally in the xy—plane

by varying both x and y simultaneously.

3.3 Subroutine IQS

The values of the initial fields that contribute to the

4

fields at a point x at time t are all on the Information—

-4 Collecting Sphere (ICS) centered at the field point x and with a

radius ct.

If the initial fields have the form given by eqs. (13) and

(14), with a <£(z) specified in (IS) or (19), the integrands that

contribute significantly to the result are confined to a finite

8
cylinder whose axis is the z —axis, whose radius is a -function of

the decay constant k of the gaussian, and which is limited by the

planes z = 0 and z = z , where z is a function of the decay

m m

constant a of the double exponential. In particular, we have set

2

k ix^ + y*") and az equal to a maximum value given by the variable

CUT in the input file. We thus have to find the regions of the

ICS that lie inside the cylinder. Depending on the radius of the

sphere and the location of its center, the whole spherical

surface may be included, one or two regions may be inside the

cylinder, or none at all. Once the nature and the limits of

these regions are determined, each is covered by patches

determined by first dividing the regions into strips by circles

of constant polar angle Q and then dividing each strip or partial

strip into patches by circles of constant azimuthal angle 0. The

sizes of these patches are determined in two different ways. A

length given through the input variable DPATCH is used to

determine the magnitudes of the sides of the patches except when

this length is comparable to the radius of the sphere. If this

is the case, there may be too few patches to assure accurate

values for the integral, and a maximum size of the patches is

determined by dividing the sphere into N1 strips and dividing the

equatorial strip into N2 patches, N1 and N2 being values supplied

also through the input file.

After the above preliminary computations are done for a given time and field point, successive calls to the subroutine

ICS return values of the components of the unit vector to the

center of a patch and the solid angle subtended by the patch.

9
given by

AO = 2 a# sin© sin(?iA©). (20)

We choose a circular cap o-f angular diameter equal to A© as the

patch at either of the poles? the correspond! ng solid angle is

Afl = 4tt sin^OiA©). (21)

n

These values might be approximated by sin©A©Atf and n(?iA©)'

w

, respectively, if the accuracy needed allows this substitution.

We nave not determined precisely when this is allowable, but we

found that, in a particular case, the switch made a significant

difference.

Subroutine ICS calls the subroutines PATCH to determine the

components of the unit vector for a patch on a partial strip,

RPATCH to do likewise for a full strip, and FILL to transfer

these values and the solid angle to COMMON. These subroutines

are listed together with ICS.

3.4 Pr og^r am BALPLT

We use this program to produce plots of the energy densi ty

the el ectromagnetic field as a function of the one or two

variables that change over a given interval. The input comes

from the files BALIN.XXX and BALOUT.XXX, and the output is

directed to a plotter or pr i nter /pi otter

When only one variable changes, the data on the dependent

and independent variables, as well as labels, are passed to the

subroutine DRAW4, a general purpose plotting subprogram described

in reference C21. When two variables change, the number of

10
points in either variable is reduced to a maximum of 100 if

necessary, and then the data are passed to subroutine PL0T3D,

which produces three-dimensional plots and is described in

reference [-33.

The nature of the Green function for the wave equation is

such that the fields propagate exactly with the speed of light.

Thus, the field in the xz —plane at a time t is confined to a ring

of radius ct and of a width equal to the size of the initial

pulse. When the grid of points where the fields are calculated

and plotted is superimposed on this ring, the appearance of the

resulting three-dimensional plot is that of a large number of

isolated peaks if the distance between the points is comparable

to the width of the ring. This problem can be minimized by

choosing a density of points that is high enough.

This program can be easily adapted to show all components of

the electromagnetic field, or to show the energy density when

only this quantity is saved in the file BALOUT.XXX.

3.5 Program BALDSK

This program was developed to prepare a raster file of

values proportional to the energy density of the electromagnetic

field that can be shown on the monitor of the image processing

facility or photographed. We assume that the time t is fixed,

and that the grid extends over a range of values of z and

positive values of x. We first interpolate in the x-direction to

256 points and then reflect this set of values about the z-axis

to take advantage of the symmetry in the field, and then w<>
interpolate to 512 such records in the z-direction. The

intensity of a single graph is scaled to values between 0 and

255, and the maximum intensity of a particular graph in a series

for varying values of t is scaled to fit a range between 150 and

255

3.6 Program BALCHK

To use BALCHK, the fields have to be computed for a fixed

value of t and for values of z and positive x that cover the

region where the fields are significantly different from zero.

BALCHK performs a numerical integration to compute the values of

the total energy and momentum of the pulse; these values allow us

to check the validity of the calculations by comparing them to

the initial values obtained from the analytical expressions (15),

(16), and (18) or (19). The x- and y-components of the momentum

vanish? the z-component of the momentum and the energy density

are proportional to

dV (E

2 1 F 1 +

8(ec

2 - 4<x/3 + ff 2

) +1 1 4k

(

a + £)

5 40°’-'

’ ( 22 )

J

-4 ? '?-+'? dV(E^ + c^B^)

o o

8 (g

2 - 4ctff + £

2

) (a + 0)

5

64

(a + $ ) 5

( 23 )

12
To compute these integrals numerically at later times, we

divide the space where we compute the -fields into tori o-f square

cross sections based on the grid in the xz —plane, and we multiply

the volume o-f each torus by the average of the energy density or

Poynting vector as computed at the four corners of the square in

the grid.

Since all the contributions to the energy are positive, if

the grid is too coarse (that is, if it misses much of the ring

where the fields are concentrated) , the agreement can be poor

(that is, the discrepancy can be larger than a few percent)

without the values of the field being incorrect. In these cases,

the conservation of momentum generally shows better agreement

than the conservation of energy.

The initial values of the z—component of the Poynting vector

(15) are positive, so that the pulse as a whole moves in the

positive z-direction.

If only the energy density of the field is saved, we can

still check for the conservation of energy.

4. Sample Outputs

We use a pri nter /pi otter to obtain the outputs from the

program BALPLT. In fig. 1 we show three—di mensi onal plots of the

energy density of the simple pulse given by eq. (19) as

-12 functions of x and z at the initial time and after 1.2x10 s,

-4 which corresponds to a distance of 3.6x10 m. The pulse moves

mainly forward in the z—direction and to the sides, and the pealj 0 20 intensity decreases from 3.3x10 to 5.9x10 in arbitrary

13
uni ts.

Figure 2 shows the same features for a modulated pulse

specified by eq. (18). Although the wave number of the

modulation is not much bigger than the size of the pulse (k =

5 —1 4 —1

1x10 m and a = 2x10 m ), the shape of this pulse changes

little over the chosen period of time, the motion is mostly in

the z-direction and the maximum intensity decreases from 0.74 to

0.46 only. The superposition of the platting grid on the ridges

that belong to the energy density function produces additional

peaks in the plot. Figure 3 shows the x— , y— , and z-components

of the fields E and cB for the unmodulated pulse. They are shown

-4 -4

at the point (1x10 , 1x10 , 0) as functions of time between 0

-12 -12

and 1.2x10 s, and the same fields at the time 1.2x10 s for

-4

y = 0 and z = 2x10 as a function of x between the values

-4 -4

—6x10 and 6x10 m (this plot exhibits the symmetry about

the z—axis: and B^ are antisymmetric while B_, is symmetric).

Figure 4 shows two photographs taken with a special camera

of an image produced from a finer grained computation of the

energy density shown in fig. 1, with the fields calculated on a

256x256 grid. This camera can produce a moving picture by

photographi ng a sequence of frames.

14
Fig, 1. Energy density of the simple pulse as a function of

^ 2

and z at the initial time and after 1.2x10 s.

cylindrical symmetry about the z-axis. The scale of distances

shown in meters and the maximum intensities in arbitrary units.

x

1 s

1 s

15
fig. 1 for a modulated

1

Fig. Same as pul se.
CO


 ZT

Fig. 3. Plots of the components of the electromagnetic

fields? the curves are labeled from 1 to 6 to identify E^, E.,,

E_., and B-., respectively. Top! functions of time at the

point (1x10 1x10 0) . Bottom: functions of x at a time

— 12 —4

t = 1.2x10 for fixed y = 0 and z = 2x10 m.

17
Fig 4. Photographs of energy densities shown in fig. 1.

The pulse is propagating mainly to the left. Symmetry is used to

extend the graph to negative values of x.

IQ
Ref erences

Cl] Marx, E. Integral Equations for Transient Electromagnetic

Fields. NBS Technical Note 1157, February 1982.

C2] Noon, T. V. User’s Manual for Modular Anal ysi s—Package

Libraries ANAPAC and TRANL. Harry Diamond Laboratories

Technical Report TR-1782, November 1976.

C3] Nelson, D. L. Perspective Plotting of Two-Dimensional

Arrays PL0T3D, a Computer Program for a Digital Plotter,

University of Maryland Technical Report 553, March 1966.

19
0£E®Odi_x

In this appendix we provide listings of the programs and

subroutines BALL, POINT, ICS, BALPLT, BALDSK, and BALCHK.

These programs use some general purpose software provided

with the Interdata computer system or developed for plotting. We

do not discuss or list these subroutines here? they are SYSIO,

PLOT, SYMBOL, DRAW4, and PL0T3D.

20
non

non

PROGRAM BALL

C

C MICRO AND OPTICAL METROLOGY GROUP

C NATIONAL BUREAU OF STANDARDS

C

C EGON MARX 11/29/83

C

C THIS PROGRAM IS USED TO COMPUTE THE FIELDS OF A SPATIALLY LOCALIZED

C PULSE AT A TIME T AND AT PRESCRIBED FIELD POINTS WHEN THE

C INITIAL FIELDS ARE GIVEN AS SOLENOIDAL FUNCTIONS OF A VECTOR

C VARIABLE (X» Y» Z).

C

C R 1 » R2, R3 = COMPONENTS OF THE UNIT VECTOR TO A PATCH

C SOLANG = SOLID ANGLES SUBTENDED BY A PATCH

C ELI , EL2, EL3 = ARRAYS TO STORE COMPONENTS OF THE ELECTRIC FIELD

C CBM1, CBM2, CBM3 = ARRAYS TO STORE COMPONENTS OF THE MAGNETIC FIELD

C CAPPA = DECAY CONSTANT FOR GAUSSIAN IN X- AND Y-DIRECTIONS

C CUT = CUTOFF PARAMETER FOR INTEGRAND

C ALPHA = RISE CONSTANT FOR DOUBLE EXPONENTIAL IN Z-DIRECTION

C BETA = DECAY CONSTANT FOR DOUBLE EXPONENTIAL IN Z-DIRECTION

C K = (REAL) WAVE NUMBER FOR MODULATION IN Z-DIRECTION

C (NO MODULATION FOR K = 0. ) C DPATCH = APPROXIMATE SIZE OF PATCH

C NI, NJ = NUMBER OF POINTS IN GRID FOR CALCULATION

C

DIMENSION ELI (512), EL2< 512) , EL3( 512) , CBM1 ( 512) , CBM2<512>, CBM3(512)

LOGICAL EX

REAL K, KZP

CHARACTER QAL*3, FL1*9, FL2#10

COMMON /RCOM/ R, R 1 , R2, R3, SOLANG

COMMON /ICOM/ N1,N2, NX

COMMON /PARAM/ ALPHA, BETA, CAPPA, CUT, DPATCH, P I , TOP I , P ITO, FOP I ,

C

INITIALIZE VALUES

NR=0

C=3. E8

P IT0=2. *ATAN< 1. ) P 1=2. *PITO

TOP 1=4. #P ITO

FOP 1=8. *PITO

FOP IM1=. 5/TOPI

READ IN VALUES OF CONSTANTS FOR THE PULSE

OPEN ( 1 , FILE= 'BALL. QAL 7

, STATUS= 'OLD 7) READ (1,3) QAL, IFL

FL1='BALIN. 7 //QAL

FL2= 'BALOUT. 7 //QAL

OPEN ( 3, F I LE=FL 1 , STATUS= 7 OLD 7 ) READ ( 3, 1 ) CAPPA, CUT, ALPHA, BETA, K, DPATCH

READ ( 3, 2 ) NI , NJ

IF(NI. GT. 512. OR. NJ. GT. 512) CALL EXIT(l)

NJ4=NJ*4

NJ44=NJ4+4

21
ooo

ooo

oooooo

ooo READ FILE QUALIFIER AND CHECK FOR RESTART

I NQU I RE ( F I LE=FL2 , EXIST=EX)

IF( IFL. EQ. 1. AND. EX) THEN

OPEN < 2* FILE=FL2, STATUS= 7 OLD ' , ACCESS= 'DIRECT '

,

1 FQRM= 'UNFORMATTED ', RECL=NJ4>

INQUIRE ( 2, SI ZE=NR > NX=NR/6

IF(NX. GE. NI) THEN

CLOSE ( 2

)

STOP

ENDIF

ELSE

OPEN ( 2, FILE=FL2, STATUS= 'RENEW ' , RECL=NJ4, BLOCKS I ZE=NJ44,

1 FORM= 'UNFORMATTED', ACCESS= 'DIRECT '

)

ENDIF

CLOSE ( 2

)

GET INITIAL PULSE DATA

CALL POINT < DUM, DUM, DUM, DUM, NI, NJ, . TRUE. ) COMPUTE THE FIELDS FOR EACH TIME AT THE GIVEN POINT

DO 350 I=NX+1 , NI

DO 300 J=l, NJ

FIND THE COORDINATES OF THE FIELD POINT AND THE TIME

CALL POINT<X, Y, Z, T, NI, NJ, . FALSE. ) IF ( T. NE. 0. . OR. Z. GE. 0. ) GO TO 110

COMPUTATIONS FOR T=0

XP=X

YP=Y

ZP=Z

XP2=XP**2

YP2=YP**2

ZP2=ZP**2

CXP2=CAPPA*XP2

CYP2=CAPPA*YP2

CX2Y2=CXP2+CYP2

AZP=ALPHA*ZP

IF < AZP. LT. -25. . OR. CX2Y2. GT. 25. ) GO TO 110

EAZ=EXP ( AZP

)

EBZ=EXP < BETA-w-ZP

)

DEXZ=EAZ—EBZ

DEXZP=ALPHA*EAZ-BETA*EBZ

IF < K. EQ. O. ) THEN

FZ=ZP2*DEXZ

FPZ=ZP2*DEXZP+2. *ZP*DEXZ

ELSE

KZP=K*ZP

22 ,
ooo

ooo

ooo

oooo»-ooo

SKZ=SIN <KZP

)

CKZ=COS(KZP)

FZ=DEXZ*SKZ

FPZ=DEXZP*SKZ+DEXZ*K*CKZ

END IF

TCXYP=2. *CAPPA*XP*YP

EC X2Y2=EXP < -C X2Y2

)

E1=YP*ECX2Y2*FPZ

E2=-XP*ECX2Y2*FPZ

E3=0.

CB1=XP*ECX2Y2*FPZ

CB2=YP*ECX2Y2*FPZ

CB3=2. *<CX2Y2-1. )*ECX2Y2*FZ

GO TO 210

110 E1=0.

E2=0

E3=0.

CB 1=0.

CB2=0.

CB3=0.

IF ( T. EQ. 0. ) GO TO 210

ADD THE CONTRIBUTIONS FROM EACH PATCH ON THE UNIT SPHERE

CALL ICS<X, Y, 1 , T, *200, *210)

COMPUTE THE COMPONENTS OF THE VECTOR

FROM THE SOURCE POINT TO THE FIELD POINT

RU1=—R

RU2=-R2

RU3=-R3

SLJ=SOLANG

RV1=R*RU1

RV2=R*RU2

RV3=R*RU3

COMPUTE THE COORDINATES OF THE SOURCE POINT

XP=X—RV1

YP=Y-RV2

ZP=Z-RV3

COMPUTE THE INTEGRANDS

XP2=XP**2

YP2=YP**2

ZP2=ZP**2

CXP2=CAPPA*XP2

CYP2=CAPPA*YP2

CX2Y2=CXP2+CYP2

AZP=ALPHA*ZP

FIND DOUBLE EXPONENTIAL AND ITS DERIVATIVES

23
oorj

ooo

non

EAZ»EXP<AZP>

EBZ«EXP(BETA*ZP>

DEXZ«EAZ-EBZ

DEXZP«ALPHA*EAZ~BETA*EB z DEXZPP=ALPHA**2*EAZ-BETA**2*EBZ

IF(K. EQ. 0. ) THEN

FZ«ZP2*DEXZ

FPZ«ZP2*DEXZP+2. #ZP*DEXZ

FPPZ»ZP2*DEXZPP+4. *ZP*DEXZP+2. *DEXZ

ELSE

KZP=K*ZP

SKZ»SIN<KZP>

CKZ»C0S(KZP)

FZ=DEXZ*SKZ

FPZ=DEXZP*SKZ+DEXZ*K*CKZ

FPPZ=<DEXZPP-K**2*DEXZ)*SKZ+2. *DEXZP*K*CKZ

END IF

COMPUTE FIELDS AT SOURCE POINT

TCXYP*2. *CAPPA*XP*YP

EC X2Y2-EXP ( -CX2Y2

)

E01®YP*ECX2Y2*FPZ

E02=-XP*ECX2Y2*FPZ

E03«0.

CB01=XP*ECX2Y2*FPZ

CB02«YP*ECX2Y2*FPZ

CB03=2. * ( C X2Y2-1 . )*ECX2Y2*FZ

COMPUTE GRADIENTS OF FIELDS AT SOURCE POINT

DE01 1=-ECX2Y2*TCXYP*FPZ

DE02 1 =EC X2Y2* ( 1 . -2. *CYP2)*FPZ

DE03 1 »EC X2Y2*YP*FPP

DE012=-ECX2Y2*( 1. -2. *CXP2)*FPZ

DE022=-DE01

DE032=-ECX2Y2*XP*FPPZ

DE013=0.

DE023-0.

DE033®0.

DB011«-DE012

DB02 1 “-DE022

DB03 1

s= -DE032

DB012~DB021

DB022«DE021

DB032«DE031

DB013=a ECX2Y2*2. *CAPPA»XP*<2. -CX2Y2>*FZ

DB023«ECX2Y2*2. *CAPPA*YP*<2. -CX2Y2)*FZ

DB033=s ECX2Y2*2. *<CX2Y2-1. >*FPZ

COMPUTE DIRECTIONAL DERIVATIVES ALONG R

R DEO 1 =R V 1 *DEO 1 1 +RV2*DE02 1 +RV3*DE03

RDEQ2=RV1*DE012+RV2*DE022+RV3*DE032

RDE03=RV1*DE013+RV2*DE023+RV3*DE033

24
ooooo

pj

oor>

ooo

R DB 0 1 =R V 1 *DBO 1 1 +RV2*DB02 1 +RV3*DB03

RDB02=RV 1 *DB012+RV2*DB022+RV3*DB032

RDB03=RV 1 *DBO 1 3+RV2*DB023+RV3*DB033

COMPUTE VECTOR PRODUCTS

EVRO 1 =E02*RU3-EG3*RU2

EVR02=E03*RU 1 -EG 1 *RU3

EVR03=E0 1 *RU2-E02*RU

CBVR01=CB02*RU3-CB03*RU2

CBVR02=CB03*RU1-CB01*RU3

CBVR03=CB01*RU2-CB02*RU1

RDEVR 1=RDE02#RU3

RDEVR2=-RDE0 1 *RU3

RDEVR3=RDE01*RU2-RDE02*RU1

RDBVR 1 =RDB02*RU3-RDB03*RU2

RDBVR2=RDB03*RU1-RDB01*RU3

RDBVR3=RDB01*RU2-RDBQ2*RU1

ADD CONTRIBUTIONS TO THE FIELDS AT THE FIELD POINT

00

DEl=SLJ*<E01-RDE01+<2. *CBVR01-

DE2=SLJ# ( E02-RDE02+ < 2. *CBVR02-

DE3=SLJ# ( E03—RDE03+ <2. *CBVR03-

DB1=SLJ# ( CB01 —RDB01 — ( 2. *EVR01-

DB2=SLJ# ( CB02-RDB02- ( 2. *EVR02-

DB3=SLJ* ( CB03—RDB03— ( 2. *EVR03-

E1=E1+DE1

E2-E2+DE2

E3=E3+DE3

CB1=CB1+DB1

CB2=CB2+DB2

CB3=CB3+DB3

GO TO 120

E1=E1*F0PIM1

E2=E2*F0PIM1

E3=E3*F0PIM1

CB1=CB1#F0PIM1

CB2=CB2*F0PZM1

CB3=CB3*F0PIM1

RDBVR 1 ) ) RDBVR2) ) RDBVR3 ) ) RDEVR1 ) ) RDEVR2) ) •RDEVR3) )

SAVE THE COMPONENTS OF THE FIELDS

IF MAGNITUDE LESS THAN l.E-39, SET EQUAL TO ZERO TO AVOID

UNDERFLOW WHEN SQUARING

210 IF(ABS<E1 ). LT. 1. E-39) E1=0.

IF(ABS(E2). LT. 1. E-39) E2=0.

IF(ABS<E3). LT. 1. E-39) E3=0.

IF<ABS(CB1 ). LT. 1. E-39) CB1=0.

IF ( ABS ( CBS ) . LT. 1. E-39) CB2=0.

IF(ABS(CB3>. LT. 1. E-39) CB3=0.

ELI < J)=E1

EL2 ( J ) =E2

EL3 ( J ) =E3

CBM1 < J)=CB1

25
CBM2< J)=CB2

CBM3< J)=CB3

300 CONTINUE

C

C SAVE ARRAYS

C

OPEN (2, FXLE=FL2, STATUS" 'OLD ' < ACCESS* "DIRECT

1 FORM* ' UNFORMATTED 7, RECL=NJ4

)

CALL WRITUN(EL1, NJ, NR+1

CALL WRITUN<EL2> NJ< NR+2)

CALL WRITUN(EL3, NJ, NR+3)

CALL WRITUN<CBM1, NJ, NR+4)

CALL WRITUN<CBM2, NJ» NR+5)

CALL WRITUN<CBM3, NJ, NR+6)

NR=NR+6

CLOSE < 2

)

350 CONTINUE

12

3

STOP

FORMAT (E12. 5)

FORMAT (414)

FORMAT <C3, /, II > END

C

SUBROUTINE WR ITUN ( A, N, IR

)

C

C

C

C

THIS SUBROUTINE IS USED TO WRITE UNFORMATTED DATA TO

A DIRECT ACCESS FILE

DIMENSION A ( N

WRITE(2, REC=IR) A

RETURN

END

26
nan

non

SUBROUTINE POINTCX, Y, Z, T, NI, NU, INF)

C

C MICRO AND OPTICAL METROLOGY GROUP

C NATIONAL BUREAU OF STANDARDS

C

C EGON MARX 11/29/03

C

C THIS SUBROUTINE RETURNS A SEQUENCE OF FIELD POINTS

C TO THE MAIN PROGRAM BALL.

C

LOGICAL INP

DATA DX, DY, DZ, DT, NT Is NT2/4*0. , 2*0/

COMMON /RCOM/ R, Rt, R2, R3, SOLANG

COMMON /ICON/ N1,N2, NX

COMMON /PARAM/ ALPHA, BETA. CAPPA, CUT, BPATCH, P I , TOP I , P ITO, FOP I , C

GO TO 200 TO READ NEW INPUT DATA IF "INP" IS TRUE.

IF (INP) GO TO 200

PASS NEW COORDINATES

X = XX

Y=YY

Z=ZZ

T=TT

NJ1=NJ1+1

C CHECK FOR END OF ROW

IF ( NJ1 . LE. NJ ) THEN

C INCREMENT VARIABLE

I F ( NT 1 . EG. 1 ) THEN

XX=XX+OX

ELSE IFCNT1. EQ. 2) THEN

YY=YY+DY

ELSE IFCNT1.EQ. 3) THEN

ZZ— ZX+DZ

ELSE

TT“TT-8”DT

END IF

ELSE

IN J 1 = 1 C RESET VARIABLE TO INITIAL VALUE

IF < NT 1 . EQ. 1 ) THEN

XX=X0

ELSE IF < NT 1. EG. 2) THEN

YY=YQ

ELSE

ZZ=ZQ

END IF

C INCREMENT OTHER VARIABLE

IFCNT2. EQ. 2) THEN

YY=YY+DY

ELSE XFCNT2. EQ. 3) THEN

ZZ=ZZ+DZ

ELSE

27
o

o

o

o

o

TT=TT+DT

END IF

END IF

RETURN

STATEMENTS TO READ IN FIELD-POINT COORDINATES,

THE INITIAL TIME, NUMBER OF THETA STRIPS,

AND MAXIMUM NUMBER OF PATCHES ON STRIP

200 READ<3,

> XX, YY, ZZ, TT

C SAVE INITIAL VALUES

XO=XX

YO=YY

Z0=ZZ

READ ( 3, 3 ) Nl, N2

C READ PULSE SIZE AND DURATION

READ(3, 1 > DELX, DELY, DELZ, DELT

NJ1 = 1 C CHECK FOR PROPER GRID VALUES

IF ( NI . LE. 0. OR. NJ. LE. 1 ) CALL EXIT(7)

C COMPUTE COORDINATE INCREMENTS

NN=NJ

IF (DELT. NE. 0. ) THEN

DT=DELT

/

( NN-1

)

NN=NI

NT 1=4

NT2=4

END IF

IF (DELZ. NE. 0. ) THEN

DZ=DELZ/ (NN-1

)

NN=NI

NT 1=3

NT2=MAX0 ( NT2, 3)

END IF

IF (DELY. NE. 0. ) THEN

DY=DELY/ (NN-1

)

NN=NI

NT 1=2

NT2=MAX0 ( NT2, 2)

END IF

IF (DELX. NE. 0. > THEN

DX=DELX/(NN-1

)

NT1 = 1 NT2=MAX0 ( NT2,

) END IF

C STOP IF ALL DELX ARE ZERO

IF ( NT1 . EG. 0 ) THEN

WRITE ( 6, * ) 'ALL INCREMENTS ARE ZERO'

CALL EXIT ( 8

)

END IF

C GIVE WARNING IF ONLY ONE INCREMENT IS GIVEN BUT NI IS NOT 1 IF ( NT1 . EG. NT2. AND. NI. GT. 1 > WRITE<6, *)

1 'WARNING: ONLY ONE INCREMENT GREATER THAN ZERO'

IF ( NX. EG. 0 ) RETURN

C

28
o

a

260

13

ADVANCE COORDINATES IF RESTARTING

IF ( NI . EQ. 1 ) CALL EXIT<9)

DO 260 1=1, NX

IF(NT2. EQ. 2) THEN

YY=YY+DY

ELSE IF(NT2. EQ. 3) THEN

ZZ=ZZ+DZ

ELSE

TT=TT+DT

END IF

CONTINUE

RETURN

FORMAT ( E12. 5)

FORMAT < 214)

END

29
SUBROUTINE ICS< X, Y, Z, T, *, *

)

C

C MICRO AND OPTICAL METROLOGY GROUP

C NATIONAL BUREAU OF STANDARDS

C

C EGON MARX i 1/29/83

C

C THE FOLLOWING GROUP OF SUBROUTINES IS USED TO COMPUTE A PATCH

C COVERING FOR THE INFORMATION-COLLECTING SPHERE, AND THEY RETURN

C TO THE MAIN PROGRAM BALL THE VALUES OF THE COMPONENTS OF THE

C UNIT VECTOR TO A PATCH AND THE SOLID ANGLE SUBTENDED BY THE PATCH

C

C DPATCH = APPROXIMATE LENGTH OF A SIDE OF A PATCH

C N1 = NUMBER OF THETA STRIPS (ALTERNATIVE)

C N2 = MAXIMUM NUMBER OF ANGLES PHI (ALTERNATIVE)

C IP - POINTER FOR PATCHES

C IRP = POINTER FOR PATCHES

C

DIMENSION PMIN ( 2 ) , PMAX ( 2 ) , PRP (2 )

,

RPMIN(2) , RPMAX (2)

COMMON /RCOM/ R, Rl, R2, R3, SOLANG

COMMON /ICOM/ N1,N2, NX

COMMON /PARAM/ ALPHA, BETA, CAPPA, CUT, DPATCH, PI, TOPI, PITO, FOP I, C

DATA LBL/O/

GO TO ( 150, 200, 250), LBL

IP=0

IRP=0

G= 1 80. /PI

C COMPUTE RADIUS OF CYLINDER WHERE INITIAL DATA ARE SIGNIFICANT

RHO=SQRT ( CUT/CAPPA

)

C COMPUTE COORDINATE OF BOTTOM OF CYLINDER

ZM*3 —CUT /ALPHA

C COMPUTE RADIUS OF INFORMATION-COLLECTING SPHERE

R=C*T

C CHECK THAT Z IS SUCH THAT THERE MAY BE AN INTERSECTION

IF ( Z. LT. R. AND. Z. GT. ZM-R ) GO TO 70

RETURN 2

C COMPUTE DISTANCE OF FIELD POINT TO Z-AXIS

70 RP=SGRT(X**2+Y**2)

C COMPUTE AZIMUTH ANGLE OF FIELD POINT

IF ( X. EG. 0. . AND. Y. EQ. 0. ) THEN

PHI0=0.

ELSE

PHI 0=ATAN2 ( -Y, -X)

END IF

C COMPUTE ANGULAR SIZE OF PATCHES

DTH=BPATCH/R

NP=P I /DTH

C COMPARE TO MINIMUM NUMBER OF STRIPS

IF ( NP. LT. N1

)

DTH=PI/N1

C COMPUTE PATCH SIZE IN PHI-DIRECTION

DL2=DPATCH

IF (2-H-NP. LT. N2) DL2=T0P I*R/N2

DTH2=DTH*. 5

C COMPUTE THETA FOT TOP AND BOTTOM OF CYLINDER

THETAT=0.

30
THETAB=P

IF ( ABS ( Z ) . LE. R> THETAT=ACOS < -Z/R ) +1 . E-5

IF ( ABS < ZM—Z ) . LE. R) THETAB=ACOS < ( ZM-Z > /R

)

C COMPUTE DISTANCE FROM FIELD POINT TO CYLINDER

RPR=ABS(RP-RHO)

C CHECK IF FIELD POINT IS OUTSIDE CYLINDER

IF(RP. GT. RHO) THEN

C CHECK FOR INTERSECTION

IF(R LE. RPR) THEN

RETURN 2

ELSE

C COMPUTE THETA FOR INTERSECTION

THETA 1=ASI N < RPR/R

)

THETA2=P I -THETA

C CHECK FOR SECOND INTERSECTION

IF(R. LE. RP+RHO) THEN

C CHECK POSITION OF TOP AND BOTTOM OF CYLINDER

IF ( THETA2. LE. THETAT. OR. THETA 1 . GE. THETAB ) THEN

RETURN 2

ELSE

C SAVE RANGE OF THETA FOR PARTIAL STRIP

IP=IP+1

PMIN( IP )=AMAX1 (THETA1, THETAT)

PMAX ( IP ) =AMIN1 < THETA2, THETAB)

PRP < IP )=RP

END IF

ELSE

C COMPUTE ANGLES FOR SECOND INTERSECTION

THETA3=AS I N ( ( RP+RHO ) /R

)

THETA4=P I—THETA3

C CHECK POSITION OF TOP AND BOTTOM OF CYLINDER

IF (THETA2. LE. THETAT. OR. THETA 1 . GE. THETAB. OR.

1 (THETA3. LE. THETAT. AND. THETA4. GE. THETAB) ) THEN

RETURN 2

ELSE

IF< THETAT. LT. THETA3) THEN

C SAVE RANGE OF THETA FOR PARTIAL STRIP

IP=IP+1

PMIN( IP >=AMAX1 (THETAT, THETA1

)

PMAX ( IP )=AMIN1 (THETAB, THETA3

)

PRP ( IP >=RP

IF ( THETAB. GT. THETA4) THEN

C SAVE RANGE OF THETA FOR PARTIAL STRIP

IP=IP+1

PMIN(IP) =THETA4

PMAX ( IP ) =AMIN1 (THETAB, THETA2

)

PRP (IP )=RP

END IF

ELSE

C SAVE RANGE OF THETA FOR PARTIAL STRIP

IP=IP+1

PMIN ( IP >=AMAX1 (THETAT, THETA4

)

PMAX ( IP) =AMIN1 (THETAB, THETA2

)

PRP ( IP ) =RP

END IF

31
END IF

END IF

END IF

C IF THE FIELD POINT IS INSIDE THE CYLINDER

ELSE

C CHECK IF CYLINDER DOES NOT INTERSECT SPHERE IN TWO CURVES

IF

(

R. LE. RP+RHO) THEN

C CHECK IF SPHERE IS INSIDE CYLINDER

IF(R.LE. RPR) THEN

C SAVE RANGE OF THETA FOR FULL STRIP

IRP=IRP+1

RPMIN< IRP )*=THETAT

RPMAX ( IRP )=THETAB

C IF SPHERE INTERSECTS CYLINDER IN ONE CURVE

ELSE

C FIND EXTREME ANGLES

THETA 1=ASIN< RPR /R

)

THETA2-P I -THETA

C CHECK POSITION OF TOP OF CYLINDER

IF< THETAT. LT. THETA 1 ) THEN

C SAVE RANGE OF THETA FOR FULL STRIP

IRP=IRP+1

RPMIN< IRP )=THETAT

RPMAX ( IRP) =AMIN1 ( THETA

, THETAB

)

END IF

IF (THETAT. LT. THETA2. AND. THETAB. GT. THETA 1 > THEN

C SAVE RANGE OF THETA FOR PARTIAL STRIP

IP=IP+1

PMIN< IP )=AMAX1 < THETA li THETAT)

PMAX( IP)=AMIN1 (THETA2, THETAB)

PRP ( IP ) ——RP

END IF

IF ( THETAB. GT. THETA2) THEN

C SAVE RANGE OF THETA FOR FULL STRIP

IRP=IRP+1

RPMIN ( IRP)=AMAX1 (THETA2, THETAT)

RPMAX ( IRP) =THETAB

END IF

END IF

C IF SPHERE INTERSECTS CYLINDER IN TWO CURVES

ELSE

C COMPUTE TWO EXTREME ANGLES

THETA3=ASIN< < RP+RHO ) /R

)

THETA4=P I —THETA3

C CHECK POSITION OF TOP AND BOTTOM OF CYLINDER

IF (THETAT. GE. THETA3. AND. THETAB. LE. THETA4 ) THEN

RETURN 2

ELSE

C COMPUTE OTHER TWO ANGLES

THETA 1=ASIN (RPR /R

)

THETA2=P I -THETA

C CHECK TOP OF CYLINDER

IF (THETAT. LT. THETA 1 ) THEN

C SAVE RANGE OF THETA FOR FULL STRIP

IRP=IRP+1

32
RPMIN < IRP ) s=THETAT

RPMAX ( IRP )=AMIN1 (THETAB, THETA1

)

END IF

IF ( THETAT. LT. THETA3. AND. THETAB. GT. THETA1 ) THEN

C SAVE RANGE OF THETA FOR PARTIAL STRIP

IP=IP+1

PMIN( IP )=AMAX1 (THETAT, THETA 1

)

PMAX< IP )=AMIN1 (THETAB, THETA3

)

PRP ( IP ) =-RP

ENDIF

IF (THETAT. LT. THETA2. AND. THETAB. GT. THETA4 ) THEN

C SAVE RANGE OF THETA FOR PARTIAL STRIP

IP=IP+1

PMIN( IP)=AMAX1 (THETAT, THETA4)

PMAX( IP >=AMIN1 (THETAB, THETA2)

PRP ( IP ) =-RP

ENDIF

IF(THETAB. GT. THETA2) THEN

C SAVE RANGE OF THETA FOR FULL STRIP

IRP=IRP+1

RPMIN( IRP >=AMAX1 (THETAT, THETA2

)

RPMAX ( IRP )=THETAB

ENDIF

ENDIF

ENDIF

ENDIF

C END OF GEOMETRICAL CONSIDERATIONS BETUEEN SPHERE AND CYLINDER

C CHECK FOR PARTIAL STRIPS

140 IF( IP. EQ. 0) GO TO 190

C INDICATE PROCESSING OF PARTIAL STRIP

LBL=1

C GET DATA FOR STRIP

THMIN=PMIN( IP)

THMAX=PMAX( IP)

PR=PRP( IP)

C DECREMENT POINTER

IP=IP-1

C COMPUTE COMPONENTS OF UNIT VECTOR TO PATCH

150 CALL PATCH(THMIN, THMAX, DTH, DL2, RHO, PR, PHIO, *140)

RETURN

C CHECK FOR FULL STRIPS

190 IF ( IRP. EQ. 0 ) GO TO 210

C INDICATE PROCESSING OF FULL STRIP

LBL=2

C GET DATA FOR STRIP

THMIN=RPMIN( IRP)

THMAX=RPMAX( IRP)

C DECREMENT POINTER

IRP=IRP—

C COMPUTE COMPONENTS OF UNIT VECTOR TO PATCH

200 CALL RPATCH(THMIN, THMAX, DTH, DL2, *190)

RETURN

C INDICATE NO MORE STRIPS

210 LEL=3

RETURN

33
o

o

oo

C ALL DONE FOR THIS FIELD POINT

250 LBL=0

RETURN 1 END

C

SUBROUTINE PATCH < THMIN, THMAX, DTH, DL2, RHQ, RP, PHIO, *)

THIS SUBROUTINE COMPUTES THE COMPONENTS OF THE UNIT VECTOR

TO A PATCH IN A PARTIAL STRIP

COMMON /RCOM/ R, Rl, R2, R3, SOLANG

COMMON /PARAM/ ALPHA, BETA, CAPPA, CUT, DPATCH, P I , TOP I , P ITQ, FOP I , C

DATA LBL/O/

GO TO (100, 110, 120), LBL

C CHECK FOR FIELD POINT ON AXIS OF CYLINDER (NO PARTIAL STRIPS)

IF ( RP. EQ. 0) RETURN 1 C COMPUTE RANGE IN THETA

DELTH=THMAX-THMIN

C COMPUTE NUMBER OF STRIPS

NN1=DELTH/DTH

DT1=DELTH/ (NN1+1

)

DT2=DT1*. 5

B = . 5/RP

A= ( RH0**2-RP**2 ) *B

TH=THM I N+DT2

C COMPUTE SOLID ANGLE SUBTENDED BY PATCH

SLG=2. *DL2*SIN(DT2)/R

100 ST=SIN(TH)

RST=R*ST

CT=COS(TH)

C COMPUTE SIZE OF PATCH IN AZIMUTHAL DIRECTION

0PHI=DL2/RST

C COMPUTE RANGE IN AZIMUTH FOR PARTIAL STRIP

DELPHI=ACQS ( B*RST-A/RST > NN2=DELPHI /DPHI

DPHI=2. *DELPHI/(2*NN2+1

)

PHI 1=PHI0

PHI2=PHI0

C SAVE COMPONENTS OF UNIT VECTOR AND SIZE OF CENTER PATCH

CALL FILL(ST, CT, PHIO, SLG)

1 = 1 IF(NN2. EQ. 0) GO TO 130

C INDICATE ADDITIONAL PATCHES

LBL=2

RETURN

C INCREMENT ANGLES

110 PHI 1=PHI 1+DPHI

PHI2=PHI2—DPHI

C SAVE COMPONENTS OF UNIT VECTOR AND SIZE OF PATCH

CALL FILL(ST, CT, PHI1, SLG)

C INDICATE SECOND ANGLE

LBL=3

RETURN

C SAVE COMPONENTS OF UNIT VECTOR AND SIZE OF PATCH

120 CALL FILL(ST, CT, PHI2, SLG)

34
o

o

r>n

1 = 1 + 1 C CHECK IF DONE WITH STRIP

IFd.LE. NN2) THEN

C INDICATE BACK TO FIRST ANGLE

LBL=2

RETURN

ENDIF

C INCREMENT THETA

130 TH=TH+DT

C CHECK IF DONE

IF(TH. LE. THMAX) THEN

C INDICATE START OF NEXT STRIP

LEL=1

RETURN

ENDIF

C INDICATE ALL DONE

LBL=0

RETURN 1 END

C

SUBROUTINE RPATCH ( THMIN, THMAX* DTH, DL2, *)

THIS SUBROUTINE COMPUTES THE COMPONENTS OF THE UNIT VECTOR

TO A PATCH IN A FULL STRIP

COMMON /RCOM/ R, R 1 , R2, R3, SOLANG

COMMON /PARAM/ ALPHA, BETA, CAPPA, CUT, DPATCH, P I , TOP I , P ITO. FOP I ,

C

DATA LEL/O/

GO TO (80, 90, 100, 110), LBL

IF ( THMIN. EQ. 0. ) THEN

C PROCESS CAP AT TOP OF SPHERE

DT2=AMIN1(DTH*. 5, THMAX

)

S1=F0PI*SIN(DT2*. 5)**2

CALL FILL(0. , 1. , 0. , SI ) IF ( DT2. EQ. THMAX) RETURN 1 THMIN=DT2

LBL=1

RETURN

ENDIF

80 IF (THMAX. EQ. PI ) THEN

C PROCESS CAP AT BOTTOM OF SPHERE

DT2=AMIN1 (DTH*. 5, PI-THMIN)

Si=F0PI*SIN(DT2*. 5>**2

CALL FILL ( 0. , -1. , 0. , SI

)

THMAX=P I—DT2

IF ( THMIN. GE. THMAX) RETURN 1 LBL=2

RETURN

ENDIF

C COMPUTE RANGE IN THETA

90 DELTH=THMAX-THM I

NN1=DELTH/DTH

DT1=DELTH/ (NN1+1

)

DT2=DT1*. 5

TH=THMIN+DT2

35
o

oo

o

100 ST=SIN<TH>

RST=R*ST

CT=COB < TH

)

C COMPUTE NUMBER OF PATCHES IN STRIP

NN2=MAX0< INT<T0PI*RST/DL2+1. ), 3)

C COMPUTE SIZE OF PATCH IN AZIMUTHAL DIRECTION

DPHI=T0PI/NN2

PHI—O.

C COMPUTE SOLID ANGLE SUBTENDED BY PATCH

SLG=2. *DPHI*ST*SIN(DT2)

1 = 1 C INDICATE NEXT PATCH

LBL=4

C SAVE COMPONENTS OF UNIT VECTOR AND SIZE OF PATCH

110 CALL FILL < ST, CT, PHI , SLG

)

C INCREMENT ANGLE

PHI=PHI+DPHI

1 = 1 + 1 IF ( I . LE. NN2 ) RETURN

TH=TH+DT1

IF(TH. LT. THMAX ) THEN

C INDICATE END OF STRIP PROCESSING

LBL=3

RETURN

ENDIF

C INDICATE ALL DONE

LBL=0

RETURN 1 END

C

SUBROUTINE FILL ( ST, CT, PHI , SLG

)

THIS SUBROUTINE COMPUTES AND STORES THE COMPONENTS OF THE

UNIT VECTOR TO THE PATCH, AND THE SIZE OF THE PATCH

COMMON /RCOM/ R, R 1 , R2, R3, SOLANG

R1=ST*C0S<PHI

R2=ST*SIN<PHI

)

R3=CT

SOLANG=SLG

RETURN

END

36
oo

o

o

oon

noo

PROGRAM BALPLT

C

C MICRO AND OPTICAL METROLOGY GROUP

C NATIONAL BUREAU OF STANDARDS

C

C EGON MARX 11/29/83

C

C THIS PROGRAM IS USED TO MAKE A THREE-DIMENSIONAL PLOT

OF THE ENERGY DENSITY FROM THE OUTPUT OF PROGRAM BALL IF TWO

COORDINATES VARY, OR A PLOT OF EACH COMPONENT OF THE ELECTROMAGNETIC

C FIELD IF ONLY ONE COORDINATE VARIES.

C

DIMENSION EB( 100, 100), XLAB<2, 4), VAR (512)

DIMENSION ELI (512), EL2(512), EL3(512), CBM1 (512), CBM2(512), CBM3(512)

INTEGER PRLINE(20)

CHARACTER QAL*3, FL1*10, FL2*9

EQU I VALENCE ( EB , VAR

)

DATA XLAB/ 7 X—COORD. ', 'Y-COORD. ', 'Z-COORD. ', ' TIME '/

INITIALIZE THE PLOTTING ROUTINES

READ ( 5, « ) IANG1, IANG2

THETA=I ANG1

PHI=I ANG2

READ ( 5, 3) QAL

FL1 = 'BALOUT. ' //QAL

FL2= 'BALIN. 7 //QAL

OPEN <4, FILE=FL2, STATUS= 'OLD

) READ ( 4, 2 ) NI,NJ

CHECK FOR VALID VALUES FOR PLOTTING ANGLES

IF ( NI . GT. 1. AND. < IANG1. EQ. 0. OR. IANG2. EQ. 0) ) THEN

WRITE (6,*) 'SUPPLY NONZERO VALUES FOR THETA AND PHI'

CALL EXIT(l)

END IF

NJ4=NJ*4

ND I=(NI— i)/i 00+

NDJ=(NJ-1 )/ 100+1

READ (4,1) XO, YO, ZO, TO

READ (4, 4) IDUM, IDUM

READ (4, 1 ) DELX, DELY, DELZ, DELT

OPEN (2, F I LE=FL 1 , STATUS= ' OLD '

)

SELECT VALUES FOR PLOTTING

READ IN DATA FILES

11=0

DO 110 1=1, NI

CALL READUN (ELI, NJ

)

CALL READUN ( EL2, NJ)

CALL READUN ( EL3, NJ)

CALL READUN ( CBM1, NJ)

CALL READUN ( CBM2, NJ)

CALL READUN ( CBM3, NJ)

37
100

110

120

C

C

C

130

140

150

160

170

180

C

C

C

IF ( NI . EQ. 1 ) GO TO 200

I F ( MOD ( I - 1 , ND I > . NE. 0 ) GO TO 110

11=11+1

JJ=0

DO 100 J=l, NJ, NDJ

JJ3* JJ+ i EB(II, JJ)=EL1 ( J ) **2+EL2 ( J ) **2+EL3 ( J ) **2+

1 CBM1 ( J) **2+CBM2 < J ) *#2+CBM3 ( J ) **2

CONTINUE

CONTINUE

UMAX=0.

DO 120 1 = 1, II

DO 120 J=l, JJ

UMAX=AMAX 1 ( UMAX , EB ( I , J )

)

CONTINUE

PREPARE A THREE-DIMENSIONAL PLOT

CALL PLOTS < 0, 0, 0)

ENCODE (PRL I NE, 13) UMAX

CALL SYMBOL ( . 5, 9. 2, . 1, PRLINE, 0. , 36)

IF(DELY. NE. 0. . AND. DELT. NE. 0. > GO TO 130

IF ( DELZ. NE. 0. . AND. DELT. NE. 0. ) GO TO 140

IF (DELY. NE. 0. . AND. DELZ. NE. 0. ) GO TO 150

IF ( DELX. NE. 0. . AND. DELZ. NE. 0. ) GO TO 160

IF ( DELX. NE. 0. . AND. DELY. NE. 0. ) GO TO 170

CALL SYMBOL ( . 5, 10. , . 15,

1 'ENERGY DENSITY AS A FUNCTION OF X AND T',0. ,39)

ENCODE (PRLINE, 7) XO, YO, ZO, TO, DELX, DELT

GO TO 180

CALL SYMBOL (. 5, 10. , . 15,

1 'ENERGY DENSITY AS A FUNCTION OF Y AND T',0. ,39)

ENCODE (PRLINE, 8) XO, YO, ZO, TO, DELY, DELT

GO TO 180

CALL SYMBOL (. 5, 10. , . 15,

1 'ENERGY DENSITY AS A FUNCTION OF Z AND T',0. ,39)

ENCODE (PRLINE, 9) XO, YO, ZO, TO, DELZ, DELT

CALL SYMBOL ( . 5, 10. , . 15,

1 'ENERGY DENSITY AS A FUNCTION OF Y AND Z',0. ,39)

ENCODE (PRLINE, 10) XO, YO, ZO, TO, DELY, DELZ

GO TO 180

CALL SYMBOL ( . 5, 10. , . 15,

1 'ENERGY DENSITY AS A FUNCTION OF X AND Z',0. ,39)

ENCODE (PRLINE, 11) XO, YO, ZO, TO, DELX, DELZ

GO TO 180

CALL SYMBOL (. 5, 10. , . 15,

1 'ENERGY DENSITY AS A FUNCTION OF X AND Y',0. ,39)

ENCODE (PRLINE, 12) XO, YO, ZO, TO, DELX, DELY

CALL SYMBOL ( . 5, 9. 5, . 1, PRLINE, O. , 75)

CALL PLQT3D(EB, 0. , 0. , 100, 100, 1, 1, II, 1, 1, JJ, 1, THETA, PHI

)

CALL PLOT ( 0. , 0. , 999)

STOP

IF ONLY ONE COORDINATE CHANGES, PREPARE MATERIAL FOR PLOT

38
0->

U

IU

O

O

O

O

200 IF ( DELX. NE. 0. ) THEN

V=X0

DV=DELX / ( NJ- 1

)

NTV=1

END IF

IF (DELV. NE. 0. ) THEN

V=Y0

DV=DELY/(NJ-1

)

NTY=2

END IF

IF (DELZ. NE. 0. ) THEN

V=Z0

DV=DELZ/(NJ-1

)

NTY=3

END IF

IF ( DELT. NE. 0. ) THEN

V=TO

DV=DELT / < NJ— 1

NTY=4

END IF

DO 210 1=1, NJ

VAR ( I

) =V

V=V+DV

210 CONTINUE

ENCODE (PRLINE, 6) XO, YO, ZO, TO

CALL DRAW4( 1, 7, 2, 4, 16, 17, XLAB ( 1, NTY) , 'FIELD AMPLITUDES '

,

1 'ELECTROMAGNETIC FIELDS, 1 : El, 2: E2, 3: E3, 4:CB1, 5: CB2,

2 ', PRLINE)

PREPARE ONE PLOT WITH ALL OF THE 6 FIELD COMPONENTS AS A

FUNCTION OF THE INDEPENDENT VARIABLE

7

8

9

10

CALL DRAW4 < 2, 7, 1, NJ, 7 1 100, VAR, ELI, 0. , 0. ) CALL DRAW4 ( 2, 7, 1, NJ, 7 2

7

, 110, VAR, EL2, 0. , 0. ) CALL DRAW4<2, 7, 1, NJ, '3', 120, VAR, EL3, 0. , 0. ) CALL DRAW4 (2, 7, 1, NJ, '4', 130, VAR, CBM1, 0. , 0. > CALL DRAW4<2, 7, 1, NJ, 7 5

7

, 140, VAR, CBM2, 0. , 0. ) CALL DRAW4<2, 7, 1, NJ, 7 6

7

, 150, VAR, CBM3, 0. , 0. ) CALL DRAW4 ( 3, 7, 0, 0, 2, -NJ, ELI, EL2, 2. , 2. ) CALL PLOT < 0. , 0. , 999)

STOP

FORMAT ( E12. 5)

FORMAT (//////, 214)

FORMAT (C3)

FORMAT < 214)

FORMAT

(

7 INITIAL VALUES: X0= , ,1PE8. 1,

1

7

, T0=', E8. 1

)

FORMAT < 7 X0= 7, 1PES. 1, 7 Y0= 7 ,E8. 1,

1

7 DELX= 7

» ES. 1, ' DELT= '

, E8. 1

)

FORMAT <

7 XG= 7

, 1PE8. 1, ' Y0= 7, E8. 1,

1

7 DELY= 7

, E8. 1 , 7 DELT= 7, E8. 1

)

FORMAT <

7 X0= 7

, 1PE8. 1, 7 Y0= 7 ,EB. 1,

1

7 DELZ= 7

* E8. 1 , 7 DELT= 7 ,E8. 1)

FORMAT (

7 XG= 7

, 1PE3. 1, 7 Y0= 7 ,E8. 1,

1

7 DELY= 7

, E8. 1 , 7 DELZ= 7

, E8. 1

)

Y0=

Z0= 7, E8. 1,

Z0= 7, E8. 1,

Z0= 7

, E8. 1,

Z0= 7, E8. 1,

E8. 1, 7

, Z0=

T0= 7

, E8. 1.

T0= 7, E8. 1,

T0= 7

, E8. 1,

T0= 7

, E8. 1,

6: CB3

, E8. 1,

39
o

o

r>

11 FORMAT (

' XO- ' > 1 PE8. 1 , ' Y0-SE8. It' Z0-SE8.li' TO- S ES. 1,

1

' DELX-SE8.lt' DELZ-SE8. 1)

12 FORMAT < 'XO-S 1PE8. 1, ' Y0-SE8.li' Z0-SE8.lt ' T0-SE8. 1,

1

' DELX-SES.lt' DELY-SE8. 1)

13 FORMAT ( 'MAXIMUM ENERGY DENSITY IS ', 1PE10. 3)

END

C

SUBROUTINE READUN(A, N>

THIS SUBROUTINE IS USED TO READ UNFORMATTED DATA

DIMENSION A(N)

READ ( 2 ) A

RETURN

END

40
PROGRAM BALDSK

C MICRO AND OPTICAL METROLOGY GROUP

C NATIONAL BUREAU OF STANDARDS

C

C EGON MARX 1 1/29/83

C

C THIS PROGRAM PREPARES FILES TO BE PROCESSED BY THE I2S

C FACILITY TO BE PUT ON A MONITOR OR PHOTOGRAPHED.

C THE INPUT IS THE FILE PRODUCED BY THE PROGRAM BALL.

THE OUTPUT IS A 512X512 SET OF INTEGERS THAT REPRESENT THE

C SCALED INTENSITY OF THE PICTURES.

C WHEN A SERIES OF OUTPUT FILES IS PROCESSED SIMULTANEOUSLY,

C THE MAXIMUM INTENSITY OF EACH OUTPUT IS SCALED

C

DIMENSION EB ( 512) , IEB ( 128) , EB1 ( 512) , EB2< 512)

DIMENSION ELI (512), EL2(512), EL3( 512) , CBM1 ( 512 ) , CBM2< 512 ) , CBM3 ( 512

DIMENSION AMX(IOOO)

LOGICAL EX

CHAR AC TER*12 FL

DATA IFL/O/

OPEN (5, FTLE= 'BALDSK. INP 7

, STATUS31 'OLD 7

) READ (5,*) N!,N2

IF(N1. GT. N2. OR. N2. GE. 10Q0) CALL EXIT<254)

IF ( N1 . EQ. N2 ) GO TO 265

C GET OLD FILE OF MAXIMA, IF IT EXISTS

INQUIRE (FILE- "BALDSK. SAV 7

, EXIST=EX > IF (EX) THEN

OPEN (3, FILE= 'BALDSK. SAV 7

, STATUS* 'OLD 7> GO TO 240

END IF

C OR CREATE A NEW ONE IF NOT

OPEN (3, FILE* 'BALDSK. SAV 7

, 3TATUS= 'NEW ',

1 FORM* ' UNFORMATTED 7

, RECL=4, BL0CKSIZE=8>

C OPEN A FILE FOR OUTPUT OF MAXIMA

OPEN ( 1 , F I LE= ' BALDSK. MAX', STATUS* 7 RENEW 7

) AMMAX=G.

ALMAX=1 . E20

C GET THE MAXIMA FROM THE DIFFERENT FILES

DO 200 1 = 1, 1000

FL= BALOUT, 7 //IT0C (I-i, IDUM)

OPEN (2, FILE=FL, STATUS= 7 OLD 7, ERR=220, SHARE= 7 SR0 7 ) C FIND THE SIZE OF THE FILES

IFdFL. EQ. 0) THEN

IFL=1

INQUIRE (2, SIZE=NI, RECL=NJ4)

NI=NI/6

NJ=NJ4/4-l

END IF

AMAX=Q.

DO 110 11=1, N

CALL READUN(EL1, NJ)

CALL READUN(EL2, NJ)

CALL READUN(EL3, NJ)

CALL READUN(CEM1, NJ)

41
CALL READUNCCBM2, NJ)

CALL READUN(CBM3» NJ)

DO 100 J=1,NJ

AMAX=AMAX 1 ( AMAX , EL 1 ( J > **2+EL2 < J ) **2+EL3 ( J ) **2+

1 CBtll ( J)#*24-CBM2( J)**2+CBM3< J)**2)

100 CONTINUE

110 CONTINUE

CL0SE<2)

AMX ( I ) -AMAX

IF (AMAX. LE. AMMAX) GO TO 120

AMMAX=AMAX

NN=I

120 IF ( AMAX. G£. ALMAX) GO TO 200

AL!iAX=AMAX

MM—

200 CONTINUE

1=1000

220 WRITE < 3 ) 1-1

DO 230 J=l, 1-1

WRITE ( 3 ) AMX(J)

WRITE (1,4) J— 1 , AMX ( J

)

230 CONTINUE

WRITE (1,5) NN-1, AMMAX

WRITE (1*6) MM-1, ALMAX

WRITE(3) AMMAX

WRITE (3) ALMAX

GO TO 260

C RETRIEVE INFORMATION FROM SAVED FILES

240 READ ( 3 ) I DO 250 J=l, I READ(3) AMX(J)

250 CONTINUE

READ(3) AMMAX

READ (3) ALMAX

260 C0EF=105. / ( AMMAX-ALMAX

)

C CREATE RASTER FILES FOR PICTURES

265 DO 700 NN=N1 + 1 # N2+1

FL= 'BALQUT. '//ITOC (NN-1, IDUM)

OPEN ( 2, FILE=FL* STATUS= ' OLD ' , SHARE='SRO'

)

C FIND INFORMATION ON FILE SIZE

IF ( IFL. EQ. 0 ) THEN

IFL=1

INQUIRE (2, SI ZE=NI » RECL=NJ4>

NI=NI/6

NJ=NJ4/4-l

ENDIF

FL= ‘'BALDSK. '//IT0C(NN-1, IDUM)

OPEN (4, FILE=FL, FORM= ' B I NARY ' » STATUS= 7 RENEW '

,

1 BLOCKS I ZE=512

)

WRITE (6, 11) NN-1

C FIND MAXIMUM FOR SINGLE PICTURE

IF ( N1 . EQ. N2 ) THEN

AMAX=Q.

DO 2S0 11=1, NI

CALL READUN ( EL 1 , NJ

)

42
270

280

C

C

C

300

C

310

320

C

330

C

c

340

CALL READUN <EL2, NJ)

CALL READUN <EL3, NJ)

CALL READUN < CBM 1, NJ)

CALL READUN(CBM2, NJ)

CALL READUN <CBM3» NJ)

DO 270 J= 1 , N J

AMAX=AMAX1 <AMAX, ELI ( J)**2+EL2< J)**2+EL3< J)**2+

1 CBM 1 < J ) **2+CBM2 < J ) **2+CBM3 < J ) **2

)

CONTINUE

CONTINUE

FAC=255. /AMAX

REWIND 2

ELSE

FAC= ( COEF* ( AMX ( NN ) -ALMAX ) + 1 50. ) /AMX < NN

)

ENDIF

READ IN FIRST RECORD

CALL READUN < EL 1 , NJ

)

CALL READUN ( EL2, NJ)

CALL READUN < EL3, NJ)

CALL READUN ( CBM1, NJ)

CALL READUN ( CBM2i NJ)

CALL READUN ( CBM3, NJ)

LIMIT THE NUMBER OF POINTS TO 256 FOR SYMMETRY

NJ1=MIN0 < NJ, 256)

SCALE INTENSITIES

DO 300 J=1,NJ1

EB ( J ) —FAC* < EL 1 ( J ) **2+EL2 < J ) **2+EL3 ( J ) **2+

1 CBM1 ( J)**2+CBM2< J)**2+CBM3< J)**2)

CONTINUE

INTERPOLATE FIRST RECORD

I J=257

EB02=EB ( 1

)

DO 320 XI=2, NJ1

EBQ 1=EB02

EB02=EB (II)

K 1 = < 5 1 2- 1 J ) / < NJ 1 + 1 — 1 1

)

FC— 1. /KI

DO 310 K=0, Kl-1

EB2( I J)=FC*(EB01*<K1-K)+EB02*K)

XJ=XJ+1

CONTINUE

CONTINUE

EB2< IJ)=EB02

FILL IN THE REFLECTED VALUES ON RECORD

DO 330 11=1, 256

EB2 (II) =EB2( 513-1 I

)

CONTINUE

PROCEED WITH OTHER RECORDS, INTERPOLATING IN OTHER DIRECTION

KJ=1

DO 430 1=2, N

SHIFT RECORD BETWEEN BUFFERS

DO 340 K=l, 512

EB1 ( K ) =EB2 ( K

)

CONTINUE

IJ=257

43
c

350

C

360

370

C

380

C

C

390

400

C

C

C

420

430

C

READ ANOTHER RECORD

CALL READUN (ELI* NJ

)

CALL READUN ( EL2, NJ > CALL READUN ( ELS, NJ)

CALL READUN ( CBM1, NJ)

CALL READUN ( CBM2, NJ)

CALL READUN (CBM3, NJ)

SCALE NEW RECORD

DO 350 J=1 , NJ1

EB ( J ) =FAC* ( EL 1 ( J ) **2+EL2 < J ) **2+EL3 < J ) #*2+

1 CBM 1 ( J ) **2+CBM2 ( J ) **2+CBM3 < J > **2

)

CONTINUE

INTERPOLATE NEW RECORD

EB02=EB ( 1

)

DO 370 11=2, NJ1

EB01=EB02

EB02=EB (II)

Kl = ( 512— I J) / (NJ1+1— I I

)

FC=1. /K1

DO 360 K=0, Kl-1

EB2( IJ)=FC*(EB01*(K1-K)+EB02*K>

IJ=IJ+1

CONTINUE

CONTINUE

EB2( I J)=EB02

FILL IN REFLECTION OF NEW RECORD

DO 380 11=1,256

EB2( II )=EB2( 513-1 I

)

CONTINUE

INTERPOLATE BETWEEN RECORDS

Ll = < 512—KJ )/(NI + l — I

)

FC=1 . /LI

DO 420 L=0, Ll-1

IK=1

DO 400 11=1, 128

IVAL=0

PACK VALUES IN BYTES

DO 390 K=1 , 4

I VAL= I VAL*256+ 1 NT ( FC* ( EB 1 < I K ) * ( L 1 -L ) +EB2 ( X K > *L > > IK=IK+1

CONTINUE

IEB (II ) =XVAL

CONTINUE

DO SPECIAL UNFORMATTED I/O TO WRITE TO RASTER FILE

CALL SYSIO( IPB, 57, 4, IEB, 512, 0)

KJ=KJ+1

CONTINUE

CONTINUE

DO LAST RECORD

I K= 1 DO 500 11=1, 128

IVAL=0

DO 450 K=l, 4

44
non

o

*-• o

ui

vi

non

IVAL=IVAL*256+INT(EB2< IK)

)

IK=IK+1

450 CONTINUE

IEB (II ) =X UAL

500 CONTINUE

DO SPECIAL UNFORMATTED I/O TO WRITE TO RASTER FILE

CALL SYSIO< IPB, 57, 4, IEB, 512, 0)

CLOSE (2)

00 CONTINUE

STOP

FORMAT (

' FILE ', 13,', MAX= ' , 1PE12. 5)

FORMAT < //, ' FILE ', 13, LARGEST MAX=', 1PE12. 5)

FORMAT < //, ' FILE ',13,', SMALLEST MAX= ', 1PE12. 5

)

1 FORMAT( IX, 13)

END

SUBROUTINE READUNCA, N)

THIS SUBROUTINE IS USED TO READ UNFORMATTED DATA

DIMENSION A(N)

READ < 2

)

A

RETURN

END

45
oooo

o

o

o

o

oooooooooo

PROGRAM BALCHK

MICRO AND OPTICAL METROLOGY GROUP

NATIONAL BUREAU OF STANDARDS

EGON MARX 11/29/83

THIS PROGRAM CHECKS THE CONSERVATION OF ENERGY AND MOMENTUM

OF THE ELECTROMAGNETIC FIELD BY COMPARING THESE QUANTITIES

AT TIME T WITH THEIR INITIAL VALUES.

DIMENSION EBH512), EB2< 512), EXCB 1(512 >,EXCB2( 512)

DIMENSION ELI (512), EL2(512>, EL3(512>, CBMK512), CBM2012), CBM3012)

REAL K, K2

CHARACTER FL1*12, FL2*12, QAL*3

DATA IFL/O/

P 1=4. *ATAN( 1. ) TOP 1=2. *PI

OPEN FILE WHERE THE QUALIFIERS OF FILES TO BE PROCESSED ARE STORED

OPEN ( 1 # FILE* 'BALCHK. QAL', STATUS* 'OLD '

)

OPEN FILE TO SAVE THE OUTPUT

OPEN (3, FILE* ' BALCHK. OUT ' » STATUS* 'RENEW '

)

READ (1,5, END=400 ) QAL

READ PULSE PARAMETERS

FL 1 = ' BALIN. 7 //QAL

OPEN (4, FILE=FL1, STATUS* 'OLD '

)

READ (4, 1 ) CAPPA, DUM, ALPHA, BETA, K, DUM

READ(4, 2) NI,NJ

READ(4, 1 ) XO, YO, ZO, TO

READ(4, 2) I DUM, I DUM

READ (4, 1 ) DELX, DELY, DELZ, DELT

DX=DELX/(NJ-1

)

DZ=DELZ/ ( NI-1

)

IF ( IFL. EQ. 0 ) THEN

IFL=1

COMPUTE INITIAL VALUES OF THE ENERGY AND MOMENTUM OF THE

FIELDS FROM THE ANALYTIC EXPRESSIONS.

C2=CAPPA*#2

A2=ALPHA**2

B2=BETA**2

AB=ALPHA*BETA

APB*ALPHA+BETA

IF ( K. EQ. O. ) THEN

ERG=TQPI*(. 25/C2* ( . 25/ALPHA**3+8. *(A2-4. *AB

1 +B2)/APB**5+. 25/BETA**3>

2 +. 375/CAPPA*( 1. /ALPHA#-* 5—64. /APB**5+1. /BETA**5> ) EXB=P I#. 25/C2*(. 25/ALPHA**3+8. *(A2-4. *AB

1 +B2)/APB**5+. 25/BETA**3)

ELSE

K2=K**2

APB2=APB**2

ERG=T0PI*K2*(. 25/C2*(. 25/ALPHA—4. * ( AB+K2 ) / ( APB* ( APB2+4. *K2))

46
non

non

250

260

270

300

310

+ . 25/BETA ) +1 . /CAPPA*(. 25/ ( ALPHA# ( A2+K2 )

)

2 -4. / ( APB* < APB2+4. *K2 ) > +. 25/ ( BETA* ( B2+K2 > ) >

)

EXB=. 25*PX*K2/C2*( . 25/ALPHA-4. *<AB+K2)/(APB*<APB2+4. *K2>

)

1 +. 25/BETA)

END IF

WR ITE < 3# 4 ) ERG# EXB

END IF

FL2='BAL0UT. V/QAL

READ IN FIELD VALUES AND COMPUTE ENERGY DENSITY AND POYNTING VECTOR

OPEN (2, FILE-FL2# STATUS- "OLD ' > DO THE NUMERICAL INTEGRATION TO FIND THE ENERGY AND MOMENTUM

ERG=0.

EXB=0.

CALL READUN ( EL 1 # NJ

)

CALL READUN ( EL2# NJ)

CALL READUN ( EL3# NJ)

CALL READUN ( CBM1 # NJ)

CALL READUN ( CBM2# NJ)

CALL READUN < CBM3, NJ)

DO 250 J=l, NJ

EB2 ( J ) =EL 1 ( J ) **2+EL2 ( J ) **2+EL3 ( J ) **2+

1 CBM1 ( J ) **2+CBM2 ( J ) **2+CBM3 ( J ) **2

EXCB2 < J ) =EL1 < J ) *CBM2 ( J ) —EL2 < J ) *CBM1 (J)

CONTINUE

DO 310 1=1, N I -

DO 260 J=l, NJ

EB1 < J)=EB2( J)

EXCB1 (J)=EXCB2<J>

CONTINUE

CALL READUN < EL 1 # NJ > CALL READUN (EL2, NJ)

CALL READUN ( EL3, NJ)

CALL READUN ( CBM1 , NJ

)

CALL READUN ( CBM2# NJ)

CALL READUN ( CBM3, NJ)

DO 270 J=l, NJ

EB2 < J ) =EL 1 < J ) **2+EL2 ( J ) **2+EL3 ( J ) **2+

1 CBM1 < J)**2+CBM2< J)**2+CBM3< J)**2

EXCB2 ( J ) =EL1 ( J ) *CBM2 ( J ) -EL2 ( J ) *CBM1 < J

)

CONTINUE

X1=X0

DO 300 J~1 # NJ—

X2=X1+DX

SURF=X2**2—X 1**2

ERG=ERG+. 25* < EB 1 ( J)+EB2( J)+EB1 ( J+l )+EB2( J+l ) )*SURF

EXB-EXB+. 25*(EXCB1 ( J)+EXCB2( J)+EXCB1 ( J+l )+EXCB2( J+l ) )*SURF

X 1=X2

CONTINUE

CONTINUE

ERG=ERG*P I *DZ

EXB=EXB*PI*DZ

47
ooo

o

cn

4*u

ju

UIRITE<3, 3) FL1 * ERG» EXB

GO TO 100

400 STOP

1 FORMAT (E12. S)

FORMAT ( 214

)

FORMAT( IX, C12, 2X, 3E15. 8)

FORMAT < //, 19X, 'ENERGY', 8X, 'MOMENTUM',/, ' INITIAL VALUES', 2E15. 8)

FORMAT (A3)

END

SUBROUTINE READUN( A, N)

THIS SUBROUTINE IS USED TO READ UNFORMATTED DATA

DIMENSION A(N)

READ ( 2 ) A

RETURN

END

48
NBS-114A (REV. 2-80

U.s. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO.

NBSIR 84-2835

2. Performing Organ. Report No, 3. Publication Date

May 1984

4. TITLE AND SUBTITLE

Free-Space Propagation of Light Pulses

5. AUTHOR(S)

Egon Marx

6. PERFORMING ORGANIZATION (If joint or other than NBS, see instructions)

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

7. Contract/Grant No.

S. Type of Report & Period Covered

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

10. SUPPLEMENTARY NOTES

I | 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,

bibliography or literature survey, mention it here)

If document includes a significant

A transient electromagnetic field in free space is completely specified when the

initial values of the electric and magnetic fields are given. Green's function

for the scalar wave equation can then be used to find the field at later times.

A group of computer programs that implement these equations and process the

output are presented in this report.

12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semi color s computer programs; conservation of energy and momentum; electromagnetic field

propagation; finite light pulses; graphic displays of energy density; transient

electromagnetic fields

13. AVAILABILITY

22 Unl imited

| | For Official Distribution. Do Not Release to NTIS

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

20402.

BTI Order From National Technical Information Service (NTIS), Springfield, VA. 22161

14. NO. OF

PRINTED PAGES

53

15. Price

$10.00

UlCOuu-OC


V-