module xcor_mod

implicit none

integer*8,public,save				:: planF,planI

contains

!--------------------------------
subroutine xcor(a,b,n,x)
real 	:: a(n),b(n),x(2*n-1),ai
integer	:: na,nb,nx,i,j,k,l,m,n
na=n
nb=n
nx=2*n-1

x=0
l=0
do i=1,na
   ai=a(i)
   m=nb
   do j=1,nb
      k=l+j
      x(k)=x(k)+ai*b(m)
      m=m-1
   end do
   l=l+1
end do
 
end subroutine
!--------------------------------
subroutine xcor1(a,b,x,n)
real 	:: a(n),b(4*n+1),x(3*n+2),y(3*n+2)
integer	:: i,n,ni,nf

ni=int(n/2)+1
nf=4*n+1-int(n/2)

do i=ni,nf
   x(i-int(n/2))=sum(b(i-int(n/2):i+int(n/2))*a)
end do

end subroutine
!-------------------------------------------------------------------------
subroutine owfftw(tfld)

complex				:: tfld(:)

integer, parameter			:: FFTW_R2HC=0,FTW_HC2R=1,FFTW_DHT=2,FFTW_REDFT00=3,FFTW_REDFT01=4
integer, parameter			:: FFTW_REDFT10=5,FFTW_REDFT11=6,FFTW_RODFT00=7,FFTW_RODFT01=8,FFTW_RODFT10=9
integer, parameter			:: FFTW_RODFT11=10,FFTW_FORWARD=-1,FFTW_BACKWARD=+1,FFTW_DESTROY_INPUT=1
integer, parameter			:: FFTW_UNALIGNED=2,FFTW_CONSERVE_MEMORY=4,FFTW_EXHAUSTIVE=8,FTW_PRESERVE_INPUT=16
integer, parameter			:: FFTW_PATIENT=32,FFTW_ESTIMATE=64,FFTW_ESTIMATE_PATIENT=128
integer, parameter			:: FFTW_BELIEVE_PCOST=256,FFTW_DFT_R2HC_ICKY=512,FFTW_NONTHREADED_ICKY=1024
integer, parameter			:: FFTW_NO_BUFFERING=2048,FFTW_NO_INDIRECT_OP=4096,FFTW_ALLOW_LARGE_GENERIC=8192
integer, parameter			:: FFTW_NO_RANK_SPLITS=16384,FFTW_NO_VRANK_SPLITS=32768,FFTW_NO_VRECURSE=65536

call sfftw_plan_dft(planF,1,size(tfld),tfld,tfld,fftW_FORWARD,fftW_PATIENT)
call sfftw_plan_dft(planI,1,size(tfld),tfld,tfld,fftw_BACKWARD,fftW_PATIENT)

end subroutine
!----------------------------------------------------------------------------
subroutine xcorw(a,b,tfld,n)
real 	            :: a(n),b(n)
complex               :: tfld(n)
complex,allocatable   :: ca(:),cb(:)
integer	            :: n,i

allocate(ca(n),cb(n))
ca=0;cb=0

a=a/maxval(abs(a))
b=b/maxval(abs(b))

tfld=cmplx(a,0)
call sfftw_execute(planF)
ca=tfld
tfld=cmplx(b,0)
call sfftw_execute(planF)
cb=tfld
tfld=ca*conjg(cb)
call sfftw_execute(planI)

cb=tfld
tfld(int(n/2)+1:n)=tfld(1:int(n/2))
do i=1,int(n/2)
  tfld(i)=cb(int(n/2)+i+1)
end do

deallocate(ca,cb)
 
end subroutine

end module
