;PRO Venus_combine_data6 ;------------------------------------------------------------------------------------------ ; Compile: .r Candace_Venus_5577_part16.pro ; execute: Candace_Venus_5577_part16 ;------------------------------------------------------------------------------------------ ; Start with 'function function_name'. The function will use the ; variables followed after the function name defined (lambda, A, F, ; and pder for the gfunct function) ;------------------------------------------------------------------------------------------ ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; **************************** 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 ;------------------------------------------------------------------------------------------ ; function OIflux ; A[0]=CO2 flux contribution, A[1]=co2red, A[2]=h2ored, A[3]=co2green, A[4]=h2ogreen ;------------------------------------------------------------------------------------------ function OIflux, lineratio, A, F; A[0]=CO2 flux contribution, A[1]=co2red, A[2]=h2ored, A[3]=co2green, A[4]=h2ogreen B = A[1] / (A[2]+A[1]) C = (lineratio*A[2] - A[4]) / (A[3] - lineratio*A[1]) F = A[0]/(B*C) return, F END ;PRO Venus_combine_data6 ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************** SECTION 2 - READING IN DATA *********************************************** ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;------------------------------------------------------------------------------------------ ; Read in info for the exposures to be combined. ;------------------------------------------------------------------------------------------ LOADCT,39 SET_PLOT,'PS' !P.CHARTHICK=3 !P.POSITION=[.19,.12,.85,.95] ;set lower left and upper right corners of plotting box !X.CHARSIZE = 1.45 !Y.CHARSIZE = 1.45 !X.THICK=4 !Y.THICK=4 angstrom = '!6!sA!r!u!9 %!6!n' ;define the angstrom character READCOL, 'combine_images_inputs.dat', index, position, month, day, year, airmass, exposure, relative_vel, FORMAT='A,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' 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 ;------------------------------------------------------------------------------------------- ; Loop over how many wavelengths there are to analzye. ;-------------------------------------------------------------------------------------------- FOR i = 0, N_ELEMENTS(labmda_true)-1 DO BEGIN ;------------------------------------------------------------------------------------------- ; 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 IF (STRMATCH(filename[k], '*flux*') EQ 1) AND (STRMATCH(filename[k], '*5577*') EQ 1) THEN BEGIN ; 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 sum_flux_err = sum_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 sum_ray_err = sum_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 sum_norm_err = sum_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) ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************** SECTION 3 - Setting up plotting area**************************************** ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% feature_upper_bound = 0.07 ; Upper wavelength bracket around emission feature; cant be too large or it will grab terrestrail feature_lower_bound = 0.07 ; Lower wavelegnth braket around emission doppler_shift = relative_vel[0] * lambda_true[i] / 3e5 Venus_lambda_guess = lambda_obs[k]+ doppler_shift[k] ; Set Veusian greenline wavelength guess Venus_lambda_upper = lambda_obs[k]+ doppler_shift[k] + feature_upper_bound ; Upper wavelength bracket around emission cant be too large or it will grab terrestrail Venus_lambda_lower = feature_lambda[k]+ doppler_shift[k] - feature_lower_bound ; Lower wavelegnth braket around emission PRINT, Venus_lambda_lower ;========================================================================================= ; FLUX PLOT ;========================================================================================= parinfo = replicate({fixed:0, limited:[1,1], limits:[1e-20,1e6]},5) parinfo[0].limits = [1e-20,1e-8] ; Venusian Flux parinfo[1].limits = [DOUBLE(Venus_lambda_lower),DOUBLE(Venus_lambda_upper)] ; Venusian wavelength parinfo[2].LIMITS = [dispersion[i]-0.001, dispersion[i]+0.001] ; Dispersion of echelle parinfo[3].limits = [1e-20,1e-8] ; Terrestrial Flux parinfo[4].limits = [DOUBLE(lambda_obs[k] - 0.1), DOUBLE(lambda_obs[k] + 0.1)] ; Terrestrial wavelength 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 P=fltarr(5,1);5 is number of parameters, 1 is number of spectra to be analyzed perror=fltarr(5,1) fluxfit = fltarr(N_ELEMENTS(flux),1) ;this used to have 2048 instead of N_ELEMENTS fluxfittell = fltarr(N_ELEMENTS(flux),1) residual = fltarr(N_ELEMENTS(flux),1) totresidual = fltarr(1) chisquared = fltarr(1) redchi = fltarr(1) averesid = fltarr(1) ew5577 = fltarr(1) ew5577_terr = fltarr(1) ewerr5577 = fltarr(1) ewerr5577_terr = fltarr(1) DEVICE, FILENAME = day[0]+month[0]+year[0]+'_'+lambda_true[k]+'_combined_flux.ps',/COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES P[*] = MPFITFUN('gfunct', lambda_flux, avg_flux[*], avg_flux_err[*]/ avg_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_true[i] - plot_lower_bound, lambda_true[i] + plot_upper_bound], /xstyle, $ yrange=[(flux_min - flux_min_err*1.1),(flux_max + flux_max_err*1.1)], /ystyle, psym=4, $ xtitle='Wavelength ('+angstrom+')', ytitle='Flux (ergs cm!u-2 !n s!u-1 !n '+angstrom+'!u-1!n)', $ title = 'Combined', 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, Venus_night_final[*,i], Venus_night_final_err[*,i], psym=4, thick = 4 ;plots the error. Thick for connecting lines ;--------------------------------------------------------------------------------------- ; 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. ;----------------------------------------------------------------------------------------- wavelength = findgen(3000)/1000 + lambda_true[i] - 1.5 ; Where the x range will start 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=82,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 oplot, wavelength, (P[0]*exp(-(wavelength-P[1])^2.0/(2.0*(P[2]^2.0)))), $ linestyle=4,COLOR=150,THICK=8 ;Terrestrail Gaussian only print, 'parameters:', P print, 'parameter errors:', perror ;------------------------------------------------------------------------------------------ ; 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_flux - P[1])^2.0 / (2.0 * (P[2]^2.0))) fluxfittell[*] = P[3] * exp(-(lambda_flux - P[4])^2.0 / (2.0 * (P[2]^2.0))) residual[*] = (fluxfit[*,i] + fluxfittell[*] - (avg_flux[*])) / (flux_err[*,i]) fiterr=replicate(0.0,N_ELEMENTS(avg_flux)) ;!!! chisquared[i] = 0.0 totresidual[i] = 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[i])/4.0 ew5577 = 2.50663 * P[0] * P[2] if sqrt(redchi) gt 1.0 then begin ewerr5577 = 2.50663 * perror[0] * P[2] * sqrt(redchi) endif else begin ewerr5577 = 2.50663 * perror[0] * P[2] endelse averesid=totresidual/N_ELEMENTS(Venus_night_final) ;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 print, '********************PLOT COMPLETE***************' ENDFOR STOP END