PRO analyse_msst, scope

;; analyse_msst
;;
;; Purpose:
;;
;;  Analysis of MSST 4m time-series spectroscopy 
;;
;; Use:
;;
;;  analyse_msst, <scope>
;;
;; Example:
;;
;;  analyse_msst, 'apo_red'
;;  analyse_msst, 'twin_blue'
;;  analyse_msst, 'twin_red'
;;
;;

;; Standard FAST common blocks
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


;; identify the dataset with an index
IF scope EQ 'apo_red'   THEN isc = 0
IF scope EQ 'twin_blue' THEN isc = 1
IF scope EQ 'twin_red'  THEN isc = 2

  restore, filename=scope+'_norm'
  specf = spec

;; compute the cross-correlations

  dm = alog10 ( lc ) 

  afits = fltarr(6,nspec)

  xcorrelate,0,1, 'MSST: '+scope

  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,1.0]
    FOR i=0,nspec-1 DO oplot, ccfx, ccf2y(*,i)
    window,2
    shade_surf, ccf2y(wfl:wfu,*), ccfx(wfl:wfu), time, zrange=[0,1.0]
  ENDIF


;; measure radial velocities from the ccfs

  xcorr_fits,5,1,0,0


;; rescale velocities from log space
;; IF wave(0) < 10. THEN afits(1,*) = (10^(-afits(1,*) - 1) * 2.99792458e5)


;; analyse the velocity shifts

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

  dft_periods, tv,rv, clf,cldft, 0
print,'V338 Ser: velocities'
  find_peaks, clf, cldft, 7, [2.08,1.99,1.89,2.10,2.74,3.99,4.06], tp, pk


;; analyse the equivalent width variations

specr=spec
FOR is = 0,nspec-1 DO specr(*,is) = spec(*,is) - specm 


  ew_extract, wave, time, specr, ew

  dft_periods, time,ew, ewf, ewdft, 0
print,'V338 Ser: equivalent widths'
  find_peaks, ewf, ewdft, 2, [2.08,2.10], tp, pk


;; analyse the light curve variations

  dft_periods, time,dm, frq,dft, 0
print,'V338 Ser: light curve'
  find_peaks, frq, dft, 4, [2.08,1.99,1.89,2.10], 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=[1,5]
    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,2
    !p.multi=[0,2,2,0,0]
    nf = n_elements(cldft)
    plot,tv  ,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