# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module build_steering3d_filters
  use sorting_mod
  use gsmooth_mod
  use steering3d
  use nhelix
  implicit none
  integer, private, save :: n(3),ns,nwind
  integer, private, save,pointer :: use_refs(:)
  real,private,save :: o(3),d(3),min_amp,max_amp,range,hwidth
  contains
  subroutine default_steering_filters_3d()
    n=1
    o=0
    d=1
    ns=3
    nwind=2
    allocate(use_refs(1))
    min_amp=.1
    max_amp=.99
    range=1
    hwidth=2.5
  end subroutine
  subroutine init_steering_filters_3d(n_in,o_in,d_in,use_refs_in,ns_in&
    &,nwind_in,min_amp_in,max_amp_in,range_in,hwidth_in)
    integer, optional :: n_in(:),use_refs_in(:),ns_in,nwind_in
    real, optional :: d_in(:),o_in(:),min_amp_in,max_amp_in,range_in&
      &,hwidth_in
    if (present(n_in)) then
      n(1:size(n_in))=n_in
    end if
    if (present(o_in)) then
      o(1:size(o_in))=o_in
    end if
    if (present(d_in)) then
      d(1:size(d_in))=d_in
    end if
    if (present(ns_in)) then
      ns=ns_in
    end if
    if (present(nwind_in)) then
      nwind=nwind_in
    end if
    if (present(min_amp_in)) then
      min_amp=min_amp_in
    end if
    if (present(max_amp_in)) then
      max_amp=max_amp_in
    end if
    if (present(range_in)) then
      range=range_in
    end if
    if (present(hwidth_in)) then
      hwidth=hwidth_in
    end if
    if (present(use_refs_in)) then
      deallocate(use_refs)
      allocate(use_refs(size(use_refs_in)))
      use_refs=use_refs_in
    end if
  end subroutine
  subroutine build_3d_steer_p_map(refs,oref,dref,pvalx,pvaly)
    real :: refs(:,:,:),pvalx(:,:,:),pvaly(:,:,:)
    real :: oref(:),dref(:)
    real , allocatable :: slopex(:,:,:),slopey(:,:,:),weight(:,:,:)
    integer :: nref(3),i1,i2,nuse,n1,n2,iuse
    real :: posx,posy,fx,fy,posz,fz,sx,sy,pos1,pos1x,pos1y,pos2x,pos2y&
      &,f
    real :: dist,dist1
    integer :: ix,iy,iz,i3,ia,loc2,loc1
    nref(1)=size(refs,1)
    nref(2)=size(refs,2)
    nref(3)=size(refs,3)
    nuse=size(use_refs)
    n1=size(refs,1)
    n2=size(refs,2)
!first find slopes along the reflectors
    allocate(slopex(n1,n2,nuse),slopey(n1,n2,nuse),weight(n(1),n(2),n(3&
      &)))
    weight=0.
    slopex=-999.
    slopey=-999.
    do iuse=1,nuse
      do i2=1,n2
        if (all(refs(1:2,i2,use_refs(iuse))>0.)) then
          slopex(1,i2,iuse)=(refs(2,i2,use_refs(iuse))-refs(1,i2,use_refs&
            &(iuse)))/dref(1)
        end if
        if (all(refs(n1-1:n1,i2,use_refs(iuse))>0.)) then
          slopex(n1,i2,iuse)=(refs(n1,i2,use_refs(iuse))-refs(n1-1,i2&
            &,use_refs(iuse)))/dref(1)
        end if
        do i1=2,n1-1
          if (all(refs(i1-1:i1+1,i2,use_refs(iuse))>0.)) then
            slopex(i1,i2,iuse)=(refs(i1+1,i2,use_refs(iuse))-refs(i1-1&
              &,i2,use_refs(iuse)))/dref(1)/2.
          end if
        end do
      end do
      do i1=1,n1
        if (all(refs(i1,1:2,use_refs(iuse))>0.)) then
          slopey(i1,1,iuse)=(refs(i1,2,use_refs(iuse))-refs(i1,1,use_refs&
            &(iuse)))/dref(2)
        end if
        if (all(refs(i1,n2-1:n2,use_refs(iuse))>0.)) then
          slopey(i1,n2,iuse)=(refs(i1,n2,use_refs(iuse))-refs(i1,n2-1&
            &,use_refs(iuse)))/dref(2)
        end if
        do i2=2,n2-1
          if (all(refs(i1,i2-1:i2+1,use_refs(iuse))>0.)) then
            slopey(i1,i2,iuse)=(refs(i1,i2+1,use_refs(iuse))-refs(i1,i2&
              &-1,use_refs(iuse)))/dref(2)/2.
          end if
        end do
      end do
    end do
!convert the slopes into velocity space coordinate system
    where(slopex.ne.-999.)
      slopex=slopex*d(2)/d(1)
end where
    where(slopey.ne.-999.)
      slopey=slopey*d(3)/d(1)
end where
!put the slopes into velocity space cube
    pvalx=0
    pvaly=0
    do i3=1,nuse
      do i2=1,n(2)
        posy=o(3)+d(3)*(i2-1)
        fy=(posy-oref(2))/dref(2)+1.
        iy=fy
        fy= fy-iy
        do i1=1,n(1)
          posx=o(2)+d(2)*(i1-1)
          fx=(posx-oref(1))/dref(1)+1.
          ix=fx
          fx= fx-ix
          if (ix<1 .or. ix >= nref(1) .or. iy<1 .or. iy>=nref(2)) then
            cycle
          end if
          if (any(slopex(ix:ix+1,iy:iy+1,i3).eq.-999.)) then
            cycle
          end if
          if (any(slopey(ix:ix+1,iy:iy+1,i3).eq.-999.)) then
            cycle
          end if
          posz= (1.-fx) * (1. -fy)*refs(ix,iy,use_refs(i3))+(fx) * (1.&
            & -fy)*refs(ix+1,iy,use_refs(i3))+(1.-fx) * (fy)*refs(ix&
            &,iy+1,use_refs(i3))+(fx) * (fy)*refs(ix+1,iy+1,use_refs(i3))
          sx= (1.-fx) * (1. -fy)*slopex(ix,iy,i3)+(fx) * (1. -fy)*slopex&
            &(ix+1,iy,i3)+(1.-fx) * (fy)*slopex(ix,iy+1,i3)+(fx) * &
            &(fy)*slopex(ix+1,iy+1,i3)
          sy= (1.-fx) * (1. -fy)*slopey(ix,iy,i3)+(fx) * (1. -fy)*slopey&
            &(ix+1,iy,i3)+(1.-fx) * (fy)*slopey(ix,iy+1,i3)+(fx) * &
            &(fy)*slopey(ix+1,iy+1,i3)
          fz=(posz-o(1))/d(1)+1.
          iz=fz
          if (iz< 1 .or. iz >= n(1)) then
            cycle
          else
            pvalx(iz:iz+1,i1,i2)=sx
            pvaly(iz:iz+1,i1,i2)=sy
            weight(iz:iz+1,i1,i2)=1.
          end if
        end do
      end do
    end do
!now interpolate down
!LOOP OVER OUTPUT INTERPOLATING SLOPE BETWEEN REFLECTORS
    do i3=1,n(3)
      do i2=1,n(2)
        loc1=n(1)
        pos1x=pvalx(1,i2,i3)
        i1=n(1)
        pos1y=pvaly(1,i2,i3)
        do while (i1>=1)
          do while (i1>=1 .and. weight(i1,i2,i3).eq.0.)
            i1= i1-1 !FIND NEXT SLOPE
          end do
          if (i1.ne.0) then
!IF WE DIDN'T REACH THE BOTTOM
            pos2y=pvaly(i1,i2,i3)
            pos2x=pvalx(i1,i2,i3)
            loc2=i1
            dist1=loc1-loc2
            do ia=loc2+1,loc1-1
!INTERPOLATE BETWEEN THE TWO POINTS
              f=real(ia-loc2)/dist1
              pvalx(ia,i2,i3)=(1.-f)*pos2x+f*pos1x
              pvaly(ia,i2,i3)=(1.-f)*pos2y+f*pos1y
            end do
            loc1=loc2
            pos1x=pos2x
            pos1y=pos2y
!RESET
          end if
          i1= i1-1
        end do
        pvalx(1:loc1,i2,i3)=pvalx(loc1,i2,i3) !COPY LAST SLOPE DOWN
        pvaly(1:loc1,i2,i3)=pvaly(loc1,i2,i3) !COPY LAST SLOPE DOWN
      end do
    end do
    call gauss_smooth(pvalx,5,9,9,1.)
    call gauss_smooth(pvaly,5,9,9,1.)
    deallocate(slopex,slopey,weight)
  end subroutine
  subroutine create_3d_steer_2filt_op(pvalx,pvaly,amp,filtx,filty)
    real :: pvalx(:,:,:),amp(:,:,:),pvaly(:,:,:)
    type(nfilter) :: filtx,filty
    integer i1,i2,i,nh,i3,nx,ny
    real :: ox,dx,oy,dy,range,a(2),h(2),p(2)
    i=0
!write not lets just set it to 1000 sampling
    nx=1000
    ny=1000
    ox=minval(pvalx)
    range=maxval(pvalx)-ox
    dx=range/(nx-1)
    oy=minval(pvaly)
    range=maxval(pvaly)-oy
    dy=range/(ny-1)
    allocate(filtx%hlx(nx+1),filtx%pch(product(n)))
    allocate(filty%hlx(ny+1),filty%pch(product(n)))
    nh=ceiling(hwidth*2.)
    h=hwidth
    a=amp(1,1,1)
    call steer3d_init(h,a,n(1:2))
    do i1=1,nx
      call allocatehelix(filtx%hlx(i1),nh)
      call allocatehelix(filty%hlx(i1),nh)
      call create_steer3d_op((/ox+dx*(i1-1),oy+dy*(i1-1)/),filtx%hlx(i1&
        &),filty%hlx(i1))
    end do
    i=0
    do i3=1,n(3)
      do i2=1,n(2)
        do i1=1,n(1)
          i= i+1
          filtx%pch(i)=min(nx,max(1,nint((pvalx(i1,i2,i3)-ox)/dx+1.)))
          filty%pch(i)=min(ny,max(1,nint((pvaly(i1,i2,i3)-oy)/dy+1.)))
        end do
      end do
    end do
    call allocatehelix(filtx%hlx(nx+1),1)
    call allocatehelix(filty%hlx(ny+1),1)
    filtx%hlx(nx+1)%flt=0.
    filtx%hlx(nx+1)%lag=0
    filty%hlx(nx+1)%lag=0
    filty%hlx(nx+1)%flt=0.
  end subroutine
end module
