program FOC2
  !! Compute a differential focusing metric
  !! Paul Sava (paul@sep.stanford.edu)
  !! 05/28/02

  !! input:  z-x (zero offset section)
  !! output: z-x (focusing semblance)

  use sep
  implicit none

  real, dimension(:,:), allocatable :: data
  real, dimension(:,:), allocatable :: fmap

  integer :: nz,  nx
  integer :: iz,  ix
  integer :: jz,  jx
  real    :: oz,  ox
  real    :: dz,  dx

  integer :: n3, i3

  integer :: xwin,zwin !! total window = 2*n+1
  real    :: sum1, sum2, stmp, srat

  integer :: style

  call sep_init()

  call from_history("n1",nz)
  call from_history("o1",oz,0.)
  call from_history("d1",dz,1.)

  call from_history("n2",nx)
  call from_history("o2",ox,0.)
  call from_history("d2",dx,1.)

  call from_history("n3",n3,1)

  call from_param("xwin",xwin,1)
  call from_param("zwin",zwin,1)

  call from_param("style",style,1)

  call sep_close()

  allocate(data(nz,nx)); data=0
  allocate(fmap(nz,nx)); fmap=0.

  do i3=1,n3
     write(0,*) i3,n3

     call sep_read(data)
     data = data/maxval(data)

     select case(style)
     case(1)
        do ix=1+xwin,nx-xwin-1
           do iz=1+zwin,nz-zwin-1
              sum1 = 0
              sum2 = 0
              do jx=ix-xwin,ix+xwin
                 do jz=iz-zwin,iz+zwin
                    sum1=sum1 + sqrt ( &
                    (data(jz  ,jx+1) - data(jz,jx))**2 + &
                    (data(jz+1,jx  ) - data(jz,jx))**2)
                    sum2=sum2 + data(jz,jx)**2
                 end do
              end do
              fmap(iz,ix) = sum1/sum2
           end do
        end do
     case(2)
        do ix=1+1,nx-1
           do iz=1,nz
              fmap(iz,ix) = data(iz,ix+1) - data(iz,ix-1)
           end do
        end do
     end select

     call sep_write(fmap)

  end do

  call exit(0)
end program FOC2













