;.r Venus_combine_data4.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 ;------------------------------------------------------------------------------------------ ; 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 ;------------------------------------------------------------------------------------------ ; 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' READCOL, 'combine_images_inputs.dat', index, position, month, day, year, airmass, exposure, relative_vel, FORMAT='A,A,A,A,I,F,F,F' ;--------------------------------- ;Read in data from data files ;--------------------------------- ; openr, lun, venus_night_name[k], /GET_LUN ;**** UPDATE VARIABLE **** ; Venus = fltarr(3,55) ; readf, lun, Venus ; FREE_LUN, lun ; lambda = double(reform(Venus[0,*])) ; wavelength value from Venus array (flux and wavelength) ; Venus_night = double(reform(Venus[1,*])+1.0) ; Flux values from Venus array ; Venus_night_err = double(reform(Venus[2,*])) ; Flux values from Venus array READCOL,'output_0015_OI_5577_flux.dat', lambda1_flux, flux1_flux, flux1_flux_err, FORMAT='F,F,F' READCOL,'output_0015_OI_5577_rayleigh.dat', lambda1_rayleigh, flux1_rayleigh, flux1_rayleigh_err, FORMAT='F,F,F' READCOL,'output_0015_OI_5577_normalized.dat', lambda1_norm, flux1_norm, flux1_norm_err, FORMAT='F,F,F' READCOL, 'output_0016_OI_5577_flux.dat', lambda2_flux, flux2_flux, flux2_flux_err, FORMAT='F,F,F' READCOL,'output_0016_OI_5577_rayleigh.dat', lambda2_rayleigh, flux2_rayleigh, flux2_rayleigh_err, FORMAT='F,F,F' READCOL,'output_0016_OI_5577_normalized.dat', lambda2_norm, flux2_norm, flux2_norm_err, FORMAT='F,F,F' READCOL, 'output_0017_OI_5577_flux.dat', lambda3_flux, flux3_flux, flux3_flux_err, FORMAT='F,F,F' READCOL,'output_0017_OI_5577_rayleigh.dat', lambda3_rayleigh, flux3_rayleigh, flux3_rayleigh_err, FORMAT='F,F,F' READCOL,'output_0017_OI_5577_normalized.dat', lambda3_norm, flux3_norm, flux3_norm_err, FORMAT='F,F,F' READCOL, 'output_0018_OI_5577_flux.dat', lambda4_flux, flux4_flux, flux4_flux_err, FORMAT='F,F,F' READCOL,'output_0018_OI_5577_rayleigh.dat', lambda4_rayleigh, flux4_rayleigh, flux4_rayleigh_err, FORMAT='F,F,F' READCOL,'output_0018_OI_5577_normalized.dat', lambda4_norm, flux4_norm, flux4_norm_err, FORMAT='F,F,F' READCOL, 'output_0019_OI_5577_flux.dat', lambda5_flux, flux5_flux, flux5_flux_err, FORMAT='F,F,F' READCOL,'output_0019_OI_5577_rayleigh.dat', lambda5_rayleigh, flux5_rayleigh, flux5_rayleigh_err, FORMAT='F,F,F' READCOL,'output_0019_OI_5577_normalized.dat', lambda5_norm, flux5_norm, flux5_norm_err, FORMAT='F,F,F' READCOL, 'output_0020_OI_5577_flux.dat', lambda6_flux, flux6_flux, flux6_flux_err, FORMAT='F,F,F' READCOL,'output_0020_OI_5577_rayleigh.dat', lambda6_rayleigh, flux6_rayleigh, flux6_rayleigh_err, FORMAT='F,F,F' READCOL,'output_0020_OI_5577_normalized.dat', lambda6_norm, flux6_norm, flux6_norm_err, FORMAT='F,F,F' READCOL, 'output_0021_OI_5577_flux.dat', lambda7_flux, flux7_flux, flux7_flux_err, FORMAT='F,F,F' READCOL,'output_0021_OI_5577_rayleigh.dat', lambda7_rayleigh, flux7_rayleigh, flux7_rayleigh_err, FORMAT='F,F,F' READCOL,'output_0021_OI_5577_normalized.dat', lambda7_norm, flux7_norm, flux7_norm_err, FORMAT='F,F,F' READCOL, 'output_0000_OI_5577_flux.dat', lambda8_flux, flux8_flux, flux8_flux_err, FORMAT='F,F,F' READCOL,'output_0000_OI_5577_rayleigh.dat', lambda8_rayleigh, flux8_rayleigh, flux8_rayleigh_err, FORMAT='F,F,F' READCOL,'output_0000_OI_5577_normalized.dat', lambda8_norm, flux8_norm, flux8_norm_err, FORMAT='F,F,F' ;----------------SKY WITH VENUS DAY SUBTRACTION--------------------- READCOL, 'output_NNNN_OI_5577_flux.dat', lambda0_flux, flux0_flux, flux0_flux_err, FORMAT='F,F,F' READCOL,'output_NNNN_OI_5577_rayleigh.dat', lambda0_rayleigh, flux0_rayleigh, flux0_rayleigh_err, FORMAT='F,F,F' READCOL,'output_NNNN_OI_5577_normalized.dat', lambda0_norm, flux0_norm, flux0_norm_err, FORMAT='F,F,F' lambda = lambda1_norm ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************** SECTION 3 - Combine Data ************************************************** ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;----------------------------------------------------------------------------------------------------------- ; Combine a straight average of the data, no weighting, as they should ; all have the same exposure time. ;----------------------------------------------------------------------------------------------------------- equator_flux = (flux1_flux + flux2_flux + flux3_flux + flux4_flux + flux5_flux + flux6_flux + flux7_flux)/7.0 equator_rayleigh = (flux1_rayleigh + flux2_rayleigh + flux3_rayleigh + flux4_rayleigh + flux5_flux + flux6_rayleigh + flux7_rayleigh)/7.0 equator_norm = (flux8_norm - flux0_norm) ;equator_norm = (flux1_norm/MAX(flux1_norm) + flux2_norm/MAX(flux2_norm))/2.0 ;----------------------------------------------------------------------------------------------------------- ; Calculate error for combined spectra, using the propagation: ; error = sqrt[ (flux1_err)^2 + (flux2_err)^2 + ...(fluxn_err)^2 ] ;;----------------------------------------------------------------------------------------------------------- equator_flux_err = sqrt((flux1_flux_err)^2.0 + (flux2_flux_err)^2.0) equator_rayleigh_err = sqrt((flux1_rayleigh_err)^2.0 + (flux2_rayleigh_err)^2.0) equator_norm_err = sqrt((flux8_norm_err)^2.0 + (flux0_norm_err)^2.0) ;----------------------------------------------------------------------------------------------------------- ; Calculate the Doppler shift and set a range for the Venusian green line ; as well an inital guess for the position ; Calculate max and min of whole data set ;----------------------------------------------------------------------------------------------------------- doppler_shift = relative_vel[0] * 5577.34 / 3e5 Venus_lambda_guess = doppler_shift + 5577.3 ; Set Veusian greenline wavelength guess Venus_lambda_upper = 5577.3 + doppler_shift + 0.08 ; Upper wavelength bracket around emission Venus_lambda_lower = 5577.3 + doppler_shift - 0.08 ; Lower wavelegnth braket around emission ;----------------------------------------------------------------------------------------------------------- ; Calculate max and min flux of all files read in. This value will be ; used to define the range of the plots ;----------------------------------------------------------------------------------------------------------- equator_flux_max = MAX(equator_flux, index) equator_flux_max_err = equator_flux_err(index) equator_rayleigh_max = MAX(equator_rayleigh, index) equator_rayleigh_max_err = equator_rayleigh_err(index) equator_norm_max = MAX(equator_norm, index) equator_norm_max_err = equator_norm_err(index) equator_flux_min = MIN(equator_flux, index) equator_flux_min_err = equator_flux_err(index) equator_rayleigh_min = MIN(equator_rayleigh, index) equator_rayleigh_min_err = equator_rayleigh_err(index) equator_norm_min = MIN(equator_norm, index) equator_norm_min_err = equator_norm_err(index) ;----------------------------------------------------------------------------------------------------------- ; Find the max value for the terrestrial and the venusian peaks in ; flux, rayleigh, and normalized units and their wavelength ;----------------------------------------------------------------------------------------------------------- linerange = WHERE((lambda GE 5577.23)AND(lambda LE 5577.4)) ;Specify wavelength range we will analyze equator_terr_flux_max = MAX(equator_flux(linerange),index) equator_terr_flux_max_err = equator_flux_err(linerange(index)) equator_terr_lambda_max = lambda1_flux(linerange(index)) equator_terr_rayleigh_max = MAX(equator_rayleigh(linerange),index) equator_terr_rayleigh_err = equator_rayleigh_err(linerange(index)) equator_terr_norm_max = MAX(equator_norm(linerange),index) equator_terr_norm_max_err = equator_norm_err(linerange(index)) linerange = WHERE((lambda GE Venus_lambda_lower) AND (lambda LE Venus_lambda_upper)) ;Specify wavelength range we will analyze equator_venus_flux_max = MAX(equator_flux(linerange),index) equator_venus_flux_max_err = equator_flux_err(linerange(index)) equator_venus_lambda_max = lambda1_flux(linerange(index)) equator_venus_rayleigh_max = MAX(equator_rayleigh(linerange),index) equator_venus_rayleigh_err = equator_rayleigh_err(linerange(index)) equator_venus_norm_max = MAX(equator_norm(linerange),index) equator_venus_norm_max_err = equator_norm_err(linerange(index)) IF Venus_lambda_guess LT 5577.3 THEN BEGIN ; Set lower wavelength bounds to be Venus if lambda_lower = Venus_lambda_guess ; Venus is blue shifted of terrestrial lambda_upper = equator_terr_lambda_max ENDIF IF Venus_lambda_guess GT 5577.3 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_lambda_guess ENDIF lambda_continuum_lower = lambda_lower - 0.15 lambda_continuum_upper = lambda_upper + 0.15 linerange_continuum = WHERE ( (lambda LE (lambda_continuum_lower)) OR (lambda GE (lambda_continuum_upper + 0.15))) mean_continuum_equator_flux = MEAN((equator_flux(linerange_continuum))) ;calculate mean of continuum mean_continuum_equator_ray = MEAN((equator_rayleigh(linerange_continuum))) mean_continuum_equator_norm = MEAN((equator_norm(linerange_continuum))) sigma_continuum_equator_flux = stddev(equator_flux(linerange_continuum)) ; calculate standard deviation of the sigma_continuum_equator_ray = stddev(equator_rayleigh(linerange_continuum)) ; mean of the continuum sigma_continuum_equator_norm = stddev(equator_norm(linerange_continuum)) ;-------------------------------------------------------------------------------------------- ; Calculate how many standard deviations the greenline is above the average of the continuum ;-------------------------------------------------------------------------------------------- Venus_greenline_confidence_equator_flux = (equator_venus_flux_max - mean_continuum_equator_flux) / sigma_continuum_equator_flux Venus_greenline_confidence_equator_ray = (equator_venus_rayleigh_max - mean_continuum_equator_ray) / sigma_continuum_equator_ray Venus_greenline_confidence_equator_norm = (equator_venus_norm_max - mean_continuum_equator_norm) / sigma_continuum_equator_norm ;-------------------------------------------------------------------------------------------- ; Calculate the upper and lower bounds of the standard deviation above the mean so as to ; plot the deviations on the plot. ;-------------------------------------------------------------------------------------------- sigma_equator_flux_continuum_upper = mean_continuum_equator_flux + sigma_continuum_equator_flux sigma_equator_flux_continuum_lower = mean_continuum_equator_flux - sigma_continuum_equator_flux sigma_equator_ray_continuum_upper = mean_continuum_equator_ray + sigma_continuum_equator_ray sigma_equator_ray_continuum_lower = mean_continuum_equator_ray - sigma_continuum_equator_ray sigma_equator_norm_continuum_upper = mean_continuum_equator_norm + sigma_continuum_equator_norm sigma_equator_norm_continuum_lower = mean_continuum_equator_norm - sigma_continuum_equator_norm ;------------------------------------------------------------------------------------------ ; Print information to file ;------------------------------------------------------------------------------------------ printf, lun_1, 'Venus_combined', 'Equator', 'Flux ', lambda[0], lambda_continuum_lower, lambda_continuum_upper, $ lambda[N_elements(lambda)-1], mean_continuum_equator_flux, sigma_continuum_equator_flux, $ equator_terr_flux_max, equator_venus_flux_max, Venus_greenline_confidence_equator_flux, $ FORMAT='(A,13X, A, 5X, 15A, 5X, F, 4X, F, X, e5.2, e5.2, 10X, e5.2, X, e5.2, X, f5.2)' printf, lun_1, 'Venus_combined', 'Equator', 'Ray ', lambda[0], lambda_continuum_lower, lambda_continuum_upper, $ lambda[N_elements(lambda)-1], mean_continuum_equator_ray, sigma_continuum_equator_ray, $ equator_terr_rayleigh_max, equator_venus_rayleigh_max, Venus_greenline_confidence_equator_ray, $ FORMAT='(A,13X, A, 5X, 15A, 5X, F,4X,F, X, e5.2, e5.2, 10X, e5.2, X, e5.2, X, f5.2)' printf, lun_1, 'Venus_combined', 'Equator', 'Norm ', lambda[0], lambda_continuum_lower, lambda_continuum_upper, $ lambda[N_elements(lambda)-1], mean_continuum_equator_norm, sigma_continuum_equator_norm, $ equator_terr_norm_max, equator_venus_norm_max, Venus_greenline_confidence_equator_norm, $ FORMAT='(A,13X, A, 5X, 15A, 5X, F, 4X, F, X, e5.2, e5.2, 10X, e5.2, X, e5.2, X, f5.2)' ;------------------------------------------------------------------------------------------ ;---------------------------------------------------------------------------------------------------------- ; Set a wavelength for all the points to be plotte again ; Calculate the average airmass for combined equator and south limb exposure. ;----------------------------------------------------------------------------------------------------------- lambda = lambda1_flux airmass_equator = (airmass[0]+airmass[1])/2.0 ;average airmass of combined exposures ;------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------ ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************************** FLUX UNITS EQUATOR ***************************************** ; ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PRINT,'' PRINT, '' ; tmpread='' ; READ,'Enter 1 for Flux units, 2 for Rayleighs, or 3 for relative scale: ',tmpread ; If(tmpread EQ '1')Then Begin PRINT,'' PRINT, '' print, 'The following calculations will be done in cgs flux units for the equator' PRINT,'' PRINT, '' ;------------------------------------------------------------------------------------------------ ; Find error due to scatter in the point around 0. Do not consider the terrestrial or Venusian ; green line peaks, only the noise. Data is written to the same file as the indiviatual spectra ;------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------ ; 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-20,1e6]},5) parinfo[0].limits = [1e-20,1e-08] ; Venusian Flux parinfo[1].limits = [DOUBLE(Venus_lambda_lower[0]), DOUBLE(Venus_lambda_upper[0])] parinfo[2].limits = [0.079,0.079] parinfo[3].limits = [1e-20,1e-8] ; Terrestrial Flux parinfo[4].limits = [5577.23, 5577.4] start = [DOUBLE(equator_venus_flux_max), DOUBLE(Venus_lambda_guess[0]), 0.079D, DOUBLE(equator_terr_flux_max), 5577.3D] ;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(lambda),1) ;this used to have 2048 instead of N_ELEMENTS fluxfittell = fltarr(N_ELEMENTS(lambda),1) residual = fltarr(N_ELEMENTS(equator_flux),1) totresidual = fltarr(1) chisquared = fltarr(1) redchi = fltarr(1) averesid = fltarr(1) ew5577 = fltarr(1) ewerr5577 = fltarr(1) ;------------------------------------------------------------------------------------------------------------ ; ------------------------------------ PLOTTING IN FLUX UNITS EQUATOR---------------------------------------- ;------------------------------------------------------------------------------------------------------------ angstrom = '!6!sA!r!u!9 %!6!n' ;define the angstrom character print, 'begin loop' 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' DEVICE,FILENAME=strcompress(string(day[0])+month[0]+'2022_SKY_'+position[0]+'_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, equator_flux[*], equator_flux_err[*]/equator_flux_err[*] * 1e-18, start, $ PARINFO=parinfo, NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1) plot, lambda1_flux, equator_flux[*], xrange=[5576.5, 5578.35], /xstyle, $ yrange=[equator_flux_min - equator_flux_min_err*1.2, equator_flux_max + equator_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])+', 2012 - '+position[0]+' Limb'), $ 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,lambda1_flux, equator_flux[*], equator_flux_err[*], psym = 4, thick = 4 ;plots the error. wavelength = findgen(3000)/1000 + 5576 ;--------------------------------------------------------------------------------------- ; 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, [5575, 5579], [sigma_equator_flux_continuum_upper, sigma_equator_flux_continuum_upper], THICK = 4, LINESTYLE = 1, COLOR = 10 oplot, [5575, 5579], [sigma_equator_flux_continuum_lower, sigma_equator_flux_continuum_lower], THICK = 4, LINESTYLE = 1, COLOR = 10 oplot, [5575, 5579], [mean_continuum_equator_flux, mean_continuum_equator_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[*] - (equator_flux[*])) / (equator_flux_err[*]) fiterr=replicate(0.0,N_ELEMENTS(flux1_flux)) ;!!! chisquared = 0.0 totresidual = 0.0 for j=0,N_ELEMENTS(equator_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(flux1_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 ; endfor ;------------------------------------------------------------------------------------------ ; Print data to file ;------------------------------------------------------------------------------------------ openw, lun, 'output_Equator_combined_flux.dat', /GET_LUN ;**** UPDATE VARIABLE **** printf, lun, '#---------------------------------------------------' printf, lun, '# Wavelength Flux Error' printf, lun, '# (A) (cgs) ' printf, lun, '#---------------------------------------------------' for i=0, N_ELEMENTS(lambda)-1 do begin printf, lun, lambda[i], equator_flux[i], equator_flux_err[i] endfor FREE_LUN, lun ;------------------------------------------------------------------------------------------ DEVICE,/CLOSE ;------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------ ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ************************************** RAYLEIGH UNITS EQUATOR**************************************** ; ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PRINT,'' PRINT, '' print, 'The following calculations will be done in Rayleigh units for the equator' PRINT,'' PRINT, '' terr_max = -1 FOR n=16, 19 DO BEGIN if equator_rayleigh[n] GT terr_max THEN terr_max = equator_rayleigh[n] ENDFOR venus_max = -1 FOR n=19, 23 DO BEGIN if equator_rayleigh[n] GT venus_max THEN venus_max = equator_rayleigh[n] ENDFOR ;------------------------------------------------------------------------------------------ ; Set limits for your parameters and give an initial guess parinfo = replicate({fixed:0, limited:[1,1], limits:[1e-20,1e6]},5) parinfo[0].limits = [0, equator_venus_rayleigh_max*2.0] ; Venusian Relative parinfo[1].limits = [DOUBLE(Venus_lambda_lower[0]), DOUBLE(Venus_lambda_upper[0])] parinfo[2].limits = [0.078,0.08] parinfo[3].limits = [0., equator_terr_rayleigh_max*2.0] ; Terrestrial Relative parinfo[4].limits = [5577.23, 5577.4] start = [DOUBLE(equator_venus_rayleigh_max), DOUBLE(Venus_lambda_guess[0]), 0.079D, $ DOUBLE(equator_terr_rayleigh_max), 5577.3D] 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(equator_rayleigh),1) ;this used to have 2048 instead of N_ELEMENTS fluxfittell = fltarr(N_ELEMENTS(equator_rayleigh),1) residual = fltarr(N_ELEMENTS(equator_rayleigh),1) totresidual = fltarr(1) chisquared = fltarr(1) redchi = fltarr(1) averesid = fltarr(1) ew5577 = fltarr(1) ewerr5577 = fltarr(1) ;------------------------------------------------------------------------------------------------------------ ; -------------------------------- PLOTTING IN RAYLEIGH UNITS EQUATOR---------------------------------------- ;------------------------------------------------------------------------------------------------------------ angstrom = '!6!sA!r!u!9 %!6!n' ;define the angstrom character print, 'begin loop' 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' DEVICE,FILENAME=strcompress(string(day[0])+month[0]+'2022_SKY_'+position[0]+'_combined_rayleigh_FINAL_errors.ps',/remove_all),$ /COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES ;------------------------------------------------------------------------------------------------------------- ;--------------------------------------- BEGIN GAUSSIAN FIT FOR PLOT ----------------------------------------- ;------------------------------------------------------------------------------------------------------------- ; for i=0, 0 do begin ;0-9 is 10, the number of spectra that are considered. For one spectra, loop goes from 0 to 0 ; P[*,i] = MPFITFUN('gfunct', lambda, Hartley_Oct1[*,i], Hartley_Oct1_err[*,i], start, PARINFO=parinfo, $ ; NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1) P[*] = MPFITFUN('gfunct', lambda, equator_rayleigh[*], equator_rayleigh_err[*]/equator_rayleigh_err[*] * 1e-18, start, PARINFO=parinfo, $ NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1) plot, lambda, equator_rayleigh[*], xrange=[5576.5, 5578.35], /xstyle, $ yrange=[equator_rayleigh_min - equator_rayleigh_min_err*1.2, equator_rayleigh_max + equator_rayleigh_max_err*1.2], $ ystyle=9, psym=4,xtitle='Wavelength ('+angstrom+' )', $ ytitle='Venusian Intesnsity (R/'+angstrom+')', title=strcompress(month[0]+' '+string(day[0])+', 2012 - '+position[0]+' Limb'), $ charsize=1.5, THICK=8 ;----------------------------------------------------------------------------------------------------- ; Set up different range on right axis for the Terrestrial intensity. Divide by airmass to account for ; Csc z, as Rayleighs are defined at an airmass of 1, we need to scale to 1. ; To get the axis to go down and not up, you need to do a XYOutS command and change the orientation by ; 90 degrees the opposite way. ;----------------------------------------------------------------------------------------------------- axis,yaxis=1,yrange=[(equator_rayleigh_min - equator_rayleigh_min_err*1.2)/airmass_equator, (equator_rayleigh_max + equator_rayleigh_max_err*1.2)/airmass_equator], $ ystyle=1, charsize=1.5 xloc = !X.Window[1] + 0.102 yloc = (0.9 - 0.15) / 2.0 + 0.15 XYOutS, xloc, yloc, 'Terrestrial Intensity (R/' +angstrom+')', Alignment=0.5, Orientation=-90.0, charsize=2.0, /Normal ;--------------------------------------------------------------------------------------- ; 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, equator_rayleigh[*], equator_rayleigh_err[*], psym = 4, thick = 4 ;plots the error. ;--------------------------------------------------------------------------------------- ; Plots Guassian of the Doppler shifted Venusian green line: ; P[0,i] * exp(-(wavelength-P[1,i])^2.0 / (2.0 * (P[2,i]^2.0)) ; ; and the terrestrial 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 + 5576 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 ; Venusian plus terrestrial ;--------------------------------------------------------------------------------------- ; 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=3, COLOR=145,THICK=8 ;Venusian Gaussian only ; oplot, wavelength, P[3]*exp(-(wavelength-P[4])^2.0/(2.0*(P[2]^2.0))) + P[0] * exp(-(wavelength-P[1])^2.0 / (2.0 * (P[2]^2.0))), linestyle=2, COLOR=83,THICK=8 ;--------------------------------------------------------------------------------------- ; Plot the upper and lower bounds to a 1 sigma stand deviation around the mean ; of the scatter of the continuum ;--------------------------------------------------------------------------------------- oplot, [5575, 5579], [sigma_equator_ray_continuum_upper, sigma_equator_ray_continuum_upper], THICK = 4, LINESTYLE = 1, COLOR = 10 oplot, [5575, 5579], [sigma_equator_ray_continuum_lower, sigma_equator_ray_continuum_lower], THICK = 4, LINESTYLE = 1, COLOR = 10 oplot, [5575, 5579], [mean_continuum_equator_ray, mean_continuum_equator_ray], THICK = 4, LINESTYLE = 0, COLOR = 10 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 - 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[*] - (equator_flux[*])) / (equator_flux_err[*]) fiterr=replicate(0.0,N_ELEMENTS(equator_rayleigh)) ;!!! chisquared = 0.0 totresidual = 0.0 for j=0,N_ELEMENTS(equator_rayleigh)-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(equator_rayleigh) ;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 ; endfor ;------------------------------------------------------------------------------------------ ; Print data to file ;------------------------------------------------------------------------------------------ openw, lun, 'output_Equator_combined_rayleigh.dat', /GET_LUN ;**** UPDATE VARIABLE **** printf, lun, '#---------------------------------------------------' printf, lun, '# Wavelength Rayleigh Error' printf, lun, '# (A) ' printf, lun, '#---------------------------------------------------' for i=0, N_ELEMENTS(lambda)-1 do begin printf, lun, lambda1_rayleigh[i], equator_rayleigh[i], equator_rayleigh_err[i] endfor FREE_LUN, lun ;------------------------------------------------------------------------------------------ DEVICE,/CLOSE ;------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------ ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ************************************** NORMALIZED UNITS EQUATOR**************************************** ; ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PRINT,'' PRINT, '' ; tmpread='' ; READ,'Enter 1 for Flux units, 2 for Rayleighs, or 3 for relative scale: ',tmpread ; If(tmpread EQ '1')Then Begin PRINT,'' PRINT, '' print, 'The following calculations will be done in normalized units for the equator' PRINT,'' PRINT, '' ; terr_max = -1 ; FOR n=10, 15 DO BEGIN ; if equator_norm[n] GT terr_max THEN terr_max = equator_norm[n] ; ENDFOR ; venus_max = -1 ; FOR n=15, 20 DO BEGIN ; if equator_norm[n] GT venus_max THEN venus_max = equator_norm[n] ; ENDFOR ;------------------------------------------------------------------------------------------ ; Set limits for your parameters and give an initial guess parinfo = replicate({fixed:0, limited:[1,1], limits:[1e-20,1e6]},5) parinfo[0].limits = [0, 5] ; Venusian Relative parinfo[1].limits = [DOUBLE(Venus_lambda_lower[0]), DOUBLE(Venus_lambda_upper[0])] parinfo[2].limits = [0.078,0.08] parinfo[3].limits = [0.0, 1.2] ; Terrestrial Relative parinfo[4].limits = [5577.23, 5577.4] start = [DOUBLE(equator_venus_norm_max), DOUBLE(Venus_lambda_guess[0]), 0.079D, DOUBLE(equator_terr_norm_max), 5577.29D] ;INITIAL GUESSES, Relat 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(equator_norm),1) ;this used to have 2048 instead of N_ELEMENTS fluxfittell = fltarr(N_ELEMENTS(equator_norm),1) residual = fltarr(N_ELEMENTS(equator_norm),1) totresidual = fltarr(1) chisquared = fltarr(1) redchi = fltarr(1) averesid = fltarr(1) ew5577 = fltarr(1) ewerr5577 = fltarr(1) ;------------------------------------------------------------------------------------------------------------ ; --------------------------------- PLOTTING IN NORMALIZED UNITS EQUATOR------------------------------------- ;------------------------------------------------------------------------------------------------------------ angstrom = '!6!sA!r!u!9 %!6!n' ;define the angstrom character print, 'begin loop' 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' DEVICE,FILENAME=strcompress(string(day[0])+month[0]+'2022_SKY_'+position[0]+'_combined_norm_FINAL_errors.ps',/remove_all),$ /COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES ;------------------------------------------------------------------------------------------------------------- ;--------------------------------------- BEGIN GAUSSIAN FIT FOR PLOT ----------------------------------------- ;------------------------------------------------------------------------------------------------------------- ; for i=0, 0 do begin ;0-9 is 10, the number of spectra that are considered. For one spectra, loop goes from 0 to 0 ; P[*,i] = MPFITFUN('gfunct', lambda, Hartley_Oct1[*,i], Hartley_Oct1_err[*,i], start, PARINFO=parinfo, $ ; NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1) P[*] = MPFITFUN('gfunct', lambda, equator_norm[*], equator_norm_err[*]/equator_norm_err[*] * 1e-18, start, PARINFO=parinfo, $ NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1) plot, lambda, equator_norm[*], xrange=[5576.5, 5578.35], /xstyle, $ yrange=[equator_norm_min * 1.1, equator_norm_max * 1.1], $ ystyle=1, psym=10,xtitle='Wavelength ('+angstrom+' )', $ ytitle='Normalized Units', title=strcompress(month[0]+' '+string(day[0])+', 2012 - '+position[0]+' Limb'), $ charsize=1.5, THICK=8 oplot, [5576.57, 5576.74], [0.8, 0.8], linestyle = 2, color = 82, thick = 8 oplot, [5576.57, 5576.74],[0.65, 0.65], linestyle = 4, color = 250, thick = 8 xyouts, 5576.77, 0.78, 'Terrestrial', color = 82,charsize = 1.7 xyouts, 5576.77, 0.73, '+ Venus', color = 82,charsize = 1.7 xyouts, 5576.77, 0.63, 'Terrestrial', color=250, charsize = 1.7 ;--------------------------------------------------------------------------------------- ; 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, equator_norm[*], equator_norm_err[*] ;plots the error. wavelength = findgen(3000)/1000 + 5576 ;--------------------------------------------------------------------------------------- ; 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, [5575, 5579], [sigma_equator_norm_continuum_upper, sigma_equator_norm_continuum_upper], THICK = 4, LINESTYLE = 1, COLOR = 10 oplot, [5575, 5579], [sigma_equator_norm_continuum_lower, sigma_equator_norm_continuum_lower], THICK = 4, LINESTYLE = 1, COLOR = 10 oplot, [5575, 5579], [mean_continuum_equator_norm, mean_continuum_equator_norm], 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[*] - (equator_flux[*])) / (equator_flux_err[*]) fiterr=replicate(0.0,N_ELEMENTS(flux1_norm)) ;!!! chisquared = 0.0 totresidual = 0.0 for j=0,N_ELEMENTS(equator_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(flux1_norm) ;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 ; endfor ;------------------------------------------------------------------------------------------ ; Print data to file ;------------------------------------------------------------------------------------------ openw, lun, 'output_Equator_combined_norm.dat', /GET_LUN ;**** UPDATE VARIABLE **** printf, lun, '#---------------------------------------------------' printf, lun, '# Wavelength Normalized Error' printf, lun, '# (A) Flux' printf, lun, '#---------------------------------------------------' for i=0, N_ELEMENTS(lambda)-1 do begin printf, lun, lambda1_norm[i], equator_norm[i], equator_norm_err[i] endfor FREE_LUN, lun ;------------------------------------------------------------------------------------------ DEVICE,/CLOSE FREE_LUN, lun_1 DEVICE,/CLOSE STOP end