* Program to simulate effect of wave motion on a sampling of * water particles. * (Uses analysis in which particles travel in ellipses.) 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, + h * Define constants (pi and gravitational acceleration) pi = acos(-1.d0) g = 9.8d0 * Define wavelength, amplitude, and depth of channel lambda = 10.d0 a = 0.5d0 h = 1.d0 * Compute wave velocity, period, and angular velocity c = sqrt(g*lambda/2.d0/pi*tanh(2.d0*pi*h/lambda)) 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 = -h ymax = 2.0d0*a * Open and initialize plot input file open(unit=1,file='wave3.kumac',status='new') write(1,1) 1 format('zone 1 2'/ + 'igset txal 23'/ + 'igset txci 1'/ + 'fort/file 66 wave3.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 / 4.d0 * Find its x and y positions at current time xp = x0p + a*cosh(2.d0*pi*(y0p+h)/lambda) + /sinh(2.d0*pi*h/lambda) + *sin(omega*(t-t0p)) yp = y0p + a*sinh(2.d0*pi*(y0p+h)/lambda) + /sinh(2.d0*pi*h/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