TECHNICAL REPORT CERC-87-16

US Army Corps COMBINED REFLECTION AND DIFFRACTION of Engineers BY A VERTICAL WEDGE by H. S. Chen

Coastal Engineering Research Center

DEPARTMENT OF THE ARMY

Waterways Experiment Station, Corps of Engineers PO Box 631, Vicksburg, Mississippi 39180-0631

8a,

December 1987 Final Report

\ NY Lif, iI Hl eC ns. Le \ Approved For Public Release, Distribution Unlimited it | \

|

|

|

|

Prepared for DEPARTMENT OF THE ARMY US Army Corps of Engineers Washington, DC 20314-1000

Under Waves at Entrances Work Unit 31673

(.S. Rees Bos ea ee ieee. Kea CEES.

TECHNICAL REPORT CERC-87-16

US Army Corps COMBINED REFLECTION AND DIFFRACTION of Engineers BY A VERTICAL WEDGE

by H. S. Chen

Coastal Engineering Research Center

DEPARTMENT OF THE ARMY Waterways Experiment Station, Corps of Engineers PO Box 631, Vicksburg, Mississippi 39180-0631

Be: I 1. 1 ; fe b; i : \ _ Woods Hole Qeeanosra sits } Instit

December 1987

EA Oss NE Ie Wi Tf one is

; a Final Report sai) We) BREAN WATE ik o//// |i\ \ is erty | fray \ Approved For Public Release, Distribution Unlimited Zi fay yey yi a IS) | / | \ a i if |i | \ Ba / } / | | off | |

Prepared for DEPARTMENT OF THE ARMY US Army Corps of Engineers Washington, DC 20314-1000

Under Waves at Entrances Work Unit 31673

Destroy this report when no longer needed. Do not return it to the originator.

The findings in this report are not to be construed as an Official Department of the Army position unless so designated by other authorized documents.

This program is furnished by the Government and is accepted and used by the recipient with the express understanding that the United States Government makes no warranties, expressed or implied, concerning the accuracy, completeness, reliability, usability, or suitability for any par- ticular purpose of the information and data contained in this program or furnished in connection therewith, and the United States shall be under no liability whatsoever to any person by reason of any use made thereof. The program belongs to the Government. Therefore, the recipient further agrees not to assert any proprietary rights therein or to represent this program to anyone as other than a Government program.

The contents of this report are not to be used for

advertising, publication, or promotional purposes.

Citation of trade names does not constitute an

official endorsement or approval of the use of such commercial products.

Unclassified SECURITY CLASSIFICATION OF THIS PAGE

Form Approved

REPORT DOCUMENTATION PAGE OMB No. 0704-0188

Exp. Date Jun 30, 1986

Ja. REPORT SECURITY CLASSIFICATION 1b. RESTRICTIVE MARKINGS Unclassified

2a. SECURITY CLASSIFICATION AUTHORITY 3. DISTRIBUTION / AVAILABILITY OF REPORT

Approved for public release; distribution

2b. DECLASSIFICATION / DOWNGRADING SCHEDULE i unlimited.

4. PERFORMING ORGANIZATION REPORT NUMBER(S) 5. MONITORING ORGANIZATION REPORT NUMBER(S)

Technical Report CERC-87-16

6a. NAME OF PERFORMING ORGANIZATION 6b. OFFICE SYMBOL 7a. NAME OF MONITORING ORGANIZATION USAEWES, Coastal Engineering (If applicable)

Research Center

6c. ADDRESS (City, State, and ZIP Code) 7b. ADDRESS (City, State, and ZIP Code) PO Box 631 Vicksburg, MS 39180-0631

8a. NAME OF FUNDING/SPONSORING 8b. OFFICE SYMBOL 9. PROCUREMENT INSTRUMENT IDENTIFICATION NUMBER ORGANIZATION (If applicable)

US Army Corps of Engineers 8c. ADDRESS (City, State, and ZIP Code) 10. SOURCE OF FUNDING NUMBERS

PROGRAM PROJECT TASK WORK UNIT ELEMENT NO. | NO NO. ACCESSION NO Washington, DC 20314-1000 31673

11. TITLE (Include Security Classification)

Combined Reflection and Diffraction by a Vertical Wedge

12. Chesenn ouanene) hen,

13a. TYPE OF REPORT 13b. TIME COVERED. 14. DATE OF REPORT (Year, Month, Day) ]15. PAGE COUNT Final report FROM December 1987 75

16. ee ee NOTATION

val rom National Technical Information Service, 5285 Port Royal Road, Springfield,

VA 22161.

18. SUBJECT TERMS (Continue on reverse if necessary and identify by block number) Breakwaters (LC) Me os Numerical analysis (LC) Rasa a | ee | Water waves 19. ABSTRACT (Continue on reverse if necessary and identify by block number)

This report presents an analytical solution for the combined wave reflection and diffraction by a vertical wedge, of arbitrary wedge angle, subject to the excitation of a plane simple harmonic wave train coming from a far field. Results of two special cases are calculated: one a thin semi-infinite breakwater and the other a wedge of 90-deg angle. The results are presented in amplification factor diagrams. A subroutine WEDGE is also documented in the appendix. The subroutine can be used to calculate the case of arbitrary wedge angle.

20. DISTRIBUTION / AVAILABILITY OF ABSTRACT 21. ABSTRACT SECURITY CLASSIFICATION Ed) UNCLASSIFIED/UNLIMITED (L] SAME AS RPT. OO otic users Unclassified

22a. NAME OF RESPONSIBLE INDIVIDUAL 22b. TELEPHONE (Include Area Code) | 22c. OFFICE SYMBOL

DD FORM 1473, 84 mar 83 APR edition may be used until exhausted. SECURITY CLASSIFICATION OF THIS PAGE All other editions are obsolete linoilacenened

WUC MMMM

0 0301 O041Leb?

SECURITY CLASSIFICATION OF THIS PAGE

eee SECURITY CLASSIFICATION OF THIS PAGE

PREFACE

The work in this report was authorized by the Office, Chief of Engineers (OCE), Coastal Engineering Functional Area of Civil Works Research and Development, under Waves at Entrances Work Unit 31673, Harbor Entrances and Coastal Channels Program, at the Coastal Engineering Research Center (CERC) of the US Army Engineer Waterways Experiment Station (WES). Messrs. John H. Lockhart, Jr., and John G. Housley were OCE Technical Monitors. Dr. Charles L. Vincent is CERC Program Manager.

This report was prepared by Dr. H. S. Chen, Coastal Oceanography Branch (CR-O), Research Division (CR). Work was performed under direct supervision of Dr. Edward F. Thompson, Chief, CR-O, and Mr. H. Lee Butler, Chief, CR; and under general supervision of Dr. James R. Houston and Mr. Charles C.

Calhoun, Jr., Chief and Assistant Chief, CERC, respectively.

This study was initiated after a discussion by Dr. Vincent and the au- thor on the possibility of implementing a scheme to redistribute wave energy behind islands in numerical models. The author acknowledges and appreciates the review and comments provided by Drs. Edward F. Thompson and Norman W. Scheffner. This report was edited by Ms. Shirley A. J. Hanshaw, Information Products Division, Information Technology Laboratory, WES.

Commander and Director of WES is COL Dwayne G. Lee, CE, and Technical

Director is Dr. Robert W. Whalin.

CONTENTS

IINIMVNGIG 50 0000000000000000000000000000000000000000000000000000000000000

ARS Telas TON EROWDIUG ILO 5 6 9000000000000000000000000000009000000000000000 PART II: BOUNDARY VALUE PROBLEM... .cccccccccccrcccrcrcccccccsccsscosees

MEVMaMAeICaL IO rMMLEVELOMNs o4660000000000000000000000000000000000000 AMAISAEHOAIL SOIMEROMs 690000000000000000000000000000000000000000000 IK Speeiwad, CASG@Scaodccood0d0 00 DDD DDO CODD DOODDD DODD ODDO DDDDD0000

PART III: CALCULATION AND RESULTS... 2... ccccccccvcccceccervccsvcccsees

Vertical Wedge of O-Deg Wedge Angle.....ccccccccccscccrcecrccccrcecs Vertical Wedge of 90-Deg Wedge Angle......cccesccccvcccccccrcccce

PART IV: (GOMGMWISIMONS 6 G50 000000000000000000000000000000000000000000000 INIDIDIDRVAMGINS 6 9600000000006000000000000000000000000000000000060000006000000 JNO (AQ SW RNOWALION, MINE 6 6 0000000000000000000000000000000000050000 JNPIEANDIDG 138 INOWMEMON 5 6900000000000000000000000000000000000000000000000

COMBINED REFLECTION AND DIFFRACTION BY A VERTICAL WEDGE

PART I: INTRODUCTION

1. The boundary value problem of linear wave reflection and diffraction by a vertical wedge of arbitrary wedge angle has been well formulated and presented by Stoker (1957) among many other investigators. The technique to obtain an analytical solution for the problem is also depicted in the cited book. However, analytical solutions are not available for the problem, except for the special case of wave diffraction by a thin semi-infinite breakwater, that is, a wedge with wedge angle equal to zero.

2. The solution of the thin semi-infinite breakwater was presented in the dimensionless diffraction diagrams by Wiegel (1962). The diagrams have been especially useful in preliminary engineering design and have been in- cluded in the Shore Protection Manual (SPM) (1984). Although equally useful, the combined reflection and diffraction diagrams are not available, perhaps because of the complexity of the diagrams which makes them difficult to create without using modern high-speed computers for computation and graphing.

3. The objectives of the present study are (a) to obtain an analytical solution for the combined wave reflection and diffraction by a vertical wedge of arbitrary wedge angle subject to excitation of a plane simple harmonic wave train coming from infinity and (b) to provide the combined reflection and dif- fraction diagrams. The diagrams included in this report have two cases: one for a thin semi-infinite breakwater and the other for a 90-deg vertical wedge. Subroutine WEDGE for computing the combined reflection and diffraction by a vertical wedge of arbitrary wedge angle is also documented in the report (Appendix A).

PART II: BOUNDARY VALUE PROBLEM Mathematical Formulation

4. In this study our primary interest is the wave reflection and diffraction by a vertical wedge of arbitrary wedge angle in a constant water depth h * subject to the excitation of monochromatic incident waves of infi- nitesimal amplitude coming from infinity. Let (r,0,z) be cylindrical coordi- nates, with z= 0 representing the undisturbed water free surface and upward direction representing the positive z-axis. The tip of the wedge is chosen to be the origin of the coordinates and two rigid walls of the wedge to coincide with 6=0 and 6= 85 » respectively, as illustrated in Figure 1. Carte- sian coordinates (x,y,z), corresponding to the cylindrical coordinates, are

also occasionally used and shown in the same figure. Therefore, the wedge

zZ

Figure 1. A vertical wedge of arbitrary wedge angle

* For convenience, symbols and abbreviations are listed in the notation (Appendix B).

angle is 27 - on » and the water region is defined by oF 2620 and 0222-h.

5. The velocity field for the wave reflection and diffraction in an ideal fluid can be represented by the velocity potential function 4(r,®,z,t) which must satisfy the Laplace equation, where t is the temporal coordinate. We assume that the waves are sinusoidal in time with radian frequency w . Water depth is constant, and the bottom is rigid and impermeable. Therefore, the vertical and temporal components of the velocity potential function, which follow from separation of variables, can be factored out and the velocity po-

tential written as

mn cosh k(z + h) iwt O(r,8,z,t) = Ay cosh kh o(r,8)e (1) where Ay = -iga,/w i= yal g = gravitational acceleration aS incident wave amplitude

k = wave number

¢ = horizontal component of the velocity potential function

6. Substituting Equation 1 into the Laplace equation and using both the kinematic and dynamic boundary conditions at the free surface, the Laplace equation is then reduced to the Helmholtz equation which is written in polar

coordinates as follows:

ES aot gine ai gh kro =0 (2)

where k must be a real number and satisfy the dispersion relationship 2 W = gk tanh kh (3)

7. The free surface displacement NN from the mean water level z= 0 can be obtained from linear wave theory and is represented as 1 ao

= Oe iwt n(r,8,t) = a 8G aj o(r, Oe (4)

8. Thus only the horizontal part of the velocity potential function 9 is needed to be determined as a solution of Equation 2 in the water region on > 620 , with the following boundary conditions at the rigid and imperme-

able walls of the wedge:

ad _ 0

ona ae @ = 0 ernal 85 (5)

9. A condition at infinity is also required to ensure a unique solu- tion. The classic approach is to use the Sommerfeld radiation condition at infinity which states that the scattered wave 0 must behave like a cylin-

drical outgoing progressing wave at infinity such that '

a lim vr (ste + sn4,) = 0 (6)

The total wave represented by 9 is the linear superposition of an incident wave %; » a reflected wave from the the ®=0 wall of the wedge o and

the scattered wave Os from the tip of the wedge.

9= 6,4 6,4 9, (7) Equation 6 can be satisfied if eo tkr od ~ at ro (8) = ise

10. The incident wave coming from a large distance from the tip of the wedge is assumed to be a plane progressive wave of amplitude ay and incident angle a to the x-axis as given by

=f etkr cos (6-a)

(9)

Consequently, the perfectly reflected wave from the y = 0 wall of the wedge

is

a etkr cos (6+0) (10)

Thus the boundary value problem (in which the governing equation is Equa- tion 2, the boundary condition is Equation 5, and the radiation condition is

Equation 6) is completely formulated.

Analytical Solution

11. Analytical solution to the problem formulated in the preceding section is obtained by following the solution technique by Stoker (1957). To obtain the solution, the water region is divided into three subregions--I, II, and III--by the incident wave ray passing through the tip of the wedge and the reflected wave ray reflected away from the tip of the wedge, as shown in Fig- ure 2. Obviously, the total wave in subregion I is the sum of the incident, reflected, and scattered waves; the total wave in subregion II, where the reflected wave does not exist, is the sum of the incident and scattered waves; and the total wave in subregion III, where the incident and reflected waves

have been shaded out, is only the scattered wave. For certain combinations of

y

Tt

Figure 2. Three subregions and the wedge

the wedge angle and incident wave angle, subregions II and III may not exist

at all. In general, the solution function can be written as

$= ¢(r,8) + 6 (r,8) (11) where 5 + om T-a>O>0 o (F298) = o; Tta>O6@>7T-a (12) 0 oe > @ > ta ©

The equation reveals that oF is the sum of the incident and reflected waves 5 and $. and is a known function. The scattered wave oe is the only un- known function to be determined in the problem. Nevertheless, the total wave ¢ instead of the scattered wave 0. is the desired solution to be obtained in this study.

12. The solution for the total wave $ is pursued. The finite cosine

transform of 9 , denoted by $ » is introduced by the formula

VT

o(kr,n) = i o(kr,®) cos ut (13)

where n= 0,1,2,... are integers, and v is related to the wedge angle as

defined by

@ = vt (14)

Applying the finite cosine transform and using the boundary condition in Equa-

tion 5, Equation 2 becomes

etree +r = + (IED - (2) ¢=0 (15)

Equation 15 is a form of the Bessel equation for which general solutions are Ja sy 6k) and Y/vike) ;

respectively. Since Ya/vik®) are singular at the origin, the solution is

the Bessel functions of the first and second kinds,

chosen to be

o(kr,n) = a Jy 6k) (16)

where a, are constants to be determined. 13. Taking the finite cosine transform of Equation 11 and using Equa-

tion 16, we have

VT VOT

a Jy (kX) - i 6 cos = (17) 0 0

~o- n Q ie} n Is Q @D i]

or

(kr) - co (18)

2 |

Son ss

Then applying the operation lim... vr(9/dr + ik) to both sides of Equa-

tion 18, and using the Sommerfeld radiation condition (Equation 6) we have

VT F 0 , ng Lim (Fe is ix) a. Tes) © i 6, cos * ae| = 0 (19) (

14. Equation 19 can be asymptotically evaluated to determine ale Firstly, the first term involving the Bessel function is evaluated. The func- tion J py Sk) at reo behaves asymptotically (Abramowitz and Stegun 1964)

as follows:

Jay (ke) ~ Vo cos (ex - = - z) (20)

Hence, we have

Lim "(Fe a, th) py : E: ei (kr-nn/2vtn/4) (21)

Secondly, the second term involving the integral of oF is evaluated. The asymptotic behavior of the integral over 9% = (0,v™) and at large distance r>© can be found by the method of stationary phase. The integral, after sub-

stituting oP from Equations 9, 10, and 12, can be written as

vt TO. e ie ° 6 i d cos ng -f E: Sees=e)) + epliee cos(842) | Coste fo) v av) 0 0 TO. i C= () + gue eerie) cos = do (22) T-a

In the integrals, there are three points of stationary phase at 9 = 4% and 6=7 +a. If the same argument as that of Stoker (1957) is followed, of the three contributions only the first one 9 =a furnishes a nonvanishing con- tribution for r*@ when the operator Y“r(9/dr + ik) is applied to it. The physical significance of this statement is that only the incident wave is

effective in determining the cosine coefficients of the solution. Therefore,

VT lin “e(% + ax) f 9, cos ne d6 ~ 2V2Tk cos = sare) (23) 0

Substituting Equations 21 and 23 into Equation 19, we obtain the unknown co-

efficients a

na eint/2v

anes ZECOS Ra (24) 15. Since the solution ¢ in the cosine series expression is ee ee 2 = nO 6(r,8) = = B(r,0) += JY F(xsn) cos & (25) n=1

10

the solution is obtained by substituting Equations 16 and 24 into Equation 25

as follows:

inn/2v a. ng

o(r,9) = = JQ (kr) + 2 » e Jay ke) cos - cos —— (26)

n=1

Equation 26 is the solution for the combined wave reflection and diffraction by a vertical wedge of arbitrary wedge angle and is considered to be extended from the solution by Stoker (1957) who only solved the problem of a thin semi- infinite breakwater. The solution in Equation 26 and the one by Stoker are not only in nonclosed form but also in terms of Bessel functions. It seems that the calculations of the solutions are very difficult without using a mod- ern high-speed computer. This is probably the reason why Stoker arrived at his solution expressed in the same cosine series but did not use it to calcu- late the result. Instead, he further transformed the expression into a very complex integral form for further approximation in calculating the result.

16. Notably, the solution at the origin point is obtained by simply

substituting r= 0 into Equation 26 to arrive at 2 $(0,8) == (27)

Therefore, wave response at the origin point depends only on the wedge angle

and does not depend on the incident wave angle.

Two Special Cases

17. The solutions for two special cases are used to verify Equation 26: one for the case of a thin semi-infinite breakwater and the other for the case of an infinite wall extending from x = -~ tom,

18. The vertical wedge should reduce to a thin semi-infinite breakwater as the wedge angle reduces to 0 deg. Therefore, solution of the combined wave reflection and diffraction by a thin semi-infinite breakwater is obtained by substituting v=2 (that is, oR = 27) into Equation 26 which then becomes

11

OGe,@)) = JQ (kr) + 2 » einm/4 5 (kr) cos = cos x (28)

n=l

n/2

Equation 28 is precisely the same one obtained by Stoker (1957).

19. The vertical wedge should also become an infinite wall extending from x = to © with the water occupying only the half plane of y20 as the wedge angle increases to 180 deg. In this situation the scattered wave is absent from the solution, and the total wave is only the sum of the incident

and reflected waves as follows:

GeO) = eikr cos (6-0) + eikr cos (8+0) ag

After expansion of the exponential functions in terms of Bessel functions

(Abramowitz and Stegun 1964), Equation 29 becomes

foe)

$(r,8) = 2 Jo (kr) ap 2 » i"J (kr) cos n@ cos n@ (30)

n=1

Equation 30 is the same equation reduced from Equation 26 by substituting

v= 1 into it.

12

PART III: CALCULATION AND RESULTS

20. Results of the combined reflection and diffraction by a wedge of arbitrary wedge angle can be calculated from Equation 26. Since the solution is not only in terms of Bessel functions but also in a nonclosed form, the computer program WEDGE is therefore written to calculate the solution.

21. In the program the subroutine BESJ for calculating Bessel function of fractional or integer order was used. The subroutine was originally writ- ten by Amos, Daniel, and Weston in 1975 (Morris 1984) and is collected in the Naval Surface Weapons Center Library of Mathematics Subroutines (Morris 1984).

22. In the calculation the summation of the infinite terms in Equa- tion 26 was carried out to the term which is preceded by eight successive terms of the absolute value of the Bessel function, all equal to or less than Oman The solution has a truncation error less than Omen and it is of the order of one.

23. In this study, results of the combined wave reflection and diffrac- tion for the wedge are calculated for two cases: one for a vertical wedge of

O-deg wedge angle and the other for a vertical wedge of 90-deg wedge angle.

Vertical Wedge of 0-Deg Wedge Angle

24. When the wedge angle is equal to zero, the wedge is actually a thin semi-infinite breakwater extending from x =0 to ». Figure 3 shows the thin semi-infinite breakwater along with the polar coordinates. In this case the diffraction results for various incident wave angles in the water region from 6 =m to 2m and r/dA < 10 , where A is the incident wave length, have already been presented by Wiegel (1962) and are shown in the SPM, Vol- ume I (1984). The present results combine reflection and diffraction effects and cover the water region from 6 =0 to 2m and r/dA < 10. Therefore, the present results for this particular case can be considered to be a comple- mentary and extended version to the ones in the SPM.

25. In this study wave response was calculated at 1,460 grid points intersected by r/\ = 0.5,(0.5),10 , which means that the values of r/) are from 0.5 to 10.0 with each value increment being 0.5. Hereafter, all similar

expressions are to be interpreted in the same way (e.g., 96 = 0,(7/36),2m for

13

an foo r

og ¢ es 7 S f\ va Ohm, <7 2 2 ~., - 2 ells elena te oy Bac eet : Pantckety nee Bae ot a7 . FoR, ne ee Oe. Uh aoe . ee . Bee a a pa er . Since ae Do 7 ae - : : ~ os oa Oi 4 . i Sete Se, x A a iad * ~2orn- een ra ~ os, + is s Deane ll Ly re me a ares f es ti Poet ' os ce + v a - 5 £ ~ 7 x we a ¢ ae 7 ea . i 5 Mo sie ne A 4 et oh Deen eFes oL So he en 24 a z ce . aa Q ~~ . os me s f fi he v Ch ese ' i ' mG me : x 4 rf ca SBS EN leet - "a i in ~~ a: x os Peas at ca a ao a ng ge 32 eee a -oe , Pe ae} ae f TP Lo SS eit t fo. 5 f uf Ae i : a ~¥L he" % a = % f ¥e vw ca a 3 ae a ae ae aa ny a: + f Oe : ~ Reon! x 5 f t ie oe 38 ea me ' oN Se oy < \ 0, ? ! - ' ' ys Ot - ? 5 f ae “Sup LOSS Sarl OS Diath 8 co Me t t i PO a el ee 7 Boe x POR, é s ;- Meo ween Oh Hh, on oo y x ope oe : t as SOR iietens he tee Se) Set, 4 aoodi 0 t t ~-. 7 ? oT > ~ . at : au w- we Un “EP Sc : HRS age Sad aia ataleteseed | Fa hEN Lyn ci¥ 4 Data T Fe Bly eee sigs piace ea CRUE ay OI ek ar LAAT UT Mga eee math (fer Eh I 1 4 1 u QPoe Ha f ened we Rpatrct wot Bet % oae 4 e i . : I --2. Dey, SHH =e Beoo? : ' I : I i I (Ped) Poae SAP eeree yy I \ 1 I I 4 I I g 1 ul t L l H Srv enced ae 4 1 I I ! t ! I bosedhosctbassd bossa ee SS ey oe Sa I 1 I c 4 et ? 1 1 1 , | s i ! TE ee esa eae ! ! I I ! | 7 x 3 -s. t ! t ! Bee See ete aye a ae 4 5 : y H i : : 4 8 eee ant) ak MBS Seg ene yom te en Cot caer, : i ! I 4 4 Seer hk x my oe ohten SDR Bos Loe cy. ty é 1 4 5 4 Brae cy : ng wen 0 Base To a Sane F) Pict i 4 I ! —W¥ * ue ee Oe EE pista Sarg Poy hit We Sef]. F noo 4 4 4 rehe aoe A £ : i} & ee Pe OS OF) F ! tp fe 4 35 ~ . . Oa . Ps £ t fot. 4 oe > oo . as w = » 4 3 x ~ f——1 + é <2 t f ? Ree Ae - eo ee hd wang ne) ee ee aie ae are a! 5 x nd x ~ T B . . a ws 4 Yee “% x a} a eos : Roo Ae, of Tata ? £ a ee of x - Sg 7+ oy Ruane FSA ft x # 5 aaa 4 “< a 8 % 1 ~ 2 ‘. 8 me XQ oN : : : oi BF ~2e : Be O, ae Sates = ! Pa /0 : es ae oe E 7 Poe ? ¢ x + mS Se tea ae oT ee f f +a * » ean SS Pars 1 H a a aoe ge ra f 4 ~ ~ ~ 7 Tou . 1 Poa x r gz a x OS MELO =e . eo ' ca A G D & = =— ~- fo. . SS ~ a . ~te 1 2 oh r x > = canon S H 0 Pete s aot Z ny oN OG a5 8 . mor Ne ra ge ¢ S, Ss “st ~~ ' - ous “y ~~ ~ Pee ee Bis ra Sa bins 2 oe OS 5 t 4 On a a 7, “6 © 5 fh c Cras pee eae ~ one Soa rice te aia Ke ere ree ~ 1 0 ie o 4 . aoe ieee ee Woes 3, a Racor oo coe ee --4--7- ~ oe Poe ' Seca : +7

eee IS oy egy a

Figure 3. Thin semi-infinite breakwater and polar coordinates

the incident wave angle a = 0,(mn/12),7) . The wave response at the origin

point is obtained substituting v = 2 into Equation 27, as follows:

o(0,6) = 1 (31)

Those calculated values were used to interpret the value for each non-

overlapping pixel of size O.lr/A by O.lr/X in the area within the 10r/A

radius from the origin. A diagram was then constructed by patching those

The wave response diagrams for each incident

pixels over the entire area. Notably, the values in the dia-

wave angle are shown in Figures 4 through 15.

grams constitute the amplification factor which is defined as the ratio of the

total wave height to the incident wave height. Therefore, in subregions II

14

89 16 .28 46 68 .88 86 .16 38

= = = © ®@® 8& 8&8 &Y &

Figure 4. Amplification factor diagram for the thin semi-infinite breakwater for incident wave angle = 0 deg

and III (as defined in Figure 2) where the reflected wave is absent, the am- plification factor is essentially the diffraction coefficient as defined in the SPM.

26. Figures 4 through 15 reveal that the amplification factors in sub- region I change very rapidly between 0 and 2.35 over the subregion, and the diagram patterns become very complex because of the interesting superposition of the incident, reflected, and scattered waves. (In the legend of Figures 4 through 15, the width of the pixel is one incident wave length, and the values are amplification factors.) Such patterns would be very difficult to con- struct without using a high-speed computer and computer graphics. In sub- regions II and III, the amplification factors change smoothly from 1.15

roughly along the reflected wave ray reflected from the origin point to nearly

15

hh eet SS BS S

Amplification factor diagram for the thin semi-infinite breakwater for incident wave angle

- 88 20 35 -85 -15 -45 #5 -85 .35

Mh = se ee SF SS &

Amplification factor diagram for the thin semi-infinite breakwater for incident wave angle

-68 25 35 -85 5 UG: -45 7S -8S ooh

Figure 7.

Amplification factor diagram for the thin semi-infinite breakwater for incident wave angle = 45 deg

18

i ee ee

-88 -29 33 -83 -15 -45 #5 -85 -35

woe hw =— se = FSF FS SF

Amplification factor diagram for the thin semi-infinite breakwater for incident wave angle

a ):) 25 35 85 1S 4S 7S BS 35

Mh == =a = FG FS &

Amplification factor diagram for the thin semi-infinite breakwater for incident wave angle

68 8 £35 85 -15 45 #5 -85 “35

ww = = = ©& SF FS ®&

Amplification factor diagram for the thin semi-infinite breakwater for incident wave angle

-68 .29 -39 -85 -15 45 75 -65 35

me = = = fF fF = @&

Amplification factor diagram for the thin semi-infinite breakwater for incident wave angle

66 28 9 85 15 45 orf 65 .35

Ny hes = =| @ O&O SF @&

Amplification factor diagram for the thin semi-infinite breakwater for incident wave angle

68 29 39 -65 15 -45 HS -65 .35

Figure 13. Amplification factor diagram for the thin semi-infinite breakwater for incident wave angle = 135 deg

24

Ne = = = ®SF @ SF @&

Amplification factor diagram for the thin semi-infinite breakwater for incident wave angle

-68 29 33 85 15 45 75 05 35

RN" - - & & * &

Amplification factor diagram for the thin semi-infinite breakwater for incident wave angle

- 68 29 33 -85 -1S 45 5 -85 -35

0.00 at the back wall of the wedge in the shadow zone. The diagram patterns are relatively smooth and simple.

27. Notably, the diagrams do not include phase information of the wave response which is usually unimportant in most engineering practice. Should the phase of the wave response need to be known, one can always use the com— puter program WEDGE to calculate it.

28. The contour diagrams of the amplification factor, similar to the ones presented by Wiegel (1962) and shown in the SPM (1984), were also plotted as typically shown in Figure 16. Examination of those contour diagrams indi- cates that, in subregion III, the results are identical to Wiegel's results. But in subregion II the contour patterns for the amplification factor (or the diffraction coefficient K" used in the SPM (1984)) equal to 1.0; thus the present results are far more complicated than Wiegel's. The author believes that Wiegel's results may lose accuracy because of insufficient resolution of the computational tools during the late fifties and early sixties. Neverthe- less, such inaccuracies are usually either tolerable or immaterial in most engineering practice.

29. The contour patterns in subregion I are very complex, and it is difficult to track specific contours. Therefore, for clarity only the patched

diagrams are presented, and the contour diagrams are omitted in this report. Vertical Wedge of 90-Deg Wedge Angle

30. When the wedge angle is equal to 7/2(6, = 31/2) , the vertical wedge occupies the entire fourth quadrant of the space as shown in Figure 17. Wave response was calculated at 1,100 grid points intersected at r/d = 0.5,(0.5),10.0 and 6 = 0,(1/36),31/2 for the incident wave angle a = 0,(1/12),7 . The wave response at the origin is obtained by substituting

v = 1.5 into Equation 27 as follows:

6(0,0) =< (32)

Those calculated results were used to construct the amplification factor dia- grams by following the same procedures described for the case of the thin semi-infinite breakwater. The diagrams are shown in Figures 18 through 27. Because of symmetry of the results, the diagrams for the incident wave angle

a > 31/4 can be obtained from those for an incident wave angle of tT-a.

2H,

Ne ee er i eee en ee eens —_

Figure 16. Amplification factor contour diagram for the thin

semi-infinite breakwater for incident wave angle = 90 deg Therefore, the diagrams for the incident wave angles omitted in this report.

Sil,

a = 1/6,(1/12),7 are

The diagrams indicate that, for each corresponding incident wave

angle, the results in subregion I are very similar to those obtained for the

vertical wedge of O0-deg wedge angle. However, the results in subregions II

and III from both wedges are discernibly different.

28

, - . , .— ? et t ra y ' HS, t Cr usin: At ies if Fo. Had : Bay “3 Ds t I ? Cs Dee Uae I t Os 5 1 ! ' I Bog Re | iL aeeneny beeen een spen abs === PRCT pBOcH DOSS poeact OSS UCG| > 4 1 | H 4 \ wee ey ee | ; | \ 5 Ys spethoOO Got See! Hts U 1 1 1 A ioe BO at? ol St tar 1 1 5 4 Zoaeo fy yao me 08 ie 1 \- a ry Ch OG, | \ Nanos \ eens OTS Geer 1 5 act 4 we Yow Pt eS a 4 = ' 4 Apts oh : et ies yee? + uA ine Oe Tae eet) 4 4 a=) : eee , iF ay Me ave ot -4 = “TO ee Se at ae Pe a ' AS ay) eee My OR Secs eae 4 Noonan \. Ds iM . ral 4 » ef Ya ' H 4 -% 1 Yo ~ ipaeeren oe vee Popa Qo a oe eee ee Sf * Ae tS we A i ats = H . _ oe ak Soe “ey, O t Desc “D. a ed x i ae a ' ty r % oe Ese a aS a 4 0 oes & aa ' ee ! ees Wae> Ss Be ee aan ~~ a ~s i y i anal bs, Baad Cea ' pate Sy ra micscs 4 -<T N aN 1 ~—.) . : we ' (> : e ' ~ ' : -_ a cer) u le ae Tinian H wet ~ . =e Pe s Saas ' y - Senet aac7 :

-=—--~

_-<--

A 90-deg wedge and polar coordinates 29

Figure 17.

Figure 18. Amplification factor diagram for the

90-deg wedge for incident wave angle

30

0 deg

—_~ = = © ® ® & FS S'

-08 -18 -26 48 -68 -88 -98 -18 -36

Figure 19, Amplification factor diagram 90-deg wedge for incident wave angle =

31

for the 15 deg

Nh ss & SO BS

-88 -25 “EES -85 1 -45 75 -85 -35

Nye = —- SF 2 @& ©

Amplification factor diagram for the 90-deg wedge for incident wave angle

-68 .25 -55 -85 15 -45 75 -85 35

Figure 21. Amplification factor diagram for the 90-deg wedge for incident wave angle = 45 deg

33

6.68 6.25 6.55 8.65 1.15 1.45 1.75 2.65 2.35

Figure 22. Amplification factor diagram for the

90-deg wedge for incident wave angle

34

60 deg

mm = = =—§ ©& & & @

.68 .25 55 65 -15 4S 7S -6S 35

NN = = = © FSF FS &

Amplification factor diagram for the 90-deg wedge for incident wave angle

-68 25 -55 .65 13 45 #5 65 .35

Figure 24. Amplification factor diagram for the 90-deg wedge for incident wave angle = 90 deg

36

8.68 8.25 6.55 6.85 1.15 1.45 1.75 2.05 2.35

Mme = = = & & S&S &

Amplification factor diagram for the 90-deg wedge for incident wave angle

-48 .25 -55 .85 13 -45 75 85 .35

Ne NM =e |= =| OO F&F &

Amplification factor diagram for the 90-deg wedge for incident wave angle

-68 25 -55 65 US 45 45 65 35

Nhe =|=- |= =| © FB F&F &

Amplification factor diagram for the 90-deg wedge for incident wave angle

-68 .25 -55 -65 15 4S 75 85 .35

PART IV: CONCLUSION

32. An analytical solution for the combined wave reflection and dif- fraction by a vertical wedge of arbitrary wedge angle is obtained and ex- pressed in Equation 26. The analytical solution is in terms of Bessel functions and in nonclosed form. The computer subroutine WEDGE, written for calculating the solution, is documented in Appendix A.

33. The amplification factor diagrams for a vertical wedge of O-deg wedge angle and a vertical wedge of 90-deg wedge angle are calculated and pre- sented. The calculated results indicate that the wave response in subre- gion I, where the incident, reflected, and scattered waves all exist, is in a very complex pattern with the amplification factor varying from 2.35 to 0.0 over the subregion. The wave response in subregions II and III is a rela- tively simple pattern with the amplification factor decreasing from 1.15 roughly along the reflected wave ray reflected from the origin point to nearly 0.00 at the back wall of the wedge in the shadow zone.

34. Diagrams of the special case of a vertical wedge of O-deg wedge angle can be considered complementary and extended versions to the ones pre-

sented in the SPM (1984).

40

REFERENCES

Abramowitz, M., and Stegun, I. A. 1964 (Jun). Handbook of Mathematical Func- tions, National Bureau of Standards, Applied Mathematics Series 55, pp 358-495.

Morris, A. H., Jr. 1984 (Jun). "NSWC Library of Mathematics Subroutines," NSWC TR 84-143, Strategic Systems Department, Naval Surface Weapons Center, Dahlgren, Va., pp 43-44.

Shore Protection Manual. 1984. 4th ed., 2 vols, US Army Engineer Waterways Experiment Station, Coastal Engineering Research Center, US Government Printing Office, Washington, DC.

Stoker, J. J. 1957. Water Waves, Interscience Publishers, Inc., New York, pp 109-133.

Wiegel, R. L. 1962. (Jan). "Diffraction of Waves by a Semi-infinite

Breaker,'"' Journal of the Hydraulics Division, American Society of Civil En- gineers, Vol 88, No. HY1, pp 27-44.

41

APPENDIX A: SUBROUTINE WEDGE

1. Subroutine WEDGE is used to calculate the value of 94 in Equa-

tion 26 which is generally a complex number. Its absolute value is the ampli- fication factor, and its phase is the phase indicator from the phase of the incident wave. As mentioned in the main text, 4 is a function of the Bessel function of either fractional or integer order, depending on the wedge angle, and is the aummaeien of a series of infinite terms. The subroutine BESJ, doc- umented in the Naval Surface Weapons Center (NSWC) Library of Mathematics Sub- routines (Morris 1984)* is used in the WEDGE subroutine. The programming of the WEDGE subroutine is very straightforward if a truncation term in the se- ries in Equation 26 is determined. The program is written in FORTRAN language

and is listed in this appendix.

Description

2. The following subroutine is available for computing 9 in

Equation 26: CALL WEDGE(F, FABS, FPHA,XRL, XTH, WEDGEA, WAVEA, IDX)

where the arguments are all real values except F which is a complex value.

Input arguments are as follows:

a. (XRL,XTH)=(r/A,60) where (r,9) are polar coordinates of the lo- cation where $ is to be computed, and A is the incident wave length. Therefore, XRL is the radius vector or radius distance normalized by the incident wave length. XTH is the vectorial angle in degree.

- WEDGEA = wedge angle in degree.

Io

c. WAVEA = incident wave angle in degree. d. IDX = an index (set to 0 in this subroutine).

Output arguments are as follows:

* References cited in the Appendix can be found in the References at the end of the main text.

Al

a. F= @ in Equation 26 (wave response normalized by incident wave amplitude).

b. FABS = amplification factor, the absolute value of 9.

ig ies] Las] Il

phase difference, the phase of $6.

Example and Test Run

3. To serve as an example as well as to ensure the subroutine is the correct one, the user should run the test program listed in Figure Al and make

sure the output is the same as that listed in Table Al.

FROGRAM WEDGE! (INPUT ,OUTPUT, TAPES=INPUT, TAFES=QUTPUT) TESTI CONPLEX F TEST2 PI=3. 141592654 TESTS WRITE (4,4) TEST4

4 FORMAT(//3%,' WEDG-ANG WAV-ANG LOCATION WAVE RESPONSE TESTS { ABS-VAL PHASE’/3X,’ (DEG) (DEB) XRL XTH(DEB) (NORMTESTS PALIZED) AMP-FAC (RAD) ‘/3X,35¢° -*)/) TEST?

WEDGEA=90, TEST WAVEA=135. TESTS XRL=2.0 TESTIO ¥TH=30. TEST11 IDY=0 TESTI2 CALL WEDGE (F,FABS,FPHA, XRL,XTH,WEDGEA,WAVEA, IDX) TESTI WRITE(6,40) WEDGEA,WAVEA,XRL,XTH,F,FARS,FPHA TESTIA 40 FORMAT(1X,8F9. 2) TESTIS STOP TESTI4 END TESTI7

Figure Al. Computer program list 1

Table Al

Sample Output of the Test Program

WEDG-ANG WAV-ANG LOCATION WAVE RESPONSE ABS-VAL FHASE (DEG) (DEG) RRL XTH(DEG) (NGRMALIZED) AMP-FAC iRAD) 90.00 135.00 2.00 40.66 a5 Ail 268 a fe Bo Mil

A2

4. The input arguments in the test programs are as follows: a. WEDGEA = 90; the wedge angle is 90 deg. b. WAVEA = 135; the incident wave angle is 135 deg.

c. XRL = 2.0 and XTH = 30; the location of the wave response $4 to be computed is at radial vector of two incident wave length dis-— tances and vectorial angle of 30 deg in polar coordinates.

d. IDX = 0; the index IDX is set to 0. This case is shown in Figure A2. The outputs are given in Table Al which is self-explanatory.

5. If the location of the wave response $ to be computed is a very large distance, for example XRL greater than 18 or so, the output might print the message that the number of terms is insufficient for summation in comput- ing $ . In this situation, the user must replace the integer 200 in PARAMETER(NN=200), listed in Card WEDGE15 in the SUBROUTINE WEDGE list, by a

larger integer to ensure the accuracy.

Figure A2. Example problem

A3

Subroutine WEDGE Listing

6. Subroutine WEDGE, which is listed in this section, calls subroutine BESJ which, in turn, calls subroutine JAIRY and function GAMLN. Subroutines BESJ and JAIRY and function GAMLN are borrowed from the NSWC Library of Mathe- matics Subroutines (Morris 1984). Including these borrowed subroutines here in the list is only for the purpose of allowing subroutine WEDGE to be self- contained and complete. The computer program of subroutine WEDGE is listed in

Figure A3.

A4

SUBROUTINE WEDGE(F,FABS,FPHA,XRL,XTH,WEDGEA,WAVEA, IDX)

THIS COMPUTER PROGRAM WAS WRITTEN BY H.S. CHEN OF CERC IN 1986 UNDER CORP’ CIVIL WORK R&D PROGRAM. NEITHER ANY OF AGENCIES NOR ANY INDIVIDUAL ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY FOR THE ACCURACY OF THE PROGRAM.

Es I 2

WAVE REFLECTION AND DIFFRACTION BY A VERTICAL WEDGE.

INPUT : XRL = R/L, (RADIUS VECTOR) / (WAVE LENGTH) XTH = VECTORIAL ANGLE IN DEGREE WEDGEA = WEDGE ANGLE IN DEGREE WAVEA = INCIDENT WAVE ANGLE IN DEGREE OUTPUT : F = FHI, WAVE RESPONSE NORMALIZED BY INCIDENT WAVE AMPLITUDE FABS = ABSOLUTE VALUE OF PHI, THE AMPLIFICATION FACTOR FPHA = PHASE OF PHI IN RADIAN.

PARAMETER (NN=200)

DIMENSION BJ (NN) ,W(NN) , XM(NN)

COMPLEX F,1M

DATA TOLR/1.E-8/,1TER/8/

PI=3.141592654

CPI=PI/180.

XKR=XRL#2. O#P I

TH=XTH#CPI

WA=WAVEA*CPI

XNU=(360.-WEDGEA) /180.

XM(1)=1.0

IF(IDX.NE.0) GOTO 14

CALL BESJ(XKR,0.0,1,W,NZ)

Bd(1)=W(1)

ICOUNT=0

DO 10 N=1,NN

NL=N+1

ICOUNT=ICOUNT+1

IF(ICOUNT.LE.ITER) GOTO 8

NNN=N

GOTO 14 8 XM(N1)=FLOAT(N) /XNU

M=INT(XM(NL) )

ALPHA=XM(NI) -M

ML=M+1

CALL BESJ(XKR,ALPHA,M1,W,NZ)

IF(N1.EQ.NN) WRITE(4,9) XKR,ALPHA,ML,W(ML) ,NZ 9 FORMAT(/'’ ##** NO. OF TERMS FOR SUMMATION IS INSUFFICIENT *#*%',

1 /’ XKR,ALPHA,M1,W(ML),NZ =',2F10.4,15,E15.6,15/)

BJ (N41) =W(ML)

Figure A3. Computer program list 2 (Sheet 1 of 25)

A5

WEGDEL

Re HE EF RE KR HE KE KK KR RR KR K€ EK KK HR HE KR KE KE RE & KF KF & E¥NOTICEL]

*¥NOTICE2 *NOTICES *NOTICE4 *¥NOTICES

# FF HK He EH EK HK KF RH HE KH HF HH HH HK KR HR HK HK HH HE HE KR KF RNOTICES

-WEDGE2 WEDGES WEDGE4 WEDGES WEDGES WEDGE7 WEDGES WEDGE? WEDGE10 WEDGE1L1 WEDGE12 WEDGELS

-WEDGE14 WEDGE1S WEDGE14 WEDGE17 WEDGE18 WEDGE19 WEDGE20 WEDGE21 WEDGE22 WEDGE23 WEDGE24 WEDGE25 WEDGEZ4 WEDGE27 WEDGE28 WEDGE29 WEDGES WEDGES1 WEDGES2 WEDGE3S3 WEDGES4 WEDGESS WEDGE364 WEDGES7 WEDGESS8 WEDGE39 WEDGE40 WEDGE41 WEDGE42 WEDGE43 WEDGE44

oa

AIAIAIYIAIGIDYIYVGYIYVWIMVWYW IIIA AIAAIIIBWIHIOIGN IHS

IF (ABS(BJ(N1)).GT.TOLR) ICOUNT=0 10 CONTINUE 14 CONTINUE F=BJ(1)/2. DO 20 N=1,NNN N1=N+1 XMN=XM (NL) TM= (0.0, 1.0) ##XMN#BJ (NI) #COS (XMN#WA) #COS (XMN*TH) F=F+TM 20 CONTINUE F=4./XNU*F FR=REAL (F) FI=AIMAG(F) FABS=SORT (FR¥FRtFI¥#FI) IF (WA.LE.1.E-8) FABS=FABS/2. IF(FABS.LT.TOLR) GOTO 30 FPHA=ATAN2(FI,FR) RETURN 30 FPHA=0.0 RETURN END SUBROUTINE BESJ(X,ALPHA,N,Y,NZ) WRITTEN BY D.E. AMOS, S.L. DANIEL AND M.K. WESTON, JANUARY, 1975. REFERENCE SAND-75-0147

ABSTRACT BESJ COMPUTES AN N MEMBER SEQUENCE OF J BESSEL FUNCTIONS J/SUB(ALPHAtK-1)/(X), K=1,...,N FOR NON-NEGATIVE ALPHA AND Xx. A COMBINATION OF THE POWER SERIES, THE ASYMPTOTIC EXPANSION FOR X TO INFINITY AND THE UNIFORM ASYMPTOTIC EXPANSION FOR NU TG INFINITY ARE APPLIED OVER SUBDIVISIONS OF THE (NU,X) PLANE. FOR VALUES OF (NU,X) NOT COVERED BY ONE GF THESE FORMULAE, THE GRDER IS INCREMENTED OR DECREMENTED BY INTEGER VALUES INTO A REGION WHERE ONE OF THE FORMULAE AFPLY. BACKWARD RECURSION IS APPLIED TO REDUCE ORDERS BY INTEGER VALUES EXCEPT WHERE THE ENTIRE SEQUENCE LIES IN THE OSCILLATORY REGION. IN THIS CASE FORWARD RECURSION IS STABLE AND VALUES FROM THE ASYMPTOTIC EXPANSION FOR X TO INFINITY START THE RECURSION WHEN IT IS EFFICIENT TO DO SO. LEADING TERMS OF THE SERIES AND UNIFORM EXPANSION ARE TESTED FOR UNDERFLOW. IF A SEQUENCE IS REQUESTED AND THE LAST MEMBER WOULD UNDERFLOW, THE RESULT IS SET TO ZERO AND THE NEXT LOWER ORDER TRIED, ETC., UNTIL A

MEMBER COMES ON SCALE OR ALL MEMRERS ARE SET TO ZERO. OVERFLOW CANNOT OCCUR. BESJ1 CALLS SUBROUTINE JAIRY AND FUNCTION GAHMLN. DESCRIPTION OF ARGUMENTS INPUT X | Hollen @

ALPHA - ORDER OF FIRST MEMBER OF THE SEQUENCE, ALPHA.GE.90

N - NUMBER OF MEMBERS IN THE SEQUENCE, N.GE.1 OUTPUT Y - A VECTOR WHOSE FIRST N COMPONENTS CONTAIN

VALUES FOR J/SUB(ALPHAtK-1)/(X), K=1,...,N

Figure A3. (Sheet 2 of 25)

A6

WEDGE45 WEDGE46 WEDGE47 WEDGE48 WEDGE49 WEDGES WEDGES1L WEDGES2 WEDGESS WEDGES4 WEDGESS WEDGES6 WEDGES7 WEDGES8 WEDGES9 WEDGE6O WEDGES1 WEDGEG2 WEDGES3S WEDGE64 WEDGESS BESJ101 BESJ102 BESJ103 BESJ104 BESJ105 BESJ106 BESJ107 BESJ108 BESJ109 BESJ110 BESJi11 BESJ112 BESJ113 BESJ114 BESJ1I15 BESJ116 BESJ117 BESJ118 BESJ119 BESJ120 BESJ1i21 BESJ1i22 BESJ1253 BESJ124 BESJ125 BESJ126 BESJ127 BESJ128 BESJ129 BESJ130 BESJi31 BESJ132 BESJ133 BESJ134

OHOonrw i ~Iianrnraianain

ise)

ia]

On

1

1

5

1

NZ - ERROR NZ=0 NZ=-1 NZ=-2 NZ=-3 NZ.GT.

ERROR CONDITIONS

IMPROPER INPUT AR UNDERFLOW

DOUBLE PRECISION DX,T DIMENSION Y(N)

DIMENSION Ciii,10) DIMENSION C1(88) DIMENSION A1(952), DIMENSION GAMA(26) ,TE DIMENSION FNULIM(2) ,P DIMENSION CR(10) ,DR(1

EQUIVALENCE EQUIVALENCE EQUIVALENCE EQUIVALENCE EQUIVALENCE EQUIVALENCE EQUIVALENCE

ROG) sc (11,9) ,€ (ALFA(L,1 (BETA(1,1 (BETA(1,5 DATA ELIM1,ELIM2, TOL

DATA PP(1)/8.72909153 PP(3)/1.24578576

TOLS=LN(1.E-3) DATA TOLS

DATA UPOL(1) ,CON1,CON 6, 6666666666667E- ol, 1.0416656566667E-01/

DATA RTWO,PDF,RTTP,PI 7. 8539816339745E-01,

DATA FNULIMN(1)/100./,

- A NON-FATAL ERROR

vALFA(26,4) “(C2 (22 A2(52)

(ALFA(1,3),

(BETA(1,3)

INDICATOR NORMAL RETURN - COMPUTATION COMPLETED X 15 LESS THAN 0.0 ALPHA IS LESS THAN 0.0 N 1S LESS THAN 1

0 LAST NZ COMPONENTS OF Y SET T0 0.9 BECAUSE OF UNDERFLOW

GUMENTS - A FATAL ERROR (NZ.6T.0)

RX ,DTM, DEN BETA (26,5)

BL (52) ,B2(52) ,B3(26)

MP (3) ,KMAX(5) ,AR(8) ,BR(10) P(4)

0)

wUPOL (10)

1(1)) 2(1)) ),A1(1)) PAN) BHA) »B21)) ),B3(1))

i Obds ; 44, f 1.E-15 PP(2)/2.6569393226503E-01/,

Q35SSE+00/, E PP(4)/7.70133574743039E-04/

3555 B6559E-O1/,

/-6.9077552799821E+00/

2,CONS ,CONS48 / L OOTOOOODODODOE SO: SB, 355333333333 3E-01, 1.4142155623731E+00, DT / 1.3483997249265E+00,

7.9788456080286E-01, 1,.5707963267949E+00/

FNULIM(2)/60./

CE=-ALOG(TOL) , TCE=-0.75*ALOG(TOL) DAA CE nee / J.45387 763749 LE+O1, 2.5904082296185E+01/ DATA INLIN / 150 /

DATA AR(1)/8.35503472

Figure A3.

AR(2)/1

282265745563 5E-

(Sheet 3 of 25)

22222E-02/,

A7

~.

BESJ135 BESJ136 BESJ137 BESJ138 BESJ1359 BESJ140 BESJ141 BESJ142 BESJ143 BESJ144 BESJ145 BESJ146 BESJ147 BESJ148 BESJ149 BESJ150 BESJ151 BESJ1S2 BESJ153 BESJ154 BESJ1595 BESJ1564 BESJ157 BESJ158 BESJ159 BESJ160 BESJ161 BESJ162 BESJ165 BESJ164 BESJ1465 BESJ146

(BESJ167

BESJ168 BESJ169 BESJ170 BESJ1i71 BESJ172 BESJi73 BESJ174 BESJ1735 BESJ176 BESJ177 BESJ178 BESJ179 BESJ180 BESJ181 BESJ192 BESJ183 BESJ184 BESJ185 BESJ186 BESJ187 RBESJ188 BESJ199A

ted Poe

Gl hoe

oOo oon bein OOO OF Sted be

con Pid he OOO OS ed he

oN Oo} Om et hoe

ARCS) /2.9194902646414E-91/, AR(4)/8.8162725744375E-01/,

AR(3)/3. AR (7) /

BRL) BR(3) BR(S) BR(7) BR(9)

Citi) Ci(3) G7)

7 /

/

OME) 7

Ci(st)/

Ci(s6)/ C1(38)/ Cic4i)/

C1(S1)/

C1(55)/ 0 CCRT) P= 7/0 C1(59)/-6. EWG) (2s

C1(63)/

£1(78)/

C1(82)/

C2¢1) C2(3) C2(35) C2(7) C2(9)

SS Ss SS

/ /

GAGta) 7

»4291918790055E+05/, C2(2) ~P98O1LSILBSSB1E+06/, C2(4) »81356352265865E+06/, C214) »3164517248456E+05/, C2(8) »4998S0481B8112E+03/, C2(10)/ ©O/, €2(12)/ 3.2844698520720E+06/, C20) POI a C2(15)/-7. EACH) [3c

214082818628E+00/, AR(6)/1.49957629869535E+01/, >B925013011S86E+01/, AR(8)/4.74451539859246E+02/

S8SSSSSSSS5SE-G1/, BRI2) /-9.8741319444444E-02/, a QOSS91IS9IOE-G1i/, BR(4) /-3.1722720267841E-01/, 242914795712E-01/, BRIS) /-3.5112030408264E+00/, 72726356203568E+01/, BR(8) /-8.2281439097186E+01/, 23555 7052367E+02/, BRIL0)/-3.35162185685480E+03/

»OB8S3S5S533335SE-01/, C1(2) / 1.2500000000000E-017, sO gy GCE) 7 WoO 5 GUMS) 7 Mo07, Gite) 7% WoO, DoOUg UNE) 7 WoO/5 GUM) # O75, ECXCLOIY Oo0/, EXEC) F/O Cia) Jaa alls) 7 Wo

O/, Ci(12)/ 3,.3420138888889E-01/7, 0104166666667E-01/, C1(14)/ 7.0312500000000E-02/, Diy (GOte i DoW, (OWN WoO, TOC Wali.

50/4 ENCAO)/ WsO7, (ENCE) Oo0/5 NCAA Wst0/ 4 (i (23) Faia Ey CA) So CCA) AO

O2581259764506E+00/, C1(24)/ 1.84646267356111E+00/, FL2Z1LO93S7I0OOOE-O1/, Cit26)/ 7. 3242187500000E-02/,

SOG ENCE Ose, LNCS )/ WoOi, Ck tsOo 7 Oot0/ 5 50/5 MKSAN/ O07 5 CLUS) DoO/. Ci(34)/ 4,

6695844254262E+00/, CL(35)/-1.1207002616223E+01/, 7TB9L2Z3S5S51562E+00/, C1137) /-2.35640869140625E+00/,

~1215209960938E-01/, C1(39)/ 0.0/7, C1(40)/ 0.0/7, -O/, C1(42)/ 6.0/, C1(43)/ 0.0/7, C1(44)/ 0.0/,

Ease) (Bq C1(47)/-9. (ECA) P= 7c

B212072558200E+01/, C1(46)/ 8.4636217674601E+017, 181824154S5240E+01/, CL(48)/ 4.25349987453588Et01/, 3687943594796E+00/, C1(5G)/ 2.2710800170899E-01/7,

»O/, C1(S52)/ 0.0/, CL1(53)/ 0.07, €1(54)/ 0.0/

-0/, C1(56)/ 2,.1257015003922E+02/,

6525246814118E+02/, Ci(58)/ 1.0599904525280E+03/, 993579627357 613E+02/, CL(60)/ 2.1819051174421E+02/, G4914350486952E+01/, C1(62)/ 5.7250142097473E-017,

»0/, €1(64)/ G.0/7, C1(65)/ 0.0/7, €1(66)/ 5.0/7, C1(67)/-1. (ah C\5))) 7 8g (3 (7/1) Yee (EN 7/35) PON 6 CHES) We

9194576623194E+03/, C1(68)/ 8.0617221817373E+03/, 3586550006434E+04/, C1(70)/ 1.165539353546864E+04/, 30564697B6134E+03/, C1(72)/ 1.2009029132164E+03/, 0809091978B840E+02/, C1(74)/ 1.7277275025845E+00/, Oly ELETO)E O.005 CLLPTIE D.O%,

»O204291330966E+04/, C1(79)/-9. 6980598 S88638E+04/, Ci(8o0)/ 1.

254700123253E+05/, C1(81)/-2.0340017728042E+05/,

»2220046498S502E+05/, C1 (835) /-4.1192654968898E+04/, C1(84)/ 7. C1(86)/ 4.

1095143024894E+03/, C€1(85)/-4.9391530477309E+02/, 0740420012735E+00/, C1(87)/ 0.0/, C1(88)/ 0.0/

1. S11763561466350E+06/, 3.7632712976564E+06/ , 1,.2683652733216E+04/, 4, 2

SS

o218768981263E+047, 43805296995 56E+01/,

970GB1911B4S2E+07/, C214) / 5.09526024926465E+07/, ALOS1L4B21153535E+07/, C2(16)/ 6.6344512274729E+07/, 7567176660763E+07/, C2(18)/ 1.35288767166422E+07/,

Figure A3. (Sheet 4 of 25)

A8

BESJ189B BESJ189C€ BESJ189D BESJ190

ESEJ1914A BSEJ191B BSEJI91C BSEJ191D BSEJI91E BESJ192

BESJ193A BESJ193B BESJ193C BESJ193D BESJ193E BESJ195F BESJ1936 BESJ193H BESJ1931 BESJ193J BESJ193K BESJ193L BESJ193M BESJ193N BESJ1930 BESJ1935P BESJ9310 BESJ1935R BESJ193S BESJ194

BESJ195A BESJ195B BESJ195C BESJ195D BESJ195E BESJ199F BESJ1956 BESJ195H BESJ1951 BESJ195d BESJ 195K BESJ19SL BESJ1935M BESJ195N BESJ1950 BESJ194

BESJ197A BESJ197B BESJ197C BESJ197D BESJL97E BESJ197F BESJ1976 BESJ197H BESJ1971

DATA

oOoOnyA oon Pol ro K— OOO OW SP he

DATA

oon eS Ww Pe

DATA

Oommen Str OO CON OO eed Pe

DATA

cd Poe

AOU) a Al (3) / ALCS) / AL(7) / ALC?) # AL(11)/

AL(L7)/ 1

A1l(37

AlL(39)/-7. A1(41)/-5. AL(43)/-4., AL1(435)/-3. AL(47)/-2.

A1(49)/-2

AL(S1)/-1. AZ(L) /- A2(3) / A2(5) / A2(7) / A2(9) 7 A2(11)/ 7

A2(19)/ 3

A2(2 A2(25

A2(3 A2(37

AZ(39)/-1. A2(41)/ 2. A2(43)/ 4. A2(45)/ 5.

» 7856181 280864E+06/, »SBB6089753717E+04/,

-4444444444444E-05/, .8489288489288E-05/, 466913572741 79E-04/, .6182429706150E-

04/7, 32721040083235E-64/,

.0073885876275E- 047, AL(13)/ 1. ANUS) We

73060773591 788E- s027477415091E-

o4/, 04/,

»5166881654592E-04/, AL(19)/ i. AL(21)/ ft. AL(23)/ 9.3 AL(25)/ 8. AL(27)/ 6. AL(29)/-1. AL(S1)/-1. AL(33)/-1. ALCS) 7-1. VID

1640527147474E-047, 0377241042299E-64/, 212051724950E-035/, 4296310571570E-05/, 93735541355 459E-047, 4198627255669E-05/, SOBNSSSGOS305E-047, 4680975564647E-04/, 1974497568425E-047, S7699S54989119E-05/7,

2937434815522E-05/7, 6961156600937E-035/, 4813486800988E-05/, 5540053297204E-05/, BS99679S90418E-05/,

»2845929716472E-05/,

8481621752767E-05/,

S4A21197145774E-04/,

O446550359494E-05/, »6747110669971E-04/, »J650142760860E-04/, »1488669202982E-04/, » 644984192 A2Z(13)/ 4, A2(15)/ 2. AZ(17)/ 1.

5090E-05/, 7439429929051E-05/, 6993971497922E-05/, 3005667479396E-05/,

~J925748581935E-06/, AZ(21)/-2, a) / 2c Wier ke A2(27)/ 3. A2(29)/-6. A2(31)/-3. A2(33)/-2. cD) Hie iee )/-6.2

6539676969794E-06/, 7273929609125E-06/, S130471509356E-067, 7819419920177E-04/7, 3793850631886E-05/7, 1091625602736E-04/, 78950273791 352E-04/, 7TSE2Z4528034085E-047, 298287282021E-035/7,

72335430236696E-05/, S139544314829E-05/, 53800612449019E-05/, 6839620854582E-05/7,

Figure A3.

C2 / C /

0) 2) AL(2)

Al(6) / AL(8) /

AL(12)/ 1

AL(20)7

AL(40)/-6.

A1(42)/

AL(48) /-

AL(52)/-

(N30) YS

A2(6) / A2Z(8) /

A2Z(10)/ 9. A2(12)/ 4. A2(14)/ 3S. A2(16)/ 1. A2(18)/ 7. A2(20)/ 1. aye A2(24)/-8.1 A2(25)/-1. A2(28)/ 2. A2(30)/-2. A2(32)/-3. A2(34)/-2. A2(36)/-1.25 A2(38)/-4.

A2(22

A2(40)/ SS. A2(42)/ 3. 5. 245952949591 1E-

A2(44)/

A2(46)/ 3.

920 222)

2.6599558954626E- 2.4873 AL(10)/ 2.1 ~8626763663754E-04/, AL(14)/ 1. A1(16)/ 1. EEO Ihe 1.0979829837271E-04/7, CAN A1L(24)/ 8. Al(26)/ 8. AL(28)/ 2. AL(30)/-1. OU CS2) Joos AL(34)/-1. AL(36)/-1. A1(38)/-8.

9 (YA) ho

-0818640461266E+05/, »1001714062

6925E+02/

07792 6592768783

207792E-04/, 245E-04/7, 04/, O43734466E-04/, 6342485712346E-04/,

4109170592902E-04/, 4O50349739127E-047, 2366744559825E-04/,

82526078369356E-035/, 8571085247871E-05/, 0249754B40779E-05/, 2224174518292E-04/, 1644493147205E-04/, 12192491B810E-04/, 3381550396749E-04/, 0618431920797E-04/, 2692304558819E-05/

4404

255772102E-05/,

-5.0473104430356E-05/, A1(44)/-3. Ai(46)/-3. 2.0495 A1(50)/-2. 1.6651933002139E-05/

9868872771760E-05/, 1741425660902E-05/, 2272063487E-05/, 0535275310648E-05/,

eo ice eee sees Uae (HME Ma 1 1

Seen ee -04/, 0222587683 a 04/,

Gisroionae vas,

4586909303469E-05/, O0757033496520E-05/, 6275751200554E-05/, 9321093824794E-05/, 8262086574450E-06/, 4404004981425E-07/7, 9134686709849E-06/, 1726937967866E-06/, 0201141879802E- 05/, 0247195276182E-04/, SBS59B230603501E-04/, 1368011524758E-04/, 2856408261914E-04/, S4406306069E-04/7, 6286073058812E-05/

606904B239460E 62642745985679E

-06/, -05/, 05/, 9A349B2039310E-05/,

(Sheet 5 of 25)

A9

BESJ197J BESJ197K BESJ198

BESJ199A RESJ199B BESJ199C BESJ199D BESJ199E BESJ199F BESJ1996 BESJ199H BESJ1991 BESJ199d BESJ199K. BESJ199L BESJ199M BESJ199N BESJ1990 BESJ199P BESJ199@ BESJ199R BESJ199S BESJ200

BESJ201A BESJ201B BESJ201C BESJ201D BESJ201E BESJ201F BESJ2016 BESJ202

BESJ203A BESJ203B BESJ203C BESJ203D BESJ203E BESJ2035F BESJ 2036 BESJ203H BESJ2051 BESJ203) BESJ 203K BESJ203L BESJ 203M BESJ205N BESJ2030 BESJ203P BESJ2030 BESJ203R BESJ2035 BESJ204

BESJ205A BESJ205B BESJ205C BESJ205D

con eS

Om NO Of eS ed ro OO NO OP ol he

ce on fs Gy ho oe

te} PO ~O CO NS & Of & Oy PO

oonoowS

boo

DATA

DATA

DATA

AZ(47)/ &. A2Z(49)/ &. A2(S1)/ 3.7

Bl (1) Bi(3) BL (3) Bi(7) Bi (9) Bi(ii)/ Bi (13)/

SSS Ss SS PI td fe SI or bo

Bi(2

si (27) #1 B1L(29)/

0647852757842E-035/, O157789453946E-035/, FABLES HI VEVESORY 5

Ger

»7998872141355E-02/, ~BBS014022351135E-03/, 22475 »1443042172729E-04/, »694310076064BE- ~S481888931850E-04/7, sd2211615 Jeol (OS) fhe BNA YY the BL CHS) / the BUEG2119) ele SY We El Cza) 7 sn

311058920E-03/, 04/,

S4957E-04/, 9754183803305E-04/, 5931689966182E-04/, 3144806811996E-04/, 1044914450460E-04/, 4199822420424E-05/, 1346626214280E-05/,

4928295321 543E-037, -5.0291654957204E- Ba Saale

04/, 7946399697078E-047,

B1(33)/-3.96141953504646E-035/, B1(35)/-1.2608973578025E-057, B1(37)/ 8.0999616541427E-06/, B1L(39)/ 1.7396412547293E-035/, BiC4i)/ 2.1446526379082E-05/, B1(43)/ 2.2896778381471E-05/, B1(45)/ 2.5032197608091E-035/, BL{47)/ 2.2500588110529E-05/, B1(49)/ 2,.1641842744810E-05/, BL(S1)/ 2.0638874978217E-035/, B2(1) / 3.9221507672129E-04/, B2(3) / 2.7952063399Z202E-04/7, 245) / 6.9327110565704E-035/, B2(7) /-1.55744996354327E-035/, B2(9) /-4.18861861459669E-035/, B2(11)/-4.8766544741379E-03/, B2(13)/-4.7473562089009E-035/ B2(15)/-4.3353509649451127E- 05/, BGS 5,9482263860322E-05/, eee 3.3779330612337E-05/, 2(21)/-2.9526956173081E-057, ie 3) /-2.5800617466688E-035/, B2(25)/-2. 2582355095 1835E-05/, B2(27)/-4.74617796559976E-04/, B2(29)/-3.2039022806704E-04/, B2(21)/-4.2577810128544E-057, B2(33)/ 7.9709268407968E-05/7, B2(35)/ 1.1246677525220E-047 , B2(37)/ 1.0865163484877E-04/, B2(39)/ 9.2929839659355E-05/, WE Pea eA SNS Wai ¢ B2(43)/ 35.92596454732320E-05/,

Figure A3.

AZ (48) / A2(50)/ A2(352)/

BL(2) / B1(4) 7 Bi(é) / Bi(8) / Bicto)/ Bi(t2)/ B1(14)/ Bi(16)/ B1(18)/ B1(20)7/ Bi(22)/ B1(24)/ Bi¢26)/ B1(28)/ B1(30)/ I C2) B1i(34)/ B1(36)/ B1(38)/

B1(40)/ Bi(42)/ B1i(44)/ B1 (46) / B1 (48) / Bi(50)/ 61(52)/

BAVA) ff B2(4) / 2Co) F H2(8) 7 2(10)/ 2Ch2) 7 B2(14)/ B2(14)/ B2(18)/ Be (20) 7 B2(22)/ B2 (24) 7 B2(26) / AL CAE})) 7 B2(30) / A CRA) / B2(34)/ H2(36)/ B2¢38) /

B2(40)/ B2(42)/ 2(44)/

A10

FI td ee FD OO

cn con c-

SOF eee pre ono on

PO rl br ho bro Po or

oo = a 2

Rr id ore BPP Po td te

r=

cn c- co

» 08023 »89199657354470E-0357,

30785 »2823607372035E-05/7, »20981015356199E- © 1150754925622 -0116924199708E

90778844E-035/,

5280437558585E-05/

-T99S491L1L064359E-03/, -9009660676105E-

O3/,

2287887657294E-04/,

.717872817897GE-047, »93235283546292E-04/, 8895 »2228058079888E- -76836855901972E-04/, »A4S479OSOLITISE- » 20245444949350E-04/, »0182877074057E ~7413054575385E »9900226964622E-05/, «782047095 46S9E-04/, 9482 -9400855046082E-04/, »12035892907610E-05/, ©42892608575735E-07/, »53690700926215E-

214849575E-04/7,

o4a/, O4/,

-04/ 5 -05/,

2135851275E-04/,

03/

»9867297884215E-05/,

23994659923246E-05/, S8981118E-05/,

O5/, Wai! 4 -05/

793258155238E-04/, ee cert EN 04/7, 25868306999E-05/,

TAO ISES WS 05/7, 4, -8701003118674E-05/, »cSOBLSOS8LS8463E »0923501931577

6900488937914E-05/, -05/, 5E-05/, S0B5716753541E-05/,

~1588856077211E-035/, ~7oF7 891 482834E-05/, - 4130983 »1147965676891E »7786456714732E-04/, »6110501611996E »F457129429497E-05/7, »OS138236 »1310364210848E -O143795159766E-04/

5676128E-05/, =05/,

-04/,

OB27E-04/, -04/,

AO2Z9S1LSS0L609E-035/, SVE SAAZ NG /a Se o5/, 2216930882699E-05/,

(Sheet 6 of 25)

BESJ205E BESJ205F BESJ2056 BESJ206

BESJ207A BESJ2078 BESJ2070 BESJ207D BESJ207E BESJ207F BESJ 2076 BESJ207H BESJ2071 BESJ207d BESJ207K BESJ207L BESJ207M BESJ207N BESJ2070 BESJ207P BESJ 2070 BESJ207R BESJ 2075 BESJ208

BESJ 209A BESJ209B BESJ209C BESJ209D BESJ209E BESJ209F BESJ 2096 BESJ210

BESJ211A BESJ211B BESJ2Z11C BESJ211D BESJ211E BESJ211F BESJ2116 BESJ211H BESJ2111 BESJ211J BESJ 211k BESJ211L BESU211M BESJ211N BESJ2110 BESJ211P BESJ2110 BESJ211R BESJ2115 BESJ2t2

BESUZ13A BESJ213B BESJ213C

moon

M4 poe ~O CON O&O OWS OP re ao cn & C4

4 Poe

Beir O ON oO On SS

B2(45)/ 4. 5853948516536E-05/, 82(46)/ 4.0144551389149E-05/, B2(47)/ 3, 5048173003133E-05/, B2(48)/ 3.0515799503435E-05/, B2(49)/ 2.64935611995052E-05/, B2(50)/ 2.2936363369100E-05/, B2(51)/ 1.9789305666402E-05/, B2(52)/ 1.7009198463541E-05/ DATA BS(1) / 7.3646581057258E-04/, B3(2) / 8.72790805144619E-04/, B3(3) / 6.22614862573514E-04/, B3(4) / 2.8599815419430E-04/, B3(S) / 3.8473757287937E-06/, B3(6) /-1.8790600363597E-04/, B3(7) /-2.9760364659456E-047, B3(8) /-3.4599812683267E-04/, BS3(9) /-3.53558247091604E-04/, B3(10)/-3.3571563577505E-04/, B3(11)/-3.04352112478904E-04/, B3(12)/-2,.6672272304761E-04/, BS(13)/-2.2755421412282E-04/, B3(14)/-1.8992261195456E-04/, B3(15)/-1.5505891859909E-04/, B3(16)/-1.2377824076187E-04/, B3(17)/-9.6292614771764E-05/, BS(18)/-7.2517832771442E-05/, B3(19)/-5. 2207002889563E-05/, #3(20)/-3.5034775051190E-G5/, B3(21)/-2.0848976103555E-05/, B3(22)/-8.7010609684977E-04/, B3(23)/ 1.1569868667510E-G6/, B3(24)/ 9.1642647412278E-06/, B3(25)/ 1.5647778542887E-05/, B3(26)/ 2.0822362948247E-05/ DATA GAMA(1L) /6.29960592494744E-01/, GAMA(2) /2.5198420997898E-01/, GAMA(S) /1.54790300415466E-01/, GAMA(4) /1.10713062414616E-01/, GAMA(S) /8.5730939552740E-02/, GAMA(6) /6.9716131695868E-02/, GAMA(7) /5.8608567199371E-02/, GAMA(8) /5.04698873535631E-02/, GAMAI9) /4.42600580689146E-02/, GAMA(10)/3.9372066154351E-02/, GANA(11)/3.5428319592446E-02/, GAMA(12)/3.2181885750210E-02/, GAMA(13)/2.9464624079116E-02/, GAMA (14) /2.7158167711293E-02/, GAMA(15)/2.5176827297386E-02/, GAMA (16) /2.3457075530608E-02/, GAMA(17)/2.1950839013491E-02/ GAMA (18) /2.0621082823565E-02/, GAMA(19)/1.94388240897988E-02/, GAMA(20)/1.838106S380068E-02/, GAMA(21)/1.7429321323196E- 02/, GAMA (22) /1.65685835778661E-02/, GAMA(23) /1.5796528598792E-02/, GAMNA(24) /1.5972950149410E- -02/, GAMA(25) /1.4419325083994E-02/,

GAMA(26)/1.

S8184805973534E-02/

TEST INPUT ARGUMENTS

NZ=

KT=1

IF(N-1) 92

Kale?

NN=N

IF (X)

IF (ALPHA)

POE

IF (N.EQ. 1)

Ie?

60 TO 118

Til

DOMit9) l=11).N (IN 0F

RETURN

CONTINUE

IF (ALPHA.LT.O.)

108,

109

93,110,120 91,114,116

RETURN

60 T0 91

Figure A3.

All

(Sheet 7 of 25)

BESJ213D BESJ213E BESJ213F BESJ2136 BESJ214 BESJ215A BESJ215B BESJ215C BESJ215D BESJ215E BESJ219F BESJ2156 BESJ215H BESJ2151 BESJ213d BESJ 215K BESJ215L BESJ215M BESJ2164 BESJ217A BESJ217B BESJ2170C BESJ217D BESJ217E BESJ217F BESJ2176 BESJ217H BESJ2171 BESJ217J BESJ 217K BESJ217L BESJ217M BESJ217N BESJ218 BESJ219 BESU220 BESJ221 BESJ222 BES ae BESJ2 Rie BESJ22 BESJ22 BESJ22 BESJ229 BESJ23 BESJ231 BESJ232 seer ee BESJ23 BESJ237 BESJ23 BESJ239

[os al oe a oe a oe

oO

(52) f3) (cal (se) fs)

850

840

880

DFN=DBLE (FLOAT (N) )+DBLE (ALPHA) -1.D+0 FNU=DFN

XO2=X*.5

SX02=X02*X02

DECISION TREE FOR REGION WHERE SERIES, ASYMPTOTIC EXPANSION FOR X TO INFINITY AND ASYMPTOTIC EXPANSION FOR NU TO INFINITY ARE

APPLIED.

IF(SXO02.LE. (FNUt+1.)) GO TO 850 TA=AMAX1 (20. ,FNU) IF(X.GT.TA) GO TO 880 IF(X.GT.12.) GO TO 840 XO2L=ALOG (X02) NS=SX02-FNU

60 TO 852

FN=FNU

FNPI=FN+1.

XO2L=AL0G (X02)

IS=KT

IF(X.LE.90.5) GO TO 134 NS=0 DFN=DFN+DBLE (FLOAT (NS) ) FN=DFN

FNPL=FN+1.

IS=KT

IF(N-1+NS.G6T.0) I5=3

60 TO 134 NS=AMAX1(36.-FNU,O.) DFN=DFN+DBLE (FLOAT (NS) ) FN=DFN

IS=KT

IF(N-1+NS.G6T.0) I5=3

60 TO 130

CONTINUE

RTX=SQRT(X) TAU=RTWO#RTX TA=TAUtFNULIM (KT) IF(FNU.LE.TA) GO TO 500 FN=FNU

IS=KT

UNIFORM ASYMPTOTIC EXPANSION FOR NU TO INFINITY

CONTINUE

XX=X/FN

W2=1.-XX#XX

ABW2=ABS (W2)

RA=SQRT (ABW2) IF(ABW2.67.9.2775) GO TO 200

CASES NEAR X=FN, ABS(1.-(X/FN)¥*#2).LE.0.2775 COEFFICIENTS OF ASYMPTOTIC EXPANSION BY SERIES

Figure A3. (Sheet 8 of 25)

Al2

PESJ240 BESJ241 BESJ242 BESJ243 BESJ244 BESJ245 BESJ246 BESJ247 BESJ248 BESJ249 BESJ290 BESJ251 BESJ252

BESJ258 BESJ209 BESJ260 BESJ261 BESJ262 BESJ263 BESJ264 BESJ265 BESJ256 BESJ267 BESJ268 BESJ269 BESJ270 BESJ271 BESJ272 BESJ2753 BESJ274 HESJ275 BESJ276 BESJ277 BESJ278 BESJ279 BESJ280 BESJ2981 BESJ282 BESJ283 BESJ284 BESJ285 BESJ286 BESJ287 BESJ288 BESJ289 BESJ270 BESJ291 BESJ292 BESJ293 BESJ294

Iona

ba bh

ioe) noonn

ra

ZETA AND TRUNCATION FOR A(ZETA) AND B(ZETA) SERIES BESJ295

BESJG296 KMAX IS TRUNCATION INDEX FOR A(ZETA) AND B(ZETA) SERIES=NAX(2,SA) BESJ297 BESJ298 SA=0. BESJ299 IF (ABW2.EQ@.G.) GO TO 21 BESJ300 SA=TOLS/ALOG (ABW2) BESJS301 SB=SA BESJ302 DO 22 1=1,5 BESJ303 KMAX (IT) =AMAX1(SA,2.) BESJ304 SA=SA+SB BESJS05 CONTINUE BESJ306 KB=KMAX (9) BESJ307 KLAST=KB-1 BESJ308 SA=GAMA (KB) BESJ309 DO 24 K=1,KLAST BESJ310 KB=KB-i BESJS11 SA=SA#W2+GAMNA (KB) BESJ312 CONTINUE BESJ313 Z=W2*SA BESJ314 AZ=ABS(Z) BESJS15 RTZ=SQRT (AZ) BESJ3516 FNIZ=FN¥*CON2 BESJ317 RTARY=RTZ#FN13 BESJ318 ARY=-RTARY*#RTARY BESJ319 AZS2=AZ#RTZ*CON1 BESJS2 ACZ=FN*¥AZ32 BESJ321 MEN oLE55)) (0) Wo) 227/ BESJ522 BESJ323 TEST FOR UNDERFLOW, 1.E-2B0=EXP(-444.), ONE WORD LENGTH BESJ324 UP FROM UNDERFLOW LIMIT OF CDC 46600 BESJ325 BESJ326 IF(ACZ.GT.ELIM2) GO TO 180 BESJ327 ARY=-ARY BESJS2 PHI=SQRT(SQRT (SA+SA+SA+SA) ) BESJ32 BESJ330 B(ZETA) FOR S=0 BESJ331 BESJ332 KB=KMAX (3) BESJ333 KLAST=KB-1 BESJ334 SB=BETA(KB, 1) BESJS39 DO 23 K=1,KLAST BESJ334 KB=KB-1 BESJ337 SB=SB#W2+KETA(KB, 1) BESJ358 CONTINUE BESJ339 KSP1=1 BESJ340 FN2=FN#FN BESJ341 RFEN2=1./FN2 BESJ342 RDEN=1. BESJ343 ASUM=1. BESJ 344 RELB=TOL*ABS (SB) BESJ345 BSUM=SB BESJ346 DO 25 KS=1,4 BESJ347 KSPI=KSP1+1 BESJ348 RDEN=RDEN#RFN2 BESJ349

Figure A3. (Sheet 9 of 25)

Al3

oa

ioe)

ann

onan

A(ZETA) AND B(ZETA) FOR S=1,2,35,4

KB=KMAX (5-KS) KLAST=KB-1 SA=ALFA(KB,KS) SB=BETA(KB,KSP1)

DO 26 K=1,KLAST KB=KB-1 SA=SA#W2+ALFA(KB,KS) SB=SB#W2+BETA(KB,KSP1) CONTINUE

TA=SA#RDEN

TB=SB*#RDEN ASUM=ASUM+TA BSUM=BSUM+TB IF(ABS(TA).LE.TOL.AND.ABS(TB).LE.RELB) GO TO 152 CONTINUE

CONTINUE

BSUM=BSUM/ (FN*#FN1S)

60 TO 400

CONTINUE TAU=1./RA T2=1./W2 IF(W2.GE.0.) 60 TO 30

CASES FOR (X/FN).GT.SQRT(1.27735)

AZ32=ABS (RA-ATAN(RA) ) ACZ=AZ32#FN

CZ=-ACZ

Z732=1,.5*AZ32 RTZ=Z52**CON2 FNIS=FN#*CON2 RTARY=RTZ#FN1S ARY=-RTARY#RTARY

60 T0 150

CONTINUE

CASES FOR (X/FN).LT.SQRT (0.7225) AZS2=ABS(ALOG((1.+RA)/XX) -RA)

TEST FOR UNDERFLOW, 1.E-280 = EXP(-644.), ONE WORD LENGTH UP FROM UNDERFLOW LIMIT OF CDC 4600

ACZ=AZ32*FN

CZ=ACZ

IF(ACZ.GT.ELIMN2) GO TO i980 Z32=1.5#AZ32

RTZ=Z232**CON2 FNIS=FN**CON2 RTARY=RTZ#FN15 ARY=RTARY*RTARY

Figure A3. (Sheet 10 of 25)

Al4

BESJ3550 BESJ351 BESJ352 BESJ353 BESJ354 BESJ355 BESJ356 BESJ397 BESJ398 BESJ359 BESJ340 BESJ341 BESJ342 BESJ563 BESJ364 BESJ365 BESUS466 BESJ367 BESJ548 BESJ349 BESJS70 BESJ371 BESJ372 BESJ373 BESJ3574 BESJ375 BESJ376 BESJ377 BESJ378 BESJ379 BESJ380 BESJ381 BESJ382 BESJ383 BESJ3584 BESJ385 BESJ386 BESJ387 BESJ388 BESJS89 BESJS90 BESJ391 BESJ392 BESJ393 BESJ394 BESJ399 BESJ396 BESJ397 BESJS98 BESJ399 BESJ460 BESJ401 BESJ402 BESJ403 BESJ404

oO

150

101

155 145

400

402

CONTINUE PHI=SQRT((RTZ+RTZ) *TAU) TB=1.

ASUM=1.

TFN=TAU/FN UPOL(2)=(C(1,1)*T24+C(2,1)) #TFN RCZ=CONI/C2Z CRZ32=CONS48*RCZ BSUM=UPOL (2) +CRZ32 RELB=TOL¥ABS (BSUM) AP=TFN

KS=0

KP1=2

RZDEN=RCZ

DO 155 LR=2,8,2

COMPUTE TWO U POLYNOMIALS FOR NEXT A(ZETA) AND B(ZETA)

LRPI=LR+1 DO 101 K=LR,LRP1

KS=KS+1

KPL=KP1+1

S1=C(1,K)

DO 102 J=2,KP1 S1=S1*T2+C(J,K)

CONTINUE

AP=AP*TFN

UPOL(KP1)=AP*S1

CR(KS) =BR(KS) #RZDEN RZDEN=RZDEN*RCZ DR(KS)=AR(KS) #RZDEN CONTINUE

SUMA=UPOL (LRP1) SUMB=UPOL (LR+2) +UPOL (LRP1) #CRZ32 JU=LRP1

DO 151 JR=1,LR

JU=JU-1 SUMA=SUMA+CR (JR) #UPOL (JU) SUMB=SUNB+DR(JR) #UPOL (JU) CONTINUE

TB=-TB

IF(W2.GT.0.) TB=ABS(TB) ASUM=ASUM+SUMA*TB BSUM=BSUM+SUMB#TB

IF (ABS(SUMA).LE. TOL. AND. ABS(SUMNB).LE.RELB) GO TO 165 CONTINUE

TB=RTARY

IF(W2.6T.0.) TB=-TB BSUM=BSUN/TH

CONTINUE CALL JAIRY(ARY,RTARY,ACZ,AI,DAI) TEMP (1S)=PHI* (AI *ASUM+DAI*BSUM) /FNIS GO TO (401,202,650), IS

TEMP (1) =TEMP (3)

Figure A3. (Sheet 11 of 25)

Al5

BESJ405 BESJ406 BESJ407 BESJ408 BESJ409 BESJ410 BESJ411 BESJ412 BESJ413 BESJ414 BESJ4195 BESJ416 BESJ417 BESJ418 BESJ419 BESJ 420 BESJ421 BESJ422 BESJ423 BESJ424 BESJ425 BESJ426 BESJ427 BESJ428 BESJ429 BESJ430 BESJ431 BESJ432 BESJ433 BESJ434 BESJ4355 BESJ434 BESJ437 BESJ438 BESJ439 BESJ440 BESJ441 BESJ442 BESJ443 BESJ444 BESJ445 BESJ446 BESJ447 BESJ448 BESJ449 BESJ450 BESJ451 BESJ452 BESJ4535 BESJ454 BESJ455 BESJ456 BESJ4397 BESJ458 BESJ499

oO

oa

ET=1 BESJ440

401 I8=2 BESJ441 DFN=DFN-1.D+0 BESJ462 FN=DEN BESJ443 60 TO 130 BESI444

BESJ445 SERIES FOR (X/2)#**2.LE.NU+1 BESJ446 BESJ447

134 CONTINUE BESJ448 GLN=GAMLN (FNP1) BESJ449 ARG=FN*XO2L-GLN BESJ470 IF(ARG.LT.-ELIM1) GO TO 123 BESJ471 EARG=EXP (ARG) BESJ472

300 CONTINUE BESJ473 §sil. BESJ474 AK=3. BESJ475 eh. BESJ474 Tain BESJ477 S1=FN BESJ478 MO 123 Reilly BESJ479 $2=T2+51 BESJ480 T=-T#5X02/52 BESJ481 S=S+T BESJ482 IF(ABS(T).LT.TOL) GO TO 127 BESJ495 T2=T2+AK BESJ484 AK=AK+2. BESJ485 S1=S1+FN BESJ486

125 CONTINUE BESJ487

7 CONTINUE BESJ488 TEMP (1S) =S*EARG BESJ499 60 TO (301,202,600), IS BESJ490

301 EARG=EARG*FN/ X02 BESJ491 DFN=DFN-1.D+0 BESJ492 FN=DFN BESJ493 1S=2 BESJ494 GO TO 360 BESJ495

BESJ496 SET UNDERFLOW VALUE AND UPDATE PARAMETERS BESJ497 BESJ498

180 Y(NN)=0, BESJ499 NN=NN-1 BESJ500 DFN=DFN-1.D+0 BESJSO1 FN=DFN BESJ502 IF (NN-1) 170,171,130 BESJ 503

17 KT=2 BESJ504 15=2 BESJ505 GO TO 130 BESJ504

123 Y(NN)=0. BESJS07 NN=NN-1 BESJSO8 FNPL=FN BESJ509 DFN=DFN-1.D+0 BESJS10 FN=DFN BESJ511 IF(NN-1) 170,172,173 BESJS12

IW2MKeT=2 BESJ5S13 1S=2 BESJS14

Figure A3. (Sheet 12 of 25)

Al6

oa

aInadnn

170

ioe) P=) i)

300

ott 512

IF(SX02.LE.FNP1) GO TO i33 60 TO 130 ARG=ARG-XO2L+ALOG (FNPI1) IF(ARG.LT.-ELIM1) GO TO 123 60 TO 134

NZ=N-NN

RETURN

BACKWARD RECURSION SECTION

CONTINUE

NZ=N-NN

IF(KT.EQ@.2) GO TO 250 CONTINUE

BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA ¥Y<(NN) =TEMP (1)

Y(NN-1) =TEMP (2) IF(NN.EQ@.2) RETURN DX=X

TRX=2.D+0/DX DTIM=DFN#TRX

TM=DTM

K=NN+1

DO 230 I[=3,NN

K=K-1 ¥(K=-2) =TM#Y (K-21) -Y¥ CK) DTM=DTM-TRX

TM=DTM

CONTINUE

RETURN

Y(1)=TEMP (2)

RETURN

ASYMPTOTIC EXPANSION FOR X TO INFINITY WITH FORWARD RECURSION IN OSCILLATORY REGION X.GT.MAX(29, NU), PROVIDED THE LAST MEMBER

OF THE SEQUENCE IS ALSO IN THE REGION.

CONTINUE IN=ALPHA-TAU+2. IF(IN.LE.O) GO TO 562

INPI=IN+1 DALPHA=ALPHA-FLOAT(INP1) KT=1

GO TO 5i1

DALPHA=ALPHA

IN=06

IS=KT

ARG=X-PIDT#DALPHA-PDF SA=SIN(ARG)

SB=COS (ARG) RA=RTTP/RTX

ETX=8. #X

DX=DALPHA

DX=DX+DX

DTM=DX*DX

Figure A3. (Sheet 13 of 25)

Al7

BESJS15 BESJS16 BESJS17 BESJ518 BESJS19 BESJ320 BESJ921 BESJG22 BESJS23 BESJS2

BESJS29 BESJ526 BESJ927 BESJ528 BESJ529 BESJS3

BESJSS1 BESJ552 BESJ933 BESJ534 BESJS35 BESJ536 BESJ537 BESJSS

BESJ539 BESJS540 BESJ541 BESJ542 BESJ343 BESJ544 BESJS45 BESJS46 BESJ347 BESJS548 BESJ549 BESJ590 BESJS51 BESJ592 BESJS53 BESJS54 BESJS95 BESJ956 BESJ597 BESJ598 BESJ399 BESJ560 BESJS61 BESJ562 BESJS563 BESJ964 BESJ969 BESJ566 BESJ567 BESJ568 BESJ569

ioe)

ao

on oa ~~

T2=DTM-1.D+0 T2=T2/ETX $2=T2 RELB=TOL*ABS(T2) TL=ETX

S1=1.

FN=1,

AK=8.

DO 504 K=1,13 TI=TI+ETX FN=FN+AK DX=FN TRX=DTM-DX AP=TRX T2=-T2*AP/T1 S1=S1+T2 TLSTI+ETX AK=AK+8. FN=FN+AK DX=FN TRX=DTH-DX AP=TRX

T2= T2*AP/T1 §2=52+T2 IF(ABS(T2).LE.RELB) GO TO 505 AK=AK+8.

4 CONTINUE 3 TEMP(1S)=RA*(S1*SB-S2*SA)

60 TO (506,507) ,1S DALPHA=DALPHA+1. 15=2

TB=SA

SA=-SB

SB=TB

60 TO 303

FORWARD RECURSION SECTION

IF(KT.E8.2) GO TO 250 S1=TEMP (1)

S2=TEMP (2)

TX=2./X

TM=DALPHA*TX IF(IN.E@.0) GO TO 520

FORWARD RECUR TO INDEX ALPHA

DO S10 I=1,1N

=52

S2=TM*52-S1

TM=TM+TX

S$1=5

CONTINUE

IF(NN.EQ@.i1) GO TO 535

$=82

Figure A3.

(Sheet 14 of 25)

A18

BESJ570 BESJS71 BESJ572 BESJS73 BESJS74 BESJS75 BESJS7

BESJS77 BESJS78 BESJS79 BESJ580 BESJ581 BESJ582 BESJS83 BESJ584 BESJS85 BESJS86 BESJS87 BESJS88 BESJS89 BESJS90 BESJS91 BESJS92 BESJS953 BESJ594 BESJS95 BESJ596 BESJS97 BESJ598 BESJS99 BESJ400 BESJ401 BESJ4602 BESJ603 BESJ604 BESJ605 BESJ4606 BESJ407 BESJG608 BESJ609 BESJ410 BESJ411 BESJé12 BESJé13 BESJ414 BESJ415 BESJ616 BESJ417 BESJ4618 BESJ419 BESJ620 BESJ421 BESJ622 BESJ623 BESJ624

aInnon

920

600

6950

677

603

S2=TH#S2-S1 TM=TH+TX 51=5 CONTINUE

FORWARD RECUR FROM INDEX ALPHA TO ALPHA+N-1

Y(t) =S1

¥(2)=52

IF(NN.E@.2) RETURN DO 530 I=3,NN Y¥(T)=TM*¥Y (1-1) (1-2) TM=TM+TX

CONTINUE

RETURN

Y¥(1)=S2

RETURN

BACKWARD RECURSION WITH NORMALIZATION BY ASYMPTOTIC EXPANSION FOR NU TO INFINITY OR POWER SERIES.

CONTINUE

COMPUTATION OF LAST ORDER FOR SERIES NORMALIZATION KM=AMAX1(S.-FN,O.)

TFN=FN+FLOAT (KM)

TA=X02L-TA

TB=-(1.-1.5/TFN)/TFN IN=CE/(-TA+SQRT (TAX TA-CE*¥TH))+1.5 IN=IN+KN

GO TO 603

CONTINUE

COMPUTATION OF LAST ORDER FOR ASYMPTOTIC EXPANSION NORMALIZATION

GLN=AZ32+RA

IF(ARY.GT.30.) GO TO 475 RDEN=(PF(4) #ARY+PP(3)) *ARY+1. RZDEN=PP (1) +PP (2) #ARY TA=RZDEN/RDEN

IF(W2.LT.G.16) GO TO 451 TB=GLN/RTARY

60 TG 677

TB=(1. 259992104940. 16798947 30#W2) /FNLS GO 10 677

CONTINUE

TA=CONL*¥TCE/ACZ

TA=((0, O4935827160#TA-O. LLLLLLLL11) #TAtG. 6666666667) €TAXARY IF(W2.LT.0.10) GO TO 451 TB=GLN/RTARY

IN=TA/TB+1.5

IFC IN. GT.INLIM) GO TO 4062 DX=FLOATC(IN)

DTM=DFN+DX

DX=X

TRX=2.D+0/Dk

DTM=DTM#TRX

Figure A3. (Sheet 15 of 25)

Al19

BESJ625 BESJ626 BESJ627 BESJ628 BESJ429 BESJ630 BESJ631 BESJ652 BESJ433 BESJ634 BESJ6355 BESJ436 BESJ637 BESJ438 BESJ639 BESJ4640 BESJ641 BESJ642 BESJ643 BESJ644 BESJ649 BESJ4646 BESJ647 BESJ648 BESJ649 BESJ4650 BESJ651 BESJ652 BESJ695 BESJ654 BESJ655 BESJ696 BESJ4&57 RESJ658 BESJ659 BESJ460 BESJ661 BESJ642 BESJ6635 BESJ664 BESJ6465 BESJ466 BESJ667 BESJ468 BESJ469 BESJ670 BESJ471 BESJ6/72 BESJ473 BESJ4674 BESJ4735 BESJ676 BESJ677 BESJ478 BESJ479

aod

a

605

601

604 615

oy

TM=DTM TA=0. TB=TOL KK=1 CONTINUE

BACKWARD RECUR UNINDEXED

DO 601 I=1,1N

S=TH

TB=TM*TE-TA

TA=S

DTM=DTM-TRX

TM=DTM

CONTINUE NORMALIZATION IF(KE.NE.1) GO TO 4604 TA=(TA/TB) #TEMP (3) TB=TEMP (3)

KK=2

IN=NS

IF(NS.NE.0) GO TO 605 Y (NN) =TB

NZ=N-NN IF(NN.E@.1) RETURN 5=TB

TB=TM#TB-TA

TA=S

DTM=DTM-TRX

TM=DTM

K=NN-1

Y(K)=TB

TF (NN.EQ@.2) RETURN KM=K-1

BACKWARD RECUR INDEXED

DO 602 I=1,KM Y(K-1) =TM#Y (K) (K#1) DTM=DTM-TRX

TM=DTM

K=K-1

CONTINUE

RETURN

CONTINUE NZ=-2 RETURN CONTINUE NZ=-3 RETURN CONTINUE NZ=-1

Figure A3.

(Sheet 16 of 25)

A20

BESJ4580 BESJ681 BESJ4682 BESJ483 BESJ4684 BESJS85 BESJ4686 BESJ487 BESJ4688 BESJ689 BESJ490 BESJ4691 BESJ492 BESJ693 BESJ4694 BESJ495 BESJ694 BESJ697 BESJ498 BESJ699 BESJ700 BESJ701 BESJ702 BESJ703 BESJ704 BESJ705 BESJ706 BESJ707 BESJ708 BESJ709 BESJ710 BESJ711 BESJ712 BESJ713 BESJ714 BESJ715 BESJ716 BESJ717 BESJ718 BESJ719 BESJ720 BESJ721 BESJ722 BESJ723 BESJ724 BESJ725 BESJ726 BESJ727 BESJ728 BESJ729 BESJ730 RBESJ731 BESJ732 BESJ7335 BESJ734

NOaonnrorinrwarIiononiwaig

ioe)

CDC 6600 ROUTINE

AND ITS DERIVATIVE DAI(X) FOR JBESS INPUT: X - ARGUMENT, COMPUTED BY JBESS, X UNRESTRICTED RX RX=SORT(ABS(X)), COMPUTED BY JBESS C - C=2.*#(ABS(X)*#1.5)/3., COMPUTED BY JBESS OUTPUT: AI - VALUE OF FUNCTION AI(X) DAI - VALUE OF THE DERIVATIVE DAI(X) WRITTEN BY D.E. AMOS, S.L. DANIEL & M.K WESTON DIMENSION AKI (14) ,AK2(23) ,AKS (14) DIMENSION AJP(19) sAUN(19) ,AC15) ,B(15) DIMENSION DAK1 (14) ,DAK2(24) ,DAK3 (14) DIMENSION DAJP(19) ,DAIN(19) ,DA(15) ,DB(15) DATA Ni,N2,N5,N4/14,23,19,15/ DATA M1,M2,M3,M4/12,21, 17,13 / DATA FP112,CON1,CON2,CONS,CON4 ,CONS/

1

5 <

Cl roe

Poe 10 CO SO ON fe C4 PO NO oS

fy PQ oe

“Nc on eS

1

RETURN END

SUBROUTINE JAIRY(X,RX,C,Al,DAI)

1.5089969389958E+00, 3.8000458986729E-01,

DATA AKI(1) AK1(2) /-1. AK1(4) / 8. AK1(6) / 1. AK1(B) /-1

DATA AK2(1) AK2(2) / AK2(4) / AK2(6) / AK2(8) / AK2(10)/ AK2(12)/ AK2(14)/ AK2(16)/ AK2(18)/ AK2(20)/ AK2(22)7

DATA AK3(1) AKS (2) AKS(4) AK3 (4) AKS (8)

AK3(14)/ DATA AJP(1) AJP (2)

98762 »8729020974102E- 06/, »21135041695591E-07/,

29245470528 265E- 097, »26863505659574E-10/, ~JOL41035901371E-11/, ©§4464450473454E-12/,

»7010603055935E- -1977479916481E-14/,

144734079E-

et et oO ee lO) Rl ed oo CO I |

HON fe cal ike Laie AK3(10)/-6. AR S(L2)/-1.

fHIkq

RASC]

JAIRY COMPUTES THE AIRY FUNCTION AI(X)

/ 2. 2042309098779E- 2929024278770E-01/7, 2284415200554E-04/, 6382428017212E-05/7,

»29621999359335E-07/, AKL(190)/ 1. AKL(12)/ 2. AK1(14)/-35.

Rese see ae VASES g O32535257423563E-12/, 16169897 78508E- Maid / 2.7436615086960E-

3979096975 F690E-03/,

2924875E- o4/, 1713 1890E-05/,

13/,

/ 2.8027 7812704284438E-03/, b3524996526900E-06/, D2 29435502291 6E-09/, 2440825180026E-11/, S63752444639351E-135/, 1270513469106E-14/,

-2.46480324312435E-16/ 4 7.7895295643758E-

B43565654568GE-01/,

Figure A3.

b. SEBO ONRIO TES Ol,

WII.

de die dr din died

01/,

AKL (3) AK1 (3) AK1(7) AK1 (9) AKI (11)/-

SS

1. 20 3. 8. 3.

AKI(13)/-1.

O1/, AK2 (3) AK2(5) AK2(7) AK 2 (9)

AK2(13)/ AK2(15)/-1. AK2(17)/-9,

AK2(19)/-4.

AK2(21)/-4, AK2(23)/-3. OLA AK3(3) / AK3(5) / ARSC) if AK3(9) / AK3(11)/ AK3(13)/

02/, AIP(S) /

fio (No JI 5 7-4. AK2(11)/-3. -2,0390244716791E-

cO cn cn ~O

9.03515471619678E+00, 8. 6602540378444E-01/

OSBBLL6335919E-02/, 34614345891 23E-04/

0690258957319E-07/, 2290815882567E-09/, 3916546561568E-11/, 1067954609788E-14/,

5733922062119E-03/, 1212491739992E-04/, 3680422537055E-06/, 7589279396229E-07/, 0924537427061E-08/, 09/~" 3670476763957E-10/, 3139829654843E- 12/, 4384026199096E-13/, 5076010450328E-14/, 1907704086507E-15/

© OS542257962900E-05/, i pivomiaeaneanes 08/, »4713840457655E-10/7, » 60477904117 3128576196625

AEST BA o E-14/, 2267976598135E-15/,

~O1412605214617E-02/7,

(Sheet 17 of 25)

A21

BESJ7395 BESJ736 AIRY1 AIRY2 AIRYS AIRY4 AIRYS AIRY6 AIRY7 AIRY8 AIRY? AIRY10 AIRY11 AIRY12 AIRY13 AIRY14 AIRY195 AIRY16 AITRY17 AIRY18 AIRY19 AIRY 20 AIRY21 AIRY22 AIRY25 AIRY 24 AIRY29 AIRY26 AIRY27 AIRY28 AIRY29 AIRY30 AIRY31 AIRY32 AIRY33 AIRY 34 AIRY35 AIRY36 AIRY37 AIRY38 AIRY39 AITRY40 AIRY41 AIRY42 AIRY43 AITRY44 AIRY45 AIRY46 AIRY 47 AIRY 48 AIRY49 AIRYSO AIRYS1 AIRYS2 AIRYS3

oOo ON Oo OF & hd by

Oo On OO Pei pe

“SM & OF f hd Pore

—“SoCe On B YI Po re

(Ea) (23) C>

Noo Se

m= 10 CO NO CN Fs ty Por

AJP(4) / 3.0534272427761E-02/, AJP(S) /-4.9542470251308E-03/, AJP (4) /-1.7274955256395E-03/, AIP(7) / 2.4313763783919E-04/, AJP(8) / 5.0456477751708E-05/, AJP(9) /-6.16316582469521E-06/, AJP (10) /-9.0398574551077E-07/, AIP(11)/ 9.7024377835588E-08/, AJP(12)/ 1.0963945330520E-08/, AJP(13)/-1.0471633058877E-09/, AJP (14) /-9.6035944134465E-11/, AJP(15)/ §.2535878945413E-12/, AJP (16)/ 6.3612343901877E-13/, AIP(17)/-4.9662961411602E- 14/, AJP (18) /-3.2981028892962E-15/, AJP(19)/ 2.3579825203110E-16/ DATA AJN(1) / 3.8049788761724E-02/, AJN(2) /-2.4531954184555E-01/, AJN(3) / 1.6582062370270E-01/, AIN(4) / 7.4933004581879E-02/, AJN(S) /-2.6347628810664E-02/, AIN(5) /-5.9253559730498E-03/, AIN(7) / 1.4474440958980E-03/, AIN(8) / 2.1831183132222 -04/, AIN(9) /-4.1065207748030E-05/, AIN(10) /-4.6687499417177E-06/, AIN(L1)/ 7.1521880727714E-07/, AIN(12)/ 6.5296477085463E-08/, AIN(13)/-8.4428402756595E-09/, AJN(14) /-6.4419615897698E-10/, AIN(15)/ 7.20802 28650528E-11/ AIN(14)/ 4.7246543171785E-12/, AJN(17) /-4.6602263254704E- 137, AIN(18) /-2.6776271038919E-14/, AJN(19)/ 2, 3616131657002E-15 DATA A(1) / 4.9027542474279E- D1/, A(2) / 1 .5764727794620E-03/, A(3) /-9.6619596314031E-05/, A(4) / 1.3591608028882E-07/, A(5) / 2.9815734265486E-07/, Ald) /-1.868247675599B8E-08/, A(7) /-1.0368573766714E-09/, A(B) / 3.28660818434335E-10/, A(9) /-2.570914106327BE-11/, A(10)/-2.32357655300468E-12/, A(11)/ 9.5752327904826E-13/, A(12)/-1.2034082864972E-13/, A(13) /-2.9090771677072E-15/, A(14)/ 4.5565445458015E-15/, A(15) /-9.9900387481026E-16/ DATA B(1) / 2.7859355280308E-01/, Bi2) /-3. Seiteey ieee 53E-03 B(3) /-2.3114967738499E-05/, B(4) / 4.7131784226356E- nee B(S) /-1.1241590793133E-07/, B(6) /-2. NOLOOSDL AREAS, B(7) / 2.6094807530219E-09/, B(8) /-3.5509813610122E-11/, B(9) /-3.508499784238B8E-11/, B(10)/ 5.8300718795420E-12/, B(11)/-2.0464482875333E-13/, H(12)/-1.1052917947674E-13/, B(13)/ 2.9772477803878E-14/, B(14)/-2.8820511100994E-15/, B(15)/-3.3265631169617E-16/ DATA NID,N2D,N3D,N4D/14,24,19,15/ DATA M1D,M2D,M3D,M4D/12,22,17,13/ DATA DAK1(1) / 2.0456784230789E-01/, DAKi(2) /-6.6132273990566E-02/, DAKI(3) /-8.4984590098929E-03/, DAKL(4) / 3.1218349155629E-03/, DAKL(S) /-2.7001648982943E-04/, DAK1(5) /-6.3563629867939E- 08/4 DAKI(7) / 3.0239771240951E-06/, DAK1(8) /-2.1831119533009E-07/, DAKI(9) /-5.3619428933293E-10/, DAKIi(10)/ 1.1309803562231E- ae DAK1 (11) /-7.4302383462907E-11/, DAK1(12)/ 4.2890417082689E-13/, DAKL{13)/ 2.2381092575454E-13/, DAKI (14) /-1, 3914013564118E- 14/ DATA DAK2(1) / 2.9333234388323E-01/, DAK2(2) /-8.0619878474311E-03 3/, DAK2(3) / 2.42540172333514E-03/, DAK2(4) /-.8229754885024E-04/, DAK2(5) / ines 8642775118E-04/, DAK2(6) /-4.9745744768406E-05/, DAK2(7) / 1.3209068123950E-05/, DAK2(8) /-3.4952824044494E-06/, DAK2(9) / 9.2436245107894E-07/, DAK2(10)/-2.4473257152187E-07/, DAK2(11)/ 6.4930783764891E-08/, DAK 2 (12) /-1.7271762150154E-09/, DAK2(13)/ 4.6072576360466E-09/, DAK2 (14) /-1.2324905529155E-09/, DAK2(15)/ 3.3062040949810E-10/, DAK2(15)/-8.8925209977240E-11/, DAK2(17)/ 2.3977351987930E-11/, DAK2 (19) /-S.4801392115345E-12/, DAK2(19)/ 1.7551013202373E-12/, DAK2(20) /-4. 7630382983344E-13/, DAK2(21)/ i. 29498241 10081E-13/,

Figure A3.

(Sheet 18 of 25)

A22

AIRYS4 AIRYS9 AIRYS6 AIRYS7 AIRYS8 AIRYS9 AIRY60 AIRY61 AIRY62 AIRY6$ AIRY 64 AIRY45 AIRY46 AIRY67 AIRY48 AIRY49 AITRY70 AIRY71 AIRY72 AIRY73 AIRY 74 AIRY735 AIRY76 AIRY77 AIRY78 AIRY79 ATRY80 AIRY81 ATRY82 AIRY83 AIRY84 AIRY85 AITRY86 AIRYS7 AIRY8E& AIRY89 AITRY99 AIRY91 AIRY92 AIRY93 AIRY?4 AIRY9S AIRY96 AIRY97 AIRY98 AIRY99 AIRY100 AIRY101 AIRY162 AIRY103 AIRY104 AIRY105 AIRY106 AIRY107 AIRY108

Ld bo

oe to a Su ee ee

Oo wo SO OF S fy Po ee

oon oO Pw Pe

SWOOP

“1 oo~ CN & te Poor

DAK2 (22) /-3.5267 -2.62786914342 / 2.8467582881135E-01/7,

DAK2(24)/ DATA DAKS(1 DAKS(2) / DAKS(4) 7 DAKS(6) / DAKS(8) / DAKS(10)/ DAK3S (12)/ DAKS(14)/ DATA DAJP(1 DAJP(2) /- DAJP(4) / DAJP(6) /- DAJP(8) / DAJP (10) /- DAJP(12)/

DAJP(14)/-3.,

DAJP(16)/ DAJP(18)/ DATA DAUN(1 DAJN(2) / DAJN(4) /- DAIN(4) / DAJN(8) /- DAJN(1G)/ DAJN(12)/- DAJN(14)/ DAJN(16) /- DAIN(18) / DATA DAI)

DA(3)

221043E-14/, 29E-15/

»9507307261908E-03/, »84907283944634E-06/, »0592563445715E-69/, ~ 563557688851 354E-11/, »8857435378444E-13/, »1737476261721E-14/, »0957477309706E-16/ ) ff 1.20262933468882E-01/, 1.6794842923050E-02/, 8.4556029509887E-04/, 2.2982786094548E-05/, 3.76343991134692E-07/, 4, 29611335 200301E-09/7, S7245881346190E-11/, »2612065309577E-13/,

ep

boo co on

2 -1.1260437448512E-15/, )

7 1.0859453963297E-02 8.553513197485709E-02/, 8.7842072529426E-02/, 9,4167406050324E-63/, 4. 1115734315683E-047, 9,8763368220840E-06/, 1.507985001S147E-07/, 1.5959991741922E-09/, 1. 24105S55050030E-11/, 7.35940097115574E-14/,

2514076235408E

4.9162732110460E-017,

DAK2(23)/ 9.

DAK3S (3) DAKS(5) DAKS(7) /-5. DAK3(9) 7-3. DAK S(11)/-8. DAK3(13)/-1.

i-4, POM a

6.952191 35151 146E-027,

DAIP(3) DAJP (3S) DAJP (7) DAJP (9) DAJP(11)/ 3. DAJP(13)/-3. DAJP(15)/ 2. DAJP(17)/-1. DAJF(19)/ 7. Af 6

DAIJN(3) /-3. DAJN(S) / 5. DAJN(7) /-3.

DAJN(9) / 1.

DAJN(11)/-1. DAJN(135)/ 2. DAJN(15)/-2. DAJN(17)/ 1.

DAJN(19)/-7. DA(2) / 3. DA(4) 7-4.

Wg foie Yo

[oH

/ / DAS) / DA(7) /- DIAG DALI) /- DA(13) /- DA(15) DATA DB(1) DB(S) DB(7) DBC?) /-

Ps Lode) [F(G. 61.5.) LF (X.6T. 1.2) T=(X+X-1, Way ey J=N1 FL=AK1 (J) F2=

60 60

2

i

i=

DB(3) /-8.4 i 8

YOR

4

-05/,

»1SL5888055463E-087, SOIT Tldor 226097,

S57 2404342042E-11/,

3. 2408911983032E-13/, »4G75524797406E-14/, »1790078647740E-16/

HAY NS

252852

SO944235E-O1/, 219009E-05/7,

»4238972021762E-07/, »3637783584458E ~40995S691G5B19E-12/, DIB, aii. DB CUS) (ot DB(15)/-7.1

-09/,

S48 7401513 702E-i13/, 7725354306791 1E-147, 179333575 79S53E-16/

TO 206

60 TO 150 2) *CON4

Figure A3.

A23

DA(&) / DA(8) /- DA(10)/- DA(12)/ DA(14)/

DB(2) / DB(4) # DE(6) /- DB(8) / DB(10)/- DB(i2)/- DB(14)/

oe ike Sie Ly

les

fi =e 6. 5. 4, ie

6200315158592E-15/,

8348115033798E-05/, 0141849117858E-07/, BS32529140058E-10/, 9088909477950E-12/, 6858825645219E-14/, 68523514651092E-15/,

»7801023526382E-03/,

97146140182135E-03/, 4288962070198E-05/,

2.2906787091599E-G6/,

456639353555956E-08/, S8673569121499E-10/, 72696091 06654E-12/, S876S205235850E-14/, SCS TAS 17 7

527706811 306E-01/, 325190697605E-02/, 3218702601900E-03/, 129732689135E-04/, 9731296981239E-06/, 3268766952539E- OB/, 076659226683B8E-10/, 39631765351 04E-12/, 3288747562750E-15/ 11164930827496-03/, 617697761721 4E- 9729580465652E-08/, 4475282664204E-10/, 9965506584722E-i2/, 6209895256874E-13/, 5938481128449E-16/,

1 sa] 0

06/,

4421 2835341992E-05/, SBO40S51841871E-04/, 2428689470978E-09/, »1699104265667E-10/, 1867422109358E-12/, 9OL9057660871E-14/, 5595061044266E-15/,

(Sheet 19 of 25)

AIRY109 AIRY116 AIRY11i1 AIRY1L12 AIRY113 AIRY114 AIRY115 AIRY116 AIRY117 AIRY118 AIRY119 AIRY120 AIRY121 AIRY122 AIRY123 AIRY124 AIRY125 AIRY126 AIRY127 AIRY128 AIRY129 AIRY130 AIRY131 BEJS132 AIRY133 AIRY134 AIRY135 AIRY1S6 AITRY137 AIRY137 AIRY139 AIRY140 AIRY141 AIRY142 AIRY143 AIRY144 AIRY145 AIRY146 ATRY147 AIRY148 AIRY149 AIRY150 AIRY151 AIRY152 AIRY1935 AIRY154 AIRY1595 AIRY196 AIRY157 AIRY158 AIRY1959 AIRY140 AIRY141 AIRY162 AIRY16$

1035

106

DO 105 1=1,h1 J=J-1

TEMPL=F 1 FI=TT¥F1-F2+AK1 (J) F2=TEMP1

CONTINUE AL=T#F1-F2+AK1 (1)

J=N1D

F1=DAK1 (J) F2=0.

DO 106 I=1,M1D J=J-1

TEMP1=F1

FL=TT¥F1-F2+DAK1 (J)

F2=TEMP1 CONTINUE

DAI=-(T#F1-F2+DAK1 (1) )

RETURN

CONTINUE T=(X+X-CON2) ¥CONS We Sul kee L

J=N2

FL=AK2(J)

F2=0.

pO 155 I=1,M2 J=d-1

TEMP1=Fl FL=TT#FI-F2+AK2 (3) F2=TEMP 1

CONTINUE RTRX=SORT (RX) EC=EXP(-C)

AI=EC¥(T#FI-F2+AK2

J=N2D FL=DAK2(d) F2=0.

DO 156 1=1,M2D J=J-1

TEMP1=F 1

(1))/RTRX

FI=TT¥FI-F24DAK2 (0)

F2=TEMP 1 CONTINUE

DAL=-EC#(T#F1-F2+DAK2(1)) #RTRX

RETURN

CONTINUE T=10. /C-1. TT=T+T

J=N1 FI=AK3(J) F2=0.

DO 205 I=1,M1 J=J-1

Figure A3.

(Sheet 20 of 25)

A24

AIRY164 AIRY165 AIRY166 AIRY147 AIRY168 AIRY149 AIRY179 AIRY171 AIRY172 AIRY173 AIRY174 AIRY175 AIRY176 AIRY177 AIRY178 AIRY179 AIRY180 AIRY181 AIRY182 AIRY183 AIRY184 AIRY185 AIRY186 AIRY187 AIRY188 AIRY189 AIRY190 AIRY191 AIRY192 AIRY193S AIRY194 AIRY195 AIRY196 AIRY197 AIRYi98 AIRY199 AIRY200 AIRY201 AIRY202 AIRY203 AIRY204 AIRY205 ALRY2Z06 AIRY 207 AIRY208 AIRY209 AITRY210 AIRY211 AIRY212 AIRY213 AIRY214 AIRY215 AIRY216 AIRY217 AIRY218

ba oS cn

300

309

306

TEMP1=F 1

FL=TT#F1-F2+AK3 (J)

F2=TEMP1

CONTINUE

RTRX=SQRT (RX)

EC=EXP(-C) AI=EC#(THFI-F2+AKS(1))/RTRK J=N1D

FL=DAKS(J)

F2=0,

DO 206 I=1,MiD

J=J-1

TEMP1=F1 FL=TT#F1-F2+DAKS (J) F2=TEMPI

CONTINUE DAI=-RTRX#EC#(T#FI-F2+DAKS(1)) RETURN

CONTINUE IF(C.G67T.5.) 60 TO 350 T=.44#C-1.

TT=T+T

J=N3

FL=AJP (J) EL=AUN(J)

F2=0.

E2=0.

DO 305 I=1,M5 J=d-1

TEMP1=F 1

TEMP2=E1 FLSTT#Fi-F2+AUP (J) EL=TT#E1-E2+AUN (J) F2=TEMPI

E2=TEMP2

CONTINUE

AL=(T#EL-E2+AIN(L))-X¥(T#FL-F2+ANP (1) )

J=N3D

FL=DAUP (J) EL=DAIN(J)

F2=0.

E2=0.

DO 306 I=1,M3D

J=J-1

TEMP1=F 1

TEMP2=E1

Fl = TT¥FL-F2+DAJP (J) El= TT#E1-E2+DAIN(J) F2=TEMP1

E2=TEMP2

CONTINUE

DAI=X#X*¥(T#FI-F2+DAJP(1))+(T#E1-E2+DAIN(1))

RETURN

Figure A3.

(Sheet 21 of 25)

A25

AIRY219 AIRY22 AIRY221

LD ms a) ~

Le)

oo

AIRY2351 AIRY232 AIRY233 AIRY234 AIRY235 AIRY23 AIRY237 AIRY23 AIRY239 AITRY240 AIRY241 AIRY242 AIRY243 AIRY244 AIRY245 AIRY246 AIRY 247 AIRY248 AIRY249 ALRY250 AIRY251 AIRY252 AIRY253 AIRY254 AIRY255 AIRY256 AIRY257 AIRY258 AIRY259 AIRY260 AIRY261 AIRY262 AIRY263 AIRY264 AIRY269 AIRY266 AIRY267 AIRY268 AIRY269 AIRY270 AIRY271 AIRY272 AIRY275

oO

(sa) (52) (se) fee) isel Ge) Ga

350

CONTINUE T=10./C-1.

TT=T+T

J=N4

FL=A(d)

E1=B(J)

F2=0,

E2=0.

DO 310 I=1,H4 J=d-1

TEMP1=F 1

TEMP2=E1 FL=TT#F1-F2+A(d) E1=TT¥E1-E2+B(J) F2=TEMP1

E2=TEMP2

CONTINUE TEMP1=T#F1-F2+A (1) TEMP2=T¥EL-E2+B (1) RTRX=SORT (RX) CV=C-FPI12 CCV=COS(CV) SCV=SIN(CV) Al=(TEMP1*¥CCV-TEMP2*SCV) /RTRX

J=N4D

FLi=DA(J) EL=DB(J)

F2=0,

E2=0.

DO 3ii 1=1,M4D J=J-1

TEMP1=F 1 TEMP2=E1

FI=TT#F1-F2+DA(J) EL=TT#E1-E2+DB(d) F2=TEMF1

E2=TEMP2

CONTINUE TEMPL=T#F1-F2+DA(1) TEMP2=T#E1-E2+DB (i) EL=COV*CONS+. S¥SCV E2=SCV*CONS-.S*CCV DAI=(TEMPL*E1-TEMP2*E2) *RTRX RETURN

END

FUNCTION GAMLN(X)

WRITTEN BY D. E. AMOS, SEPTEMBER, 1977. REFERENCES * SAND-77-1518

* COMPUTER APPROXIMATIONS BY J.F.HART, ET.AL., SIAM SERIES IN APPLIED MATHEMATICS, WILEY, 19468, P.135-134.

* NBS HANDBOOK OF MATHEMATICAL FUNCTIONS, AMS S55, BY M. ABRAMOWITZ AND I.A. STEGUN, DECEMBER. 1955, P.257

ABSTRACT Figure A3. (Sheet 22 of 25)

A26

AIRY274 ATRY2735 AIRY276 AIRY277 ALRY27 AIRY279 ATRY280 AIRY281 AIRY282 ATRY283 AIRY284 AIRY285 AIRY286 ATRY287 AIRY288 AIRY289 AIRY290 ATRY291 AITRY292 AIRY293 AIRY294 AIRY295 AIRY296 AIRY297 AITRY298 ATRY299 AIRY300 AIRY3O01 AIRY302 AIRY303 AIRY304 ATRY305 AIRY306 AIRY307 AIRY308 AITRY309 ATRY210 ATRYS11 AIRYS12 ATRY313 AIRY314 AIRYS15 AIRYS14 AIRYS17 AIRY318 GLNI GLN2 GLNS GLN4 GLNS GLN6 GLN7 GLNS GLN9 GLN1LO

NIOOIAINIAAMAAWAAAABAAHANHM

GAMLN COMPUTES THE NATURAL LOG OF THE GAMMA FUNCTION FOR

X.GT.G. A RATIONAL CHEBYSHEY APPROXIMATION IS USED ON B.LT.X.LT.1000., THE ASYMPTOTIC EXPANSION FOR X.SE.1000. A RATIONAL CHEBYSHEV APPROXIMATION ON 2.LT.X.LT.3. FOR O.LT.X.LT.8. AND X NON-INTEGRAL, FORWARD OR BACKWARD RECURSION FILLS IN THE INTERVALS 0.LT.X.LT.2 AND S.LT.X.LT.8. FOR X=1.,2 3100., GAMLN IS SET TO NATURAL LOGS OF FACTORIALS. DESCRIPTION OF ARGUMENTS INPUT X SeXGh 0 OUTPUT GAMLN - NATURAL LOG OF THE GAMMA FUNCTION AT xX DIMENSION GLN(100) ,P(S) ,Q(2) ,PCOE(9) ,Q@COE(4) DATA XLIM1,XLIM2,RTWPIL/ EG 9 aoe » 97.189385322047E- DATA P(1)/7.663451880000E-04/, P(2)/-5.940956105200E-04/, CSS) / To PESTS NOSE TSER O47, PCA) /-2.7777777546577E-03/,

i

2 DATA DATA

PW he

ted Poo

oon oon Pine Oona one

DATA

Cn ee tel Pa oe

RCo ies Q(1)/- PCOE(1)/2.

PCOE(3)/1. PCOE(S)/2, PCOE(7)/2. PCOE(9) /6.

OCOE(1)/1. QCOE(3)/9.

GLN(1) /0. GLN(4) /1., GLN(4) /4, GLN(8) /8.52 GLN(10)/1. GENCI2) 7107 GLN(14)/2.2 GLN(16)/2, GLN(18)/3.3 GLN(20)/3., 2)/4,.5 GLN(24)/5. GLN(26)/5. WAG 4 GLN(30)/7. GLN(32)/7., GLN(34)/8. GLN(34)/9, Nao Fe GLN(40)/1.,

GLN(2

GLN(28

SLN(38

BLN(42)/1.

SLN(44)/1.

GLN(45)/1.

GLN(48)/1., "445657439463E+02/,

(Sheet 23 of 25)

BLN(S0) /1

edrdend

ITIIIIIIT

-02/

7777777777786 - O3/,

973786644810E-03/ O931159546710E-01/ 159943128451E+00/ 0782472535179E+01/ 2003583800713E+01/

DODDDDD00NDOE+00/,

Q(2)/8. 3333535355333 3E- 2581943959028E- 980671310204E-01/7, 338067999387E+00/, 503677233002E+01/,

»PCOE(2)/9. ,PCOE(4)/3. »PCOE(4) /6. »PCOE(8)/3.

QCOE(2

02/ 03/,

)/-8.906016659498E+00/,

822521104714E+00/ ,QCOE (4) 76. 20038380071 3E+O1/

O/7, GLN(2) /0.0/, TILTS94G9228E+00/, TE7AIIT4A27E2E+O0/, es ere OLB2748008E+01/, O2307B84587E+01/, 216385512E+01/, 9927138384E+01/, SSOS0734S014E+01/, FS359BB41B8720E+01/, 38013889848E+01/, 160667554776E+01/, BOG3IE0S22298E+01/, JI7SSB62701E+01/, 12570289671 7E+01/, B09222555332E+01/, JOS44S701758E+01/, 213617560369E+01/, P33051245479E+017, 066317502506E+02/ 14G342117815E+02/, 219330815 1S4E+02/, 29 ee eae

3480272243 73E+

ee ee ies Ie coc Gy

Figure A3.

GLN (3) GLN(3)

GLN(7) /6. GLN(9) /1. GLN(I1)/1. GLN(13)/1. GLN(15)/2. GLN(17)/3.,

GLN(19)/

3.63590445

GLN(21) 7/4, GLN(23)/4.,

GLN (25)

GLN(27)/6. GLN(29)/6, GLN(31)/7.

GLN(33)/8. GLN(35)/8. GLN(37)/9,

GLN{39)/1. GLN(41)/1,

GLN(43)/1. GLN(45)/1.

GLN(47)/1. GLN(49)/1. GLNCS1)/1.

A27

75.931471805599E-01/, /3.178053830348E+00/,

S79251212010E+00/, 050460290274E+01/, S10441257308E+01/, 9987214495466EtO1/, S19122118274E+01/, DS71BS01060BE+91/, Z20B03E+01/, 255561646075E+01/, B471181351B84E+0i/,

/S.A7TB4A7T2Z9S9BLIE+O1/,

261701761 00E+01/, 7TEB974AZLSTLISE+Oi/, Fees Oe e eS Eun.

LSS7ISPASSLZE+O1/, 898082754220E+01/, S71G69454214Et+01/, O29681996145E+02/ LOSZ06397148E+02/, 177718813997E+92/, 205172711 494E+02/, S295Z5750356E+02/, 40673923648 2E+02/, 48477766951 8E+02/,

AND

o1/

GLNi1 GLN12 GLNIS GLN14 GLNIS GLN14 GLN17 GLN18 GLNI9 GLN20 GLN21 GLN22 GLN23 GLN24 GLN25 GLN26 GLN27 5LN28 GLN29 GLNSO GLN3S1 GLNS2 GLN33 GLN34 GLN3S GLNS6 GLN37 GLN38 GLN39 GLN40 GLN41 GLN42 GLN43 GLN44 GLN45 GLN46 GLN47 SLN48 GLN49 GLNSO GLNGi GLNS2 BLNSS GLNS4 6GLNSS GLNS6 SLNS7 GLNS8 GLNO? GLNSO GLN61 GLNé2 GLN6S GLN64 GLNGS

a

Lap)

(val (sa) (o2)

an

10

15

~0 CO “I cr

L4 poe OO coONM co Cn & 4 Po or

oon or cn

DATA

DATA

GLN(S2)/1. GLN(S4)/1. BLN(S4) /1 GLN(S8)/1. GLN(60)/1, GLN(61)/1. GLN(S6S)/1. GLN(65)/2. GLN(67)/2. GLN(69) 7/2. GLNI71)/2. GLN(73)/2. GLN(75)/2. GLN(77)/2. GLN(79)/2. GLN(81)/2. GLN(83)/2. GLN(85)/2. GLN(S7)/3., GLN(89)/3. GLN(91)/3. GLN(93)/3, GLN(95)/3. GLN(97)/3. GLN(99)/3.

S24095925945E+02/, 6035511 2B215SEt+02/,

~68S274454484E+62/,

7639584840 70E+02/, 8453358288614 E+G2/ 8862917542357E+02/, 968661814572Z9E+02/, O51681994826E+02/, LS5322414945E+02/, 219564418191 E+92/, SO439043546S58E+02/, SB97B3895618E+02/, 4§73729140962E+02/, 562211355500E+02/, 649216497985E+02/, 736735124285 7E+02/, B24742926876Et+02/, FL S2Z39500943E+027, OO2Z2074B6470E+02/, U91L641935802E+02/, 18152463596202E+02/, 271852877038E+02/, So2511819792E+02/, 4337940706235E+02/, J45359085S5194E+02/,

NDX=X T=X-F POT

LOAT(NDX) £@.6.0)

DX=XLIMI-X%

IF (DX.

RATIONAL CHEBYSHEY APPROXIMATION ON 2.

LT.0.0)

NXM=NDX-2 PX=PCOE (1)

DO 10 1=2,9 PX=T#PX+PCOE (1) QX=OCDE(1)

DO 15 I=2,4 QX=T#OX+QCOE (1)

DGAM=

IF (NX IF (NX

BACKWARD RECURSION

PX/QX M.GT.9)

MEG. 0)

60 TO 2 60 T0 2

60 TO 51

GO TO 40

DGAM=DGAN/(1.+T)

IF (NX

DGAM=

M.EQ.-1) DGAN/T

60 T

Ac aud

GAMLN=ALOG (DGARM) RETURN

FORWARD RECURSION FOR 3.LT.X.LT.9

Figure A3.

A28

PON Wolbiokalbiok

GLN

GLN(62)/1.

GLN(54)/2

GLN(646)/2. GLN(48)/2. GLN(70)/2,

GLN(72)/2.

GLN(74) 7/2. GLN(76)/2. GLN(78)/2. GLN(S0)/2.

GLN(82)/2. GLN(84)/2. GLN(84)/2,

SLN(88)/3. GLN(90)/3. GLN(92)/3. GLN(94)/3.

GLN(96)/3. GLN(98)/3. GLN(100)/3.

Cag) fin GLN(S55)/1. GLN(S7)/1. GLN(S9)/1.

Va Wolbia

t=)

le

Yh ts of

ied id od

ca

cn cn po ce a- ha om

hI ~J = cO

92739047 287B8E+02/, OLOO9S163993E+02/, 0934258675 25E+02/, L77Sh9S41140E+02/, 2619054983235 7E+02/, BATOL7234428E+02/, AS26884 900S0E+92/ : S18904022097E+02/, SOS649F09719E+027, 9291097 651L0E+02/, 7B0675734404E+02/, B48951552954E+02/, 9576660135508E+02/, O46868SS74657E+02/, 136528299499E+02/, 26534991 267E+02/, 5171788719695 +02 bf 40B8150588708E+02/, 499541180408E+02/ 5913542053695E+02/

3 FOR GAMMA(X)

(Sheet 24 of 25)

GLN66 GLN67 GLN468 GLN69 GLN7O GLN7i GLN72 GLN73 GLN74 GLN75 GLN76 6LN77 65LN78 GLN79 GLN8O GBLN81 GLN82 GLN83S GLNB4 GLNS85 GLN86 GLN87 GLN8B8 GLN89 GLNIO GLNI1 GLN92 GLN93 GLN94 GLNIS GLN9S GLN97 GLN9S GLN99 GLNLOG GLNIO1 GLNI62 GLNIOS GLNIO4 GLN1IOS GLN1OS GLN107 GLN168 GLN109 GLN110 GLNi1i GLNI12 GLNi13 BLNI14 GLN115 GLN116 SLN117 GLN118 GLN119 GLNi20

on]

lao)

oO

40

af

XX=2.4T

DO 24 1=1,NkR DGAMN=DGAN*X Xx XX=KX#1. GANLN=ALOGS(DGAN) RETURN

X.GT.xXLIM1

RX=1./%

RXX=RX#RX

IF((X-XLIM2).LT.0.) GO TO 41 PX=O(1) #RXX+O(2) GAMLN=PX#RX+(X-.5) #ALOG(X) -K+RTWPIL RETURN

X.LT.XLIM2

PX=P (1) SUM=(X-.5) #ALOG (X) -X DO 42 1=2,5 PX=PX#RXX+P (1) CONTINUE GAMLN=PX#RX+SUM+RTWEIL RETURN

TABLE LOOK UP FOR INTEGER ARGUMENTS LESS THAN OR EQUAL 100.

IF(NDX.GT.100) GO TO 40 GAMLN=GLN(NDX)

RETURN

END

Figure A3. (Sheet 25 of 25)

A29

GLNL21 GLNL22 GLNi23 GLNL24 GLNI2S GLNI2 GLNi27 GLN128 GLNi29 GLN1 20 BLN13i GLN1S2 GLNI33 GSLN134 GLNISS GLH136 GLN137 GLN138 GLNI29 GLN140 SLN141 GLNI42 GLN143 GLNi44 GLNI45 GLN146 GLN147 GLN148 GLNi49 GLNLSO GLNIS1 GLNIS2

deh oa ; oe wig ih aa rT: e 2a pare me

a ant iaa ta ee

O50, Kal i ; racy i", a O: Lh al ft +“ ae

oar cfhis a ¥

tet “Ate 4 EF mney nnbhiok unite o) tants

2 oe eigen

eee a ua en ee 16 7? AR

| Acre oe

28 - P| a 5 : :

a us : Shih) shyb3' #e nay’ np avemniin wpe dot Ue 08 sate es . - 5 ,

; = Ww y ae 4 a mite es reLeolorie Le ee fut 7. a3. aaa, cana etd =

PRONE I Pe

a | teen ae es ae ee _ Risto sgtetienn . ae 4) ub j ; 7 ; : aie : a : - Fr : Fonsio| ' Ee

7 - grits”

7 et et we - wig ~ cas LF acwie HAS Ty.

- Fie ©

: ity |

ie BHR ie

a) we

ire, i a io >.

APPENDIX B: NOTATION

Incident wave amplitude Softy

Gravitational acceleration

Water depth

v=1

Bessel function of the first kind

Wave number

Non-negative integers

Polar coordinate

Temporal coordinate

Horizontal coordinate

Horizontal coordinate

Bessel function of the second kind Vertical coordinate

Incident wave angle

Free surface displacement

Polar coordinate

Angle related to wedge angle

Value related to wedge angle

Velocity potential function

Horizontal component of the velocity potential function Defined in Equation 12

Incident wave velocity potential function Reflected wave velocity potential function

Scattered wave velocity potential function

Finite cosine transform of ¢

Wave radian frequency

Bl