Newer
Older
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.
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
;
; 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))
Loading full blame...