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) 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 contrast_gaussfit=parms_g[0] sigma_contrast_gaussfit=sigma_parameters_g[0] end ; ********************************************************************************* ; ********************************************************************************* function gaussian, ax, dp, xaxisf=x, ccfnf=y, measure_errors=err pcon=ax[0] rvc=ax[1] fwhm=ax[2] a=ax[3] xaxisf=x left=where(xaxisf gt rvc) right=where(xaxisf le rvc) yfit=dblarr(n_elements(xaxisf)) argl0=2.0d0*alog(2.0d0)*((xaxisf[left]-rvc)^2) argl=argl0/((fwhm*(1.0d0+a))^2) yfit[left]=1.0d0-pcon*exp(-argl) argr0=2.0d0*alog(2.0d0)*((xaxisf[right]-rvc)^2) argr=argr0/((fwhm*(1.0d0-a))^2) yfit[right]=1.0d0-pcon*exp(-argr) if (n_params() gt 1) then begin requested=dp dp=make_array(n_elements(xaxisf), n_elements(ax), /double, value=0.0d0) dp[*,0]=0.0d0 dp[*,1]=0.0d0 dp[*,2]=1.0d0/err[*] dp[*,3]=1.0d0/err[*] endif return, (y-yfit)/err end ; ********************************************************************************* function bi_gaussian, ax, dp, xaxisf=x, ccfnf=y, measure_errors=err pcon=ax[0] rvc=ax[1] fwhm=ax[2] a=ax[3] xaxisf=x left=where(xaxisf ge rvc) right=where(xaxisf le rvc) yfit=dblarr(n_elements(xaxisf)) argl0=2.0d0*alog(2.0d0)*((xaxisf[left]-rvc)^2) argl=argl0/((fwhm*(1.0d0+a))^2) yfit[left]=1.0d0-pcon*exp(-argl) argr0=2.0d0*alog(2.0d0)*((xaxisf[right]-rvc)^2) argr=argr0/((fwhm*(1.0d0-a))^2) yfit[right]=1.0d0-pcon*exp(-argr) if (n_params() gt 1) then begin requested=dp dp=make_array(n_elements(xaxisf), n_elements(ax), /double, value=0.0d0) dp[*,0]=0.0d0 dp[*,1]=0.0d0 dp[*,2]=1.0d0/err[*] dp[*,3]=1.0d0/err[*] endif return, (y-yfit)/err end ;************************************************************************************** ; ; NAME: ; READFITS ; PURPOSE: ; Read a FITS file into IDL data and header variables. ; EXPLANATION: ; READFITS() can read FITS files compressed with gzip or Unix (.Z) ; compression. FPACK ( http://heasarc.gsfc.nasa.gov/fitsio/fpack/ ) ; compressed FITS files can also be read provided that the FPACK software ; is installed. ; See http://idlastro.gsfc.nasa.gov/fitsio.html for other ways of ; reading FITS files with IDL. ; ; CALLING SEQUENCE: ; Result = READFITS( Filename/Fileunit,[ Header, heap, /NOSCALE, EXTEN_NO=, ; NSLICE=, /SILENT , STARTROW =, NUMROW = , HBUFFER=, ; /CHECKSUM, /COMPRESS, /FPACK, /No_Unsigned, NaNVALUE = ] ; ; INPUTS: ; Filename = Scalar string containing the name of the FITS file ; (including extension) to be read. If the filename has ; a *.gz extension, it will be treated as a gzip compressed ; file. If it has a .Z extension, it will be treated as a ; Unix compressed file. If Filename is an empty string then ; the user will be queried for the file name. ; OR ; Fileunit - A scalar integer specifying the unit of an already opened ; FITS file. The unit will remain open after exiting ; READFITS(). There are two possible reasons for choosing ; to specify a unit number rather than a file name: ; (1) For a FITS file with many extensions, one can move to the ; desired extensions with FXPOSIT() and then use READFITS(). This ; is more efficient than repeatedly starting at the beginning of ; the file. ; (2) For reading a FITS file across a Web http: address after opening ; the unit with the SOCKET procedure ; ; OUTPUTS: ; Result = FITS data array constructed from designated record. ; If the specified file was not found, then Result = -1 ; ; OPTIONAL OUTPUT: ; Header = String array containing the header from the FITS file. ; If you don't need the header, then the speed may be improved by ; not supplying this parameter. Note however, that omitting ; the header can imply /NOSCALE, i.e. BSCALE and BZERO values ; may not be applied. ; heap = For extensions, the optional heap area following the main ; data array (e.g. for variable length binary extensions). ; ; OPTIONAL INPUT KEYWORDS: ; /CHECKSUM - If set, then READFITS() will call FITS_TEST_CHECKSUM to ; verify the data integrity if CHECKSUM keywords are present ; in the FITS header. Cannot be used with the NSLICE, NUMROW ; or STARTROW keywords, since verifying the checksum requires ; that all the data be read. See FITS_TEST_CHECKSUM() for more ; information. ; ; /COMPRESS - Signal that the file is gzip compressed. By default, ; READFITS will assume that if the file name extension ends in ; '.gz' then the file is gzip compressed. The /COMPRESS keyword ; is required only if the the gzip compressed file name does not ; end in '.gz' or .ftz ; ; ; EXTEN_NO - non-negative scalar integer specifying the FITS extension to ; read. For example, specify EXTEN = 1 or /EXTEN to read the ; first FITS extension. ; ; /FPACK - Signal that the file is compressed with the FPACK software. ; http://heasarc.gsfc.nasa.gov/fitsio/fpack/ ) By default, ; (READFITS will assume that if the file name extension ends in ; .fz that it is fpack compressed. The FPACK software must ; be installed on the system ; ; HBUFFER - Number of lines in the header, set this to slightly larger ; than the expected number of lines in the FITS header, to ; improve performance when reading very large FITS headers. ; Should be a multiple of 36 -- otherwise it will be modified ; to the next higher multiple of 36. Default is 180 ; ; /NOSCALE - If present and non-zero, then the output data will not be ; scaled using the optional BSCALE and BZERO keywords in the ; FITS header. Default is to scale. ; ; /NO_UNSIGNED - By default, if the header indicates an unsigned integer ; (BITPIX = 16, BZERO=2^15, BSCALE=1) then READFITS() will output ; an IDL unsigned integer data type (UINT). But if /NO_UNSIGNED ; is set, then the data is converted to type LONG. ; ; NSLICE - An integer scalar specifying which N-1 dimensional slice of a ; N-dimensional array to read. For example, if the primary ; image of a file 'wfpc.fits' contains a 800 x 800 x 4 array, ; then ; ; IDL> im = readfits('wfpc.fits',h, nslice=2) ; is equivalent to ; IDL> im = readfits('wfpc.fits',h) ; IDL> im = im[*,*,2] ; but the use of the NSLICE keyword is much more efficient. ; Note that any degenerate dimensions are ignored, so that the ; above code would also work with a 800 x 800 x 4 x 1 array. ; ; NUMROW - Scalar non-negative integer specifying the number of rows ; of the image or table extension to read. Useful when one ; does not want to read the entire image or table. ; ; POINT_LUN - Position (in bytes) in the FITS file at which to start ; reading. Useful if READFITS is called by another procedure ; which needs to directly read a FITS extension. Should ; always be a multiple of 2880, and not be used with EXTEN_NO ; keyword. ; ; /SILENT - Normally, READFITS will display the size the array at the ; terminal. The SILENT keyword will suppress this ; ; STARTROW - Non-negative integer scalar specifying the row ; of the image or extension table at which to begin reading. ; Useful when one does not want to read the entire table. ; ; NaNVALUE - This keyword is included only for backwards compatibility ; with routines that require IEEE "not a number" values to be ; converted to a regular value. ; ; /UNIXPIPE - When a FileUnit is supplied to READFITS(), then /UNIXPIPE ; indicates that the unit is to a Unix pipe, and that ; no automatic byte swapping is performed. ; ; EXAMPLE: ; Read a FITS file test.fits into an IDL image array, IM and FITS ; header array, H. Do not scale the data with BSCALE and BZERO. ; ; IDL> im = READFITS( 'test.fits', h, /NOSCALE) ; ; If the file contains a FITS extension, it could be read with ; ; IDL> tab = READFITS( 'test.fits', htab, /EXTEN ) ; ; The function TBGET() can be used for further processing of a binary ; table, and FTGET() for an ASCII table. ; To read only rows 100-149 of the FITS extension, ; ; IDL> tab = READFITS( 'test.fits', htab, /EXTEN, ; STARTR=100, NUMR = 50 ) ; ; To read in a file that has been compressed: ; ; IDL> tab = READFITS('test.fits.gz',h) ; ; ERROR HANDLING: ; If an error is encountered reading the FITS file, then ; (1) the system variable !ERROR_STATE.CODE is set negative ; (via the MESSAGE facility) ; (2) the error message is displayed (unless /SILENT is set), ; and the message is also stored in !!ERROR_STATE.MSG ; (3) READFITS returns with a value of -1 ; RESTRICTIONS: ; (1) Cannot handle random group FITS ; ; NOTES: ; (1) If data is stored as integer (BITPIX = 16 or 32), and BSCALE ; and/or BZERO keywords are present, then the output array is scaled to ; floating point (unless /NOSCALE is present) using the values of BSCALE ; and BZERO. In the header, the values of BSCALE and BZERO are then ; reset to 1. and 0., while the original values are written into the ; new keywords O_BSCALE and O_BZERO. If the BLANK keyword was ; present (giving the value of undefined elements *prior* to the ; application of BZERO and BSCALE) then the *keyword* value will be ; updated with the values of BZERO and BSCALE. ; ; (2) The use of the NSLICE keyword is incompatible with the NUMROW ; or STARTROW keywords. ; ; (3) On some Unix shells, one may get a "Broken pipe" message if reading ; a Unix compressed (.Z) file, and not reading to the end of the file ; (i.e. the decompression has not gone to completion). This is an ; informative message only, and should not affect the output of READFITS. ; PROCEDURES USED: ; Functions: SXPAR() ; Procedures: MRD_SKIP, SXADDPAR, SXDELPAR ; ; MODIFICATION HISTORY: ; Original Version written in 1988, W.B. Landsman Raytheon STX ; Revision History prior to October 1998 removed ; Major rewrite to eliminate recursive calls when reading extensions ; W.B. Landsman Raytheon STX October 1998 ; Add /binary modifier needed for Windows W. Landsman April 1999 ; Read unsigned datatypes, added /no_unsigned W. Landsman December 1999 ; Output BZERO = 0 for unsigned data types W. Landsman January 2000 ; Update to V5.3 (see notes) W. Landsman February 2000 ; Fixed logic error in use of NSLICE keyword W. Landsman March 2000 ; Fixed byte swapping for Unix compress files on little endian machines ; W. Landsman April 2000 ; Added COMPRESS keyword, catch IO errors W. Landsman September 2000 ; Option to read a unit number rather than file name W.L October 2001 ; Fix undefined variable problem if unit number supplied W.L. August 2002 ; Don't read entire header unless needed W. Landsman Jan. 2003 ; Added HBUFFER keyword W. Landsman Feb. 2003 ; Added CHECKSUM keyword W. Landsman May 2003 ; Restored NaNVALUE keyword for backwards compatibility, ; William Thompson, 16-Aug-2004, GSFC ; Recognize .ftz extension as compressed W. Landsman September 2004 ; Fix unsigned integer problem introduced Sep 2004 W. Landsman Feb 2005 ; Don't modify header for unsigned integers, preserve double precision ; BSCALE value W. Landsman March 2006 ; Use gzip instead of compress for Unix compress files W.Landsman Sep 2006 ; Call MRD_SKIP to skip bytes on different file types W. Landsman Oct 2006 ; Make ndata 64bit for very large files E. Hivon/W. Landsman May 2007 ; Fixed bug introduced March 2006 in applying Bzero C. Magri/W.L. Aug 2007 ; Check possible 32bit overflow when using NSKIP W. Landsman Mar 2008 ; Always reset BSCALE, BZERO even for unsigned integers W. Landsman May 2008 ; Make ndata 64bit for very large extensions J. Schou/W. Landsman Jan 2009 ; Use PRODUCT() to compute # of data points W. Landsman May 2009 ; Read FPACK compressed file via UNIX pipe. W. Landsman May 2009 ; Fix error using NUMROW,STARTROW with non-byte data, allow these ; keywords to be used with primary array W. Landsman July 2009 ; Ignore degenerate trailing dimensions with NSLICE keyword W.L. Oct 2009 ; Add DIALOG_PICKFILE() if filename is an empty string W.L. Apr 2010 ; Set BLANK values *before* applying BSCALE,BZERO, use short-circuit ; operators W.L. May 2010 ; Skip extra SPAWN with FPACK decompress J. Eastman, W.L. July 2010 ; Fix possible problem when startrow=0 supplied J. Eastman/W.L. Aug 2010 ; First header is not necessarily primary if unit supplied WL Jan 2011 ;- function READFITS, filename, header, heap, CHECKSUM=checksum, $ COMPRESS = compress, HBUFFER=hbuf, EXTEN_NO = exten_no, $ NOSCALE = noscale, NSLICE = nslice, $ NO_UNSIGNED = no_unsigned, NUMROW = numrow, $ POINTLUN = pointlun, SILENT = silent, STARTROW = startrow, $ NaNvalue = NaNvalue, FPACK = fpack, UNIXpipe=unixpipe On_error,2 ;Return to user compile_opt idl2 On_IOerror, BAD ; Check for filename input if N_params() LT 1 then begin print,'Syntax - im = READFITS( filename, [ h, heap, /NOSCALE, /SILENT,' print,' EXTEN_NO =, STARTROW = , NUMROW=, NSLICE = ,' print,' HBUFFER = ,/NO_UNSIGNED, /CHECKSUM, /COMPRESS]' return, -1 endif unitsupplied = size(filename,/TNAME) NE 'STRING' ; Set default keyword values silent = keyword_set( SILENT ) do_checksum = keyword_set( CHECKSUM ) if N_elements(exten_no) EQ 0 then exten_no = 0 ; Check if this is a Unix compressed file. (gzip files are handled ; separately using the /compress keyword to OPENR). if N_elements(unixpipe) EQ 0 then unixpipe = 0 if unitsupplied then unit = filename else begin len = strlen(filename) if len EQ 0 then begin filename =dialog_pickfile(filter=['*.fit*;*.fts*;*.img*'], $ title='Please select a FITS file',/must_exist) len = strlen(filename) endif ext = strlowcase(strmid(filename,len-3,3)) gzip = (ext EQ '.gz') || (ext EQ 'ftz') compress = keyword_set(compress) || gzip[0] unixZ = (strmid(filename, len-2, 2) EQ '.Z') fcompress = keyword_set(fpack) || ( ext EQ '.fz') unixpipe = unixZ || fcompress ; Go to the start of the file. openr, unit, filename, ERROR=error,/get_lun, $ COMPRESS = compress, /swap_if_little_endian if error NE 0 then begin message,/con,' ERROR - Unable to locate file ' + filename return, -1 endif ; Handle Unix or Fpack compressed files which will be opened via a pipe using ; the SPAWN command. if unixZ then begin free_lun, unit spawn, 'gzip -cd '+filename, unit=unit gzip = 1b endif else if fcompress then begin free_lun, unit spawn,'funpack -S ' + filename, unit=unit,/sh if eof(unit) then begin message,'Error spawning FPACK decompression',/CON free_lun,unit return,-1 endif endif endelse if keyword_set(POINTLUN) then mrd_skip, unit, pointlun doheader = arg_present(header) or do_checksum if doheader then begin if N_elements(hbuf) EQ 0 then hbuf = 180 else begin remain = hbuf mod 36 if remain GT 0 then hbuf = hbuf + 36-remain endelse endif else hbuf = 36 for ext = 0L, exten_no do begin ; Read the next header, and get the number of bytes taken up by the data. block = string(replicate(32b,80,36)) w = [-1] if ((ext EQ exten_no) && (doheader)) then header = strarr(hbuf) $ else header = strarr(36) headerblock = 0L i = 0L while w[0] EQ -1 do begin if EOF(unit) then begin message,/ CON, $ 'EOF encountered attempting to read extension ' + strtrim(ext,2) if ~unitsupplied then free_lun,unit return,-1 endif readu, unit, block headerblock = headerblock + 1 w = where(strlen(block) NE 80, Nbad) if (Nbad GT 0) then begin message,'Warning-Invalid characters in header',/INF,NoPrint=Silent block[w] = string(replicate(32b, 80)) endif w = where(strcmp(block,'END ',8), Nend) if (headerblock EQ 1) || ((ext EQ exten_no) && (doheader)) then begin if Nend GT 0 then begin if headerblock EQ 1 then header = block[0:w[0]] $ else header = [header[0:i-1],block[0:w[0]]] endif else begin header[i] = block i = i+36 if i mod hbuf EQ 0 then $ header = [header,strarr(hbuf)] endelse endif endwhile if (ext EQ 0 ) && ~(keyword_set(pointlun) || unitsupplied ) then $ if strmid( header[0], 0, 8) NE 'SIMPLE ' then begin message,/CON, $ 'ERROR - Header does not contain required SIMPLE keyword' if ~unitsupplied then free_lun, unit return, -1 endif ; Get parameters that determine size of data region. bitpix = sxpar(header,'BITPIX') byte_elem = abs(bitpix)/8 ;Bytes per element naxis = sxpar(header,'NAXIS') gcount = sxpar(header,'GCOUNT') > 1 pcount = sxpar(header,'PCOUNT') if naxis GT 0 then begin dims = sxpar( header,'NAXIS*') ;Read dimensions ndata = product(dims,/integer) endif else ndata = 0 nbytes = byte_elem * gcount * (pcount + ndata) ; Move to the next extension header in the file. Use MRD_SKIP to skip with ; fastest available method (POINT_LUN or readu) for different file ; types (regular, compressed, Unix pipe, socket) if ext LT exten_no then begin nrec = long((nbytes + 2879) / 2880) if nrec GT 0 then mrd_skip, unit, nrec*2880L endif endfor case BITPIX of 8: IDL_type = 1 ; Byte 16: IDL_type = 2 ; Integer*2 32: IDL_type = 3 ; Integer*4 64: IDL_type = 14 ; Integer*8 -32: IDL_type = 4 ; Real*4 -64: IDL_type = 5 ; Real*8 else: begin message,/CON, 'ERROR - Illegal value of BITPIX (= ' + $ strtrim(bitpix,2) + ') in FITS header' if ~unitsupplied then free_lun,unit return, -1 end endcase if nbytes EQ 0 then begin if ~SILENT then message, $ "FITS header has NAXIS or NAXISi = 0, no data array read",/CON if do_checksum then begin result = FITS_TEST_CHECKSUM(header, data, ERRMSG = errmsg) if ~SILENT then begin case result of 1: message,/INF,'CHECKSUM keyword in header is verified' -1: message,/CON, errmsg else: endcase endif endif if ~unitsupplied then free_lun, unit return,-1 endif ; Check for FITS extensions, GROUPS groups = sxpar( header, 'GROUPS' ) if groups then message,NoPrint=Silent, $ 'WARNING - FITS file contains random GROUPS', /INF ; If an extension, did user specify row to start reading, or number of rows ; to read? if N_elements(STARTROW) EQ 0 then startrow = 0 ;updated Aug 2010 if naxis GE 2 then nrow = dims[1] else nrow = ndata if ~keyword_set(NUMROW) then numrow = nrow if do_checksum then if ((startrow GT 0) || $ (numrow LT nrow) || (N_elements(nslice) GT 0)) then begin message,/CON, $ 'Warning - CHECKSUM not applied when STARTROW, NUMROW or NSLICE is set' do_checksum = 0 endif if exten_no GT 0 then begin xtension = strtrim( sxpar( header, 'XTENSION' , Count = N_ext),2) if N_ext EQ 0 then message, /INF, NoPRINT = Silent, $ 'WARNING - Header missing XTENSION keyword' endif if ((startrow NE 0) || (numrow NE nrow)) then begin if startrow GE dims[1] then begin message,'ERROR - Specified starting row ' + strtrim(startrow,2) + $ ' but only ' + strtrim(dims[1],2) + ' rows in extension',/CON if ~unitsupplied then free_lun,unit return,-1 endif dims[1] = dims[1] - startrow dims[1] = dims[1] < numrow sxaddpar, header, 'NAXIS2', dims[1] if startrow GT 0 then mrd_skip, unit, byte_elem*startrow*dims[0] endif else if (N_elements(NSLICE) EQ 1) then begin ldim = naxis-1 lastdim = dims[ldim] while lastdim EQ 1 do begin ldim = ldim-1 lastdim = dims[ldim] endwhile if nslice GE lastdim then begin message,/CON, $ 'ERROR - Value of NSLICE must be less than ' + strtrim(lastdim,2) if ~unitsupplied then free_lun, unit return, -1 endif dims = dims[0:ldim-1] for i = ldim,naxis-1 do sxdelpar,header,'NAXIS' + strtrim(i+1,2) naxis = ldim sxaddpar,header,'NAXIS' + strtrim(ldim,2),1 ndata = ndata/lastdim nskip = long64(nslice)*ndata*byte_elem if Ndata GT 0 then mrd_skip, unit, nskip endif if ~SILENT then begin ;Print size of array being read if exten_no GT 0 then message, $ 'Reading FITS extension of type ' + xtension, /INF if N_elements(dims) EQ 1 then $ st = 'Now reading ' + strtrim(dims,2) + ' element vector' else $ st = 'Now reading ' + strjoin(strtrim(dims,2),' by ') + ' array' if (exten_no GT 0) && (pcount GT 0) then st = st + ' + heap area' message,/INF,st endif ; Read Data in a single I/O call. Only need byteswapping for data read with ; bidirectional pipe. data = make_array( DIM = dims, TYPE = IDL_type, /NOZERO) readu, unit, data if unixpipe then swap_endian_inplace,data,/swap_if_little if (exten_no GT 0) && (pcount GT 0) then begin theap = sxpar(header,'THEAP') skip = theap - N_elements(data) if skip GT 0 then begin temp = bytarr(skip,/nozero) readu, unit, skip endif heap = bytarr(pcount*gcount*byte_elem) readu, unit, heap if do_checksum then $ result = fits_test_checksum(header,[data,heap],ERRMSG=errmsg) endif else if do_checksum then $ result = fits_test_checksum(header, data, ERRMSG = errmsg) if ~unitsupplied then free_lun, unit if do_checksum then if ~SILENT then begin case result of 1: message,/INF,'CHECKSUM keyword in header is verified' -1: message,/CON, 'CHECKSUM ERROR! ' + errmsg else: endcase endif ; Scale data unless it is an extension, or /NOSCALE is set ; Use "TEMPORARY" function to speed processing. do_scale = ~keyword_set( NOSCALE ) if (do_scale && (exten_no GT 0)) then do_scale = xtension EQ 'IMAGE' if do_scale then begin if bitpix GT 0 then $ blank = sxpar( header, 'BLANK', Count = N_blank) $ else N_blank = 0 Bscale = sxpar( header, 'BSCALE' , Count = N_bscale) Bzero = sxpar(header, 'BZERO', Count = N_Bzero ) if (N_blank GT 0) && ((N_bscale GT 0) || (N_Bzero GT 0)) then $ sxaddpar,header,'O_BLANK',blank,' Original BLANK value' ; Check for unsigned integer (BZERO = 2^15) or unsigned long (BZERO = 2^31) if ~keyword_set(No_Unsigned) then begin no_bscale = (Bscale EQ 1) || (N_bscale EQ 0) unsgn_int = (bitpix EQ 16) && (Bzero EQ 32768) && no_bscale unsgn_lng = (bitpix EQ 32) && (Bzero EQ 2147483648) && no_bscale unsgn = unsgn_int || unsgn_lng endif else unsgn = 0 if unsgn then begin if unsgn_int then begin data = uint(data) - 32768US if N_blank then blank = uint(blank) - 32768US endif else begin data = ulong(data) - 2147483648UL if N_blank then blank = ulong(blank) - 2147483648UL endelse if N_blank then sxaddpar,header,'BLANK',blank sxaddpar, header, 'BZERO', 0 sxaddpar, header, 'O_BZERO', Bzero,' Original BZERO Value' endif else begin if N_Bscale GT 0 then $ if ( Bscale NE 1. ) then begin if size(Bscale,/TNAME) NE 'DOUBLE' then $ data *= float(Bscale) else $ data *= Bscale if N_blank then blank *= bscale sxaddpar, header, 'BSCALE', 1. sxaddpar, header, 'O_BSCALE', Bscale,' Original BSCALE Value' endif if N_Bzero GT 0 then $ if (Bzero NE 0) then begin if size(Bzero,/TNAME) NE 'DOUBLE' then $ data += float(Bzero) else $ ;Fixed Aug 07 data += Bzero if N_blank then blank += bzero sxaddpar, header, 'BZERO', 0. sxaddpar, header, 'O_BZERO', Bzero,' Original BZERO Value' endif endelse if N_blank then sxaddpar,header,'BLANK',blank endif ; Return array. If necessary, first convert NaN values. if n_elements(nanvalue) eq 1 then begin w = where(finite(data,/nan),count) if count gt 0 then data[w] = nanvalue endif return, data ; Come here if there was an IO_ERROR BAD: print,!ERROR_STATE.MSG if (~unitsupplied) && (N_elements(unit) GT 0) then free_lun, unit if N_elements(data) GT 0 then return,data else return, -1 end