                  COUPLE, November 20, 1997 Version 
         November 20, 1997 Documentation by Richard B. Evans 
         
               Woods Hole Oceanographic Institution
         404 Bigelow Bldg., MS #9, Woods Hole, MA 02543, USA
                          (508) 289-3513

                                and

             Science Applications International Corporation
         21 Montauk Avenue, Suite 201, New London, CT 06320, USA
                          (860) 442-0564

     INTRODUCTION:

       COUPLE is a program that performs a stepwise coupled mode [1] 
     calculation of the two dimensional underwater acoustic field due a 
     vertical array of harmonic point or line sources in cylindrical or 
     plane geometry respectively. The environment may be range dependent, 
     consisting of a series of horizontally stratified fluid filled 
     regions. Each of the regions consists of a water layer and a 
     sediment or bottom layer that are assumed to have a piecewise linear 
     depth variation in the index of refraction squared and a piecewise 
     exponential depth variation in the density. The index of refracton 
     squared and density are continuous in the water layer and in sediment 
     layer, but they may be discontinuous at the water sediment interface. 
     The index of refraction squared may be complex in both the water and 
     bottom sediment layers to allow for attenuation.

       The modal calculation in each of the regions is accomplished with 
     a search or Galerkin procedure. If the problem consist of a 
     homogeneous water layer over a homogeneous bottom sediment layer 
     and only the water depth varies with range, then the complex 
     eigenvalues can be found using a search procedure based on Newton's 
     method. 

       In the more general case with inhomogeneous water and bottom 
     sediment layers, range dependent material properties and range 
     dependent bathymetry, the Galerkin matrix method is used. The 
     Galerkin solution is constructed on a basis set from a simple 
     background model consisting of a homogeneous water layer over a 
     homogeneous bottom sediment layer, which retains the density 
     discontinuity at the water sediment interface.

       The Galerkin procedure, with a finite basis set, yields a complex 
     matrix eigenvalue problem. The complex eigenvalues are found using 
     a Jacobi type algorithm designed for complex symmetric matrices. 
     This routine does not order the eigenvalues by increaing real part.
     A description of the Galerkin method for the constant density water 
     layer and homogeneous sediment layer case can be found in [2]. When 
     density variations are present, the eigenvalue problem becomes a 
     generalized eigenvalue problem and the Cholesky factorizaton is 
     applied before the Jacobi type algorithm. 

       The range dependence is described by a series of ranges and water
     depths along with the associated sound speed, density and attenuation
     profiles. A pressure release boundary is imposed at the surface of
     the water and at the greatest depth in the bottom sediment. Although
     the water depth varies with range, the total depth of the water and 
     sediment is fixed. 

       A decoupling algorithm [3] is used to solve the matrix two-point
     boundary value problem in range that occurs in the stepwise coupled 
     mode formulation. This eliminates the numerical instability that 
     occurs with other methods such as shooting.

       The purpose of this document is to describe the inputs and outputs 
     of COUPLE. There is also a brief description of how to compile, link 
     and run the FORTRAN source code that comprises the program.



     I. INPUTS:

     Block 0 (1 record): TITLE
     -------
     Block 1 (1 record): HB,FREQ,ROHW,ROHB,DBPWL
     -------
     HB = Total thickness of water and sediment layer (m).
     FREQ = Frequency (Hz).
     ROHW = Density of the water layer (g/cm**3).
     ROHB = Density of the sediment layer (g/cm**3).
     DBPWL = Attenuation in the sediment layer (dB/wavelength).
     
     Block 2 (1 record): CB,CW
     -------
     CB = Sound speed in the sediment (m/s).
     CW = Sound speed in the water (m/s).
     
     Block 3 (1 record):  M,NATEN,IEIG,NDPTH,IGEOM,IREFL
     -------
     M = Number of contributing modes (Maximum=MX=190).
     NATEN = Number of incremental parts of DBPWL to use in the search for
             eigenvalues.
     IEIG = Flag used to indicate the source of the estimates of complex
            eigenvalues in region 2 through N in the search mode.
               IEIG=0 - The real eigenvalue problem is used and the 
               attenuation is introduced, based on NATEN.
               IEIG=1 - The complex eigenvalues of the previous region 
               are used to start.
     NDPTH = Number of incremental parts of the depth step to use with 
             IEIG=1 in the search for complex eigenvalues.
     IGEOM = Flag used to indicate cylindrical or plane geometry.
                IGEOM=0 - Cylindrical geometry is used with point sources.
                IGEOM=1 - Plane geometry is used with line sources.
     IREFL = Flag to control the source condition at range zero.
                IREFL=0 - Perfectly cylindrically symmetric geometry is 
                          assumed.
                IREFL=1 - The source is backed by an absorber at negative
                           ranges in cylindrical or plane geometry.

     Block 4 (1 record): NSDEP,ZSMIN,ZSINC
     -------
     NSDEP = Number of source depths (Maximum=MD=400).
     ZSMIN = Minimum source depth (m).
     ZSINC = Increment in source depth (m)

     Block 5 (1 record): NDEP,ZMIN,ZINC
     -------
     NDEP= Number of receiver depths in the pressure calculation 
           (Maximum=MD=400).
     ZMIN = Minimum receiver depth (m).
     ZINC = Increment between the receiver depths (m).

     Block 6 (1 record): RMIN,RMAX,RINC
     -------
     RMIN = Minimum range in the pressure calculation (km).
     RMAX = Maximum range in the pressure calculation (km).
     RINC = Increment between ranges in the pressure calculation (km).

     Block 7 (1 record): N,IPRT,IOUT,MAMP
     -------
     N = Number of environmental regions, headed by range depth pairs, to 
         follow (Maximum=NRX=2002).
     IPRT = Print flag for details of the calculation such as eigenvalues,
            coupling matrices and coefficients of outgoing and ingoing 
            waves.
     IOUT = Flag that controls the matching condition used at the vertical
            interfaces in the range marching scheme.
               IOUT=0 - The program computes the fully elliptic two way 
                        solution that matches both the pressure and radial 
                        particle velocity.
               IOUT=1 - The one way pressure matching solution is found.

               IOUT=-1 - A single scatter approximation is used.
          MAMP = Flag to control generation of a file of mode amplitudes. 

     ********  The following input block is repeated for J=1,N   ********

     Block 8(J) (1+NPW(J)+NPB(J) records):
     ----------

     First record: RANG(J),IRLIN(J),DPTH(J),NPW(J),NPB(J),IREIG(J)

     RANG(J) = Range to the end of the J th region (km).
     IRLIN (J) = Number of additional regions generated in subroutine  
                 BATHY to represent a linearly sloping bottom between 
                 RANG(J) and RANG(J+1).
     DPTH (J) = Depth of the water in the J th region (m).
     NPW(J) = Number of points in the water sound speed, attenuations, 
              density profile in the J th region.
             NPW(J)=0 - The water is assumed to be homogeneous with sound 
                        speed CW and the search procedure is used.
             NPW(J)>1 - The water column may be nonhomogeneous, a sound 
                        speed profile must be provided and the Galerkin 
                        procedure is used.
     NPB(J) = Number of points in the bottom sound speed, attenuation, 
              density profile in the J th region. NPB(J) must follow the 
              same two cases as NPW(J). 
              NPB(J)=0 - The bottom is homogeneous with sound speed CB, 
                         density RHOB and attenuation DBPWL.
              NPB(J)>1 - The bottom may be nonhomogeneous and a bottom 
                         sound speed, attenuation, density  profile must 
                         be provided.
     IREIG (J) = Sequence number for storage or retrieval of eigenvalues 
                 on a direct access file for use with perodic media only.
                      IREIG(J)=0 - Omit storage on direct access file
                      IREIG(J)>0 - Write eigenvalues on direct access file
                                   at record IREIG(J) 
                      IREIG(J)<0 - Read eigenvalues from direct access 
                                   file at record ABS(IREIG(J))
                              
     Next NPW(J) records: DEPW(J,K),SVPW(J,K),DBPWLW(J,K),RHOW(J,K), 
                          K=1,NPW(J) 

     DEPW(J,K) = Depth (m) of sound speed in the water (DEP(J,1)=0.0).
     SVPW(J,K) = Sound speed in the water (m/s).
     DBPWLW(J,K) = Attenuation in the water (dB/wavelength).
     RHOW(J,K) = Density in the water (g/cm**3).

     Last NPB(J) records: DEPB(J,L),SVPB(J,L),DBPWLB(J,L),RHOB(J,L), 
                          L=1,NPB(J)

     DEPB(J,L) = Depth (m) of sound speed in the bottom measured from 
                the water surface (DEPB(J,1)=DEP(J,NPW(J))).
     SVPB(J,L) = Sound speed in the bottom (m/s).
     DBPWLB(J,L) = Attenuation in the bottom (dB/wavelength).
     RHOB(J,L) = Density in the bottom (g/cm**3).


     II. SUGGESTIONS AND COMMENTS:

     Block 1:
     -------
     HB should be large enough to provide space for an absorber in the 
     deepest part of the problem. The artificial absorber, where the 
     attenuation increases significantly with depth, is often included 
     at the bottom of the physical sediment.
     ROHW, ROHB and DBPWL are used in search mode and ignored in the 
     Galerkin mode.

     Block 2: 
     -------
     CB is always nonzero. It is used to compute the index of refraction 
     everywhere. CB is the bottom sound speed in the search mode and it 
     is the bottom sound speed for the background two layer problem in 
     every region in the Galerkin procedure. CW is used in the search 
     mode and is ignored in the Galerkin mode.

     Block 3:
     -------
     M should be chosen to be significantly larger than the number of 
     propagating modes to represent high angles and to help in the 
     expansion of the modes in the Galerkin procedure. Convergence 
     should be checked by increasing M until the change in the 
     calculated field is insignificant.

     NATEN, IEIG and NDPTH are only used in the search mode. IEIG=0 
     is always assumed in the first region. If IEIG=1, then NDPTH is 
     used when redundant eigenvalues are found, which means the one 
     was missed (See subroutine SEARCH). NATEN, IEIG and NDPTH can be 
     set to zero, in the Galerkin mode, to avoid confusion.

     IREFL=0 yields the source condition described in [1 and 2] for 
     perfectly cylindrically symmetric geometry. In this case, there 
     exists standing waves, in range, in the region containing the source. 
     This is a multiple scattering effect. In perfectly cylindrically 
     symmetric geometry, the alternative IREFL=1 is nonphysical, but it 
     is  consistent with the single scatter approximation. IREFL=1 is 
     always used in plane geometry.

     Block 4
     -------
     The sources may be in both the water layer and/or the bottom  
     sediment layer. If NSDEP>1, then a vertical array of in-phase 
     point or line sources projecting a horizontal beam is used. The 
     field is found by a superpostion of the fields due to the NSDEP 
     sources. The sum of sources is divided by NSDEP. See the main 
     program COUPLE.

     Block 6
     -------
     RMIN can be zero, with IOUT=0, in the calculation of the scattered 
     field that is regular at the source. The other fields are irregular 
     at range zero.
     RMAX can be larger than the end of the last environmental region 
     RNG(N).

     Block 7 
     -------
     If IPRT>0, then details of the calculation are printed for the first, 
     last and every IPRT region in between. 
     If IOUT=-1, an approximation to the single scatter approximation due 
     to Porter, Jensen and Ferla [4] is used. See subroutine PROPM for 
     the implementation. If MAMP>0, then the mode amplitudes are written 
     out in CPRES every MAMP range steps. This file can get very big so 
     pick MAMP carefully.

     Block 8(J) :
     ----------
     In long flat regions with IOUT=0, it may be necessary to include 
     zero height bathymetry changes to activate the decoupling algorithm 
     to normalize the solution and accurately compute ingoing and outgoing 
     waves.

     RNG(J) is the range to the end of the J th region so RNG(1)>0. 
     RNG(N) has no significance, since the last region extends 
     indefinitely.

     If IRLIN(J)>0 then interpolation and extrapolation procedures are 
     used. As the water depth increases (decreases) the sound speed, 
     attenuation, density profile in the water will be extrapolated 
     (interpolated). As the water depth increases (decreases),the bottom 
     sound speed, attenuation, density profile is shifted down (up) and 
     the interpolation (extrapolated) takes place in deepest part of 
     sediment profile at HB.

     NPW(J) controls the procedure used to find eigenvalues.

     If NPW(J)=0, then the water is assumed to be homogeneous with 
     sound speed CW. No depth sound speed pairs are input. The eigenvalues 
     are found in subroutine SEARCH based on NATEN, NDPTH and IEIG. The 
     analytic forms of the depth functions are used.

     If NPW(J)=2 or more, then the water column may be nonhomogeneous. 
     A sound speed, attenuation, density profile must be input. The 
     eigenvalues are found by the Galerkin method in subroutine GALRKN 
     based on the Jacobi type matrix eigenvalue algorithm. The depth 
     functions are constructed from the matrix eigenvectors.

     NPB(J) follows the same two cases as NPW(J). If NPW(J)=2 or more, 
     then NPB(J)=2 or more and a bottom sound speed, attenuation, density  
     profile must be provided. If NPW(J)=0 is then it is necesary that 
     NPB(J)=0 and no bottom profile may be input.

     IREIG(J) is used with perodic media to avoid recomputation of 
     identical sets of eigenvalues using subroutine SEARCH or GALRKN. 
     It is not valid when IRLIN(J) is nonzero because it is too hard 
     to anticipate the sequence numbers needed. It is not valid with a 
     mixture of calls to subroutines SEARCH and GALRKN to make it easy 
     to keep track of the records written on direct access file.

     RHOW(J,K) overrides the value of ROHW (Note the different spelling).

     DBPWLB(J,L) and RHOB(J,L) override the values of ROHB and DBPWL.


     III. EXAMPLES:

       The following four examples are chosen because the calculated 
     transmission loss appears in the references [1-5] and they may be 
     used to test the code. The Galerkin procedure is less prone to 
     problems, but input for the search procedure is given, in some 
     examples, for contrast. The results from both procedures should be 
     essentially the same.

     Example 1. Kutschale Problem [1].

     Kutschale Problem, Search procedure, 70 Modes
     5000.  25.0   1.0   1.15   0.5         !HB,FREQ,ROHW,ROHB,DBPWL
     1704.5    1500.                        !CB,CW
     70  20  1  0  0  0                     !M,NATEN,IEIG,NDPTH,IGEOM,IREFL
     1   25.0    0.0                        !NSDEP,ZSMIN,ZSINC
     1   25.0    0.0                        !NDEP,ZMIN,ZINC
     0.05   5.1   .05                       !RMIN,RMAX,RINC
     2   1   0   0                          !N,IPRT,IOUT,MAMP
     2.50     0    32.5   0  0   0          !RANGE,IRLIN,DPTH,NPW,NPB,IREIG
     5.10     0    32.5   0  0   0

     Kutschale Problem, Galerkin procedure, 25 Modes
     1000.  25.0   0.0   0.00   0.0         !HB,FREQ,ROHW,ROHB,DBPWL
     1704.5    00.00                        !CB,CW
     25  0  0  0  0  0                      !M,NATEN,IEIG,NDPTH,IGEOM,IREFL
     1   25.0    0.0                        !NSDEP,ZSMIN,ZSINC
     1   25.0    0.0                        !NDEP,ZMIN,ZINC
     0.05   5.1   .05                       !RMIN,RMAX,RINC
     2   1   0   0                          !N,IPRT,IOUT,MAMP
     2.50     0    32.5   2  3   0          !RANGE,IRLIN,DPTH,NP,NPB,IREIG
      0.0   1500.0   0.0    1.00            !DEPW,SVPW,DBPWLW,RHOW
     32.5   1500.0   0.0    1.00
     32.5   1704.5   0.5    1.15            !DEPB,SVPB,DBPWLB,RHOB
     250.0  1704.5   0.5    1.15
     1000.0 1704.5   2.5    1.15
     5.10     0    32.5   2  3   0
      0.0   1500.0   0.0    1.00            !DEPW,SVPW,DBPWLW,RHOW
     32.5   1500.0   0.0    1.00
     32.5   1704.5   0.5    1.15            !DEPB,SVPB,DBPWLB,RHOB
     250.0  1704.5   0.5    1.15
     1000.0 1704.5   2.5    1.15


     Example 2. Square wave 5-10 km patch of roughness [2 and 5]. The 
     published transmission loss plots are under sampled, at a range 
     increment of .05 km, in both references. The .01 km range increment 
     used here shows a good deal more structure.

     Patch of roughness, Square wave corrugation, Search procedure, 40 modes
     1500.  25.   1.0   2.5   0.5           !HB,FREQ,ROHW,ROHB,DBPWL
     1704.5 1500.0                          !CB,CW
     40  20  0  0  0  0                     !M,NATEN,IEIG,NDPTH,IGEOM,IREFL
     1  18.0   0.0                          !NSDEP,ZSMIN,ZSINC
     1  50.0   0.0                          !NDEP,ZMIN,ZINC
     0.00   18.1   .010                     !RMIN,RMAX,RINC
     101   10   0   0                       !N,IPRT,IOUT,MAMP
     5.00     0    100.   0  0   1          !RANGE,IRLIN,DPTH,NPW,NPB,IREIG
     5.05     0     90.   0  0   2       
     5.10     0    100.   0  0   -1
     5.15     0    90.    0  0   -2 

     (Repeat the last two records 47 times increasing both ranges by .1 km
     in each cycle. The -1,-2 sequence eliminates recomputation of 
     eigevalues.)

     9.90     0    100.   0  0   -1
     9.95     0    90.    0  0   -2 
     20.00     0    100.   0  0  -1             


     Example 3. PE Workshop II, Test Case 3b [3].

     PE Workshop II, Test Case 3b , Galerkin procedure, 28 modes
     400.  250.0   0.0   0.0   0.0          !HB,FREQ,ROHW,ROHB,DBPWL
     1590.0   00.00                         !CB,CW
     28  0  0  0  0  0                      !M,NATEN,IEIG,NDPTH,IGEOM,IREFL
     1   99.5    0.0                        !NSDEP,ZSMIN,ZSINC
     1   99.5    0.0                        !NDEP,ZMIN,ZINC
     4.9  10.1   .01                        !RMIN,RMAX,RINC
     2   1   0   0                          !N,IPRT,IOUT,MAMP
     5.00     0    100.0   2  2   0         !RANGE,IRLIN,DPTH,NPW,NPB,IREIG
       0.0    1500.0    0.0   1.0           !DEPW,SVPW,DBPWLW,RHOW 
     100.0    1500.0    0.0   1.0
     100.0    1590.0    0.5   1.2           !DEPB,SVPB,DBPWLB,RHOB 
     400.0    1590.0    0.5   1.2
     5.10     0    100.0   2  2   0
       0.0    1500.0    0.0   1.0           !DEPW,SVPW,DBPWLW,RHOW 
     100.0    1500.0    0.0   1.0
     100.0    1590.0    0.5   1.2           !DEPB,SVPB,DBPWLB,RHOB
     400.0    1590.0    0.5   1.2


     Example 4. Acoustical Society of America (ASA) Benchmark Wedge [4].

     ASA Benchmark Wedge, Galerkin procedure, 45 modes, 200 steps
     1500.  25.0   1.0   1.50   0.0          !HB,FREQ,ROHW,ROHB,DBPWL
     1700.0    00.                           !CB,CW
     45  0  0  0  0  0                       !M,NATEN,IEIG,NDPTH,IGEOM,IREFL
     1  100.0    0.0                         !NSDEP,ZSMIN,ZSINC
     2   30.0    120.0                       !NDEP,ZMIN,ZINC
     0.01   4.01   .01                       !RMIN,RMAX,RINC
     2   50   0   0                          !N,IPRT,IOUT,MAMP
     0.01     200   200.0   2  3   0         !RANGE,IRLIN,DPTH,NPW,NPB,IREIG
     0.0      1500.0   .0    1.00            !DEPW,SVPW,DBPWLW,RHOW
     200.0    1500.0   .0    1.00
     200.0    1700.0   .5    1.50            !DEPB,SVPB,DBPWLB,RHOB
     1000.0   1700.0   .5    1.50
     1500.0   1700.0  2.5    1.50
     4.000      0     0.01   2  3   0
     0.0      1500.0   .0    1.00            !DEPW,SVPW,DBPWLW,RHOW
     0.01     1500.0   .0    1.00
     0.01     1700.0   .5    1.50            !DEPB,SVPB,DBPWLB,RHOB
     1000.0   1700.0   .5    1.50
     1500.0   1700.0  2.5    1.50
     
     ASA Benchmark Wedge, Search procedure, 90 modes, 200 steps
     3000.  25.0   1.0   1.50   0.5          !HB,FREQ,ROHW,ROHB,DBPWL
     1700.0    1500.0                        !CB,CW
     90  30  1  20  0  0                     !M,NATEN,IEIG,NDPTH,IGEOM,IREFL
     1  100.0    0.0                         !NSDEP,ZSMIN,ZSINC
     2   30.0    120.0                       !NDEP,ZMIN,ZINC
     0.01   4.01   .01                       !RMIN,RMAX,RINC
     2   50   0   0                          !N,IPRT,IOUT,MAMP
     0.01     200   200.0   0  0   0         !RANGE,IRLIN,DPTH,NPW,NPB,IREIG
     4.000      0    0.01   0  0   0


     IV. INPUT AND OUTPUT FILES:

     COUPLE.DAT = Free format input data as described above.

     COUPLE.LOG  = Formatted log file (may appear on screen).
     COUPLE.PRT  = Formatted print file.
     COUPLE.TL   = Formatted transmission loss.
     COUPLEO.TL  = Formatted outgoing transmission loss if IOUT=0
     COUPLEI.TL  = Formatted ingoing transmission loss if IOUT=0
     COUPLES.TL  = Formatted scattered transmission loss in region 1 
                   if IOUT=0  

     COUPLE.CPR  = Unformatted real and imaginary parts of the complex 
                   pressure.
     COUPLEO.CPR = Unformatted outgoing complex pressure if IOUT=0
     COUPLEI.CPR = Unformatted ingoing complex pressure if IOUT=0
     COUPLES.CPR = Unformatted scattered complex pressure in region 1 
                   if IOUT=0  
     COUPLE.AMP  = Unformatted mode amplitudes when MAMP>0


     V. COMPILING, LINKING AND RUNNING:

       Windows NT or 95, Command window (DOS): There is a batch file  
     called "FALLSUBS.BAT" that compiles all the subroutines (70) using 
     the Lahey Fortran 90 compiler. The Fortran source code has names of 
     the form *.FOR. The compiler option 'vax' is used along with 
     targeting for your your processor. To compile the subroutines type 
     "FALLSUBS.BAT" at the DOS prompt in command window. There is a batch 
     file "FLCOUPLE.BAT" that compiles the main program COUPLE and links 
     it with the subroutines that are listed in the file "COUPLE.LNK" 
     read by the linker. To compile and link the main program COUPLE, 
     type "FLCOUPLE.BAT" at the DOS prompt in the command window.

       There is a batch file called "COUPLE.BAT" that imports the name 
     of the *.dat input file, runs the program and handles the output 
     file naming. To run COUPLE create an input file called "name.dat" 
     and type "COUPLE.BAT name" at the DOS prompt in the command widow.  
     When the program finishes, some or all of the output files listed 
     above are created with "name" in place of COUPLE and with the 
     various extensions. The string "name" should be no more than 7 
     characters for backward compatibiitly with DOS, where the maximun 
     is 8. The format of the output files is established in the 
     subroutines HEADER and CPRES.

     Unix: There is a unix make file called "makefile" that lists all 
     the subroutines. To compile and link type "make couple" at the unix 
     prompt. There is a executable script file called "couple.bat" that 
     runs the program and handles the input and output file naming. To run 
     couple perform the following steps. First create an input file, as 
     described above, with a unique name like "newrun.dat", Next use the 
     vi editor to replace the "oldrun" from the last run with the "newrun" 
     everywhere in couple.bat. This can be done with the vi command 
     "%s/oldrun/newrun/g" at the ":" prompt in vi. Exit vi. Lastly run 
     couple in the batch mode by typing the three lines "batch <cr>" and  
     "couple.bat <cr>" and "control d <cr>" at the unix prompt.
     
       It is necessary to open a file in the main program, say 
     COUPLE.LOG, to except the output written with WRITE(* that normally 
     comes to the screen in a DOS window. It is also be necessary to 
     change timming calls in the main program COUPLE and the date call 
     in HEADER, which are system dependent.


     REFERENCES:

     1. R. B. Evans, "A coupled mode solution for acoustic propagation 
     in a waveguide with stepwise depth variations of a penetrable 
     bottom," J. Acoust. Soc. Am. 74, 188-195 (1983).

     2. R. B. Evans and K. E. Gilbert, "Acoustic propagation in an 
     refracting ocean waveguide with an irregular interface," Comp. 
     Maths. Appl. 11, 795-805 (1985). 

     3. R. B. Evans, "The decoupling of stepwise coupled modes," J. 
     Acoust. Soc. Am. 80, 1414-1418 (1986).

     4. M. B. Porter, F. B. Jensen and C. M. Ferla, "The problem of 
     energy conservation on one way models," J. Acoust. Soc. Am. 89, 
     1058-1067 (1991).

     5. R. B. Evans and K. E. Gilbert, "The perodic extension of 
     stepwise coupled modes," J. Acoust. Soc. Am. 77, 983-988 (1985).

