#<
#semblance
#
#Usage:
#semblance.x  <in.H >out.H 
#
#Input Parameters:
#
#
#Output Parameters:
#
#
#
#Description:
#
#
#
#>
#%end of self-documentaiton
#-------------------------------------------------
#
#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
  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 init_sep3d("in",in,"INPUT")
	naxes=getch("gather_axes","d",gather_axes)
	if(naxes == 0) { gather_axes=3; naxes=1;}
	else{
		if(any(gather_axes(1:naxes) < 2) || any(gather_axes(1:naxes) > in%ndims)){
			write(0,*)  "gather_axes",gather_axes(1:naxes)
			call erexit("invalid gather axes")
		}
	}

	if(1==getch("param_axis","d",param_axis(1))){
		if(param_axis(1) > 1 &&  param_axis(1) <= in%ndims) {
			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)
		}
	}
	else call erexit("no param_axis specified")

	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+=1
	}

		

	if(type(1:5) == "parab"){
		label(2:2-1+nkey)="moveout"
		#INITILIZE PARAB
	}
	else call erexit("invalid type")


	call par_init(out,"OUTPUT" ,"FLOAT","REGULAR",n,o,d)

	where(nwind !=1) nwind=out%n

	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){ #we need to initilize
			fin(gather_axes(2))=i2-1
			fout(nkey+3)=i2-1
		}
		do i1=1,nloop(1){
			
			fin(gather_axes(1))=i1-1
			fout(3)=i1-1

			if(verb) write(0,*) "working on" ,i1,i2," of ",nloop


			#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)))
			
##			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) == "parab"){
				#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.)  delta=delta*-1.
								fract=(tau+delta-in%o(1))/in%d(1)+1.
								it=fract; fract-=it
								if(it > 0 && it < in%n(1)){
            			tempr=(input(it,itr)*(1-fract)+fract*input(it+1,itr))
								}
								else if(it==in%n(1)) tempr=input(it,itr)
								else tempr=0
          			if (tempr!=0){
      			     numer(iz,ipar)  += tempr
         			   denom(iz,ipar)+= tempr**2
         			   count(iz,ipar)+= 1
							}
          	}
			    	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+=numer(i0,ipar)**2
        				densum+=denom(i0,ipar)*count(i0,ipar)
      				}
      				if(densum==0)  semb(ism,ipar)=0.
      				else semb(ism,ipar)=numsum/densum
						}
					}
!!$OMP END DO
!!$OMP END PARALLEL

				}

				#WRITE OUT GATHER
				call srite_window("out",out%ndims,out%n,nwind,fout,jout,4,semb)


			}
		}
	}

				deallocate(input)
		




}