hee ah se LS, Ay Coat Ge Fe Ba af es Les. Army Eo: C F— Coastal Engineering Research Center DEVELOPMENT OF A METHOD FOR NUMERICAL CALCULATION OF WAVE REFRACTION TECHNICAL MEMORANDUM NO.6 DEPARTMENT OF THE ARMY CORPS OF ENGINEERS wii chbbLeoo MIN TOEO oO AMAT IOHM/Taw EMCO DEVELOPMENT OF A METHOD FOR NUMERICAL CALCULATION OF WAVE REFRACTION by W. Harrison and W. S. Wilson TECHNICAL MEMORANDUM NO.6 October 1964 Material contained herein is public property and not subject to copyright. Reprint or re-publication of any of this material shall give appropriate credit to U.S.Army Coastal Engineering Research Center HIS PUBLICATION WITHIN THE UNITED STATES IS MADE BY THE U, S.ARMY LIMITED FREE DISTRIBUTION OF T N. W., WASHINGTON D.C. 20016 COASTAL ENGINEERING RESEARCH CENTER 5201 LITTLE FALLS ROAD, FOREWORD A knowledge of the energy influx along a shoreline is of consider- able importance to an understanding of many beach processes and long- term trends in beach responses. This study represents a step toward the development of a rapid method for routine determinations of wave power along a shoreline, using observed or hindcast deep-water wave character- istics and high-speed computer programs for the calculation of wave refraction. Specifically treated here is a computer program and pro- cedure for the construction of wave rays. An example of the method is presented in which wave rays are brought from deep water into the Atlantic shoreline of the City of Virginia Beach, Virginia. Detailed explanations of the computer programs used, instructions for their operation, and sample inputs and outputs are given in appendices. A series of sugges- tions is also given for improvement of the method presented. This report was prepared by Dr. Wyman Harrison, Associate Marine Scientist, Virginia Institute of Marine Science, in pursuance of Contract DA-49-055-—CIV-ENG-64-5 with the Coastal Engineering Research Center, and in collaboration with Mr. W.Stanley Wilson, formerly a graduate student at the Institute. The study was supported by the Coastal Engineering Research Center (formerly the Beach Erosion Board of the Corps of Engineers). Computing Centers at Northwestern University, the College of William and Mary, and NASA (Langley Field, Virginia) extended every cooperation. Lieutenant G. Griswold, Oregon State University, first suggested to the authors the use of numerical methods for the calculation of wave refraction. The computer program used in this report for calculating wave refraction was extensively modified from a program under development in 1963 by Lieu- tenant G. Griswold and Mr. F. Nagle of the U. S. Navy Weather Research Facility, Norfolk, Virginia, and Mr. E. Mehr of New York University. Mr. J. Gaskin of IBM and LCDR C. Barteau of the U. S. Navy Weather Research Facility aided in adapting the Griswold-Mehr program for use on the 7094 computer, while Mr. J. Curran and Mr. R. Libutti of IBM, aided in attempts to adapt the program to the 1620 computer. Both programs were considerably revised by Wilson prior to their presentation in the appendices of this study. The authors are indebted to Professor W. C. Krumbein of Northwestern University, Consultant to the Coastal Engineering Research Center, and to Mrs. Betty Benson of the Northwestern University Computing Center, for making available the surface-fitting program that is used in the PRMAT subprogram and SURFCE subroutine of the wave-refraction program. Pro- fessor B. R. Cato of the Mathematical Department of the College of William and Mary, provided helpful advice at various stages of the study. This report is published under authority of Public Law 166, 79th Congress, approved July 31, 1945, as supplemented by Public Law 172, 88th Congress, approved 7 November 1963. CONTENTS Page No.. ANDSTRAGIE Sm SAN 25 ene ain ee eet RON cece ee eed Be ere ai ees cen Be INTRODUCTION ea Sy eae eg SN Oe et ee eas el ee Ve ee ae oe ACI UiNG CF PrOjece = Shae ao So Se) a Sa Je) SS 6 Earlier Studies of Wave Refraction - - - - - - = = = = Scone Of Presemt Study =o SoS Ss 6s 5S Ss Ss 4 952 S 6 & NETO 1 Se cy Sh Sa) ea eel Gee es) dy I eee PrevtousmWomk y= =I Hi = y= =) a ae a Construction of Wave Rays - - - - - = -=- = = = = = = Computer Operations = = = = = = = = = = = = = = -= Representation of Output - = - - = - = = = = = = = PROGRAM DEVIRLORMENT So oS os Ss Ss oso Ss Ss Sos Ss s&s Ss ss oo Imeeroolation SUREACCS S 2 =o S Ss Ss Ss 6 Ga Ss is Ss « 9 Ray CicVvAtUSe ApororaintleOm = = |e Sse oes Ss Ss = 56 6 & 12 Geidéd Comeiderationg © Ss = Ss eis 6 S&S s'S4 6 & SBS =& Wave Ray Representation PN NSE SS) oA GE i AEE UAE SC ices Ba a a TROTIMS OF Ray CONtTRUGEIONS =& = SS Ss & Ss S95 Ss = =! 19 CONCLUSION: (6. “e 6. Stee He aR ee ee Es! ey eS: GE iS aS RERINENCGES S68 co Se SoS aoe SoS Ss os So 6 oa SS S&S Soe Ss 6G 14 0 wmaoOKRAHR WWNHD F 0 1 hh W dH | pa Ww TABLE 1 - Computer Input Specifications for Waves of 4-Sec. Period 16 TABLE 2 - Computer Input Specifications for Waves of 6-Sec. Period 7 LIST OF FIGURES Be Se a a eae ie ei sai 18 CURE Seine Ole Py eee yo eee Cay ae hy Totes Sco By ih f eNO lac O APPENDICES A - Computer program for compilation of wave velocities as a function of wave period and water depth Sats Se ee 30 B - Computer program for distribution of Wave Velocity values over a grid of depth values as a function of wave period, - 34 C - Computer program to produce matrices for deriving equations of linear surfaces - - - - = = = = = = = = = = 38 D - Computer program for wave refraction (linear-interpolation) using the IBM 1620 - - - - = - = = = = = = = = 40 E - Computer program for wave refraction (linear-interpolation) sHineerehe! LB Ma OOF = Mae MeN Te demu wea Se) cia Se Ma SS Pie | 5's ; : ; OS 5 ad OG ios od F - Derivations relating — with — and —= with —- - - - - 59 eax aK 3Y SY G - Method used in obtaining most frequent combinations of deep- water height, period, and direction of waves capable of striking Virginia Beach Sey has (epee eek et (ies eri ior i 60 er ) 4 Fee ie: Lane ee SSP EIP! erred ‘ . tt ei a Ti aS Mah Be BRS HS e wiley 2 ¢ pat es dss i asi: vee ily r 1 a rey i 4 . i DEVELOPMENT OF A METHOD FOR NUMERICAL CALCULATION OF WAVE REFRACTION by W. Harrison U. S. Coast and Geodetic Survey, Washington, D. Ca and Virginia Institute of Marine Science and W. S. Wilson The Johns Hopkins University, Baltimore, Maryland? ABSTRACT Steps in wave-ray construction are as follows: 1. Select wave periods and approach angles for each series of rays to be constructed. 2. Prepare a grid of depth values for the area of investigation. 3. Use a computer program to obtain a table of water depths and related wave velocities for each wave. period Sellefetedmin Ge) 4, Make a grid of wave-velocity values for each wave period selected in (1), for the area of investigation, using a second computer program which takes as input the depth grid of (2) and the appropriate depth-velocity table of (3). 5. Derive, using a third computer program, matrices for use in the linear interpolation scheme of the computer program of (6). 6. Calculate, for each wave period specified in (1), the points along a wave ray using a computer program which takes for input: a. The appropriate velocity grid of (4). b. The matrices of (5). c. The origin points and approach angles of (1) for given wave rays. lpresent address. Study completed while at the Virginia Institute of Marine Science. A linear-interpolation scheme (using the least-squares method) is used in determination of wave velocity at a given point along a ray. Ray curvature is then calculated at this point and an iteration procedure is solved to obtain the position of the next point. The ray terminates at the shore or grid border. An example of usage of the programs is presented in which 52 rays for waves of 4- and 6-second period, from six different deep-water origin points, are brought into the Atlantic shore- line of the City of Virginia Beach, Virginia. The procedure outlined is in the developmental stage, and suggestions for improvements are given that should offer a quick, accurate, and objective method of constructing wave rays. INTRODUCTION Background of Project An ultimate understanding of the changes in topography of a given shoreline will be derived from a considerably fuller knowledge of the input and transformation of energy along the strand than is currently available. Because by far the greatest amount of energy expended in the beach-ocean-atmosphere system is associated with breaking waves in the surf zone, it becomes of the greatest importance to evaluate the long- term distribution of wave energy along the shoreline. This energy dis- tribution depends mainly upon the energy of the original deep-water waves, as modified by refraction when the waves move shoreward through shallow water. It was with these general considerations in mind that the authors embarked upon a study of wave refraction at Virginia Beach, Virginia, where previous studies (cf. U. S. Congress, 1953; Harrison and Wagner, 1964) suggested the desirability of determining wave power distributions. The authors were armed with an “operational"' computer program said to be capable of rapid and accurate calculation of wave-ray paths. As in the case of many "operational" computer programs (cf. the cogent comments of Adams, 1964), the authors soon discovered that they were either in diffi- culty simply trying to make the program operate, or that they were not in agreement with certain of the logical steps involved. It soon became apparent that the entire method needed to be reviewed and revised. The emphasis here on development of a numerical method for calcu- lation of wave refraction does not mean that the long-range goal of assessing wave-power distributions along shorelines has beén set aside. The authors feel, moreover, that the high-speed method ultimately adopted for calculating wave refraction is of such basic importance to the larger problem that a detailed presentation is warranted at this time. Earlier Studies of Wave Refraction The investigations of Krumbein (1944) and Munk and Traylor (1947) were among the first in which steps were taken to assess the link between wave refraction, the energy distribution along the shore, and beach erosion or deposition. These workers constructed wave-refraction diagrams for selected deep-water wave periods by hand-drawn methods, a time-consuming process that at best permits only partial assessment of the effect of a spectrum of waves on shore erosion. About ten years ago, Pierson, Neumann, and James (1953) considered wave refraction effects in some detail and concluded (1953, p. 186) that to use only significant height and the average "period" of the waves in deep water and to refract the waves with these two numbers will lead to totally unrealistic results. These authors went on to outline a method (1953, p. 197) for forecasting wave heights and characteristics at a point in shallow water near a coast. The method involves construction of a ‘large number of ray diagrams by graphical methods. Except for the great amount of work involved, the method holds promise for a realistic assessment of energy distributions along a shore- line, assuming the deep-water spectra can be adequately approximated. Although more "practical" methods for determining wave refraction (cf. Silvester, 1963) and for assessing wave energy along coasts (cf. Walsh, Reid, and Bader, 1962) have been proposed, the present authors feel that the adequate treatment of the problem of energy inflow must of necessity include large-scale computations. The treatment of wave refraction, for example, should be as thorough as that outlined by Pierson, Neumann, and James (1952, p. 197), and this requires many man hours. The advent of high-speed computers offers the possibility of significant reductions in the time required for construction of refraction diagrams, wave spectral forecasts, and, resultant distributions of energy at the shore. The machine methods given below for constructing refraction diagrams represent a step in the attempt to develop a comprehensive computer program for the routine determination of wave energy along shorelines. Scope of Present Study It will be assumed that the deep-water wave spectral periods, heights, and directions are known, and that they may serve as input to a computer program for determining wave refraction as part of an overall method such as proposed by Pierson, Neumann, and James (1953). Examples of the machine refraction of six separate combinations of wave period and direction off Virginia Beach, Virginia, are given to illustrate the method of refracting wave rays. The assumptions and methodology used in constructing the wave rays are fully treated, and suggestions for im- proving and amplifying the program are then presented. METHOD Previous Work The presentiy acceptable method for constructing wave-ray diagrams involves the hand method of graphical construction first proposed by Arthur, Munk, and Isaacs (1952) and later repeated in Pierson, Neumann, and James (1953) and U. S. Army, Corps of Engineers (1961) Technical Report No. 4, of the Beach Erosion Board. This procedure for graphical construction of wave rays involves the use of Snell's Law, sin Oy sin oP) Cy where for two successive points ona ray, 1 and 2, C = wave velocity, and a = angle between the wave crest and the bottom contours. In the hand- refracting methods, Snell's Law requires the assumption that between points 1 and 2 there exist straight and parallel contours. The computer program for-construction of wave rays (that is to be described later) involves the fitting of interpolative surfaces to points of a wave-velocity grid. This program, although requiring certain assumptions of its own, avoids the assumption of straight and parallel contours between successive points on a ray. With this program an expression for ray curvature is solved quite independently of the use of Snell's Law. The authors' introduction to the possibilities of computer programs for the numerical construction of wave rays was gained originally through discussions with Lieutenant Gale M. Griswold, then of the U. S. Navy Weather Research Facility, Norfolk, Virginia. Ideas on the subject had appeared in mimeographed reports (Mehr, 1961, 1962a, 1962b; Griswold and Nagle, 1962) and a published paper (Griswold, LOGS). A computer program, although not actually operational, was provided by Griswold for develop- ment in this study. (This program will be referred to henceforth as the Griswold-Mehr program. ) Construction of Wave Rays Selection of Input. The first step in wave-ray determination begins with the construction of a grid of depth values for the area of interest. Geographical limits for a grid of depth values take into account the probable origin points for waves of all deep-water (d/Lo>>0.5) approach angles of interest. By convention, the grid origin is located in the southwestern corner of the area to be covered, with the X axis extending eastward along the southern border and the Y axis extending northward along the western border. The grid interval is selected such that, in a given cell, bottom contours can be represented by straight and parallel lines. (This representation will be discussed in detail later.) Actual or interpolated depths at grid intersections are recorded to the nearest foot, and all real depths are made positive. Contours are then drawn of depth values extending several grid units out from shore. At this point, contours that are symmetrical reflections of nearshore bottom contours are drawn on the land, over a perpendicular distance from the shore of two grid units. "Depths" at the grid points on land are then derived in the region of the symmetry contours on land; these "depth" values are made negative. All other land "depths" (that is, those farther than two grid units from the shore) are made zero. (The procedure of assigning negative values for nearshore "depths" on land is required for the fitting of wave-velocity surfaces.) The depth values so obtained are prepared as input to the DISTV program (Appendix B). Example of Input. An area of the mid-Atlantic Bight (Figure 1) is chosen to illustrate the method of data input and the computation of wave refraction. A specific wave-ray target area if the region covered by the depth grid is the shoreline of the Borough of Virginia Beach in the City of Virginia Beach, Virginia (Figure 1). The wave input data used here are not representative of any given wave spectrum observed or hindcast. They approximate the 15 most commonly observed combinations of wave height, period, and direction observed at the Chesapeake Lightship (Figure 1). These combinations were determined and condensed for input to the computer at the request of the Coastal Engineering Research Center. The method used in determination of the combinations is described in Appendix G. The authors have used the six wave period -- direction combinations (Tables 1 and 2), extracted from the 15 height-period-direction combinations of Appendix G, merely to illustrate the method of ray construction. The depth grid (Figure 1) of the example was chosen so that the origin points of six second waves approaching Virginia Beach from 60° and 90° T would be positioned in water depths of greater than 92 feet (d/Lo>0.5 for T = 6 sec.). Outlines of the 99 X 81 unit depth grid used in the example are shown on Figure 1; the grid origin is located at 7691.9'W, 36°39.5'N, and the grid interval is 3,040 feet. U. S. Coast and Geodetic Survey (boat- sheet) charts 5988, 5990, 5991, 5992, 5993, and 6595 were used in obtaining depth values. Where these charts did not offer coverage, depths were picked off charts 1222 and 1227. Smoothed contours of the depth values at grid points appear on Figure 2. Deep-water starting points and directions for the 52 wave rays of the example are given in Tables 1 and 2, and are shown plotted on Figure 3. An example of the coding of thewave-ray and wave-velocity input data is given under the heading "INPUT" in Appendix D. Computer Operations Programs Used. For each wave period selected for the investigation, a table of depths and their associated wave velocities is prepared. This is accomplished by solving the following equation (U. S. Army, Corps of Engineers, 1942; U. S. Navy, Hydrographic Office, 1944; Mason, 1950) by an iteration process: _ 2d C= TT tanh Tic) where C = wave velocity, g = acceleration due to gravity, T = wave period, and d = water depth. This equation has been programmed, and a sample depth-velocity table has been prepared for T = 4 seconds. (See COMPV in Appendix A.) As in other wave refraction studies (Pocinki, 1950; Pierson, 1951; Pierson, Neumann, and James, 1953), it is assumed that wave velocity is a function only of water depth and wave period, as expressed by the above equation. Various factors such as bottom friction, percolation, reflection, currents, and winds are considered as having no effect on the refracting waves. Given the absolute value of a water depth, it is possible to check the appropriate table, which has been previously prepared, for the associated wave velocity. Preparation of an entire velocity grid for each wave period to be studied is then carried out. This procedure has also been programmed, and a portion (from X =.0 to X = 19, and from Y = 0 to Y = 2) of the depth grid described above has been derived for T = 4 seconds. (See DISTV in Appendix B.) The procedure for determining the wave-velocity value at an arbitrary point in a grid cell involves fitting a plane to the velocity values at the four grid intersections that bound the grid cell in question. This linear surface is fit by the least-squares method, using an equation of the form: Cr Et BOX eB, where C = wave velocity, E's = coefficients, and X,Y are the grid coor- dinates for the arbitrary point. With the least-squares method of surface fitting, it is possible to obtain certain matrices which are used each time the equation of a plane is derived for a grid cell. These matrices have been derived in PRMAT. (See Appendix C.) That portion of the surface-fitting procedure which must be carried out each time a linear equation is to be derived is given in SURFCE subroutine of MAIN 1620 (Appendix D) and MAIN 7094 (Appendix E). Determination of each desired ray for a given period is accomplished by first specifying origin coordinates and an angle of approach. The actual ray is constructed by plotting and connecting the series of successive calculation points (computed by MAIN 1620, Appendix D, or MAIN 7094, Appendix E) which range across the velocity grid until striking the beach or grid margin. (As discussed by Griswold (1963, p. 1722) the rays may be run from the shore seaward, if the construction of a refraction diagram at a point is desired.) Velocity at each point is calculated as mentioned above; ray curvature (FK) is calculated by using the following expression (Munk and Arthur, 1952): FK = = | sin A (&) - cos A () | where A = approach angle, as defined in Appendix D. In order to determine Xpn+1, Yn+i4, An+1, and FKyn+ 1 for calculation point n+1 (with those similar values known for point n), the following equations (Griswold and Nagle, 1962; Griswold, 1963) are solved by an iteration procedure: AA = (FE + FK 40/2 A =A + AA n+1 n A = (AL + Aled /2 Xie =X + DcosA vee = Y +DsinA where D equals the incremented distance between points n and n+l. (See MOVE subroutine of MAIN 1620, Appendix D, and of MAIN 7094, Appendix E.) The computers to be used for the programs cited in the appendices are the IBM 1620 (with floating-point hardware and 60K memory) and the IBM 7094. FORTRAN II is the language used for the 1620 programs; FORTRAN Iv is used for the 7094 programs. The computer and language for a given program are specified on the first comment card of each program print-out. The number of digits carried by a computer during internal computations for floating-point numbers (designated by f) and the number of digits carried for fixed-point numbers (designated by k) vary with the different programs. For the 1620 programs, f is given in Columns 2 and 3 and k is given in Columns 4 and 5 of the first card of each program print-out. For the 7094 program, f = 8 and k = 5. Example of Sample Output. Sample input and output listings are given in Appendix D for three wave rays. These three rays, numbered 1, 2, and 3 in the OUTPUT listings, correspond to rays numbered 33, 32, and 31, respec-— tively, of Table 1 and Figures 2 and 6. Representation of Output Possible Methods. Output from the computer programs could be presented in a variety of ways, depending upon the requirements of the investigator. Precision work on caustics, for example, might entail skillful plotting of the values at MAX points with precision drafting techniques. If a rapid analysis of ray paths from deep water to the shore is all that is required, however, the investigator might consider an X-Y plotter attachment for rapid printing of the computer output. For many shore engineering studies, only the final few hundred yards of the wave rays will be of practical value, while for studies of energy or wave-power distribution along a fixed segment of the shore, it is possible that only the terminal points of the rays at some specified distance from shore would be of value for presentation. Example. The results of the computer refraction of the 52 wave rays of Figure 2 are presented in Figures 3-8, which show only the last portions of the wave rays relatively close to shore. Successive calculation points were plotted on a grid and then connected by a smooth curve. The rays themselves vary from the almost unmodified ones’for 4-second waves that approach the beach relatively head on (Figure 5), to the 6-second ones that actually cross (Figure 8) within the limits of the figure. Attention is drawn to the terminal points of the wave rays; some of these points (such as those of rays Nos. 14, 37, and 52 as shown in Figures 4, 7, and 8, respectively) are closer to the shoreline than others. These variations are due to the fact that a ray terminates because either the next calculation point is in an area of zero or negative velocity, or curvature approximations are not converging. Thus, with a constant in- cremented distance (D), all rays do not reach a similar distance from the shoreline. Ray No. 24 (Figure 5), on the other hand, makes a pronounced curve near its terminus that appears out of context when compared to the adjacent wave rays. This is due to a combination of (1) a poor "fit" by the linear-interpolation surface at this point, and (2) poor grid control (i.e., poor representation of depth values recorded on the initial depth grid). Suggestions for elimination of these factors are discussed in the section "Grid Considerations." PROGRAM DEVELOPMENT The computer programs presented in the appendices have by no means been refined to the fullest extent. At the moment the following factors may be considered of most importance to their continued development: (1) improvement of the surface-fitting scheme for wave velocity interpolation, (2) improvement of the method of ray-curvature approximation, (3) addition of a provision for changes in grid scale and incremented distance as a ray moves shoreward, (4) refinement of ray-plotting techniques, and (5) testing of the revised program at a suitable place in nature. Interpolation Surfaces "Forced-cubic" Interpolation Surface. First mentioned by way of review is the surface-fitting scheme for interpolation of wave velocity used in the Griswold-Mehr program. Their scheme involves fitting a cubic surface that (1) passes exactly through the velocity values at the grid intersections of the given grid cell, and (2) is the cubic surface of best fit (by the least-squares method) to the velocity values at the eight grid intersections closest to and surrounding the four grid inter- sections of the given cell. This cubic surface is called a "forced-—cubic" surface because it is "forced" to pass through the four innermost velocity values. Because it permits the taking of first and second derivatives for use in wave intensity calculations, as explained by Griswold (1963, p. 1720), it was the first surface-fitting program used in this study. An example of the results of its use appears in Figure 9C. It is obvious from the figure that this method of interpolation is invalid. The forced feature of the surface creates undesired ridges and/or troughs in the velocity surface. Thus, depending upon the location of a given calcu- lation point in a grid cell, the ray can be deflected erratically from a "normal" path. Results such as those exemplified in Figure 9C necessitated a search for an interpolation surface that could adequately portray the general trend of wave velocity change in a given cell. Quadratic-interpolation Surface. The Griswold-Mehr program was altered by insertion of a quadratic surface of best fit (by the least- Squares method) subroutine for preparation of the interpolation surface. In order to derive a quadratic surface, at least six data points are necessary. It was decided to use the velocity values at the closest 12 grid intersections. The results (Figure 9B), however, did not yield a satisfactory interpolation surface; this was evident when calculation points for a ray had moved within two grid units of the shore. This is because the negative "land" velocity values, used in derivation of the surface, made the general tilt of the surface steeper than indicated by the velocity values at the central four grid intersections. This is shown in Figure 9B where the rays bend, after moving within two grid units of the shore, more than do the rays produced by the linear-interpola- tion program (Figure 9A) or the rays produced by the graphical-construction method (Figure 9D). Linear-interpolation Surface. The linear-interpolation programs (MAIN 1620, Appendix D; and MAIN 7094, Appendix E) were developed in order to remedy the excessive tilt created by the quadratic-interpolation program. The advantage offered by the linear-interpolation program is that only the velocity values at the central four grid intersections are needed in order to derive the surface. The assumption, presented under the heading "Selec- tion of Input," stating that the grid interval be selected such that bottom contours in any grid cell be represented by straight and parallel lines, is not actually valid. Although bottom topography may be represented by a plane (i.e., it changes linearly) in a given cell, the associated velocity values may be changing exponentially (i.e., velocity is an exponential function of depth) in the same cell. However, for purposes of this program, it is assumed that the velocity values in a given grid cell can be adequately represented by a plane. The use of the variable PCTDIF in MAIN 1620 (card number RAYN 19, Appendix D) and MAIN 7094 (card number RAYN 20, Appendix E) serves to give an indication of the relative "fit" of a surface to a given set of velocity values. The output values for PCTDIF (see OUTPUT, Appendix D) give the percent difference between only one of the four velocity values at the nearest four grid intersections and its related value on the plane fit to the same four values. PCTDIF by no means represents the degree of "fit" of the surface to the four values because these percent differences may vary for each of the four values to which a plane is fit. This is especially true in cells where there are great differences between the four velocity values (such as near the shore). PCIDIF does, however, give an estimate of the relative error encountered in the interpolation in a given grid cell. Just as the quadratic-interpolation program yielded an excessive tilt when the given grid cell in use fell within two grid units of the shore, the linear-interpolation program yields an excessive tilt when the given cell falls within one grid unit of the shore. Therefore, it is expected that PCTDIF will be large in value when the grid cell being interpolated is located near the shore. In view of this fact, the rays run with the linear surface-fitting program should be considered as rough approximations only, when closer than one grid unit to the shore. Interpolation Surface Versus Graphical Method of Ray Construction. Consideration of the above-mentioned surface-fitting programs led to the choice of the linear surface of best fit as the most suitable one for use in construction of wave rays. In general, it was found that rays run with the linear surface of best fit compared most favorably with rays constructed by graphical methods, and this is the case in the plotting example (Figure 9 A and D). An advantage of the least-squares linear-interpolation method over the graphical method lies in the fact that ray curvature can be com- puted at a number of points between contours. Aside from the absolute validity of the two methods of wave-ray construction, the computer method is estimated to be many times faster than the hand method when only the terminal points are desired of a large number of rays. Where only a few rays are desired, the hand method is clearly the most rapid and economical. Because the practice of refracting only a few rays for a few wave periods yields totally inadequate information upon which to assess wave heights and energies at the shore, it seems clear that the real considerations are not so much the man-hours involved in depth-grid and program-deck preparation for the computer as they are in the necessity for the more realistic results that can be afforded by the computer construction of wave rays for entire wave spectra. ? Future Interpolation Schemes. A more valid interpolation scheme would be obtained if a grid of depth values were input to MAIN 1620 or MAIN 7094 instead of a grid of velocity values. The assumption, that the grid interval be selected such that depth contours in any grid cell be represented by straight and parallel contours, would then be perfectly valid. In this case, a depth value would be obtained from the interpola- tion surface at a given calculation point. Then, using the procedure used by COMPV (Appendix A), the depth value could be converted into a wave- velocity value. It is noted that, in order to obtain curvature (FK) at a calculation point, && and 35 are needed. If a linear-interpolation program Ge used, with an input depth grid, expressions would be available for x and a: on oor ou tne. ne Lact onsntp (derived in Appendix F) could then e used to obtain 2+ andl: x oy oc = od ZA le = od Z 0x Ox TOV Ove | whe re BM sotphen lym 0 adherens nani! Ck" Ck" " Se a St ut -1 1-k"C TSENCMASIEIC GH We cue : In this expression k' = T/4T and k" = 21/gT. If a quadratic surface were to be used for interpolation similar relations could be derived to relate 2 2 2 2 ue and ome with oad and ole , respectively. ae a" ae aye A reason for using a quadratic (or a higher-order polynomial) surface is that the second partial derivatives can be obtained. This fact was used in the Griswold-Mehr program to obtain values for Beta (Munk and Arthur, 1952), the coefficient of refraction, and H/Ho. (As mentioned, however, their interpolation scheme rendered their results invalid.) It is apparent that certain problems with the quadratic (ar higher-order polynomial) interpolation schemes must be resolved before such sophisticated parameters as Beta or H/Ho can be estimated. It should be possible to derive a quadratic surface by heavily weighting the velocity values at the central four grid intersections while still using the surrounding eight velocity values. This procedure may produce a good interpolation surface; if so, it would be quite worth while to calculate estimates of the above parameters. Ray Curvature Approximation An entire grid of velocity values is not a smooth and continuous surface when observed as a set of planes in which each grid cell is represented by a specific plane of best fit to its four bounding velocity values. With this in mind, it is not surprising that all curvature approxi- mations (see MOVE subroutine of MAIN 1620, Appendix D; or MOVE subroutine of MAIN 7094, Appendix E) fail to converge to a single value when determining a new point along a ray path. This is especially apparent when adjacent planes present a large discontinuity in wave-velocity values at their common edges. This fact is the reason for the variable MIT included in MOVE subroutine. In case the curvature approximations oscillate among three or more values after 20 iterations, the ray is terminated, as no valid curvature approximation can be made. Although infrequent, sometimes (not shown in OUTPUT, Appendix D) the curvature approximations oscillate between two values. In this case, the message "CURVATURE APPROXIMATED" is included with the output in order that the operator note that the curvature used for the given point was an average of two values. The Griswold-Mehr program did not have such a check on the curvature approximations; the average of the last two approximations after 10 iterations was used as the new curvature value, if the values did not converge. This is another reason for the erratic ray behavior near the shore in Figure 9C. Grid Considerations For reasons outlined in a previous section ("Future Interpolation Schemes") assume that a grid of depth values will be input and used for derivation of the interpolation surfaces. It will be necessary in the future program to transfer a ray to a larger scale grid when approaching the shore in order to allow better grid control (i.e., better representation of depth values) in an area where the velocity values are rapidly changing. There is, however, a limit to the maximal scale of a grid as beach and nearshore topography are constantly changing, especially during each storm. Thus the selection of grid interval is a question to be pursued in additional studies. It will also be necessary to vary the distance (D) incremented between successive calculation points along a ray in order that the ray may proach as close to the shoreline as possible. As an example, the value of D could be assigned such that D = 0.5 when d/Lg>0.5 and D = d/Ly when d/L .<0.5. Wave Ray Representation At present, output must be plotted by hand; and, if there are a large number of rays to be constructed, this can be a tedious and time- consuming process. The use of an X-Y plotter, if adaptable to changing grid scales, would be an ideal scheme for rapid plotting of wave rays. Another scheme might involve the hand or machine plotting of the (numbered) terminal points of the wave rays, along the shore or grid margin. Testing of Ray Constructions The absolute validity of the ray constructions (Figures 3-8) is impossible to evaluate without a test in nature. It would be desirable to test the ray constructions in an area where significant refraction of relatively constant-period wave trains occur and where the bathymetry is well known. The use of aerial photography would be an invaluable aid in conducting such a test. It is possible to test ray constructions by comparison with analytic solutions for wave rays passing over algebraically-described surfaces (Pocinki, 1950). Griswold (1953), in testing the Griswold-Mehr program, used both a straight uniformly-sloping beach and a parabolic bay and found little variation between computed and analytic rays. However, such a theoretical test is not sufficient reason for acceptance of a computer program utilizing an interpolation scheme. In these theoretical tests, an algebraically-described surface of velocity values is input to the computer. The close agreement then noted between computed and analytic rays is due to the fact that the given interpolation scheme can easily and accurately represent the velocity surface when interpolating within any grid cell. In order to adequately test such a computer program, an irregular surface of velocity values must be provided. This necessity supports the need for a test in nature, as mentioned in the previous paragraph. CONCLUSION The procedure described in this report, when fully improved using the accompanying suggestions, will constitute a rapid and accurate method of wave ray construction. At the present stage of development, however, the procedure must be accepted with reservation. REFERENCES Adams, B. B., 1964, Program pool debated: Geotimes, v. 9, p.4. Arthur, R. S.; Munk, W. H., and Isaacs, J. D., 1952, The direct construction of wave rays: Trans. Amer. Geophys. Union, v. 33, no. Griswold, G. M., and Nagle, F. W., 1962, Wave refraction by numerical methods: Mimeo. Rept., U. S. Navy Weather Research Facility, 19 p. Griswold, G. M., 1963, Numerical calculation of wave refraction: Jour. Geophys. ROS 5 Wo ES. wo IVISAL( 23 Harrison, W., and Wagner, K., 1964, Beach changes at Virginia Beach, Virginia: U. S. Army, Corps of Eng., Coastal Eng. Res. Center Tech. Memo. (in press). Johnson, J. W.; O'Brein, M. P.; and Isaacs, J. D., 1948, Graphical construc- tion of wave refraction diagrams: U. S. Navy Hydro. Off. Pub. No. 605, SS 50) Kaplan, Wilfred, 1952, Advanced calculus: Addison-Wesley Pub. Co., Reading, Mass., 679 p. Mason, M. A., 1950, The transformation of waves in shallow water: in Proc. First Conf. Coastal Engr., Counc. Wave Res., Berkeley, Calif., p. 22-32. Mehr, E., 1961, Surf refraction: Prog. Rept. No. 1 to U. S. Navy Weather Res. Fac., Norfolk (Contract No. N189-188-5411A, mimeo. rept.), 4 p. , 1962a, Surf refraction: Prog. Rept. No. 2 (see citation immediately above), wl p. 9 1962b, Surf refractions; computer refraction program: Final Rept. of N. Y. Univ. Statistical Lab. (Coll. Eng., Res. Div.) to U. S. Navy Weather Res. Fac., Norfolk, on Contract No. N189-188-5411A, 21 p. Munk, W. H., and Traylor, M. A., 1947, Refraction of ocean waves: a process linking underwater topography to beach erosion: Jour. Geology, v- 5), Wa Ie BSs Munk, W. H., and Arthur, R. S., 1952, Wave intensity along a refracted ray: in Gravity waves: U. S. Dept. Comm., Nat. Bur. Standards Circ. 521, pe 95-108. Pierson, W. J., Jr., 1951, The interpretation of crossed orthogonals in wave refraction phenomena: U. S. Army, Corps of Engineers, Beach Erosion Board, Tech. Memo. 21, 83 p. Pierson, W. J., Jr.; Neumann, G.; and James, R. W., 1953, Observering and forecasting ocean waves: U. S. Navy Hydro. Off. Pub. No. 603, 284 p. Pocinki, L. S., 1950, The application of conformal transformations to ocean wave refraction problems: Trans. Amer. Geophys. Union, v. 31, p. 856- 866. Silvester, R., 1963, Design waves for littoral models: Jour. Waterways and Hanbor (Dav .,5 Amer. SOc. Civ. Ene. vj. 69, WNS. p. Sr—4i U. S. Army, Corps of Engineers, 1942, A summary of the theory of oscillatory waves: Beach Erosion Board Tech. Rept. No. 2, 43 p. , 1961, Shore protection, planning and design: Beach Erosion Board Tech. Rept. No. 4, 21 p. U. S. Congress, 1953, Virginia Beach, Va., beach erosion control study: 83rd Congress, 1st Session, House Doc. 186, 45 p. U. S. Navy, Hydrographic Office, 1944, Breakers and surf, principles in forecasting: Pub. No. 234, 54 p. Walsh, D. E.; Reid, R. 0.; Bader, R. G., 1962, Wave refraction and wave energy on Cayo Arenas, Campeche Bank: Texas A. and M, Dept. of Oceano- graphy Ref. 62-6T, 62 p. Wilson, W. S., 1964, Improving a procedure for numerical construction of wave rays: Masters Thesis, School of Marine Science, College of William and Mary, Williamsburg, Virginia. TABLE 1.-Computer input specifications for waves of 4~sec period Ray Ray origin Grid origin Direction from which no. angle for ray travels computer (A) Cos) (Cr) i 20.0 Bo 50) Go 50 30 2 20.0 RsAs (O00 30 3 20.0 39.25 75.50 30 4 240.0 52 7/5100 30 5 2h0 .0 LOO Who 5O 30 6 20.0 MLa i FAC 30 i 20.0 ne) ih Ges} 0) 30 8 210.0 2659239. 23 60 9 210.0 PMS 3S38 60 10 210.0 BO Si 52 60 afin 210.0 28.45 36.66 60 12 210.0 28.96 35.80 60 ali 210.0 29.47 34.94 60 14 210.0 29.98 34.08 60 15 210.0 BOgMe) 3-22 60 16 210.0 31.00) 32.36 60 iY 210.0 SL G0) Subs SO 60 18 180.0 250 26.150 90 19 180.0 PAO 25550 90 20 180.0 250) 2.50 90 Pill 180.0 21.50 23.50 90 22 180.0 Pilie50 22,50 90 23 180.0 PAL 50) 2, SO) 90 pn 180.0 PS Om ee2 On © 90 25 180.0 P50). 1G. 50 90 26 180.0 250 UIo50 90 oT 180.0 Pio 50) Io 50 90 28 120.0 siscil OHC6 150 29 120.0 oO Oi 150 30 120.0 Oks) O84 Si/ 150 31 120.0 1622 OSES 150 32 120.0 15.36 02.99 150 33 120.0 14.50 02.50 150 Table 2.-Computer input specifications for waves of 6-sec period. Ray Ray origin Grid origin Direction from which no. angle for ray travels computer (A) (x) (x) (ME) 34 210.0 85.50 73.50 60 35 210.0 SSO, 7Aseh 60 36 210.0 G52 Wil 60 37 210.0 87.03 70.92 60 38 210.0 87.54 70.06 60 39 210.0 88.05 69.20 60 ho 210.0 88.56 68.34 60 Ta 210.0 89.07 67.48 60 he 210.0 89.58 66.62 60 43 180.0 Ci, 50 2S.50 90 yh SOO Ci, 50 25050 90 45 180.0 CisSO) Des 50 90 h6 180.0 91.50 23.50 90 LT 180.0 91.50 22.50 90 48 180.0 GiLe5O ilo 50 90 hg 180.0 91.50 20.50 90 50 180.0 91.50 19.50 90 bil 180.0 Co 50) Msig 59 90 52 180.0 91.50 17.50 90 Gl. FIGURES Text Map of portion of mid-Atlantic bight showing generalized bathymetry and areas covered by depth grid and detailed ray diagrams. Map showing contours of depth values associated with intersection points of the primary grid, and the starting points and directions of the 52 wave rays run with the computer programs. Wave-ray diagram for rays numbered 1-7 (fig. 2). Wave-ray diagram for rays numbered 8-17 (fig. 2). Wave-ray diagram for rays numbered 18-27 (fig. 2). Wave-ray diagram for rays numbered 28-33 (fig. 2). Wave-ray diagram for rays numbered 34-2 (fig. 2). Wave-ray diagram for rays numbered 43-52 (fig. 2). Comparisons of wave rays 31-33 (fig. 2) constructed by three different programs for fitting velocity surfaces (the linear-surface program, A; the quadratic-surface program, B; and the forced-cubic-surface program, C) and by conventional graphical methods, D. Appendix G Map showing former location of U. S. Navy wave gage off Cape Henry, Virginia, in 20 feet of water. G28 Taree ATLANTIC OCEAN AREA OF DEPTH GRID § 500 \O is} =x m 3S 2m — CHESAPEAKE - —37°00' LIGHTSHIP a E> f £ : iS f) iS AREA OF — 36°30! RAY DIAGRAMS g | 75°00 %, FIGURE |. GENERALIZED DEPTH CONTOURS IN FATHOMS EXPLANATION so DEPTH CONTOUR IN FEET Sa CLOSED DEPRESSION 28 SMOOTHED SHORELINE ! WAVE RAY ORIGIN POINTS 2 DIRECTIONS OF TRAVEL 3 AND RAY NUMBERS 30—| DEPTH GRID TICK \_| LATITUDE OR 37°00 “LonGITUDE. TICK CAPE CHARLES CHESAPEAKE LIGHTSHIP. _____Match Line ATLA NT IE p ) OCEAN "4 Ces FIGURE 2. MAP SHOWING CONTOURS OF DEPTH VALUES ASSOCIATED o— WITH INTERSECTION POINTS OF THE PRIMARY GRID, AND THE STARTING POINTS AND DIRECTIONS OF THE 52 WAVE RAYS RUN WITH THE COMPUTER PROGRAMS Match Line Va pe -. North Property iD Camp Pendleton \S5th Street Fishing Pier Y=15 aly Dam Neck Mills e Sandbridge FIGURE 3 WAVE RAYS FOR 4 SEC. WAVES FROM 30° TRUE NORTH (2 (2 ISth Street Fishing Pier North Property oD Camp Pendleton FIGURE 4 veis WAVE RAYS FOR 4 SEC. WAVES Dam Neck Mills FROM 60° TRUE NORTH 7) Sandbridge 23 FIGURE 5 WAVE RAYS FOR 4 SEC. WAVES FROM) 902 TRUE SNORT 18 19 20 2 15th Street 22. Fishing Pier 23 25 North Property Line — 26 Camp Pendleton ig ll Dam Neck Mills Sandbridge 24 FIGURE 6 WAVE RAYS FOR 4 SEC. WAVES BROMESIS OL. aE RUE SNORE I5th Street Fishing Pier North Property (> Camp Pendleton Y=15 af Dam Neck Mills CY) Sandbridge 25 4o 15th Street Fishing Pier North Property eo Camp Pendleton FIGURE WAVE RAYS FOR 6 SEC. WAVES FROM 60° TRUE NORTH eee Dam Neck Mills = Sandbridge 26 FIGURE 8 WAVE RAYS FOR 6 SEC. WAVES FROMM IOS RUE SNORT I5th Street Fishing Pier North Property mo Camp Pendleton Y=15 as Dam Neck Mills 2 Sandbridge Cnt ») Lia dille) SES Elet IO J0vsYNS DILVYAGVNO HLIM NOILV10d44S1NI SAVY JAYM GQ3SLONYLSNOD SNIHOVW ol 423) U! sunojuoD yidag apis 0 uo 432) Ob0E OF 1189 plo +91UT 3apny ols sis 02 E 4a1ul we aapny Br 02 ol Ss ~—— X A ¢ WS 483) Ul sinojuoD yydag \ apis 0 uo [| 429) ObOE 1199 iia Oo a o fo} a P. *s A ila 1838 JO S39VSHYNS YVANII HLIM NOILV10d44LNI ps SAVY JAAVM G3LOINYLSNOD ANIHOVW 28 SAOVAUNS ALIDOTSA ONILLIS YO4 SWVY9OUd LN3YS4dIG SSYHL Ad G3LONYLSNOD (2 aNbi4) E¢-1€ SAVY JAVM 4O SNOSIYVdWOD 6 AYNDI4 | I [ T [T I | og SI é ol S 02 il Ss . 4X =<~——— x A \ A [== y Gh = ee 49a) Ul Sinoju0D yydaqg | oS aaj Ul SinojuoD yydag \ 1) Sete \ rs sce 4393 ObOE 1139 Bea 1139 roe > b L 9 os a = \ si K lhe 491ul oer Fe aapny 3 a % ° Pe ®. P \ \\ Ne / / \ = \\ sz— Le ~ 5 © Z30v4yNS DIGND GjAdYO4S HLIM NOILV10dY43SLNI SAVY JAVM Sim N oc = GSLONYLSNOD ATIVOIHd VHS SAVY JAVM quoi G3LONYLSNOD ANIHOVW AiuaH Dd ® l Sos l 1 ZY) APPENDIX A Computer Program for Computation of Wave Velocities as a Function of Wave Period and Water Depth PROGRAM TITLE: COMPV. Variables Used in Program: INCHES So d5000 --- Number of different wave periods for which velocity computa- tions are to be made. Mente eae npeciodm (seconds) KADEP eerie Wateridep thm cnecty re CXX(K)........ Wave velocity (feet/second) for water of K depth. XL.......2.2-- 1/2 deep-water wave length (feet) for a wave of TT period. L..eeeee0e-ee-- XL YOUNnded down to the nearest integer. Summary of Program: NOTT and the first TT are the first variables input to the computer. For water of K depth, where K ranges from 1 to L in one-foot increments, CXX's are then computed by an iteration pro¢edure so that the third digit to the right of the decimal is significant. After all CXX's have been computed for the first TT, they are rounded to the nearest hundredth and output. Computations then proceed using the next TT. After computations have been made using NOTT different TT's, the program stops. Remarks: If wave periods greater than 20 seconds are to be used with this program, CXX will require larger dimensions. The output from this program consists of a table of water depths and asso ciated wave velocities for each specified wave period. The CXX values of these serve as input for the DISTV program in Appendix B. 30 The COMPV program, as well as several of the following programs, utilizes a subroutine for the internal rounding of values before they are output. This subroutine (ROUND) may profitably be described at this point. SUBROUTINE TITLE: ROUND .« Variables Used in Subroutine: VALUE.......-. Corresponds to the value in the calling program which is to be output. WIRGooooGCOK0R0 10% where N is the number of digits which are to be output to the right of the decimal. Summary of Subroutine: If there are N digits to the right of the decimal in the output specification for VALUE, ROUND tests the N+lst digit to see if it is equal to or greater than 5. If it is, the Nth digit is rounded up. Otherwise, the Nth digit remains unchanged. | Remarks: Mode specifications, f and k, for ROUND must be idential with those for the calling program; k for ROUND must be one greater than the maximum number of digits to be output for any VALUE specified by the calling program. 3| XXKXXXXXXXXKXXKKXXKXKXXKXKXKXKXXKXXKXKX SOURCE PROGRAM XXXXXXXXXXXXXXXKXXXKK KKK KXKXAMK AX * 8 8 € 1620s FORTRAN IIs COMPUTATION OF WAVE VELOCITIES(FEET/SECOND) AS COMPV O1 C A FUNCTION OF WATER DEPTH(FEET) AND WAVE PERIOD( SECONDS) e COMPV 02 E THEORY FROM HeOcoPUBe234(1944)s PROGRAMED BY WeSeWILSON» COMPV 03 € JUNE 179 19646 COMPV 04 DIMENSION CXX (1025) COMPV 05 TANHF(X) = (EXPF(X)-EXPF(-X)) /(EXPF(X)+EXPF(=X)) COMPV 06 P = 361415927 COMPV 07 G = 3202 COMPV 08 READ 10sNOTT COMPV 09 10 FORMAT (13) COMPV 10 DO 1000 NOT=1sNOTT COMPV 11 READ 209TT COMPV 12 20 FORMAT (F5el) COMPV 13 XL = Oe5#G*(TT##200)/(200*P) COMPV 14 Le xe COMPV 15 CXXO = TT*G/(2e0*P) COMPV 16 EEE =| 5ie5 COMPV 17 BAR = 2e0*P/TT COMPV 18 DO 2000 K=leL COMPy 19 DEP = K COMPV 20 DO 3000 M=1590 COMPyV 21 CXX(K) = CXXO*TANHF( (BAR¥DEP)/CCC) COMPV 22 IF (ABSF(CXX(K)-CCC)-e0005) 55300053000 COMPyV 23 3000 CCC = (CXX(K)+CCC)/200 COMPY 24 5 IF (SENSE SWITCH 1) 493 COMPy 25 4 TYPE 900sKsM COMPV 26 900 FORMAT (2HK=91523H»M=913) COMPyV 27 3 VALUE = CXX(K) COMPV 28 CALL ROUND (VALUE9100¢) COMPV 29 2000 CXX(K) = VALUE COMPV 30 PUNCH 1009 TTosL COMPV 31 100 FORMAT (8HPERIOD =sF5e1ls25H SECONDS» MAXIMUM DEPTH =9I597He FEETe)COMPV 32 PUNCH 200 COMPV 33 200 FORMAT (/5(5HDEPTHs1X »6HVELCTY 93X)/) COMPV 34 PUNCH 3009(K»oCXX(K) s9K=1l9L) COMPV 35 300 FORMAT (5(159F70e293X)) COMPV 36 PUNCH 700 COMPV 37 700 FORMAT (///) COMPV 38 1000 CONTINUE COMPV 39 TYPE 800 COMPV 40 800 FORMAT (16HTHIS IS THE ENDe) COMPV 41 END COMPV 42 * 8 8 SUBROUTINE ROUND (VALUE »DEC) ROUND 01 ¢ PROGRAMED BY WeSeWILSON»s JULY 189 19646 ROUND 02 IVALUE = VALUE * DEC * 106 ROUND 03 IF (IVALUE) 10051045100 ROUND 04 104 VALUE = 060 ROUND 05 32 100 101 102 103 GO TO 103 JVALUE * VALUE * DEC XVALUE = IVALUE XVALUE = XVALUE / 1006 YVALUE = JVALUE IF ((XVALUE = YVALUE) = 005) VALUE = (YVALUE + le) GO TO 103 VALUE = YVALUE / DEC RETURN END / DEC 99.9, 0.0.9.9.0.0.9.9.0.099.9. 0.0.9.2. 0.0.0,0.0.0.0.0.0.00.0.0.04 1 400 XXX XXKXK KKK KK AX XK XK XK XX KK XK KX KKK KXXKXX OUTPUT PERIOD = DEPTH VELCTY 1 5060 6 12283 ll 16017 16 18610 21 19022 26 19084 31 20017 36 20234 420 SECONDS» MAXIMUM DEPTH = DEPTH VELCTY 7082 13.67 16064 18237 19637 19293 20022 20236 DEPTH 10291019101 INPUT ROUND ROUND ROUND ROUND ROUND ROUND ROUND ROUND ROUND ROUND ROUND KK KH KK HK HK IK KH HK HK MH MM HK IKK HMM HK KK HM KK NOTT UU XM HK MH KKK MA MK KH NHK HH KKK NM HH KK KK KKH KK 40.0 FEETo VELCTY DEPTH 904§ & 14040 9 17.07 14 18062 19 19.51 24 20200 29 20026 34 20038 39 33 VELCTY 10077 15006 17045 18284 19064 20007 20029 20040 DEPTH VELCTY 11287 15065 17079 192004 19675 20012 20032 20041 APPENDIX B Computer Program for Distribution of Wave Velocity Values Over a Grid of Depth Values as a Function of Wave Period PROGRAM TITLE: DISTV. Variables Used in Program: Maiveyetieerien Naver Deriodms (seconds). Dodie ceeecee ses NOM a) Sven 4 pOsition) Onvar era GW) aXe lie Vee eaees where the grid origin is at X = 0, ¥Y = 0. CMAT(I,J)..... For position I,J on the grid, CMAT represents, in the input, water depth (feet or fathoms). After conversion, CMAT re- presents, in the output, wave velocity (grid units/second). MM,NN......... Maximum I and J, respectively, for the grid. FMOP.......... Allows conversion of CMAT to feet, if input is in fathoms. GRED. se ossee se Grid interval (feet): Ingo6 6G OOOO OO 6t 1/2 deep-water wave length (feet) rounded down to the nearest integer. (CXX(K),K=1,L) For a wave of TT period, CXX is the array of wave velocities (feet/second) where K ranges from 1 to L in one-foot incre- ments. JX,JY...--2--- On each output card JX and JY represents the X and Y grid co- ordinates of the first velocity value on that same card. Summary Of Program: MM, NN, FMOP,GRID,TT, and (CXX(K),K=1,L) are the initial input. Then the first row of CMAT depth values is read into the computer. (Fathoms are converted if specified by FMOP.) For each depth value the computer "looks up" the associated CXX and divides this value by GRID. When all values in 34 the first row have been converted from a depth in feet to a velocity in grid units/second, they are rounded to the eighth place to the right of the decimal and output. When all rows for a given grid have been similarly input, converted, and output, the computer pauses in order that another set of data may be input. Remarks: Each output card includes (at the extreme right-hand side) TT, JX, and JY. Depths greater and L are assigned wave velocities equal to those assigned for L. The output from this program serves as input for MAIN 1620 and MAIN 7094 (Appendices D and £), MM must be a multiple of 10. 35 XX XK XX KK XXX KK X KKK AXXR XK KAKA KKK AKKX SOURCE PROGRAM XXAXKXKXKKKXAXKAKAAKKKAK AXKAAKAARK AK * 8 8 Cc 16205 FORTRAN IIs DISTRIBUTION OF WAVE VELOCITIES(GRID UNITS/ DISTV 01 ig SECOND) OVER A GRID OF DEPTHS(FEET OR FATHOMS) AS A FUNCTION DISTV 02 G OF WAVE PERIOD(SECe) »PROGRAMED BY WeSeWILSON® MAY 19919646 DISTV 03 DIMENSION CXX(1025) »CMAT (200) DISTV 04 10 READ 409 MM»oNNoFMOP »GRIDsTT DIstv 05 40 FORMAT (I149T49F2e09F7e09F501) DISTV 06 XL = 0¢5*50118%( TT*#*200) DISTV O7 SG Sede DISTV 08 READ 209 (CXX(K) sK=loeL) DISTV 09 20 FORMAT (5(5X9F70293X)) DISTV 10 PUNCH 100sTT»sGRID DISTv 11 100 FORMAT (8HPERIOD =9F5e1ls17H SECesGRID SIZE =eF7e0s6H FEETe/) DISTV 12 DO 3000 J=1sNN DISTV 13 READ 3059 (CMAT(1I)oI=19MM) DISTV 14 30 FORMAT (10F4e1) OISTv 15 IF (FMOP) 29291 DISTV 16 1 DO 1000 1=19MM DISTV 17 1000 CMAT(I) = 6e0*CMAT(T) DIsTv 18 2 DO 2000 I=1 MM DISTv 19 NVALUE = 1 DISTV 20 IF (CMAT(1T)) 39200094 DISTV 21 3 CMAT(T) = -(CMAT(TI)) DISTV 22 NVALUE = 2 DISTV 23 4 IF (CMAT(TI)=XL) 69695 DISTV 24 5 CMAT(I) = XL DISTV 25 6 XK = CMAT(T) DISTV 26 K = XK DISTV 27 CMAT(I) = CXX(K)/GRID DISTV 28 CALL ROUND (CMAT(I)910¢E8) DISTV 29 GO TO (200097) sNVALUE DISTV 30 7 CMAT(I) = =(CMAT(T)) DISTV 31 2000 CONTINUE DISTV 32 MM5 = MM/5 DISTv 33 if al DISTV 34 BDO 5000 JM=1>5MM5 DISTv 35 JX = Il-1 DISTV 36 JY Oo dA DISTV 37 III = 11+4 DISTV 38 PUNCH 200s(CMAT(I) sIT2IIs III) sTT»IXsJY DISTV 39 200 FORMAT (5(F100893X) 92X9FS5elsI40I4) DISTv 40 Mit ES Shieh teal DISTV 41 5000 CONTINUE DISTV 42 3000 CONTINUE : DISTV 43 TYPE 4009TT DISTV 44 400 FORMAT (F50e1920H SECeGRID COMPLETEDe) DISTV 45 PAUSE DISTV 46 PUNCH 300 DISTV 47 36 300 FORMAT (///) GO TO 10 END XX XKK KKK MK KKK KK HM IK KKK KICK XK OOK XK OOKCOK 20 3 0 1 5260 6 120683 Te LiGrenl 7 16 18.610 Aa | WO 26 19284 31. 20017 36 20034 -00 =-00 -00 -32 -24 -09 -00 -00 -00 -32 -22 008 -00 -00 -00 -30 -15 020 XK KKK KKH KKK HK XXX KKK KK KKK KK XK KKK XKXKKX OUTPUT 3040. 460 2 7082 7 13267 12 16264 UY UWISsy 22 19437 27 19293 32 20022 37 20636 =-00 -00 -00 -00 027 030 040 043 -00 -00 -00 -00 025 030 039 043 -00 -00 -00 -00 030 036 042 045 -00 045 -00 045 -00 049 PERIOD = 420 SECesGRID SIZE = 3040 0200000000 0200000000 0200000000 0200000000 0200000000 0200000000 -00665132 — 000646053 - 000495395 © 00671382 000671382 000671382 0-00000000 0200000000 0200000000 0200000000 0200000000 0200000000 - 200665132 -200637171 000473684 e 00671053 000671382 000671382 0200000000 0200000000 0200000000 0.00000000 0200000000 0200000000 —.00661842 -000585197 000626316 © 00671382 000671382 000671382 IN 9 14 17 18 19 20 20 20 41 055 43 056 -43 057 e F PUT 045 4 040 9 e207 14 062 19 051 24 200 29 026 34 238 39 EETe 0200000000 0200000000 000655592 ©00671382 0200000000 0.00000000 090649671 °00671382 02000009000 - 200671382 ©00661842 000671382 37 10677 15206 17045 18284 19064 20¢07 20029 20¢40 0200000000 -000671382 ©00661842 ©00671382 0200000000 -000671382 000661842 000671382 0200000000 - 200671382 ©00669079 ©00671382 DISTV 48 DISTV 49 DISTV 50 KKK KKM KKK NK KK KK KH KK HK KH KK HI HK IK MK HK KKK MM »9NNoFMOP sGRIDsTT 5 11¢87 10 15265 NEY yoy 20 19204 AS sos 30 20¢12 35 20032 40 20041 0000 1000 0001 1001 0002 1002 420 400 HK KKK KKK KKM HK HK HK HICH HN IH HHI IK ION IK HK NNUNNPRPRP RR OOOO APPENDIX C Computer Program to Produce Matrices for Deriving Equations of Linear Surfaces PROGRAM TITLE: PRMAT. Variables Used in Program: (X(I),I=1,4).. The X co-ordinates of the four data points. (Y(I),I=1,4).. The ¥Y co-ordinates of the four data points. SoHMoo600c060d . Matrices used in determiningrthe coefficients for an equation of a linear surface of best-fit to data values at four points. Summary of Program: Let C represent wave velocity at any grid position X,Y. By specifying the X,Y values for four positions, this program outputs S and EM. When later manipulated with a particular set of C values for these same four positions, S and EM produce the coefficients (E's) of linear equation of the form, Ca E, + EX + Det This equation represents the least-squares plane of best-fit to this particular set of C values. Remarks: This program was extracted from a general program provided by W. C. Krumbein that fits a first, second, or third-order polynomial surface to specified regularly-spaced data points by the least-squares method. The output from this program serves as spate Oh MAIN 1620 and MAIN 7094 (Appendices D and £). 38 MX KK KKK KKK KKK KKK XX KK XK KX XXX XMKX SOURCE PROGRAM XXXKXXXKXKKX KKK KK XK MM HMM KK HHH ¥*16 6 Cc 16209 FORTRAN II» PRODUCE MATRICES FOR DERIVING EQUATIONS OF PRMAT O01 C LINEAR SURFACESe PRMAT 02 < THEORY FROM WeCeKRUMBEIN» PROGRAMED BY BETTY BENSON» PRMAT 03 G MODIFIED BY WeSeWILSONs JULY 175 19640 PRMAT 04 DIMENSION X(4)9 Y(4)9 EM(45s3)9 S(3093) PRMAT 05 READ 109 (X(I)s9T#194) PRMAT 06 READ 109 (Y(I)soT=194) PRMAT 07 10 FORMAT (4F5e0) PRMAT 08 DO 1000 L=1s4 PRMAT 09 EM(Lo1) = le PRMAT 10 EM(L92) = X(L) PRMAT 11 1000 EM(L»93) = Y(L) PRMAT 12 DO 2000 [=1s3 PRMAT 13 DO 2000 J=1s3 PRMAT 14 S(IsJ) = 000 PRMAT 15 DO 3000 L=194 PRMAT 16 3000 S(IsJ) = S(IsJ) +EM(LoT) #EM(L9J) PRMAT 17 2000 CONTINUE PRMAT 18 DO 4000 K=193 PRMAT 19 DIV = S(KsK) PRMAT 20 S(K9K) = 160 PRMAT 21 DO 5000 J=153 PRMAT 22 5000 S(KsJ) = S(KsJ)/DIV PRMAT 23 DO 4000 I=193 PRMAT 24 IF (T=K) 19400091 PRMAT 25 1 DIV = S(ToK) PRMAT 26 S{I9K) = 000 PRMAT 27 DO 6000 J=193 PRMAT 28 6000 S(IsJ) = S(IsJ)I-DIV¥S(K9J) PRMAT 29 4000 CONTINUE PRMAT 30 PUNCH 500 PRMAT 31 500 FORMAT (SOHMATRICES FOR DERIVING EQUATIONS OF LINEAR SURFACES) PRMAT 32 PUNCH 1005 ((S(TsJ) 9J=193) 912193) PRMAT 33 100 FORMAT (6F 128) PRMAT 34 PUNCH 2009 (CEM(LoI) oL=194)sI=193) PRMAT 35 200 FORMAT (12F6e02) PRMAT 36 END PRMAT 37 XH XX KK KK XK KKK KKK KKK MK RK XX KK KXKXKXKKK TNPUT KXMKKX KKK KKK MH HK KM HN HK KHIM KKM KR KK KK HX 0.0 1.0 120 0.0 (X( 1) oT=154) 060 060 160 160 (YCT) oT2=194) XX KX XXX KKK KX XXX KK AKK KK KK XK KKKKXKKKKK OUTPUT NXRXXKKKKKKK KK KK MRK KKK HNN NK HK KK HK KKK KK MATRICES FOR DERIVING EQUATIONS OF LINEAR SURFACES 075000000 =e50000000 =-.50000000 =.50000000 1.00000060 0.00000000 =250000000 0.00000000 14200000000 1600 1600 1600 100 0200 1-00 1200 0600 0200 0200 1600 1-00 39 APPENDIX D Computer Program for Wave Refraction (Linear-Interpolation) Using the IBM 1620 PROGRAM TITLE: MAIN 1620. Input Variables: S,HM.......... Matrices used in determining coefficients for the plane of best-fit to four data points. XLABL......... Arbitrary title used to designate each batch of input. MM NN. sseece . Maximum values of I and J, respectively, for the CMAT grid of vellocity values; T=X +1, J =Y + 1, with the grid origin lee@ecech G65 26 = On Mis Oc CHECK... .ccces Allows either a two-dimensional CMAT to be input, or a single- rowed CMAT (extending from shore to deep water) to be input and all other rows to be made identical to the first. (This procedure creates a set of grid values which can be character- ized by straight and parallel contours. ) TBs Sooo .sss.. Wave period (seconds). NOJ.......---- Number of wave rays to be run in each batch. D..eeeeeceeese- Distance incremented between successive points along a ray path. CMAT(I,J)..... Wave velocity (grid fiateyeccen a) at grid position I,J. A..ecccseeeees Angle (degrees) measured from the direction of increasing X along the X-axis moving counter-clockwise to the dire¢tion of travel of a wave ray. X,Y........-.- Grid origin position for a wave ray. 40 Output Variables: WIWABIS Goo G0000 Defined previously. TR h9 OltiOIOO.G Defined previously. INOIWs 6660 see... Davch number. Woosodana00000 Ray number for batch NOT. MIX cocaod0000 - Number assigned to a point along a ray path where calculations are made. MAX = 1 at the origin point of ray. FOOL, WAR 506000 X,Y co-ordinates for a given MAX. ANGLE......... Angle (degrees) between X-axis and ray, as previously defined for A. PCTDIF.....+-. Percent difference between the value of one of the four points to which a plane is fit and that corresponding value of the plane (see text, "Interpolation Surfaces"). Variables in Common: S, HM.......-. Defined previously. Heesceeeeeeees Coefficients of the equation of a plane of best-fit to four grid points. YVW......-.... Matrix used in determining E's. CWE gc00G0006 Defined previously. Crccccceeeeeee Values of the four CMAT values to which a plane is fit. XLABL......... Defined previously. D.....--2----- Defined previously. IM6 60000000 -.-. Defined previously. QiMooac0o00000 Wave-velocity for a given MAX. Io aga00000000 Number of iterations used in obtaining the curvature of a wave ray for a given MAX. NGO........... Designates whether a wave ray has moved within one grid unit 4\ of the edge of the grid. AMM, ANN....... Used in determining NGO. MAX......-.-.. Defined previously. MUEPG 6 od0000 -.. Designates whether the last two curvature estimates for a given MAX are less than 0.00009/D, whether the 18th and the 20th estimates are less than this value, or whether neither of the above is true. Summary of Program: MAIN reads S and EM and sets NOT = 1 to denote batch one. It then reads the information (XLABL, MM, NN, CHECK, TT, NOT, D, and CMAT) for processing the first batch of rays. For N = 1 the program reads the information (Gi, 35° ekavel A) for processing the first wave ray. MAIN converts A from degrees to radians and then punches the title information (XLABL, TT, NOT, N, and the column headings). Control is then transferred to the RAYN subroutine. When RAYN has determined the path of the first ray, control is transferred back to MAIN. MAIN then reads the information for the second ray (N = 2) and proceeds as before. When NOJ ray paths have been determined for NOT = 1, the computer pauses allowing time to load the necessary information for the second batch. After pressing start, MAIN sets NOT = 2, and the information for this batch is read. MAIN then continues as before. The program is terminated when the desired number of batches has been processed. Remarks: 5; This program calls four subroutine/RAYN, SURFCE, MOVE, and ROUND. Descrip- tions of these subroutines follow. (ROUND has been described previously. ) 42 SUBROUTINE TITLE: RAYN. Summary of Subroutine: RAYN calls SURFCE to obtain FK (ray curvature) for MAX = 1. MAX is then set equal to two. RAYN calls MOVE in order that X, Y, A, and FK be obtained for MAX = 2. ANGLE (equal to A but expressed in degrees instead of radians) and PCTDIF are then computed. XXX and YYY are set equal to X and Y, respectively, so that significance is not lost when the values to be output are rounded. ROUND is called for XXX, YYY, ANGLE, and PCTDIF before they are punched. MAX is incremented by one, and the entire process continues until 1) MAX = 500 , 2). CXY is' equal to or less than zero, 3) MIT = 3, or 4) NGO = 3. The first condition guards against a ray going in an endless circle. The second pre- vents a ray from leaving the water and going over land. The third prevents invalid curvature approximations. The fourth condition stops a ray if it reaches the edge of the grid. Any of these four conditions causes the message, "RAY STOPPED" to be punched and transfers control back to MAIN. If MIT = 2, the message "CURVATURE APPROXIMATED" is punched and computations proceed as usual. Remarks: The use of Sense Switch 3 keeps the operator informed of the current value of MAX. SUBROUTINE TITLE: SURFCE. Variables Used in Subroutine: I,FI......+.-- X + 1 rounded down to the nearest integer. J,FJe.e+eee--- Y + 1 rounded down to the nearest integer. Wbooaocco00000 o& we dh o Wilko Wabooccaccc0000 er do Wo 43 We UN 00000000 Values of FI and FJ, respectively, the last time SURFCE was called. (Not valid for MAX = 1.) IMK&5 g00000G0000 Curvature of the wave ray at X,Y. Summary of Subroutine: For a specified A and X,Y position on the CMAT grid, SURFCE first determines I, J, XL, and YL. If MAX does not equal one, SURFCE checks to see if ZI = FI and ZJ =FJ. If both equalities are true, the # values from the previous operation of SURFCE are still valid, and the CXY and FK computations can be made directly. Otherwise, it is necessary to derive new E values, using C, IM, and S, before computing CXY and FK. Remarks: Program steps corresponding to SURFCE2O0 through SURCH2/7 in the listing that follows were obtained from a program provided by Dr. W. C. Krumbein. (This program is described with the PRMAT program in Appendix Gs) SUBROUTINE TITLE: MOVE. Variables Used in Subroutine: FKBAR......-.-- Curvature used to obtain DELA for a given IT. TMM 6 obcGoad MINS yor ITU = AUS). ll FKKP.......... FKBAR for IT = n-l where current FKBAR is for IT =n. (For IT = 1, FKKP equals the FKBAR used in determining XX and YY the previous time MOVE was called.) PKK .cc-ceee-- X, Y, AS and HK wailues, respectively.) for) Max snr iewinercce MAX = n when MOVE was called. DELX. 2. ccee ~. XX = X. 44 Summary of Subroutine: If MAX equals three or more when MOVE is called, FKBAR for the previous MAX is used in obtaining approximations of XX, YY, and AA. (For MAX = 2, FKBAR is set equal to FK.) SURFCE is called and returns FKK for this approximation at XX,YY. KBAR is then redefined as (FK+FKK)/2. If the difference between FKBAR and FKKP is less than 0..00009/D, the current XX, YY, AA, and FKK values are accepted for the new point. If the difference is greater, FKKP is set equal to FKBAR, and the current FKBAR is used to obtain another set of XX, YY, AA, and FKK values. The difference between FKBAR and FKKP is again tested. This cycle may repeat a maximum of 20 times before termination. If the cycle stops before IT = 20, MET is set equal to one. If the cycle stops: at IT = 20, and if the difference between FKBAR and FKKPP is less than 0.00009/D, then MIT is set equal to two, and FKBAR is defined as (FKBAR+FKKP)/2 for obtaining XX, YY, AAA, and FKK. If IT = 20, and this difference is greater than 0.00009/D, then MIT is set equal to three. When MIT = 3, control is transferred back to RAYN immediately. When MIT = 1 or 2, XX and YY are tested to see if the new point has reached the édge of the grid. NGO = 2 if the ray has reached the edge of the grid, and NGO = 1 if it has not. Remarks: MIT = 1 when the curvature approximations are converging to a single vaiiue. MIT = 2 when the approximations are oscillating between two values. In this case the average of the two values is taken as the new curvature. MIT = 3 when the approximations are oscillating among three or more values. No valid curvature approximation can be made in this case; this is one basis for the terminiation of a ray. The use of Sense Switch 2 allows the operator to observe successive IT and 45 FKBAR values. The maximum difference (0.00009/D) in the curvature approximations is such that the output ANGLE is significant in the hundredths digit. 46 KX KH KKK MK KKM HHH HMM HK KK HK HK KKK HK HK KK NA * 403 FORMAT (///12A4/8HPERIOD =9F50e1s6H SECes910H BATCH NOceI399H» RAY 1391He//4X 9 3HMAX 9 6X 9 LHX 9 8X 9 LHY 9 8X9 SHANGLE 9 4X 9 6GHPCTDIF 9 2//1792F9e29Fl1lle2) 15 3 INOec» 1620s FORTRAN IIs LINEAR-INTERPOLATION WAVE REFRACTION PROGRAMe ORIGINAL PROGRAM BY GRISWOLDsNAGLEsAND MEHRe THIS PROGRAM ADAPTED FROM ORIGINAL BY WILSON»HARRISON »sKRUMBEINs AND BENSONe DIMENSION S(393) 5EM( 493) 5E(3) »¥YVW(3) sCMAT (50951) 9C(4) »XLABL(12) COMMON SsEMsEsYVWsCMAT sC oXLABL »DoTT »CXY 9IT 9NGO 9 AMMo ANN 9 MAX o MIT READ 59(¢€S(I9J) »J=193) s1=193) FORMAT (6F1208) READ 7o((CEM(LoI) oL=1904) »I=193) FORMAT (12F 6602) NOT = 1 READ 4009 XLABL FORMAT( 12A4) READ 402 9MMoNN» CHECK 9 TT oNOJ 9D FORMAT (2149F3ce097X9F50e19159F40l) AMM = MM=1 ANN = NNe1l IF (CHECK) 10910920 READ 119((CMAT(I 9J) »IT=19MM) »J=1 9NN) FORMAT (5(F100893X)) GO TO 14 Jzil READ lls (CMAT(I09J) 9I=159MM) DO 77 J=2o9NN DO 77 I=1l9MM CMAT(T9J) = CMAT(T91) DO 15 N=l»5 NOJ READ 69A0XsY FORMAT (F70e292F6e02) MAX = 1 PUNCH 4039XLABL»oTT »NOT oN oMAXoX9YOA A=A¥ 00174532925 CALL RAYN (X9YsA) CONTINUE PAUSE NOT = NOT + 1 GO TO 9998 END SUBROUTINE RAYN (XsY9A) DIMENSION S(393) 9EM( 493) 9E(3) sYVW(3) »CMAT(50951) 9C(4) » XLABL(12) COMMON SoEMsEsYVWsCMAT 9C »XLABL sDoTT sCXYsIT #»NGO» AMM 9 ANN 9 MAX o MIT CALL SURFCE (XsYsAsFK) MAX=1+MAX 47 MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN RAYN RAYN RAYN RAYN RAY¥N RAYN SOURCE PROGRAM XXXXXKKXXXKXXKKMAKK MAKKAH KM KK KX HA KK NK 100 103 101 399 396 395 200 397 12 13 8 6 MPN 318 IF (SENSE SWITCH 3) 1002101 TYPE 103 9MAX FORMAT (4HMAX=914) IF (MAX=500) 399915915 CALL MOVE (XsYoAsFK) TF -CEXY) 2159159396 GO TO (3979395515)s MIT WRITE 2009 MAX FORMAT (32HCURVATURE APPROXIMATED FOR MAX =914) ANGLE=A*57e29577951 XXX = X WOAL Evg PGTDINE = ABSE CUG(3) -E CRE C2). =F 3) G03) 21 0 Oe CALL ROUND (XXX9100e) CALL ROUND (YYY:100e) CALL ROUND (ANGLE »1000e) CALL ROUND (PCTDIF 9100) PUNCH 129MAX9XXX9YYY o9ANGLE sPCTDIF FORMAT (I1792F9e29Flle29F10e0e1) GO TO (3915) »NGO PUNCH 13 FORMAT (12HRAY STOPPEDe) RETURN END SUBROUTINE SURFCE (XsYsAo9FK) DIMENSION S(393) 2EM(493) E(3) sYVW(3) sCMAT (50951) 9C(4) sXLABL(12) COMMON SoEMoEsYVWsCMAT »CsXLABL »DoTT »CXY9IT»NGOsAMM oANN 9 MAX o MIT I1=X+lo J=Y+1le FI=I FJ=J XL=X+1le-FI YL=Y+1le-FJ IF (MAX=1) 19194 VE CZT—-F1) Lsi2is1 WE (ZI=F)) 15391 Zao Zj) SB Fd) C(1)=CMAT(I9J) C(2)=CMAT(I4+1sJ) C(3)=CMAT(I4+19J+1) C(4)=CMAT( I 9J+1) DO 318 Il=193 YVW(TT) = Oo DO 318 L=1s4 YVW(IT) = YVWOITT)4C(L)*EM( Lo IT) DO 319 II=193 ECII) = Oc 48 RAYN O7 RAYN 08 RAYN 09 RAYN 10 RAYN 11 RAYN 12 RAYN 13 RAYN 14 RAYN 15 RAYN 16 RAYN 17 RAYN 18 RAYN 19 RAYN 20 RAYN 21 RAYN 22 RAYN 23 RAYN 24 RAYN 25 RAYN 26 RAYN 27 RAYN 28 RAYN 29 RAYN 30 SURFCEO1 SURFCE02 SURF CE03 SURFCE04 SURFCE05 SURFCE06 SURFCEO7 SURFCEO8 SURFCEO9 SURFCE10 SURFCE11 SURFCE12 SURFCE13 SURFCE14 SURFCE15 SURFCE16 SURFCE17 SURFCE18 SURFCE19 SURFCE20 SURFCE21 SURFCE22 SURFCE23 SURFCE24 SURFCE25 319 102 104 39 On Ww 38 DO 319 JJ=193 ECII) = ECIIV+S(TIsJJ)*#YVW( JJ) CXY = E(1) + E(2)¥XL + EU3)*YL FK =(E(2)*SINF(A)-E(3)*COSF(A)I/CXY RETURN END SUBROUTINE MOVE (XsYsAo0FK) DIMENSION S(393) 0EM( 4093) 9E(3) sYVW(3) sCMAT(50951) 0€{4) oXLABL(12) COMMON SoEMoEsYVWoCMAT oC sXLABL sDs TT sCXY oI TsNGOoAMMoANR oMAX oMIT IF (MAX = 2) 10291029104 FKBAR=FK MIT = 1 DO 20 IT#1920 DELA=FKBAR*D AA=A+DELA ABAR=A+e5*DELA DELX=D#COSF (ABAR) DELY=D#SINF(ABAR) XX=X+DELX YY=Y+DELY GO TO (10196) MIT CALL SURFCE (XX9YYs2AAs9FKK) IF (CXY) 38938910 FKBAR = 00e5 *® (FK + FKK) IF (SENSE SWITCH 2) 8989899 TYPE 9009 ITsFKBAR FORMAT (149E14-8) IF (IT = 18) 593799 FKKPP = FKBAR TF (MAX = 2) 79709 IF (IT = 1) 2092009 IF (ABSF(FKKP=FKBAR) = (0-00009/D)) 696920 FKKP = FKBAR IF (ABSF(FKKPP = FKBAR) = (0200009/0)) 18518917 MIT = 3 GO TO 38 FKBAR = 0065 * (FKBAR + FKKP) MIT = 2 GO TO 39 NGO = 1 IF ((XX—100)*( (AMM—10e00)—XX))29293 IF (€YY—1¢0)*( (ANN=100)-YY) 129298 NGO = 2 X = XxX Y 2 YY A = AA FK = FKK RETURN END 49 SURFCE26 SURFCE27 SURFCE28 SURFCE29 SURFCE30 SURFCE31 MOVE 01 MOVE 02 MOVE 03 MOVE 04 MOVE 05 MOVE 06 MOVE 07 MOVE 08 MOVE 09 MOVE 10 MOVE 11 MOVE 12 MOVE 13 MOVE 14 MOVE 15 MOVE 16 MOVE 17 MOVE 18 MOVE 19 MOVE 20 MOVE 21 MOVE 22 MOVE 23 MOVE 24 MOVE 25 MOVE 26 MOVE 27 MOVE 28 MOVE 29 MOVE 30 MOVE 31 MOVE 32 MOVE 33 MOVE 34 MOVE 35 MOVE 36 MOVE 37 MOVE 38 MOVE 39 MOVE 40 MOVE 41 MOVE 42 MOVE 43 MOVE 44 XXXXXXXXXXXXXXXXXXXK KX KK KK XX KK KK KKK 275000000 =.50000000 =250000000 000000000 NoGO WMEOO AGOO ” W5OO 1620s GRID OFF VAeCAPES» Ao) oe © 400 0200000000 0200000000 0200000000 0200000000 — 200665132 - 200646053 © 00671382 000671382 0200000000 0200000000 0200000000 0200000000 —.00665132 -200637171 ©00671053 000671382 0-00000000 0200000000 0.00000000 0.00000000 - 00661842 -.00585197 2° 00671382 000671382 0200000000 0200000000 0200000000 0200090000 = 000646053 000354276 200671382 000671382 0200000000 0200000000 0200000000 0200000000 -2 00612500 ©00514803 © 00671382 ©00671382 0200000000 0200000000 0200000000 0200000000 = 000495395 000637171 200671382 200671382 0200000000 0200000000 0200000000 -200671382 200422039 ©00660197 ©00671382 ©00671382 0200000000 0200000000 0200000000 -200669737 © 00531908 200666447 ©00671382 ©00671382 0200000000 0200000000 0200000000 - 200669737 © 00655592 000666447 © 00671382 200671382 000000000 0200000000 -.00669737 000666447 200663487 000668421 200671382 000671382 0200000000 0200000000 =. 00669079 —000666447 © 00665132 000668421 © 00668421 000671382 INPUT -~.50000000 -.50000000 1200000000 0.00 1200 1-00 0-00 LINEAR INTERPos 8/5/646 3 05 0e00000000 0-00000000 0600000000 0.00000000 = 200495395 200655592 200671382 200671382 0e00000000 0600000000 0200000000 0.00000000 000473684 2©00649671 200671382 200671382 0600000000 0.00000000 0200000000 -.00671382 © 00626316 200661842 200671382 200671382 9200000000 0.00000000 -.00671382 -.00671382 000657895 200669079 200671382 2©00671382 0e00000000 0.00000000 -000671382 -e00670395 200665132 200669737 000671382 200671382 0e00000000 0200000000 =200671382 -—.00670395 2000661842 200671382 200671382 200671382 0200000000 0.00000000 —000670395 200666447 200669079 °00671382 200671382 °00671382 02000000000 0.00000000 =000669737 —000666447 200670395 200671382 200671382 200671382 0200000000 0.00000000 ~.00668421 -.e00665132 200669737 200671382 200671382 200671382 0e00000000 0200000900 -000665132 000652632 200669737 200671382 200671382 200671382 0e00000000 0.00090000 =000665132 -000619737 2000670395 200671382 000671382 200671382 50 0200 1200000000 0200 0e00000000 -°00671382 ©00661842 °©00671382 0200000000 -—2©00671382 ©00661842 200671382 0e00000000 - 200671382 °00669079 ©00671382 0200000000 -— 200669079 000671382 ©00671382 0200000000 — 000666447 ©00671382 ©00671382 0200000000 —200661842 ©00671382 000671382 000000000 -200649671 ¢00671382 ©00671382 0200000000 -—e00612500 ©00671382 ©00671382 0200000000 — 000495395 °00671382 000671382 0200000000 °00390461 ©00671382 ©00671382 0200000000 ©00531908 ©00671382 ©00671382 0.00000000 1-00 1-600 40 FPP ELE eee 8 @ © e@ Co 9o00O0 00000 > e XXXKKXKXKKXKK RX KKK KX XX KX KX KKK KKK KAKA KKK KK Ll 0-00000000 — 200667434 © 00665132 e 00671382 0200000000 - 000666447 © 00665132 © 00669737 0200000000 — 200665132 © 00667434 e 00669737 0200000000 - 000655592 © 00665132 e00669079 0-00000000 -—2°00649671 e00665132 °00671382 0-00000000 — 200646053 © 00661842 © 00671382 0-00000000 -—-00649671 © 00665132 e 00671382 0.00000000 -—200649671 © 00667434 ©00671382 0200000000 -—¢00649671 ° 00665132 © 00671382 0200000000 -— 200641776 e 00663487 © 00671382 0200000000 -200637171 e 00665132 © 00671382 12060 14650 120060 15036 120¢0 16622 0200000000 0200000000 0200000000 0200000000 4.0 0 11 — 200666447 -000652632 -000547368 ©00649671 400 Bait 000666447 000669737 ©00671382 ©00669079 4-0 10 11 ©00671382 ©00671382 e00671382 e00671382 400 15 11 0200000000 0e00000000 000000000 - 00669079 400 ahr - 000665132 -000649671 ©00257237 ©00657895 400 Sel 000667434 ©00671053 ©00671382 ©00669079 400 10 12 ©00671382 ©00671382 ©00671382 ©00671382 GeO) eo 2 0200000000 0e00000000 0200000000 —000666447 40 @). - 342} -200661842 -— 000626316 ©00514803 ©00661842 400 Gy AN) 000666447 ©00669079 © 00671382 ©00669737 4e0 10 13 000671053 ©00671382 © 00671382 ©00671382 4¢0 15 13 0200000000 02e00000000 0200000000 — 000669737 40 Oo 14 -—000641776 -.00585197 ©00637171 000661842 4.0 5 14 000667434 ©00670395 000671382 000671382 4e¢0 10 14 000671053 ©00671382 © 00671382 ©00671382 4¢0 15 14 0200000000 0200000000 -—200668421 -000661842 420 @) 3G) -200637171 ~-2©00257237 e 00646053 ©00657895 400 5S 000667434 °00669079 «00671053 ©00671382 460 10 15 ©00671382 000671382 200671382 ©00671382 40.0 15 15 0200000000 0200000000 - 200665132 -—000657895 4.0 Oo 16 - 000626316 000514803 © 00655592 °00657895 400 3 Ue 000667434 ©00669079 © 00671382 ©00671382 4e-0 10 16 ©00671382 e00671382 ©00671382 °00671382 4-0 15 16 0200000000 0200000000 -— 200663487 —000657895 420 Oo 17 — 000626316 ©00547368 °©00660197 °00660197 400 Sli \ 200669079 ©00670395 °00671382 ©00671382 4e0, 10 17 000671382 00671382 °00671382 °00671382 4e0\ 15 17 0200000000 000000000 -— 200660197 -—200657895 4.0 @) alts) -\.00585197 ©00637171 © 00652632 ©00661842 \4e0 Sees 000668421 000669737 ©00670395 ©00671382 400 10 18 600671382 ©00671382 e00671382 ©00671382 ‘460 \15 18 0200000000 0200000000 -200657895 =000655592 4e0 @) oO) = 200422039 ©00649671 °00652632 °00660197 420 By 1S) 000666447 ©00671053 ©00671382 °00671382 42.0 \10 19 000671382 000671382 °00671382 ©00671382 400 15 | 19 0200000000 — 000663487 - 200655592 = 000652632 4-0 0 \20 0200000000 000655592 ©00655592 ©00655592 40 5 \20 e00665132 °00669079 °00671382 ©00671382 400 10 (20 ©00671382 000671382 ©00671382 °00671382 400 15 20 0200000000 -—200660197 —200655592 200649671 420 Oo 21 000547368 000655592 °00655592 °00655592 40 5) 2.1! 000667434 ©00670395 ©00671382 000671382 400 10 21 000671382 °00671382 ©00671382 ©00671382 4e0 15 21 0250 RAY 1 02099 RAY 2 03049 RAY 3 S| XXXXXXXXXXKXXXXXKXXXXXXXKXKXXXXXXKXKXXXKXKX OUTPUT XXXXXXXXXXXKXKXKXKXKXKXKXKXKXKKK KKK KKK KKK KK 16205 GRID OFF VAeCAPESs LINEAR INTERPe.98/16/64e PERIOD = 4¢0 SECes BATCH NOw 1s RAY NOe le MAX X Y ANGLE PCTNIF ] 14650 2050 120-00 2 14025 2093 120207 el 3} 14200 3037 120614 00 4 13275 3080 120022 020 5 13250 4023 12028 el 6 13024 4066 120.23 el U 12.99 5099 120¢51 3 8 12¢74 Dei2 12080 °3 9 1248 5295 121.210 03 10 12022 6028 121228 0-0 11 11-96 608) 121048 02 2 11270 Te23 121.270 020 13 11643 7266 121283 0.0 14 Wiloahy 808 121294 ol WS) 19290 8251 122019 02 16 10.64 8293 122059 e2 Wy NOSS7 9035 12288 el 18 19.09 GeTl 123295 el 19 9281 10.18 125027 404 20 9051 10.58 129256 44 21 9-18 10095 133285 4o4 22 8.75 11222 162047 3023 RAY STOPPED. 1620s GRID OFF VAeCAPESs LINEAR INTERPe 98/16/64e PERIOD = 4e0 SECes BATCH NOe ls RAY NOco 26 MAX x Y ANGLE PCTDIF 1 1536 2099 120200 2 MB) aba 3042 120-00 000 3 146286 3086 120-00 0.0 4 14.61 4079 120-00 0-0 5 14036 4072 120200 000 6 14.11 52016 120200 000 U 13286 5059 120-90 0e0 8 13661 6092 12000 020 9 136026 6045 120-00 0560 — 10 13611 6089 120200 00 11 12286 7032 120.02 000 2. 1261 T0775 120-06 0-0 13 1236 8019 120e11 020 52 14 12611 8e62 120617 020 15 11.286 9205 120224 020 16 11.60 9248 120230 0e0 W7/ 1135 9.91 120237 0e0 18 11.10 1034 120044 000 19 10.85 1078 12051 el 20 10.59 11.21 120057 000 2 10034 110264 12065 020 22 10208 12207 120671 el 23 9.83 12250 120.88 ol 24 9057 1292 121019 el 25 9e31 13035 121641 el 26 9205 12278 121055 ol 27 8.78 1420 121-98 05 28 8052 14063 122672 05 29 8024 15604 123027 04 30 7091 15042 139233 2961 3)il 7046 15263 W/o als} 29] RAY STOPPEDe 1620s GRID OFF VAeCAPESs LINEAR INTERPe 98/16/64e PERIOD = 4ce0 SECes BATCH NOs 1s RAY NOc 36 MAX x Y ANGLE PCTDIF 1 16022 3049 120.200 Q 15297 2092 120200 000 3 15072 4236 120200 020 4 15047 419 120200 0-0 5 15022 5022 120.00 020 6 14-97 5266 120200 000 v 14.72 6009 120-00 000 8 14047 6052 120200 000 9 14622 6095 120200 0-0 10 13097 72039 120200 0-0 i 13672 7082 120-00 020 UZ 13047 8225 120200 0-0 LES 13622 8-69 120-00 000 14 12697 9012 120-03 00 15 12272 9255 12008 000 16 12047 9298 120014 0-0 17 12022 10642 120019 020 18 11.97 10285 120024 00 19 1le71l 1128 12035 020 20 11046 11.71 120651 000 21 11621 12014 12062 000 22 10095 12057 120.68 el 23 10.70 13-00 120671 el 24 10044 13243 12072 ol 53 25) 10019 13086 26 9293 1429 2at 9067 14672 28 9e41 15015 29 9216 15258 30 8290 16200 Bil 8.63 16043 32 837 16086 3)3) Bell 1728 34 7284 17270 35 7056 18.11 36 7026 18.252 37 6091 18.87 RAY STOPPED. 12073 120-81 120-96 121212 12129 12144 121-56 121268 121.79 123651 125042 125092 144.59 54 APPENDIX E Computer Program for Wave Refraction (Linear-Interpolation) Using the IBM 7094 PROGRAM TITLE: MAIN 7094. Summary of Program: This program is essentially the same as that of MAIN 1620; both programs give identical results if processing the same input data. The chief differ- ences between the programs are as follows: (1) MAIN 7094 is written in FORTRAN IV, while MAIN 1620 is written in FORTRAN IT. (2) NOTT (the maximum number of batches to be run) is specified for MAIN 7091. (3) Sense Switches are not used in MAIN 709}. (4) ROUND subroutine is not used in MAIN 7094 since the 7094 computer auto- matically rounds off output data. (5) Due to the large storage capacity of the 7094, dimensions for CMAT in MAIN 7094 may greatly exceed those in MAIN 1620. (6) Since the 7094 computer outputs to a printer, MAIN 7094 has been pro- grammed to print title information at the top of each output page. (7) MAIN 7094 operates approximately 200 times faster than MAIN 1620. Remarks: MAIN 1620 is best suited for experimental work connected with further pro- gram refinements and development, while MAIN 7094 is best suited for pro- cessing large numbers of rays. 55 XXXMXKXKXKXKX KKK KX XXKXKKKXKXKXKKXXKXXX SOURCE PROGRAM XXXXXKXXKXKK KKK KM KKK HK KKK KKK KKK KK $IBFTC MAIN MAIN O1 G 70945 FORTRAN IVs LINEAR-INTERPOLATION WAVE REFRACTION PROGRAMe MAIN 02 C ORIGINAL PROGRAM BY GRISWOLD»sNAGLE»AND MEHRe THIS PROGRAM ADAPTED MAIN 03 Cc FROM ORIGINAL BY WILSON»HARRISON»KRUMBEINS AND BENSONe 8/4/640 MAIN 04 DIMENSION S(393)9EM( 493) 9E(3) sYVW(3) »CMAT(100982) 9C(4) »XLABL(12) MAIN 05 COMMON 5 » EM o — » YVW » CMAT iG MAIN 06 COMMON XLABL » D 5 Tt » CXY OT » NGO MAIN O07 COMMON AMM » ANN » MAX » NOT » N »9 MIT MAIN 08 READ (595) ((S(Is9J)9J=193)s9l=193) MAIN 09 5 FORMAT(6F1268) MAIN 10 READ (597) ((EM(LsI)9L=194) sl=193) MAIN 11 7 FORMAT(12F602) MAIN 12 READ (59500) NOTT MAIN 13 500 FORMAT (13) MAIN 14 DO 399 NOT=1leNOTT MAIN 15 READ (59400) XLABL MAIN 16 400 FORMAT (12A6) MAIN 17 READ (59402) MMsNN»CHECK »TTsNOJsD MAIN 18 402 FORMAT (2149F30007XsF5elsI59F4ol) MAIN 19 AMM = MM=1 MAIN 20 ANN = NN-1 MAIN 21 ITF (CHECK) 10910920 MAIN 22 10 READ (5911) ((CMAT(I9J)9T=19MM) 9J=19NN) MAIN 23 11 FORMAT (5(F100¢8s3X)) MAIN 24 | GO TO 14 MAIN 25 AG) Je i MAIN 26 READ (5911) (CMAT(19J)9I=19MM) MAIN 27 DO 77 J=29NN MAIN 28 DO 77 I=19MM MAIN 29 77 CMAT(T9J) = CMAT(I91) MAIN 30 14 DO 15 N=l»s NOJ MAIN 31 READ (596) AsXsY MAIN 32 6 FORMAT (FT70e292F6e02) MAIN 33 MAX = 1 MAIN 34 WRITE (69403) XLABL»oTTsNOToNoMAX 9X 9A MAIN 35 403 FORMAT (1H191246/9H PERIOD =9F5016H SECess1l0H BATCH NOcosl399H» RAMAIN 36 LY NOes1391He//4X 9 3HMAX 9 6X 9 LHX 98X 9 LHY 9 8X 9 SHANGLE 94X 96HPCTDIF// MAIN 37 21792F9e29Flle2) MAIN 38 A=A% 00174532925 MAIN 39 CALL RAYN (X9YsA) MAIN 40 15 CONTINUE MAIN 41 399 CONTINUE MAIN 42 WRITE (699999) MAIN 43 9999 FORMAT (18H1 THIS IS THE ENDe) MAIN 44 CALL EXIT MAIN 45 STOP MAIN 46 END MAIN 47 56 C SIBFTC RAYN G SUBROUTINE RAYN (X9YosA) DIMENSION S(393)9EM( 493) 9E(3) sYVW(3) sCMAT( 100982) 9C(4) sXLABL(12) COMMON S 9 EM 9 |e » YVW 9 CMAT COMMON XLABL » D » TT » CXY OH) 0 COMMON AMM » ANN » MAX » NOT » N CALL SURFCE (X9YsAsFK) MAX=1+MAX IF (MAX = 800) 399915915 CALL MOVE (XeYsAsFK) IF (CXY) 159159396 GO. TO (3979395915)s5 MIT PUNCH 2005 MAX FORMAT(33H CURVATURE APPROXIMATED FOR MAX =914) ANGLE=A*57029577951 IF (MOD(MAXs40)) 2025920 WRITE (697) XLABLsTT»NOT9N FORMAT (1H1912A6/9H PERIOD =9F5e1l96H SECess10H BATCH NOes1399H»s 1Y NOesI391He//4X 9 3HMAX 9 6X 9 LHX 9 8X 9 LHY 98X 9 SHANGLE 9 4X s6HPCTDIF//) PCTDIF = ABSF((C(3)-E(1)-E(2)-E(3))/C(3))*1006 WRITE (6912) MAX sX»9Y ANGLE sPCTDIF FORMAT (1792F9e2sFlle2sFl10el) GO TO (3915) »NGO WRITE (6913) FORMAT (13H RAY STOPPEDe) RETURN END S$IBFTC SURFCE PN SUBROUTINE SURFCE (XsYsAsFK) DIMENSION S$(333)59EM( 493) 9E (3) sYVW(3) »CMAT(100982) 9C(4) »XLABL(12) COMMON §S » EM 9 E » YVW » CMAT COMMON XLABL » D 9 TT » CXY » IT COMMON AMM » ANN » MAX » NOT » N I=X+1l. J=Y+1le FI=I FJ=J XL=X+1le-FI YL=Y+1le-FJ IF (MAX=-1) lslo4 Ife (Zur) oeow IF (ZJ-FJ) 19391 Zu S [py ZJ = FJ C(1)=CMAT(T J) C€(2)=CMAT(I+19J) C(3)=CMAT(I+19J+1) C(4)=CMAT( I 9J+1) 57 9 9 9 9 > 9 C NGO MIT C NGO MIT RAYN Ol RAYN 02 RAYN 03 RAYN 04 RAYN O05 RAYN O06 RAYN O7 RAYN 08 RAYN O09 RAYN 10 RAYN 11 RAYN 12 RAYN 13 RAYN 14 RAYN 15 RAYN 16 RAYN 17 RARAYN 18 RAYN 19 RAYN 20 RAYN 21 RAYN 22 RAYN 23 RAYN 24 RAYN 25 RAYN 26 RAYN 27 SURFCEO1 SURFCEO2 SURF CE03 SURFCE0O4 SURFCEO5 SURFCE06 SURFCEO7 SURF CE0O8 SURFCE09 SURFCE10 SURFCE11 SURFCE12 SURFCE13 SURFCE14 SURFCE15 SURFCE16 SURFCE17 SURFCE18 SURFCE19 SURFCE20 SURFCE21 318 DO 318 II=193 YVW(TT) = Oo DO 318 L=194 YVWCIT) = YVWCTI)+C(L)*EM(LoTT) DO 319 If=1s3 E(II) = Oc DO 319 JJ=193 E(TI) = ECIII+S(TIsJI)*#YVW( JJ) Oar 2 esl} se ECLQIIe se (EUS) AME FK = (E(2)*SIN(A)-E(3)*®COS(A))/CXKY RETURN END S$IBFTC MOVE 102 104 39 an Ww 38 SUBROUTINE MOVE (XsYoAsFK) DIMENSION S(393) 9EM( 493) 9E(3) sYVW(3) sCMAT( 100082) 9C(4&) oXLABL(12) COMMON S » EM 90 f » YVW » CMAT mG COMMON XLABL » D 5 Wi 9 CXY > IT » NGO COMMON AMM » ANN » MAX 9 NOT oN » MIT IF (MAX = 2) 10291025104 FKBAR=FK MIT = 1 DO 20 IT=1920 DELA=FKBAR*®D AA=A4+DELA ABAR=A+e5*DELA DELX = D*COS(ABAR) DELY = D*SIN(ABAR) XX=X+DELX YY=Y+DELY GO TO (10196) MIT CALL SURFCE (XX9YYsAAsFKK) IF (CXY) 38938910 FKBAR = 0¢5 * (FK + FKK) IF (IT — 18) 5293799 FKKPP = FKBAR IF (MAX = 2) 79799 IF (IT = 1) 2092009 TF (ABS (FKKP-FKBAR) -— (0000009/D)) 696920 FKKP = FKBAR IF (ABS (FKKPP = FKBAR) -— (02¢00009/D)) 18518917 MIT = 3 GO TO 38 FKBAR = O0e5 * (FKBAR + FKKP) MIT = 2 GO TO 39 NGO = 1 = IF ((XX—100) #( (AMM—100)—-XX) 129293 IF ((YY—-100)*( (ANN=100)-YY))29298 NGO = 2 X = XX V6 BAA A= AA FK = FKK RETURN END 58 SURFCE22 SURFCE23 SURFCE24 SURFCE25 SURFCE26 SURFCE27 SURFCE28 SURFCE29 SURF CE30 SURFCE31 SURFCE32 SURFCE33 MOVE 01 MOVE 02 MOVE 03 MOVE 04 MOVE 05 MOVE 06 MOVE 07 MOVE 08 MOVE 09 MOVE 10 MOVE 11 MOVE 12 MOVE 13 MOVE 14 MOVE 15 MOVE 16 MOVE 17 MOVE 18 MOVE 19 MOVE 20 MOVE 21 MOVE 22 MOVE 23 MOVE 24 MOVE 25 MOVE 26 MOVE 27 MOVE 28 MOVE 29 MOVE 30 MOVE 31 MOVE 32 MOVE 33 MOVE 34 MOVE 35 MOVE 36 MOVE 37 MOVE 38 MOVE 39 MOVE 40 MOVE 41 MOVE 42 MOVE 43 APPENDIX F Qu Derivations Relating QC with 2d, and 3C with 9 Given: C= tanh | eTld ” E ] Ox ayo MC Then: tanh} 2d} _ enc Cay - Ee efd _ tanh “1 fone Ge = 2 d= ORs =o [iia iil Ave) 6 dia jl > Aue TI gt gr Y K Let: Teo cy TOG, eimgl KY Ss Piuyeae Then: @ = Gk? fin (@ & RG) sc tm (Go x"c)] Given: Cla) 6 (a) andardil— 21x59) However,it may also be condidered that: = iG) gael © S GOK, 99) Therefore (Kaplan, 1952, p.86): 2d = d(d) . dC anddd _ d(a) . ac Tx 6a ox Ou) dot, os Sel — tlle re, - In(1-k"C)\9c +fin(1+k"C) - In(1-k"c)\oc Ox de Ox gs eG) a ie! ae y - (-k")\aC + fan(1+k"C) - In(1-k"c)\ac rp KC 1-k"CJQx ox 2 -9d 1 OX. “Soe ok! Similary: ac -gd. 1 iL QY oY &k' a 4 In(1+k"C) - In(1-k"C eC 59 APPENDIX G Method Used in Obtaining Most-Frequent Combinations of Deep-Water Height, Period, And Direction of Waves Capable of Striking Virginia Beach A rough estimate of the 15 most-frequent combinations of wave height, period, and direction of deep-water waves capable of striking Virginia Beach is gained by an analysis of five representative years of wave observations at the Chesapeake Lightship (fig. 1). Results of the analysis appear in table Gl and indicate that only six combinations of period and direction are involved in the 15 combinations. These six combinations are rough approxima- tions for input data, not only because of the crude methods involved in wave observation and recording, but also because it is assumed that one can pre- serve wave-direction angles for the waves observed at Chesapeake Lightship when searching for deep-water origin points that will yield rays capable of striking the target area. This is in error for some of the 4-second and all of the 6-second waves because they have already undergone refraction by the time they have reached Chesapeake Lightship; it is, therefore, invalid to select deep-water starting points for most waves by assuming that deep- water wave directions in water depths greater than 64 feet (Lightship water depth) are parallel to those observed at the Lightship. The 4 and 6-second waves that dominate the Lightship observations of table Gl are a product of the convention used here in interpreting the Ship's International Code. The convention is based upon an analysis of wave gage records that indicate a dominant 4-second period (table Ge) at a wave recording point off Cape Henry (fig. G1). (The Cape Henry wave gage was operated by the U. S. Navy between 1951 and 1956. It consisted 60 of a stepped-resistance wave staff; strip-chart recording were analyzed the significant wave method at 4-hour intervals.) Thus, code figure 2 Code Table 17 of the Ship's International Code, which stands for waves of 5-second period or less, was taken arbitrarily to represent waves of h4egecond period, while code figure 3 (5 to 7 seconds) was taken as 6 seconds. 61 Table Gl.-The fifteen most-frequently observed combinations of wave period, height, and direction (from 30°-150°T) as observed at Chesapeake Lightship during five representative years Combina- Period* Height Direction from No. of tion which wave front observa- (sec) (£t) travels (°T) tions** a 4 1 150 392 2 4 3 60 286 3 m 1 60 279 4 4 1 90 265 5 4 3 150 208 6 4 3 90 166 T 4 5 60 103 8 6 3 60 100 9 4 3 30 Qh 10 6 3 90 93 I. 4 1 30 We 12 6 5 60 12 13 6 al 90 50 14 iF 6 60 ha i) in 5 30 39 *Waves coded 2 {table 17, Ship's International Code) were considered as 4-sec waves; those coded 3 were considered 6-sec waves. See explanation in text and table G2. **TOtal observations | of waves capable of striking Virginia Beach (traveling from 30° through 150° True) = 3016 "08S €°G ST sucTZeATESqO FO Jos ageTduoo ayy soz spotued [Te jo ueoW oUTy,y SSCTO TePWy. —————— ee SNE ESE 6°9 +°S Gr Canis De Ditty Gr Oe LoS 62 Hot O°h G*g Qrg 0's ot G*tm oe G° Z°9 or€ Go H°? me G*g Ge ORS g°2 Gag G*z G°T Gul ord OH Gt 9°€ eT 6"t Die Sica JENS Gt ere Ge EPE GPE Sou, IES, EME q°? 6H g°e ECS Ow, Tene 6°2 gre OS GS. 1's G°9 ett Ord Er€ o'r 6°€ 62 G°9 TusiG OE cre Aor Tir 0°9 Sou ES) 9°9 weg 0°6 6°9 6°€ 0°9 q°9 Ort q°S o°h €°G L°S G°9 €°6 9°6 Lee 8°g 1°9 O's G°s 9°OL Me Mp q°S ae) ToUE oH ATE SAS SSE OTL G°8 SOE 6°S o°s 6°€I IE OR O°HL LST Gite OE) Ry CIE ate | OSE EOS Dire Gor AG MeL HOST E°OT EEE Fae I SL ME KOPN ear HIE Orge ayoSe Ory APN TZO0 SOE xSe en 4° QT Gack So) G°6 6°TT c°9OT 4°6 €°OL eve gre Tg €°G G°9 G°g 1° OL 9°OT Q°s 6°9 9°98 L°8 6°t 6°9 UG Org Cae HE 1°S Wes Eg et cy Tene a°€ LoS g°€ Gg et G*g 8° ge G* H* CoE g° Vea G*T earn tor Ie. o'2 ie ie Ee c° i’ Gore yO] (%) fouenberg (088) A a ard SR a ae ee ea AS RO ee oe DE ee eee oe ee Se EQEG GLE Ely Q6E €Th 92S €Ge 6Sh LLG HES PATA RSE GES (93 2°O H) spmooey on LOL oud AON LOO dus ONV Tar Nar AVN udv UVIN qa NVL >4quOW (sowry + pequeserdex ATaqetduoo yyuow yoes) saeak oat yequeserdet Mmof JO potted @ Fox aBeS anem Aruey odeg ay, 9e peATesqo spotsed anen ZOZ eTqeq Aouenbezz e@ JO uoTAAOg-"2y sTaqelL 63 CHESAPEAKE BAY BEACH RUDEE INLET *., LOCATION OF “+. WAVE GAGE INSTALLATION LOCATION OF AREA ce BATTERY Ba WORCESTER - ADMINISTRATION BUILDING 400 800 1200 1600 2000 FEE fe} 100 200 300 400 500 600 METERS FIGURE GI FORMER LOCATION OF U.S.NAVY WAVE GAGE OFF CAPE HENRY VIRGINIA, IN 20 FEET OF WATER 64 “SUOT}BUTWIAZ}9p 2UTJNOI IOF uoT}eTNITeD a},eINd9e pue ptdez sastwoid jnq a9e}s [Tej,uaudoTaAap ay} uT ST poyiow eyuL ‘“etursItA ‘ydeag eTUTSITA je aTOYS ay} 0} Ue|IDO IT {URTV 94} UT 1a}zeM daap woiz jyYysnorq aie sAeI JAM YITYM UT pajUasaid ST poyzow 294} JO atdwexa uy ‘suwerso1d rz3azndwod paads — ysty pue ‘SOT}STI9zIeIeYD JABM IajeMdaap yseoputTY IO paAiasqo Sutsn UOTLIVIFII 9AM FO UOT}LETNITYD IOF paqriosap st ainpaooid y STIFL III dadIdISSVIONN 9 “ON WNGNVYOWSW TVOINHOSL “S "M ‘UOSTIM IT “mM ‘uostizey I “saotpuadde g pue “p soinstz 6 ‘satqe} g Sutpnzout ‘dd 69 ‘po6T 31aq0190 sueiso1d rtayndwod “¢ UOsTTM “S‘"M pue uosTIzeH “M Aq NOILOVUSAY AAVM JO eTUTSITA NOILVINOTVO IVOI YAWNN YOd GOHLAW V SO INAWdOTSA ga “yoeag eTUTSITA °2¢ UOT}IVIFIY 2PACM “T “O° ‘“HSVM “HO ‘YAINSO “SHU “OYONT IVLSVOO AWUV’S‘A “SUOTJBUTWI9}9p JUT4{NOT TOF uoT}eTNITeD 92}eIND0e pue prtdei sastwoid jnq a3ejs [ejuawdoTaAap ay}, ut st poyiaw PyL “eTUTSITA ‘Yoeeg eTUTSITA 4e aTOYS 94} 0 UID IT JURTIV ay} ut Ia}eM daap worz yYysno1q are skeI aAeM YITYM UT pajuasaid ST poyjow a4} JO atTduexa uy ‘“sweiso1d 13}ndwod pasds — ysty pue “soT}StiazIe1eYyD aAeM I39}eMdaap 4seopuTY IO paAiasqo Sutsn UOT}IVIFOI QAVM FO UOT}RTNITeD IOF paqtio2sap st ainpasoid y eTITL III ddaIdyISSVDNN 9 “ON WNGNVYOWSN TVOINHOEL “S °M ‘UOSTTM IT “M ‘uostiiey I “saotpuodde g pue “p sain3tz 6 ‘setqe} @ Sutpnzout ‘dd pg ‘po6T 19qG0}90 suweiso0id rtayndwog ‘“¢ UOSTTM “S “M puke uostizeH *M Aq NOILOVUNAY AJAVM JO BTUTSITA NOILVINOTVO IVOINSWNN YO GOHLAW Y JO LNAWdOTSA gC “yoeag erursItA “z UoTzIeIFAY 2AeM “T “O’d ‘"HSWM ‘HO ‘UHINSO “SHU “DYONT IVLSvOO AWuv’s‘n “suotj}euTWIa}9ap aut NOI JOF UOT}ETNITeD 9}eIND90e pue ptdez sastwoid j4nq 98e}s [Te},uewdoTaAap ay} ut st poyjow eyuL “erursitA ‘yoeag etursitA je arOYS ay} 0} Uk|9DQ IT\URTIY 24} UT 1a}eM daep worzZ yYysnoIq are SABI 9AeM YOTYM UT pajUuasamd ST poyjow 94} Jo atdwexa uy ‘sweis0id 133ndwod paads — yaty pue ‘soT}STIa}IeIeYD 2AeM Iaj,eMdaap yseopuTY IO paAirasqo 3utsn UOTLIVIFII BACEM JO UOT}ETNITeD IOF paqtiosep st ainpadoid y eTITL III “S$ "M ‘UOSTIM II “mM ‘uostizey I “Y sweriso0id reyndwog ‘¢ BIUTSITA “yoeag eTUTSIIA “2 UOTJIVIFIY DACM “T aad IdISSVIONA 9 “ON WNGNVYOWSW IVOINHOAL “sodtpuadde g pue sainstzZ 6 ‘satqe} ¢g Butrpntout ‘dd pg “p96T 129q0;}90 UOSTTM “S "M pue uostizeH “mM Aq NOTLOVUAY JAVM SO NOILVINOTIVO TVOIYSWNN YOd GOHLYW V JO LNAWdOTSA ad “O° ‘*HSVM “dO ‘USINHO “SHU “OUONA IVLSVOO AWUV’S‘A “SUOTJEUTWIA}9p VUT}NOI IOF uOoT}eTNITeO |a}eINd0~e pue ptdez sastword 3nq a8e3s TeJuawdoTaAep ay} UT ST poyjaMw eyuL “erursaitA ‘yoeog etuTsItA je sTOYS dy} O} Ue|ADO ITJURTIV 24} ut 19}3eM daap worIzJ jyYysnoiq are skeI 3AeM YOTYM UT paj}Uaseaid ST poyjew a4}, FO aTdwexa uy ‘“sweid01d rta}ndwod psads — ysty pue ‘SOT}STIa}yOeIeYD 2AeM IayeMdaep yseopuTY IO paATasqo 3utsn UOT}IVIFAI 9AM FO UOT}RTNITeD IOF paqti9sep st sinpadsoid y eTItL III “S "M ‘UOSTIM II “mM ‘uostTiiey I “Y sweiso1d rajndwod “¢ eTUTSITA “yoeog eTUTSITA “2 UOT}LIVIFIY QACM “T dad IdISSVIONN 9 “ON WNGNVYOWSW TVOINHOSL “saotpuedde g pue Sainsty 6 ‘satTqe} zg Sutpntout ‘dd 69 ‘po6T 19qG03590 UOSTTM “S “M pue uostiieH “mM Aq NOILOVUdAY FAVM AO NOILVINOTVO IVOIYAWAN YOe GOHLEW V AO LNAWdOTSA AG O° ‘°HSWM “dO ‘YHINSO “Sdu “OYONT IVLSVOD AWUY'S'A “SUOT}BUTWII{9P 9UTJNOI IOF uot}eTNOTe. a}zeINdde pue ptdei sastwoird jnq 98e}s [ejuswdoTaAap 94} UT St poyiew eyL “etuTsItA ‘yoeog eTuTSITA 4e azT0YS ay} 0} ue|dD0 IT4URTIY ey} UT 1a}eM daap woIz jysnorq are sAeI JAeM YOTYM UT pajuasaid ST poyyew 34} FO atdwexa uy ‘“swerso1d 1azndwos paads — y3ty pue ‘SOT}STIazIeIeYD 2AM I3a,eMdaap jseopuTY IO paAiasqo Sutsn UOTIIVIFII 9ACM FO UOT}LETNITYI IOF paqriosap st ainpaosoid y eTtL III @HIAISSYIONA 9 “ON WNGNVYOWSW IVOINHOSL “SM ‘UOSTIM IT “mM ‘uostizey I “sadtpuadde g pue “p seinstz 6 ‘Satqe} ¢g Sutpntout ‘dd 9 ‘post 3129q0190 swerso0id tayndwog ‘¢ UuOSTTM “S*°M pue® uostizeyH “mM AQ NOTLOVUIATY AAVM JO eTUTSITA NOITLVINOTVO TVOI YAWNN YO GOHLIW Y JO LNAWdOTSAgC “yoeeg eTUTSITA ‘2 UOT}PIVIFIY SAPM “T “O°d ‘“HSVM “SO ‘USINAO “SHY “OYONA IVLSVOO AWUV'S’A “SUOTJBUTWI9Z9p JUTJNOI TOF uoT}eTNITeD a3eInNd0e pue ptdez sastwoid 3nq 93ejs [Te}UawdoTaAap 94} UT ST poyiaM eyuL “eSTUTSITA ‘Yyoeag eTUTSITA 1e aTOYS ay} 0} Uv|ADO IT iURTIYV 24} UT 19}emM daap worz jyYysnorq are SABI SAM YOTYM UT pajuasaid ST poyjew 94, FO aTdwexa uy ‘“sweis0id 1aj,ndwod pasds — ysty pue ‘soTystrazoerey aAeM 1a},eMdaap }seopuTY IO paAzrasqo Sutsn UOT}LIVIFOI 9AM FO UOT}RTNITeYD IOJ paqti9sap st aInpadoid y eTITL III dad1IdISSVINN 9 “ON WACNVYOWSW TIVOINHOL “SM ‘UOSTTM II “M ‘uostTiz1ey I “saotpuedde g pue “p sain3tz 6 ‘setqe} Z@ Sutpnpout “dd pg ‘po6T 19qG0}90 suerdozd rajndwop “¢ UOsTTM “S “M pue uosTzzeH “mM Aq NOILOVWHAY JAWM JO eTUTSITA NOILVINOITVO TVOINAWAN YOd GOHLEW V JO INAWdOTaAga “yoeag eTursItA °Z UOT}IVIFIY SACM “T “O° ‘*HSVM ‘HO ‘USINEO “SHU “OYONA IVLSVOO AWUV'S‘A “SUOT}eUTWIA}9p JUTJNOI IOF UOT}eTNITeD 9}eIND0e pue ptdez sastword 4nq a8ejs ~TejuewdotTaAap ay} ur st poyjzow ayL “eSTUTSITA ‘ydeag BTUTSITA }e 9IOYS 94} 0} UBZDQ IT\URTIY ay} UT 19}zeM daap woIF yYysnorq are sAkeI |BAeM YOTYM UT pajUasaid ST poyjow 24} FO aTdwexe uy ‘“sweis0id 1a,ndwod paads — ysry pue ‘soT}sti9}Ie1eYyD aAeM Iaj,emMdaap yseoputYy Io paAzasqo Sursn UOT}IVIFII BACM JO UOT}ETNITeD IOF paqtissep st ainpaooid y Cees eH dadIdISSVIONN 9 “ON WNGNVYOWSW TVOINHOAL “S "M ‘UOSTIM II “mM ‘uostTizey I “sodtpuadde g pue “p sainstz 6 ‘satqe, ~@ Butpntout “dd pg “p96T 19qQ0390 sweisoid rezndwop ‘¢ UOSTIM “S °M pue® uostizeH *M Aq NOTLOVYNAY AAVM JO eTUTSITA NOTLVINOIVO TVOIUXWNN YOd GOHLAW V JO LNAWdOTSA SG “yoeag BTUTSITA “2 UOTPIVIFIY 9ACM “T “O°d ‘"“HSVM “dO “UAINHO “SHU “SUONA IVLSVOO AWUV’S ‘A “SUOT}JEUTWIA}9pP JUT}NOI TOF uoT}eTNITed a4yeind9e pue ptdez sastwoid 3nq 93e}4s [Tej,uawdoTaAap ay}, ut st poyiaw ayL “eTUTSITA ‘ydeag eTUTSITA 1e 2TOYS ay} 0} UkADO ITRURTLY 24} UT 19,eM daap woIy }yYysnoI1q are SABI SAeM YITYM UT pajUuasaid ST poOyjJew 34} Fo atTdwexs uy ‘sweisoid 1tayzndwoos psads — ysty » pue ‘SOT}STI9}DeIeYD SAeM I3a}eMdaap ySeopUTY IO paAZasqo 3utsn UOTJIVIFII BACM FO UOT}JETNITeD IOF paqridsap st sinpeadsoid y eT4TL III a ° ‘ = ON WHW IVOINHOSL “So °M ‘uos—tt™M II GHIAISSVIONA 9 WNdNVuol “mM ‘uostiiey Te “saotpuedde g pue “py Sainsty 6 ‘satqe} zg Butpnqout ‘dd ~9 ‘poet 12q0790 sweisoid rajndwop “¢ UOSTTM °S “M pue uostTizeH “M 4q NOTLOVYSAY FAVM JO eTUTSITA NOILVINOTVO IVOIYSWNN YOed GOHLAW VY JO LNAWdOTSA ad “yoeeg eTUTSITA “Z UOT}IVIFOY SACM “T “O°d ‘“HSVM “HO ‘USINTO “SHY “SYONA IVLSVOO AWUV'’S"A aT x ie i i. ae % i A oF