PRO xcorr, x1, y1, x2, y2, ccfx, ccfy

; Compute cross-correlation function for functions y1(x1) and y2(x2)
; The result is returned as ccfy(ccfx)

; NB.  Currently assumes that x1 and x2 are identical.

nwave=n_elements(y1)
nccf  = 131072L
IF nwave LT 65537L THEN nccf = 65536L
IF nwave LT 32769L THEN nccf = 32768L
IF nwave LT 16385L THEN nccf = 16384L
IF nwave LT  8193L THEN nccf =  8192L
IF nwave LT  4097L THEN nccf =  4096L
IF nwave LT  2049L THEN nccf =  2048L
IF nwave LT  1025L THEN nccf =  1024L
ccfy  = fltarr(nccf)

;    construct fast fourier transforms
testspec = fltarr(nccf)
testspec(0:nwave-1) = y1
ff1 = fft(testspec)
template = fltarr(nccf)
template(0:nwave-1) = y2
ff2 = fft(template)

;    convolve transforms
     ffxff = ff2 * conj(ff1) * nccf/8

;    ccf is real part of inverse transform
     ccfy = float (fft ( ffxff, +1 ))

;    swap order around
     work = ccfy (nccf/2:nccf-1)
     ccfy(nccf/2:nccf-1) = ccfy(0:nccf/2-1)
     ccfy[0]=work

;    compute ordinate
     dw = (x1(nwave-1) - x1(0))  / (nwave-1)
     ccfx = (findgen(nccf) - nccf/2) * dw

;    scaling factors to produce identity with DIPSO
     ccfx = ccfx 
     ccfy = ccfy / 2

END