PRO analyse_bruce, label, colour, mode, xfiles, xruns

;; analyse_bruce
;;
;; Purpose:
;;
;;  Analysis of simulated time-series spectroscopy of pulsating star
;;  spectra. 
;;
;; Use:
;;
;;  analyse_bruce, 'star', 'mode', 1, nruns
;;  analyse_bruce, 'star', 'mode', 0
;;
;; Example:
;;
;;  analyse_bruce, 'kpd2109', 'A', 1,  [1500]
;;  analyse_bruce, 'kpd2109', 'A', 0
;;
;;


COMMON fileblock,  nfiles, folder, flabel, nruns, nbomit, nwomit
COMMON fitsblock,  time, wave, spec, nwave, nspec
COMMON procblock,  specf, specr, specm
COMMON pspecblock, freq, twod_dft
COMMON xcorblock,  ccfx, ccfy, ccf2y, ccf2n, nccf, afits, template


;; read in a previous dataset 

IF xfiles EQ 0 THEN BEGIN


  restore, filename='idl_'+label+'.sav'



;; read in a new dataset

ENDIF ELSE BEGIN

;; define new dataset parameters

  nfiles = xfiles
  folder = strarr(nfiles)
  folder = ['./'+label]
  flabel = strarr(nfiles)
  flabel = [mode]
  nruns  = intarr(nfiles)
  nruns  = xruns
  nspec  = total(nruns)
  nwave  = 2751
  nwomit = 0                ;; omit rubbish at start of spectrum
  nbomit = 0                ;; miss first nbomit files from each run....

  ndft   = 5545             ;; dimension of dft


;; setup dataset array dimensions

  time   = fltarr(nspec)
  wave   = fltarr(nwave)
  spec   = fltarr(nwave,nspec)
  spec   = spec + 1         ;; dummy values are continuum

  specm  = fltarr(nwave)
  specf  = fltarr(nwave,nspec)
  specr  = fltarr(nwave,nspec)

  freq   = dblarr(2*ndft)
  twod_dft = fltarr(nwave,ndft)


;; read in and clean the selected data

  brucedata

  IF (strupcase(!d.name) EQ 'X') THEN BEGIN
    window,1
    shade_surf,spec,wave,time
  ENDIF


;; compute and analyse the light variations
  
  lightcurve, wave, time, spec, lc, dm
   
  IF (strupcase(!d.name) EQ 'X') THEN BEGIN
    window,0
    plot,time,dm,xtitle='day',ytitle='dmag'
  ENDIF


;; normalise the data

IF colour EQ 'blue' THEN BEGIN
  normalise, [3855, 3860, 3930, 3945, 4010, 4013, 4219, 4228, 4399, 4400, 4555, 4570, 4720, 4740, 4945, 4950]
;  normalise, [3855, 3860, 3930, 3945, 4010, 4013, 4219, 4228, 4399, 4400, 4555, 4570, 4720, 4740, 4945, 4950]
ENDIF

IF colour EQ 'red' THEN BEGIN
  normalise, [6500, 6502, 6615, 6625, 6690, 6700]
;  normalise, [6500, 6510, 6615, 6625, 6690, 6700]
ENDIF



  IF (strupcase(!d.name) EQ 'X') THEN BEGIN
    window,1
    shade_surf,spec,wave,time
  ENDIF


;; clean the data

  cleandata

  IF (strupcase(!d.name) EQ 'X') THEN BEGIN
    window,0
    plot,wave,specm,xtitle='Angstrom',ytitle='F/Fc'
    window,1
    shade_surf,specr,wave,time
  ENDIF



;; save the data for reanalysis

  save, filename='idl_'+label+'.sav', /all

ENDELSE


;; compute the cross-correlations

  dm = alog10 ( lc ) 

  afits = fltarr(6,nspec)

  xcorrelate,0,1, 'BRUCE '+'mode'

  IF (strupcase(!d.name) EQ 'X') THEN BEGIN
    wfl=nccf/2-nccf/6
    wfu=nccf/2+nccf/6
    window,0
    plot, ccfx, ccf2y(*,10), xr = [-1000,1000], yr = [0,3.0]
    FOR i=0,nspec-1 DO oplot, ccfx, ccf2y(*,i)
    window,1
    shade_surf, ccf2y(wfl:wfu,*), ccfx(wfl:wfu), time
  ENDIF


;; measure radial velocities from the ccfs

  xcorr_fits,5,1,0,0


;; analyse the velocity shifts

  xvclean,time,afits(1,*),1,0, tclean, rv

  dft_periods, time,rv, clf,cldft, 0
  find_peaks, clf, cldft, 1, [2.0822], tp, pk

;; analyse the equivalent width variations

  ew_extract, wave, time, specr, ew

  dft_periods, time,ew, ewf, ewdft, 0
  find_peaks, ewf, ewdft, 1, [2.0822], tp, pk


;; analyse the light curve variations

  dft_periods, time,dm, frq,dft, 0
  find_peaks, frq, dft, 1, [2.0822], tp, pk



  set_viewport

;; plot all the periodograms together
  IF (strupcase(!d.name) EQ 'X') THEN BEGIN
    window,0
    !p.multi=[0,2,2,0,0]
    !x.range=[0,10.]
    nf = n_elements(cldft)
    plot,clf(1:nf),cldft,xtitle='mHz',ytitle='km/s'
    plot,frq(1:nf),dft,xtitle='mHz',ytitle='dmag'
    plot,ewf(1:nf),ewdft*1000,xtitle='mHz',ytitle='mA'
    !p.multi=0
    !x.range=0
  ENDIF


;; plot all the variables together
  IF (strupcase(!d.name) EQ 'X') THEN BEGIN
    window,1
    !p.multi=[0,2,2,0,0]
    nf = n_elements(cldft)
    plot,time,rv,xtitle='day',ytitle='km/s'
    plot,time,dm,xtitle='day',ytitle='dmag'
    plot,time,ew,xtitle='day',ytitle='mA'
    !p.multi=0
    !x.range=0
  ENDIF


;; save selected results
;
;  out2dip,'BRUCE',  wave, template,  clf,cldft,  time, afits(1,*), bg_t, bg_v
;


END