# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!<
!scan_it
!
!Usage:
!call scan_it (arg1, arg2 arg3)
!
!Input Parameters:
!
!
!Output Parameters:
!
!
!
!Description:
!
!
!
!>
!-------------------------------------------------
!
!Author: Robert Clapp, ESMB 463, 7230253
!
!Date Created:Thu Jan 30 10:58:22 2003
!
!Purpose:
!
!
module scan_mod
  use sep3d_struct_mod
  use sep
  use sep3d_struct_mod
  use put_3d_mod
  use map_pts_reg_mod
  implicit none
  real, pointer,private,save :: tantheta(:),numer(:),denom(:),count(:)&
    &,cut(:,:,:),myc(:),dips(:,:,:),maps(:,:,:)
  integer,private,save :: n(4),ngam,ndip
  integer, private, pointer, save :: idips(:)
  real,private,save :: o(4),d(4),ogam,dgam,odip,ddip
  integer, private,save :: nwind1,nwind2,nwind3
  type(sep3d), private, save :: cut_s,dip_s,cut_f,dip_f
  logical :: no_cut,no_dip
  contains
  logical function init_scan(mig_s)
    type(sep3d) :: vel_s,mig_s
    integer :: ni(4),i1,na(2)
    real :: di(4),x,y,z,ocut(3),dcut(3)
    real, pointer :: cquad(:,:)
    integer :: ncut(3),i
    integer, external :: auxin
    init_scan=.false.
    n=1
    o=0
    d=1.
    if (sep3d_ndims(mig_s).eq.3) then
      n(1:3)=mig_s%n(1:3)
      o(1:3)=mig_s%o(1:3)
      d(1:3)=mig_s%d(1:3)
    else if (sep3d_ndims(mig_s).eq.4) then
      n=mig_s%n
      o=mig_s%o
      d=mig_s%d
    end if
    allocate(tantheta(n(2)))
    do i1=1,n(2)
      tantheta(i1)=tan(atan(1.)/45.*(o(2)+d(2)*(i1-1)))**2
    end do
    allocate(numer(n(1)),denom(n(1)),count(n(1)))
    call from_param("nwind1",nwind1,min(5,n(1)))
    call from_param("nwind2",nwind2,min(5,n(3)))
    call from_param("nwind3",nwind3,min(5,n(4)))
    call from_param("ocut1",ocut(1),o(1))
    call from_param("ocut2",ocut(2),o(3))
    call from_param("ocut3",ocut(3),o(4))
    call from_param("dcut1",dcut(1),10.*d(1))
    call from_param("dcut2",dcut(2),10.*d(3))
    call from_param("dcut3",dcut(3),10.*d(4))
    call from_param("ncut1",ncut(1),min(n(1),nint( (o(1)+d(1)*(n(1)-1.&
      &)-ocut(1))/dcut(1)+1.5)))
    call from_param("ncut2",ncut(2),min(n(3),nint( (o(3)+d(3)*(n(3)-1.&
      &)-ocut(2))/dcut(2)+1.5)))
    call from_param("ncut3",ncut(3),min(n(4),nint( (o(4)+d(4)*(n(4)-1.&
      &)-ocut(3))/dcut(3)+1.5)))
    allocate(cut(ncut(1),ncut(2),ncut(3)))
    call init_sep3d(cut_s,"SCRATCH","FLOAT","REGULAR",ncut,ocut,dcut)
    allocate(cquad(4,1))
    cquad(1:3,1)=ocut(1:3)+dcut(1:3)*.5 * ncut(1:3)
    cquad(4,1) =mig_s%o(2)+mig_s%d(2)*(mig_s%n(2)-1.)
    no_cut=.true.
    call interp_ireg_to_reg(cut_s,cquad(1:3,:),cquad(4,:),cut,(/1.,1.&
      &,1./))
    if (-1.ne.auxin("dips")) then
      call init_sep3d("dips",dip_s,"INPUT")
      if (dip_s%n(1).ne.n(1)) then
        call seperr("depth axis of dips and semblance must be the same"&
          &)
      end if
      if (sep3d_ndims(dip_s).eq.2) then
        allocate(dips(dip_s%n(1),dip_s%n(2),1))
      else
        allocate(dips(dip_s%n(1),dip_s%n(2),dip_s%n(3)))
      end if
      call sep3d_grab_headers("dips",dip_s,i)
      if (.not. sep3d_read_data("dips",dip_s,dips)) then
        call seperr("trouble reading dips")
      end if
      dips=atan(dips)
      no_dip=.false.
    else
      no_dip=.true.
    end if
    if (.not. make_map()) then
      write(0,*) "trouble making map"
      return
    end if
    allocate(idips(n(1)))
    idips=(0.-odip)/ddip+1.5
    call init_coord(cut_s%n(1:3),cut_s%o(1:3),cut_s%d(1:3))
    allocate(myc(mig_s%n(1)))
    init_scan=.true.
  end function
  logical function make_map()
    real :: mymin,mymax
    integer :: igam,idip,itheta
    real :: gamma,mydip,theta
    call from_param("ogam",ogam)
    call from_param("dgam",dgam)
    call from_param("ngam",ngam)
    call from_param("maxdip",mymax,60.)
    call from_param("ddip",ddip,.5)
    ndip=mymax/ddip*2+1
    odip=-mymax*atan(1.)/45.
    ddip=ddip*atan(1.)/45.
    allocate(maps(ndip,n(2),ngam))
    do igam=1,ngam
      gamma=(igam-1)*dgam+ogam
      do itheta=1,n(2)
        theta=((itheta-1)*d(2)+o(2))*atan(1.)/45.
        do idip=1,ndip
          mydip=(idip-1)*ddip+odip
! mydip=atan(mydip)
! mydip=0.
          maps(idip,itheta,igam)= (1-gamma)/cos(mydip)*(sin(theta)**2&
            &)/ (cos(mydip)**2-sin(theta)**2)
! maps(idip,itheta,igam)=
! (1-gamma)*
! (sin(theta)**2)/
! (cos(mydip))**2-sin(theta)**2
! maps(idip,itheta,igam)=
! (1.-gamma)*tan(theta)**2
!if(abs(mydip) < .001)
! fz=(sc*tantheta(i2)*z-o(1))/d(1)+1.
        end do
      end do
    end do
! fz=(sc*tantheta(i2)*z-o(1))/d(1)+1.+i1
    make_map=.true.
  end function
  logical function grab_cut(pos)
    real :: pos(2),myp(3)
    integer :: i
    myp(2:3)=pos(1:2)
    if (.not. no_cut) then
      do i=1,n(1)
        myp(1)=o(1)+d(1)*(i-1)
        myc(i)=get_val_lin_pos(cut,myp(1),myp(2),myp(3))
      end do
    else
      myc=o(2)+d(2)*(n(2)-1)
    end if
    grab_cut=.true.
  end function
  logical function grab_dip(pos)
    real :: pos(:)
    integer :: i2,i3,i1
    grab_dip=.true.
    if (no_dip) then
      return
    end if
    i3=1
    i2=(pos(1)-dip_s%o(2))/dip_s%d(2)+1.5
    if (sep3d_ndims(dip_s).eq.3) then
      i3=(pos(2)-dip_s%o(3))/dip_s%d(3)+1.5
    end if
    if (i2< 1 .or. i2 > size(dips,2) .or. i3<1 .or. i3>size(dips,3))&
      & then
      write(0,*) "pos",pos
      write(0,*) "position out of range of dips",i2,i3
      grab_dip=.false.
      return
    end if
    do i1=1,size(dips,1)
      idips(i1)=max(1,min(ndip,nint((dips(i1,i2,i3)-odip)/ddip+1.)))
    end do
!idips=(0.-odip)/ddip+1.5
!idips=1
  end function
  logical function compute_semb(pos,gamma,data,semb)
    real :: data(:,:),semb(:),gamma
    integer ::ix,iy,iz,i1,i2,ia,idip,igamma
    real :: fz,sc,amp,densum,numsum,pos(2),a,z,b
    compute_semb=.false.
    if (.not. grab_cut(pos)) then
      write(0,*) "POSITION",pos
      call seperr("trouble grabbing cut function")
    end if
    if (.not. grab_dip(pos)) then
      write(0,*) "POSITION",pos
      call seperr("trouble grabbing dip function")
    end if
    igamma=min(ngam,max(1,nint((gamma-ogam)/dgam+1.)))
    semb=0.
    sc=1.-gamma
    numer=0
    denom=0
    count=0
    ia=(5-o(2))/d(2)+1.
    do i2=1,n(2)
      a=o(2)+d(2)*(i2-1)
      a=a/45.*atan(1.)
      z=o(1)-d(1)
      do i1=1,n(1)
        z= z+d(1)
        if (myc(i1)>=a) then
!cut conditional would be here
          fz=(sc*tantheta(i2)*z-o(1))/d(1)+i1
!if(i1==100 )write(0,*) "BBB",a,z,gamma,fz,sc*tantheta(i2)*z
! fz=((1-gamma)/(1-gamma*(1-cos(dip(i1))))*(sin(a)**2)/
! (cos(dip(i1))**2-sin(a)**2)*z-o(1))/d(1)
          iz=fz
          fz= fz-iz
          if (iz > 0 .and. iz < n(1)) then
            amp=(1.-fz)*data(iz,i2)+fz*data(iz+1,i2)
            if (amp>0) then
              numer(i1)= numer(i1)+amp
              denom(i1)= denom(i1)+amp**2
              count(i1)= count(i1)+1
            end if
          else
            count(i1)= count(i1)+1
          end if
        end if
      end do
    end do
! call srite("denom.H",denom*count,size(denom)*4)
! call srite("numer.H",numer*numer,size(numer)*4)
! nwind1=0
    do i1=1,n(1)
      numsum=0
      densum=0
      do i2=max(1,i1-nwind1),min(i1+nwind1,n(1))
        if (count(i2) > 0.) then
          numsum= numsum+numer(i2)**2
          densum= densum+denom(i2)*count(i2)
        end if
      end do
      if (densum>0.) then
        semb(i1)=numsum/densum
      end if
    end do
    compute_semb=.true.
  end function
end module
