	program CircleWave
	
	real*8 grid(120,120,500)
	integer allow(120,120)
	real*8 load(120,120)
	integer shape(2)
	real*8 dx,dy,dt,Lx,Ly,Lto,t
	real*8 e
	integer i,j,k,nx,ny,nt,rank
	
	e=2.718d0
	rank=2
	
	
	Lx=15.0d0
	Ly=15.0d0
	Lto=20.0d0
	
	nx=120
	ny=120
	nt=500
	
	shape(1)=nx
	shape(2)=ny
	
	
	dx=Lx/nx
	dy=Ly/ny
	dt=Lto/nt

C loading up allow array

	do i=1,nx
	  do j=1,ny
	    x=i*dx
	    y=j*dy
	    if((x-Lx/2)**2+(y-Ly/2)**2 .lt. (Lx**2)/4) then
	      allow(i,j)=1
	      write(81,*) x,y
	    endif
	  enddo
	enddo	
	
C done	
	
C loading initial waveform with zero initial velocity	
	
	do i=1,nx
	   do j=1,ny
	     if(allow(i,j) .eq. 1) then
	      grid(i,j,1)=e**-((i*dx-Lx/3)**2+(j*dy-Lx/3)**2)
     &	+e**-((i*dx-2*Lx/3)**2+(j*dy-2*Lx/3)**2)
	      grid(i,j,2)=grid(i,j,1)
	     endif
	   enddo
	enddo
	
C done loading initial waveform

C boundary conditions are all zero

C now for main loop

	do i=2,nt-1
	   do j=2,nx-1
	     do k=2,ny-1
	       if(allow(j,k) .eq. 1) then
	        grid(j,k,i+1)=2*grid(j,k,i)-grid(j,k,i-1)
     &	+((dt/dx)**2)*(grid(j+1,k,i)-2*grid(j,k,i)+grid(j-1,k,i))
     &	+((dt/dy)**2)*(grid(j,k+1,i)-2*grid(j,k,i)+grid(j,k-1,i))
	       endif
	     enddo
	   enddo
	enddo

	write(*,*) 'halfway'

C main body done.  Now load to file


	do i=1,nt
	  t=(i-1)*dt
	   do j=1,nx
	      do k=1,ny
	        if(allow(j,k) .eq. 1) then
	          load(j,k)=grid(j,k,i)
	        endif
	      enddo
	   enddo
	  call gft_out('Ripple',t,shape,rank,load)
	enddo
	
	write(*,*) (((dt/dx)**2)+((dt/dy)**2))**(0.5)
	
	stop
	end
