Program Refraction

	real*8 grid(300,100,1000)
	real*8 load(300,100)
	real*8 coords(4)
	integer shape(2)
	real*8 dx,dy,dt,Lx,Ly,Lto,x,y,t,L1,L2
	real*8 e
	integer i,j,k,nx,ny,nt,rank

	e=2.718d0
	rank=2


	Lx=20.0d0
	Ly=10.0d0
	Lto=10.0d0

	coords(1)=dx
	coords(2)=Lx
	coords(3)=dy
	coords(4)=Ly

	nx=300
	ny=100
	nt=1000

	shape(1)=nx
	shape(2)=ny


	dx=Lx/nx
	dy=Ly/ny
	dt=Lto/nt

	L1=(dt/dx)**2
	L2=(dt/dy)**2

C loading initial waveform

	do i=1,nx
	   do j=1,ny
	     grid(i,j,1)=exp(-10*(i*dx-15)**2)
	     grid(i,j,2)=grid(i,j,1)-dt*2*10*(i*dx-15)*grid(i,j,1)
	   enddo
	enddo

C done loading initial waveform

C boundary conditions are all zero and are satisfied automatically
C as a real*8 array comes from the store full of zeros.

C now for main loop

	do i=2,nt-1
	   do j=2,nx-1
	   x=j*dx

C make periodic in Y.  At y=0
	grid(j,1,i+1)=2*grid(j,1,i)-grid(j,1,i-1)
     &	+L1*v(x,dy)*(grid(j+1,1,i)-2*grid(j,1,i)+grid(j-1,1,i))
     &	+L2*v(x,dy)*(grid(j,2,i)-2*grid(j,1,i)+grid(j,ny,i))
C done


	     do k=2,ny-1
	     y=k*dy
	        grid(j,k,i+1)=2*grid(j,k,i)-grid(j,k,i-1)
     &	+L1*v(x,y)*(grid(j+1,k,i)-2*grid(j,k,i)+grid(j-1,k,i))
     &	+L2*v(x,y)*(grid(j,k+1,i)-2*grid(j,k,i)+grid(j,k-1,i))

	     enddo


C make periodic in Y.  At y=Ly
	y=ny*dy
	grid(j,ny,i+1)=2*grid(j,ny,i)-grid(j,ny,i-1)
     &	+L1*v(x,y)*(grid(j+1,ny,i)-2*grid(j,ny,i)+grid(j-1,ny,i))
     &	+L2*v(x,y)*(grid(j,1,i)-2*grid(j,ny,i)+grid(j,ny-1,i))
C done

	   enddo
	   write(*,*) i
	enddo

C main body done.  Now load to file


	do i=1,nt
	  t=(i-1)*dt
	   do j=1,nx
	      do k=1,ny
	        load(j,k)=grid(j,k,i)
	      enddo
	   enddo
	  call gft_out_bbox('refraction',t,shape,rank,coords,load)
	  write(*,*) i
	enddo

	write(*,*) ((dt/dx)**2)+((dt/dy)**2)

	stop
	end

C Here is our velocity function.  I call it velocity but it really represents
C velocity squared.

	function v(x,y)

	   real*8 x,y

	   if(y .gt. 2*x-10) then
		v=5
	   else
	        v=1
	   endif
	return
	end