asym_ind.pro

                       pro asymmetry_analysis
                        
                       common path, path, n_file 

;
; This IDL procedure computes the CCF asymmetry indexes, its FWHM, and its contrast using files provided by the HARPS Data Reduction Software, 
; see Lanza et al. 2018, A&A, in press. It is recommended that you read that paper and in particular its Appendix A before compiling and using 
; this IDL procedure. Please cite that paper if you use this procedure to produce data appearing in a scientific publication or a thesis. 
; 
; To compile this procedure, first you have to compile the procedure file MPFIT.PRO and then you have to compile this file twice 
; because its first compilation will report an apparent error. MPFIT.PRO is accessible at http://cow.physics.wisc.edu/~craigm/idl/fitting.html
; 
; Input information to be specified 
; 
; path: specify the path of the subdirectory containing the *_ccf_*.fits and *_bis_*.fits files on your own computer; 
; these file are produced by the HARPS Data Reduction Software (DRS) that is described in the DRS User Manual and HARPS-N User Manual 
; (see Sect. 8 - Data Products) on the web site: http://www.tng.iac.es/instruments/harps/
; 
; file_list_ccf [string]: a file in the path subdirectory listing all the *_ccf_* files to be analysed
;
; file_list_bis [string]: a file in the path subdirectory listing all the *_bis_* files to be analysed
;
; operation [string]: a string specifying interactive operation with plots on the X11 terminal and warnings for the user; 
;            set it to 'interactive' for this option, otherwise leave it blank to have an automated processing [not recommended for the first run];
;            even when the processing is automated, the procedure may still ask the user for confirmation when some deviation from the expected data accuracy 
;            is detected. Therefore, pay attention on a possible request from the procedure even if you put operation=' '; 
;            
; spectrogr_string [string]: specifies which spectrograph has acquired the data; it is required to read the keywords fields 
;            in the fits files; it can be 'TNG' for HARPS-N@TNG or 'ESO' for HARPS@ESO La Silla; please check the dimensions of the matrixes extracted
;            from the FITS files of the DRS of HARPS@ESO because it can be different from those of the files of the DRS of HARPS-N@TNG used in this version.  
;
; file_out_ind [string]: the name of the output file where the output will be written. It is an ASCII file that contains the following columns from left to right 
;                        [they are all floating point numbers, except for maxcpp that is an integer]: 
;                        
;  - time: the barycentric julian date as given by the DRS keyword DRS BJD;
;  - rv_ccf: the radial velocity as given by the DRS keyword CCF RVC in km/s;
;  - rv_ccf_err: the error on the radial velocity as given by the DRS keyword CCF NOISE in km/s;
;  - maxcpp: the maximum number of counts per pixel as given by the DRS keyword CCF MAXCPP;
;  - bisector: the BIS calculated as specified in Lanza et al. 2018 in km/s;
;  - sterr_bis: the standard error on the BIS calculated as specified in Lanza et al. 2018 in km/s;
;  - delta_V_nardetto: the \Delta V indicator calculated as specified in Lanza et al. 2018 in km/s;
;  - sigma_delta_V: the standard deviation of the \Delta V indicator in km/s calculated by mpfit.pro by means of the covariance matrix and the best fit residuals [not recommended for use, see Lanza et al. 2017];
;  - v_asy: the old V_asy indicator as defined by Figueira et al. 2013, A&A 557, A93; it is given only for 
;           reference, but it is not recommended for use because of its systematic correlation with the RV;
;  - stdev_v_asy: the standard deviation of the old V_asy indicator calculated as specified in Lanza et al. 2018;
;  - v_asy_mod: the new CCF asymmetry indicator V_asy(mod) defined and calculated as specified in Lanza et al. 2018;
;  - stdev_v_asy_mod: the standard deviation of the new CCF asymmetry indicator v_asy_mod calculated as specified in Lanza et al. 2018;
;  - fwhm_drs: the full width at half maximum of the CCF as given by the DRS keyword CCF FWHM in km/s;
;  - fwhm_gaussfit: the full width at half maximum (FWHM) of the CCF as computed by our Gaussian best fit with mpfit.pro in km/s;
;  - sigma_fwhm_gaussfit: the standard deviation of the FWHM in km/s as computed by mpfit.pro by means of the covariance matrix and the residuals of the best fit [not recommended for use, see Lanza et al. 2017];
;  - contrast_drs: the fractional depth of the CCF with respect to its continuum (CCF contrast) as given by the DRS keyword CCF CONTRAST;
;  - contrast_gaussfit: the CCF contrast as given by our Gaussian best fit with mpfit.pro;
;  - sigma_contrast_gaussfit: the standard deviation of the CCF contrast as given by our Gaussian best fit using the covariance matrix and the residuals of the best fit calculated by mpfit.pro [not recommended for use, see Lanza et al. 2017];
;  - fwhm_set: the FWHM in km/s obtained as the median of the 200 realizations with random and correlated noise [recommended for use, see Lanza et al. 2018];
;  - sigma_fwhm_set: the standard deviation of the FWHM in km/s as obtained from the 200 realizations with random and correlated noise [recommended for use];
;  - delta_v_nardetto_set: the \Delta V indicator in km/s obtained as the median of the 200 realizations with random and correlated noise [recommended for use];
;  - sigma_delta_v_set: the standard  deviation of the \Delta V indicator in km/s obtained from the 200 realizations with random and correlated noise [recommended for use];
;  - contrast_set: the fractional depth of the CCF with respect to its continuum (CCF contrast) obtained as the median of the 200 realizations with random and correlated noise [recommended for use];
;  - sigma_contrast_set: the standard deviation of the CCF contrast as obtained from the 200 realizations with random and correlated noise [recommended for use]. 
;
;
;   VERSION 1.0 - 10 April 2018 - Author: A. F. Lanza with contributions by L. Malavolta, S. Desidera, A. Bignamini 
;   e-mail: nuccio.lanza@oact.inaf.it 
;
;
; the path to the subdirectory containing the *_ccf_* and *_bis_* DRS files

   path='/Users/nlanza/home/GAPS/WP5000/Asymmetry_indicators/paper/IDL_macro_for_distribution/distribution_cds/input/'
   
; specify the file with the list of the names of *_ccf_* files 
; note an example file name for a *_ccf_*.fits file: 'HARPN.2012-09-29T05-46-30.718_ccf_K5_A.fits' ; 
; files not to be read, e.g., because of low S/N, can be flagged with a # that comments that line.

 file_list_ccf='flistccf.lis' 

; specify the file with the list of the names of the *_bis_* files
 file_list_bis='flistbis.lis'
 
; specify the operation mode: 'iteractive' or blank ' ' 
 operation='interactive' ; or ' ' for automated operation without plots and user checks    
 
; specify which spectrograph have acquired the data; for HARPS@ESO use 'ESO'; for HARPS-N@TNG use 'TNG'
 spectrogr_string='TNG' 

; reading the names of the *_ccf_* and *_bis_* files 
 
 data_file=path+file_list_ccf
 readcol, data_file, filenames, format='a', comment='#'
 n_spectra=n_elements(filenames)
 
 path_bisect_files=path
 list_bis_files=path_bisect_files+file_list_bis
 readcol, list_bis_files, filenameb0, format='a', comment='#'
 n_spectra_bis=n_elements(filenameb0)
 
 if(n_spectra_bis ne n_spectra) then begin 
  print, 'The number of CCF files is not the same as the number of BIS files: stop'
  stop
 endif 

 filenameb=path_bisect_files+filenameb0

 
 close, 1,2,3,4 ; closing possibly open units 

; specify the name of the output file where the indicators are listed (it will be written in the same directory as the present procedure file) 
 
 file_out_ind='asymmetry_analysis_HD108874.dat'
 openw, 1, file_out_ind 
  
; analysing spectra by extracting information from the *_ccf_* and *_bis_* files one after the other
     
     for j=0, n_spectra-1 do begin 

      n_file=j      
                             
      filename_ccf=filenames[j]
      filename_bis=filenameb[j] 

          asymmetry_indexes, filename_ccf, filename_bis, operation, spectrogr_string, time, rv_ccf, rv_ccf_err, maxcpp, bisector, sterr_bis, delta_V_nardetto, sigma_deltaV, v_asy, stdev_v_asy, $
                v_asy_mod, stdev_v_asy_mod, fwhm_drs, fwhm_gaussfit, sigma_fwhm_gaussfit, contrast_drs, contrast_gaussfit, sigma_contrast_gaussfit, $ 
                fwhm_set, sigma_fwhm_set, delta_v_nardetto_set, sigma_delta_v_set, contrast_set, sigma_contrast_set



                  
     printf, 1, time, rv_ccf, rv_ccf_err, maxcpp, bisector, sterr_bis, delta_v_nardetto, sigma_deltaV, v_asy, stdev_v_asy, $
               v_asy_mod, stdev_v_asy_mod, fwhm_drs, fwhm_gaussfit, sigma_fwhm_gaussfit, contrast_drs, contrast_gaussfit, sigma_contrast_gaussfit, $
               fwhm_set, sigma_fwhm_set, delta_v_nardetto_set, sigma_delta_v_set, contrast_set, sigma_contrast_set, format='(1x,3f20.8,i9,20f20.8)'
     
  endfor 
  close, 1    
     
                             stop 
                             end 
;*************************************************************************

                    
                       
pro asymmetry_indexes, filename_i, filename_bi, oper, sp_string, time, rv_ccf, rv_ccf_err, maxcpp, bisector, sterr_bis, delta_V_nardetto, sigma_deltaV, v_asy, stdev_v_asy, $
  v_asy_mod, stdev_v_asy_mod, fwhm_drs, fwhm_gaussfit, sigma_fwhm_gaussfit, contrast_drs, contrast_gaussfit, sigma_contrast_gaussfit, $ 
   fwhm_set, sigma_fwhm_set, delta_v_nardetto_set, sigma_delta_v_set, contrast_set, sigma_contrast_set


                             
                           common path, path, number_file    
                            
                             set_plot, 'x' 
                             !p.thick=1.0
                             !p.charsize=1.5
                             !p.charthick=1.0
                             
                             if(oper eq 'interactive' or oper eq 'INTERACTIVE') then begin 
                                operflag=1
                             endif else begin
                                operflag=0 
                             endelse
                             
                             filename=path+filename_i
                             result0=readfits(filename, header)
  
                                           
                             crval1=sxpar(header, 'CRVAL1')
                             cdelt1=sxpar(header, 'CDELT1')
                             nx=sxpar(header, 'NAXIS1')
                             ny=sxpar(header, 'NAXIS2')
                             
                             if(ny ne 70) then begin 
                               print, 'The second dimension of the CCF matrix is not 70: STOP'
                               stop
                             endif 
                             
                             snr=dblarr(ny)
                             
; reading the S/N ratios of the different orders 
                        
                               keyword='HIERARCH '+sp_string+' DRS SPE EXT SN'
                             
                             for i=0, ny-2 do begin ; note ny-2 instead of ny-1 because there is no SNR value for ny-1
                               if(i le 9) then begin
                                 fmt='(i1)'
                                 keys=keyword+string(i,format=fmt)
                                 x=where(strmid(header,0,28) eq keys)
                                 snr[i]=float(strmid(header[x[0]], 31, 4))
                               endif else begin 
                                 fmt='(i2)'
                                 keys=keyword+string(i,format=fmt)
                                 x=where(strmid(header,0,29) eq keys)
                                 snr[i]=float(strmid(header[x[0]], 32, 4))
                               endelse                              
                             endfor 


; reading MAXCPP parameter 
                           
 
                             keyword='HIERARCH '+sp_string+' DRS CCF MAXCPP'
                             x=where(strmid(header,0,27) eq keyword)
            
                             maxcpp=double(strmid(header[x[0]], 30, 6))
                             print, 'Maximum number of count per pixel in the CCF (e) ', maxcpp

; expressing the CCF in photoelectron counts times the line mask weights per pixel extracting it from the last line of the matrix from the CCF file

                             result=result0[*,69]
                                    
                             ccf=result ; total(result, 2) 
                             stdev_ccf=sqrt(ccf) ; assuming Poisson photon shot noise to maximize the error (see Lanza et al. 2018) 
                              
; defining the radial velocity axis values
                              
                             xaxis=cdelt1*dindgen(nx)+crval1
                             
; estimating the continuum by a Gaussian fit to the CCF                         
                         
                         nterms=4
                         result_gaussfit=gaussfit(xaxis, ccf, a_gauss_fit, nterms=nterms)
                         continuum_level_gfit=a_gauss_fit[3]
  
                             
; plotting the CCF with errorbars (one standard deviation) 
                  if(operflag eq 1) then begin 
                             plot, xaxis, ccf, linestyle=0, /ynozero, title=' Cross-Correlation Function (CCF)', $ 
                               xtitle='RV (km/s)', ytitle='CCF counts ', pos=[0.15, 0.15, 0.9, 0.9] 
                             errplot, xaxis, ccf-stdev_ccf, ccf+stdev_ccf 
                             oplot, [min(xaxis), max(xaxis)], [continuum_level_gfit, continuum_level_gfit], linestyle=1
                             
                             print, 'Plotting the CCF for data file no. ', number_file
                             print, 'The dotted line is the continuum level estimated from the Gaussian best fit'
                             print, 'Do you want to procede ? (press any key to continue) '
                             ans=get_kbrd()  
                  
                  endif        
                             data_file=sxpar(header, 'FILENAME')
                             print, 'Name of the data file ', data_file
                             
                             
                             time_str_des='HIERARCH '+sp_string+' DRS BJD' 
                             x=where(strmid(header,0,20) eq time_str_des)
                             time=double(strmid(header[x[0]], 23, 16))
                             
                             print, 'Modified BJD of the observation (i.e., BJD-2400000) ', time-2400000.0d0 

; reading the RV from the data file 

                             
                             keyword_rv='HIERARCH '+sp_string+' DRS CCF RVC'
                        
                             
                             x=where(strmid(header,0,24) eq keyword_rv)
                             rv_ccf=double(strmid(header[x[0]], 27, 16))
                             
                             print, 'RV from the DRS ', rv_ccf, ' (km/s)' 

; reading the error of the RV from the data file

                       
                       keyword_rv_err='HIERARCH '+sp_string+' DRS CCF NOISE'
                       x=where(strmid(header,0,26) eq keyword_rv_err)
                       rv_ccf_err=double(strmid(header[x[0]], 29, 19))
                             
                       print, 'RV from the DRS ', rv_ccf, ' (km/s)' 
                       print, 'Error on RV (photon shot noise only) ',rv_ccf_err, ' (km/s)'

;                       
; reading the FWHM of the CCF from the data file (in km/s) 
;                     
                      
                      keyword_fwhm_drs='HIERARCH '+sp_string+' DRS CCF FWHM'
                      x_fwhm=where(strmid(header,0,25) eq keyword_fwhm_drs)
                      fwhm_drs=double(strmid(header[x_fwhm[0]], 28, 17))
                      
; reading the contrast of the CCF from the data file (in percent, then converted to fraction of the unity)
                      
                      keyword_contrast_drs='HIERARCH '+sp_string+' DRS CCF CONTRAST'
                      x_contrast=where(strmid(header,0,29) eq keyword_contrast_drs)
                      contrast_drs=double(strmid(header[x_contrast[0]], 32, 17))/100.0d0
                                      
; defining the normalized CCF  
                         
                         ccfn=ccf/continuum_level_gfit
                         min_ccfn=min(ccfn, imin)
                         
; computing the BIS from the normalized CCF                                                  

                         ndepth=100 ; number of intervals considered along the full normalized CCF depth
                         delta_int=(1.0d0-min_ccfn)/double(ndepth)
                         
                         wave_bis=dblarr(ndepth)
                         bis_xaxis0=dblarr(ndepth)
                         
                         for k=1, ndepth-1 do begin
                            bis_xaxis0[k]=(1.0d0/double(ndepth))*double(k)
                            flux_level=1.0d0-delta_int*double(k)
                            wave1=interpol(xaxis[0:imin], ccfn[0:imin], flux_level, /QUADRATIC)
                            wave2=interpol(xaxis[imin:(nx-1)], ccfn[imin:(nx-1)], flux_level, /QUADRATIC)
                            wave_bis[k]=0.5d0*(wave1+wave2)                         
                         endfor 
                         wave_bis[0]=wave_bis[1]
                         bis_xaxis0[0]=bis_xaxis0[1]
                         
                         bisector= mean(wave_bis[10:40])-mean(wave_bis[60:90])
                         print, 'BIS = ', bisector, ' (km/s)'

; a posteriori estimate of the standard error of the BIS 
 
                         sterr_bis=stdev(wave_bis[10:40]-wave_bis[60:90])/sqrt(10.0d0) 
                         
                         print, 'BIS standard deviation  = ', sterr_bis, ' (km/s)'
                         print, 'BIS error from RV error = ', (2.5d0)*rv_ccf_err, ' (km/s)'

; if the BIS error estimated from the RV error is larger, then we use it                         

                         if (((2.5d0)*rv_ccf_err) ge sterr_bis) then sterr_bis=(2.5d0)*rv_ccf_err


; reading the BIS file; the DRS gives it from 5 to 95 percent relative depth in steps of 0.01 of relative depth 
; (see the header of the bisector file) 

                        
                             filename0=filename_i
      
                             
                             comp_bis=readfits(filename_bi, header_bis)
                             
                             cdelt1_bis=sxpar(header_bis, 'CDELT1')
                             crval1_bis=sxpar(header_bis, 'CRVAL1')
                             nbis=sxpar(header_bis, 'NAXIS1') 
                             
                             bis_xaxis=crval1_bis+cdelt1_bis*dindgen(nbis)

                             
                             
                             vindex=indgen(100)
                             vindex1=[vindex[10:40],vindex[60:90]]
                             vindex2=[vindex[5:35],vindex[55:85]]
                             
             if(operflag eq 1) then begin                       
                             plot, bis_xaxis0[vindex1], wave_bis[vindex1], psym=4, /ynozero, xtitle='Normalized depth', $ 
                                   ytitle='RV (km/s)', title='CCF line bisector'  
                             oplot, bis_xaxis[vindex2], comp_bis[vindex2], linestyle=0 
                             
                             print, 'Plotting CCF bisectors of file no. = ', number_file 
                             print, 'Diamonds - presently computed bisector; solid line - bisector from DRS _bis_ file'
                             print, 'Should we continue ? (press any key to continue) '
                             anss=get_kbrd()
             endif 

; checking the closeness of the BIS computed here to that given by the HARPS DRS 

                             
                             delta_bis0=comp_bis[5:35]-comp_bis[55:85]
                             delta_bis1=wave_bis[10:40]-wave_bis[60:90]
                             delta_bis_m=abs(mean(delta_bis0)-mean(delta_bis1))/max([mean(delta_bis0),mean(delta_bis1)])
                             
                             if(delta_bis_m gt 0.01d0) then begin 
                               print, 'For ', filename0, ' the computed BIS differs from that of the HARPS-N DRS by more than 1 percent '
                               print, delta_bis_m
                               print, 'MAXCPP ', maxcpp
                               print, 'Do you want to procede ? (Y/N) '
                               ans=get_kbrd()
                               if(ans eq 'N' or ans eq 'n') then stop
                             endif 
                              
;
; computing the bi-Gaussian best fit to the normalized CCF below some threshold level;  
; the IDL function MPFIT.PRO, that applies Levenberg-Marquardt minimization, is used. It should be
; compiled in advance from a separate file: mpfit.pro [see the header of this procedure] 

;
; now preparing to fit the CCF with a Gaussian or a bi-Gaussian function by means of a
; Levenberg-Marquardt constrained minimization algorithm as implemented in the IDL procedure
; MPFIT.PRO 
;

; only the part of the normalized CCF below THRESHOLD is considered for these best fits 

                             threshold=0.95d0*continuum_level_gfit 
                             
                             
                             flim=where(ccfn le threshold) 
                             
                             if(flim[0] eq -1) then begin
                               print, 'Not enough points along the CCF below the threshold: STOP'
                               print, filename
                               stop 
                             endif 
                             
                             ccfnf=ccfn[flim]
                             xaxisf=xaxis[flim]
                             measure_errors=stdev_ccf[flim]/max(ccf)
                            
                             nxf=n_elements(flim)
             

                             
             fa={xaxisf:xaxisf, ccfnf:ccfnf, measure_errors:measure_errors}

             parinfo=replicate({value:0.D, fixed:0, limited:[0,0], limits:[0.D,0.D]},4)
             
; estimating the parameters and fixing their allowed ranges of variations for a constrained minimization of the chi square

                           pcon_e=1.0d0-min(ccfnf,iminf) ; central depth of the Gaussian in relative units
                              
                              if(pcon_e lt 0.0d0) then begin 
                                print, 'Error in the depth of the CCF: STOP'
                                stop
                              endif 
                              
                            rvc_e=xaxisf[iminf] ; central RV of the CCF 
                            
                            flux_level=1.0d0-(delta_int*0.5d0*ndepth)
                            wave1=interpol(xaxisf[0:iminf], ccfnf[0:iminf], flux_level, /QUADRATIC)
                            wave2=interpol(xaxisf[iminf:(nxf-1)], ccfnf[iminf:(nxf-1)], flux_level, /QUADRATIC)
                              
                              fwhm_e=double(abs(wave2-wave1)) ; full width at half maximum of the CCF
                              a_e=0.01d0 ; asymmetry parameter in the bi-Gaussian function
           
               parinfo(0).fixed=0
               parinfo(0).limited(0)=1
               parinfo(0).limited(1)=1
               parinfo(0).limits(0)=0.25d0*pcon_e
               parinfo(0).limits(1)=2.0d0*pcon_e

               parinfo(1).fixed=0
               parinfo(1).limited(0)=1
               parinfo(1).limited(1)=1
               
               if(rvc_e gt 0.0d0) then begin 
                    rv1=0.25d0*rvc_e
                    rv2=2.0d0*rvc_e
               endif else begin 
                    rv1=2.0d0*rvc_e
                    rv2=0.25d0*rvc_e
               endelse 
               
               parinfo(1).limits(0)=rv1
               parinfo(1).limits(1)=rv2
               
               parinfo(2).fixed=0
               parinfo(2).limited(0)=1
               parinfo(2).limited(1)=1
               parinfo(2).limits(0)=0.25d0*fwhm_e 
               parinfo(2).limits(1)=2.0d0*fwhm_e 
               
               parinfo(3).fixed=0
               parinfo(3).limited(0)=1
               parinfo(3).limited(1)=1
               parinfo(3).limits(0)=-0.2d0 
               parinfo(3).limits(1)=0.2d0 
                
               param=[pcon_e, rvc_e, fwhm_e, a_e]
                             dp=dblarr(4)
                             dp[*]=1.0d0
                              
                             
; Fixing some parameters controlling convergence and operation of MPFIT

          quiet=1
          maxiter=500


          ftol=1.0d-16 
          xtol=1.0d-16
          gtol=1.0d-16
          autoder=1 ; numerical derivatives are computed by MPFIT.PRO 

; Computing the best fit model with the bi-Gaussian function of Nardetto et al. (2006)

          myfunct1='bi_gaussian'
          parms = MPFIT(MYFUNCT1, param, FUNCTARGS=fa, NFEV=nfev, $
                 MAXITER=maxiter, ERRMSG=errmsg, NPRINT=nprint, QUIET=quiet, $
                 FTOL=ftol, XTOL=xtol, GTOL=gtol, NITER=niter, $
                 STATUS=status, ITERPROC=iterproc, ITERARGS=iterargs, $
                 COVAR=covar, PERROR=perror, BESTNORM=bestnorm, $
                 PARINFO=parinfo,autoderivative=autoder,nocatch=1)

                 if(status eq 0) then begin 
                   print, 'Improper input parameters '
                   print, errmsg
                 endif   
                  
                 print,  ' MPFIT exit STATUS = ', status, ' No. of iterations ', niter 
                 print,  ' Bi-Gaussian fitting: chi square = ', bestnorm
                

; computing the best fit to the normalized CCF 

                 dev=bi_gaussian(parms, dp, xaxisf=xaxisf, ccfnf=ccfnf, measure_errors=measure_errors)
                 ccfn_bigaussian_fit=-dev*measure_errors+ccfnf

; Estimating the errors on the parameters from the covariance matrix and the residuals of the best fits (see description inside mpfit.pro file)

                 sigma_parameters=dblarr(4) 
                 dof=n_elements(ccfnf)-n_elements(parms)
                 sigma_parameters[*]=perror[*]*sqrt(bestnorm/dof)

                 print, 'Reduced chi square ', bestnorm/dof 
                 print, 'PARAMS = ', parms
                 print, 'SIGMA PARAMS = ', sigma_parameters 

                 erase 
                 
             if(operflag eq 1) then begin 
                 set_plot, 'x' 
                 plot, xaxis, ccfn, psym=4, /ynozero, pos=[0.1,0.35,0.9,0.9], /noerase, $ 
                        ytitle='Normalized CCF', title='Best fit to the CCF with a Bi-Gaussian function'
                   
                 err_p=stdev_ccf/max(ccf)
                 errplot, xaxis, ccfn-err_p, ccfn+err_p
                 oplot, xaxisf, ccfn_bigaussian_fit, linestyle=0
                 
                 ydiff=ccfn[flim]-ccfn_bigaussian_fit
                 
                 plot, xaxisf,ydiff, /ynozero, pos=[0.1,0.1,0.9,0.3], linestyle=0, /noerase, $ 
                       xtitle='RV (km/s)', ytitle='Residuals'
                 print, 'Press any key to procede '
                 vzzv=get_kbrd() 
                 erase 
             endif 

; computing the Gaussian (symmetric gaussian) best fit by fixing the asymmetry to zero

          myfunct2='gaussian'
          
                    parinfo(3).fixed=1
                    a_ef=0.0d0
                                        
                    paramg=[pcon_e, rvc_e, fwhm_e, a_ef] 

;                    paramg[3]=a_ef 
                   
          parms_g = MPFIT(MYFUNCT2, paramg, FUNCTARGS=fa, NFEV=nfev, $
                 MAXITER=maxiter, ERRMSG=errmsg, NPRINT=nprint, QUIET=quiet, $
                 FTOL=ftol, XTOL=xtol, GTOL=gtol, NITER=niter, $
                 STATUS=status, ITERPROC=iterproc, ITERARGS=iterargs, $
                 COVAR=covar, PERROR=perror, BESTNORM=bestnorm, $
                 PARINFO=parinfo,autoderivative=autoder,nocatch=1)

                 sigma_parameters_g=dblarr(4) 
                 dof=n_elements(ccfnf)-n_elements(parms_g)
                 sigma_parameters_g[*]=perror[*]*sqrt(bestnorm/dof)

          print, 'Parameters of the Gaussian best fit ', parms_g
          print, 'Symmetric Gaussian fitting: chi square = ', bestnorm
          delta_V_nardetto=parms_g[1]-parms[1]
          print, 'Delta V = ', delta_V_nardetto, ' (km/s)'
          sigma_deltaV=sqrt(sigma_parameters_g[1]^2 + $
                            sigma_parameters[1]^2)
          print, 'Standard deviation of Delta V = ', sigma_deltaV, ' (km/s)'     
          
          fwhm_gaussfit=sqrt(2.0d0)*parms_g[2]
          sigma_fwhm_gaussfit=sqrt(2.0d0)*sigma_parameters_g[2] ; standard deviation of the FWHM of the CCF in km/s from the covariance matrix and the residuals of the best fit          
          
          contrast_gaussfit=parms_g[0]
          sigma_contrast_gaussfit=sigma_parameters_g[0]
          
          print, 'FWHM of the CCF (from the Gaussian best fit) ', fwhm_gaussfit, ' (km/s) with STDEV = ', sigma_fwhm_gaussfit, ' (km/s)' 
          print, 'FWHM of the CCF as given by the DRS ', fwhm_drs, ' (km/s)'
          print, ' '
          print, 'Contrast of the CCF (from the Gaussian best fit) ', contrast_gaussfit, ' with STDEV ', sigma_contrast_gaussfit
          print, 'Contrast of the CCF as given by the DRS ', contrast_drs 
  
;
; numerically evaluating the derivatives of the normalized CCF function using a Savitzy-Golay smoothing filter  
; (see Press et al. 2007, Numerical Recipes, 3rd Ed., Sect. 14.9) 
;

; computing the coefficients of the Savitzy-Golay filter to be passed to v_asymmetry 

          nleft=10
          nright=10
          order_der=1 ; order of the derivative of the CCF vs. RV 
          degree=4 ; degree < 4 for derivatives is not recommended 
          
          sav_coeff=savgol(nleft,nright,order_der,degree,/double)*(factorial(order_der)/(cdelt1^order_der))

; computing the first derivative of the CCF           

         ccfn_deriv=convol(ccfn, sav_coeff, /edge_truncate) 

; computing the Savitzy-Golay coefficients to plot the fit to the CCF 
          
          order=0
          
          sav_coeff_bf=savgol(nleft,nright,order,degree,/double)

; computing the best fit with the S-G fit          

         ccfn_bf_savgol=convol(ccfn, sav_coeff_bf, /edge_truncate) 
          
          nder=n_elements(ccfn_deriv)
          if (nder ne nx ) then begin
            print, 'The number of elements of the derivative vector is not correct: stop'
            stop 
          endif 
  if (operflag eq 1) then begin 
          plot, xaxis, (ccfn_bf_savgol-ccfn), linestyle=0, /ynozero, pos=[0.2, 0.1, 0.9, 0.9], $ 
                xtitle='RV (km/s)', ytitle='CCF fit residuals (normalized)', $ 
                title='Residuals of the fit with the Savitzy-Golay filter to the normalized CCF'
          print, 'Now plotting the difference between CCF and its S-G best fit for spectrum no. ', number_file
          print, 'Press any key to continue '
          vzzv=get_kbrd()

      endif 
          check_sg=max(abs(ccfn_bf_savgol-ccfn)/max(ccfn)) 
          if(check_sg gt 1.0d-03) then begin 
             print, 'The best fit to the CCF based on the Savitzy-Golay filter is not good: ' 
             print, 'the relative difference is ', check_sg
             print, 'Press any key to continue '
             fff=get_kbrd() 
          endif 
          
          res_sg_bf=ccfn-ccfn_bf_savgol ; computing residuals with the best fit of Savitzy-Golay filter 
          
 
;  
; computing V_asy as defined by Figueira et al. (2013, A&A 557, A93) and V_asy_mod and estimating their standard deviations in the case of 
; a Gaussian noise and in the case of a correlated noise by means of the prayer-bead method (see Cowan et al. 2012, ApJ 747, 82; Sect. 4.2);  
; the final standard deviations are obtained by adding the single variances in quadrature, respectively. The same procedure is applied to 
; evaluate the FWHM of the CCF and the \Delta V indicator of Nardetto et al. (2006, A&A 453, 309) and their standard deviations. 
;    

         v_asymmetry, xaxis, ccfn, sav_coeff, sav_coeff_bf, v_asy0_no_noise, v_asy0_mod_no_noise 
     
            nsam=100
            v_asy_trial=dblarr(nsam)
            v_asy_mod_trial=dblarr(nsam)
            v_asy_trial_1=dblarr(nsam) 
            v_asy_mod_trial_1=dblarr(nsam)
            
            delta_v_trial=dblarr(nsam)
            fwhm_trial=dblarr(nsam)
            contrast_trial=dblarr(nsam)
            delta_v_trial_1=dblarr(nsam)
            fwhm_trial_1=dblarr(nsam)
            contrast_trial_1=dblarr(nsam)
                     
            sigma_ccf=sqrt(ccf)/continuum_level_gfit
             
            nccf_elem=n_elements(ccfn)
            seed=1003L ; seed for the random number generator of IDL 
                                 
            for j=0, nsam-1 do begin
              ccfn_trial=ccfn+sigma_ccf*randomn(seed, nccf_elem) ; generating a noisy normalized CCF

              v_asymmetry, xaxis, ccfn_trial, sav_coeff, sav_coeff_bf, v_asy0, v_asy_mod0
              asym_gauss_parameters, xaxis, ccfn_trial, stdev_ccf, delta_v_nardetto0, fwhm_gaussfit0, contrast_gaussfit0
              v_asy_trial[j]=v_asy0
              v_asy_mod_trial[j]=v_asy_mod0
              delta_v_trial[j]=delta_v_nardetto0
              fwhm_trial[j]=fwhm_gaussfit0
              contrast_trial[j]=contrast_gaussfit0

            endfor
            
            nperr=n_elements(res_sg_bf)
            if(nsam gt nperr) then begin 
               print, 'No. of sampling exceeding the number of allowed prayer-bead noise realizations: stop
               stop 
            endif 
            for j=0, nsam-1 do begin
              if(j eq 0) then verror=2.0d0*res_sg_bf
              if(j ge 1 and j lt nperr-1) then begin
                verror=2.0d0*[res_sg_bf[j:(nperr-1)], res_sg_bf[0:(j-1)]]
              endif
              if(j ge nperr-1) then verror=2.0d0*res_sg_bf
              ccfn_trial=ccfn+verror
              v_asymmetry, xaxis, ccfn_trial, sav_coeff, sav_coeff_bf, v_asy0, v_asy_mod0
              asym_gauss_parameters, xaxis, ccfn_trial, stdev_ccf, delta_v_nardetto0, fwhm_gaussfit0, contrast_gaussfit0
              v_asy_trial_1[j]=v_asy0
              v_asy_mod_trial_1[j]=v_asy_mod0
              delta_v_trial_1[j]=delta_v_nardetto0
              fwhm_trial_1[j]=fwhm_gaussfit0
              contrast_trial_1[j]=contrast_gaussfit0
            endfor

             
            
            v_asy=median([v_asy_trial, v_asy_trial_1])
            v_asy_mod=median([v_asy_mod_trial, v_asy_mod_trial_1])
            dev_v_asy=abs(v_asy0_no_noise - v_asy)/(0.5d0*(v_asy+v_asy0_no_noise))
            dev_v_asy_mod=abs(v_asy0_mod_no_noise-v_asy_mod)/(0.5d0*(v_asy_mod+v_asy0_mod_no_noise))
            if(dev_v_asy gt 0.01d0) then begin 
             print, 'The value of V_asy computed from the median of 200 noisy realizations'
             print, 'has a relative deviation of ', dev_v_asy, ' with respect to that computed'
             print, 'directly from the CCF '
             print, 'Do you want to continue ? (Y/N) '
             ftyf=get_kbrd()
             if(ftyf eq 'N') then stop
            endif
            
             if(dev_v_asy_mod gt 0.01d0) then begin 
             print, 'The value of V_asy_mod computed from the median of 200 noisy realizations'
             print, 'has a relative deviation of ', dev_v_asy_mod, ' with respect to that computed'
             print, 'directly from the CCF '
             print, 'Do you want to continue ? (Y/N) '
             ftyf=get_kbrd()
             if(ftyf eq 'N') then stop
            endif
            
;            
; excluding outliers from the sets used to compute the standard deviations of v_asy and v_asy_mod 
;             
            thres=10.0d0
            ccxtrial=where(abs(v_asy_trial) le abs(thres*v_asy))
            ccxtrial_mod=where(abs(v_asy_mod_trial) le abs(thres*v_asy_mod))

             if(n_elements(ccxtrial) lt 0.3d0*nsam) then begin 
               print, 'It is difficult to estimate v_asy because of the noise: stop'
               stop              
             endif
              if(n_elements(ccxtrial_mod) lt 0.3d0*nsam) then begin 
               print, 'It is difficult to estimate v_asy_mod because of the noise: stop'
               stop              
             endif
            
            stdev_v_asy_0=stdev(v_asy_trial[ccxtrial]) 
            stdev_v_asy_1=stdev(v_asy_trial_1)
            
            stdev_v_asy=sqrt(stdev_v_asy_0^2 + stdev_v_asy_1^2) 

            stdev_v_asy_mod_0=stdev(v_asy_mod_trial[ccxtrial_mod]) 
            stdev_v_asy_mod_1=stdev(v_asy_mod_trial_1)
            
            stdev_v_asy_mod=sqrt(stdev_v_asy_mod_0^2 + stdev_v_asy_mod_1^2) 

   print, 'V_asy = ', v_asy , ' (V_asy is non-dimensional for a normalized CCF)'    
   print, 'Standard deviation of V_asy ', stdev_v_asy 
   print, 'V_asy_mod = ',v_asy_mod
   print, 'Standard deviation of V_asy_mod ', stdev_v_asy_mod                     

         
      fwhm_set=median([fwhm_trial, fwhm_trial_1])
      sigma_fwhm_set=sqrt(stdev(fwhm_trial)^2 + stdev(fwhm_trial_1)^2)

      delta_v_nardetto_set= median([delta_v_trial, delta_v_trial_1])
      sigma_delta_v_set=sqrt(stdev(delta_v_trial)^2 + stdev(delta_v_trial_1)^2)

      contrast_set=median([contrast_trial, contrast_trial_1])
      sigma_contrast_set=sqrt(stdev(contrast_trial)^2 + stdev(contrast_trial_1)^2)



                             
                             end 
                             
;*********************************************************************
;
;
;
;
; ********************************************************************************
;

            pro v_asymmetry, xaxis, ccfn_in, sav_coeff, sav_coeff_bf, v_asy, v_asy_mod

;            
; This procedure computes the old V_asy index of Figueira et al. (2013, A&A 557, A93), not recommended for applications, 
; and the new V_asy_mod that is the recommended index. 
; 
; computing the filtered CCF and assuming it as the new CCF 

             ccfn_bf_savgol=convol(ccfn_in, sav_coeff_bf, /edge_truncate) 
             ccfn=ccfn_bf_savgol 

; computing the first derivative of the CCF           

         ccfn_deriv=convol(ccfn_in, sav_coeff, /edge_truncate) 

; computing the v_asy index as defined in Sect. 6 of Figueira et al. (2013, A&A 557, A93) and the 
; new V_asy_mod as defined in Lanza et al. (2018, A&A in press)  

                         ndepth=100 ; number of intervals considered along the full normalized CCF depth
                         
                         min_ccfn=min(ccfn, imin)
                         delta_int=(1.0d0-min_ccfn)/float(ndepth)
                         
                         wei_red=dblarr(ndepth)
                         wei_blue=dblarr(ndepth)
                         
                         wei_red_mod=dblarr(ndepth)
                         wei_blue_mod=dblarr(ndepth)
                         
                         nx=n_elements(xaxis) 
                         
                         for k=1, ndepth-1 do begin
                           
                            flux_level=1.0d0-delta_int*double(k)
                            rv_blue=interpol(xaxis[0:imin], ccfn[0:imin], flux_level, /QUADRATIC)
                            deriv_blue=interpol(ccfn_deriv[0:imin], ccfn[0:imin], flux_level, /quadratic)  
                            wei_blue[k]=((rv_blue*deriv_blue)^2)/flux_level
                            wei_blue_mod[k]=(deriv_blue^2)/flux_level
                            
                            rv_red=interpol(xaxis[imin:(nx-1)], ccfn[imin:(nx-1)], flux_level, /QUADRATIC)
                            deriv_red=interpol(ccfn_deriv[imin:(nx-1)], ccfn[imin:(nx-1)], flux_level, /quadratic)  
                            wei_red[k]=((rv_red*deriv_red)^2)/flux_level
                            wei_red_mod[k]=(deriv_red^2)/flux_level
                                                
                         endfor 
                         
                         wei_mean=0.5d0*(wei_red+wei_blue)
                         wei_mean_mod=0.5d0*(wei_red_mod+wei_blue_mod)
                         nl=5   ; selecting the lower index of the vector
                         nu=95   ; selecting the upper index of the vector 
                         
                         v_asy=total((wei_red[nl:nu]-wei_blue[nl:nu])*wei_mean[nl:nu])/total(wei_mean[nl:nu])
                         v_asy_mod=total((wei_red_mod[nl:nu]-wei_blue_mod[nl:nu])*wei_mean_mod[nl:nu])/total((wei_mean_mod[nl:nu])^2)
                      

 
            end 
;
;
; **************************************************************************************************************
;
;
;
; **************************************************************************************************************
;
pro asym_gauss_parameters, xaxis, ccf, stdev_ccf, delta_v_nardetto, fwhm_gaussfit, contrast_gaussfit

  ; computing a Gaussian and a bi-Gaussian best fit to the normalized CCF below some threshold level;
  ; the IDL function MPFIT.PRO, that applies Levenberg-Marquardt minimization, is used. It should be
  ; compiled in advance from the separate file mpfit.pro

 
  ; estimating the continuum by a Gaussian fit to the CCF 

  nterms=4
  result_gaussfit=gaussfit(xaxis, ccf, a_gauss_fit, nterms=nterms)
  continuum_level_gfit=a_gauss_fit[3]
  ccfn=ccf/continuum_level_gfit

  ; only the part of the normalized CCF below THRESHOLD is considered for these best fits

  threshold=0.95d0*continuum_level_gfit 


  flim=where(ccfn le threshold)

  if(flim[0] eq -1) then begin
    print, 'Not enough points along the CCF below the threshold: STOP'
    print, filename
    stop
  endif

  ccfnf=ccfn[flim]
  xaxisf=xaxis[flim]
  measure_errors=stdev_ccf[flim]/max(ccf)

  nxf=n_elements(flim)



  fa={xaxisf:xaxisf, ccfnf:ccfnf, measure_errors:measure_errors}

  parinfo=replicate({value:0.D, fixed:0, limited:[0,0], limits:[0.D,0.D]},4)

  ; estimating the parameters and fixing their allowed ranges of
  ; variations for a constrained minimization of the chi square

  pcon_e=1.0d0-min(ccfnf,iminf) ; central depth of the Gaussian in relative units

  if(pcon_e lt 0.0d0) then begin
    print, 'Error in the depth of the CCF: STOP'
    stop
  endif

  rvc_e=xaxisf[iminf] ; central RV of the CCF

  ndepth=100 ; number of intervals considered along the full normalized CCF depth

  min_ccfn=min(ccfn, imin)
  delta_int=(1.0d0-min_ccfn)/float(ndepth)

  flux_level=1.0d0-(delta_int*0.5d0*ndepth)
  wave1=interpol(xaxisf[0:iminf], ccfnf[0:iminf], flux_level, /QUADRATIC)
  wave2=interpol(xaxisf[iminf:(nxf-1)], ccfnf[iminf:(nxf-1)], flux_level, /QUADRATIC)

  fwhm_e=double(abs(wave2-wave1)) ; full width at half maximum of the CCF
  a_e=0.01d0 ; asymmetry parameter in the bi-Gaussian function

  parinfo(0).fixed=0
  parinfo(0).limited(0)=1
  parinfo(0).limited(1)=1
  parinfo(0).limits(0)=0.25d0*pcon_e
  parinfo(0).limits(1)=2.0d0*pcon_e

  parinfo(1).fixed=0
  parinfo(1).limited(0)=1
  parinfo(1).limited(1)=1

  if(rvc_e gt 0.0d0) then begin
    rv1=0.25d0*rvc_e
    rv2=2.0d0*rvc_e
  endif else begin
    rv1=2.0d0*rvc_e
    rv2=0.25d0*rvc_e
  endelse

  parinfo(1).limits(0)=rv1
  parinfo(1).limits(1)=rv2

  parinfo(2).fixed=0
  parinfo(2).limited(0)=1
  parinfo(2).limited(1)=1
  parinfo(2).limits(0)=0.25d0*fwhm_e 
  parinfo(2).limits(1)=2.0d0*fwhm_e 

  parinfo(3).fixed=0
  parinfo(3).limited(0)=1
  parinfo(3).limited(1)=1
  parinfo(3).limits(0)=-0.2d0
  parinfo(3).limits(1)=0.2d0

  param=[pcon_e, rvc_e, fwhm_e, a_e]
  dp=dblarr(4)
  dp[*]=1.0d0

  ; Fixing some parameters controlling convergence and operation of MPFIT

  quiet=1
  maxiter=500

  ftol=1.0d-16 
  xtol=1.0d-16
  gtol=1.0d-16
  autoder=1 ; numerical derivatives are computed by MPFIT.PRO

  ; Computing the best fit model with the bi-Gaussian function of Nardetto et al. (2006)

  myfunct1='bi_gaussian'
  parms = MPFIT(MYFUNCT1, param, FUNCTARGS=fa, NFEV=nfev, $
    MAXITER=maxiter, ERRMSG=errmsg, NPRINT=nprint, QUIET=quiet, $
    FTOL=ftol, XTOL=xtol, GTOL=gtol, NITER=niter, $
    STATUS=status, ITERPROC=iterproc, ITERARGS=iterargs, $
    COVAR=covar, PERROR=perror, BESTNORM=bestnorm, $
    PARINFO=parinfo,autoderivative=autoder,nocatch=1)

  if(status eq 0) then begin
    print, 'Improper input parameters '
    print, errmsg
  endif

  ; computing the best fit to the normalized CCF

  dev=bi_gaussian(parms, dp, xaxisf=xaxisf, ccfnf=ccfnf, measure_errors=measure_errors)
  ccfn_bigaussian_fit=-dev*measure_errors+ccfnf

  ; Estimating the errors on the parameters (see description inside mpfit.pro file)

  sigma_parameters=dblarr(4)
  dof=n_elements(ccfnf)-n_elements(parms)
  sigma_parameters[*]=perror[*]*sqrt(bestnorm/dof)

  ;    print, 'Reduced chi square ', bestnorm/dof
  ;    print, 'PARAMS = ', parms
  ;    print, 'SIGMA PARAMS = ', sigma_parameters

  ;    erase
  ;    plot, xaxis, ccfn, psym=4, /ynozero, pos=[0.1,0.35,0.9,0.9], /noerase, $
  ;      ytitle='Normalized CCF', title='Best fit to the CCF with a Bi-Gaussian function'

  err_p=stdev_ccf/max(ccf)
  ;      errplot, xaxis, ccfn-err_p, ccfn+err_p
  ;      oplot, xaxisf, ccfn_bigaussian_fit, linestyle=0

  ;      ydiff=ccfn[flim]-ccfn_bigaussian_fit

  ;      plot, xaxisf,ydiff, /ynozero, pos=[0.1,0.1,0.9,0.3], linestyle=0, /noerase, $
  ;        xtitle='RV (km/s)', ytitle='Residuals'


  ; Computing the Gaussian (symmetric gaussian) best fit by fixing the asymmetry to zero

  myfunct2='gaussian'

  parinfo(3).fixed=1
  a_ef=0.0d0

  paramg=[pcon_e, rvc_e, fwhm_e, a_ef] 
  
;  paramg[3]=a_ef

  parms_g = MPFIT(MYFUNCT2, paramg, FUNCTARGS=fa, NFEV=nfev, $
    MAXITER=maxiter, ERRMSG=errmsg, NPRINT=nprint, QUIET=quiet, $
    FTOL=ftol, XTOL=xtol, GTOL=gtol, NITER=niter, $
    STATUS=status, ITERPROC=iterproc, ITERARGS=iterargs, $
    COVAR=covar, PERROR=perror, BESTNORM=bestnorm, $
    PARINFO=parinfo,autoderivative=autoder,nocatch=1)

  sigma_parameters_g=dblarr(4)
  dof=n_elements(ccfnf)-n_elements(parms_g)
  sigma_parameters_g[*]=perror[*]*sqrt(bestnorm/dof)
  
  delta_V_nardetto=parms_g[1]-parms[1]
 
  sigma_deltaV=sqrt(sigma_parameters_g[1]^2 + $
    sigma_parameters[1]^2)