# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module build_steering2d_filters
  use sorting_mod
  use gsmooth_mod
  use steering
  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_2d()
    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_2d(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:min(3,size(n_in)))=n_in
    end if
    if (present(o_in)) then
      o(1:min(3,size(o_in)))=o_in
    end if
    if (present(d_in)) then
      d(1:min(3,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_2d_steer_p_map(refs,pval)
    integer n1,n2,i1,i2,ierr,getch,nuse,bi,ei,ia,ib,loc1,loc2
    real :: o1,o2,d1,d2,wt1,wt2
    real, allocatable :: slopes(:,:),tempr(:),temp2r(:),weight(:,:)
    complex :: refs(:,:)
    real :: pval(:,:)
    real :: pos1,pos2,dist,smooth,f,f2,f1,dist1
    integer :: npts,nref,nmid,izero
    integer :: n_t(3),n_z(3)
    real :: o_z(3),d_z(3),o_t(3),d_t(3),x,z
    real,pointer :: vel(:,:),map(:,:),tt(:,:,:)
    npts=size(refs,1)
    nref=size(refs,2)
    nuse=size(use_refs)
    if (any(use_refs > nref ) .or. any(use_refs < 1)) then
      call erexit("steering filter reflectors outside bounds")
    end if
    o1=o(1)
    o2=o(2)
    d1=d(1)
    d2=d(2)
    n1=n(1)
    n2=n(2)
    allocate(slopes(npts,nuse))
    allocate(tempr(nwind*2+1),temp2r(npts))
    do i2=1,nuse
      slopes(1,i2)=(aimag(refs(2,use_refs(i2)))-aimag(refs(1,use_refs&
        &(i2)))+.0001)/ (real(refs(2,use_refs(i2)))-real(refs(1&
        &,use_refs(i2))+.0001))
      do i1=2,npts-1
        slopes(i1,i2)=((aimag(refs(i1+1,use_refs(i2)))-aimag(refs(i1-1&
          &,use_refs(i2)))+.0001)/ (real(refs(i1+1,use_refs(i2)))&
          &-real(refs(i1-1,use_refs(i2)))+.0001))
      end do
      slopes(npts,i2)=(aimag(refs(npts,use_refs(i2)))-aimag(refs(npts-1&
        &,use_refs(i2))+.0001))/ (real(refs(npts,use_refs(i2)))-real&
        &(refs(npts-1,use_refs(i2))+.0001))
!median smooth it
      temp2r=slopes(:,i2)
      where(temp2r>3.)
        slopes(:,i2)=3.
end where
      where(temp2r<-3.)
        slopes(:,i2)=-3.
end where
      do i1=1,npts
        bi=max(1,i1-nwind)
        ei=min(npts,i1+nwind)
        nmid=ei-bi+1
        tempr(1:nmid)=slopes(bi:ei,i2)
        call quick_sort(tempr(1:nmid),nmid)
        temp2r(i1)=tempr(ceiling(nmid/2.))
      end do
      slopes(:,i2)=temp2r
    end do
    slopes=slopes*d2/d1
!now we interpolae
    allocate(weight(n1,n2))
    pval=0
    weight=0
! IMPOSE A PYRAMID SMOOTHER OVER REFLECTOR POSITION
    smooth=ns
    do i2=1,nuse
      do i1=1,npts
        pos1=(aimag(refs(i1,use_refs(i2)))-o1)/d1+1.
        pos2=(real(refs(i1,use_refs(i2)))-o2)/d2+1.
        loc1=floor(pos1)
        f1=pos1-loc1
        loc2=floor(pos2)
        f2=pos2-loc2
        do ib=max(loc2-ns,1),min(loc2+ns,n2)
          wt2=max(smooth-abs(pos2-ib),0.)
          do ia=max(loc1-ns,1),min(loc1+ns,n1)
            wt1=max(smooth-abs(pos1-ia),0.)
            weight(ia,ib)= weight(ia,ib)+wt1*wt2
            pval(ia,ib)= pval(ia,ib)+wt1*wt2*slopes(i1,i2)
          end do
        end do
      end do
    end do
    where(weight.ne.0.)
      pval=pval/weight !NORMALIZE BY WEIGHT
end where
    izero=max(nint(-o1/d1)+1,1)
    pval(izero,:)=0.
    weight(izero,:)=1.
!LOOP OVER OUTPUT INTERPOLATING SLOPE BETWEEN REFLECTORS
    do i2=1,n2
      loc1=1
      pos1=pval(1,i2)
      i1=1
      do while (i1<=n1)
        do while (i1<=n1 .and. weight(i1,i2).eq.0.)
          i1= i1+1 !FIND NEXT SLOPE
        end do
        if (i1.ne.n1+1) then
!IF WE DIDN'T REACH THE BOTTOM
          pos2=pval(i1,i2)
          loc2=i1
          dist1=loc2-loc1
          do ia=loc1+1,loc2-1
!INTERPOLATE BETWEEN THE TWO POINTS
            f=real(ia-loc1)/dist1
            pval(ia,i2)=(1.-f)*pos1+f*pos2
          end do
          loc1=loc2
          pos1=pos2 !RESET
        end if
        i1= i1+1
      end do
      pval(loc1+1:,i2)=pval(loc1,i2) !COPY LAST SLOPE DOWN
    end do
    allocate(tt(n1,n2,1))
    tt(:,:,1)=pval
    tempr(1)=7
    call gauss_smooth(tt,3,5,1,tempr(1))
    pval=tt(:,:,1)
    deallocate(slopes,tempr,temp2r,weight)
  end subroutine
  subroutine build_2d_steer_amp_map(fault,amp)
    real :: x,z
    integer :: i1,i2
    real :: min_d,fract
    real :: amp(:,:)
    complex :: fault(:,:)
    do i2=1,n(2)
      x=o(2)+d(2)*(i2-1)
      do i1=1,n(1)
        z=o(1)+d(1)*(i1-1)
        min_d=minval(sqrt((z-aimag(fault(:,:)))**2+(x-real(fault(:,:))&
          &)**2))
        if (min_d > range) then
          amp(i1,i2)=max_amp
        else
          fract=min_d/range
          amp(i1,i2)=max_amp*fract+min_amp*(1.-fract)
        end if
      end do
    end do
  end subroutine
  subroutine create_2d_steer_op(pval,amp,filt)
    real :: pval(:,:),amp(:,:)
    type(nfilter) :: filt
    integer i1,i2,i,nh
    i=0
    allocate(filt%hlx(product(n)+1),filt%pch(product(n)))
    nh=ceiling(hwidth*2.)
    do i2=1,n(2)
      do i1=1,n(1)
        i= i+1
        call allocatehelix(filt%hlx(i),nh)
        call steer_init(hwidth,amp(i1,i2),n(1))
        call steering2(pval(i1,i2),filt%hlx(i))
        filt%pch(i)=i
      end do
    end do
    allocate(filt%hlx(size(filt%hlx))%flt(1))
    allocate(filt%hlx(size(filt%hlx))%lag(1))
    filt%hlx(size(filt%hlx))%lag=0
    filt%hlx(size(filt%hlx))%flt=0.
  end subroutine
end module
