; ****************************************************************************
; 		PV-WAVE CROSS-CORRELATION ROUTINES
;			T. Statler, Dec. 1994
;			Public version Aug. 1996
;	Look for !!!! marks to find things that must be set by the user
;
;	-- Uses PVWave FFT routines. Minimally tested under IDL.
; ****************************************************************************
;
pro corrfit,seed,nshift,galaxy,template,param,chisq,covariance

; --Shift spectrum so that correlation peak will come out around pixel 80.
nshift=fix(nshift)
galaxy=shift(galaxy,nshift)

; !!!!!!!!!
; This is a prescription for estimating the noise in the galactic spectrum.
; It may or may not work for your application
gerr=sqrt(smooth(smooth((galaxy-smooth(galaxy,7))^2,99),99))
; !!!!!!!!!

; !!!!!!!!!
; This is the number of km/s per pixel in the (log-binned) spectra.
; NOTE: This must also be set in another place, below.
delta=46.6892                           ; velocity sampling interval
; !!!!!!!!!

corrpeak,galaxy,gerr,template,xcuse,xceuse,i1,i2
corrauto,template,ac

; !!!!!!!!!
; --Use these three lines for fitting gaussian only.
; scale=[.1,.1,100.]
; start=[(i1+i2)*delta/2.,200.,1.]
; del=[1.,1.,.01]

; --Use these three lines for fitting Gauss-Hermite series up to h4.
 scale=[.1,.1,100.,100.,100.]
 start=[(i1+i2)*delta/2.,200.,1.,0.,0.]
 del=[1.,1.,.01,.01,.01]
; !!!!!!!!!

; --If you want a plot of the fitted part of the cross-correlation peak
; plot,(indgen(i2-i1+1)+i1)*delta,xcuse,psym=4,xra=[i1-.5,i2+.5]*delta, $
;	yra=[min(xcuse),max(xcuse)]
 empty

nxl=1
nx=1
while (nx gt 0 or nxl gt 0) do begin

nxl=nx
corramoeba,start,scale,finish,xcuse,xceuse,ac,i1,i2,nx,chisq,seed
start=finish

endwhile

 makeb,finish,b
 xcm=float(fft(fft(b,-1)*ac,1))
 xcm=xcm(i1:i2)

; --If you want to see the result of the fit overplotted in the XC peak
; oplot,(indgen(i2-i1+1)+i1)*delta,xcm

param=finish

; Calculate covariance matrix

n=n_elements(param)
alpha=fltarr(n,n)
deriv=fltarr(i2-i1+1,n)

for i=0,n-1 do begin
	x1=param
	x1(i)=x1(i)-del(i)/2.
	par=x1
	makeb,par,b
	xcm=float(fft(fft(b,-1)*ac,1))
	xcmu1=xcm(i1:i2)
	x2=param
	x2(i)=x2(i)+del(i)/2.
	par=x2
	makeb,par,b
	xcm=float(fft(fft(b,-1)*ac,1))
	xcmu2=xcm(i1:i2)
	deriv(*,i)=(xcmu2-xcmu1)/del(i)
endfor

for i=0,n-1 do for j=0,n-1 do begin
	di=deriv(*,i)
	dj=deriv(*,j)
	alpha(i,j)= (di # xceuse # dj + dj # xceuse # di)/2.
endfor

covariance=invert(alpha)

; --Unshift spectrum and velocity
param(0)=param(0)-nshift*delta
galaxy=shift(galaxy,-nshift)

end

pro makeb,pp,b
 
; !!!!!!!!!
; --Use this line for fitting gaussian
; p=[pp,0.,0]
 
; --Use this line for fitting up to h4
 p=pp
; !!!!!!!!!

; !!!!!!!!!
; This is the number of km/s per pixel in the (log-binned) spectra.
; NOTE: This must agree with the assignment above
delta=46.6892                           ; velocity sampling interval
; !!!!!!!!!

; !!!!!!!!!
n=1300				; full length of padded spectra in pixels
; !!!!!!!!!

fac=0.3989422
rt3=1.732051
rt24=4.8989795

j=findgen(n)
jj=p(0)/delta
s=abs(p(1)/delta)
gamma=p(2)
 
jnear=fix(jj+.5)
jdel=fix(6*s) > 1
jev1=(jnear-jdel) > 0
jev2=(jnear+jdel) < n-1
jeval=j(jev1:jev2)
argue=(jeval-jj)/s
arg2=argue^2
b0=fac*gamma/s*exp(-arg2/2.)            ; line profile
bigH3=(2.*arg2-3.)*argue/rt3
bigH4=((arg2-3.)*4.*arg2+3.)/rt24
b=fltarr(n)
b(jev1:jev2)=b0*(1.+p(3)*bigH3+p(4)*bigH4)
b = b > 0.

; Make sure broadening function is normalized
; (This is necessary if the b.f. is not well sampled, i.e. for small
; dispersions. However, if your results change when you remove these two
; lines it means that your resolution is too small, your template is lousy,
; or both, and you are courting danger.)
tb=total(b)
if (tb gt 0.) then b=b/total(b)*gamma

end

pro correval,pp,xcuse,xceuse,ac,i1,i2,f

makeb,pp,b

xcmodel=float(fft(fft(b,-1)*ac,1))
xcmodel=xcmodel(i1:i2)
diffs=xcmodel-xcuse
f=diffs # xceuse # diffs /(float(i2-i1+1-n_elements(pp))) ;reduced chi square
f=f(0)

end

pro corramoeba,start,scale,finish,xcuse,xceuse,ac,i1,i2,nx,chisq,seed

n=n_elements(start)
dir=randomu(seed,n)-.5
move=fltarr(n)+dir/abs(dir)
deltareflect=2.
deltacontract=.5
origsize=1.
for i=0,n-1 do origsize=origsize*abs(move(i))
tol=(1.e-6)^n
size=origsize

; Set up initial simplex.

tart=start*scale
points=fltarr(n,n+1)
for i=0,n do points(*,i)=tart
for i=1,n do points(i-1,i)=points(i-1,i)+move(i-1)

; Evaluate the function at the vertices.

fpoints=fltarr(n+1)
for i=0,n do begin
	p=points(*,i)
	 correval,p/scale,xcuse,xceuse,ac,i1,i2,f
	fpoints(i)=f
endfor

niter=0
nx=0
while (niter le 100 and size gt tol) do begin

; Sort function values from lowest to highest; reorder simplex and f values.

so=sort(fpoints)
points=points(*,so)
fpoints=fpoints(so)

; Get center of simplex face opposite high point.

face=points(*,0:n-1)
center=avg(face,1)

; Get reflection vector, reflected point; evaluate f at reflected point.

reflect=center-points(*,n)
new=center+reflect
correval,new/scale,xcuse,xceuse,ac,i1,i2,fnew

; Is the new point at least better than the second-worst previous point?

if (fnew lt fpoints(n-1)) then begin

; If it is, then at least we can replace the old worst point

	points(*,n)=new
	fpoints(n)=fnew
;						print,'Reflected'

; But is it also better than ALL previous points?

	if (fnew lt fpoints(0)) then begin

; If so, check whether further extrapolation will do even better.
; (Remember that the nth point in the list is now the best!)

		new=center+deltareflect*reflect
		 correval,new/scale,xcuse,xceuse,ac,i1,i2,fnew
		if (fnew lt fpoints(n)) then begin
			points(*,n)=new
			fpoints(n)=fnew
			nx=nx+1
;						print,'Extrapolated'
;						print,fnew
		endif
	endif

endif else begin

; OK, so it wasn't better than the second-worst. If it's at least
; better than the worst point, replace the worst before looking for
; a contraction.

	if (fnew lt fpoints(n)) then begin
		points(*,n)=new
		fpoints(n)=fnew
	endif

; Either way, the nth point is still the worst. 
; Now find a point between the current worst point and the face center.

	reflect=center-points(*,n)
	new=center-deltacontract*reflect
	 correval,new/scale,xcuse,xceuse,ac,i1,i2,fnew

; If this is better than the worst, accept it and move on.

	if (fnew lt fpoints(n)) then begin
		points(*,n)=new
		fpoints(n)=fnew
;						print,'Contracted'
	endif else begin

; Otherwise contract toward the low point.

		directions=fltarr(n,n+1)
		for i=1,n do directions(*,i)=points(*,0)-points(*,i)
		points=points+deltacontract*directions

; And we need to evaluate the functions at the moved points

		for i=1,n do begin
			p=points(*,i)
			 correval,p/scale,xcuse,xceuse,ac,i1,i2,f
			fpoints(i)=f
		endfor
;						print,'Frustrated'

	endelse
	
endelse

niter=niter+1
size=1.
for i=0,n-1 do size=size*(max(points(i,*))-min(points(i,*)))

endwhile

; Finish up-- resort; solution is best point.

so=sort(fpoints)
finish=points(*,so(0))/scale
chisq=fpoints(so(0))

; print,'Done!'
 print,finish,fpoints(so(0))

end

pro corrpeak,galaxy,gerr,template,xcuse,xceuse,i1,i2

n=n_elements(template)

; spectra should be pretapered if you want to, before calling this.

; Correlation:

xc=float(fft(fft(galaxy,-1)*conj(fft(template,-1)),1))

; Find peak:

n=n_elements(xc)
sxc=smooth(xc,3)

; !!!!!!!!!
; It is assumed here that the primary correlation peak is going to be found
; within 20 pixels of pixel 80, i.e. that the redshift of the galaxy is
; about 80*delta. You can fiddle with this a bit, but DO NOT change
; it to zero if you have a low-redshift object. See discussion in the
; documentation.
wh=where(sxc(60:100) eq max(sxc(60:100)))
ic=wh(0)+60
; !!!!!!!!!

; Find first minima either side of peak

splus=shift(xc,-1)
sminus=shift(xc,1)
diff=(sminus < xc) - (xc < splus)
wh=where(diff eq 0)
if (wh(0) eq -1) then print,'Problem finding first local minimum'

; !!!!!!!!!
; The 5's below mean that the boundaries of the fitted region will not
; be closer than 5 pixels to the XC maximum. You can change this.
i1=ic-wh
i1(where (i1 le 5))=10000
i1=wh(where(i1 eq min(i1)))
i2=wh-ic
i2(where (i2 le 5))=10000
i2=wh(where(i2 eq min(i2)))
; !!!!!!!!!

i1=fix(i1(0))
i2=fix(i2(0))

; !!!!!!!!!
; Similarly, the 20's below set the boundaries no further than 20 pixels
; either side of peak.
i1=i1>(ic-20)
i2=i2<(ic+20)
; !!!!!!!!!

xcuse=xc(i1:i2)
m=n_elements(xcuse)

; Error bars:

st=dblarr(n,m)
for j=i1,i2 do st(*,j-i1)=double(shift(template,j))*gerr
st=st/float(n)
s=transpose(st)
xceuse=invert(s # st)
xceuse=float(xceuse)

end

pro corrauto,template,ac

; returns the complex conjugate of the fourier transform of the template's
; autocorrelation function.

ac=fft(template,-1)
ac=ac*conj(ac)*n_elements(template)

end

pro filter,spectrum,taperwidth,specfiltered

spe=fft(spectrum,-1)			; transformed spectrum

n=n_elements(spectrum)
taper=fltarr(n)+1.			; build taper for fourier transforms
taperwidth=fix(taperwidth)
theta=findgen(taperwidth+1)/float(taperwidth)*!pi
right=(cos(theta)+1.)/2.
left=1.-right
taper(0:taperwidth)=0.
taper(taperwidth+1:2*taperwidth+1)=left
taper(n-1-taperwidth:n-1)=0.
taper(n-1-2*taperwidth:n-1-taperwidth)=right

spe=spe*taper				; taper the transforms
specfiltered=float(fft(spe,1))

end
;
; $Id: sum.pro,v 1.3 1994/02/24 14:22:24 erickson Exp $
;
FUNCTION SUM,ARRAY,DIMENSION

; NAME:
;	SUM
; PURPOSE:
;	Total up an array over one of its dimensions.
; HISTORY
;	MODIFIED from SUM 91/09/04 by John Black
;	This is a bug fix of the supplied routine SUM, which fails in some
;	circumstances due to 16 bit arithmatic relating to array indices. The
;	main bug I suspect of the origonal routine is the line
;
;	XK = FIX( XIK / NI ) 
;
;       here I've replaced the 'FIX' by a 'LONG'. I've also made alot of the 
;	numbers involved in the for loops long as well for extra safety, but
;       these might not be needed.
;       
;
IF (N_PARAMS() LT 2) THEN BEGIN
	PRINT,'*** Function SUM must be called with two parameters:'
	PRINT,'                   ARRAY , DIMENSION'
	RETURN,ARRAY
ENDIF
;
S = SIZE(ARRAY)
N_DIM = S(0)
IF N_DIM EQ 0 THEN BEGIN
	PRINT,'*** Variable must be an array, name= ARRAY, routine SUM.'
	RETURN,ARRAY
END ELSE IF (DIMENSION GE N_DIM) OR (DIMENSION LT 0) THEN BEGIN
	PRINT,'*** Dimension out of range, name= ARRAY, routine SUM.'
	RETURN,ARRAY
END ELSE IF N_DIM EQ 1 THEN BEGIN	;Trivial case, equivalent to TOTAL.
	F = TOTAL(ARRAY)
	RETURN,F
ENDIF
;
S2 = S(1+WHERE(INDGEN(N_DIM) NE FIX(DIMENSION)));Set up array for output variable.
F = MAKE_ARRAY(DIMENSION=S2, TYPE=S(S(0)+1))
;
;  Calculate product of sizes of dimensions lower than, equal to, and higher
;  than DIMENSION (NI,NJ,NK respectively).
;
NI = 1L ; Make sure that NI is a long integer and result from using it is too
IF DIMENSION GT 0 THEN FOR M = 1,DIMENSION DO NI = NI * S(M)
NJ = S(DIMENSION+1)
NK = 1L ; Make sure that NK is a long integer and result from it is too
IF DIMENSION LT N_DIM-1 THEN FOR M = DIMENSION+2,N_DIM DO NK = NK * S(M)
;
;  Set up index arrays.
;
XIK = LINDGEN( NI * NK )
XJ = LINDGEN( NJ )
NIJ = NI*NJ
;
;  Choose whether it is more efficient to loop over NI and NK ...
;
IF NI*NK LT NJ THEN BEGIN
	FOR I = 0L,NI-1 DO FOR K = 0L,NK-1 DO $
		F(I+NI*K) = TOTAL( ARRAY(I + NI*XJ + NIJ*K) )
;
;  ... or over NJ.
;
END ELSE BEGIN
	XI = XIK MOD NI
	XK = LONG(XIK / NI) ;replaces  XK = FIX( XIK / NI )
	FOR J = 0L,NJ-1 DO F = F + ARRAY(XI + NI*J + NIJ*XK)
ENDELSE
RETURN,F
END
;
; $Id: avg.pro,v 1.1 1991/05/22 16:53:05 jeffry Exp $
;
FUNCTION AVG,ARRAY,DIMENSION

; NAME:
;	AVG
; PURPOSE:
;	Calculate the average value of an array, or calculate the average
;	value over one dimension of an array as a function of all the other
;	dimensions.

ON_ERROR,2
S = SIZE(ARRAY)
IF S(0) EQ 0 THEN BEGIN
	PRINT,'*** Variable must be an array, name= ARRAY, routine AVG.'
	RETURN,ARRAY
ENDIF
IF N_PARAMS() EQ 1 THEN BEGIN
	AVERAGE = TOTAL(ARRAY) / N_ELEMENTS(ARRAY)
END ELSE BEGIN
	IF ((DIMENSION GE 0) AND (DIMENSION LT S(0))) THEN BEGIN
		AVERAGE = SUM(ARRAY,DIMENSION) / S(DIMENSION+1)
	END ELSE BEGIN
		PRINT,'*** Dimension out of range, name= ARRAY, routine AVG.'
		RETURN,ARRAY
	ENDELSE
ENDELSE
;
RETURN,AVERAGE
END