PRO xcorrelate, star, temp_file, temp_v, temp_b, starname, wavwind

;; cross-correlates individual spectra with means

;; star ; select type of template (self, model, self-model)
;; temp_file ; file containing model spectrum
;; temp_v ; shift model by this velocity
;; temp_b ; broaden model by this width 



COMMON fitsblock
COMMON procblock
COMMON xcorblock, ccfx, ccfy, ccf2y, ccf2n, nccf, afits, template


;; Announce programme and options

print, 'XCORRELATE, star, gplot'
print
print, 'Cross-correlate a series of spectra with a template'
print 
IF (star eq 0) THEN $
    print, 'star  =', star, '; template = mean spectrum'
IF (star eq 1) THEN $
     print, 'star  =', star, '; template = model'
IF (star eq 2) THEN $
    print, 'star  =', star, '; template = model - mean'


;; Define holders for CCFs

nccf  = 4096
; IF nwave < 1025 THEN nccf = 1024
ccf2n = fltarr(nccf, nspec)
ccf2y = fltarr(nccf, nspec)
ccfy  = fltarr(nccf)
afits = fltarr(6,nspec)

; x-axis for ccf
;     mode;  all;  Hg;  Hd;  Fe;  all
ix1 =  0;   53; 800; 150; 330;   53
ix2 =  nwave-1; 1003; 920; 320; 800; 1003

IF (n_elements(wavwind) EQ 2 ) THEN BEGIN
ax1 = where(wave GE wavwind(0))
ax2 = where(wave LE wavwind(1))
ix1 = ax1(0)
ix2 = ax2(n_elements(ax2)-1)
ENDIF



; wavelength
IF ( wave(0) gt 10. ) THEN BEGIN
  dw = (wave(ix2) - wave(ix1))  / (ix2-ix1-1)
  mw =  wave(nwave/2)
  ccfx = (findgen(nccf) - nccf/2) * dw / mw * 2.99792458e5 

; log wavelength
ENDIF ELSE BEGIN
  dw = (wave(ix2) - wave(ix1))  / (ix2-ix1-1)
  ccfx = (findgen(nccf) - nccf/2) * dw 
  ccfx = (10.^(-ccfx) - 1.) * 2.99792458e5 

ENDELSE



;;  TEMPLATE definition

;  Extract the template from the mean spectrum (specm)
xcorr_template, specm, nccf, ix1, ix2, template, wave, starname

colors = GetColor(/Load)
IF (strupcase(!d.name) EQ 'X') THEN BEGIN
     window,3,title="xcorrelate"
    !p.multi=[0,0,2,0,0]
     plot, wave,template; ,yrange=[-0.6,0.2]
ENDIF

; do we want to redefine the template ?
IF (star gt 0) THEN BEGIN

; read model as template
  bluestar = fltarr(nccf)
  alasrd,temp_file,tw,tf
  vcorr,tw,tf,temp_v/2.9979E5,1
  dw = (wave(nwave-1)-wave(0))/(nwave-1)
  tf = qqsm ( tw,tf,temp_b/2.534 )
  resample,tw,tf,wave(0),wave(nwave-1),dw
  bluestar(0:nwave-1) = tf - 1.0

; define red star spectrum as difference
  redstar  = template - bluestar


; redefine template
  IF (star eq 1) THEN template = bluestar
  IF (star eq 2) THEN template = redstar
  IF (star eq 3) THEN template = helium

ENDIF

IF (strupcase(!d.name) EQ 'X') THEN BEGIN
IF (star gt 0) THEN BEGIN
    oplot,wave,bluestar-0.1,color=colors.blue
    oplot,wave,redstar+0.1,color=colors.red
ENDIF
oplot,wave,template,color=colors.green

    plot, wave,template; ,yrange=[-0.6,0.2]
ENDIF



;;  AAUTOCORRELATION FUNCTION

; Construct acf for subtracting from ccf
fftemp   = fft(template)
ffxff    = fftemp * conj(fftemp) * (nccf/8)
acfy     = float (fft ( ffxff, +1 ))

; Swap order around
work = acfy (nccf/2:nccf-1)
acfy(nccf/2:nccf-1) = acfy(0:nccf/2-1)
acfy(0:nccf/2-1)=work



testspec = fltarr(nccf)
IF (ix1 gt 1) THEN testspec (0:ix1-1) = 0
IF (ix2 lt nccf-2) THEN testspec (ix2+1:nccf-1) = 0

; Loop over spectra
FOR it = 0,nspec-1 DO BEGIN

;  construct test spectrum
   testspec(ix1:ix2) = specf(ix1:ix2, it) - 1.0

IF (strupcase(!d.name) EQ 'X') THEN  plot, wave,testspec,yrange=[-0.6,0.2]

;  check data is valid 
   ccfy = 0
   IF ( stddev(testspec(ix1:ix2)) gt 1.e-6 ) THEN BEGIN

;    construct fast fourier transform
     fftest = fft(testspec)

;    convolve transforms
;     ffxff = fftemp * conj(fftest) * nccf/8
     ffxff = fftest * conj(fftemp) * 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:nccf/2-1)=work

;    save residual ccf in 2d array 
     ccf2n(*,it) = ccfy - acfy 

   ENDIF

;  save full ccf in 2d array
   ccf2y(*,it) = ccfy

IF (strupcase(!d.name) EQ 'X') THEN  plot,ccfy

ENDFOR


!p.multi=0

END