# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!<
!create3d_nfilter
!
!Usage:
!nfilter_3d.x <in.H >out.H [pars]
!
!
program create_nfilter3d
  use sep
  use nbound_bob_mod
  use interp_filter_mod
  use io_nfilt
  implicit none
  integer :: n(3),j(3),center(3),na(3),ncoef,n123,niter,iloc1,iloc2&
    &,nold(3)
  integer :: nuse(3),n_t(3),n_tau,use_faults(100)
  integer :: beg1,end1,beg2,end2,nbig(3),nc,i1,i2,ntau,npts,nuse_ref&
    &,beg3,end3,i3
  real :: obig(3),dbig(3),d_tau
  integer :: nbig123
  real :: o(3),d(3),hwidth,eps,o_t(3),d_t(3),ouse(3),duse(3),min_amp
  real :: max_amp,range,oold(3),dold(3)
  character(len=128) :: label,tag_in,tag_ref
  real, pointer :: data(:),covar(:),temp(:,:,:),model(:)
  real, pointer :: map(:,:,:),velocity(:,:,:),tempv(:,:,:),dip_xz(:,:,:&
    &),dip_yz(:,:,:)
  logical, pointer :: mask(:),msk(:)
  logical :: test_it,rf,is_2d
  type(nfilter) :: filt1,filt2
  real :: scale_lap
  logical :: tau,isotropic,non_stationary,verb
  integer :: use_refs(100),ns,nwind,nref(3),ierr
  real :: oref(3),dref(3)
  real :: dot1(2),dot2(2)
  integer,external :: getch,auxin
  real, allocatable :: refs(:,:,:),j1(:)
  real, pointer :: field(:,:,:),tt(:,:,:),py(:,:,:)
  real :: x,z
  logical :: have_faults,from_initial
  integer :: n1f,n2f,n1r,n2r,n1rn,n2rn,i,smooth_wind
  integer :: noff,nparts,ifilt,ilag,nin,nout,bin,bout
  real :: diprange,angrange,ddip,daz
!read in parameters
  call sep_get_data_axis_par("in",1,n(1),o(1),d(1),label)
  call sep_get_data_axis_par("in",2,n(2),o(2),d(2),label)
  call sep_get_data_axis_par("in",3,n(3),o(3),d(3),label)
  call from_param("verb",verb,.false.)
  nbig=n
  obig=o
  dbig=d
  nuse=nbig
  ouse=obig
  duse=dbig
  nbig123=product(nuse)
  n123=product(n)
  if (verb) then
    write(0,*) "padding axis 1",beg1,end1," axis2 ",beg2,end2," axis3 "&
      &,beg3,end3
  end if
  is_2d=.false.
  if (nuse(3).eq.1) then
    is_2d=.true.
  end if
  call from_param("read_filt",rf,.false.)
  allocate(dip_xz(n(1),n(2),n(3)))
  call sreed("in",dip_xz,size(dip_xz)*4)
  dip_xz=d(2)/d(1)*dip_xz!*.7
  dip_yz=d(3)/d(1)*dip_yz
  call from_param("noff",noff,7)
  call from_param("diprange",diprange,30.)
  call from_param("angrange",angrange,30.)
  call from_param('nparts',nparts,100)
  call from_param("ddip",ddip,1.)
  call from_param("daz",daz,1.)
  if (.not. init_filt_interp(nbig(1),noff,diprange,angrange,nbig(2)&
    &,nparts,ddip,daz,.true.)) then
    call seperr("trouble initializing filter interpolation")
  end if
  if (.not. create_steer_integrate_op_2d(filt1,dip_xz)) then
    call seperr("trouble creating filter")
  end if
  call nbound_bob(nbig,filt1)
  deallocate(dip_xz)
  if (.not. rite_nfilt("out",filt1,nbig,o, d) ) then
    call seperr("trouble riting filter")
  end if
end program
