# Program acoustic 2D coded in Julia
# nx, nz, nt = 10256, 10256, 2001
    using Gadfly

function acoustic(n)



    X = Int32[]
    Y = Int32[]
    Z = Float32[]

    mxnx = 203
    mxny = 203
    mxnz = 203
    mxnt = 2003


    nx = mxnx - 2
    ny = mxny - 2
    nz = mxnz - 2
    nt = mxnt - 2


    dx   = 0.02               # (km)
    dz   = 0.02               # (km)     
    dt   = 0.002              # (s)


    fmax = 10                 # (Hz)
    xmax = (nx - 1) * dx      # (km)
    zmax = (nz - 1) * dz      # (km)
    tmax = (nt - 1) * dt      # (s)


    isnap = 500
 
    u1 = zeros(Float32,mxnx,mxnz)                 
    u2 = zeros(Float32,mxnx,mxnz)                 
    u3 = zeros(Float32,mxnx,mxnz)                 
    vp = zeros(Float32,mxnx,mxnz)                 
     f = zeros(Float32,mxnx,mxnz)                 
     w = zeros(Float32,nt)

     rzloc = 1

     for iz = 2:nz+1
         for ix = 2:nx+1		
            vp[ix,iz] = 1.5
         end
      end


#     Gaussian Source Time Function to set w
#      w = fdgaus (w, fmax, dt, nt)
#      println("w=$w")
#      function fdgaus(w,cutoff,dt,nt)
#      println("cutoff=$cutoff, dt=$dt, nt=$nt")

       cutoff = fmax
       phi = 4 * atan(1.0)
       a = phi * (5.0 * cutoff / 8.)^2
       amp = sqrt(a/phi)

       for i = 1:nt
           t = (i - 1) * dt
           arg = -a * t^2
           if arg < -32.0 arg = -32.0 end
           w[i] = amp * exp(arg)
       end

       t0 = 0
#      for i = 1:nt  
#       if w[i] < 0.001 * w[1]]
#        icut = i
#        t0 = (icut - 1) * dt 
#        break
#       end
#      end 

        i = 1
        while true
            if w[i] < 0.001 * w[1]
                icut = i
                t0 = (icut - 1) * dt
                break
            end
            i += 1
        end

        for i = 1:nt
            t = (i-1) * dt
            t = t - t0
            arg = -a * t^2
            if arg < -32.0 arg = -32.0 end
            w[i] = -2. * sqrt(a) * a * t * exp(arg) / sqrt(phi)
        end

        smax = 0.0

        for i = 1:nt
            if abs(w[i]) > smax smax = abs(w[i]) end
        end

        for i = 1:nt
            w[i] = w[i] / smax
            if abs(w[i]) < 0.0001 w[i] = 0.0 end
        end

#       println("smax=$smax")
#       println("w=$w")
#       p0 = plot(x = 1:mxnt,y = w,Geom.rectbin)
#       display(p0)
#
#      End of fdgaus
#

      sxloc = trunc(mxnx / 2)               # Source location in x-direction
      szloc = trunc(mxnz / 2)               # Source location in z-direction (depth)          
      f[sxloc,szloc] = 1

tic()
     for it = 1:nt              # Marching over Time
       println("time it=$it")
         for iz = 2:nz+1
             for ix = 2:nx+1
                 u3[ix,iz] = 2 * u2[ix,iz] - u1[ix,iz] + ((u2[ix+1,iz] -2.* u2[ix,iz] + u2[ix-1,iz]) / dx^2  + (u2[ix,iz+1] -2. * u2[ix,iz] + u2[ix,iz-1]) / dz^2 ) * (vp[ix,iz] * dt )^2
             end
         end

         ix = sxloc
         iz = szloc
         u3[ix,iz] = u3[ix,iz] + f[ix,iz] * w[it] * (vp[ix,iz] * dt )^2



         for iz = 1:mxnz
             for ix = 1:mxnx
               u1[ix,iz] = u2[ix,iz]
             end
         end
         
         for iz = 1:mxnz
             for ix = 1:mxnx
               u2[ix,iz] = u3[ix,iz]
             end
         end


#         uncomment the following to display pictures
#         if it % isnap == 1 
#             for x in linspace(1,mxnx,mxnx), y in linspace(1,mxnz,mxnz)
#                  push!(X,x)
#                  push!(Y,y)
#                  push!(Z,u3[x,y])
#              end
#              p = plot(x=X,y=Y,color=Z,Geom.rectbin)
#              display(p)    

#         end

     end                  # time loop
toc()


end # acoustic end