* Program to simulate effect of wave motion on a sampling of
* water particles. (Uses deep-water approximation in which
* particles travel in circles.)

      implicit none

* Declare variable types

      integer it,ix,isign,iy,nt,nx,ny
      double precision xmin,xmax,ymin,ymax,lambda,f,a,x0,pi,
     +                 t,t0p,x0p,y0p,xp,yp,g,c,period,omega

* Define constants (pi and gravitational acceleration)

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

* Define wavelength and amplitude

      lambda = 10.d0
      a = 0.5d0

* Compute wave velocity, period, and angular velocity

      c = sqrt(g*lambda/2.d0/pi)
      period = lambda/c
      omega = 2.d0*pi/period

* Define x position for which the crest occurs at t=0

      x0 = 0.d0

* Define boundaries of plot

      f = 1.1d0
      xmin = x0 - f * lambda
      xmax = x0 + f * lambda
      ymin = -5.d0 * a
      ymax = 2.0d0*a

* Open and initialize plot input file

      open(unit=1,file='wave1.kumac',status='new')
      write(1,1)
 1    format('zone 1 2'/
     +       'igset txal 23'/
     +       'igset txci 1'/
     +       'fort/file 66 wave1.ps'/
     +       'meta 66 -111')


* Loop over intervals of time, starting at t=0 
* (1 picture per time interval)

      nt = 3
      do it = 0, nt

* Define current time

        t = it * period / 4

* Define border for this plot

        write(1,5) xmin, xmax, ymin, ymax, 
     +             x0, ymax*0.85d0, t/period
 5      format('null',4(1x,f14.6)/
     +         'igset txci 1'/
     +         'itx',1x,f14.6,1x,f14.6,1x,'''t = ',f5.3,'*T''')

* Loop over intervals of x (horizontal position)

        nx = 16
        do ix = 0, nx

* Look to both left and right of x=x0 (center of picture)

          do isign = -1, 1, 2

* Define nominal still-water x position of particle

            x0p = x0 + (ix * lambda / nx) * isign

* Define color of particles with common x0p 
* (vertical column when t=t0p)

            write(1,8) mod(ix,7) + 1
 8          format('igset txci ',i1)
        

* Find the time t0' when this particle would be at the
* top of its orbit (t0' > 0 to right of x=x0, t0' < 0 to left)

            t0p =  (x0p-x0) / c

* Loop over intervals of y (depth of still-water position)

            ny = 20
            do iy = 0, ny

* Define nominal still-water y position of particle

              y0p = - iy * a / 2.d0

* Find its x and y positions at current time

              xp = x0p + a*exp(2.d0*pi*y0p/lambda)
     +                    *sin(omega*(t-t0p))
              yp = y0p + a*exp(2.d0*pi*y0p/lambda)
     +                    *cos(omega*(t-t0p))

* Put marker on picture for this point 
* (but only if it's within bounds of frame)

              if(xp.gt.xmin.and.xp.lt.xmax.and.yp.gt.ymin)then
                write(1,10) xp,yp,mod(iy,10)
 10             format('itx',1x,f14.6,1x,f14.6,3x,i1)
              endif

            enddo
          enddo
        enddo

* If not last picture, put in wait statement
 
        if(it.ne.nt) write(1,20)
 20     format('wait')

      enddo

      stop
      end