* Program to simulate "standing" wave response to obstacle at
* bottom of steady-flowing stream of water in trough. 

* Formulas are based on Lighthill's "Waves in Fluids", chap. 3

* In this program the obstacle is a "cosine bump" of length L, i.e.,
* half of a single period of a cosine wave of wavelength 2*L,
* centered on x=0. The maximum height of the bump is a.
* The depth z of the bottom is therefore

*               |  -h                  if |x| >  L/2
*         z =   |
*               |  a*cos(pi*x/L)       if |x| <= L/2

      program stand1

      implicit none

      integer iv,ivlo,ivhi,niter,icol(50)
      double precision h,l,a,xmin,xmax,zmin,zmax,g,pi,vmax,
     +                 vtry(50),v,f,alpha,diff,k0,kh,kl,zeta,
     +                 x1,y1,y2,y3,y4,yoff,xv,yv,gamma,factor

* Define nominal depth of trough (in feet)

      h = 2.d0

* Define length of "cosine bump" (in feet)

      l = 0.5d0

* Define maximum height of "cosine bump" (in feet)

      a = 0.4d0

* Define distances upstream and downstream to plot

      xmin = -2.d0
      xmax = +6.d0

* Define minimum and maximum z values of plot

      zmin = - h * 1.2d0
      zmax = + h * 1.2d0

* Define constants

      g = 32.d0
      pi = acos(-1.d0)

* Find maximum allowed velocity for which standing waves
* can be created

      vmax = sqrt(g*h)

* Define the stream velocity values to try (as fractions of vmax)

      ivlo = 1
      ivhi = 5
      vtry(1) = 0.30d0
      vtry(2) = 0.50d0
      vtry(3) = 0.75d0
      vtry(4) = 0.90d0
      vtry(5) = 0.95d0
*      vtry(6) = 0.99d0

* Select colors for plotting

      icol(1) = 1
      icol(2) = 2
      icol(3) = 3
      icol(4) = 4
      icol(5) = 6
*      icol(6) = 6

* Initialize plot file

      x1 = xmax * 0.4d0
      y1 = zmax * 0.9d0
      yoff = zmax*0.15d0
      y2 = y1 - yoff
      y3 = y1 - 2.d0*yoff
      y4 = y1 - 3.d0*yoff
      open(1,file='stand1.kumac',status='new')
      write(1,5) xmin,xmax,zmin,zmax,
     +           pi,a,l,h,
     +           xmax,xmin,
     +           -l/2.d0,l/2.d0,
     +           xmin,-h,-l/2.d0,-h,
     +           l/2.d0,-h,xmax,-h,
     +           x1,y1,h,
     +           x1,y2,vmax,
     +           x1,y3,l,
     +           x1,y4,a
 5    format('fort/file 66 stand1.ps'/
     +       'meta 66 -111'/
     +       'zone 1 1'/
     +       'null',4(1x,f10.6)/
     +       'igset plci 1'/
     +       'igset txci 1'/
     +       'igset chhe 0.35'/
     +       'nxd = 100'/
     +       'nxu = 100'/
     +       'nxb = 50'/
     +       'pi = ',f10.8/
     +       'a = ',f10.6/
     +       'l = ',f10.6/
     +       'h = ',f10.6/
     +       'sigma xd = array([nxd],0#',f10.6,')'/
     +       'sigma xu = array([nxu],',f10.6,'#0)'/
     +       'sigma xb = array([nxb],',f10.6,'#',f10.6,')'/
     +       'sigma bump = -[h] + [a]*cos([pi]*xb/[l])'/
     +       'graph [nxb] xb bump l'/
     +       2('line',4(1x,f10.6)/)/
     +       'itx ',f10.6,1x,f10.6,' Depth = ',f4.2,' ft'/
     +       'itx ',f10.6,1x,f10.6,' Max V = ',f5.2,' ft/s ("5#)'/
     +       'itx ',f10.6,1x,f10.6,' Bump length = ',f4.2,' ft'/
     +       'itx ',f10.6,1x,f10.6,' Bump height = ',f4.2,' ft')

* Loop over trial values of the stream velocity

      do 100 iv = ivlo, ivhi
      
        v = vtry(iv)*vmax

* Determine the gravity wave number k0 consistent with the depth
* and current trial value of v
*
* Must numerically solve  v**2/gh = tanh(k0*h)/(k0*h)
*
* or  f = tanh(alpha)/alpha where f = vtry(iv)**2

        f = vtry(iv)**2

        alpha = 1.d0/f
        
        niter = 0
 10     niter = niter + 1
        diff = f*alpha - tanh(alpha)
        if(abs(diff).lt.1.d-10)then
          k0 = alpha/h
          write(6,15) v,vmax,f,alpha,niter
 15       format(' For (v=',f10.6,'/vmax=',f10.6,')**2=',f8.6,
     +           ', k0*h=',f10.6,'(niter=',i2,')')
          goto 20
        else 
          alpha = alpha - diff/(f-1./cosh(alpha)**2)
        endif
        if(niter.gt.20)then
          write(6,*)'Error: Failed to converge on alpha'
          write(6,*)'       Trying to solve for f=',f
          stop
        endif
        goto 10

* Compute amplitude of downstream standing wave

 20     kl = k0*l
        kh = k0*h
        zeta = -16.d0*pi*a * kl*cos(kl/2.d0)/(pi**2-kl**2)
     +         *sinh(kh)/(sinh(2.d0*kh)-2.d0*kh)

        if(abs(zeta).gt.h)then
          write(6,*)'Warning: Zeta=',zeta,' skipping this V'
          goto 100
        endif        
 
        write(6,22) kl,zeta
 22     format(' -> k0*L=',e14.6,' zeta=',e14.6)
        
* Write downstream sine-wave info to plot file

        xv = xmin * 0.5d0
        yv = y1 - yoff*(iv-1)
        write(1,25) zeta,k0,icol(iv),icol(iv),xv,yv,vtry(iv)
 25     format('zeta = ',e14.6/
     +         'k0 = ',e14.6/
     +         'sigma zd = [zeta]*sin([k0]*xd)'/
     +         'igset plci ',i1/
     +         'graph [nxd] xd zd l'/
     +         'igset txci ',i1/
     +         'itx ',f10.6,1x,f10.6,' V/V?MAX! = ',f4.2)

* Determine the tension wave number k0 consistent with the 
* trial value of v, using the approximate formula:

*          v = sqrt(k0*T/rho)

* where T is the surface tension and rho is water density:

*          gamma = T/rho = 0.0026 ft**3/s**2

        gamma = 2.6d-3
        k0 = v**2/gamma

* Compute amplitude of upstream standing wave 

        kl = k0*l
        kh = k0*h
        if(kh.lt.50.d0)then
          factor = sinh(kh)/(sinh(2.d0*kh)-2.d0*kh)
        else
          factor = 1.d0/(1.d0-2.d0*kh)
        endif
        zeta = +16.d0*pi*a * kl*cos(kl/2.d0)/(pi**2-kl**2)
     +         *factor

        write(6,22) kl,zeta
        
* Write upstream sine-wave info to plot file

        write(1,30) zeta,k0
 30     format('zeta = ',e14.6/
     +         'k0 = ',e14.6/
     +         'sigma zu = [zeta]*sin([k0]*xu)'/
     +         'graph [nxu] xu zu l')

 100  continue

      stop

      end