!<
!semblance
!
!Usage:
!semblance.x  <in.H >out.H
!
!Input Parameters:
!
!
!Output Parameters:
!
!
!
!Description:
!
!
!
!>
!-------------------------------------------------
!
!Author: Robert Clapp, ESMB 463, 7230253
!
!Date Created:Tue Jan 18 21:41:44 PST 2000
!
!Purpose:
!
program semblance 
  use sep3d_struct_mod
  use sep
  implicit none
  integer n1,n2,i1,i2,n3,nkey
  real, allocatable :: input(:,:),semb(:,:)
  type(sep3d) :: in
  type(sep3d) :: out
  integer, external :: getch
  integer  :: gather_axes(10),naxes,kindex(2),param_axis(2)
  character(len=128) ::  type
  logical :: verb
  integer,allocatable:: n(:), fin(:),nin(:),fout(:),nwind(:),jout(:),jin&
    &(:)
  real, allocatable :: o(:),d(:)
  character(len=128),allocatable :: label(:)
  character(len=128),allocatable :: unit(:)
  real  :: opar,dpar,ratio
  real, allocatable :: numer(:,:),denom(:,:),key_vals(:),wt(:)
  integer, allocatable :: count(:,:)
  integer :: npar,ndim_out,nout,nloop(2),smooth,iz,i0,it,begi,endi,ism
  integer :: itr,ipar,ikey,nsz,iaxes
  real :: delta,fract,tau,tempr,numsum,densum,scale
  call initpar()
  call doc('/homes/sep/bob/papers/thesis/Src/main/semblance.rs90')
  call init_sep3d("in",in,"INPUT")
  naxes=getch("gather_axes","d",gather_axes)
  if (naxes .eq. 0) then
    gather_axes=3
    naxes=1
  else
    if (any(gather_axes(1:naxes) < 2) .or. any(gather_axes(1:naxes) >&
      & in%ndims)) then
      write(0,*)  "gather_axes",gather_axes(1:naxes)
      call erexit("invalid gather axes")
    end if
  end if 
  if (1.eq.getch("param_axis","d",param_axis(1))) then
    if (param_axis(1) > 1 .and.  param_axis(1) <= in%ndims) then
      allocate(key_vals(in%n(param_axis(1))))	
      do i1=1,in%n(param_axis(1))
        key_vals(i1)=in%o(param_axis(1))+in%d(param_axis(1))*(i1-1)
      end do
    end if
  else
    call erexit("no param_axis specified")
  end if 
  call from_param("verb",verb,.false.)
  call from_param("type",type,"parab")
  call from_param("opar",opar)
  call from_param("dpar",dpar)
  call from_param("npar",npar)
  call from_param("nout",nout,in%n(1))
  call from_param("smooth",smooth,in%n(1)*2/nout)
  ratio=real(in%n(1))/real(nout)
  ndim_out=2+naxes
  allocate(n(ndim_out),o(ndim_out),d(ndim_out),label(ndim_out),unit&
    &(ndim_out))
!ALLOCATE OUR WINDOWING ROUTINES ASIGN INITIAL VALUES
  allocate(nin(in%ndims),fin(in%ndims))
  allocate(nwind(ndim_out),fout(ndim_out),jout(ndim_out),jin(in%ndims)&
    &)
  jin=1
  jout=1
  fin=0
  nin=in%n
  fout=0
  n(1)=nout
  o(1)=in%o(1)
  d(1)=in%d(1)*real(nout)/real(in%n(1))
  n(2)=npar
  o(2)=opar
  d(2)=dpar
  iaxes=3
  nloop=1
  do i1=1,naxes 
    n(iaxes)=in%n(gather_axes(i1))
    o(iaxes)=in%o(gather_axes(i1))
    d(iaxes)=in%d(gather_axes(i1))
    nwind(iaxes)=1
    nin(gather_axes(i1))=1
!		label(iaxes)=in%label(gather_axes(i1))
!		unit(iaxes)=in%unit(gather_axes(i1))
    nloop(i1)= n(iaxes)
    iaxes=		iaxes+1
  end do 
  if (type(1:5) .eq. "parab") then
    label(2:2-1+nkey)="moveout"
!INITILIZE PARAB
  else
    call erexit("invalid type")
  end if 
  call par_init(out,"OUTPUT" ,"FLOAT","REGULAR",n,o,d)
  where(nwind .ne.1)
    nwind=out%n
end where
  call sep3d_write_description("out",out)
  call hclose()
  nsz=size(key_vals)
  allocate(input(in%n(1),nsz),wt(nsz))
  allocate(numer(in%n(1),npar),denom(in%n(1),npar),count(in%n(1),npar)&
    &,semb(nout,npar))
!LOOP OVER AXES
  do i2=1,nloop(2)
    if (naxes > 1) then
!we need to initilize
      fin(gather_axes(2))=i2-1
      fout(nkey+3)=i2-1
    end if
    do i1=1,nloop(1)
      fin(gather_axes(1))=i1-1
      fout(3)=i1-1
      if (verb) then
        write(0,*) "working on" ,i1,i2," of ",nloop
      end if
!READ IN GATHER
      call sreed_window("in",in%ndims,in%n,nin,fin,jin,4,input)
      do itr=1,nsz
        wt(itr)=sum(abs(input(:,itr)))
      end do 
!#			wt=wt/(sqrt(wt(1))+wt)
!			wt=1/wt
!
!			do itr=1,nsz
!				input(:,itr)=input(:,itr)*wt(itr)
!
!			input(:,23:)=0.
      call srite("check.H",input,size(input)*4)
      semb=0
      numer=0
      count=0
      denom=0
!CORRECT GATHER
      if (type(1:5) .eq. "parab") then
!PARALIZED THIS
!!$OMP PARALLEL private(iz,tau,delta,fract,it,tempr,begi,endi,i0,ism,numsum,densum) default(shared)
!      !$OMP DO
        do 	ipar=1,npar 
          scale=opar+(ipar-1)*dpar
          do itr=1,nsz 
            do iz=1,in%n(1)
              tau=in%o(1)+in%d(1)*(iz-1)
              delta=scale * key_vals(itr)**2
              if (key_vals(itr)<0.) then
                delta=delta*-1.
              end if
              fract=(tau+delta-in%o(1))/in%d(1)+1.
              it=fract
              fract= fract-it
              if (it > 0 .and. it < in%n(1)) then
                tempr=(input(it,itr)*(1-fract)+fract*input(it+1,itr))
              else if (it.eq.in%n(1)) then
                tempr=input(it,itr)
              else
                tempr=0
              end if 
              if (tempr.ne.0) then
                numer(iz,ipar)  =      			     numer(iz,ipar)  + tempr
                denom(iz,ipar)=         			   denom(iz,ipar)+ tempr**2
                count(iz,ipar)=         			   count(iz,ipar)+ 1
              end if
            end do 
            do ism=1,nout 
              numsum=0
              densum=0
              begi=max(nint( ratio*ism -smooth/2.),1)
              endi=min(nint( ratio*ism +smooth/2.),in%n(1))
              do i0=begi,endi     
                numsum=        				numsum+numer(i0,ipar)**2
                densum=        				densum+denom(i0,ipar)*count(i0,ipar&
                  &)
              end do 
              if (densum.eq.0) then
                semb(ism,ipar)=0.
              else
                semb(ism,ipar)=numsum/densum
              end if
            end do
          end do 
!!$OMP END DO
!!$OMP END PARALLEL
        end do 
!WRITE OUT GATHER
        call srite_window("out",out%ndims,out%n,nwind,fout,jout,4,semb&
          &)
      end if
    end do
  end do 
  deallocate(input)
end program