nx = 100 & ny = 100 & xrange = [0,20]
mx = 300 & my = 300 & yrange = [0,20]
nt = 100 & tstep = 2

;; Boxcar shape
rj = beselj(findgen(mx)/7,0)-beselj(findgen(mx)/5,0)
zj = exp(-(findgen(my)/7)^2/(2*5.^2))

kr = 0.05*findgen(mx)
kz = kr
rr = findgen(nx)*(xrange(1)-xrange(0))/nx
zz = findgen(ny)*(yrange(1)-yrange(0))/ny

rd = fltarr(nx,mx)
zd = fltarr(ny,my)
for i = 0L, mx-1 do rd(*,i) = rj(i)*beselj(rr*kr(i),1)
for i = 0L, my-1 do zd(*,i) = zj(i)*cos(zz*kz(i))

ones = fltarr(nx)+1
t0 = systime(1)
images = fltarr(nx, ny, nt)
for k = 0, nt-1 do begin 
    t = k*tstep/5. ;; Time step

    t1 = systime(1)
    print, 'Frame '+strtrim(k,2), format='(A,$)'
    bb = fltarr(nx,ny)
    for j = 0, ny-1 do begin
        b = fltarr(nx)
        for i = 0, my-1 do begin
            cgrid = ones # cos(sqrt(kz(i)^2+kr^2)*t)
            b = b + zd(j,i)*total(rd*cgrid,2)
        endfor
        bb(*,j) = b
        plot, rr, abs(b), yrange=[.01,1000], /ylog, $
          title='Time = '+strtrim(t,2)+'   z = '+strtrim(zz(j),2)
        wait, 0.05
        if j MOD 5 EQ 0 then print, '.', format='(A,$)'
    endfor
    minute = (systime(1)-t1)/60.
    print, 'done ('+strtrim(minute,2)+' min)'
    images(*,*,k) = bb
    
;    plot, rr, b
;    oplot, xx, y*0+2, linestyle=1
;    oplot, xx, y*0-2, linestyle=1
;    oplot, (findgen(20)-10)*2.5 + t, fltarr(20)+2, psym=8
;    oplot, (findgen(20)-10)*2.5 + t, fltarr(20)-2, psym=8
    
endfor

print, 'TOTAL TIME: ', (systime(1)-t0)/3600., ' hours'

end