! module for SU subroutines converted to Fortran
module auxsu
	use constants

	contains
	!-------------------------------------------------------------------
	subroutine convsu (lx,ifx,x,ly,ify,y,lz,ifz,z)
	!-------------------------------------------------------------------
	! convolution subroutine with arbitrary first lag in inputs and 
	! output
	!-------------------------------------------------------------------

		! input parameters
		integer, intent(in)	::	lx	! length of input 1-D array x
		integer, intent(in)	::	ifx	! index of first sample of input 1-D array x
		integer, intent(in)	::	ly	! length of input 1-D array y
		integer, intent(in)	::	ify	! index of first sample of input 1-D array y
		integer, intent(in)	::	lz	! length of output 1-D array z=x*y
		real,dimension(lx),intent(inout)	:: x	! first input array 
		real,dimension(ly),intent(inout)	:: y	! second input array 
	
		! output
		real,dimension(lz),intent(inout)	:: z	! output array z=x*y

		! internal variables
		integer				:: ilx	! index of last sample in x
		integer				:: ily	! index of last sample in y
		integer				:: ilz	! index of last sample in z
		integer				:: jlow,jhigh,i,j	! aux variables
		real				:: sum

		! compute index of last samples
		ilx=ifx+lx-1
		ily=ify+ly-1
		ilz=ifz+lz-1
		!write(0,*)'In conv: ilx=',ilx,' ily=',ily,' ilz=',ilz

		! do deconvolution
		do i=ifz,ilz+1
			jlow=i-ily+1
			if (jlow<ifx) jlow=ifx
			jhigh=i-ify+1
			if (jhigh>ilx) jhigh=ilx
			!write(0,*)'i=',i,' jlow=',jlow,' jhigh=',jhigh
			
			sum=0.0
			do j=jlow,jhigh
				sum=sum+x(j-ifx+1)*y(i-j-ify+2)
			!write(0,*)'i=',i,' j=',j,' x=',x(j-ifx+1),' y=',y(i-j-ify+2),' sum=',sum
			end do
			z(i-ifz+1)=sum
		end do
	end subroutine convsu
			

	!-----------------------------------------------------------------------
	subroutine rickersu (fpeak,dt,lw,iw,wavelet)
	!-----------------------------------------------------------------------
	! make a Ricker wavelet
	!-----------------------------------------------------------------------

		! input parameters
		real,intent(in)		:: fpeak	! peak frequency (Hz)
		real,intent(in)		:: dt		! time sampling interval (s)

		! output parameters
		integer,intent(out)	:: lw		! length of Ricker wavelet
		integer,intent(out)	:: iw		! index of first Ricker wavelet sample
		real,dimension(:),pointer :: wavelet	! Ricker wavelet

		! internal variables
		integer				:: it,jt	! loop counters 
		real				:: t,x		! auxiliary variables
		
		! compute first index and length of Ricker wavelet
		iw=INT(-(1+1.0/(fpeak*dt)))+1	
		lw=1-2*iw
		write(0,*)'in Ricker: fpeak=',fpeak,' dt=',dt,' iw=',iw,' lw=',lw	

		! allocate wavelet
		allocate(wavelet(lw))

		! compute wavelet
		it=iw-1
		t=(it+1)*dt
		do jt=1,lw
			x=PI*fpeak*t
			x=x*x
			wavelet(jt)=EXP(-x)*(1.0-2.0*x)
!write(0,*)'it=',it,' wavelet(',jt,')=',wavelet(jt)
			it=it+1
			t=t+dt
		end do
	end subroutine rickersu

	!---------------------------------------------------------------------
	subroutine ints8r(nxin,dxin,fxin,yin,yinl,yinr,nxout,xout,yout)
	!---------------------------------------------------------------------
	! Interpoaltion of a uniformly-sampled real function y(x) via a 
	! table of 8-coefficient sinc approximations
	!---------------------------------------------------------------------

		! input parameters
		integer,intent(in)	:: nxin	! number of samples in input array y(x)	
		real,intent(in)		:: dxin	! sampling interval in input array y(x)	
		real,intent(in)		:: fxin	! x-value of first sample in input array 
		real,dimension(nxin),intent(in)	:: yin	! values at which the y array is input
		real,intent(in)		:: yinl	! value used to extrapolate y to the left of y(1)
		real,intent(in)		:: yinr	! value used to extrapolate y to the right of y(nxin)
		integer,intent(in)	:: nxout ! number of x-values in output array
		real,dimension(nxout),intent(in) :: xout ! x-values at which y is output 
		
		! output
		real,dimension(nxout),intent(out) :: yout ! array (nxou) of sinc-interpolated y

		! internal variables
		real,dimension(NTABLE,LTABLE),save	:: table	! table of sinc coeffcients
		logical,save		:: tabled=.FALSE.
		integer				:: jtable
		real				:: frac
		integer				:: i,j

		! tabulate sinc interpolation coefficients if not already tabulated
		if (.NOT. tabled) then
			do jtable=2,NTABLE-1
				frac=REAL(jtable-1)/REAL(NTABLE-1)
				call mksinc(frac,LTABLE,table(jtable,:))
!write(0,*)'frac=',frac,' table(',jtable,',1)=',table(jtable,1)
			end do
			do jtable=1,LTABLE
				table(1,jtable)=0.0
				table(NTABLE,jtable)=0.0
			end do
			table(1,LTABLE/2)=1.0
			table(NTABLE,LTABLE/2+1)=1.0
			tabled=.TRUE.
		end if

!do i=1,NTABLE
!	do j=1,8
!	write(0,*)'table(',i,',',j,')=',table(i,j)
!	end do
!end do

		! interpolate using tabulated coefficients
		call intt8r(NTABLE,table,nxin,dxin,fxin,yin,yinl,yinr,nxout,xout,yout)

	end subroutine ints8r

	!--------------------------------------------------------------------------
	subroutine intt8r (mtable,table,nxin,dxin,fxin,yin,yinl,yinr,nxout,xout,yout)
	!--------------------------------------------------------------------------
	! Interpoaltion of a uniformly-sampled real function y(x) via table of
	! 8-coefficient interpolators
	!--------------------------------------------------------------------------
	
		!input parameters
		integer,intent(in)	:: mtable	! table length
		real,dimension(NTABLE,LTABLE),intent(in)	:: table ! of sinc coeffs
		integer,intent(in)	:: nxin	! number of x-values in input y(x)
		real,intent(in)		:: dxin	! sampling interval in input array 	
		real,intent(in)		:: fxin	! first sample in input array 
		real,dimension(nxin),intent(in)	:: yin	! values at which the y array is input
		real,intent(in)		:: yinl	! value used to extrapolate y to the left
		real,intent(in)		:: yinr	! value used to extrapolate y to the left
		integer,intent(in)	:: nxout ! number of samples in output array
		real,dimension(nxout),intent(in) :: xout ! x-values at which y is output 
		
		! output
		real,dimension(nxout),intent(out) :: yout ! sinc-interpolated y

		! internal variables
		integer		:: ioutb,nxinm8,ixout,ixoutn,kyin,ktable,itable
		real		:: xoutb,xouts,xoutn,frac,fntablem1,yini,sum
		integer		:: i,j

		! compute constants
		ioutb=-3-8
		xoutf=fxin
		xouts=1.0/dxin
		xoutb=8.0-xoutf*xouts
		fntablem1=REAL(ntable-1)
		nxinm8=nxin-8
!write(0,*)'ntable=',ntable,' xoutf=',xoutf,' xouts=',xouts,' xoutb=',xoutb, &
!	&' fntablem1=',fntablem1,' nxinm8=',nxinm8

		! loop over output samples
		do ixout=1,nxout
			! determine pointers into the table
			xoutn=xoutb+xout(ixout)*xouts
			ixoutn=INT(xoutn)
			kyin=ioutb+ixoutn+1
			frac=xoutn-REAL(ixoutn)
			if (frac>0.0) then
				ktable=frac*fntablem1+1.5
			else
				ktable=(frac+1.)*fntablem1+0.5
			end if
!write(0,*)'ixout=',ixout,' xoutn=',xoutn,' ixoutn=',ixoutn,' kyin=',kyin,&
!&' frac=',frac,' ktable=',ktable

			! if totally within input array, use fast method
			if (kyin>=0 .AND. kyin <=nxinm8) then
				yout(ixout)= 						&
					yin(kyin+1)*table(ktable,1)+	& 
					yin(kyin+2)*table(ktable,2)+	& 
					yin(kyin+3)*table(ktable,3)+	& 
					yin(kyin+4)*table(ktable,4)+	& 
					yin(kyin+5)*table(ktable,5)+	& 
					yin(kyin+6)*table(ktable,6)+	& 
					yin(kyin+7)*table(ktable,7)+	& 
					yin(kyin+8)*table(ktable,8)

!do i=1,8
!write(0,*)'i=',i,' yin=',yin,' table=',table(i,ktable)
!end do
			else
				! sum over 8 tabulated coeffcients
				sum=0.0
				do itable=1,8
					if (kyin<1) then
						yini=yinl
					else if (kyin>nxin) then
						yini=yinr
					else
						yini=yin(kyin)
					end if
!write(0,*)'itable=',itable,' kyin=',kyin,' yini=',yini
					sum=sum+yini*table(ktable,itable)
!write(0,*)'itable=',itable,' table=',table(ktable,itable),' sum=',sum
					kyin=kyin+1
				end do		
				yout(ixout)=sum
			end if
		end do
	end subroutine intt8r

	!----------------------------------------------------------------------
	subroutine mksinc(d,lsinc,sinc)
	!----------------------------------------------------------------------
	! compute least-squares optimal sinc interpolation coefficients
	!----------------------------------------------------------------------
	
		! input parameters
		real,intent(in)		:: d		! fractional dustance to interpolation point
		integer,intent(in)	:: lsinc	! length of sinc approximation lsinc%2=0
										! and lsinc<=20
		
		! output
		real,dimension(lsinc),intent(out)	:: sinc ! array(lsinc) of interpolation
													! coefficients

		! internal variables
		integer				:: j
		real(kind=DOUBLE),dimension(20) :: s,a,c,work
		real(kind=DOUBLE)				:: fmax

		! compute auto-correlation and cros-correlation arrays
		fmax=0.066+0.265*LOG(REAL(lsinc,kind=DOUBLE))
		if (fmax>=1.0) fmax=1.0
!write(0,*)'fmax=',fmax
		do j=1,lsinc
			a(j)=dsinc(fmax*(j-1))
			c(j)=dsinc(fmax*(lsinc/2-j+d))
!write(0,*)'j=',j,' a=',a(j),' c=',c(j)
		end do

		! solve symmetric Toeplitz system for the sinc approximation
		call stoepd(lsinc,a,c,s,work)
		do j=1,lsinc
			sinc(j)=s(j)
!write(0,*)' j=',j,'sinc=',sinc(j)
		end do
	end subroutine mksinc

	!----------------------------------------------------------------------
	function dsinc(x)
	!----------------------------------------------------------------------
	! evaluate double precission sinc(x) for double precission x
	!----------------------------------------------------------------------
	
		! input
		real(kind=DOUBLE),intent(in)	:: x	
	
		! output
		real(kind=DOUBLE)				:: dsinc

		! internal variables
		real(kind=DOUBLE)				:: pix

		if (x .EQ. 0.0) then
			dsinc=1.0
		else
			pix=PI*x
			dsinc=sin(pix)/pix
		end if
	end function dsinc

	!----------------------------------------------------------------------
	subroutine stoepd(n,r,g,f,a)
	!----------------------------------------------------------------------
	! solve a symmetric Toeplitz linear system of equations Rf=g for f
	! with double precission arithmetic
	!----------------------------------------------------------------------
	
		! input
		integer,intent(in)							:: n ! dimension of system
		real(kind=DOUBLE),dimension(n),intent(in)	:: r ! top row of Toeplitz matrix
		real(kind=DOUBLE),dimension(n),intent(in)	:: g ! RHS column vector

		! output
		real(kind=DOUBLE),dimension(n),intent(out)	:: f ! solution: LHS column vector
		real(kind=DOUBLE),dimension(n),intent(out)	:: a ! Solution to Ra=v		

		! internal variables
		integer							:: i,j
		real(kind=DOUBLE)				:: v,e,c,w,bot
			
		if (r(1) .EQ. 0.0) return
		a(1)=1.0
		v=r(1)
		f(1)=g(1)/r(1)
	
		do j=2,n
			! solve Ra=v as in Claerbout, FGDP, p. 57
			a(j)=0.0
			f(j)=0.0
			e=0.0
			do i=1,j
				e=e+a(i)*r(j-i+1)
			end do
			c=e/v
			v=v-c*e
!write(0,*)'j=',j,' e=',e,' c=',c,' v=',v

			do i=1,(j+1)/2
				bot=a(j-i+1)-c*a(i)
				a(i)=a(i)-c*a(j-i+1)
				a(j-i+1)=bot
			end do

			! use a and v above to get f(i), i=1,2,...j+1
			w=0.0
			do i=1,j
				w=w+f(i)*r(j-i+1)
			end do
			c=(w-g(j))/v
			
			do i=1,j+1
				f(i)=f(i)-c*a(j-i+1)
			end do	
		end do
	end subroutine stoepd
end module auxsu
