.r Venus_combine_data.pro ; ; ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; **************************** SECTION 1 - GAUSSIAN FUNCITONS *********************************************** ; ; Gaussian functions defined and called on in the code for fitting. ; Two functions are calculated, one for a blended Gaussian ; consisting of two Gaussians, and another of a single Gaussian ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;------------------------------------------------------------------------------------------ ; function gfunct ; define Gaussian function for deblend case. ; This will be called later on in the code and the variable F will be ; returned. ; A[] array is defined near the bottom of the code in: ; parinfo[*].value = [2.31e-13, 5577.01, 0.079, 2.11e-12, 5577.232]; **** UPDATE VALUES ***** ; ; A[0] = Amplitude of Venusian green line ; A[1] = central wavelength of of Venusian green line ; A[2] = sigma (spread in the data, defined by the dispersion of the instrument, same for both lines) ; A[3] = Amplitude of the Terrestrial green line ; A[4] = central wavelength of the Terrestrial green line ; ;------------------------------------------------------------------------------------------ ; function gfunct,lambda, A, F, pder ;define Gaussian function for deblend case ; F=A[0]*exp(-(lambda-A[1])^2.0/(2.0*(A[2]^2.0)))+A[3]*exp(-(lambda-A[4])^2.0/(2.0*(A[2]^2.0))) ; IF N_PARAMS() GE 4 THEN $ ; pder = [[exp(-(lambda-A[1])^2/(2*(A[2]^2)))], [(lambda-A[1])*A[0]/A[2]^2*exp(-(lambda-A[1])^2/(2*(A[2]^2)))], $ ; [(lambda-A[1])*A[0]/A[2]^3*exp(-(lambda-A[1])^2/(2*(A[2]^2)))]] ; return, F ;END ;------------------------------------------------------------------------------------------ ; function gfunct2 ; define Gaussian function for individual line ;------------------------------------------------------------------------------------------ ; function gfunct2,lambda, A, F, pder ;define Gaussian function for individual line ; F=A[0]*exp(-(lambda-A[1])^2.0/(2.0*(A[2]^2.0))) ; IF N_PARAMS() GE 4 THEN $ ; pder = [[exp(-(lambda-A[1])^2/(2*(A[2]^2)))], [(lambda-A[1])*A[0]/A[2]^2*exp(-(lambda-A[1])^2/(2*(A[2]^2)))],$ ; [(lambda-A[1])*A[0]/A[2]^3*exp(-(lambda-A[1])^2/(2*(A[2]^2)))]] ; return, F ; END ;------------------------------------------------------------------------------------------ ; Open file using openu, for update, so the data from the previous files are not overwritten. ; Data written will include the peak of the Venusian green line and how many standard ; deviations it is away from the mean in the scatter around the cointuum. Idealy, the peak ; of the Venus green line would be 3-sigma above the mean, but 1 is ok.Print data to file ;------------------------------------------------------------------------------------------ ;openw, lun_1, 'output_continuum_mean_error.dat', /GET_LUN, width = 250, /append ;append to data is added to end of file ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************** SECTION 2 - READING IN DATA *********************************************** ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;READCOL, 'combined_files.dat', filename, index, position, month, day, year, airmass, exposure, relative_vel, $ ; FORMAT='A,A,A,A,I,F,F,F' free_lun, lun READCOL, 'combine_images_inputs.dat', index, position, month, day, year, airmass, exposure, relative_vel, $ FORMAT='A,A,A,A,I,F,F,F' READCOL, 'combine_images_files.dat', filename, FORMAT ='A' READCOL, 'wavelength_combined_input.dat', lambda_true, lambda_obs, dispersion, FORMAT = 'F,F,F' READCOL, 'gaussian_input.dat', gaussian_flux_terr, gaussian_flux_venus, gaussian_ray_terr, $ gaussian_ray_venus, gaussian_norm_terr,gaussian_norm_venus, $ FORMAT='F,F,F,F,F,F' sum_flux = 0 sum_ray = 0 sum_norm = 0 sum_flux_err = 0 sum_ray_err = 0 sum_norm_err = 0 angstrom = '!6!sA!r!u!9 %!6!n' ;define the angstrom character LOADCT,39,/SILENT !P.CHARTHICK=2.5 !P.POSITION=[.17,.12,.87,.95] ;set lower left and upper right corners of plotting box !X.CHARSIZE = 1.3 !Y.CHARSIZE = 1.3 !X.THICK=4 !Y.THICK=4 SET_PLOT,'PS' lambda_plot_max = lambda_true + 0.7 lambda_plot_min = lambda_true - 0.7 print, '---------------------start loop---------------' ;------------------------------------------------------------------------------------------- ; Loop over how many wavelengths there are to analzye. ;-------------------------------------------------------------------------------------------- FOR i = 0, N_ELEMENTS(labmda_true)-1 DO BEGIN print, i ;------------------------------------------------------------------------------------------- ; Going through all the output files, group the flux, rayleigh, and norm files for each ; emission line. ;-------------------------------------------------------------------------------------------- FOR k = 0, N_ELEMENTS(filename)-1 DO BEGIN print, k IF (STRMATCH(filename[k], '*flux*') EQ 1) AND (STRMATCH(filename[k], '*5577*') EQ 1) THEN BEGIN print, filename[k] ; IF (STRMATCH(filename[k], '*5577*') EQ 1) THEN BEGIN ;FOR j = 0, N_ELEMENTS(index)-1 DO BEGIN ; IF (STRMATCH(filename[k], file) EQ 1) THEN BEGIN ; READCOL, filename[k], lamda, flux, error, FORMAT = 'F,F,F' ; print, filename[k] ; print, j, k ; ENDIF ; ENDFOR READCOL, filename[k], lamda, flux, error, FORMAT = 'F,F,F' sum_flux = sum_flux + flux flux_err = flux_err + error^2.0 lambda_flux = lambda ; ENDIF ENDIF IF (STRMATCH(filename[k], '*rayl*') EQ 1) AND (STRMATCH(filename[k], '*5577*') EQ 1) THEN BEGIN READCOL, filename[k], lamda, flux, error, FORMAT = 'F,F,F' sum_ray = sum_ray + flux ray_err = ray_err + error^2.0 lambda_ray = lambda ENDIF IF (STRMATCH(filename[k], '*norm*') EQ 1) AND (STRMATCH(filename[k], '*5577*') EQ 1) THEN BEGIN READCOL, filename[k], lamda, flux, error, FORMAT = 'F,F,F' sum_norm = sum_norm + flux norm_err = norm_err + error^2.0 lambda_norm = lambda ENDIF ENDFOR avg_flux = sum_flux/(n_elements(index)) avg_ray = sum_flux/(n_elements(index)) avg_norm = sum_flux/(n_elements(index)) flux_min = MIN(avg_flux) flux_max = MAX(avg_flux) ray_min = MIN(avg_ray) ray_max = MAX(avg_ray) norm_min = MIN(avg_norm) norm_max = MAX(avg_norm) flux_min_err = MIN(flux_err) flux_max_err = MAX(flux_err) ray_min_err = MIN(ray_err) ray_max_err = MAX(ray_flux_err) norm_min_err = MIN(norm_flux_err) norm_max_err = MAX(norm_flux_err) avg_airmass = MEAN(airmass) ;----------------------------------------------------------------------------- ; Upper and lower wavelength bracket (in Angstroms) around emission feature; ; cant be too large or it will grab terrestrail ; Calculate doppler shift of emission line based on OBSERVED terrestrial ; emission feature ; Set the Veusian emission line wavelength with Venus_lambda_guess ; Set upper and lower wavelength bracket around emission line (cant be too large or ; it will grab terrestrail) ;----------------------------------------------------------------------------- feature_upper_bound = 0.07 feature_lower_bound = 0.07 doppler_shift = relative_vel[0] * lambda_true[i] / 3e5 Venus_lambda_guess = lambda_obs[i]+ doppler_shift Venus_lambda_upper = lambda_obs[i]+ doppler_shift + feature_upper_bound Venus_lambda_lower = lambda_obs[i]+ doppler_shift - feature_lower_bound terr_lambda_upper = lambda_obs[i] + feature_upper_bound terr_lambda_lower = lambda_obs[i] - feature_lower_bound ;----------------------------------------------------------------------------------------------------------- ; Find the max value for the terrestrial and the venusian peaks in ; flux, rayleigh, and normalized units and their wavelength ;----------------------------------------------------------------------------------------------------------- linerange = WHERE((lambda GE terr_lambda_lower) AND (lambda LE terr_lambda_upper)) ;Specify wavelength range we will analyze terr_flux_max = MAX(avg_flux(linerange), pt) terr_flux_max_err = flux_err(linerange(pt)) terr_flux_lambda_max = lambda_flux(linerange(pt)) linerange = WHERE((lambda GE Venus_lambda_lower) AND (lambda LE Venus_lambda_upper)) ;Specify wavelength range we will analyze Venus_flux_max = MAX(avg_flux(linerange), pt) Venus_flux_max_err = flux_err(linerange(pt)) Venus_flux_lambda_max = lambda_flux(linerange(pt)) IF Venus_lambda_guess LT lambda_obs[i] THEN BEGIN ; Set lower wavelength bounds to be Venus if lambda_lower = Venus_lambda_guess ; Venus is blue shifted of terrestrial lambda_upper = terr_flux_lambda_max ENDIF IF Venus_lambda_guess GT lambda_obs[i] THEN BEGIN ; Set lower wavelength bounds to be terr if lambda_lower = equator_terr_lambda_max ; Venus is red shifted of terrestrial lambda_upper = Venus_flux_lambda_max ENDIF lambda_continuum_lower = lambda_lower - 0.15 lambda_continuum_upper = lambda_upper + 0.15 linerange_continuum = WHERE ( (lambda_flux LE (lambda_continuum_lower)) OR (lambda_flux GE (lambda_continuum_upper))) mean_continuum_flux = MEAN((avg_flux(linerange_continuum))) ;calculate mean of continuum mean_continuum_ray = MEAN((avg_ray(linerange_continuum))) mean_continuum_norm = MEAN((avg_norm(linerange_continuum))) sigma_continuum_flux = stddev(avg_flux(linerange_continuum)) ; calculate standard deviation of the sigma_continuum_ray = stddev(avg_rayleigh(linerange_continuum)) ; mean of the continuum sigma_continuum_norm = stddev(avg_norm(linerange_continuum)) ;-------------------------------------------------------------------------------------------- ; Calculate how many standard deviations the greenline is above the average of the continuum ;-------------------------------------------------------------------------------------------- Venus_greenline_confidence_flux = (Venus_flux_max - mean_continuum_flux) / sigma_continuum_flux Venus_greenline_confidence_ray = (Venus_ray_max - mean_continuum_ray) / sigma_continuum_ray Venus_greenline_confidence_norm = (Venus_norm_max - mean_continuum_norm) / sigma_continuum_norm ;-------------------------------------------------------------------------------------------- ; Calculate the upper and lower bounds of the standard deviation above the mean so as to ; plot the deviations on the plot. ;-------------------------------------------------------------------------------------------- sigma_flux_continuum_upper = mean_continuum_flux + sigma_continuum_flux sigma_flux_continuum_lower = mean_continuum_flux - sigma_continuum_flux sigma_ray_continuum_upper = mean_continuum_ray + sigma_continuum_ray sigma_ray_continuum_lower = mean_continuum_ray - sigma_continuum_ray sigma_norm_continuum_upper = mean_continuum_norm + sigma_continuum_norm sigma_norm_continuum_lower = mean_continuum_norm - sigma_continuum_norm ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************************** PLOT DATA ***************************************** ; ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;============================================================================================================ ; ---------------------------------------- FLUX PLOT -------------------------------------------------------- ;============================================================================================================ PRINT,'' PRINT, '' print, 'The following calculations will be done in cgs flux units for the equator' PRINT,'' PRINT, '' ;------------------------------------------------------------------------------------------ ; define parameters to be fit THIS WAS CHANGED A LOT!!!! ; in the last parinfo, set the guessed flux, Venusian greenline ; position, dispersion of the instrument at that wavelength (in Angstroms), guessed ; telluric flux, position, and dispersion at that wavelength ; ; In the updated version of this code that Kenza altered, parinfo was ; changed from fixed values to "limits" and an initial guess was ; made. This means that for each parameter (there are 5, the flux of ; Venusian, wavelength of Venusian, dispersion, flux of telluric, and ; wavelength of telluric), a rough guess is made on each, just to ; start somewhere, and then the fit is run across a range which is ; defined in you parinfo[#].limits. ; ; For example, the range of the first parameter (peak flux of Venusian ; line) is set to 2.13e-13, and the range is set to 1e-20 through ; 1e-10. The numbers need to be defined as doubles (IDL defaults to ; float) or else not enough values past the decimal will be considered. ; parinfo = replicate({value:0.D, fixed:0, limited:[1,1], limits:[0.,1e6]},5) ; parinfo = replicate({value:0.D, fixed:1, limited:[0,0], limits:[0.0,0.0]},6) ; parinfo[0].fixed = 0 Venusian guessed flux ; parinfo[1].fixed = 0 Venusian line position ; parinfo[3].fixed = 0 Dispersion at that order (0.0709 for 47th order of echelle) ; parinfo[3].fixed = 0 Terrestrial wavelength ; parinfo[4].fixed = 0 Terrestrial flux ; parinfo[*].value = [2.31e-13, 5577.012, 0.079, 2.11e-12, 5577.232]; ; parinfo[*].value = [4.5e-13, 5577.01, 0.0709, 2.3e-12, 5577.23, 0.0709]; ; parinfo[*].value = [2e-14, 5577.13, 0.0709, 9e-14, 5577.31, 0.0709]; ;------------------------------------------------------------------------------------------ ;Set limits for your parameters and give an initial guess parinfo = replicate({fixed:0, limited:[1,1], limits:[1e-25,1e6]},5) parinfo[0].limits = [1e-25,1e-08] ; Venusian Flux parinfo[1].limits = [DOUBLE(Venus_lambda_lower[0]), DOUBLE(Venus_lambda_upper[0])] parinfo[2].limits = [dispersion[i]-0.001, dispersion[i]+0.001] parinfo[3].limits = [1e-25,1e-8] ; Terrestrial Flux parinfo[4].limits = [DOUBLE(lambda_obs[k] - 0.1), DOUBLE(lambda_obs[k] + 0.1)] start = [DOUBLE(gaussian_flux_venus[0]), DOUBLE(Venus_lambda_guess), DOUBLE(dispersion[i]), DOUBLE(gaussian_flux_terr[0]), DOUBLE(lambda_obs[i])] ;INITIAL GUESSES, Flux ; start = [DOUBLE(equator_venus_flux_max), DOUBLE(Venus_lambda_guess[0]), 0.079D, DOUBLE(equator_terr_flux_max), 5577.45D] ;INITIAL GUESSES P=fltarr(5,1);5 is number of parameters, 1 is number of spectra to be analyzed perror=fltarr(5,1) ;------------------------------------------------------------------------------------------- ; this used to have 2048 instead of N_ELEMENTS but since the array was ; changed to just be the wavelength range of the Doppler shifted ; Venusian and telluric green line, and I may slightly alter that range ; values for other data sets, the number of elements was just written ; as the same number of Venus_solsub (remember that the number of ; element for Venus_solsub was changed right before this section ;------------------------------------------------------------------------------------------- fluxfit = fltarr(N_ELEMENTS(avg_flux),1) ;this used to have 2048 instead of N_ELEMENTS fluxfittell = fltarr(N_ELEMENTS(avg_flux),1) residual = fltarr(N_ELEMENTS(avg_flux),1) totresidual = fltarr(1) chisquared = fltarr(1) redchi = fltarr(1) averesid = fltarr(1) ew5577 = fltarr(1) ewerr5577 = fltarr(1) DEVICE,FILENAME=strcompress(string(day[0])+month[0]+year[0]+lambda_true[i]+'_combined_flux_FINAL_errors.ps',/remove_all),$ /COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES ;------------------------------------------------------------------------------------------------------------- ;--------------------------------------- BEGIN GAUSSIAN FIT FOR PLOT ----------------------------------------- ;------------------------------------------------------------------------------------------------------------- P[*] = MPFITFUN('gfunct', lambda_flux, avg_flux[*], flux_err[*]/flux_err[*] * 1e-18, start, $ PARINFO=parinfo, NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1) plot, lambda_flux, avg_flux[*], xrange=[lambda_plot_min, lambda_plot_max], /xstyle, $ yrange=[flux_min - flux_min_err*1.2, flux_max + flux_max_err*1.2], $ ystyle=1, psym=4, xtitle='Wavelength ('+angstrom+' )', $ ytitle='Flux (ergs cm!u-2 !n s!u-1 !n '+angstrom+'!u-1!n)', title=strcompress(month[0]+' '+string(day[0])+', '+year[0]), $ charsize=1.5, THICK=8 ;--------------------------------------------------------------------------------------- ; Plot error on each point. ; If the errors are off the scale, multiple by a small number ; just to compare relative sizes of errors on each point. ;-------------------------------------------------------------------------------------- oploterror,lambda_flux, avg_flux[*], flux_err[*], psym = 4, thick = 4 ;plots the error. wavelength = findgen(3000)/1000 + lambda_plot_max ;--------------------------------------------------------------------------------------- ; Plots Guassian of terrestrail green line: ; P[0,i] * exp(-(wavelength-P[1,i])^2.0 / (2.0 * (P[2,i]^2.0)) ; ; and the Doppler shifted green line: ; P[3,i] * exp(-(wavelength-P[4,i])^2.0 / (2.0 * (P[5,i]^2.0))) ; ; added together. It's drawn as a dashed line. If it's ; not matched in flux to the observed emission, something went ; wrong somewhere. ;----------------------------------------------------------------------------------------- oplot, wavelength, ((P[0] * exp(-(wavelength-P[1])^2.0 / (2.0 * (P[2]^2.0))) + $ P[3] * exp(-(wavelength-P[4])^2.0 / (2.0 * (P[2]^2.0))))), $ linestyle=2, COLOR=83,THICK=8 ;--------------------------------------------------------------------------------------- ; Guassian plot of just the terrestral green line, so you can check ; how it compares to the observed line. If this line fits better ; than the blended line, then you don't have green line emission. ;--------------------------------------------------------------------------------------- oplot, wavelength, (P[3]*exp(-(wavelength-P[4])^2.0/(2.0*(P[2]^2.0)))), $ linestyle=4,COLOR=250,THICK=8 ;Terrestrail Gaussian only print, 'parameters:', P print, 'parameter errors:', perror ;--------------------------------------------------------------------------------------- ; Plot the upper and lower bounds to a 1 sigma stand deviation around the mean ; of the scatter of the continuum ;--------------------------------------------------------------------------------------- oplot, [lambda_plot_min, lambda_plot_max], [sigma_flux_continuum_upper, sigma_flux_continuum_upper], THICK = 4, LINESTYLE = 1, COLOR = 10 oplot, [lambda_plot_min, lambda_plot_max], [sigma_flux_continuum_lower, sigma_flux_continuum_lower], THICK = 4, LINESTYLE = 1, COLOR = 10 oplot, [lambda_plot_min, lambda_plot_max], [mean_continuum_flux, mean_continuum_flux], THICK = 4, LINESTYLE = 0, COLOR = 10 ;------------------------------------------------------------------------------------------ ; Integrate over fit and calculate errors ; ; fluxfit - fits the gaussian of the Venusian line flux. ; This will be zero everywhere except near the emission. ; The max value it should have is that of the guessed emission ; strength, which is written in ; parinfo[*].value = [8e-14, 5577.07, 0.0709, 1e-12, 5577.31, 0.0709]; ; ; fluxfittell - fits Gaussian flux of telluric green line flux ; Same hold true here for the zero values everywhere except ; around the line. ; ; "The residual printed is in units of the 1-sigma error bars, so the residual ; should be around 1 or less for most points." ;------------------------------------------------------------------------------------------ print, "" print, "start flux fit" print, "" fluxfit[*] = P[0] * exp(-(lambda - P[1])^2.0 / (2.0 * (P[2]^2.0))) fluxfittell[*] = P[3] * exp(-(lambda - P[4])^2.0 / (2.0 * (P[2]^2.0))) residual[*] = (fluxfit[*] + fluxfittell[*] - (avg_flux[*])) / (flux_err[*]) fiterr=replicate(0.0,N_ELEMENTS(avg_flux)) ;!!! chisquared = 0.0 totresidual = 0.0 for j=0,N_ELEMENTS(avg_flux)-1 do begin ;range of point in lambda to match greenline region totresidual = totresidual + abs(residual[j]) chisquared = chisquared + residual[j]^2.0 endfor redchi = (chisquared)/4.0 ew5577 = 2.50663 * P[0] * P[2] if sqrt(redchi) gt 1.0 then begin ewerr5577 = 2.50663 * perror * P[2] * sqrt(redchi) endif else begin ewerr5577 = 2.50663 * perror * P[2] endelse averesid=totresidual/N_ELEMENTS(avg_flux) ;divided by the number of points you're sampling over (N_elements..) print, 'average residual:', averesid print, 'reduced chi-squared:',redchi print, 'integrated flux:', ew5577 print, 'error in flux:', ewerr5577 DEVICE,/CLOSE ENDFOR FREE_LUN, lun_1 STOP