# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module puck_mod
  use tri_smooth_3d_mod
implicit none
  contains
!SMOOTH IT
  subroutine puck_yz (uu, coh, pp, res )
    integer i1, i2, n1, n2
    real uu(:,:,:), xx, xt, tt, coh, pp
    real, optional :: res(:,:,:)
    real dx(size(uu,1)-1,size(uu,2),size(uu,3)-1)
    real dt(size(uu,1)-1,size(uu,2),size(uu,3)-1)
    n1=size(uu,1)
    n2=size(uu,2)
    call wavekill_13( 1., 0., uu, dx)
    call wavekill_13( 0., 1., uu, dt)
    xx = 1.e-30
    tt = 1.e-30
    xt = 0.
    xt=xt+sum(dt*dx)
    xx=xx+sum(dx*dx)
    tt=tt+sum(dx*dx)
    coh = sqrt((xt/tt)*(xt/xx))
    pp = - xt/tt
    if (present(res)) then
      call wavekill_13( 1., pp, uu, res)
    end if
  end subroutine
  logical function puck_2d (uu, coh, pp,nbox )
    integer i1, i2, n1, n2,nbox(2),ib1,ib2,ie1,ie2
    real uu(:,:), xx, xt, tt, coh(:,:), pp(:,:)
    real dx(size(uu,1),size(uu,2))
    real dt(size(uu,1),size(uu,2))
    real :: a,eps
    puck_2d=.false.
    n1=size(uu,1)
    n2=size(uu,2)
    call wavekill( 1., 0., uu, dx)
    call wavekill( 0., 1., uu, dt)
    eps=epsilon(eps)
    eps = 1.e-30
!xt+=sum(dt*dx)
!xx+=sum(dx*dx)
!tt+=sum(dt*dt)
!coh = sqrt((xt/tt)*(xt/xx)
!if(.not. triangle_smooth_2d(dx,nbox*4)){
! write(0,*) "ERROR:      Trouble smoothing dx"
! return
!}
!
!if(.not. triangle_smooth_2d(dt,nbox*4)){
! write(0,*) "ERROR:      Trouble smoothing dt"
! return
!}
    do i2=1,size(uu,2)
      ib2=max(1,i2-nbox(2))
      ie2=min(size(uu,2),i2+nbox(2))
      do i1=1,size(uu,1)
        ib1=max(1,i1-nbox(1))
        ie1=min(size(uu,1),i1+nbox(1))
        xt=sum(dt(ib1:ie1,ib2:ie2)*dx(ib1:ie1,ib2:ie2))
        xx=eps+sum(dx(ib1:ie1,ib2:ie2)*dx(ib1:ie1,ib2:ie2))
        tt=eps+sum(dt(ib1:ie1,ib2:ie2)*dt(ib1:ie1,ib2:ie2))
        coh(i1,i2)=sqrt((xt/tt)*(xt/xx))
        pp(i1,i2)=-xt/tt
      end do
    end do
    puck_2d=.true.
  end function
  subroutine puck_xz (uu, coh, pp, res )
    integer i1, i2, n1, n2
    real uu(:,:,:), xx, xt, tt, coh, pp
    real, optional :: res(:,:,:)
    real dx(size(uu,1)-1,size(uu,2)-1,size(uu,3))
    real dt(size(uu,1)-1,size(uu,2)-1,size(uu,3))
    n1=size(uu,1)
    n2=size(uu,2)
    call wavekill_12( 1., 0., uu, dx)
    call wavekill_12( 0., 1., uu, dt)
    xx = 1.e-30
    tt = 1.e-30
    xt = 0.
    xt=xt+sum(dt*dx)
    xx=xx+sum(dx*dx)
    tt=tt+sum(dx*dx)
    coh = sqrt((xt/tt)*(xt/xx))
    pp = - xt/tt
    if (present(res)) then
      call wavekill_12( 1., pp, uu, res)
    end if
  end subroutine
! vv = (aa Dx + pp Dt) uu
!
  subroutine wavekill( aa, pp, uu, vv )
    integer i1,i2, n1,n2,n3
    real aa, pp, s11, s12, s21, s22, uu(:,:), vv(:,:)
    s11 = -aa-pp
    s12 = aa-pp
    s21 = -aa+pp
    s22 = aa+pp
    vv=0
    n1=size(uu,1)
    n2=size(uu,2)
    do i1= 1, n1-1
! vv is one point shorter than uu on both axes.
      do i2= 1, n2-1
        vv(i1,i2) = vv(i1,i2) + uu(i1 ,i2) * s11 + uu(i1 ,i2&
          &+1) * s12 +uu(i1+1,i2) * s21 + uu(i1+1,i2+1) * s22
      end do
    end do
    vv(n1,:)=vv(n1-1,:)
    vv(:,n2)=vv(:,n2-1)
  end subroutine
  subroutine wavekill_12( aa, pp, uu, vv )
    integer i1,i2, n1,n2,n3
    real aa, pp, s11, s12, s21, s22, uu(:,:,:), vv(:,:,:)
    s11 = -aa-pp
    s12 = aa-pp
    s21 = -aa+pp
    s22 = aa+pp
    vv=0
    n1=size(uu,1)
    n2=size(uu,2)
    n3=size(uu,3)
    do i1= 1, n1-1
! vv is one point shorter than uu on both axes.
      do i2= 1, n2-1
        vv(i1,i2,:) = vv(i1,i2,:) + uu(i1 ,i2,:) * s11 + uu(i1&
          & ,i2+1,:) * s12 +uu(i1+1,i2,:) * s21 + uu(i1+1,i2+1,:) * s22
      end do
    end do
  end subroutine
  subroutine wavekill_13( aa, pp, uu, vv )
    integer i1,i2, n1,n2,n3,i3
    real aa, pp, s11, s13, s31, s33, uu(:,:,:), vv(:,:,:)
    s11 = -aa-pp
    s13 = aa-pp
    s31 = -aa+pp
    s33 = aa+pp
    vv=0
    n1=size(uu,1)
    n2=size(uu,2)
    n3=size(uu,3)
    do i1= 1, n1-1
! vv is one point shorter than uu on both axes.
      do i3= 1, n3-1
        vv(i1,:,i3) = vv(i1,:,i3) + uu(i1 ,:,i3) * s11 + uu(i1&
          & ,:,i3+1) * s13 +uu(i1+1,:,i3) * s31 + uu(i1+1,:,i3+1) * s33
      end do
    end do
  end subroutine
end module
