# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!<;
!calc_semblance
!
!Usage:
!call calc_semblance (arg1, arg2 arg3)
!
!Input Parameters:
!
!
!Output Parameters:
!
!
!
!Description:
!
!
!
!>
!-------------------------------------------------
!
!Author: Robert Clapp, ESMB 463, 7230253
!
!Date Created:Tue May 22 11:04:56 PDT 2001
!
!Purpose:
!
!
module calc_semblance_mod
  use sep
  use basic_util_mod
  use parse_region_mod
  use sep3d_struct_mod
  use dip_calc_mod
  implicit none
  real, private, save :: ovel(3),dvel(3),max_err,ogather(2),dgather(2)&
    &,max_rms_err_pct,pct_keep_semb,min_keep_semb
  integer,private,save:: nvel(3),nrms,ngather(2),ix_pos,iy_pos,nsmooth_semb&
    &,n_d_rms_semb
  real,allocatable,save,private :: slow0_sq(:,:,:),var_axis(:),cut(:)&
    &,aprox_off(:)
  real,private,allocatable,save :: min_rms(:),max_rms(:),rms(:),slowsq&
    &(:)
  logical, private, save :: kirchoff
  integer :: ihave
  contains
  subroutine init_semb(ov,dv,velocity,max_r,nrms, cut_field,nmig,omig&
    &,dmig)
    real :: ov(3),dv(3),omig(2),dmig(2),max_r
    real :: velocity(:,:,:)
    integer :: nrms,nang,nmig(2),i
    real :: cut_field(:),z,f
    integer :: iz
    ovel=ov
    dvel=dv
    nvel=(/size(velocity,1),size(velocity,2),size(velocity,3)/)
    allocate(slow0_sq(nvel(1),nvel(2),nvel(3)))
    max_err=max_r/50.
    nrms=nrms
    ngather=nmig
    ogather=omig
    dgather=dmig
    slow0_sq=1./(velocity*velocity)
    allocate(cut(size(cut_field)))
    allocate(var_axis(ngather(2)))
    allocate(aprox_off(ngather(2)))
    do i=1,ngather(2)
      var_axis(i)=ogather(2)+dgather(2)*(i-1)
    end do
    if (.not. kirchoff) then
      var_axis=atan(1.)/45.*var_axis
      cut=cut_field*atan(1.)/45.
    else
      cut=cut_field
    end if
    allocate(min_rms(ngather(1)),max_rms(ngather(1)),rms(ngather(1))&
      &,slowsq(ngather(1)))
  end subroutine
  logical function change_semb_loc(new_x)
    real :: new_x(2)
    real :: z,f
    integer :: iz,iv
    change_semb_loc=.false.
    if (new_x(1) < ovel(2) .or. new_x(1) >ovel(2)+ (nvel(2)-1)*dvel(2)&
      &) then
      write(0,*) "invalid x location",new_x(1),"min_max",ovel(2),ovel(2&
        &)+(nvel(2)-1)*dvel(2)
      return
    end if
    ix_pos=nint(new_x(1)-ovel(2))/dvel(2)+1.
    if (new_x(2) < ovel(3)-dvel(3)/2. .or. new_x(2) >ovel(3)+ (nvel(3)&
      &-.5)*dvel(3)) then
      write(0,*) "ovel",ovel(3),dvel(3)
      write(0,*) "invalid y location",new_x(2),"min_max",ovel(3),ovel(3&
        &)+(nvel(3)-1)*dvel(3)
      return
    end if
    iy_pos=max(1,min(nvel(3),nint((new_x(2)-ovel(3))/dvel(3)+1.)))
!extract the interval velocity at this location and
!map it to gather depth coordinates
    do iz=1,ngather(1)
      z=(iz-1)*dgather(1)+ogather(1)
      f=(z-ovel(1))/dvel(1)+1.
      iv=f
      f= f-iv
      if (iv<1) then
        iv=1
        f=0.
      end if
      if (iv>=nvel(1)) then
        iv=nvel(1)-1
        f=1.
      end if
      slowsq(iz)=(1.-f)*sqrt(1./slow0_sq(iv,ix_pos,iy_pos))+f*sqrt(1.&
        &/slow0_sq(iv+1,ix_pos,iy_pos))
    end do
!calculate the minimum and maximum rms corresponding to a given change
    call calc_rms(slowsq,rms)
    call calc_rms((1.-max_err)*slowsq,min_rms)
    call calc_rms((1.+max_err)*slowsq,max_rms)
    slowsq=1./(slowsq*slowsq)
    change_semb_loc=.true.
  end function
  subroutine sem_clean()
    deallocate(slow0_sq,min_rms,max_rms,rms,slowsq,cut,var_axis,aprox_off&
      &)
  end subroutine
  subroutine calc_rms(intvel,rms_vel)
    real :: intvel(:)
    real :: rms_vel(:)
    integer :: i1
    real :: mys,myt
    rms_vel(1)=intvel(1)
    mys=dgather(1)*intvel(1)*intvel(1)
    myt=dgather(1)
    do i1=2,ngather(1)
      mys= mys+dgather(1)*intvel(i1)*intvel(i1)
      myt= myt+dgather(1)
      rms_vel(i1)=sqrt(mys/myt)
    end do
  end subroutine
  subroutine calc_semblance(location, gather, semblance,min_range,max_range&
    &)
    real :: location(:) , gather(:,:)
    real :: semblance(:)
    real :: min_range,max_range
    integer :: ivel(3),igath,it
    integer :: ipar,bz,ez,i2,ipt,iz
    real :: numer(size(gather,1)),denom(size(gather,1)),sumn(size(gather&
      &,1))
    real :: par,tempr,delta,f,numsum,densum,z,ft
    real :: dsemb,maxoff
    ivel=coord_transform(location)
    igath=ivel(1)
    bz=min(size(gather,1),max(1,igath-nsmooth_semb))
    ez=min(size(gather,1),max(1,igath+nsmooth_semb))
    min_range=(1./(rms(ivel(1))*(1.+max_err))**2)-1./rms(ivel(1))**2
    max_range=(1./(rms(ivel(1))*(1.-max_err))**2)-1./rms(ivel(1))**2
    min_range=-max_range
    dsemb=(max_range-min_range)/(nrms-1)
    if (.not. kirchoff) then
      aprox_off=2.*location(3)*tan(var_axis)
    else
      aprox_off=var_axis
    end if
!write(0,*) "APPROX 2",minval(aprox_off),maxval(aprox_off)
    i2=1
    do while (i2<size(var_axis) .and. var_axis(i2) < cut(ivel(1)))
      i2= i2+1
    end do
!write(0,*) "check range",min_range,max_range,max_rms(ivel(1)),rms(ivel(1)),min_rms(ivel(1))
!
!write(0,*) "MAX DELTA",location(3),sqrt(abs(aprox_off(i2)**2*(/min_range,max_range/)))
!write(0,*) "so slow",bz,ez,size(gather,2),size(semblance)
    do ipar=1,size(semblance)
      numer(bz:ez)=0
      denom(bz:ez)=0
      sumn(bz:ez)=0
      par=(ipar-1)*dsemb+min_range
      do i2=1,size(gather,2)
        do ipt=bz,ez
          if (var_axis(i2) >cut(ipt)) then
            cycle
          end if
          delta=par*aprox_off(i2)*aprox_off(i2)
          f=delta+(ogather(1)+dgather(1)*(ipt-1))**2
! par*off^2+zpos^2
! write(0,*) "check",i2,ipar,delta,cut(ipt),par,sqrt(f)
          if (f>0) then
            f=sqrt(f)
! sqrt(par*off^2+zpos^2)
            f=(f-ogather(1))/dgather(1)+1.
            it=f
! write(0,*) "check2",i2,ipar,f-ipt
            f= f-it
            if (it>0 .and. it < ngather(1)-1) then
              tempr=gather(it,i2)*(1.-f) + gather(it+1,i2)*f
              if (tempr.ne.0.) then
                numer(ipt) = numer(ipt) + tempr
                denom(ipt)= denom(ipt)+ tempr**2
                sumn(ipt)= sumn(ipt)+ 1
              end if
            end if
          end if
        end do
      end do
      numsum=sum(numer(bz:ez)**2)
      densum=sum(denom(bz:ez)*sumn(bz:ez))
      if (densum.eq.0) then
        semblance(ipar)=0
      else
        semblance(ipar)=numsum/densum
      end if
    end do
  end subroutine
  function coord_transform(location) result(ivel)
    integer ivel(3)
    real :: location(3)
    ivel(2)=max(1,min(nvel(2),nint((location(1)-ovel(2))/dvel(2)+1.5))&
      &)
    ivel(3)=max(1,min(nvel(3),nint((location(2)-ovel(3))/dvel(3)+1.5))&
      &)
    ivel(1)=(location(3)-ogather(1))/dgather(1)
    ivel(1)=max(1,(min(ngather(1),ivel(1))))
  end function
  logical function moveout_dt(pos,semb,myang,delta,ang,depth)
    real :: pos(:),semb,myang,delta
    real :: dz,f,posx,posy,posz,myoff,tempr
    integer :: iloc(3)
    logical, optional :: ang,depth
    logical :: ang_l,d_l
    ang_l=.true.
    d_l=.false.
    if (present(ang)) then
      ang_l=ang
    end if
    if (present(depth) ) then
      d_l=depth
    end if
    posx=pos(1)
    if (size(pos).eq.2) then
      posy=posx
      posz=pos(2)
    else
      posy=pos(2)
      posz=pos(3)
    end if
    iloc=coord_transform((/posx,posx,posz/))
    iloc(3)=1
    moveout_dt=.false.
    if (kirchoff .and. .not. ang_l) then
      myoff=myang
    else
      myoff=2.*posz*tan(myang*atan(1.)/45.)
    end if
    delta=0
    if (kirchoff) then
      if (myoff>cut(iloc(1))) then
        return
      end if
    else
      if (myang*atan(1.)/45.>cut(iloc(1))) then
        return
      end if
    end if
!moveout will be in time
!
! dz = semb * 2.*tan(angle)
! dt = dz * slow
!
    f=(posz-ovel(1))/dvel(1)+1.
    if (f<1. ) then
      iloc(1)=1
      f=0.
    else if (f >= nvel(1)) then
      iloc(1)=nvel(1)-1
      f=1
    else
      iloc(1)=f
      f= f-iloc(1)
    end if
    tempr=semb*myoff
    if (tempr.eq.0.) then
      dz=0.
    else
      dz=sqrt(abs(tempr)) * tempr/abs(tempr)
    end if
    delta=(semb*myoff**2+posz**2)
    if (delta<0.) then
      return
    end if
    delta=sqrt(delta)-posz
    if (.not. d_l) then
      delta=sqrt(slow0_sq(iloc(1),iloc(2),iloc(3)))*delta
! delta=delta
    end if
    moveout_dt=.true.
  end function
  subroutine finish_moveout()
    call sep_put_data_axis_par("move_panel",1,ngather(1),ogather(1)&
      &,dgather(1),"Depth")
    call sep_put_data_axis_par("move_panel",2,ngather(2),ogather(2)&
      &,dgather(2),"Moveout")
    call sep_put_data_axis_par("move_panel",3,ihave,1.,1.,"instance")
    call sep_put_data_axis_par("move_line",1,ngather(2),ogather(2),dgather&
      &(2),"Moveout")
    call sep_put_data_axis_par("move_line",2,ihave,1.,1.,"instance")
    call auxputch("esize","d",8,"move_line")
  end subroutine
  subroutine add_moveout(gather,pos,moveout)
    real :: gather(:,:),pos(3),moveout,delta,x,z
    integer :: i
    call srite("move_panel",gather,size(gather)*4)
    do i=1,size(gather,2)
      if (moveout_dt(pos,moveout,ogather(2)+dgather(2)*(i-1),delta,.false.&
        &,.true.)) then
        x=ogather(2)+dgather(2)*(i-1)
        z=pos(3)+delta
      end if
      call srite("move_line",x,4)
      call srite("move_line",z,4)
    end do
    ihave=ihave+1
  end subroutine
  subroutine moveout_check_init()
    ihave=0
  end subroutine
  logical function get_vertical_coherence(vel_s,vel,err,coher,pts_field&
    &)
    type(sep3d) :: mig_s,vel_s,cut_s
    real, pointer :: err(:),coher(:),pts_field(:,:),temp(:,:)
    real, pointer :: gather(:,:),semb_field(:),cutf(:),temp1(:)
    integer :: i3,i2,iout,npts,ix,iz,iy,i4,i1
    integer, external:: auxin,getch
    logical :: is_2d
    real :: min_semb,max_semb,maxme,mypct,dsemb
    real :: vel(:,:,:)
    integer :: mymax(1),nv(3),ic
    integer,external :: sseek_block
    real :: ov(3),dv(3),pos
    logical :: moveout_check
    get_vertical_coherence=.false.
    call set_n_o_d_3d(vel_s,nv,ov,dv)
    allocate(semb_field(3*n_d_rms_semb))
    iout=0
    i3=-1
    i4=-1
    if (1.eq.getch("moveout","d",i3)) then
      moveout_check=.true.
    else
      moveout_check=.false.
    end if
    npts=size(pts_field,2)
    if (-1.eq.auxin("mig")) then
      write(0,*) "must supply migration 'mig' tag"
      return
    end if
    call auxclose("mig")
    call init_sep3d("mig",mig_s,"INPUT")
    allocate(cutf(mig_s%n(1)))
    if (-1.eq.auxin("cut")) then
      cutf=mig_s%o(1)+mig_s%d(1)*(mig_s%n(1)-1)
    else
      call auxclose("cut")
      call init_sep3d("cut",cut_s,"INPUT")
      allocate(temp1(cut_s%n(1)))
      if (.not. sreed_field("cut",temp1)) then
        write(0,*) "trouble reading in cut function"
        return
      end if
      do i1=1,mig_s%n(1)
        pos=(i1-1)*mig_s%d(1)+mig_s%o(1)
        pos=floor((pos-cut_s%o(1))/cut_s%d(1)+1.)
        if (pos < 1) then
          cutf(i1)=temp1(1)
        else if (pos >= cut_s%n(1)) then
          cutf(i1)=temp1(cut_s%n(1))
        else
          iz=pos
          pos= pos-iz
          cutf(i1)=(1.-pos)*temp1(iz)+pos*temp1(iz)
        end if
      end do
      deallocate(temp1)
    end if
    call init_semb(ov,dv,vel,3*max_rms_err_pct,n_d_rms_semb*3,cutf,mig_s%n&
      &(1:2),mig_s%o(1:2),mig_s%d(1:2))
    allocate(gather(mig_s%n(1),mig_s%n(2)))
    call moveout_check_init()
    pct_keep_semb=pct_keep_semb/100.
    is_2d=sep3d_ndims(mig_s).eq.3
    ic=0
    do i2=1,npts
      ix=(pts_field(1,i2)-mig_s%o(3))/mig_s%d(3)+1.5
      if (sep3d_ndims(mig_s).eq.4) then
        iy=(pts_field(2,i2)-mig_s%o(4))/mig_s%d(4)+1.5
      end if
      iy=1
      iz=(pts_field(3,i2)-mig_s%o(1))/mig_s%d(1)+1.5
      if (is_2d) then
        iy=1
      end if
      if (i3.ne.ix .or. i4.ne. iy) then
        i1=ix+(iy-1)*mig_s%n(3)
        if (i1.ne. sseek_block("mig",i1,size(gather)*4,0)) then
          write(0,*) "trouble seeking tag=mig to ",i1
          return
        end if
        if (.not. sreed_field("mig",gather)) then
          write(0,*) "trouble reading in migration block",ix,iy
          return
        end if
        if (.not. change_semb_loc(pts_field(1:2,i2))) then
          write(0,*) "Trouble changing semblance location"
          return
        end if
        i3=ix
        i4=iy
      end if
      call calc_semblance(pts_field(:,i2),gather,semb_field,min_semb&
        &,max_semb)
      dsemb=(max_semb-min_semb)/n_d_rms_semb/3
      maxme=maxval(semb_field)
      mymax=maxloc(semb_field)
      call calc_semblance(pts_field(:,i2),gather,semb_field,min_semb&
        &,max_semb)
      dsemb=(max_semb-min_semb)/n_d_rms_semb/3
      maxme=maxval(semb_field)
      mymax=maxloc(semb_field)
      mypct=(maxme-sum(semb_field)/size(semb_field))/(sum(semb_field)&
        &/size(semb_field))
      if (mypct > pct_keep_semb .and. maxme>min_keep_semb .and. mymax&
        &(1)>2.*size(semb_field)/6 .and. mymax(1)<4.*size(semb_field)/6) then
        ic= ic+1
        err(ic)=(mymax(1)-1)*dsemb+min_semb
        coher(ic)=maxme
        call add_moveout(gather,pts_field(:,i2), coher(ic))
        pts_field(:,ic)=pts_field(:,i2)
      end if
    end do
    call finish_moveout()
!if(.not.finish_semb_wt_moveout(nout,.true.)){
! write(0,*) "trouble finishing semb_wt_moveout"
! return
!}
    allocate(temp(3,ic))
    temp=pts_field(:,1:ic)
    deallocate(pts_field)
    allocate(pts_field(3,ic))
    pts_field=temp
    deallocate(temp)
    call sep3df_delete(mig_s,"mig")
    deallocate(semb_field,gather)
    get_vertical_coherence=.true.
  end function
  logical function check_pars_calc_semb(file,norerun)
    character(len=*) :: file
    logical :: stat_ok(10),stat_rerun(10),norerun
    check_pars_calc_semb=.false.
    stat_ok=.true.
    stat_rerun=.true.
    stat_ok(1)=check_pars_parse_region(file,stat_rerun(1))
    stat_ok(2)=check_pars_dip_calc(file,stat_rerun(2))
    stat_rerun(3)= check_par_file("n_d_rms_semb" ,n_d_rms_semb &
      &,file,101)
    stat_rerun(4)= check_par_file("min_keep_semb" ,min_keep_semb &
      &,file,.1)
    stat_rerun(5)= check_par_file("max_rms_err_pct" ,max_rms_err_pct&
      & ,file,300.)
    stat_rerun(6)= check_par_file("nsmooth_semb" ,nsmooth_semb &
      &,file,4)
    stat_rerun(7)= check_par_file("pct_keep_semb" ,pct_keep_semb &
      &,file,300.)
    stat_rerun(8)= check_par_file("kirchoff" ,kirchoff ,file,.false.&
      &)
    if (all(stat_ok)) then
      check_pars_calc_semb=.true.
    end if
    if (all(stat_rerun)) then
      norerun=.true.
    else
      norerun=.false.
    end if
  end function
end module
