NBSIR 83-2788

Piece-Wise Analytic Evaluation of

the Radiative Tail From Elastic

and Inelastic Electron Scattering

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Center for Radiation Research

Washington, DC 20234

Nuclear Physics Laboratory

Department of Physics

University of Illinois at Urbana-Champaign

Urbana, Illinois 61801

December 1983

er> w* u s department of commerce

1 3 J NATIONAL BUREAU OF STANDARDS

,U56

33-2733

1933

L 9-

NBSIR 83-2788

PIECE-WISE ANALYTIC EVALUATION OF

THE RADIATIVE TAIL FROM ELASTIC

AND INELASTIC ELECTRON SCATTERING

UATIOIIAL BU7K'. OF STAMDARDS

LlLliAitY L ' /'C . O c

/ 6<y

U

! C/K^

Leonard C. Maximon

U S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Center for Radiation Research

Washington, DC 20234

Steven E. Williamson

Nuclear Physics Laboratory

Department of Physics

University of lllionis at Urbana-Champaign

Urbana, Illinois 61801

December 1983

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

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director
I
PIECE-WISE ANALYTIC EVALUATION OF THE RADIATIVE TAIL

FROM ELASTIC AND INELASTIC ELECTRON SCATTERING

L. C. Maximon

National Bureau of Standards

Washington, D.C. 20234

S. E. Williamson

Nuclear Physics Laboratory, Department of Physics

University of Illinois at Urbana-Champaign

Urbana, II

1

1 nois 61801

ABSTRACT

We consider here the calculation of the radiative tail from the elastic

peak in medium and high energy electron scattering as well as from a discrete

inelastic level of the recoiling nucleus. We examine the method generally

used for this calculation, viz., a numerical integration of the differential

cross section over the angles of the unobserved photon, and discuss the

difficulties inherent in this numerical integration due to the sharp peaking

of the integrand. We present an alternative method for calculating the

radiative tail, in which the region of integration is divided into an arbitrary number of subintervals, the structure functions are fitted by cubic

spline functions in each subinterval, and the integrations are then performed

analytically in closed form. This method has the advantages of greatly

increased accuracy and a reduction of the computation time by a factor which

can vary between 10 and 10 3

, depending on the kinematics.

Key words: electron scattering cross sections; elastic electron scattering;

elastic radiative tail; inelastic electron scattering; inelastic

radiative tail; intermediate and high energy electron scattering.
ACKNOWLEDGEMENTS

We wish to express our gratitude to Gloria Wiersma, to whom we are

indebted for the skillful and expeditious typing of this report. We thank

David K. Kahaner of the Scientific Computing Division at the National Bureau

of Standards for suggesting the use of the Kronrod rules, embodied in the

subroutines DQK61 and D1MACH. We also thank Claude Schuhl for the gracious

hospitality accorded at the Centre d' Etudes Nucl£aires de Saclay, where much

of this work was done. This work was supported in part by the National

Science Foundation under Grant NSF PHY 80-23603.
CONTENTS

Abstract..... i

Acknowledgements i i

I. Introduction. 1

II. Preliminaries for the Piece-Wise Analytical Evaluation of the

Radiative Tail......... 6

III. Details of the Piece-Wise Analytical Evaluation of the Radiative Tail 16

Appendix (Algebraic Details)..., 31

References 55

Timing Table (Comparison of Piece-Wise Analytic Integration with Numerical

Integration) 56

Computer Program Listing 59

I. Introduction

The analysis of medium and high energy electron scattering experiments

involves, as is well known, the calculation of electromagnetic processes

(e.g., radiative corrections, Landau straggling, radiative tails, etc.) which

inevitably accompany and mask in some manner the electron-nucleus scattering

of primary interest. The methods for evaluating the cross sections for these

processes has been given in numerous articles and reviews familiar to

experimenters, and are generally incorporated into the analysis of electron

scattering experiments. We consider here the calculation of the radiative

tail from the elastic peak, as well as from a discrete inelastic level of the

recoiling nucleus. 1 We examine the method which is generally used for this

calculation, namely a numerical integration of the differential cross section

over the angles of the unobserved photon, and discuss the difficulties

inherent in this numerical integration. Finally, we present an alternative

method for calculating the radiative tail having the advantage of greatly

increased accuracy and a reduction of the computation time by a factor which

can vary between 10 and 10 3

, depending on the kinematics (see Timing Table,

p. 56).

The formalism used for the calculation of the radiative tail in almost

all analyses of high-energy electron scattering experiments is that presented

in a detailed series of papers by Tsai [1-3] 2 and later refined by Miller [4],

In all of these papers the cross section for the scattering of an electron of

initial momentum Ip^l, final momentum |£2

| at an angle 6, with emission of a

photon, is given by an integral of the form

^he assumption that the final state of the recoiling nucleus is a discrete

level does not constitute a significant restriction in practice, since in the

analysis a final continuum state is generally treated as a sum of many

discrete levels.

2 Figures in brackets indicate literature references at the end of this paper.
( 1 ) dnd|£2

l = ^ {w

2f

2 + w

i fi } dn

k

*

Here the functions f

2 and are known analytic functions of the initial and

final electron momenta and scattering angle (which remain fixed) and the

photon angles 0^,^, over which one integrates. The structure functions W

2

and W

1 are, for elastic scattering, a function of the four -momentum transfer

q

2 < 0. They are, in general, qiven as the result of a numerical calculation

in which a charge distribution has been entered into a phase shift

calculation. (That is, although a Born approximation formalism is used to

derive the cross section given in (1), one uses, in the integrand in that

expression, "experimental" structure functions, defined either by the

experimental cross section divided by the Mott cross section or by a phase

shift calculation.) It is expected, therefore, that the calculation of the

radiative tail from (1) requires a numerical integration. However, the

numerical integration of the expression appearing in the integrand in (1) is

not completely trivial. Depending on the kinematics, the functions f

2 and f1 can be very sharply peaked, — a reflection of the fact that the photons are

emitted preferentially in the direction of either the initial or the final

electron. In addition, since the functions f

2 and fx have an overall factor —, the integrand may rise sharply for q

2 near its minimum value, achieved (in

q

4

the coordinate system defined following eqs (9)) when cos 0^ = 1. These peaks

are discussed in detail both in Maximon and Isabelle [5] and in Mo and Tsai

[2], In general, the height of the peaks is of order (E/mc 2

) relative to the

background and the width is of order (mc 2/E) relative to the width of the

region over which one integrates. The contribution from the peaks is thus of

the same order of magnitude as that from the background, and for high energies

2
great care must be taken to ensure that the numerical integration procedure

has indeed converged and gives the desired accuracy. The result of the

peaking is thus that the numerical integration may be quite lengthy, or of

dubious accuracy, or both. Further, given the large number of intervals that

must be chosen for the integration in the region of the peaks, the numerical

evaluation of the structure functions W

2 and W1 in each integration interval

would be prohibitive in computation time, and rather useless since the

functions W

2 and Wx vary little over even the entire region of a peak. The

procedure generally followed 3 is therefore to divide the entire region of

integration into a fixed number of subintervals (1 < i < N) (say with N of the

order of 100) and to fit W

2 and VJ

1 by spline functions in each of these

subintervals, e.g., by a cubic function of -q 2

:

W

2 = a

i + bi (-q 2

) + c

i (-q 2

) 2 + d^-q 2) 3

Wi = a

i + bi (-q

2

) + c

i (-q

2

) 2 + d^-q 2) 3

.

( 2 )

With the overall region of integration divided into these subintervals (some

of which contain the peaks) and W

2 and WL represented by these simple analytic

expressions in each subinterval, the numerical integration is then carried out

as just discussed, the peak rendering the integration costly in time and

questionable in accuracy.

3 This is, essentially, the procedure followed by R. Altemus in the analysis of

the data on 56 Fe taken at Bates (reported in Phys. Rev. Lett. 44, 965

(1980)), by M. Deady in the analysis of the data on 40, 48Ca taken at Bates

(reported in Phys. Rev. C28, 631 (1983)), and by Z. Meziani in the analysis

of the data on 40 ’ 48Ca and 56 Fe taken at Sacla.y (to be published shortly).

3
We note that the integration procedure just described consists in fact of

two parts which are quite independent. In the first, each of the structure

functions is fitted by a spline (or other useful function) in each of the N

subintervals into which the integration region is subdivided for this purpose.

(For example, for a point charge with spin zero we have W

2 = Z

2 and W

1 = 0, so

that we may choose N = 1.) In the second part, the numerical integration of

(1) is performed in each of the N subintervals, with W

2 and W

L now given by

(2)

, using some suitable integration technique (e.g., Simpson's rule or

Gaussian integration). In the event that a peak occurs in a given subinterval, then clearly this subinterval must be further subdivided, and the

convergence of the integration procedure must be tested.

The essence of the present report is contained in the observation that

once we have carried out the fitting of W

2 and W

L by functions of the form

given in (2) (the first part of the procedure just described), there is no

longer any need to perform the integration numerically. With W

2 and W

x given

by finite polynomials in q

2 the functions f

2 and f1 are sufficiently simple

that the integration can be performed in each subinterval in closed form.

(The resulting functions involve only algebraic functions and logarithms.)

The integration over the entire region then consists merely of summing the

contribution from each of the N subintervals. The problems related to these

peaks are thereby avoided completely. The presence of sharp peaks does not

necessitate a further division of the subintervals. It is thus by eliminating

the numerical integration, the second part of the procedure usually followed,

that we reduce the time for the entire calculation by a factor which varies

between 10 and 10 3

, depending on the specific kinematics. The accuracy of the

entire calculation is also greatly increased, since it is now given by the

precision of the computer rather than the accuracy of the particular numerical

integration method.

4
As a last general comment before entering into the details, we observe

that in this note the method, which we call "piece-wise analytic" evaluation,

is applied to the integration of the doubly differential cross section over

the angles of the photon. However, it is equally applicable (with, of course,

different resulting expressions) to the situation of interest for the analysis

of experiments using photons as the electromagnetic probe. Then one integrates this same doubly differential cross section over the angles of the

electron: The expressions which result in this latter case will be presented

in a future report, in which we consider in particular the cross section for

polarized photons (integrated over the angles of the final electron).

We next present the details of the calculation of the radiative tail

given in this report. We divide this presentation into three parts: In the

first we give preliminaries to the method presented here. Notation is

defined, the differential cross section which forms the basis of our work is

presented, and a few comments are made. In the second we give the details of

the calculation as reflected directly in the computer program, a copy of which

is included with this report. 4 Finally, in the Appendix we give details of

the algebra involved in obtaining our final expression for the cross section,

given by eq (38) and the subsidiary equations defining the quantities in it.

4 The Fortran program is available from the authors. Requests should be

accompanied by a magnetic tape and a specification of the desired format.

5
II. Preliminaries for the Piece-Wise Analytical Evaluation of the Radiative

Tail

The cross section under consideration here is the radiative tail to both

elastic and inelastic scattering of electrons from a nucleus. In this report

(as in the work of Tsai) we consider emission of photons only by the electron.

Furthermore, the formalism used is that of the Born approximation. Thus in

terms of diagrams we have

We do not consider here the contribution due to radiation by the recoiling

nucleus. That has been treated in an approximate manner by Miller [4] and

will be examined by us in a report to appear later.

The notation used here is as follows: Energies, momenta, and masses are

all understood to be in energy units. When four-vectors are written, the

metric is as in Bjorken and Drel 1 [6], viz.,

V Y

( 3 )

6
Pi = ( £ 1\P l)

:

P 2 = ( e 2 >P 2 )

:

k = ( w,J<

)

:

q = (q 0 »^) :

P

! l.il) =

* (M,0)

P

2

=

( E

2 ’£.2

>

:

= (M + q 0 ,q)

m =

M = *

M =

four-momentum of the incident electron

four -moment urn of the scattered electron

four-momentum of the emitted (real) photon

four -momentum of the virtual photon exchanged between

the electron and the target particle

four -momentum of the target particle

in the lab system

four -momentum of the recoiling nucleus

in the lab system

rest mass of the electron

rest mass of the target particle.

rest mass of the final, recoiling nucleus

(= (M* - M) Cl + *4^) ” M* - M for 144* « 1). (4)

Then from energy-momentum conservation

p i + P

i = p 2 + P 2 + k • (5)

The momentum transferred to the nucleus is then

q

= P

2 -P 1 = P 1 -P 2 -k . ( 6 )

Further, it is useful to define the four-vector u by

u = (u

0

,_u) = p 1 + P

1 -p 2 = P

2 +k (7)

7
Then in the lab system

U

0 = M + Ej - c

2 (8a)

and

£ = £1 - £2 (8b)

We note that u

0 and £

are independent of both photon energy and angles.

0 = scattering angle of the electron (0 = $ (jP 1 *P. 2

)

)

0^ = angle between £

and k_

0^ = angle between £

and £x

As in Mo and Tsai [2] the integration over photon angles is performed in a

coordinate system such that the z axis is in the direction of £ = £: - £2

. As

observed there, this coordinate system has the particular advantage that q

2 is

then independent of the photon azimuthal angle . The structure functions

W

2 (-q 2

) and W^-q 2) can then be taken outside of the integral over which

can then be performed analytically. We show this explicitely in the following

paragraph.

From the notation just introduced, since the recoiling nucleus is assumed

(9)

*• to be a free particle of mass M , we have

8
( 10 )

M*

2 = P

2

2 = (u - k) 2

= u

2 - 2k *u

from which

0) - IjJ | cos 6

k

J

‘ (ID

Thus the photon energy, w, depends on the photon direction, but in the

particular coordinate system we have chosen it is a function only of 0^, and

not of (Note, however, from the expression for w given below, viz.,

Q

2 Q2 w = - e

£

- y + that in the limit 0, we have w = e

i - e2

* v; that

is, in the absence of recoil the photon energy is independent of photon

direction.)

Similarly, we can write

M*

2 = P

2

2 = (p t + q)

2

+ q

2 + 2q-Pj .

Thus in the lab system

u .

( 12 )

(13)

9
On the other hand, substituting q

= - p 2

- k in 2q*P

1 above gives

( 14 )

so that in the lab system

(15)

Thus in the coordinate system with z-axis in the direction of u., both u> and

q

2 are independent of Only the products p^k and p 2 *k depend on <a .

Our starting point is then the doubly differential cross section for

scattering of an electron with the emission of a photon, essentially as qiven

by Tsai [1] and Mo and Tsai [2], We make two minor changes in the expressions

presented in these references, discussed immediately following the cross

section written below.

d

2 o a

3

|^_2 1 2 f

+1 2Mco d (cos 0^ ) r 2ir

dfld|£2 | (2") 2 ejpj J -1

q

4 (u 0 - |_u |cos 6

k

) J°

where
I } * W

2 (-q 2

) [- [2«

z (e -») + %q2] - [2e,(«2+») + iq 2

] L

< p l'

k > ( p 2- k ) 2

+ ipj-lofivkT

Cm

z (p 1 *p 2

-

2

) + (Pi

, P 2 )[2 e iVPi‘P2t “( £ r

£2)])

- (p

'l

-k- [2(e

1 e r

c 2W e

1 2

) + %q

2

- P^Pj - m

2

]

+

(p -k) [ 2 ( e iVe l"fe

2

2

) + *p2 - p i* p 2

- m2 ]

+ W,(-q 2

) |m 2 (2m 2+q 2) (

1 + 1 )

L

V ( p ! -k ) 2

( p 2 *k ) 2 /

4p 1 *P 2 (p 1 *p

2 -2m 2

) +

(p 1 *k)(p 2 *k)

- (t^FT * T^)( 2p i’ P2+2m2- q2)

+ qj

(i7)

The first modification we have made in the expression of Tsai concerns the

structure functions. We use W

2 and Wx

, defined by the elastic scattering

cross section, which at high energies is given in terms of these functions by

da a

2cos 2 ^6[W 9 + 2tan 2^0 W,] —^ = . --9

- — . (18)

dft 4e 2 sin 4 ^0 [l + £_1 s i

n

2% 0

]

M

11
(Note that for a point charge with spin zero, W

2 (0) = Z

2

, W^O) = 0.) This is

the same as the change of notation introduced by Stein et al. [7] in

Appendix A of their paper. That appendix is a useful reference in that it is

essentially a summary of Tsai's SLAC-PUR 848 [3]. The structure functions F.

J and G. used by Tsai are related to W

2 and Wx by

j

F

j *

G. = 4M 2W

1

. (19)

J

The second modification made in the expression of Tsai concerns the term

ip 2

i 2

. , which appears as a factor to the integral just written. In the

e

2 l£.i I articles of Tsai this factor has been replaced by its high energy limit, viz.,

e—2

. While this replacement is clearly of little significance numerically for

e

i ip 2

i 2

energies of experimental interest, we retain the factor —r—-t out of a sense

e

2 1J1 1

of consistency, since no high energy approximations are made anywhere else in

the expression for the cross section. Our expression is then identical to

that given by Nguyen-Ngoc and Perez-y-Jorba [8], apart from a trivial difference in the normalization of the structure functions: they use F and G, which

are related to our W

2 and Wx by

F M

2

G = - W

1

. ( 20 )

12
In any event, this term is closely related to the physics: Writing

I P.2

1

2 /IP 2 I\ l£2

l IP 2

I —;—

r

= —r, we note that the factor comes from the fact that we

£

2 Ip.i i V

e 2 / Ipj ; a have written the cross section as

d

2

j- rather than . The term |p dftd L

I P 2i ^ “

2 comes from the incident flux, and the final factor |p 2

| comes from the density

of final states for the scattered electron: (Note that in Stein et al . [7]

d

2 o -2^

they write

d ^

- e

-, so that only should appear, but they replace this with

e—2

. ) This factor without high energy approximation has also been retained to

e

i facilitate comparison with the Bethe-Heitler cross section in the limit M

Finally, we note that in Stein et al. [7] the structure functions W

2 and

W x have been replaced by W

2 and Wx

, defined by (A42) and (A43) of Stein

et al.:

W_( q

2

) = F(q 2 )W,(q 2)

WJq 2) = F(q 2 )M

1 (q

2

) (21)

to account for both thick target bremsstrahlung (equivalent radiator) and that

part of the radiative correction which does not depend on the resolution AE.

In addition, in the work of Miller [4] these factors are further modified to

account (in an approximate manner) for the multiple emission of soft photons

from both the electron and the target particle. The important point to retain

here is that all these factors depend on the photon direction only through the

four-momentum transfer q

2 (or via the photon energy, which is again a linear

function of q

2

, as we have just seen). All these effects can therefore be

accounted for simply by an appropriate modification of the structure functions

which leaves them independent of .

13
We next perform the integration over as in Mo and Tsai [2], Tsai [3]

and Stein et al . [7], and obtain

d

2 o a

3 |£ 2

I 2 f

+1 2Ma)d(cos0u )

dnd |p /: 2

| 2v J ~ 1 q

4

(u

0

- McoseJ

( 1 ( 22 )

where

f 1 * W 2 (-q 2

) am 2 rn _ , _ . \ . i a'm2

[2 £

i( e2+ “) + -

, [2 £

2

< £ i-“) +

+ 2v

(i - y){m 2 (p 1 *p

2

- t*32

) + (p 1 *P2)[2 e

1 e2-Pi , P2+aj( e i- e 2)H

+ e

i e2+ei a>fe 2

2

) + - P 1 *P 2

- m2 ]

’ ^[2( e

i e

2

- e

2

a>fe

i 2

) + ^

2

" P 1 *P 2

' m2 l

-]

+ W^-q 2

) j^(j

3 + p) m

2 (2m 2

+q 2

) + 4v (£ - -)(p 1 *p

2 ) (p 1 *p2 -2m 2

)

+ (j - '“)(2p 1 *p

2+2m 2-q 2

) + ^

(23)

14
in which 5 * 6

a = w( e

2

- |£2 |cos0

p cos0

k

)

a'= o)( e

x - |£ 1 |cos0

s cos0

(< )

b = - w|£

2 |sin0

p sin0

k

v = (a' - a)- 1

x = (a 2 - b

2

) 1/2

y = (a 1 2 - b

2

) 1/2

|£x

| COS 0 - |£2

| COs9

p iuj

cos 0„ =

s

l£ x

l - I P.2

1

eos 0

(24)

5The symbol x which is defined in (24) and used in (23) is employed at this

point out of consistency with the notation of Tsai [3] and Stein et al . [7].

It is used nowhere else in this report as defined in (24). Elsewhere

throughout this report we define x

= -q 2

.

6The expressions given here in (24) are those given in Tsai's SLAC-PUB 848 [3]

and Stein et al . [7], They represent a considerable simplification of the

expressions given originally in Mo and Tsai [2],

15
III. Details of the Piece-Wise Analytical Evaluation of the Radiative Tail

At this point we depart from the procedure followed in these references.

Since we perform the remaining integration analytically, the most appropriate

variable of integration is q

2 rather than cos 9^, in view of the functions

W

2 (-q 2

) and W^-q 2 ). Indeed, the form of the integrand itself is then

somewhat simpler. To give just one example, referring to the overall factor

in the integrand in above expressions for the radiative tail, we now have,

with q

2 as integration variable.

We will not, however, encumber this section of the report with these algebraic

manipulations. They are presented at the end in the Appendix. Here we give

only the final form of our expression for the radiative tail, which lends

itself directly to the computer implementation of the calculation.

As variable of integration we choose finally a positive quantity, viz..

2M a) d ( cos 0^) 1 d(q 2

) (25)

q

4 (u

Q - |jj | cos

9

k

) \u\ q

4

x -q 2 > 0 . (26)

The cross section for the radiative tail is then given by

d

2 a a

3

|jd 2

1 2

dftd|p_

2

| tt e^pj

(27)

16
where

1 1= - [

6

0 + SjX + 6

2x

2

] x

2D^ 3/2 (x)

x

2D

2

3 /2 ( x

)

[«o + 5

1X + {2x2]

|u|x 2 (2A - x) VDj l/2 (x) 0

2

1/2(x),

x

2Dj 1/2 (x) x2D

2

1/2 (x)

xD 1/2 (x) xD

2

1/2 (x)

(28)

and where j }' is obtained from { } merely by replacing the parameters

3

0 ,3 1 ,... by the corresponding parameters with primes:

{3 0

, S

2

, 6

0

, 6

1# 52

, n, Y, T p

t

2

, p p p

2 }

f V » V > e

2

' , 6

0

‘ , «i' , 62

‘ , n' , Y' , Ti' , t2

' , p

i' , p

2

'

}

In (28) the quadratic functions D^(x) and ( x) are related to the quantities

a, a', and b of (24) by eqs (A19) and (A25). Referring to the terms in the

1 integrand with factor (2X - x)

indeed lies in the range of integration:

, we note in the appendix that the point x = 2

A

x < 2A < x M , m

( 29 )

17
but that D|(2A) = (2 X) , so that there is in fact no singularity in the

integrand at the point x = 2A. (This is shown explicitely later in this

report, in the Appendix.)

The parameters appearing in the integral in (27) and (28) are defined

shortly below. However, we first outline our piece-wise analytic evaluation

of the cross section. We divide the range of integration (x iX^) into N

subintervals (which can be of different and arbitrary length) writing

(30)

with

x

0

= x

rn

- X

N

= X M

' (31)

In each subinterval ( x

-j_i» x -j) we represent W

2 (x) and WL (x) by a third order

polynomial in x: The coefficients of the polynomial will, in general, be

different in each subinterval.

!n (x

i _ 1 ,x

i

)

W

2 (x) = a- + b-x + c.jX 2 + d-x 3

I I W

x (x) = a- + b-x + c-x 2 + d-x 3 (32)

We may then write

18
/* x m n r

x

i

!j dx[w 2 (x){ } + Wjx)) } ] = M dx[W 2 (x){ } + W

x (x) { } ] *m i=1 Vl

N f

X

i r vv* dx

[

(a

i 1-1 x i-l

b

i x + C-X d

i x

; ) {

+ (a- + brx + c^x 2 + d^x 3

) {

Considering the expressions just given for { } and { } (the functions D^x)

and (x ) are quadratic functions of x), the integrals can now all be

performed analytically. The integrations are all straightforward — none

result in functions any more complicated than the logarithm. However, care

must be taken in arranging the terms, and in particular the arguments of the

logarithms, so that at high energies there are no cancellations between almost

equal terms, leaving a remainder of relative order (m/ej) 2

, with a consequent

loss in precision.

Referring to the integrand, the integrals to be performed analytically

are thus of the the form

(1) * P U Jv

dx

x i-lx2

^i

3/2

( x )

[a^ + b^x + c-x 2 + d^

x

3

] [e

0 + $

x x + 62x

2

] (34a)

'i 2> = / - h + v C iX 2 + d, x

3

] (34b)

x i-lx

is -r

dx

Xi.j x

2Dj 1 / 2 (x)

[a. + b^x + c-x 2 + d -

x

3

] (34c)

19
(4) . f

1 l.i ' l

dx

x i-l

xD

i 1/2

( x )

[a- + + ^x 2 + d-x 3

] (34d)

and

f

dx

x i_l ( 2X - x

) \p 1 1/2

( x ) D

2

I/2

( x)>

[a- + b-x + c^ 2 + d^x 3 ]. (34e)

The integrals iP], iP], and il

4

] follow from iP], iP], and l{

4

! merely by

U.9I L. ^

I l j

I J. ^I J. ^I 1 )

I replacing D^(x) by D^(x) in the integrands (and replacing 3

0 ,$ 1 ,62 in I jP by

6

0 * 51 » 5 2

^ • We then write

I- =

1 . t(1) . t(1) . _J]_ r(2) + T

r ( 3

)

! l,iA 2,i |u| M Vl.i xT i(3) 2

1 2, i

P I (4) + p I (4)

- yI (5)

p

i

1 l,i ^2,1 Y1

i (35)

and define 1^ by substituting throughout the corresponding parameters with

primes:

(3 0

, Bp ...

}

+ {Bq,B^,...}

and till

{a

i ,bi ,ci ,di

} {a

i ,bi ,c

i ,d

i }

. (36)

The integral given for the cross section is then

N

,2i

[I

- hi (37)

20
and our piece-wise analytic expression for the cross section is

d

2 a

dnd|p 2

l

a

3 |p 2

| 2

* e2 l£il

N

,-

r

i

CI

’

( 38 )

Again we stress that the integrations are all performed analytically. The sum

above reflects only the subdivision such that W

2 (x) and W

x (x) can be

approximated with sufficient accuracy by third order polynomials in each

subi nterval

.

Referring to the integrals

, we see that the important

T (l) t(D t (2), t (3) T (3) t

(

4)

n.i* ! 2,i * Ai n.i’ 2,i* U.i*

integrals to be evaluated are

and

m

dx

x i-l

D

l 1/2

( x )

!,m .pJ

m

dx

x i _ i D

2

1/2 (x)

m = -2, -1,0, 1,2 (39)

1 ,m = JPv

m

dx

x i-!

D

l 3/2

( x )

.p- m

dx

x i-l

D

2 3/2

( x )

m = -2, -1,0, 1,2 (40)

21
and

--P

-i dx

X

i

_1

(2A - x) yDj 1/2 (x) D

? 1/2 (x),

(41)

In terms of these integrals we have

(i)

.

1 1, i = a

i 0

oK

l,-2 + (bi 6

o + a

i e

i) K

l,-l

+ < C

i 6

0 + b

i 6

l + a

i 6 2) K

1,0

+ (d

i S

0 + c

1 e

1 b

1 6

2 )K

1>1

+

C<«i®i + c

l®a>K

l.z + d

i B*

3 (42a)

1 2 , i " a

i 6

oK

2,-2

+ * b

i 5

o + a

i {

l* K 2,-l

+ ( Cl «

0 + + a

i «

2 )K2>0

+ (d

i

<S 0 + c^! + b

( 6

2 )K 2jl

+ (d

i 5

i + c

i 42 ) K 2,2

+ d

i 62

K 2,3

(42b)

*

(

i 2) ' a

1 1

1 fc' V

b ’*n lvTr Cit><1 " X i- l)

+ W i (*i

' X ’-l } (42c)

22
( 42d

)

I

I

I

I

I

(3)

l,i

(3)

2,i

(4)

l.i

(4)

2,i

(5)

1

a

i J l,-2

+ b

i J l,-l

a

i J2,-2

+ b

i J2,-l

a

i J l,-l

+ b

i J l,0

a

i J2,-l

+ b

i J 2,0

+ c

i J

1<0 + d.J

1#1

+ c

i J 2.0

+ d

i J 2,l

+ Ci J

lj]L + d- J^2

+ c

i J

2#1 + d

i J2>2

[a

i + 2 Xbi + (2 A) 2c. + ( 2 X) 3d . ] L

t b

i + 2Xci + (2X) 2d

i ](J

li0

- J 2t0 )

[c

i + 2Xdi ](J

lwl

- J 2jl )

(42e)

(42f)

(42g)

(42h)

As just mentioned, the integrals Ip? , Ip? , lj

2

? , Ip- , Ip? , lj

4

? , (4)' (5)'

\\ , and I

v ' follow from the above expressions by substituting throughout

,

i i the corresponding parameters with primes:

{^o»^i****} ' *

* * * f

and

fa

i ,b

i ,c

i ,d

i } > (a

i ,bi ,ci ,di

}

. (36)

Finally, we give the explicit expressions for the parameters which appear

in the integral in terms of the input variables for the analysis, viz., Ip^l,

I

p

2

1 , 6, m, M, and M , as well as the expressions for the integrals J, ,

J. y 1

1

1

23
Jo K i K 0 , and L in terms of these variables. Again we leave the

2, m i,m 2,m

details of the evaluation to the Appendix, giving here only the final results.

It is to be noted that in each case, the arqument of the loqarithm has been

written so that there are no cancellations — all terms in the argument are

positive.

For ease of writing we now write and p 2 for |£ 1

1 and |£ 2

I*

respectively. Then with

e, = /p , 2 + m

2

•*2 = 4> 2

2 + m2

u

=

2 2 M - M

2M

(43)

(44)

we define

X = £^2 - p 1 p 2cos0 - m

2

|ui| = (p L 2 - 2p 1 p 2cos6 + p 2 2

) 1/2

K = £

1 - e2

- W -

ft

2 [X + (e r

e 2-u)(e

1 -e2-|_u|)]

x m e !- c

2

- 1 u.1

2 [X + (e

1 .e2-p)(e

1 -e2+|j u|)]

X M

"

e r

e

2+| u | 1 +

(45)

(46)

24
A, =

Bn =

e

i p 2

’ e

2p i cos 9 \

2 /mp

x si n e \2

p 2 +

M )

+ {—~)

- r

e

i p 2’ e

2p l cos6

'

2p ! [ P 1.P 2

- (eiVm

2

) cos0] p 2 +

~2 y 2

( e

1 P 2

" e2P iCos 0) +

4 [p - y(p 1 e

2

-e

1 p 2cos0) ] 2 + 4(ump 2sin0) 2

M

+ -r- 2m 2Pj 2

( e

1 -e2 )sin 2 0

C, =

x , =

>1

1 A,

A„ =

< (e

2p r

e lP2cos0)\2 /mp

2sin0\ 2

vp i m )

+ v-—

;

Bo = 2p 2[p L p 2

- (e^-m^cos©]

h-

(47)

( e

2P 1 ” e

1 P 2cos 0)"j 2m 2p 2

2

( n

2 0 +

M 1

C. =

+2yjp 1 (£2p 1 -£

1 p 2cos0) +

+ 2,qi 2

) j

4 [p 2 X + y(p 2

£

1 -£2p 1 cos0) ] 2 + 4(ymp1 sin0) 2

M

x

2

"

A,

z, - 2mp

1 p 2 < sin0

n

i = f~A

1

^2 = a"m

2

*2 v x 7 = «

2 UX-x 2 j- + i) 2 1 = m2

(48)

(49a)

(49b)

(49c)

2 - 2B

2x + Cj (50a)

2 - 2B

2

x + C

2 (50b)

25
B

0 = 8m 2 e

1 (e

1 -y)x[p

1 (p 1 -p

2cos6) - (e

1 +e2 )u]

= -2m 2 2Ci( e r y) ^p 1 (p 1 -p 2cos0) - (e

1 +e

2 )(icfy))

+ X^p

1 (p rp 2COS0) - ( e!+ e2 ) vj ^1 + —j5p)

2e,

B 2

= m

2

[p 1 (p 1 -p 2cos0) - (ej+e^ficfyjJfl + -rp)

6

0 = -8m 2 e

2

( e2+p) X[p

2 (p 2

- Plcos0) + (epe^y]

<5

1 = 2m

2 2e

2

( e 2+w)(p 2 (P 2 -Picos0) + (e

1 +e

2

)( ,cfw))

+ x(p

2 (p 2 .p lC os0) + (e

1 +e2 )y)(l - -~j

6 2

= -m 2

[p 2

(

p

2-p pos 0) + ( e!+e

2 ) ( K+y) ] (l - -pp)

T

1 =

"X t*

2 " 4m 2 e

1 e 2

- 2un 2

( e f

e 2-y) - 2 pe

2 x]

t

2 = \ l xZ - 4m 2 e

l e 2

- 2wn 2

( £p£2-y) + 2y£ L x]

Pi = JL

X

2

( X+m 2) (2 £j e 2

- X) + Xk2 + ( X+m 2

) (( X/M) + y)<

+ (m 2/M 2

) X

2 + ^X2 + (

X

2 e

1 /M)

( X+m 2) (2 £^2

” X) + X<2 + ( X+m 2

) (( X/M) + y)<

P-, =

Y

=

X

2 L+ (m 2/M 2

) X

2 + %X2 - (X2 e

2 /M)

~ [( X-Hn 2

) (2 e

1 e 2

- X) + Xk2 + ( X+m 2

) (( X/M) + y)<]

A

n = 1 (51)

26
3

0

' = -8m 4 x[p

1 (p 1 -p 2cos0) - (e1 +e

2 )y]

Sj' = 4m 2 [m 2

(p 1 (p 1 -p

2cose) - ( ex + e

2

) ( y) ) + ^(p x

( P X -P 2cos 9) - (e

1 +e2 )y)]

0

2

' = -2m 2

[p 1 (p rp 2cos0) - ( e

1 +e

2

)( ic+y) ]

V = 8m 4 x[p2 (p 2

- Plcos0) + (e

1+e

2 )y]

= -4m 2 [m 2 (p 2 (p 2-p

1cos0) + ( ( ic+y) ) + *(p 2 (P 2"P i cos e) + (e

1 +e2 )y)]

V = 2m 2

[p 2 (p 2-p

1 cos0) + ( e

1 +e

2 )(icfy) ]

T T = -2 fx 2 - 2m 2

( X + m

2

) 1 2 X

. = -21X 2 - m

4

)

X

2

n' = -2 . (52)

The important integrals are then

°i,o

= s 9

n CA

i x

i - B l)

10

|A

1 xi -8,1*

sgn fA,x , - B, ) £n

/a7 1 ’- 1 1 ^l x i-l

~ B

lJ *

( 53 )

and

27
, /IB.x, - c ,| + nhxJ/c.X

J l,-1

s«"(B

lx

i - C

l) *> (—— )

/

1

B,

x

- , - c. | + D*(x.

, )/C, \

- “i: s 9

n (B

1x

i _ 1 - Cj) In I "

) (54)

/C x i-l e

where

sgn(x) =

1

0

-1

x > 0

x = 0

x < 0

(55)

Then from

q

and

^ we have

J l,l

= t°l(x i)

- D l(

x i-1^ + A^ J l,0 (56)

3B

1,2

= 2A

X

t x

i D

1

C x

i ^ +2A

x

J l,l ‘ 2A

1

J l,0 (57)

1,-2

r °iM D

i

( x i-i }

1L A

i x i-l J °1

Si + C7

J l,-1 (58)

Then the functions

m

are given by

K = UO

(AjX.j-Bj) ^ A

l x i-l"

B

l

^

°l(x i ) Dlfx i-l)

(59)

1,1

1_

(~Bjx

ix

i _ l +C

1

3

0]_ ( x

-j ) ^1 (x

i -1

)

(60)

28
K (61) l,2

= A

x

t J l,0

+ 2B

1 K

1,1

’ C

l K l,oJ

K

- Jl,3 A

x

x

-

i-1

D

l(

x

i _1

)

3B, + "A[ K

l,2

2C

1 Aj

K

l,l (62)

K

l,-1 * t

Dl (*, ) D

l(

x i-l)

+

C,

J

l,-1

+

C,

K

l,0 (63)

1,-2

x

i D

l(

x i)

3B. 2A, +

Cj

K

l,-1

" K

l,0 (64)

The integrals J 0 and K 0 can all be obtained from J, and K 0 merely by

2,m 2,m l,m 2,m

replacing all subscripts 1 by the subscript 2 in the expressions above.

Finally, we have the integral L. In order to write this without

cancellations, we define 7

3 = D^(2 X) = D^(2X) = 2k/x(X + 2m 2

) (65)

5

l,i

= (2XA

1 ‘ B

l )x

i ’ 2XB

1 + C1

h i =

< 2XA

2

" B

2^ x

i " 2XB

2 + C

2

?l,i_l

= (2XA

x - B

1 )x i-1 - 2XB

x + C

1

?2 ,i_i

= ( 2xa

2

- B

2

)

x i-i

- 2xb

2 + C2 (66)

7 See pages 51 and 52 of the appendix, leading to eg. (A. 59).

29
Then

/|C | + BD^fx- )'

- i

s 9

n ( ? i,i-i^

ln

1 ^l.i-1 1 *

- X sgn(t

2j1 ) in

U2 i

I + 6D

2

(x

i

)

1 r \ 0^2 i-1 ^ +

^ x i-l

+ i sgn(E

2>1 . 1

) *n (—‘--1 ‘ L-

- \ ) - s 9n ( 52,i)j tn l x

i ' 2X1

+ \ (

s 9"( el,i-l)

- s 9

n < 52,i-l ))*n|x

i-l

- 2X1 (67)

30
APPENDIX

Algebraic Details

In this appendix we present some of the detailed algebra leading to our

final expression for the cross section, in particular to eqs (27), (28), (38),

and ( 44) - ( 67 )

,

Our starting point is the cross section integrated over the photon

azimuthal angle given here by eqs (22) and (23), As we have noted, this

is the expression given by Mo and Tsai [2], and Tsai [3], and Stein et al

.

[7], apart from the minor modifications and typographical errors which we have

already discussed. Since we perform the remaining integration analytically,

the most appropriate variable of integration is q

2 rather than cose^ in view

of the functions W

2 (-q 2

) and W^-q 2

). Further, as we have already noted, the

form of the integrand itself is then also somewhat simpler. Referring to

eqs (22), (23), and (24), we see that our first task is then to express the

quantities listed in (24), as well as the photon energy, u>, in terms of the

variable of integration, q

2

, and to derive the differential change of variable

given by (25). We see that in (24) the variable cose^ enters only either as u>

or as cocos 0^. (The quantity b enters only as

b

2 = p 2 2 sin 2 0pCo2 sin 2

0|

<

= p 2

2 sin 2 e p

[co2 - (cocose

k

) 2

]

Thus again we have either co or cocosQ^.)

Further, we note that the expressions for cos0

p and cos0

s listed in (24)

follow from

31
£i*“ * £i‘(£i - £2

)

P; 2 - P]P 2C0S6

P]Mcoses ( Pl = Ipjl) (A.l)

and

£2

*U. - £2 *(£j - £2

)

PiP 2cos6 - P 2 2

* P 2 l£|cose

p (p 2 *|£2 |). (A. 2)

In addition, it is useful to note from u Pj - p 2 and the figure below,

that

! 32
and

Thus

p 2 p x |u|

sine “ sine " sine * s p

(A. 4)

p 2sine

p

= PjSine

s (A. 5)

and

lujsine = PjSine (A. 6)

r

|u^|sine

s = p 2sine . (A. 7)

The Jacobian expressed by eq (25) follows directly from eqs (11) and

(15) . From eq (11), viz..

O)

2

u

2 - M

2(u 0 - |uJcose

k

) *

we have

(u 2 - M*

2

) |u|d (cose.)

du> = = —

2(u

0 - h|cose

k

) 2

o)|_u| d(cose

k

) "

(u o - |uj cose

k

) *

33
and from eq (15), viz..

w e

l -

»

we have

du) =

"2M d ( q2 ) *

so that

u)|uj d (cose.) ,

(u

0 - \u\ cose

k

) = 2M d ^ q2 ^

or

2Ma) d (cos0

k

) 1 d(q 2

)

q

4 (u

0 - |_u|cos0

k

) " Tuf ~q^ *

which is eq (25).

Next we express w and ucos©^ in terms of q

2

. In eqs ( 12) -( 15) we have

shown that

(i) =

£

j (A. 8)

in which the quantities e

lt e

2 ,u, and M are defined in (4). Next, from eq (6)

we have (four-vectors are implied in the next two lines)

34
q

2 = (pj - p 2

- k) 2

= (P ! - P 2 ) 2 - 2k*(p

1 - p 2

)

= - 2X - 2^^ - e 2

) + 2u>|ujcos0

k (A. 9)

where, as in (45),

X = e

x e2 - p^.£2

- m

2

.

Thus, using (A. 8),

0 - 0

acos0

k = + ( e

i- e2M e r

e2"*0 + %q 2

Ci + • (a.io)

From (A. 8) and (A. 10) we note that both co and wcos

0^ are linear functions of

q

2

. Thus, referring to (24), we see that a and a' are also linear functions

of q

2

, and that a

2 - b2 and a' 2 - b

2 are quadratic functions of q

2

.

The algebra involved in our calculation is reduced considerably by noting

that upon making the substitutions

subscripts (1 -*• 2, 2 1), y - y, M - M (A. 11)

we have, from (A. 8) and (A. 10),

to - 00

and

wcos 0^ cocos©., ( i . e . , <cos0

k does not change).

35
and from (A.l) and (A. 2)

cosQp - cos 0^ , cos6

s * " cos0p*

Thus when a and a' in (24) are expressed in terms of q

2 as the sole variable,

we may obtain a' from a simply by making the substitutions indicated in (A. 11)

and introducing an overall minus sign. That is, under (A. 11)

a - a' , a

1 -a .

In similar fashion we see from (A. 5) that b

2 as defined in (24) remains

unchanged under the substitutions of (A. 11). Thus under (A. 11) we also have

a

2 - b2

a' 2 - b

2

.

We thus need only write a and a

2 - b

2 in terms of q

2

, and then use ( A.ll)

(A. 13) to obtain a' and a' 2 - b

2

.

We start then with a as given in (24), viz.,

a = e

2

o) - p 2cos6

p wcos9

k

In this expression we then substitute u as given by (A. 8), cosQp as given by

(24) and wcose^ as given by (A. 10). After some tedious but straightforward

algebra we find

a

l [p 1 (p 1 -p 2cos9) - (e

1 +e2 ) M ]x

2 - + ^q 2 [p 1 (p 1 -p

2COS0) - (e

1 +e2

)( ic+p)

]

(A. 14)

36
Then, applying (A. 11) and (A. 12) we have

a = -

[p 2 (p 2-p

1 cos0) + (e

1 +e

2 )y]X

+ ^q 2 [p 2 (p 2-p

1 cos9) + (e1 +e

2 )(icfii)]

(A. 15)

where X and < are given by eg (43) in the body of this report, viz.,

X = £^2 - pjPgCosQ - m

:

K = e

l - e2

- ** ' M

(A. 16)

We next consider a

2 - b2

, which, as we have noted, is a quadratic

function of q

2

. The quantity b is given in (24), viz.,

b = - t»p 2

s i n 6pS i n 6^

and from (A. 6) this may be written as

^lPz . . b = - —|^j|- sin0sin0

k

,

from which

b

2 = (p xp 2sin0)

[to2 (ocos0

k

) 2]

.

We then substitute in this expression u and ocos0

k as given by (A. 8) and

(A. 10), and obtain b

2 as a quadratic function of q2

. With a as given by

(A. 14) we then have then have, after much tedious but straightforward algebra,

the following expression for a

2 - b

2

:

37
a (A. 17) 2 - b

2 = {A, (q

2

) 2 + 2B,q 2 + C,

}

4u 2

1 11

where Ap Bp and are qiven by eqs (47). It should be noted from eqs (47)

(in particular for the integrations that follow) that

A

: > 0 , C1 > 0 . (A. 18)

With

x

= -q 2 > 0

as defined earlier in eq (26), we then write

a

2 - b

2 = — D,(x)

4u 2

1 (A. 19)

where

D^(x) = A^x 2 - 2B^x +

= A

1 [(x - x

x

) 2 + n

x 2

] (A. 20)

with

x

l = (A. 21)

and

38
V (A. 22) = a r . r

2 Vl B

l A

2

1

With Ap B^, and as given by (47) we find, after considerable

A iC i - Bj 2 = (2mp

1 p 2 < sin0) 2 = c

2 > 0 .

Again making the substitutions (A. 11), we obtain

a' 2 - b

2 = (A 9 (q2

) 2 + 2B 9q

2 + C 0 ] 4u 2

2 2 2

with A^, B^, and as given by eqs (48). We now write

where

with

a' 2 - b

2 = D 9 (x)

4u 2

2

D^(x) = A^x 2 - 2B^x +

= A 9 [(x - x

2 ) 2 + n 2

2

]

x

2

= A,

and

39

algebra,

(A. 23)

(A. 24)

(A. 25)

(A. 26)

(A. 27)

I
(A. 28)

From (A. 23) we see that A^C^ - B^ 2

(A. 11). Thus

A

2 C2

- B

2

2 = A

1 C

1

- B

x 2 = c 2

. (A. 29)

Next we consider the terms in (23) with factor v. From eqs (24), (A. 14),

and (A. 15) we have

V? - B?

2

is unchanged under the substitutions

v =

' I

2 Uq 2 + X) [p 2 (p 2 -p

1 cose) + p 1 (p 1 -p

2cose) ]

-2

q

2 + 2 X (A. 30)

since, from eq (8b)

,

\j 2 = Pi

2 _ 2p 1 p 2cos0 + p 2 2

.

Thus the terms with factor v in the integral (22) have an overall factor

1

q

4 (q

2 + 2x)

40
As may be seen by referring to these terms in (23), they also have a

coefficient of the form c

0 + c^ 2 + c

2 (q

2

) 2

. (For the terms multiplying

W^-q 2

) we have c

x = c2 = 0.) Then, to obtain our expression (28), we write

c

o + Cjq 2 + c2 (q

2

) 2

q“(q 2 + 2 X)

c

g - 2XCj + (2X)2

(2X)2 (q

2 + 2X)

1

c

0 - 2Xc

(2X) 2 (A. 31)

Finally we substitute (A. 14), (A. 15), (A. 19), (A. 25), and (A. 30) with

x = - q

2

, in the cross section given by (22) and (23), using (A. 31) for the

terms with factor v. Collecting the terms with similar denominators then

gives the expression (28).

We derive next the expression given in (46) for the limits of

integration, x m

and x^. From (11), viz..

we see that <*> is an increasing function of cose^, and from (15), viz..

it follows that q

2 is also an increasing function of cos0^, and hence that

x = -q 2 is a decreasing function of cose^. Thus the minimum value of x,

viz.

,

o> (ID

q

2 “ = e

l - c

2

- M +

2M • (15)

corresponds to cos0

k = + 1

41
and the maximum value of x, viz..

x = x^ = -q

M

2 corresponds to cose^ = - 1

We now compare eq (15) given just above with eq (A. 10), viz.,

G — G “COS 9. = -At[x + (e

1 -e2 )(£ 1 -c2-u) + *q 2 (l +

When cos

9^

= +1 these two equations give

’m

2

1 e

le 2- y + 2M "

| u

|

[x + ( Eje

2 ) ( ei" e2" y ) + +

m

e , -

1

or

2 [X + (e

1 -e2-|u> |)(er

e

2-|i)] m

1 + WW

as in (46). When cos0

k = -1 these two equations give

e

l’

e

2

-li + 2M?m 7T7T[

x + + Hi^1

or

X M

' 2 [X + ( e r

£ 2+|_u|)( e re

2-u) ]

1 +

e r

e

2+

'J£'

as in (46).

(A. 10)

-)]

42
We next show that the point x = -q 2 = 2A, where, as in (45),

A = £^2 - p x p 2cos 0 - m

2

,

lies within the range of integration:

x < 2A < x M , m

as noted in (29). Now from (A. 9) we see that x = -q 2 = 2 A implies

e

l"

e

2 cos0

k

=

|ui|

* (A. 32)

Now as we show immediately below,

£ - G 0 < -ijy2 < 1 . (A. 33)

Thus since the integration region, when we choose cos0

k as integration

variable, is the entire range of possible values, - 1 < cos0

k < 1, we see from

(A. 32) and (A. 33) that the point x = 2A corresponds to a point in the range of

integration with x as integration variable (recall that x is a uniformly

decreasing function of cos0

k

), i.e., that

x < 2A < x M . m

Further, since the photon energy, given by (15), must be positive for cose

k within the range of integration, i.e..

43
a) = e, - e - y + > 0 , 1 2

- - 2M

it follows that this is also true for the point - q

2 = 2X just considered.

Thus

<=ei" e2" w "M >0 * (A. 34)

Returning now to the proof of (A. 33), we note from e

2 < e

1 (and p 2 < p^ that

: l ’ p l =

e, +m

2

„ m

2

l + Pi e 2 + P 2

= e„ - 2

- P 2

Thus

0 < e r

e2 < p rp 2

< |p rp 2l

from which (A. 33) follows.

Finally, we conclude our presentation of algebraic details with a

consideration of the expressions given in eqs (53) - ( 67) for the integrals

defined in eqs ( 39) -

( 41 ) . We begin with g, defined by eq (39), viz..

wf dx

x M D

1 l/2(x)

(A. 35)

where from eq (50a)

D^(x) = A^x 2 - 2BjX + C

44
A straighforward integration gives [9]

_L £n /

A

l x i-1

" B

1 * D1 (x i-l

)v^l

/A.

(A. 36)

where c

0 is an arbitrary constant. Now for reasons that will be apparent

momentarily, we choose

c o = AjCj - B

^ 2

= ? = 2mp 1 p 2ic sin0 > 0 (A. 37)

from (A. 23) and (A. 34). The arguments of the logarithms in (A. 36) are then

dimensionless and, from eqs. (47) and (A. 23), are of order e/m if A^x - > 0

(where x = x. or x-_j). 0n the other hand, since

D

1 (x)A1 - (Axx - Bx

) 2 = (A

xx

2 - 2B

xx + C^Aj - (AjX - B^ 2

= A

1 C1 - Bx 2 (A. 38)

we have

A

x x - B1 + D

1 ^(x)/A

1 -(Ajx - B

x

) + D

1

?5 (x)/A

1

(A. 39)

45
Thus if A^x - < 0 we can write the arguments of the logarithms in (A. 36) in

the form given by the right hand side of (A. 39), from which we see that the

argument is then of order m/e. Thus if A^x - < 0 there is a cancellation,

leaving terms of relative order (m/e) 2

. We will therefore write the argument

so that there is no cancellation, i.e., as given by the left hand side of

(A. 39) if AjX - Bj > 0 and by the right hand side of (A. 39) if A^x - B^ < 0.

These can be combined into a single expression by writing

|A

l xi - Bi! *

1— sgn(A

1x

i . 1

- Bj) An

1

' i-l ; 1 /A

1

where

1 x > 0

sgn(x) = 0 x = 0 .

-1 x < 0

This is the result given in (53).

The integral

^ is also defined by (39), viz..

(A. 40)

Again, a straightforward integration gives [9]

46
+ — ln

r

B

i x i-i + c

i + Di^ x j-i)

/c

i /c

(A. 41)

1 x i-l

c

o'

where c

0

' is an arbitrary constant. In analogy with (A. 38) we note that

D

1 (x)C1 - (-B lX + C^ 2 = (AjX 2 - 2Bxx + - (-B^ + C^ 2

= x

2 (A

1 C1 - Bx 2

)

= x

2

?

2

. (A. 42)

Thus we again choose c

0

' = c, and note that

-B^x + C-^ + D-j^(x) xc

x C -(-BjX - Cj) + DjVj/Cj

(A. 43)

Thus if -B^x + Cj < 0 we can write the arguments of the logarithms in (A. 41)

in the form given by the right hand side of (A. 43), thus avoiding the cancellation. As before, we now write

^ in a form which avoids the possible

cancellation in the arguments of the logarithms, viz..
1 /l-B.x, + C,| + D.^xJ/C,

J i-i

'

' *= ‘•"'•V* * c .’ ” (

1 , .

1 , + C. | + D, ^(x

. ,)/C,

+ ^

*9nC-8

1*

1 . 1

C

x

) In- ^ 1 ^

x m'

which is equivalent to the result given in (54).

The integrals

g

and defined by (40) and given expl icitely by

(59) and (60), are elementary and present no particular difficulty. The

expression of

m

and min terms of

g, o*

and K

1 1* 9^ ven

( 56 ) - ( 64) , is also elementary and may be found, for example, in the integral

tables of Grffbner and Hofreiter [9].

Lastly, we have the integral L, defined in (41), viz.,

L =

1

(2X - x) .D 1 l/2(x) D

2

1 / 2 (x)

dx (A. 44)

We first note that

D

1 (2 X) = D

2

(2 X) , (A. 45)

and thus the integrand has in fact no singularity at x = 2X. This may be seen

most easily as follows: From (A. 30)

a' - a

= -^(q 2 + 2X)

= Mx - 2 X) (A. 46)

48
and from (A. 19) and (A. 25)

Thus

and

We define 8

D

x (x) - D

2

(x) = 4u_2 (a 2 - a' 2

)

.

D

x (x) - D

2

(x) = -2jj 2 (a + a')(x - 2 \)

D

x (2X) = D2

(2X)

3 = D

1 1/2 (2X) = D

2

1/2 (2X) .

In analogy with the procedure followed for

q

and

^ we first

integration straightforwardly (see [9]) and write

L

= i In N.i * BD

l tx

i

] 6 *" \

|x r

2X| c

1 tn

/

sl,j-l

+ gD itx j-i)

» \

|x i-l _2X| 5

,

1 >n

S \|x

i -2X|c

,

8 See the end of this appendix for a derivation of the value of 3.

(A. 47)

(A. 48)

(A. 49)

perform the

(A. 50)

49
where

^(x) = (2 XA

1 - B

t )x - 2XB

x + Cl

C

2 (x) = (2 AA

2

- B

2 )x - 2XB2 + C2 (A. 51)

5l,i ~ Mx

i^ * 5l,i-l

= Mx i-1^

^2,i

= Vx

i^ ^2,i-l

= ^2 ^ x

i - 1

^ (A. 52)

Then noting that

^(x) + BD^(x )

I x -2 X | ?

|x-2A|c -^(x) + SO^(x)

(A. 53)

(and a similar equation in which all subscripts 1 are replaced by the

subscript 2), we start by writing L in a form similar to that just used for

q

and so that there is no cancellation when any of the terms ^

52,i- or 52,i-l is native:

1 ( \ 1

1 ^l.i I + ( x

j

) L

= 6 sgn

(^l,i)

ln

\

| x 1 -2 X | c

1 /ISi i J + eofy. ,)

- 1

*n

(

' ix^^-zxu

'

" B sgn (^2 ,i ) *n \ | x - -2 X | c j

i /1 ,

1

1 + eoJ(x. , ) + 1 ^2,1.1)

n ( tT-pFirr^- (A. 54)

50
However, in this form of the expression for L there can be a large

cancellation for or x-

^ sufficiently close to 2A, since, as we have noted,

the integrand in the integral for L, (A. 44), has no singularity for x = 2A.

We therefore write L in a form such that this cancellation is manifest.

Equation (A. 54) may simply be written in the form

and a similar expression with i replaced by i-1. This is the result given in

eq (67). Now from (A. 51) and eqs (50a) and (50b) we have

(A. 55)

where

(A. 56)

5

: (2X) = Dx (2X)

S

2 (2\) = D2

(2A) (A. 57)

and hence from (A. 48) and eqs (50a) and (50b)

51
^(2X) = ^( 2 \) > 0 (A. 58)

Therefore, for x. sufficiently close to 2X we have sgn( 5 . .) = +1,

t

* sqn(£? ,) * +1 and thus L. in (A. 56) can be written more simply as

Here (i.e., when sgn( 5-, J = sgn(C7

• )) there is clearly no singularity for

x- = 2X. These observations can be used to avoid any possible problem arising

in a computer program from x^ or x-_^ being chosen too close to 2X. We may

simply define L by (A. 55 ), with L- (and L-^) given by (A. 59 ) if both ^

. > 0

and . > 0, and otherwise given by (A. 56 ). (It can be shown that we can not

have both < 0 and C? 1 < 0.)

Finally, we require the value of 6, defined in (A. 49), given in the body

of the report in eq (65), and appearing in the equations for L, (A. 50) and

(A.54)-(A.59) . From the definition, (A. 49),

it is clear that one can evaluate 6 simply by substituting x = 2X in the

equations defining D^(x) and D^x), (50a), and (50b). This is a rather

cumbersome procedure, however. We therefore give the following simpler

derivation:

For x = 2X we have

(A. 59)

6 = dJ(2A) = d|{2X) ,

a a

D

1 (2X) = D2

(2X) ,

52
from (A. 46) and (A. 47), and

(D = <

from eqs (15) and (A. 16). Further, also for x = 2X, we have

cos0

k =

from (A. 32), from which

sin6

k

/u 2 -

'( e r

e2 ) a

lul

2X

'

111

(A. 60)

Substituting these values for w, cose^, and sine^ in the expressions (24) for

a' and b we have, for x = 2X,

a

' = ^ I e

iii 2 - Pi( e1 - e

2 )(PrP2cos9) ] (A. 61)

u 2

and

<p ,p 2 si n0 /2X

b = - . (A. 62)

u

2

In writing (A. 61) we have used the expression for cos6

s given in (24), and in

writing (A. 62) we have used the expression for sin0

p given in (A. 6).

Recalling that u

2 = p x 2 - 2p 1 p 2cos6 + p 2 2

, a little algebra then gives

53
a' = — (e

1 +e2 )(e1 e2

- p^cosQ - m

2

)

<X(ei+e2

) (A. 63)

Substituting (A. 63) and (A. 62) in (A. 25) we then have

D

2 (2 X) = [X(e

1 +e2 ) 2 - 2p 1 2 p 2 2

si

n

2 e] . ^ 1

1

2

Again recalling that X = - p^cosQ - m

2

, a little algebra gives

M e

i +e2 ) 2 * 2p 1 2p 2

2sin 2 0 = (e

x e

2

- Pl p 2 + m

2 )^2 - 2p 1 p 2cos6 + p 2 2

)

so that

= ( X + 2m 2

) u

2

,

D

2

( 2 X) = 4

k

2 X( X + 2m 2

) (A. 64)

and

e = D^(2 X) = D^(2X) = 2</x( X + 2m 2

)

. (A. 65)

54
REFERENCES

[1] Y. S. Tsai, in Proceedings of the International Conference on Nuclear

Structure, 1963 edited by Hofstadter and Schiff (Stanford Umv. Press,

Stanford, California, 1964), p. 221.

[2] L. W. Mo and Y. S. Tsai, Rev. Mod. Phys. 41, 205 (1969).

[3] Y. S. Tsai, SLAC-PUB 848, January 1971.

[4] G. Miller, et al. Phys. Rev. D5, 528 (1972).

[5] L. C. Maximon and D. B. Isabelle, Phys. Rev 133 , B1344 (1964).

[6] J. D. Bjorken and S. D. Drell, Relativistic Quantum Mechanics (McGraw

Hill Book Co., New York, 1964), see Chap. 7.

[7] Stein et al., Phys. Rev. £12, 1884 (1975). In eq (A24) of this

reference, the authors have corrected two misprints that appear in

eq (A. 24) of SLAC-PUB 848, reference 3 above.

[8] H. Nguyen-Ngoc and J. P. Perez-y-Jorba, Phys. Rev. 136 , B1036 (1964).

[9] W. GrO’bner and N. Hofreiter, Integraltafel , Part I (Springer, Vienna,

1961).

55
The form factors used in these timing runs were derived from a DWBA

calculation for 12C at 500 MeV (see subroutine PWAFFI in the computer

program listing, pp. 98-100).

PWA denotes the piece-wise analytic evaluation of the integral as presented in this report.

All times listed are in seconds. Numerical integrals were

calculated with an accuracy of 0.0001 H.

p ^250 MeV

1 deg

.

P 2 (MeV) Simpson Gauss <n=16) Kronrod (n:61) PWA

10 21.07 5.42 1.57 . 43

125 62.03 16. 12 2.36 .22

240 3.15 1.22 . 05 .02

5 deg

.

P 2 ( MeV ) Simpson Gauss (n:16) Kronrod <n=61) PWA

10 11.66 2.06 .96 .42

125 4.37 1.30 .96 .22

240 1.86 .59 .61 .03

45 deg.

P 2 (MeV) Simpson Gauss <n=16) Kronrod (n=61) PWA

10 2.46 .84 . 87 .43

125 2.33 1.00 .78 .23

240 2.21 .07 .80 .02

90 deg.

P 2 (MeV) Simpson Gauss (n=16) Kronrod (n=61) PWA

10 3.01 1.24 . 93 .44

125 2.61 1.17 .84 .23

240 2.46 1.23 .09 .02

135 deg.

P 2 (MeV) Simpson Gauss (n:16) Kronrod <n=61) PWA

10 5.35 1.78 1.17 .43

125 3.15 1.37 .92 .23

240 3.30 1.31 .97 .01

175 deg.

P 2 (MeV) Simpson Gauss (n=16) Kronrod (n=61) PWA

10 11.71 2.43 1.06 . 44

125 4.47 1.73 1 . 18 .23

56
All times are in seconds. Numerical integrals uiere

calculated with an accuracy of 0.0001%.

Pl =500 MeV

1 deg

p 2 (MeV) Simpson Gauss (n=16) Kronrod (n:61) PWA

10 48.45 8.53 2.39 .99

250 70.24 16.29 2.90 . 44

490 3.66 1.20 1.01 .03

5 deg.

P 2 (MeV) Simpson Gauss <n=16) Kronrod (n=61) PWA

10 33. 19 3.06 2.21 .99

250 6.57 2.16 1.32 . 44

490 3.35 1.24 1.04 .03

45 deg.

P 2 (MeV) Simpson Gauss (n=16) Kronrod <n=61) PWA

10 3.26 1.31 .75 .98

250 3.76 1.64 1.33 .46

490 5.01 2.32 1.96 . 02

90 deg.

P 2 (MeV) Simpson Gauss (n=16) Kronrod (n=61) PWA

10 6.40 2.33 1.29 .99

250 4.31 1.75 1.55 .48

135 deg.

P 2 (MeV) Simpson Gauss (n=16) Kronrod (n=61) PWA

10 20.92 3.61 1.97 .99

250 6.40 2.58 1.44 .53

175 deg.

P 2 (MeV) Simpson Gauss (n=16) Kronrod (n=61) PWA

10 23.49 5.90 1.86 1.00

250 8.70 3.46 1.89 .56

57
All times listed ere in seconds. Numerical integrals were

calculated with an accuracy of 0.0001*.

Pj-750 MeV

1 deg

.

P 2 (MeV) Simpson Gauss <n=16) Kronrod <n=61) PWA

10 93.26 16.85 2.52 1.63

375 46.09 18.82 2.96 .70

710 5.72 1.86 1.23 .08

5 deg.

P 2 (MeV) Simpson Gauss ( n : 16

)

Kronrod <n:61) PWA

10 59.96 15.38 3.94 1.63

375 7.26 2.33 1.36 .70

710 4.49 1.89 1.27 .08

45 deg.

P 2 (MeV) . Simpson Gauss (n=16) Kronrod <ns61) PWA

10 5.66 2.33 1.02 1.62

375 5.00 2.09 1.75 .75

710 5.60 2.57 1.30 .07

90 deg.

P 2 (MeV) Simpson Gauss (n=16) Kronrod (n=61) PWA

10 11.47 4.40 1.92 1.62

375 6.68 3.44 1.94 1.21

135 deg.

P 2 (MeV) Simpson Gauss (n:16) Kronrod (n=61) PWA

10 26. 18 7.16 2.36 1.62

375 9.84 2.39 2.14 .98

175 deg.

P 2 (MeV) S i mpson Gauss (n=16) Kronrod <n=61) PWA

10 57. 47 15. 11 3.70 1.62

375 10.22 3.56 1.47 .88

58
COMPUTER PROGRAM LISTING

In the pages which follow we give the listing of the computer program.

However, the subroutines PROMPT, TOD, and CLK (which are called in the main

program MSPLTIM, as well as in the timing subroutines STARTT and STOPT) are

system information routines which are unique to the Perkin-Elmer system on

which the program was run. They are not, therefore, included in the program

listing.
I

1 STITL MSPLTIM — timing for Tsai and PWA

2

' C This program compares the formula for the radiative tail

3 C given by Tsai and the Piece-Wise Analytic method ( PWA

)

4 C

5 IMPLICIT DOUBLE PRECISION (A-H,0-Z)

6 c Maximum number of angles • 7 PARAMETER (MNTslO)

e c Maximum number of momenta

9 PARAMETER (MNP=10)

10 LOGICAL OPENED, FAST

n DOUBLE PRECISION TIM( 4 ) , THETA( MNT ) , PF IN ( MNP

)

12 INTEGER ICLKBC4)

13 CHARACTER FILE*40, T IME*8 , DATE*8 , YESNO*SO

14 COMMON /EUICOM/ NSTEP,EPS

15 COMMON /'CONST / EM, EM2, AM, HBARC, ALPH, PIE, ALPHA, RTD , DTR

16 COMMON /KIN/ EEX, TH, AMO, S3, P3, CTH, STH, S32, ES, P32, EP, AMF,SDP3, SDP,

17 + U32, U3, UO, U2, CTHS, CTHP, STHP, XMU

18 COMMON /CORREC/IWTILD, WTIL1, WTIL2, IS0FT,S0FT1

19 COMMON /TIME/ T I M, I CLKB

20 DATA FAST/. FALSE./

21 CALL TOD ( T IME , DATE

)

22 CALL INICON

23 INQUIRE <1,0PENED=0PENED)

24 IF ( .NOT. OPENED) OPEN ( 1 , FILEs * CON : *

)

25 2 IF ( .NOT. OPENED) THEN

26 WRITE (1,901) ' ', TIME, DATE

27 901 FORMAT (A, 'MSPLTIM version 1A ',A,1X,A,

28 + ' (last edit 09/06/83)')

29 CALL PROMPT( 1, 'Output file')

30 END I

F

31 READ (1,80) FILE

32 IF ( FILE . EQ . ' ') STOP. 'All done'

33 00 FORMAT (A)

34 OPEN (6, FILErFILE, ACCESS =' SEQUENT I AL ' ) 35 IF ( .NOT. OPENED) CALL PROMPT ( 1 ,' Exc i t at i on energy ' ) 36 READ (1,*) EEX

37 EEXrEEX/HBARC

38 IF ( .NOT. OPENED) CALL PROMPT ( 1 ,' Inc i dent momentum ' ) 39 READ (1,*) S3

40 S3: S3/HBARC

41 IF ( .NOT. OPENED) CALL PR0MPT(1,'A of target')

42 READ (1,*) AMO

43 „• AMO: AMO*AM

44 IF ( .NOT. OPENED) WRITE (1,910)

45 910 FORMAT (' Enter values for final momentum (end with negative):')

46 NP : 0

47 99 READ (1,*) P3

48 IF (P3. GE. 0. ) THEN

49 NP : NP+

1

50 PFIN(NP) :P3/HBARC

51 IF (NP.EQ.MNP) GO TO 88

52 GO TO 99

53 END I F 54 IF (NP.EQ.O) STOP

55 88 IF ( .NOT. OPENED) WRITE (1,909)

56 909 FORMAT (' Enter values for THETA (end with negative):')

57 NT : 0

60
5B

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

103

109

110

1

1 112

113

,114

.9 READ (1,*) TH

IF (TH. GE. 0. ) THEN

NTrNT+1

THETA ( NT )=TH*DTR

IF (NT.EQ.HNT) GO TO 9

GO TO 8

ENDIF

IF (NT.EO.O) STOP

9 IF ( .NOT. OPENED) CALL PROMPT ( 1

,

+ 'Starting number of steps for the integral')

READ (1,*) NSTEP

IF ( .NOT. OPENED) CALL PROMPTt 1, 'Accur acy of the integral')

READ (1,*) ACC

IF ( .NOT. OPENED) CALL PROMPT ( 1

, 'W-TILDA correction')

READ (1/80) YESNO

IF ( INDX ( YESNO/ ' Y ' ) .GT.O) THEN

IWTILD=1

ELSE

IWTILDrO

ENDIF

ISOFTrO

CALL PWAFFI

DO 502 1=1, NT

TH= THETA ( I

)

WRITE (6/901) ' 1 TIME/ DATE

WRITE (6/902) TH*RTD, S3*HBARC

,

+ AMO/AM/ EEX*HBARC/ ACC, YESNO

902 FORMAT ('0 Scattering angle =',F7.2,' degrees'/

+ T62/ ' Inc ident momentum =',F8.3, ' MeO/c'/

+ ' Mass of nucleus ='/F8.3/' AMU'/

+ T61 ,

' Exc i t at i on energy = ',F8.3,' MeU'/

+ ' Accuracy of numeric integral =',1PE8.1,

+ T6l/'W-Ti)da correction made? '/A5)

WRITE (6/903)

903 FORMAT ( ' 0 ' , T6 . 'F i n a 1 ' / T17, ' Ts a

i

' / T26/ ' For m factor'/

+ T39, ' Integr al ' , T52/ 'Total '/T65, 'PWA',T76/

+ 'Search '/ T37/ 'Integral '»T100/ 'Total ',T110/ '* Diff'

+ T126/ 'Speed '

/

+ T3/'momentum'/T15/'result'/T27/ 'time'/ T39/ 'time'/

+ T51 / 'time', T63, ' resu 1 t

' , T75, 'time', T87, 'time',

+ T99, ' t i me ' , T124, ' incre ase ' ) CALL CLK(ICLK,Y'2B', .TRUE.

)

DO 503 J=1,NP

P3=PF IN ( J

)

C Initialize some commons

CALL INIKIN

T I M ( 1 )=0. DO

TIM(2)=0. DO

TIM(3)=0. DO

TIM(4)=0.D0

IF ( J . LE . 1 . OR . . NOT . FAST ) THEN

EPS=-ACC

ELSE

EPS=ACC*SIGTS

ENDIF

CALL ST ARTT ( 2

)

CALL TSAINT (SIGTS)

61
115 CALL STOPT ( 2

)

116 • CALL ST ARTT ( 4

)

117 CALL PWASPL(SIGPWA)

118 CALL STOPT ( 4

)

119 WRITE (6/905) P3*HBARC *SIGTS»TIM(1)»TIM(2>-TIM(1)»

120 TIM(2)»SIGPWA»TIM(3)*TIM(4)-TIM(3)»TIM(4)>

121 + 200. *<SIGPWA-SIGTS>/'<SIGPWA+ SIGTS) , T I M < 2 ) TIM ( 4

)

122 905 FORMAT (1X#F8.2*T12» 1 PE 1 2 . 5 » T24 , 0PF9 . 3 , T36, F9. 3,T4e,

123 F9. 3,T60, 1PE12. 5/T72, 0PF9.3, T84,F9. 3, T96/F9.3, T108,

124 + 1PE10. 3/ T 1 20 » 0PF10. 2)

125 503 CONTINUE

126 502 CONTINUE

127 CLOSE (6)

128 GO TO 2

129 END

62
12

3

4

5

6

7

8

9

• $TI TL STARTT -- Start a tune interval

SUBROUTINE STARTT(I)

IMPLICIT DOUBLE PRECISION <A-H,0-Z)

DOUBLE PRECISION T I M f 4

)

INTEGER I CLKB ( 4

)

COMMON . /TIME/ TIM, ICLKB

CALL CLK(ICLKB(I),Y'2B', .FALSE.

)

RETURN

END

63
12

3

4

5

6

7

8

9

10

11

ST1TL STOPT — Stop a time interval

SUBROUTINE STOPT(I)

C Note that the constant is for optimized Fortran.

IMPLICIT DOUBLE PRECISION <A-H,0-Z)

DOUBLE PRECISION TIM<4)

INTEGER ICLKB ( 4

)

COMMON /TIME/ TIM, ICLKB

CALL CLKdCLKE, Y'20', .FALSE. ) TIMC I)=TIM(I )+( ICLKE-ICLKBC I ))/ 1000 0009374

RETURN

END

64
12

3

45

6

7

8

9

10

11

12

13

14

15

16

17

18

STITL INICON — Initialize constants

SUBROUTINE INICON

C Subroutine to initialize the common CONST. All energies in

C units of fm-l. Constants from RMP 52* No. 2 (April 1980) as

C reprinted by the CERN particle properties data booklet.

IMPLICIT DOUBLE PRECISION (A-H*0-Z>

COMMON /CONST/ EM* EM2* AM* HBARC* ALPH* PIE* ALPHA, RTD* DTR

HBARC: 197. 32858D0

ALPH=137. 03604D0

AMr 931 . 5016D0/HBARC

EM= . 51 10034D0/HBARC

EM2:EM**2

PIE=3. 1415926535897D0

ALPHA= 1 .ODO/ALPH

DTR= PIE/ 180. DO

RTD:180. DO/PIE

RETURN

END

65
12

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

no

on

$TITL INIKIN — Initialize kinematic variables

SUBROUTINE INIKIN

C Subroutine to initialize the commons KIN and CORREC.

C In general this must be called if any of the first

C five elements of KIN change.

IMPLICIT DOUBLE PRECISION (A-H/O-Z)

COMMON /CONST/' EM, EM2, AM, HBARC, ALPH, PIE/ ALPHA, RTD, DTR

COMMON /KIN/ EEX, TH, AMO, S3, P3, CTH, STH, S32, ES, P32 , EP , AMF , SDP3 , SDP.

+ U32,U3,U0,U2,CTHS,CTHP, STHP,XMU

COMMON /CORREC/IWTILD, WTIL1, WTIL2, IS0FT,S0FT1

C sin's and cosin's of angles

CTHrDCOS(TH)

STH=DSIN(TH)

C Incident momentum squared

S32=S3*S3

C Incident energy of the electron

ES=DSQRT ( S32+EM2

)

C Final momentum squared

P32=P3*P3

C Final energy of electron

EPrDSQRT ( P32+EM2

)

C Final rest mass of the nucleus

AMF = AMO+EEX

C This is ( AMF**2 - AM0**2 ) /AMO

XMU=EEX+. 5D0*EEX**2/AM0

C 3-vector dot product of s and p SDP3=S3*P3*CTH

C 4-vector dot product of s and p SDP =ES*EP-SDP3

C Momentum part of u=s+t-p (momentum loss of the electron)

U32=S32+P32-2.D0*SDP3

U3=DSQRT(U32)

C Energy part of u=s+t-p

UO=ES+AMO-EP

C Length of u=s+t-p (squared) (missing mass squared)

U2=U0**2-U32

C cosin of angle between s and u CTHS=(S3-P3*CTH)/U3

C cosin of angle between p and u CTHP=(S3*CTH-P3 ) /U3

C sin of angle between p and u STHP=S3*STH/U3

Initialize coefficients for change of W to W-TILDA.

WTIL1=1 . D0+28.D0*ALPHA/(9.D0*PIE)-ALPHA*DL0G(ES/EP)**2/(2. DO*PIE)+

+ ALPHA* (PIE*P IE/6. DO-SPENCE( . 5D0*( 1 . DO+CTH) ) )/PIE

WTIL2=13. D0*ALPHA/(PIE*6. DO)

Initialize coefficients for soft photon correction.

WS = ES-EP/( 1 . DO— EP* ( 1 . DO-C TH ) /AMO

)

WPr ES/ ( 1 . DO+ES* ( 1 . DO-CTH ) /AMO ) -EP

SOFTlrWS*WP/(ES*(EP+WP)

)

RETURN

END

66
12

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

4 8

49

50

51

STITL TSAINT — Evaluate Tsai by doing numerical integral

SUBROUTINE TSAINT ( UAL

)

C The formula given by Tsai in Slac-pub 848 (as we 1 1 as in

C Mo and Tsai RMP 41, 205 and Stein et . al. Phys Rev D 12,1884)

C is evaluated by doing the integral numerically. This

C is done in up to 7 separate intervals. Cal.l INIKIN before

C calling this routine.

IMPLICIT DOUBLE PRECISION CA-H,0-Z>

EXTERNAL TSAI

LOGICAL DEBUG

INTEGER NS ( 7

)

DOUBLE PRECISION EDGE ( 8

)

COMMON /CONST/ EM, EM2, AM, HBARC, ALPH, PIE, ALPHA, RTD, DTR

COMMON /KIN/ EEX , TH, AMO, S3, P3, CTH, STH, S32, ES, P32, EP, AMF, SDP3, SDP,

+ U32,U3,U0,U2,CTHS, CTHP, STHP,XMU

DATA DEBUG/. FALSE.

/

DATA EDGE C 1 ) /-I . DO/, EDGE C 8 ) /I . DO/

UAL=0. DO

THL IM=AM0*CES-EP)/C2. DO*ES*EP)

IF (THLIM.LE. 1. DO) THEN

THLIM=2. DO*DASIN( DSQRT(THLIM)

)

IF (TH.GT.THLIM) RETURN

ENDIF

THS = DATAN2 C DSQRT C 1 . D0-CTHS**2) , CTHS)

DTHS= DSQRT CEM/ES)*2.0D0

THP=DATAN2 (STHP, CTHP)

DTHP=DSQRTCEM/EP)*2.0D0

EDGE (2) =DC0S(DMIN1 (PIE, THP+DTHP)

)

EDGE ( 3 ) =CTHP

EDGE (4) rDCOSCDMAXl (0. DO, THP-DTHP )

)

EDGE (5) = DCOS C DMIN1 CPIE, THS + DTHS) > EDGE ( 6 ) = CTHS

EDGE (7 ) = DCOS C DMAX1 (0. DO, THS-DTHS)

)

IF CEDGEC4).GT.EDGEC5) ) THEN

EDGE C 4 ) = CEDGEC3)+EDGE(6) )/2. DO

EDGEC5)=EDGEC4)

ENDIF

DO 504 1=1,7

NSC I ) =0

IF (EDGE(I).LT.EDGECI+1) ) THEN

DSIGrEUALI ( EDGE ( I ) , EDGE ( 1 + 1 ) , TSAI, NSC I )

)

IF C DEBUG ) WRITE Cl,*) CH AR C I CHAR C ' A ' ) +1 - 1 )

,

+ EDGE C I ) , EDGE C 1 + 1 ) , DSIG

UAL=UAL+DSIG

ENDIF

504 CONTINUE

UAL = UAL*ALPHA**3*P32*AM0/ C PIE*EP*S3)

IF C DEBUG ) WRITE Cl, 31) NS, UAL

31 FORMAT C IX, 715, E10. 3)

RETURN

END

67
12

3

4

5

6

7e

9

10

i l 12

13

14

15

16

17

16

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

nn

STITL TSAI -- Evaluate the integrand of Tsai's formula

SUBROUTINE TSAI < CTHK, SIG

)

C Subroutine to calculate the integrand of the formula for

C the brems. cross section given by Tsai in SLAC-PUB-B49 (Jan

C 1971) (Formula A. 24). The constants in the commons CONST

C and KIN must be set. Angles- in radians# energies and

C momenta in fm-1

IMPLICIT DOUBLE PRECISION (A-H,0-Z)

COMMON /CONST/ EM# EM2# AM# HBARC# ALPH# PIE# ALPHA, RTD# DTR

COMMON /KIN/ EEX# TH, AMO# S3# P3# CTH# STH# S32# ES# P32# EP# AMF, SDP3. SDP#

+ U32#U3#U0# U2# CTHS# CTHP# STHP# XMU

COMMON /CORREC/IWTILD# WTIL1# WTIL2# ISOFT# S0FT1

C sin of photon angle

STHK rDSQRT ( 1 . D0-CTHK**2)

Energy of the real photon (avoid the cancelation of mass of target

squared

)

W = ( EM2-SDP+AM0* ( ES-EP-XMU ) ) / ( U0-U3*CTHK

)

C Momentum part of s-p-k

Q32=U32+W**2-2. D0*U3*W*CTHK

Q3=DSQRT(Q32)

C Length of q=s-p-k (4-momentum transfer to nucleus)

Q2: (ES-EP-W)**2-Q32

C

A=W*(EP-P3*CTHP*CTHK)

AA=W*(ES-S3*CTHS*CTHK

)

B2= ( W*P3*STHP*STHK )**2

X=DSQRT(A**2-B2)

YrDSQRT (AA**2-B2)

C The constant NU*(l/x-l/y) of Tsai

XYNU=(AA+A)/(X*Y*(X+Y)

)

C Calculate the integrand

SIG=W*( F0RM2C Q2)*C -A*EM2*(2. DO*ES* ( EP + W ) + . 5D0*Q2 ) /X**3 -

+ AA*EM2*(2. D0*EP*(ES-W)+. 5D0*Q2 )/Y#*3 - 2. DO +

+ 2. D0*XYNU*( EM2* (SDP-W*>f<2 ) + SDP* ( 2 . D0*ES*EP-SDP + U* ( ES-EP ) ) ) + + (2. D0*(ES*EP+ES*W+EP**2)+. 5D0*Q2-SDP-EM2)/X - + (2. D0*(ES*EP-EP*W+ES**2)+.5D0*Q2-SDP-EM2)/Y) +

+ F0RM1 (Q2)*( ( A/X**3+AA/Y**3 ) *EM2* ( 2. D0*EM2+02) + 4 . DO +

+ 4 . DO*XYNU*SDP* ( SDP — 2 . D0*EM2 ) +

+ ( 1 . DO/X-1 . D0/Y)*(2. D0*SDP+2. D0*EM2-Q2) ) ) / ( 02*Q2* ( U0-U3*CTHK )

)

IF (ISOFT. EQ.O) RETURN

SIG=SIG*S0FT1**( ALPHA* ( DLOG ( -Q2/EM2 ) -1 . DO) /PIE)

RETURN

END

68
12

3

45

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

•$T I TL PWASPL — Piece-wise analytic integration of Tsai

SUBROUTINE PWASPL(UAL)

C We assume that the COMMON KIN has been setup already

C by a call to INIKIN and that the common PWAFF has been

C initialized with a call to PWAFFI.

‘ C Size of arrays in PWAFF

PARAMETER (MNPT=200)

IMPLICIT DOUBLE PRECISION (A-H/O-Z)

DOUBLE PRECISION X ( MNPT ) , A ( MNPT ) , B ( MNPT ) , C ( MNPT ) , D ( MNPT

)

DOUBLE PRECISION AP ( MNPT ), BP ( MNPT ), CP ( MNPT ), DP ( MNPT

)

COMMON /PWAFF / NPT , X, A , B , C , D , AP , BP , CP/ DP

COMMON /KIN/ EEX,TH,AM0,P1,P2, CTH, STH,P12,El,P22, E2, AMF/SDP3, SDP,

+ U32# U3»U0/U2»CTHS/CTHP / STHP, XMU

COMMON /CONST/ EM, EM2> AM/ HBARC, ALPH, PIE/ ALPHA, RTD, DTR

UAL=0. DO

CALL PWAINICXl/XMAX)

CALL ST ARTT ( 3

)

DO 30 1=1, NPT

IF (Xl.LT.XCI) ) GO TO 40

30 CONTINUE

I =NPT

40 CALL ST0PTC3)

50 X2 = DMIN1 (X(I ) , XMAX)

CALL PWAEUACX1, X2, ACI ) , B( I ) , C< I >, D( I )

,

+ AP< I ) , BP( IJ , CP( I) , DP( I ) , XINT)

UAL =UAL + X I NT

1 = 1 + 1 XI =X2

IF (X2.LT.XMAX) GO TO 50

UAL = UAL*ALPHA*>K3*P22/( PIE*E2*P1 ) RETURN

END

69
12

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

non

$TITL PWAINI* PWAEUA — Analytic Tsai* cubic form factor

SUBROUTINE PWAI N I < XM I N , XMAX

)

PWAINI should be called once* after the common KIN has

been initialized in order to calculate constants which

do not change when only THETA-k Cor q-MU) changes.

IMPLICIT DOUBLE PRECISION CA-Z)

COMMON /PWAC0M/K1M2, K1M1, K10*K1P1* K1P2* K1P3* K2M2* K2M1, K20*

+ K2P1, K2P2.K2P3, J1M2*J1M1*J10,J1P1*J1P2* J2M2, J2M1 * J20* J2P1*

+ J2P2* L » LAMBDA

COMMON /CONST/ EM , EM2 * AM * HBARC . ALPH * P I E * ALPHA * RTD * DTR

COMMON /KIN/ EEX*TH, AMO* PI* P2* CTH* STH, PI 2* E 1 * P22 * E2* AMF*SDP3, SDP*

+ U32 * U3* UO* U2, CTHS » CTHP, STHP*XMU

LAMBDA=SDP-EM2

KAPPA: E1-E2-XMU-LAMBDA/AM0

CALL PWACOEC LAMBDA* KAPPA, P 1 , El * P2* E2* 1 . DO* Al* B1 * Cl* ETA 1 * + XMI N* BET0/BET1* BET2* TAU1*RH01* BETOP* BET IP, BET2P

)

CALL PWACOE (LAMBDA, KAPPA* P2* E2 * P 1 * E 1 * - 1 . DO * A2 * B2 * C2 * ETA2*

+ XMAX* DELO* DELI * DEL2* TAU2,RH02, DELOP* DELIP, DEL2P

)

GAMMA: ( ( L AMBD A + EM2 ) * ( 2 . DO*E 1 *E2-LAMBDA ) +LAMBDA*K APPA**2+

+ ( LAMBDA+ EM2 ) *< LAMBDA/ AMO + XMU ) *KAPPA ) /LAMBDA**2

ETA:1 . DO

TAUP:-2. D0*(LAMBDA**2-2. D0*EM2* ( LAMBDA+EM2 ) ) /LAMBDA

GAMMAP:2. D0*(LAMBDA**2-EM2*#2)/LAMBDA#*2

ETAP :-2 . DO

RHOP : C LAMBDA **2-2 . D0*EM2**2 )/LAMBDA**2

RETURN

C

C PWAEUA calculates the integral of Tsai from XI to XF assuming

C that the longitudinal form factor (W2) is a cubic defined

C by A* B* C* and D and that the transverse form factor (Wl)

C is a cubic given by AP, BP* CP* and DP.

ENTRY PWAEUA (XI , XF , A , B * C * D , AP , BP * CP* DP * UAL)

CALL PW AK I N ( X I * XF, LAMBDA* KAPPA* P1*E1*P2*E2*ETA1*A1*B1*C1*

+ Dll, DIF, J10* J1M1, J1P1*J1M2* J1P2*K10*K1M1*K1P1*K1 M2 , K1P2, K1P3)

CALL PW AK I N ( X I * XF* LAMBDA* KAPPA* P2* E2* PI* El * ETA2* A2* B2* C2*

+ D2 1 * D2F i J20* J2M1, J2P1 , J2M2* J2P2* K20 * K2M 1 * K2P 1 , K2M2* K2P2 , K2P3)

BETA:2. D0*KAPPA*DSQRT( LAMBDA* (LAMBDA+2. D0*EM2)

)

ZETA1I: (2. 0D0*LAMBDA*A1-B1 )*XI-2. D0*LAMBDA*B1 +C1

ZETA1F: (2. 0D0*LAMBDA*A1-B1 )*XF-2. DO*LAMBDA*B 1 +C

1

ZETA2I: (2. 0D0*LAMBDA*A2-B2)*XI-2. DO*L AMBDA*B2+C2

ZETA2F: (2. 0D0*LAMBDA*A2-B2) *XF-2 . D0*LAMBDA*B2+C2

C This is the formula as originally written in Maximon's notes

C L:-DLOG( ( ZET A2F+BETA*D2F ) / ( ZETA 1 F+BETA*D IF )

C + (ZETA1I+BETA*D1I)/(ZETA2I+BETA*D2I) )/BETA

C The following version of L is rewritten to avoid cancelation

S 1 1 : 1 . DO

EPS1 1:0. DO

IF (ZETA1I.LT. 0. DO ) THEN

5 1 I

— 1 .DO

EPS 1 1 :2 . DO

END I

F

S2I : 1 . DO

EPS2 I :0 . DO

IF (ZETA2I .LT.O. DO) THEN

52 I :- 1 . DO

EPS2I : 2 . DO

END I

F

70
5.8 59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

sir

i . do

EPS 1 F = 0 . DO

IF ( ZETA1F . LT . 0 . DO ) THEN

S1F=-1 .DO

EPS1F = 2. DO

ENDIF

S2F=1 . DO

EPS2F = 0 . DO

IF CZETA2F.LT. 0. DO) THEN

S2F=-1 .DO

EPS2F = 2 . DO

ENDIF

DEN=2. D0*EM*P1*P2*KAPPA*STH

L = C S IF *DLOG ( (DABSCZETA1F )+BETA*DlF )/DEN)-

+ SI I *DLOG ( ( DABS C ZETA1 1 ) +BETA*D1 1 ) /DEN )

+ S2F *DLOG ( C DABS CZETA2F)+BETA*D2F)/DEN)+

+ S2I*DL0GC ( DABS C ZET A2 1 > + BETA*D2I ) /DEN ) +

+ EPS IF #DLOG ( DABS C XF-2 . DO*LAMBDA) )-

+ EPS1 I*DLOGC DABSCXI-2. DO*LAMBDA) )-

+ EPS2F*DL0G( DABSCXF-2 . DO^LAMBDA) )+

+ EPS2I*DL0G(DABS(XI-2.D0*LAMBDA)) ) /BET A CALL PWAINTCXI, XF, BETO, BET1 , BET2, DELO/ DELI / DEL 2/ ETA/

+ TAU1/ T AU2 / RH01 / RH02/ GAMMA/ A/ B/ C/ D/ I

)

CALL PWAINTCXI, XF, BETOP, BET IP, BET2P, DELOP, DELIP, DEL2P, ETAP,

+ TAUP, TAUP/RHOP, RHOP, GAMMAP, AP, BP, CP, DP, IP)

OAL=I+IP

RETURN

END

71
12

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

$TITL PWACOE — Initialize constants constant with q-MU

SUBROUTINE PWACOE ( LAMBDA , KAPPA , PA , EA , PB , EB , SON , A , B, C » ETA,

+ X, BETO , BET1 * BET 2, T AU , RH0,BET0P,BET1P,BET2P)

C Note that the argument SGN will be either +1 or -1. This

C permits one to use the same formulae to calculate two sets

C of variables necessary for the PWA calculation.

IMPLICIT DOUBLE PRECISION (A-Z)

COMMON /CONST / EM, EM2, AM/ HBARC, ALPH, PIE, ALPHA, RTD, DTR

COMMON /KIN/ EEX, TH, AMO, PI, P2, CTH, STH, PI 2, El, P22, E2, AMF , SDP3, SDP,

+ U32, U3, U0,U2, CTHS, CTHP, STHP,XMU

Ar (PB+SGN*(EA*PB-EB*PA*CTH)/AM0)**2+(EM*PA*STH/AM0)**2 ' B=2. D0*PA*(PA*PB-(EA*EB-EM2)*CTH)*(PB+SGN*(EA*PB-EB*PA*CTH)/

+ AMO ) +SGN*2 . D0*EM2*PA**2*(EA-EB)*STH**2/AM0+2. DO*XMU*

+ ( -SGN*PB* (EA*P8-EB*PA*CTH) -LAMBDA* (LAMBDA+2.D0*EM2) / AMO

)

Cs4 . DO* ( PA*LAMBDA-SGN*XMU* ( PA*EB-EA*PB*CTH) )**2 +

+ (2. D0*EM*XMU*PB*STH)**2

ETA= ( (2. D0*EM*PA*PB*KAPPA*STH)/A)**2

X=2. DO* (LAMBDA* ( EA-EB-SGN*XMU ) * ( E A-EB-U3 ) )/

+ (1 . DO+SGN*(EA-EB-U3)/AMO)

PAB=SGN*PA*(PA-PB*CTH)

EAMU=EA*(EA-SGN*XMU)

EAAMOrl . D0+SGN*2 . DO*E A/AMO

BETO =8. D0*EM2*EAMU*LAMBDA*( PAB- ( EA+EB )*XMU)

BET 1 z-EM2* ( 4 . DO*EAMU* C PAB- ( EA+EB ) * ( KAPPA + XMU ) ) +

+ 2. DO*L AMBDA* CPAB-C EA+EB )*XMU)*EA AMO)

BET2=EM2*CPAB-( EA+EB)* (KAPPA+XMU) ) *EAAMO

TAU= (LAMBDA**2-4 . D0*EM2*EA*EB-SGN*2 . DO*XMU*

+ (EM2*(EA-EB-SGN*XMU)+EB*LAMBDA) ) /LAMBDA

RHO= ( (LAMBDA+EM2 )*(2. DO*EA*EB-L AMBDA ) +LAMBDA*KAPPA**2+

+ (LAMBDA+EM2)*(LAMBDA/AM0+XMU)*KAPPA)/LAMBDA**2+

+ EM2/AM0**2+. 5D0+SGN*E A/AMO

BET0P=-8. D0*EM2**2*LAMBDA* (PAB- (EA+EB )*XMU)

BET IP: 4 . D0*EM2* ( EM2* ( PAB- ( EA + EB )*( KAPPA + XMU ) ) +

+ LAMBDA*(PAB-(EA+EB)*XMU)

)

BET2P: -2 . D0*EM2*( PAB- ( EA+EB )*( KAPPA+XMU)

)

RETURN

END
12

3

45

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

$TITL PWAKIN — Evaluate constants which change with q-MU

SUBROUTINE PW AK I N ( X I , XF , LAMBDA , KAPPA , PA , EA , PB , EB , ETA, A, B, C, + DI , DF, JO, JM1, JP1, JM2, JP2,K0, KM1, KP1, KM2, KP2, KP3

)

IMPLICIT DOUBLE PRECISION (A-Z)

COMMON /CONST/ EM, EM2 , AM , HB ARC , ALPH , P IE , ALPHA , RTD , DTR

COMMON /KIN/ EEX, TH, AMO, PI, P2, CTH, STH, PI 2, El , P22, E2, AMF , SDP3 , SDP

,

+ U32»U3,U0» U2 , CTHS , CTHP, STHP » XMU

CHI=B/A

D I : DSQRT ( A* ( ( XI -CH I ) **2+ETA )

)

DF r DSQRT ( A* ( ( XF-CH I )**2+ETA)

)

SQTAr DSQRT ( A

)

SQTCzDSQRT(C)

C The following is the formula as originally written

C JO=DLOG( ( A*XF-B+DF*SQTA)/(A*XI-B+DI#SQTA) )/SQTA

C This is the version as corrected to reduce "cancelation"

SI = 1.. DO

IF (A*XI-B.LT.O. DO) SI=-1.D0

SF;1 .DO

IF ( A*XF-B.LT.O. DO) SF=-1.D0

DEN=2. DO*EM*PA*PB*KAPPA*STH

J0= (SF*DLOG( ( DABS ( A*XF-B)+DF *SQTA ) /DEN ) -

+ SI*DLOG( (DABS( A*XI-B)+DI*SQTA)/BEN) )/SQTA

C The following is the formula as originally written

C JM1 =-DLOGCXI/XF*(-B*XF+C+SQTC*DF)/

C + (-8*XI+C+SQTC*DI) )/SQTC

C This version is corrected to avoid cancelation

SI=1.D0

IF (B*XI-C.LT.O. DO) SI=-1.D0

SF= 1 . DO

IF (B*XF-C.LT.O.DO) SF=-1.D0

JM1 z ( SF *DLOG ( CDABSC B*XF-C)+SQTC*DF)/(XF*DEN) )-

+ SI *DLOG ( ( DABS( B*XI-C ) +SQTC*DI )/<XI*DEN) ) )/SQTC

JP1= ( DF-DI ) / A + CH I * J 0

JM2= (-DF/XF+DI/XI+B*JM1 )/C

JP2 - (XF*DF-XI*DI+3. DO*B* JP1-C* JO ) / ( 2 . DO* A

)

KOr C CA*XF-B)/DF-(A*XI-B)/DI )/(ETA*A**2)

KPl=-( (-B>i<XF + C)/DF-(-B*XI+C)/DI ) / ( ETA*A**2 ) KP2= ( JO+2. D0*B*KP1-C*K0)/A

KP3: (XF* + 2/DF-XI**2/DI+3. D0*B*KP2-2. D0*C>*<KP1 ) /A

KM1= ( 1 . DO/DF-1 . DO/D 1+ JM1 +B*KO)/C

KM2= (-1 . D0/(XF*DF) + 1 . DO/(XI*DI )+3. D0*B*KMl-2. DO*A*KO)/C

RETURN

END

73
12

3

4

5

6

7

8

9

10

1

1

12

13

14

15

16

17

18

19

20

21

22

23

STITL PWA1NT — Determine integral of long, or trans term.

SUBROUTINE PWAINT ( XI , XF, BETO, BET1 , BET2, DELO, DELI , DEL2, ETA,

+ T AU1 »TAU2,RH01,RH02, GAMMA, A, B, C, D, I

)

IMPLICIT DOUBLE PRECISION (A-Z)

COMMON /PWAC0M/K1M2,K1M1, K10,K1P1, K1P2, K1P3, K2M2 , K2M1 , K20

,

+ K2P1 , K2P2, K2P3, J1M2, J1M1 , J 1 0, J1P1#J1P2, J2M2 » J2M 1 * J20, J2P 1 , . + J2P2, L, LAMBDA

COMMON /KIN/ EEX, TH, AMO, PI , P2, CTH, STH, P12, El, P22, E2, AMF , SDP3, SDP,

+ U32,U3,U0,U2, CTHS , CTHP, STHP,XMU

CALL PWAE INC BETO, BET1 , BET2, K1M2, K1M1 ,K10,K1P1,K1P2,K1P3,

+ J1 M2, J1M1,J10,J1P1,J1P2, A, B,C,D, 111,131,141)

CALL PWAE INC DEL 0, DELI , DEL2, K2M2, K2M1 , K20, K2P1 , K2P2, K2P3,

+ J2M2, J2M1 , J20, J2P1 , J2P2, A, B,C, D, 112, 132, 142)

I2=A*< 1 . DO/XI-1 . DO/XF)+B*DLOG(XF/XI )+C*(XF-XI )+

+ . 5D0*D*(XF**2-XI**2)

15: ( A+B*2. D0*LAMBDA+C*4. D0*LAMBDA**2+D*8 . D0*LAMBDA**3 ) *L-

+ ( B + 2 . D0*LAMBDA*C+4.D0#LAMBDA*#2#D)*(J10-J20)-

+ CC+2. DO*LAMBDA*D)*( J1P1-J2P1 )-

+ D*( J1P2-J2P2)

I:-I11-I12-ETA#I2/U3+TAU1*I31-TAU2*I32-

+ RH01*I41+RH02*I42-GAMMA*I5

RETURN

END

74
1 STITL PWAEIN -- Evaluate some terms in the analytic integral

2 SUBROUTINE PWAEINC BETO, BET1 , BET2, KM2, KM1 / KO, KP1 , KP2, KP3»

3 + JM2, JM1. JO, JP1, JP2, A, B, C, D, 11/ 13, 14

)

4 IMPLICIT DOUBLE PRECISION (A-Z)

5 Il=BET0*A*KM2+( BETO*B+BET 1 *A )*KM1+

6 + (BET0*C+BET1*B+BET2*A)*K0+(BET0*D+BET1*C+BET2#B)*KP1+

7 + <BET1*D+BET2*C)*KP2+BET2*D*KP3

8 I3=A*JM2+B*JM1+C*J0+D*JP1

9 I4=A*JM1+B*J0+C*JP1+D*JP2

10 RETURN

11 END

75
12

3

4

5

6

7

6

9

10

11

12

13

14

15

16

17

non

ST1TL SPCE — Series expansion of Spence function

FUNCTION SPCE ( X

)

Series expansion of the Spence function for range 0

Note that the number of terms used in the series is

to auo i d floating underflow.

IMPLICIT DOUBLE PRECISION (A-H,0-Z>

DOUBLE PRECISION SPCE

SPCEsO. ODO

IF (X.LE.O.DO) RETURN

TOP = 1 . ODO

N=1 . DO-25. DO/DLOGCX)

DO 10 1=1,

N

TOP=TOP*X

SPCE=SPCE+T0P/I**2

10 CONTINUE

RETURN

END

to 0.5

chosen

76

I
1 STITL PWAEIN — Evaluate some terms in the analytic integral

2 SUBROUTINE PWAEINC BETO, BET1 , BET2, KM2, KM1 , KO, KP1 , KP2, KP3,

3 + JM2» JM1*J0>JP1*JP2/A»B»C,D» 1 1 > I3» 14)

4 IMPLICIT DOUBLE PRECISION CA-Z)

5 Il=BETO*A*KM2+(BETO*B+BETl*A)*KMl+

6 + <BET0*C+BET1*B+BET2*A)*K0+(BET0*D+BET1*C+BET2*B)*KP1+

7 + <BET1*D+BET2*C)*KP2+BET2*D*KP3

8 I3rA*JM2+B*JMl+C*J0+D*JPl

9 I4=A*JM1+B*J0+C*JP14D*JP2

10 RETURN

11 END

75
12

3

4

5

6

7

6

9

10

11

12

13

14

15

16

17

o

n

o

f

STITL

10

SPCE -- Series expansion of Spence function

FUNCTION SPCE C X

)

Series expansion of the Spence function for range 0

Note that the number of terms used in the series is

to avoid floating underflow.

IMPLICIT DOUBLE PRECISION (A-H,0-Z)

DOUBLE PRECISION SPCE

SPCE=0. ODO

IF (X.LE.O.DO) RETURN

TOP r 1 . ODO

N=1 . DO-25. DO/DLOGCX)

DO 10 1=1,

N

TOP=TOP*X

SPCE=SPCE+T0P/I**2

CONTINUE

RETURN

END

to 0.5

chosen
1 $TITL SPENCE — The Spence function

2 FUNCTION SPENCE C X

)

3 C Program for the computation of the Spence function

4 C values of the argument x 5 IMPLICIT DOUBLE PRECISION (A-H,0-Z)

6 DOUBLE PRECISION SPENCE

7 DATA PI23, PI26/3. 289868134D0, 1 . 644934067D0/

a IF (X.GT.2.0D0) GO TO 1 9 IF (X. GT. 1 .ODO) GO TO 2

10 IF (X.GT.0.5D0) GO TO 3

n IF <X. LT. -1 . DO) GO TO 4

12 IF (X.LT.O.ODO) GO TO 5

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

SPENCE=SPCE<X)

RETURN

Y = 1 . ODO/X

SPENCE = PI 23-SPCE < Y ) - . 5D0*DL0G ( X ) **2

RETURN

Y=(X-1 . ODO)/X

Yl: (X-l .. ODO ) * Y SPENCE=PI26+SPCE(Y)-0. 5D0*DL0G <X)*DL0G(Y1)

RETURN

Y=1 .ODO-X

SPENCE = PI 26-SPCE ( Y ) -DLOG < X ) *DLOG( Y

)

RETURN

Y=1 . ODO-X

Yl=l .ODO/Y

SPENCE =-PI 26+SPCE< Yl )-0. 5D0*DL0G < Y ) *DLOG < X*X/Y

)

RETURN

Y zX/ (X— 1 . ODO)

Yl=-X+1 . ODO

SPENCE = -SPCE C Y ) -0 . 5 DO* DLOG ( Yl

)

RETURN

END

for all real

77
12

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

16

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

*TITL EUALI — Simpson's rule integration

FUNCTION EUALI ( A / B / FUNC / NS)

C Function to evaluate the integral of the function FUNC

C betuieen 6 and B. The number of steps starts at NSTEP

C as defined in EUALI and then goes up by a factor of two

C until the integral stops changing by more than a factor

C of EPS. In this version/ which uses Simpson's rule/ the

C new values needed when NSTEP increases are not all reC calculated. Ualues computed during the calculation of

C the old (inaccurate) value of the integral are re-used.

IMPLICIT DOUBLE PRECISION (A-H/O-Z)

EXTERNAL FUNC

DOUBLE PRECISION EUALI

COMMON /EU I COM/ NSTEP/ EPS

CALL FUNC( A/ EUALA)

CALL FUNC ( B/ EUALB

)

NSrNSTEP-1

STEPO= ( B- A ) /"NSTEP

EUAL2=0. DO

DO 900 1=1 /NS

CALL FUNC (A+I*STEPO/ U)

EUAL2 = EUAL2 + U

900 CONTINUE

NS=NSTEP

C

C Start of loop

4 STEP: . 5D0*STEP0

EUAL4=0. DO

DO 901 1=1/ NS

CALL FUNC (A+STEP+C 1-1 ) *STEPO/ U

)

EUAL4=EUAL4+U

901 CONTINUE

EU= (EUALA + 2. D0*EUAL2+ 4 . D0*EUAL4+EUALB ) *STEP0/'6 . DO

IF (NS. NE. NSTEP) THEN

IF (EPS. LE. 0. DO) THEN

IF (DABS(EU-EUALIO) . LE . DABS ( EU#EPS ) ) GO TO 6

CD WRITE (1,*) A/ B/ DABS (EU-EUALIO) / DABS( EU*EPS

)

ELSE

IF ( D ABS ( EU-EUAL 1 0 ) . LE . EPS ) GO TO 6

CD WRITE (1,*) A/ B/ DABS ( EU-EUAL 1 0 ) / EPS

END IF

END I

F

. EUALIOrEU

EUAL2 = EUAL4 + E VAL2

NS = 2*NS

STEPO= STEP

GO TO 4 6 EUAL I =EU

RETURN

END

78
12

3

A5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

3e

39

40

41

42

JTITL EUALI -- Gaussian integration

FUNCTION EUALI(A,B,FUNC,NS)

C Function to evaluate the integral of the function FUNC

C between A and B. The number of intervals starts at NSTEP

C as defined in EUICOM and then goes up by a factor of two

C until the integral stops changing by more than a factor

C of EPS. In this version, NGAUSS-point Gaussian integration

C is used (for particularly tough functions).

PARAMETER (NGAUSSr 16)

IMPLICIT DOUBLE PRECISION (A-H,0-2)

EXTERNAL FUNC

DOUBLE PRECISION EUALI , 2 ( NGAUSS ) > W ( NGAUSS

)

COMMON /EUICOM/ NSTEP, EPS

NS=NSTEP

C

C Start of loop

4 STEP = ( B—A ) /NS

EUrO.DO

DO 902 1=1, NS

AX=A+( 1-1 )*STEP

BX= AX+STEP

CALL GSET ( AX, BX, NGAUSS , 2, W

)

DO 903 J=l, NGAUSS

CALL FUNC (2(J),U)

EU=EU+U*W( J)

903 CONTINUE

902 CONTINUE

IF (NS. NE. NSTEP) THEN

IF (EPS. LE. 0. DO) THEN

IF ( DABS ( EU-EUALIO ) . LE . DABS ( EU*EPS ) ) GO TO 6

CD WRITE (1,*) A, B, DABS(EU-EUALIO) , DABS(EU*EPS)

ELSE

IF (DABS(EU-EUALIO) . LE . EPS) GO TO 6

CD WRITE (1,*) A, B, DABS ( EU-EUAL I 0 ) , EPS

ENDIF

END IF

EUALIOrEU

NS= 2*NS

GO TO 4

6 EUALI=EU

RETURN

END

79
1 ST1TL GSET — N-point Gauss zeros and weights

z SUBROUTINE GSET ( AX , BX , NX , Z , W

)

3 C N-point Gauss zeros and weights for the

4 c stored in arrays Z and W respectively.

5 c 6 c This version converted from the version

7 c for the CDC-7600 to a double precision v 8 c mach ine

.

9 c 10 IMPLICIT DOUBLE PRECISION CA-H,0-Z>

11 DOUBLE PRECISION Z(21»W(21>A(2731/X(2731

12 INTEGER KTABC96)

13 c 14 c Table of initial subscripts for N=2(l)16

15 DATA KTAEK21/1/

16 DATA KTABOl/2/

17 DATA KTAB C 4 ) /4/

18 DATA KTAB ( 5 ) /6/

19 DATA KT AB C 6 ) /9/

20 DATA KTAB(71/12/

21 DATA KTAB ( 8 1 / 16/

22 DATA KTAB ( 9 ) /20/

23 DATA KTAB( 101/25/

24 DATA KTAB( 11 1/30/

25 DATA KTAB( 121/36/

26 DATA KTAB( 131/42/

27 DATA KTAB(141/49/

28 DATA KTAB ( 151/56/

29 DATA KTAB C 161/64/

30 DATA KTAB ( 20 ) /72/

31 DATA KTAB ( 24 ) /82/

32 DATA KTAB ( 28 1 /B2/

33 DATA KTABC321/94/

34 DATA KTAB ( 36 ) /94/

35 DATA KTAB (401/110/

36 DATA KTAB(441/110/

37 DATA KTAB(481/130/

38 DATA KT AB ( 52 1/130/

39 DATA KTAB(56)/130/

40 DATA KTAB(601/130/

41 DATA KTAB (64 1/154/

42 DATA KTAB(681/154/

43 DATA KTABC721/154/

44 DATA KTAB(761/154/

45 DATA KTAB(801/ie6/

46 DATA KTAB(841/186/

47 DATA KTABC881/186/

48 DATA KTABC921/186/

49 DATA KT AB ( 96 1/226/

50 c 51 c Table of abscissae (XI and weights (A! fi 52 c 53 c N = 2

54 DATA X( 1 1 /O. 577350269189626D0 /, A(ll/1

55 c N = 3

56 DATA X(21/0. 774596669241483D0 /, A(21/0

57 DATA X(31/0. OOOOOOOOOOOOOOODO /, A(3)/0( 310. 888888888888B89D0

80
58 C N = 4 . 59 DATA X <4)/0. 861 13631 1594053D0 A ( 4

)

SO

60 DATA X (5)/0. 3399810435B4856D0 A (

) SO

61 C N = 5

62 DATA X ( 6 ) /"0 . 9061 798 45938 664 DO / , A ( 6

)

so

63 DATA X <7)/0. 5384693101 05683D0 / , A ( 7

)

so

64 DATA X ( 8 ) / O . 000 000 000 000 000DO / , A ( 8

)

SO

65 C N = 6

66 DATA X ( 9) so. 9324695 1 4203 152D0 / , A ( 9

)

SO . 67 DATA X ( 10)/0. 661209386466265D0 / , A ( 1 0

)

SO

68 DATA X (11 )/0. 23861 9186083 197D0 / , A( 11 )sO

69 C N = 7

70 DATA X ( 12)/0. 949107912342759D0 / , A ( 1 2

)

SO

71 DATA X < 13 )S0. 741531 185599394 DO / , A ( 1 3 ) /O

72 DATA X ( 14 )/0. 405845 151 377397D0 / , A ( 141/0

73 DATA X ( 15 ) SO

000000000 000000 DO / , A( 15)/0

74 C N = 8

75 DATA X ( 1

) /O . 96 028 9 85 6 49 7 53 6 DO / , A( 16)/0

76 DATA X ( 17)/0. 7966664774 13627D0 /, A ( 17 ) SO

77 DATA X( 18) SO. 52553240991 6329D0 / , A ( 18) SO

78 DATA XC19)/0. 183434642495650D0 / , A ( 1

) SO

79 C N = 9

80 DATA X C 20 ) SO . 9681 60239507 626D0 / , A ( 20 ) SO

81 DATA X (21 1/0. 836031 107326636D0 / , A ( 2

) SO

82 DATA X ( 22 ) /O . 613371 4 3 270 0590 DO s

,

A ( 22 ) SO

83 DATA X ( 23 ) /O . 324 25 34 2 3403 809 DO / , A ( 23) SO

84 DATA X ( 24 ) /Q . OOOOOOOOOOOOOOODO st A ( 24 ) SO

85 C N=10

86 DATA X (25)/0. 97390652851 7 172D0 s

,

A(25)/0

87 DATA X (26)/0. 8650633566S8985D0 St A ( 26 ) SO

88 DATA X (27)/0. 679409568299024 DO St A ( 27 ) SO

89 DATA X (28)/0. 433395394 129247 DO s

,

A ( 28 ) SO

90 DATA X (29)/0. 148874338981631D0 Sr A(29)/0

91 c N = 1

92 DATA X ( 30 )/0. 978228658 146057 DO S , A ( 30 ) SO

93 DATA X (31 )/0. 88706259976S095D0 St A ( 31 )S0

94 DATA X (32)/0. 730152005574049D0 St A ( 32 ) SO

95 DATA X (33)/0. 51 9096 12920681 2D0 s

,

A( 33) SO

96 DATA X (34)/0. 269543155952345D0 St A(34)/0

97 DATA X ( 35 ) /O . OOOOOOOOOOOOOOODO s

,

A ( 35 ) SO

98 c N = 12

99 DATA X (36)/0. 981 5606342467 1 9D0 s

,

A ( 36 ) SO

100 DATA X (37)/0. 90411 7 25 6 37 04 75 DO St A ( 37 ) SO

101 DATA X (38)/0. 7699026741 94 305D0 St A ( 38 ) SO

102 DATA X (39)/0. 58731 79542866 1 7 DO s

,

A ( 39 ) SO

103 DATA X (40)/0. 36783 14989981 80 DO s

,

A ( 40 ) SO

104 DATA X (41 )/0. 12523340851 1469D0 s

,

A ( 4

) SO

105 c N = 13

106 DATA X ( 42 ) /O . 984 18305471 8588D0 s

,

A(42)/0

107 DATA X (43)/0. 91 7 598 399 22 2 97 8 DO St A ( 43 ) SO

108 DATA X (44 )/0. 80 1 57809073331 ODO s

,

A ( 44 ) SO

109 DATA X (45)/0. 642349339440340D0 s

,

A ( 45

)

SO

110 DATA X (46)/0. 44 849275 1036447 DO s

,

A ( 46 ) SO

111 DATA X ( 47 ) /O . 2304583 15955 135D0 s

,

A ( 47 ) SO

112 DATA X ( 48 ) SO . 00 000 000 0000 000 DO s

,

A ( 48 ) SO

113 c N = 14

114 DATA X ( 49 )/0. 9862838086968 12D0 s

,

A ( 49 ) SO

, 101228536290376D0

222381 034453374 DO

313706645877887D0

3626 8 3783 37 83 62 DO

08 1 27438836 1574D0

06667 1344308688D0

2692667 19309996D0

272925086777901 DO

1 069393259953 18D0

160078328543346D0

203167426723066D0

233492536538355D0

/

/

/

/

/

/

/

/

/

S

s

s

s

s

ss

040484004765316DO /

092121499837728D0 /

138873510219787D0 /

1781459B0761946D0

20781 6047536889D0

2262831 80262897 DO

23255 1553230874 DO

035119460331752D0

81
115 DATA X ( 50 ) /O . 928 4 34 08 366 3 57 4 DO / , A(50)/0.

116 . DATA X ( 5 1 )/0. 827201 31 5069765D0 / , A( 51 ) /O.

117 DATA XC52)/0. 68729290481 1685D0 / , A(52)/0.

118 DATA X ( 53 )/0. 5 15248636358 154 DO / , A(53)/0.

119 DATA X(54)/0. 3191 1 2 36 8927 890 DO / , A(54)/0.

120 DATA X ( 55 ) /O . 108 05 4 94 8707 34 4 DO / , A(55)/0.

121 C N = 1

5

122 DATA XC56)/0. 98799251 8020485D0 /

,

A(56)/0.

123 DATA X ( 57 ) /O . 937 27 3392 400706 DO / , A(57)/0.

124 DATA X ( 58 )/0. 8482065834 10427D0 /, A(58)/0.

125 DATA X(59)/0. 724417731360170D0 /, A ( 59 ) /O

.

126 DATA X(60)/0. 570972 17260B539D0 // A( 60) /O.

127 DATA X ( 6 1 )/0. 394151347077563D0 / , AC61 )/0.

128 DATA X ( 62 ) /O . 201 194093997 435 DO /

,

A(62)/0.

129 DATA X ( 63 ) /O . OOOOOOOOOOOOOOODO / , A(63)/0.

130 C N=16

131 DATA X ( 64 ) /O . 98940093499 1650D0 / , A(64)/0.

132 DATA X ( 65 ) /O . 94 4 57 5 02 307 323 3 DO /, A(65)/0.

133 DATA X(66)/0. 86563 1202387832D0 // A ( 66 ) /O

.

134 DATA X(67)/0. 755404408355003D0 /, A ( 67 ) /O

.

135 DATA X ( 68 )/0. 6 1787624440264 4D0 /, A ( 68 ) /O

.

136 DATA X ( 69 ) /O . 45801 677 7 657 22 7 DO / , A ( 69 ) /O

.

137 DATA X(70)/0. 281 60355077 9259D0 /

,

A<70)/0.

138 DATA X ( 7 1 >/0. 09501 2509837637D0 / , A ( 7 1 )/0.

139 c N:20

140 DATA X < 72 )/0. 993 128599 185094 DO / , A ( 72 ) /O

.

141 DATA X(73)/0. 963971 9272779 13D0 / , A ( 73 ) /O

.

142 DATA X ( 74 ) /O . 912234428251 325D0 / , A(74)/0.

143 DATA X(75)/0. 839 11 6971 8222 18D0 / , A(75)/0.

144 DATA X ( 76 )/0. 74 633 1906460 150D0 /, A ( 76 ) /O

.

145 DATA XC77)/0. 6360536807265 15D0 /, A ( 77 ) /O

.

146 DATA X(79)/0. 51 0867001 950827D0 /, A(78)/0.

147 DATA X(79)/0. 37 37060887 154 19D0 // A(79)/0.

148 DATA X(80)/0. 227785851 14 1645D0 /, AC 80/0.

149 DATA XC81 )/0. 076526521 133497D0 // AC81 )/0.

150 c N r 24

151 DATA X ( 82 ) /O. 995 187 2 1999702 IDO /, A ( 82 ) /O

.

152 DATA XC83)/0. 97472855597 1309D0 // AC 83) /O.

153 DATA X(84)/0. 93827 4552002732D0 /, A C 84 ) /O

.

154 DATA X(85)/0. 8864 1552700440 IDO /, AC 85)/0.

155 DATA X(86)/0. 82000 1985973902D0 /, A C 86 ) /O

.

156 DATA X ( 87 )/0. 740124 191 578554 DO /, A C 87 ) /O

.

157 DATA X(88)/0. 64809365 1936975D0 // AC 88) /O.

158 DATA X(89)/0. 545421 47 1388839D0 /, A C 89 ) /O

.

159 DATA XC90)/0. 4 337 93507 62604 5D0 /, A C 90 ) /O

.

160 DATA X(91 )/0. 315042679696163D0 /, AC91 )/0.

161 DATA X ( 92 ) /O . 19111 886747 36 1 6D0 /, A C 92 ) /O

.

162 DATA X(93)/0. 064056892862605D0 /, AC93)/0.

163 c N = 32

164 DATA X(94)/0. 99726386 184 948 IDO /, AC94)/0.

165 DATA X(95)/0. 98561 151 1545268D0 /, A C 95 ) /O

.

166 DATA X(96)/0. 964762255587506D0 /, A C 96 ) /O

.

167 DATA X(97)/0. 934 906075937 739D0 /, A C 97 ) /O

.

168 DATA X(98)/0. 896321 155766052D0 /, AC 98)/C.

169 DATA X ( 99 ) /O . 84936761 3 73 2 56 9 DO / , AC99)/0.

170 DATA XC 1 00 )/0. 7 944 8379596794 2D0/, AC 100) /O

171 DATA X( 101 )/0. 7321821 18740289D0// AC 101 )/0

000158007 159 7 60 DO /

121518570687903D0 /

1572031 67 1581 94 DO /

1 85 5 38 397 4 77 9 38 DO /

205198463721296D0 /

2152638534631 58 DO /

0307532419961 17 DO /

070366047488 108D0 /

107 159220467 172D0 /

1 39570677926 154D0 /

1 662692058 16994 DO /

186161 0000 155 62 DO /

1984314853271 11D0 /

20257824 1925561 DO /

027 152459411754D0 /

062253523938648D0 /

095158511682493D0 /

124628971 255534 DO /

1 495959888 16577D0 /

1691 565 1 9395003D0 /

1 826034 15044924D0 /

1 894506 10455069D0 /

0176140071 39 1 52 DO /

040601 4 29 8 00 3 86 DO /

0626720483341 09 DO /

083276741 576704 DO /

1019301 19817240D0 /

1 18194531961518D0 /

131 6886384491 76D0 /

142096109318382D0 /

149172986472603D0 /

1 52753387 130725D0 /

01234 1229799987D0 /

028531 3B8628933D0 /

0442774388174 19D0 /

059298584915436D0 /

07334648141 1080D0 /

086190161531953D0 /

0976186521041 13D0 /

1074442701 15965D0 /

1 15505668053725D0 /

121670472927803D0 X

1 25837 4 56346828D0 /

1 27930 195346752D0 /

00701 8610009470D0 /

016274394730905D0 /

025 3 92 065 3 09 2 62 DO /

034273862913021D0 /

0428 35 8980222 26 DO /

050998059262376DO /

. 058684 09 347 8 535 DO^

. 065822222776361 DO/

82
172 DATA Xt 102

173 DATA X( 103

174 DATA Xt 104

175 DATA Xt 105

176 DATA Xt 106

177 . DATA Xt 107

178 DATA Xt 108

179 DATA Xt 109

180 C N = 40

181 DATA Xt 110

182 DATA Xt 111

183 DATA Xt 112

184 DATA Xt 113

185 DATA Xt 114

186 DATA Xt 115

187 DATA Xt 116

188 DATA Xt 117

189 DATA Xt 118

190 DATA Xt 1.1-9 191 DATA Xt 120

192 DATA Xt 121

193 DATA Xt 122

194 DATA Xt 123

195 DATA Xt 124

196 DATA Xt 125

197 DATA Xt 126

198 DATA Xt 127

199 DATA Xt 128

200 DATA Xt 129

201 c N = 4B

202 DATA Xt 130

203 DATA X t 1 3

204 DATA Xt 132

205 DATA Xt 133

206 DATA Xt 134

207 DATA Xt 135

208 DATA Xt 136

209 DATA Xt 137

210 DATA Xt 138

211 DATA Xt 139

212 DATA Xt 140

213 DATA Xt 141

214 DATA Xt 142

215 DATA Xt 143

216 DATA Xt 144

217 DATA Xt 145

218 DATA Xt 146

219 DATA Xt 147

220 DATA Xt 148

221 DATA Xt 149

222 DATA Xt 150

223 DATA Xt 151

224 DATA Xt 152

225 DATA Xt 153

226 c N = 64

227 DATA Xt 154

228 DATA Xt 155

)/0. 663044266930215D0/, At )/0. 587715757240762D0/» At )/0. 506899908932229D0// At )/0. 4213512761 30635D0// At )/0. 331868602282127D0/, At

) / O . 2392873622521 3700// At

>/0. 144471961582796D0/, At

)/0.048307665687738D0/, At

)/0. 998237709710559D0// At

)/0. 99072623B699457D0// At

) / O . 97 72 5 99 4 99 8 37 7’ 4 DO// At )/0. 95791 68 19213791D0// At )/0. 932812808278676D0// At

)/0. 902098806968874D0// At

)/0. 86595950321 225 9D0/» At )/0. 82461223083331 IDO// At )/0. 778305651426519D0// At )/0. 727318255189927D0// At )/0. 671956684614179D0/. At )/0. 612553889667980D0/, At )/0. 549467125095123D0// At

)/0. 483075801686178D0// At

)/0, 413779204371605D0// At )/0. 341994090825758D0// At

)/0.268152185007253D0// At

)/0. 192697580701371D0// At )/0. 1 16084O70675255D0// At

)/0. 03877241750605OD0// At

)/0. 998771007252426D0// At )/0. 993530172266350D0/, At )/0. 984124583722826D0/, At

) / O . 97059 1 59254 62 47D0// At

)/0. 952987703160430D0// At

)/0.931386690706554D0/, At

)/0. 905879136715569D0// At )/0. e76572020274247D0// At

) /O . 84358826 162 43 93D0// At )/0. 807066204029442D0/, At )/0. 767159032515740D0// At

)/0.724034130923B14D0// At

)/0. 677872379632663D0// At

)/0. 628867396776513D0// At

)/0. 577224726083972D0// At 1/0.52316097 4 722233D0// At 1/0.466902904750958D0// At )/0. 408686481990716D0// At

) / O . 34 8755886292 1 6 ODO// At

I/O. 2873624B7355455D0// At

1/0.224763790394689D0// At

I/O. 161222356068891D0// At

)/0. 097004699209462D0// At )/0. 032380170962869D0// At

I/O. 999305041735772D0// At 1/0.9963401 1677 1955D0/, At

1021/0. 072345794 108B48D0/

1031/0.07819389578707000/

104 1/0. 08331 1924226946D0/

1051/0. 087652093004403D0/

1061/0. 091 173878695763D0/

1071/0. 093844399080804DO/

1081/0. 095638720079274D0/

1091/0. 0965400885 14727 DO/

1101/0. 00452 127709e533D0/

111 1/0. 010498284531 152D0/

1 121/0. 01 642105838 1907D0/

1131/0. 022245849 194 166D0/

1 14 1/0. 027937006980023D0/

1151/0. 033460 195282547D0/

1161/0. 038782 167974472D0/

1171/0. 0438709C8 1 85673D0/

1181/0. 048695807635072D0/

1191/0. 053227846983936D0/

1201/0. 05743976909939 IDO/

1211/0. 061 30624 2492928D0/

1221/0. 06480401 345660 IDO/

1231/0. 067912045815233D0/

124 1/0. 07061 1 64739 1286D0/

1251/0. 072886582395804D0/

1261/0. 074723169057968D0/

1271/0. 076 11 036 1900626D0/

1281/0. 07703981 8 164247 DO/

1291/0. 07 7 50594 7 97 8424 DO/

1301/0. 003 1 5334 6052305D0/

131 1/0. 007327553901 276D0/

1321/0. 01 1477234579234 DO/

1331/0. 01 55793 1572294 3D0/

1341/0. 0196161 60 4 57 35 5 DO/

1351/0. 023570760839324 DO/

1361/0. 027426509708356D0/

137 1 /O. 031 167227832798D0/

1381/0. 034777222564770D0/

1391/0. 03824 1351 065B30D0/

1401/0. 041545082943464D0/

141 1/0. 044674560856694D0/

1421/0. 04761 6658492490D0/

1431/0. 050359035553854 DO/

1441/0. 052890 189485 193D0/

1451/0. 055 199503699984D0/

1461/0. 057277292 100403D0/

1471/0. 059 11 4839698395D0/

1481/0. 060704439 165893D0/

1491/0. 062039423 159892D0/

1501/0. 0631 14 192286254 DO/

151 1/0. 063924238584648D0/

1521/0. 064466 164 435950D0/

1531/0. 06473769681 2683D0/

154 1/0. 00 178328072 1696D0/

1551/0. 004 14 7033260562D0/

83
229 DATA XC 156

230 DATA XC 157

231 DATA XC 158

232 DATA XC 159

233 DATA XC 160

234 DATA XC 161

235 DATA XC 162

236 DATA XC 163

237 DATA XC 164

238 DATA XC 165

239 DATA XC 166

240 DATA XC 167

241 DATA XC 168

242 DATA XC 169

243 DATA XC 170

244 DATA XC 171

245 DATA XC 172

246 DATA XC 173

247 DATA XC 174

248 DATA XC 175

249 DATA XC 176

250 DATA XC 177

251 DATA XC 178

252 DATA XC 179

253 DATA XC 180

254 DATA XC 181

255 DATA XC 182

256 DATA XC 183

257 DATA XC 184

258 DATA XC 185

259 C N = 80

260 DATA XC 186

261 DATA XC 187

262 DATA XC 188

263 DATA XC 189

264 DATA XC 190

265 DATA XC 191

266 DATA XC 192

267 DATA XC 193

268 DATA XC 194

269 DATA XC 195

270 DATA XC 196

271 DATA XC 197

272 DATA XC 198

273 DATA XC 199

274 DATA X C 200

275 DATA X C 201

276 DATA X C 202

277 DATA X C 203

278 DATA X C 204

279 DATA X C 205

280 DATA X C 206

281 DATA XC 207

282 DATA X ( 208

283 DATA X C 209

284 DATA XC210

285 DATA X C 21 1

)/0. 991013371476744D0/

)/0. 983336253884625D0/

) / 0 . 9 7 33268 2 77 8 991 ODO/

>/0. 961008799652053D0/

)/0. 94641137 4858402D0/

)/0. 929569 1721 3 1939D0/

)/0. 9 105221 37078502D0/

)/0. 8893154459951 14D0/

I/O. 865999398 154092D0/

)/0. 840629296252580D0/

)/0 . 8132653 15122797D0/

)/0. 78 39723589 4 334 IDO/

)/0. 7528 1990726053 IDO/

)/0. 71 98818501 7 1610D0/

)/0. 68523631 3054233D0/

)/0 . 64896547 1254657D0/

)/0.611155355172393D0/

)/0.571B95646202634DO/

)/0. 53 12794640 19894D0/

)/0. 4894031 45707052D0/

)/0. 446366017253464D0/

)/0. 4022701 5796399 IDO/

)/0. 357220 158337668D0/

)/0. 3 1132287 199021 ODO/

)/0,264687162208767D0/

I/O. 2 1742364 3740007D0/

)/0. 169644420423992D0/

)/0. 1214628 19296 120D0/

)/0. 0729931 2 1787799D0/

)/0. 024350292663424D0/

)/0. 99955382265 1630D0/

)/0 . 99764S864398237D0/

)/0.994227540S65688D0/

)/0. 98929 1302499755D0/

)/0.982848572738629D0/

)/0. 974909 140585727D0/

)/0. 965485089043799D0/

)/0. 95459076634363400/

)/0.942242761309872D0/

)/0. 928459877 172445D0/

)/0 . 9132631 02571757D0/

)/0. 896675579438770D0/

)/0. 87872256767821 3D0/

)/0. 85943 140666311 IDO/

)/0. 83883 147 3580255D0/

)/0. 8 169541 3868 1463D0/

)/0. 7938327 17504605D0/

)/0. 7695024 20 13504 IDO/

)/0. 744000297583597D0/

)/0. 7173651B5362099D0/

)/0. 68963764434 2027D0/

)/0. 6608598989861 19D0/

)/0. 63 107577304 687 IDO/

)/0. 60033062282975 IDO/

)/0. 56867 1 2681 22709D0/

)/0. 5361 4 59208971 31 DO/

AC 156)/0.

A ( 1 57 ) /O

.

A ( 1 58 ) /O

.

A < 1 59 ) /O

.

AC 1601/0.

AC 161 1/0.

AC 1621/0.

AC 1631/0.

AC 1641/0.

AC 1651/0.

AC 1661/0.

AC 1671/0.

AC 1681/0.

AC 1691/0.

AC 1701/0.

AC 171 1/0.

AC 1721/0.

AC 1731/0.

AC1741/0.

AC 1751/0.

AC 1761/0.

AC 1771/0.

AC 1781/0.

A C 179 )/0.

A C 180 )/0.

A C 181 /A \O

A C 182 )/0.

A C 183 )/0.

A C 184 I/O.

A C 185 )/0.

A C 186 )/0.

A C 187 )/0.

A c ie8 )/0.

A C 189 )/0.

A C 190 )/0.

A C 191 )/0.

A C 192 )/0.

A C 193 >/0.

A C 194 )/0.

A C 195 )/0.

A C 196 )/0.

A C 197 )/0.

A C 198 )/0.

A C 199 )/0.

A C 200 )/0.

A C 201 )/0.

A C 202 )/0.

A C 203 )/ 0.

A C 204 >/0.

A C 205 )/0.

A C 206 )/0.

A C 207 )/0.

A C 208 )/0.

A C 209 )/0.

A C 210 )/0.

A C 2 1

1

)/0.

006504 457 96897 8D0/

00884 6759826363D0/

011168139460131D0/

01 34630478967 18D0/

01 572603047 6024D0/

0179517 1 5 77 5697 DO/

020 1 34823 153530D0/

022270 173808383D0/

0243527025687 10D0/

026 37 7 4697 15 054 DO/

02833967261 4259D0/

030234657072402D0/

03205792835485 IDO/

033805161 837 141D0/

035472213256882D0/

037055 128540240D0/

038550 1 53 1 7861 5D0/

03995374 1132720D0/

041 26 2 56 324 2 62 3 DO/

042473515123653D0/

04 358 372 4 52 9 32 3 DO/

044590558 1 63756D0/

045491 6279274 18D0/

046284796581 31 4D0/

046968 1828 1 621 ODO/

047540 1 657 14B30D0/

047999388596458D0/

04 8 34 4 76 2 23 4 80 2 DO/

04857546744 1503D0/

048690957009 139D0/

001 1 44950003 186D0/

00266353358951 2D0/

004 180313124694D0/

00569092245 1403D0/

0071 92904768 117D0/

00 8683 94 526 9 26 ODO/

010161766041 103D0/

011624114120797D0/

01 306876 159240 IDO/

014493508040509D0/

01 5896 183583725 DO/

017274652056269D0/

01 862681 4208299D0/

019950610878141D0/

021 244026 115782D0/

022505090246332D0/

02373 1B82865930D0/

0249225357641 15D0/

026075235767565D0/

027188227500486DO/

02825981 6057276D0/

029288369583267D0/

03027232 1 759557D0/

031210174188114D0/

0321 00498673487 DO/

03294 1939397645D0/
286 DATA X ( 21

287 DATA X ( 21

288 DATA X ( 214

289 DATA X ( 215

290 DATA X ( 216

291 DATA X ( 217

292 DATA X( 218

293 DATA X ( 219

294 DATA X ( 220

295 DATA X ( 221

296 DATA X ( 222

297 DATA X ( 223

298 DATA X ( 224

299 DATA X ( 225

300 C N = 96

301 DATA X ( 226

302 DATA X ( 227

303 DATA X ( 228

304 DATA X ( 229

305 DATA X ( 230

306 DATA X ( 231

307 DATA X ( 232

308 DATA X ( 233

309 DATA X ( 234

310 DATA X ( 235

311 DATA X ( 236

312 DATA X ( 237

313 DATA X ( 238

314 DATA X ( 239

315 DATA X ( 240

316 DATA X ( 24

317 DATA X ( 242

318 DATA X ( 243

319 DATA X ( 244

320 DATA X ( 245

321 DATA X ( 246

322 DATA X ( 247

323 DATA X ( 248

324 DATA X ( 249

325 DATA X ( 250

326 DATA X ( 251

327 DATA X ( 252

328 DATA X ( 253

329 DATA X ( 254

330 DATA X ( 255

331 DATA X ( 256

332 DATA X ( 257

333 DATA X ( 258

334 DATA X ( 259

335 DATA X ( 260

336 DATA X ( 261

337 DATA X ( 262

338 DATA X ( 263

339 DATA X ( 264

340 DATA X ( 265

341 DATA X ( 266

342 DATA X ( 267

)/0. 5028041 118887B4D0/

)/0. 4686966 151 70544D0/

)/0 . 4 3387537 0831 7 56D0/

>/0.39B393405881969DO/

>/0. 362304753499487D0/

)/0. 32566437074770100/

)/0. 28852805488451 100/

) /0 . 25 095235839 227 2D0/

) /0 . 212994502857666D0/

)/0. 174712291832646D0/

)/0. 136164022809143D0/

)/0. 09740839844 1584D0/

)/0. 058504437152420DO/

)/0. 01951 13B3256793D0/

)/0 . 99 96895038 83230D0/

)/0. 9983643758631 81 DO/

)/0. 99598 1842987209D0/

) /0 . 9 92 5 4 39003 2 37 62D0/

)/0. 988054 126329623D0/

)/0 . 9 8 25 1726 356 30 14 DO/

)/0. 9759391 745851 36D0/

)/0. 96832682846326400/

)/0. 95968829 1448742D0/

)/0. 950032717784437D0/

)/0 . 9 3937 033975 27 5 5D0/

)/0. 9277 124567 223 08D0/

)/0. 9 1507 14231 20898D0/

)/0. 901460635315852D0/

)/0. 88689451740242000/

)/0. 87138850590929600/

)/0. 8549590334 3460 IDO/

)/0. 83762351 12281 87D0/

)/0. 81 94003 1073793 IDO/

)/0. 8003087441 391 40D0/

) / O . 7 8036904 38674 3 3D0/

)/0. 75960234 11 76647D0/

) / O . 7 3803064 374 44 OODO/

)/0. 7 1567681 234 8967D0/

)/0. 69256453664217100/

)/0. 6687 183 1 00439 16D0/

) /O . 644 1 634 037 8496 7D0/

)/0. 6 189258401 254 68D0/

)/0. 593032364777572D0/

)/0. 5665 104 18561 397D0/

)/0. 5393881 08324357D0/

) /O. 51 16941 771 54667D0/

)/0. 483457973920596D0/

)/0. 454709422167743D0/

)/0. 425478988407300D0/

)/0. 395797649828908D0/

)/0. 365696861472313D0/

)/0. 335208522892625D0/

)/0. 304364S4435449600/

)/0.273198B12591049D0/

)/0. 24 17431 561 63840D0/

)/0. 21 00313 10460567D0/

A ( 212 )/0. 03373321490461 IDO/

A<213)/0. 034473 120451 75 3D0/

A( 214 ) /O . 035 16 0529 04 4 74 7 DO/

A ( 215 )/0. 0357943939534 16D0/

A(216)/0. 03637374 990 5835 DO/

A(217)/0. 0368977 14638276D0/

A(218)/0. 037365490238730D0/

A<219)/0. 037776364 362001 DO/

A < 220 )/0. 038 12971 1314 477D0/

A(221 )/0. 038424 99 300 695 9 DO/

A(222)/0. 038 66 175 977 4 07 6 DO/

A ( 223 )/0. 03883965 105905 IDO/

A<224)/0. 038958395962769D0/

A ( 225 ) /O . 03901781 3656306D0/

A ( 226 ) /O . 000796792065552D0/

A(227 )/0. 00 185 3960 78894 6D0/

A( 228 )/0. 00291 073 181 7934D0/

A( 229 )/0. 003964554338444 DO/

A ( 230 ) /O . 00501 4 20 274 2 927 DO/

A( 231 )/0. 006058545504235D0/

A( 232 )/0. 007096470791 153D0/

A ( 233 ) /O . 00812687 6925698D0/

A<234)/0. 009 1 4067 1230783D0/

A ( 235 ) /O .010 1 6 077 053 5008 DO/

A ( 236 ) /O . 01 1 162102099838D0/

A(237)/0. 012151604671088D0/

A ( 238 )/0. 01 3 12822956696 IDO/

A < 239 )/0. 01 409094 17723 14D0/

A< 240) /O. 01 503872 1026994D0/

A ( 241 ) /O . 015970562902562D0/

A ( 242 )/0. 01 688547986424 5D0/

A ( 243 )/0. 01 77825023 16045D0/

A ( 244 )/0. 01866067962741 IDO/

A ( 245 )/0. 019519081 1401 45D0/

A ( 246 )/0. 020356797 154 333D0/

A(247)/0. 02 1172S39892191D0/

A(248)/0. 02 196664 4438744D0/

A ( 249

)

/O . 022737069658329D0/

A ( 250

)

/O . 0234 8 3 39 908 5 926 DO/

A ( 25 1 )/0. 02420484 1792364 DO/

A(252)/0. 02 4 90063 3 22248 3D0/

A( 253 )/0. 02557003600534 9D0/

A( 254 )/0. 02621 2340735672D0/

A C 255 )/0. 02682686672559 IDO/

A C 256

)

/O . 027 412962726029D0/

A ( 257 ) /G . 02797000761 684 8D0/

A ( 258 )/0. 02849741 1065085D0/

A ( 259 ) /O . 02899461 4 150555D0/

A( 260) /O. 02946 1 089958 167D0/

A (261 )/0. 029896344 136328D0/

A(262)/0. 0302999 15420827D0/

A ( 263 )/0. 03067 1376 123669D0/

A ( 264 )/0. 03 101 033258631 3D0/

A ( 265 ) /0 . 031316425596861 DO/

A(266)/0. 031 589330 77 0727 DO/

A ( 267

)

/O . 03182875889441 IDO/

85
343 DATA X(26B)/0. 17809688236761800/# A ( 268 )/0. 032034 45623 1992D0/

344 DATA XC269)/0. 1 459737 1 4654896D0/# AC 269 )/0. 032206204 794 030D0/

345 DATA XC270)/0. 1 1 36958501 1 0665D0/# AC 270 )/0. 03234 382256B575D0/

346 DATA XC271 )/0.081297495464425D0/# AC271 )/0. 032447 1 637 14064D0/

347 DATA X<272)/0. 048812985136049D0/# AC272)/0. 03251611871 3 86 8 DO/

348 DATA X(273)/0. 016276744849602D0/# AC273)/0. 03255061 4492363D0/

349 C

350 C Test N

351 N = NX

352 ALPHA=0 . 5 DO* C BX+AX

)

353 BETArO. 5D0*CBX-AX)

354 IF CN.EQ.l) THEN

355 Z( 1 ) = ALPHA

356 W( 1 ) rBX-AX

357 ELSE IF ( (N.GT. 1 . AND. N. LE . 16) .OR.

N

. EQ.20.0R.N.EQ.24.0R.N.EQ.32.0R.

358 + N. EQ. 40. OR. N. EQ. 48. OR. N. EQ. 64 . OR .N.EQ.80.0R.N.EQ.96) THEN

359 c Set K equal to initial subscript and store results

360 KrKTABCN)

361 M=N/2

362 DO 904 J=1,M

363 JTAB=K-1+J

364 WTEMP=BETA*A( JTAB)

365 DELTArBETA*X( JTAB)

366 Z(J)= ALPHA -DELTA

367 W( J) =WTEMP

368 JPrN+l-J

369 Z( JP)=ALPHA+DELTA

370 W( JP)=WTEMP

371 904 CONTINUE

372 IF ( (N-M-M) . EQ. 0) RETURN

373 Z ( M+ 1 ) = ALFHA

374 JMID=K+M

375 UKM + l 5 =BETA*A( JMID)

376 ELSE

377 ZN = N 378 WRITE ( 1 » 50 ) ZN

379 50 FORMAT( 'OError in GSET . N has the nonperm i ss ab

1

e value'#

380 + 1 PEI 1.3)

381 ENDIF

382 RETURN

383 END

86
1 STITL EUALI -- 61-point Gauss-Kronrod integration

2 FUNCTION EUALI (A, B, FUNC, NS)

3 C Function to evaluate the integral of the function FUNC

4 C between A and B. The number of intervals starts at NSTEP

5 C as defined -in EUICOM and then goes up by a factor of two

6 C until the error in the integral is smaller than EPS.

7 C In this version/ 61-point Gauss-Kronrod integration

8 C using the subroutine DQK61 from the NBS library CMLIB

9 C is done.

10 IMPLICIT DOUBLE PRECISION (A-H.O-Z)

11 C The maximum number of times that the integral is evaluated

12 PARAMETER (IDIUMX=10)

13 EXTERNAL FUNC

14 COMMON /EUICOM/ NSTEP/ EPS

15 NSrNSTEP

16 IDIUrO

17 c 18 c Start of loop

19 * 4 STEP: (B-A)/NS

20 EU=0. DO

21 ERR;0. DO

22 DO 902 1:1/ NS

23 AX : A+ ( 1-1 ) *STEP

24 BX : AX+STEP

25 CALL DQK61 (FUNC/ AX/ BX/ RESULT/ ABSERR / RESABS / RESASC

)

26 EU : EU + RESUL T 27 ERR:ERR+ABSERR

28 902 CONTINUE

29 IF (EPS.LE.O. DO) THEN

30 CD WRITE (1,*) A, B, ERR, DABS(EU*EPS)

31 IF (ERR. LE. DABS(EU*EPS) ) GO TO 6

32 ELSE

33 CD WRITE (1/*) A, B, ERR, EPS

34 IF (ERR.LE.EPS) GO TO 6

35 ENDIF

36 IF (NS. NE. NSTEP) THEN

37 IF (EPS. LE. 0. DO) THEN

38 IF ( DABS ( EU-EUAL 10) . LE . DABS ( EU* EPS ) ) GO TO 6

39 CD WRITE (1,*) A , B, DABS ( EU-EUAL 10 ) , DABS ( EU*EPS

)

40 ELSE

41 IF (DABS(EU-EUALIO) . LE. EPS) GO TO 6

42 CD WRITE (1/*) A, B, DABS (EU-EUALIO) , EPS

43 ENDIF

44 ENDIF

45 EUAL 10; EU

46 NS:2*NS

47 I DI U = I D IU+

48 IF ( IDIU. LE. IDIUMX) GO TO 4 49 NS: — 1 50 6 EUAL I : EU

51 RETURN

52 END

87
12

3

4

5

6

7

Q

9

10

1

1

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

DIU. - K.U.LEUUEN

DIO. - K . U . LEUUEN

SUBROUTINE DQK6 1 C F, A, B, RESULT, ABSERR, RESABS , RESA SC

)

C***BEG I N PROLOGUE DQK61

C***DATE WRITTEN 800101 (YYMMDD)

C***REUISION DATE 810101 (YYMMDD)

C*** CATEGORY NO. D1

C* <K<<KEYWORDS 61-POINT GAUSS-KRONROD RULES

C***AUTHOR PIESSENS, ROBERT, APPL. MATH. & PROGR.

C DE DONCKER, ELISE,APPL. MATH. & PROGR

C D.K. KAHANER ( NBS

)

C**#PURPOSE TO COMPUTE I r INTEGRAL OF F OUER <A,B) WITH ERROR

C ESTIMATE

C J = INTEGRAL OF DABS ( F ) OUER ( A , B

)

c***description

c

c

c

c

cc

ccc

c

c

c

ccc

c

c

c

c

c

c

cc

c

c

c

c

cc

c

c

c

cc

c

cc

c

c

c

c

cc

c

INTEGRATION rule

STANDARD FORTRAN SUBROUTINE

DOUDLE PRECISION UERSION

PARAMETERS

ON ENTRY

F DOUBLE PRECISION

SUBROUTINE DEFINING THE INTEGRAND F(X) AT X: CALL F(X,F(X>). THE ACTUAL NAME FOR F NEEDS TO BE

DECLARED EXTERNAL IN THE CALLING PROGRAM.

DOUBLE PRECISION

LOWER LIMIT OF INTEGRATION

DOUBLE PRECISION • UPPER LIMIT OF INTEGRATION

DOUBLE PRECISION

APPROXIMATION TO THE INTEGRAL I RESULT IS COMPUTED BY APPLYING THE 61-POINT

KRONROD RULE (RESK) OBTAINED BY OPTIMAL ADDITION OF

ABSCISSAE TO THE 30-POINT GAUSS RULE CRESG).

ABSERR - DOUBLE PRECISION

ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,

WHICH SHOULD EQUAL OR EXCEED DABS ( I -RESULT

)

RESABS - DOUBLE PRECISION

APPROXIMATION TO THE INTEGRAL J

ON RETURN

RESULT -

RESASC - DOUBLE PRECISION

APPROXIMATION TO THE INTEGRAL OF DABS ( F- 1 / ( B-A )

)

REFERENCES (NONE)

ROUTINES CALLED D 1 MACH , DABS, DMAX 1 , DMI N 1 , F ( USER-PROUI DED FUNCTION)

END PROLOGUE DQK61

88
58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

DOUBLE PRECISION A, DABSC, ABSERR, B, CENTR, DABS, DHLGTH, DMAX1 , DM INI

* D1MACH, EPMACH, FC, FSUM, FUAL1, FUAL2, FU1 , FU2, HLGTH, RESABS, RESASC-

* RESG, RESK, RESKH, RESULT , UFLOW, WG, WGK, XGK

INTEGER J , JTW » JTWM

DIMENSION FU1 ( 30 ) , FU2 ( 30 ) , XGK< 31 ) , WGK ( 31 ) , WG C 1 5

)

THE ABSCISSAE AND WEIGHTS ARE GIUEN FOR THE

INTERUAL (-1,1). BECAUSE OF SYMMETRY ONLY THE POSITIVE

ABSCISSAE AND THEIR CORRESPONDING WEIGHTS ARE GIUEN.

XGK - ABSCISSAE OF THE 61-POINT KRONROD RULE

XGK ( 2 ) , XGK ( 4 ) ... ABSCISSAE OF THE 30-POINT

GAUSS RULE

XGK ( 1 ) , XGK ( 3 ) ... OPTIMALLY ADDED ABSCISSAE

TO THE 30-POINT GAUSS RULE

WGK - WEIGHTS OF THE 61-POINT KRONROD RULE

WG - WEIGTHS OF THE 30-POINT GAUSS RULE

GAUSS QUADRATURE WEIGHTS AND KRONRON QUADRATURE ABSCISSAE AND WEIGHTS

81 C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETI C BY L. W. FULLERTON

82 C BELL LABS, NOV. 1981 . 83 C

84 DATA WG ( 1 ) 0.0079681924 9616660561 5465883474 674 DO /

85 DATA WG ( 2) / 0.0184664683 1109095914 2302131912 047 DO /

86 DATA WG ( 3) s 0.0287847078 8332336934 9719179611 292 DO /

87 DATA WG ( 4) / 0.0387991925 6962704959 6801936446 348 DO /

88 DATA WG ( 5) s 0.0484026728 3059405290 2938140422 808 DO /

89 DATA WG ( 6) / 0. 057.4931562 1761906648 1721689402 056 DO /

90 DATA WG < 7) / 0.0659742298 8218049512 8128515115 962 DO /

91 DATA WG ( 8) / 0.0737559747 3770520626 8243850022 191 DO /

92 DATA WG ( 9) / 0.0807558952 2942021535 4694938460 530 DO /

93 DATA WG ( 10) s 0.0868997872 0108297980 2387530715 126 DO /

94 DATA WG ( 11 ) / 0. 0921225222 3778612871 7632707087 619 DO /

95 DATA WG ( 12) / 0.0963687371 7464425963 9468626351 810 DO /

96 DATA WG ( 13) s 0.0995934205 8679526706 2780282103 569 DO /

97 DATA WG ( 14) s 0. 1017623897 4840550459 6428952168 554 DO

98 DATA WG ( 15) / 0. 1028526528 9355884034 1285636705 415 DO /

99 C

100 DATA XGK ( 1 ) / 0.9994844100 5049063757 1 3258S5705 811 DO /

101 DATA XGK ( 2) / 0.9968934840 7464954027 1630050918 695 DO /

102 DATA XGK < 3) / 0.9916309968 7040459485 8628366109 486 DO /

103 DATA XGK < 4) s 0. 9836681232 7974720997 0032581605 663 DO /

104 DATA XGK < 5) / 0.9731163225 0112626837 4693868423 707 DO

105 DATA XGK ( 6) / 0.9600218649 6830751221 6871025581 798 DO

106 DATA XGK ( 7) / 0. 9443744447 4855997941 5831324037 439 DO /

107 DATA XGK ( 8) s 0.9262000474 2927432587 9324277080 474 DO /

108 DATA XGK ( 9) / 0.9055733076 9990779854 6522558925 958 DO /

109 DATA XGK ( 10) / 0. 8825605357 9205268154 31 16462530 226 DO /

110 DATA XGK ( 11 ) / 0.8572052335 4606109895 8658510658 944 DO /

111 DATA XGK ( 12) / 0.8295657623 8276839744 28981 19732 502 DO /

112 DATA XGK < 13) s 0. 7997278358 2183908301 3668942322 683 DO /

1 1 3 DATA XGK ( 14 ) / 0. 7677774321 0482619491 797734097503 DO /

114 DATA XGK ( 15) / 0.7337900624 5322680472 6171131369 528 DO /

89
115 DATA XGK ( 16) / 0.6970504947 9331579693 2292388026 640 DO /

116 DATA XGK ( 17) / 0. 6600610641 2662696137 0053668149 271 DO /

117 DATA XGK < 18) 0.6205261029 8924286114 0477556431 109 DO /

118 DATA XGK ( 19) / 0.5793452358 2636169175 6024932172 540 DO /

119 DATA XGK ( 20) / 0.5366241481 4201989926 4169793311 073 DO /

120 DATA XGK ( 21

)

/ 0.4924804678 6177857499 3693061207 709 DO /

121 DATA XGK ( 22) / 0.4470337695 3808917678 0609900322 054 DO /

122 DATA XGK < 23) / 0.4004012548 3039439253 5476211542 661 DO /

123 DATA XGK < 24) / 0.3527047255 3087011347 1037207089 374 DO /

124 DATA XGK ( 25) / 0.3040732022 7362507737 2677107199 257 DO /

125 DATA XGK ( 26) / 0.2546369261 6788984643 9805129817 805 DO /

126 DATA XGK ( 27) / 0.2045251166 8230909143 8957671002 025 DO /

127 DATA XGK ( 28) / 0. 1530699136 0858354696 3794672743 256 DO /

120 DATA XGK ( 29) / 0. 1028069379 6673703014 7096751318 001 DO /

129 DATA XGK ( 30) / 0.0514718425 5531769583 3025213166 723 DO /

130 DATA XGK ( 31 ) / 0.0000000000 0000000000 0000000000 000 DO /

131 C

132 DATA WGK ( 1) / 0.0013890136 9867700762 4551591226 760 DO /

133 DATA WGK ( 2) / 0.0038904611 2709988405 1267201844 516 DO /

134 DATA WGK < 3) / 0.0066307039 1593129217 3319826369 750 DO /

135 DATA WGK ( 4) / 0.0092732796 5951776342 8441146892 024 DO /

136 DATA WGK ( 5) / 0.0118230152 5349634174 2232898853 251 DO /

137 DATA WGK < 6) / 0.0143697295 0704580481 2451432443 580 DO /

138 DATA WGK ( 7) / 0.0169208891 8905327262 7572289420 322 DO /

139 DATA WGK ( 8) / 0.0194141411 9394238117 3408951050 128 DO /

140 DATA WGK < 9) / 0.0210280358 2160919229 7167485738 339 DO /

141 DATA WGK ( 10) / 0.0241911620 7808060136 5686370725 232 DO /

142 DATA WGK ( 11 ) / 0.0265099540 8233310161 0601709335 075 DO /

143 DATA WGK ( 12) / 0.0287540487 6504129284 3978785354 334 DO /

144 DATA WGK ( 13) / 0.0309072575 6238776247 2884252943 092 DO /

145 DATA WGK < 14) / 0.0329814470 5740372603 1814191016 854 DO /

146 DATA WGK ( 15) / 0.0349793380 2806002413 7499670731 468 DO /

147 DATA WGK ( 16) / 0.0368823646 5182122922 3911065617 136 DO /

148 DATA WGK C 17) / 0.0386789456 2472759295 0348651532 281 DO /

149 DATA WGK < 18) / 0.0403745389 515359591

1

1995279752 468 DO /

150 DATA WGK ( 19) / 0.0419698102 1516424614 7147541205 970 DO /

151 DATA WGK ( 20) / 0. 0434525397 0135606931 6831728117 073 DO /

152 DATA WGK < 21 ) / 0.0448148001 3316266319 2355551616 723 DO /

153 DATA WGK ( 22) / 0.0460592382 7100698811 6271735559 374 DO /

154 DATA WGK ( 23) / 0.0471855465 6929915394 5261470181 099 DO /

155 DATA WGK ( 24) / 0.0401858617 5708712914 0779492298 305 DO /

156 DATA WGK ( 25) / 0.0490554345 5502977888 7528165367 238 DO /

157 DATA WGK c 26) / 0.0497956834 2707420635 7811569379 942 DO /

158 DATA WGK ( 27) / 0. 0504059214 0278234684 0893085653 585 DO /

159 DATA WGK < 28) / 0.0508817958 9874960649 2297473049 805 DO /

160 DATA WGK ( 29) / 0.0512215478 4925877217 0656282604 944 DO /

161 DATA WGK ( 30) / 0.0514261285 3745902593 3862879215 781 DO /

162 DATA WGK ( 31 ) / 0.0514947294 2945156755 8340433647 099 DO /

163 C

164 C LIST OF MAJOR VARIABLES

165 C

166 C

167 C CENTR MID POINT OF THE INTERVAL

168 C HLGTH1 - HALF-LENGTH OF THE INTERVAL

169 c DABSC ABSCISSA

170 c FVAL* FUNCTION VALUE

171 c RESG - RESULT OF THE 30-POINT GAUSS RULE

90
172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

RESK - RESULT OF THE 61-POINT KRONROD RULE

RESKH - APPROXIMATION TO THE MEAN VALUE OF F OVER < A , B ) , I.E. TO I/(B-A)

MACHINE DEPENDENT CONSTANTS

EPMACH IS THE LARGEST RELATIUE SPACING.

UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.

EPMACH = D1MACH ( 4

)

UFLOW = DlMACH(l)

C

CENTR = O. 5D+00*<B+A)

HLGTH = 0. 5D+00*(B-A)

DHLGTH = DABS ( HLGTH

)

C

C COMPUTE THE 61-POINT KRONROD APPROXIMATION TO THE

C INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.

C

C***FIRST EXECUTABLE STATEMENT DQK61

193 RESG = O.OD+OO

194 CALL F < CENTR, FC ) 195 RESK r WGK(31)*FC

1 196 RESABS = DABS ( RESK

)

197 DO 10 J=l, 15

198 JTW = J*2

199 DABSC = HLGTH*XGK( JTW)

200 CALL F ( CENTR-DABSC, FUAL1

)

201 CALL F ( CENTR+ DABSC, FOAL 2

)

202 FOl(JTW) = F0AL1

203 F02 (JTW ) = F0AL2

204 FSUM = FOAL 1 +F0AL2

205 RESG = RESG+WG(J) *FSUM

206 RESK = RESK+WGK(JTW) *FSUM

207 RESABS = RESABS+WGK (JTW)*(DABS( F0AL1 ) +DABS ( F0AL2 )

)

208 10 CONTINUE

209 DO 15 J=l, 15

210 JTWM1 = J*2-l

211 DABSC = HLGTH*XGK( JTWM1

)

212 CALL F (CENTR-DABSC, F0AL1

)

213 CALL F (CENTR+ DABSC, FOAL 2

)

214 FOl(JTWMl) = FOAL

215 F02CJTWM1) = F0AL2

216 FSUM = FOAL 1+F0AL2

217 RESK : RESK+WGKCJTWMl ) *FSUM

218 RESABS = RESABS+WGK C JTWM1 )*(DABS(F0AL1 )+DABS(F0AL2)

219 15 CONTINUE

220 RESKH = RESK*0. 5D+00

221 RESASC = WGK (31 ) *DABS ( FC- RESKH

)

222 DO 20 J = 1 , 30

223 RESASC = RESASC+WGK ( J ) * ( DABS ( FO 1 ( J ) -RESKH )+DABSCF02

224 20 CONTINUE

225 RESULT r RESK*HLGTH

226 RESABS : RESABS* DHLGTH

227 RESASC : RESASC *DHLGTH

, 228 ABSERR = DABS ( ( RESK -RESG ) *HLGTH

)

91

I
229

230

231

232

233

234

I F ( RESASC . NE . 0 . 0D+00 . AND .ABSERR.NE.O. OD+OO)

* ABSERR = RESASC*DMIN1 (0. 1D+01 * ( 0. 2D+03*ABSERR/RESASC ) ** 1 . 5D+00)

IF(RESABS.GT.UFL0U/(0.5D + 02*EPriACH) > ABSERR = DMAX1

* C (EPMACH*0. 5D+02)*RESABS, ABSERR)

RETURN

END

92
DOUBLE PRECISION FUNCTION D1MACH( I

)

-

12

3

4

5

6

7

e

9

10

n

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

53

54

55

56

57

C

C

C

C

C

Cc

c

cc

c

cc

c

cc

c

c

c

cc

DOUBLE-PRECISION MACHINE CONSTANTS

D1MACHC 1) = B**<EMIN-1), THE SMALLEST POSITIUE MAGNITUDE.

D1MACH ( 2) = B**EMAX* ( 1 - B**(-T>), THE LARGEST MAGNITUDE.

D1MACH ( 3) = B**(-T), THE SMALLEST RELATIUE SPACING.

D1MACH ( 4) = B** ( 1 -T ) / THE LARGEST RELATIUE SPACING.

D1MACH ( 5) = LOG 1 0(B)

TO ALTER THIS FUNCTION FOR A PARTICULAR ENUIRONMENT

THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACT I UATED BY

REMOUING THE C FROM COLUMN 1. WHERE POSSIBLE, OCTAL OR HEXADECIMAL CONSTANTS HAUE BEEN USED

TO SPECIFY THE CONSTANTS EXACTLY WHICH HAS IN SOME CASES

REQUIRED THE USE OF EQUIUALENT INTEGER ARRAYS.

INTEGER SMALL ( 4

)

INTEGER LARGE ( 4

INTEGER RIGHTC4)

INTEGER DIUER ( 4

)

INTEGER L0G10C4

)

DOUBLE PRECISION DMACH ( 5

)

EQUI UALENCE C DMACH < 1 ), SMALL < 1 )

)

EQUIUALENCE < DMACH ( 2 ), LARGE ( 1 )

)

EQUI UALENCE ( DMACH ( 3 ) , RI GHT ( 1 )

)

EQUIUALENCE ( DMACH ( 4 ), DIUER ( 1 )

)

EQUIUALENCE ( DMACH ( 5 ) , LOG 1 0 ( 1 )

)

MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.

38 C

39 C DATA SMALL ( 1 ) / ZC00800000

40 C DATA SMALL ( 2

)

/ 2000000000 /

41 C

42 C DATA LARGE < 1

)

/ ZDFFFFFFFF /

43 C DATA LARGE C 2 ) X ZFFFFFFFFF /

44 C

45 C DATA RIGHT ( 1

)

/ ZCC5800000 /

46 C DATA RIGHT (2 ) / ZOOOOOOOOO /

47 c 48 c DATA DIUERC 1

)

/ ZCC6800000 /

49 c DATA DIUER(2> ZOOOOOOOOO /

50 c 51 c DATA L0G1

( 1

)

/ ZD00E730E7 /

52 c DATA L0G1

( 2 ) / ZC77800DC0 /

MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.

DATA SMALL ( 1 ) 01771000000000000 /

DATA SMALL ( 2 ) / 00000000000000000 /

93
50 . C

59 C DATA LARGE ( 1

)

/ 00777777777777777 /

60 c DATA LARGE (2) / 00007777777777777 /

61 c 62 c DATA RIGHT(1

)

/ 01461000000000000 /

63 c DATA RIGHT ( 2

)

/ 00000000000000000 /

64 c 65 c DATA DIUER < 1

)

/ 01^51000000000000 /

66 c DATA DIUERC2) / ooooooooooooooooo /

67 c 60 c DATA L0G1 0 ( 1

)

/ 01157163034761674 /

69 c DATA L0G1 0 ( 2

)

/ 00006677466732724 /

70 c 71 c MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700

72 c 73 c DATA SMALL ( 1

)

/ 01771000000000000 /

74 c DATA SMALL ( 2

)

/ 07770000000000000 /

75 c 76 c DATA LARGE ( 1

)

/ 00777777777777777 /

77 c DATA LARGE ( 2

)

/ 07777777777777777 /

70 c 79 c DATA RIGHT (

1

) 01461000000000000 /

00 c DATA R I GHT (

2

) / OOOOOOOOOOOOOOOOO /

01 c 02 c DATA D I UER ( 1 ) / 01451000000000000 /

03 c DATA DIUERC2) / ooooooooooooooooo /

04 c 05 c DATA L0G1 0(1 ) 01157163034761674 /

06 c DATA L0G1

0

( 2 ) / 00006677466732724 /

87 c 80 c MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES

89 c 90 c DATA SMALL ( 1 ) / 006 04 0000 00 0000 00 000 B /

91 c DATA SMALL ( 2 ) / OOOOOOOOOOOOOOOOOOOOB /

92 c 93 c DATA LARGE ( 1 ) / 37767777777777777777B /

94 c DATA LARGE ( 2 ) / 371 67777777777777777 B /

95 c 96 c DATA R I GHT ( 1

)

/ 15604 000000000 0000 00 B /

97 c DATA R I GHT ( 2 ) / 15 000000000000000000 B /

98 c 99 c DATA D I UER ( 1 ) / 156 1 4 0000 00 000 000 000 B /

100 c DATA D I UER ( 2 ) / 1 50 1 OOOOOOOOOOOOOOOOB /

101 c 102 c DATA L0G10 ( 1 ) / 17164642023241 175717B /

103 c DATA L0G1

0

( 2 ) / 16367571421742254654B /

104 c 105 c MACHINE CONSTANTS FOR THE CRAY 1 106 c 107 c DATA SMALL (1 ) / 200004 000000 00 OOOOOOOB /

108 c DATA SMALL (2) / OOOOOOOOOOOOOOOOOOOOB /

109 c 110 c DATA LARGE ( 1 ) / 577777777777777 777777 B /

111 c DATA LARGE (

2

) / 000007 777 77 7 77 7 77 7 77 7

B

/

112 c 113 c DATA RIGHTC

1

) / 37721 4 00 000 000 000 000 OB /

,114 c DATA R I GHT ( 2 ) / OOOOOOOOOOOOOOOOOOOOOB /

94
115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

138

139

140

141

142

143

144

145

146

147

155

156

157

158

159

DATA DIUER(l) / 377224000000000000000B /

DATA DIUER ( 2 ) / OOOOOOOOOOOOOOOOOOOOOB /

DATA LOGlO(l) / 377774642023241 175717B /

DATA L0G1 0(2) / 00000757 1 42 1 742254654B /

MACHINE CONSTANTS TOR THE DATA GENERAL ECLIPSE S/200

NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD

STATIC DMACH ( 5

)

DATA SMALL/20K, 3*0/ , LARGE/7777 7K, 3*177777K/

DATA RIGHT/31420K, 3*0/ , D I UER/32020K, 3*0/

DATA L0G10/40423K/ 42023K, 50237K, 74776K/

MACHINE CONSTANTS FOR THE HARRIS 220

133 C DATA SMALL ( 1 ) , SMALL ( 2

)

/ '20000000, '00000201 /

134 C DATA LARGE ( 1 ) , LARGE (2) / '37777777, '37777577 /

135 C DATA R I GHT (1), RIGHTO) / '20000000, '00000333 /

136 C DATA D I UER ( 1 ) , DIUERO) / '20000000, '00000334 /

137 C DATA L0G10 < 1 ) , LOG 1 0(2) / '23210115, ' 10237777 /

MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.

DATA SMALL( 1 ) , SMALL ( 2 ) / 0402400000000, 0000000000000 /

DATA LARGE( 1 ) , L ARGE ( 2 ) / 0376777777777, 0777777777777 /

DATA RIGHT ( 1 ) , RIGHT (2) / 0604400000000, OOOOOOOOOOOCO /

DATA DIUER( 1 ) , DIUERO) / 0606400000000, 0000000000000 /

DATA L0G10( 1 ) ,L0G10(2) / 0776464202324, 0117571775714 /

MACHINE CONSTANTS FOR THE HP 2100

148 C 3 WORD DOUBLE PRECISION OPTION WITH FTN4

149 C

150 C DATA SMALL ( 1 ), SMALL ( 2 )

,

SMALL ( 3

)

/ 40000B, 0, 1 /

151 C DATA LARGE (1 ), LARGE ( 2 )

,

LARGE ( 3

)

/ 77777B

,

177777B, 177776B /

152 C DATA R I GHT ( 1 ) , RIGHT O )

,

RIGHTO) / 40000B, 0, 265B /

153 C DATA D I UER ( 1 )

,

DIUERO), DIUERO) / 40000B

,

0, 276B /

154 c DATA L0G1 0 ( 1 )

,

LOG 10(2), LOG 10(3) / 46420B, 46502B

,

77777B /

MACHINE CONSTANTS FOR THE HP 2100

4 WORD DOUBLE PRECISION OPTION WITH F'TN4

160 C DATA SMALL (1 ), SMALL ( 2

)

/ 40000B, 0 /

161 C DATA SMALL C3)

,

SMALL ( 4

)

/ 0, 1 /

162 c DATA LARGE ( 1 )

,

LARGE C 2

)

/ 77777B

,

177777B /

163 c DATA LARGE ( 3 )

,

LARGE ( 4 ) / 177777B, 177776B /

164 c DATA R I GHT ( 1 ) , RIGHT ( 2

)

/ 40000B

,

0 /

165 c DATA RIGHTO) , RIGHT ( 4 ) / 0, 225B /

166 c DATA DIUERC 1 )

,

D IUER ( 2

)

/ 40000B

,

0 /

167 c DATA DIUERO) , DIUER ( 4

)

/ 0, 227B /

168 c DATA L0G1 0(1 )

,

LOG 10(2) / 46420B

,

46502B /

169 c DATA LOGIOO) , LOG 10(4) / 76747B, 176377B /

170 c 171 c MACH I NE CONSTANTS FOR THE IBM 360/370 SERIES

95
172

173

174

173

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

21 1 212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

non THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86 AND

THE PERKIN ELMER (INTERDATA) 7/32.

C

C

C

C

C

C

C

C

C

C

C

C

Cc

cc

c

c

c

DATA SMALL (1

DATA LARGE (

1

DATA RIGHTd

DATA DIUER (

1

DATA L0G1 0 (

) , SMALL ( 2

)

) , LARGE ( 2

)

) ,RIGHT(2)

) , DIUER (2

)

) , LOG 1 0(2)

/ Z00100000,

/ Z7FFFFFFF,

/ Z33 100000

/

/ Z34100000,

/ Z4 1 1 3441 3 ,

ZOOOOOOOO /

ZFFFFFFFF /

ZOOOOOOOO /

ZOOOOOOOO /

Z509F79FF /

MACHINE CONSTANTS FOR THE PDP-10 ( KA PROCESSOR).

DATA SMALL (1),SMALL(2) / "033400000000,

DATA LARGE(1),LARGE(2) / "377777777777,

DATA RIGHT(1),RI GHT ( 2 ) / "113400000000,

DATA DIUERd), DIUER(2) / "114400000000,

DATA L0G10( 1 ) ,L0G10(2) / "177464202324,

"000000000000 /

"344777777777 /

"000000000000 /

"000000000000 /

"144117571776 /

MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).

DATA SMALL (1), SMALL (2) / "000400000000,

DATA LARGE (1),LARGE(2) / "377777777777,

DATA RIGHT( 1 ) , RIGHT(2) / "103400000000,

DATA DIUERd ),DIUER(2) / "104400000000,

DATA LOGlOd ) ,L0G10(2) / "177464202324,

"000000000000 /

"377777777777 /

"000000000000 /

"000000000000 /

"476747767461 /

MACHINE CONSTANTS FOR PDP-11 FORTRAN'S SUPPORTING

32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).

c DATA SMALL (1

)

, SMALL (2) / 8388608, 0 /

c DATA LARGEd) , LARGE ( 2

)

/ 2147483647, -1 /

c DATA RIGHT ( 1

)

, R I GHT ( 2

)

/ 612368384, 0 /

c DATA DIUERd ) , D I UER ( 2

)

/ 620756992, 0 /

c

c DATA LOGlOd ) , L0G1 0 ( 2

)

/ 1067065498, -2063872008 /

Smt c DATA SMALL (1

)

, SMALL ( 2

)

/ 000040000000, 000000000000 /

c DATA LARGEd ) , LARGE ( 2

)

/ 017777777777, 037777777777 /

c DATA RIGHT (1

)

, RIGHT (2

)

/ 004440000000, 000000000000 / c DATA DIUERd ) , D I UER ( 2

)

/ 004500000000, 000000000000 /

c

r* DATA LOGlOd ) , L0G1 0 ( 2

)

/ 007746420232, 020476747770 /

c MACHINE CONSTANTS FOR PDP -11 FORTRAN'S SUPPORTING

c

c 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).

VC DATA SMALL ( 1

)

, SMALL ( 2

)

/ 128, 0 /

Cr DATA SMALL ( 3

)

, SMALL ( 4

)

/ 0, 0 /

LC DATA LARGEd ) , LARGE ( 2

)

/ 32767, -1 /

C

r DATA LARGE ( 3

)

, LARGE ( 4

)

/ -1, -1 /

VC DATA RIGHTd ) , R I GHT ( 2

)

/ 9344, 0 /

C DATA RIGHT (3) , R I GHT ( 4

)

/ 0, 0 /

c DATA DIUERd ) , DIUER(2) / 9472, 0 /

c

r DATA DIUER ( 3

)

, DIUER(4) / 0, 0 /

c DATA LOGlOd ) , LOG 1 0(2) / 16282, 8346 /

c DATA L0G1 0 ( 3

)

, L0G10 ( 4

)

/ -31493, -12296 /

96
230 C DATA SMALL ( 1 ) , SMALL ( 2

)

/ 0000200. 0000000 /

231 C DATA SMALL ( 3 ) » SMALL ( 4

)

0000000. 0000000 /

232 C

233 C DATA LARGE (1 ) » LARGE ( 2

)

0077777. 0177777 /

234 C DATA LARGE(3), LARGE ( 4

)

/ 0177777. 0177777

235 C

236 C DATA RIGHT(1),RIGHT(2) / 0022200. 0000000 /

237 c DATA RIGHTO) ,RIGHT(4) / 0000000. 0000000 /

238 c 239 c DATA DIUER d ) , DIUER ( 2

)

/ 0022400. 0000000 /

240 c DATA DIUER ( 3 ) / DIUER ( 4

)

/ 0000000. ooooooo /

241 c 242 c DATA LOGlOd ) » L0G1 0 ( 2 ) / 0037632, 0020232

243 c DATA L0G10O)»L0G10(4) / 0102373, 0147770 /

244 c 245 c MACHINE CONSTANTS FOR THE UNIUAC 1100 SERIES. FTN COMPILER

246 c 247 c DATA SMALL < 1 ) * SMALL ( 2

)

/ 0000040000000, 0000000000000 /

248 c DATA LARGE ( 1 ) , LARGE ( 2

)

/ 0377777777777, 0777777777777 /

249 c DATA RIGHT d ) , RIGHT (2

)

/ 0170540000000, 0000000000000 /

250 c DATA DIUERd ) , DIUER ( 2

)

/ 0170640000000, 0000000000000 /

251 c DATA L0G10< 1 ) , LOG 10 (2

)

/ 0177746420232, 0411757177572 /

252 c 253 c 254 c MACHINE CONSTANTS FOR THE UNIUAC 1100 SERIES. FOR COMPILER

255 c 256 c DATA SMALL( 1 ) ,SMALL(2) 0000040000000, 0000000000000 /

257 c DATA LARGEd ) , LARGE ( 2 ) / 0377777777777, 0777777777777 /

258 c DATA RIGHT d ) .RIGHTO) / 0170540000000, 0000000000000 /

259 c DATA DIUERd ) , DIUERO) / 0170640000000, 0000000000000 /

260 c DATA LOGlOd ) ,L0G10(2) / 0177746420232, 0411757177572 /

261 c 262 c MACHINE CONSTANTS FOR UAX -11/780

263 c (EXPRESSED IN INTEGER AND OCTAL)

264 c 265 c DATA SMALL(l).SMALLO) / 128, 0 /

266 c DATA LARGEd ) » LARGE ( 2 ) / -32769, -1 /

267 c DATA RIGHT ( 1

)

» R I GHT ( 2

)

/ 9344, 0

268 c DATA DIUER ( 1 ) , DIUERO) / 9472, 0 /

269 c DATA LOGlOd ) .LOGIOO) / 546979738, -805665541 /

270 c 271 c DATA SMALL ( 1 )

SMALL ( 2 ) / 000000000200, 000000000000 /

272 c DATA LARGEd ) » LARGE ( 2 ) / 037777677777, 037777777777 /

273 c DATA RIGHTd ), RIGHTO) / 000000022200, 000000000000 /

274 c DATA DIUERd ) , DIUERO) / 000000022400, 000000000000 /

275 c DATA LOGlOd ). LOGIOO) / 004046437632, 031776502373 /

276 c 277 c IF (I .LT. 1 .OR. I .GT . 5)

278 c 1 CALL SETERR( 24HD1MACH - I OUT OF BOUNDS, 24, 1, 2)

280

281

282

283

D1MACH = DMACH(I)

RETURN

END

97
12

3

4

5

6

7

e

9

10

1

1

12

13

14

15

16

17

IB

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

n

n

n

o

STITL PWAFFI — Initialize form-factor table for C-12

SUBROUTINE PWAFFI

IMPLICIT DOUBLE PRECISION (A-H,0-Z)

C The size of the arrays in the common PWAFF

PARAMETER (MNPT=200)

C The number of points in the interpolation table

PARAMETER <NTBM1=179)

C The number of points + 1 PARAMETER (NTB=NTBM1+1

)

DOUBLE PRECISION XC 1 ( 1 00 ) , XC2 ( 79 ) , FC 1 C 1 00 ) , FC2 < 79 ) , X ( MNPT

)

DOUBLE PRECISION F ( NTBM1 ), DF < NTBM1 ), A ( MNPT ), B < MNPT ), C ( MNPT

)

DOUBLE PRECISION D ( MNPT ), AP ( MNPT ), BP ( MNPT ), CP ( MNPT ), DP C MNPT

)

DOUBLE PRECISION XC ( NTBM 1 ) , FC C NTBM 1

)

COMMON /CONST/ EM, EM2 , AM , HB ARC , ALPH, P IE , ALPHA , RTD , DTR

COMMON /PWAFF/ NPT, X, A, B, C, D, AP, BP, CP, DP

COMMON /CORREC/IWTILD, WTIL1 , WTIL2, I SOFT, S0FT1

EGUIUALENCE ( FC ( 1 ) , FC1 ( 1 ) ) , ( FC ( 101 ) , FC2 ( 1 )

)

EQUIUALENCE ( XC ( 1 ) , XC 1 ( 1 ) ) , ( XC < 1 01 ) , XC2( 1 )

)

Table of q and corresponding form factor from DWBA code. In principle this table could be read in rather than defined in

data statements.

DATA XC1/5. 3068D-02, 1 . 06 1 3D-01 , 1 . 5918D-01, 2. 1 222D-0 1 , 2 . 6523D-01

+ 3. 1 822D-01 , 3. 71 18D-01, 4. 2410D-01, 4. 7697D-01 , 5. 2980D-01

,

+ 5. 8258D-01, 6.3530D-01,6. 8795D-01 , 7 . 4054D-01 , 7.9305D-01

,

+ 8 . 4549D-0 1 »8.97B4D-01,9.5010D-01, 1.0023D+00, 1.0543D+00,

+ 1 . 1 063D+00, 1 . 1582D + 00, 1 . 2099D+00, 1 . 2615D+ 00, 1 . 3130D + 00,

+ 1 . 3644D+00, 1 . 4156D+00, 1 . 4667D+00, 1 . 5177D+00, 1 . 5685D+00,

+ 1 . 6191 D + 00, 1 . 6696D+00, 1 . 7199D + 00, 1 . 7701 D+00, 1 . 8201 D+00,

+ 1 . 8699D+ 00, 1 . 9195D+00, 1 . 9690D+00, 2. 0182D+ 00, 2. 0673D+00,

+ 2.1 162D + 00, 2. 1648D+00, 2 . 2 1 33D + 00 , 2 . 26 1 5D+ 00 , 2 . 3096D + 00

,

+ 2. 3574D+ 00, 2. 4050D + 00, 2. 4524D+00, 2. 4 995D+00, 2 . 54 64 D + 00,

+ 2. 5931 D+00, 2. 6395D + 00, 2. 6856D + 00, 2. 7316D+00, 2. 7772D + 00,

+ 2. 8226D+00, 2. 8678 D + 00, 2. 9127D + 00, 2. 9573D+00, 3. 0017D + 00,

+ 3. 0457 D + 00, 3. 0895D + 00, 3. 1331 D + 00, 3. 1763D + 00, 3. 2193D + 00,

+ 3. 2619D + 00, 3 . 3043D+ 00, 3. 34 64D + 00, 3. 3881 D + 00, 3. 4296D+ 00,

+ 3. 4708D + 00, 3. 51 16D+00, 3. 5522D+ 00, 3. 5924D + 00, 3. 6324 D + 00,

+ 3. 6720D+ 00, 3. 71 13D+00, 3 . 7502D + 00, 3. 7889D+ 00, 3. 8272D + 00,

+ 3. 8652D + 00, 3. 9029D+ 00, 3. 9402D + 00, 3. 9772D+00, 4 . 0138D+ 00,

+ 4. 0501 D+00, 4. 0861 D+00, 4. 1218D+00,4. 1570D+00,4. 1920D+00,

+ 4. 2266D+00, 4 . 2608D+00, 4 . 2947D+00, 4. 3282D+00, 4. 36 14 D+00,

+ 4 . 3942D+00, 4 . 4267D+00, 4 . 4588D+00, 4. 4906D+00, 4 . 5219D+00/

DATA XC2/4. 5529D+00,4. 5836D400, 4. 6139D+00, 4 . 6438D+00, 4 . 6733D+00

+ 4. 7025D+00, 4 . 7313D+00, 4 . 7598D+00, 4 . 7878D+00, 4 . 8155D+00,

+ 4. 8428D+00, 4 . 8698D+00, 4 . B963D+00, 4. 9225D+00, 4. 9483D+00,

+ 4. 9737D+00, 4 . 9987D + 00, 5 . 0234D+00, 5. 0477D+00, 5 . 0715D+ 00,

+ 5. 0950D+00, 5. 1182D+00,5. 1409D+00,5. 1632D+00,5. 1852D+00,

+ 5. 2068D+00, 5. 2280D+00, 5. 24 88D + 00, 5. 2692D+ 00, 5. 2892D + 00,

+ 5. 3088D+00, 5. 3280D+00, 5. 3469D+00, 5 . 3653D+00, 5. 3834D+00,

+ 5. 4011 D+00, 5 . 4183D+00, 5. 4352D+00, 5. 4517D+00, 5. 4678D+00,

+ 5. 4935D+00, 5 . 4988D+00, 5 . 51 37D+00, 5. 5283D+00, 5. 5424 D+00,

+ 5. 5561 D + 00, 5. 5695 D + 00, 5. 5824D + 00, 5. 5949D+00, 5. 6071 D + 00,

+ 5. 61B8D + 00, 5. 6302D + 00, 5. 641 1 D + 00, 5. 6517D+00, 5. 6618D + 00,

+ 5. 6716D+00, 5. 6810D+00, 5. 6899D+00, 5. 6985D+00, 5 . 7067D+00,

+ 5. 7 145D+00, 5.72 18D+00, 5 . 7288D+00, 5. 7354 D+00, 5. 7416D+00,

+ 5. 7474 D+00, 5 . 7527D + 00 , 5 . 7577D + 00, 5. 7623D+ 00, 5 . 7665D+00,

+ 5. 7703D+00, 5. 7737D + 00,5. 7767D + 00, 5. 7793D+00, 5. 7815D + 00,
58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

o

oo

+ 5. 7833D+00, 5. 7847D+00, 5 . 7857D+00 , 5 . 7862D+00/

DATA FC1/1 . 0000D + 00 > 9 . 8087D-01, 9. 5239D-01, 9. 1556D-01 / 8. 697 5D-01

,

+ 8. 1646D-01, 7. 5736D-01 , 6. 9420D-01, 6. 2873D-01, 5. 6261 D-01,

+ 4. 9734D-01, 4.3427D-01, 3. 7450D-01 , 3. 1890D-01 , 2 . 6807 D-01 , + 2 . 2240D-0 1 , 1 . 8203D-01 , 1 . 4692D-01, 1 . 1688D-01, 9. 1572D-02,

+ 7. 0602D-02, 5. 3507D-02,3. 9805D-02, 2. 9015D-02, 2. 0675D-02,

+ 1. 4357D-02, 9. 6756D-03, 6. 2920D-03, 3 . 916 1 D-03, 2 . 3048D-03,

+ 1 . 2588D-03, 6. 1895D-04, 2 . 6065D-04, 8. 91 12D-05, 3 . 4238D-05,

+ 4 . 5899D-05, 8.9708D-05, 1 . 4332D-04, 1 . 9333D-04, 2. 3265D-04,

+ 2. 5852D-04, 2. 7093D-04, 2. 71 43D-04 , 2 . 6234D-04, 2 . 4618D-04

,

+ 2. 2533D-04, 2. 0186D-04, 1 . 7744D-04, 1 . 5334 D-04, 1 . 3047 D-04,

+ 1 . 0941 D-04, 9. 0497D-05,7. 3864D-05 , 5 . 951 2D-05 , 4. 7338D-05,

+ 3. 7177D-05, 2.8825D-05/2. 2064D-05, 1 . 6672D-05, 1.2435D-05,

+ 9. 1554D-06, 6.6544D-06, 4 . 7753D-06, 3 . 3842D-06, 2 . 3691 D-06

,

+ 1. 6387D-06, 1 . 1202D-06»7. 56 6 2D -07 , 5. 0471D-07, 3. 3205D-07 » + 2. 1493D-07, 1 . 3630D-07, 8. 4 1 39D-08 , 5 . 0045D-08 , 2 . 825 1 D-08

,

+ 1 . 4815D—08, 7.0487D-09, 3. 1135D-09, 1 . 7327D-09, 2 . 0021 D-09

,

+ 3. 2652D-09, 5.0351D-09, 6. 9478D-09, 8 . 7363D-09, 1 . 0215D-08,

+ 1 . 1272D-08, 1 . 1 855D-08 » 1 . 1969D-08, 1 . 1660D-08, 1 . 1001D-08,

+ 1 . 0085D-08, 9. 0060D-09, 7 . 8529D-09, 6. 7016D-09, 5 . 6109D-09,

+ 4. 6210D-09, 3. 7547D-09, 3. 0203D-09, 2. 4145D-09, 1 . 9268D-09/

DATA FC2/1 . 5422D-09, 1 . 2443D-09, 1 . 0168D-09, 8 . 4 509D-1 0 , 7 . 1668D-10,

+ 6. 2124D-10, 5.5047D-10/4. 9786D-10, 4. 5831D-10, 4. 2782D-10,

+ 4. 033 ID- 10, 3. 8236D-10, 3. 6319D-10, 3. 4448D-10, 3. 2538D-10,

+ 3. 054 ID- 10, 2. 8440D-10, 2. 6246D-10, 2. 3986D-10, 2. 1698D-10,

+ 1 . 9427D-10, 1 . 7217D— 10, 1.5108D-10, 1.3133D-10, 1. 1315D-10,

+ 9.6708D-11, 8. 206 ID- 11, 6. 9205D-11, 5.8070D-11, 4.B544D-11,

+ 4. 0484D-11, 3. 3732D-1 1 ,2. 8123D-11 , 2. 3498D-1 1, 1 . 9707D-1 1

,

+ 1.6614D-11, 1 .4098D-11, 1 .2055D-11, 1 . 0398D- 1 1 , 9 . 05 1 4 D- 1 2

,

+ 7. 9547D-12, 7.0577D-12, 6. 31 95D- 1 2 , 5 . 707 1 D- 1 2 , 5 . 1943D-12,

+ 4. 760 ID- 12, 4. 3877 D- 12, 4 . 0639D-12, 3. 7782D-12, 3. 5222D-12,

+ 3. 2896D-12, 3. 0756D-12, 2. 8765D-12, 2. 6900D-12, 2. 5143D-12,

+ 2. 3433D-12, 2. 1 91 1 D-12, 2 . 0422D-12, 1 . 9009D-12, 1 . 7673D-12,

+ 1 . 641 2D- 12, 1 .5229D-12, 1 . 4126D-12, 1 . 3105D-12, 1 . 2164D-12,

+ 1. 1 302D-1 2, 1.051 3D- 12, 9. 792 ID- 13, 9. 1 362D- 1 3, 8 . 5462D- 1 3

,

+ 8. 026 ID- 13, 7 .5814 D- 13, 7 . 21 58D- 1 3 , 6 . 9264D- 1 3 , 6 . 6982D- 1 3

,

+ 6. 4928D-13, 6. 2152D-13, 5. 5475D-13, 3. 2798D-13/

Convert table from form factor versus q to form factor

times 2**2 versus q**2

DO 10 1=1, NTBM1

X ( I ) =XC( I ) * *2

F ( I )=FC( I ) *36 . DO

IF ( IWTILD . NE . 0 ) F C I ) =F < I ) * ( WT I L 1 +WT IL2*DL0G ( X ( I ) /EM2 )

)

10 CONTINUE

C

C Find the cubic coefficients for spline interpolation of

C the table. A good reference for spline interpolation is

C Elementary numerical analysis by S.D. Conte and Carl deBoor.

A ( 1 ) =36 . DO

B( 1 ) =0. DO

CC 1 ) =0. DO

D( 1 ) =0. DO

X(NTB) =1 . D74

A ( NTB ) = 0 . DO

BCNTB)=0. DO

C(NTB) =0. DO

99
115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

DCNTB) =0. DO

C Calculate the derivatives of the form factor ( DF ) according

C to the spline method at each of the tabulated points.

P=0. DO

QrO. DO

H2r 1 .D0/(X(2)-X( 1)

)

DO 908 I z 2 / NTBM

1

H1=H2

H12=H1*H1

H2= 1 . DO/(X( 1+1 ) -X ( I )

)

H22zH2*H2

Rz2. D0*(H1+H2)-P

B ( I ) :H2/R

DF< I ) z-3. DO* ( F ( 1-1 )*H12+F( I ) * ( H22-H12 ) -F ( I +1 )*H22)

DF(I)=(DF(I)-Q)/R

PrB( I )*H2

0= DF ( I )*H2

908 CONTINUE

C Now back substitution

PrO. DO

DO 909 I r 1 / NTBM1

JrNTBMl+1-I

DFC J)zDF( J)-P*B( J)

PrDFC J)

909 CONTINUE

C Using the spline approximation values for the derivatives

C at the points in the table/ calculate the coefficients

C for a cubic which exactly fits these derivatives and the

C values of the function at the end points of each interval.

DO 20 I z2/ NTBM1

Hr 1 . DO/(X( I )-X( 1-1 )

)

H2r H*H

AlrDF (I)-2. DO*F ( I ) *H

Cl r F ( I )-Al*X( I ) A2zDF ( 1-1 ) +2 . DO*F ( I-i )*H

C2=F( 1-1 ) -A2*X ( 1-1

)

D(I)r(Al+A2)*H2

C( I >r (Cl -2. D0*A1*X(I-1 )+C2-2. D0*A2*X< I ) )*H2

B ( I ) z ( A 1 *X ( I - 1 ) **2-2 . D0*C1 *X( 1-1 ) +A2*X( I )**2-

+ 2. D0*C2*X( I ) )*H2

A(I)r(Cl*X(I-l) **2+C2*X(

I

)**2)*H2

20 CONTINUE

C For this model/ we assume that Wl (the transverse form factor)

C is 0.

DO 30 Irl,NTB

AP( I ) rO. DO

BP( I ) rO. DO

CP( I ) rO. DO

DPC I ) rO. DO

30 CONTINUE

NPTzNTB

RETURN

END

100
12

3

45

6

7

8

9

10

11

12

13

14

15

16

17

18

STITL F0RM1 — Transverse form factor by interpolation

FUNCTION F0RMKQ2)

C Size of arrays in PWAFF

PARAMETER (MNPT=200)

IMPLICIT DOUBLE PRECISION (A-H/O-Z)

DOUBLE PRECISION X ( MNPT) / A ( MNPT ) / B ( MNPT ) » C ( MNPT )<D( MNPT

)

DOUBLE PRECISION AP ( MNPT ), BP (MNPT ), CP < MNPT ), DP ( MNPT

)

COMMON /PWAFF NPT/X»A»B#C»D»AP»BP# CP» DP

C CALL STARTT(l)

0Q=-Q2

DO 10 1=1/ NPT—

IF (QQ.LT.X(I)) GO TO 20

10 CONTINUE

I =NPT

20 F0RM1 = (

(

DP ( I ) *QQ+CP ( I ) )*QQ+BP( I ) )*QQ+AP( I

)

C CALL STOPT(l)

RETURN

END

101
12

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

*TITL F0Rr^2 — Longitudinal form factor by interpolation

FUNCTION F0RM2C02)

C Size of arrays in PWAFF

PARAMETER (MNPT=200)

IMPLICIT DOUBLE PRECISION (A-H,0-Z)

. DOUBLE PRECISION X ( MNPT ) , A ( MNPT ) , B ( MNPT ) , C ( MNPT ) , D < MNPT

)

DOUBLE PRECISION AP ( MNPT ), BP C MNPT ), CP ( MNPT ), DP ( MNPT

)

COMMON /"PWAFF / NPT*X*A,B/C»D* AP / BP / CP / DP

C CALL STARTT(l)

QQz-QZ

DO 10 1=1* NPT —

1

IF ( QQ . LT . X ( I ) ) GO TO 20

10 CONTINUE

I =NPT

20 F0RM2= ( (DC I)*QQ + C< I ) >*QQ + B< X ) >#QQ + A( I

)

C CALL STOPT(l)

RETURN

END

Manuscrit raqu le 7 novembre 1983

102
NBS-1 14A (rev. 2 »eo) U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPC'’'r NO.

NBSIR 83-2788

2. Performing Organ. Report No. 3. Publication

DECEMBER 1983

4. TITLE AND SUBTITLE

PIECE-WISE ANALYTIC EVALUATION OF THE RADIATIVE TAIL

FROM ELASTIC AND INELASTIC ELECTRON SCATTERING

5. AUTHOR(S)

Leonard C. Maximon and Steven E. Williamson

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.

8. Type of Report & Period Covered

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

10. SUPPLEMENTARY NOTES

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

bi bl iography or literature survey, mention it here)

If document includes a significant

We consider here the calculation of the radiative tail from the elastic

peak in medium and high energy electron scattering as well as from a discrete

inelastic level of the recoiling nucleus. We examine the method generally

used for this calculation, viz., a numerical integration of the differential

cross section over the angles of the unobserved photon, and discuss the

difficulties inherent in this numerical integration due to the sharp peaking

of the integrand. We present an alternative method for calculating the

radiative tail, in which the region of integration is divided into an arbitrary number of subintervals, the structure functions are fitted by cubic

spline functions in each subinterval, and the integrations are then performed

analytically in closed form. This method has the advantages of greatly

increased accuracy and a reduction of the computation time by a factor which

can vary between 10 and 10 3

, depending on the kinematics.

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

Key words: electron scattering cross sections; elastic electron scattering;

elastic radiative tail; inelastic electron scattering; inelastic

radiative tail; intermediate and high energy electron scattering.

13. AVAILABILITY

[)Q Unlimited

1 | For Official Distribution. Do Not Release to NTIS

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

20402.

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

14. NO. OF

PRINTED PAGES

107

15. Price

$ i3.CC

USCOMM-DC 6043-P80