	program doubleslit

	real*8 u(100,100,700)
	real*8 load(100,100)

	real*8 coords(4)

	integer allow(100,100)
	integer shape(2)
	real*8 dx,dy,dt,Lx,Ly,Lto,t,x,y
	real*8 e
	integer i,j,k,nx,ny,nt,rank
	
	e=2.718d0
	rank=2


	Lx=10.0d0
	Ly=10.0d0
	Lto=7d0

	coords(1)=0.0d0
	coords(2)=Lx
	coords(3)=0.0d0
	coords(4)=Ly

	nx=100
	ny=100
	nt=700

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


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

C loading logical array

	do i=1,nx

	  x=i*dx

	  do j=1,ny

	    y=j*dy

	    if((y .lt. 6.2 .and. y .gt. 5.8 .and. x .lt. 4.1) .or.
     &	(y .lt. 6.2 .and. y .gt. 5.8 .and. x .gt. 4.7 .and.
     &	x .lt. 5.3) .or.
     &	(y .lt. 6.2 .and. y .gt. 5.8 .and. x .gt. 5.9)) then
	        allow(i,j)=1
		write(91,*) x,y
	    endif
	  enddo
	enddo
C done loading logical

C loading initial waveform with zero initial velocity
	
	do i=1,nx

	   do j=1,ny

	    if(allow(i,j) .eq. 0) then
	     u(i,j,1)=exp(-20*(j*dy-1)**2)+exp(-20*(j*dy-1.5)**2)
     &	+exp(-20*(j*dy-2)**2)+exp(-20*(j*dy-2.5)**2)
     &	+exp(-20*(j*dy-3)**2)+exp(-20*(j*dy-3.5)**2)
     &	+exp(-20*(j*dy-4)**2)

	     u(i,j,2)=u(i,j,1)+dt*2*20*(j*dy-1)*exp(-20*(j*dy-1)**2)
	     u(i,j,2)=u(i,j,2)+dt*2*20*(j*dy-1.5)*exp(-20*(j*dy-1.5)**2)
	     u(i,j,2)=u(i,j,2)+dt*2*20*(j*dy-2)*exp(-20*(j*dy-2)**2)
	     u(i,j,2)=u(i,j,2)+dt*2*20*(j*dy-2.5)*exp(-20*(j*dy-2.5)**2)
	     u(i,j,2)=u(i,j,2)+dt*2*20*(j*dy-3)*exp(-20*(j*dy-3)**2)
	     u(i,j,2)=u(i,j,2)+dt*2*20*(j*dy-3.5)*exp(-20*(j*dy-3.5)**2)
	     u(i,j,2)=u(i,j,2)+dt*2*20*(j*dy-4)*exp(-20*(j*dy-4)**2)
	    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 k=2,ny-1


	if(allow(1,k) .eq. 0) then
	        u(1,k,i+1)=2*u(1,k,i)-u(1,k,i-1)
     &	+((dt/dx)**2)*(u(2,k,i)-2*u(1,k,i)+u(nx,k,i))
     &	+((dt/dy)**2)*(u(1,k+1,i)-2*u(1,k,i)+u(1,k-1,i))
	      endif


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


	if(allow(nx,k) .eq. 0) then
	        u(nx,k,i+1)=2*u(nx,k,i)-u(nx,k,i-1)
     &	+((dt/dx)**2)*(u(1,k,i)-2*u(nx,k,i)+u(nx-1,k,i))
     &	+((dt/dy)**2)*(u(nx,k+1,i)-2*u(nx,k,i)+u(nx,k-1,i))
	      endif


	   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)=u(j,k,i)
	      enddo
	   enddo
	  call gft_out_bbox('doubleslit',t,shape,rank,coords,load)
	  write(*,*) i
	enddo

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

	stop
	end
