;.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 ;------------------------------------------------------------------------------------------ ; 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 ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************** SECTION 2 - READING IN DATA *********************************************** ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% READCOL, 'combine_images_inputs.dat', index, position, month, day, airmass, FORMAT='A,A,A,I,F' ; FOR k=0, N_ELEMENTS(index)-1 DO BEGIN ; READCOL,'output_'+index[k]+'_flux.dat', lambda[k]_flux, flux[k]_flux, FORMATE='F,F' ; READCOL,'output_'+index[k]+'_rayleigh.dat', lambda[k]_rayleigh, flux[k]_rayleigh, FORMATE='F,F' ; READCOL,'output_'+index[k]+'_normalized.dat', lambda[k]_norm, flux[k]_norm, FORMATE='F,F' ;ENDFOR ;stop ;------------------- ;Equator exposures ;------------------- READCOL,'output_0004_flux.dat', lambda1_flux, flux1_flux, FORMAT='F,F' READCOL,'output_0004_rayleigh.dat', lambda1_rayleigh, flux1_rayleigh, FORMAT='F,F' READCOL,'output_0004_normalized.dat', lambda1_norm, flux1_norm, FORMAT='F,F' READCOL,'output_0005_flux.dat', lambda2_flux, flux2_flux, FORMAT='F,F' READCOL,'output_0005_rayleigh.dat', lambda2_rayleigh, flux2_rayleigh, FORMAT='F,F' READCOL,'output_0005_normalized.dat', lambda2_norm, flux2_norm, FORMAT='F,F' ;------------------- ;South exposures ;------------------- READCOL,'output_0006_flux.dat', lambda3_flux, flux3_flux, FORMAT='F,F' READCOL,'output_0006_rayleigh.dat', lambda3_rayleigh, flux3_rayleigh, FORMAT='F,F' READCOL,'output_0006_normalized.dat', lambda3_norm, flux3_norm, FORMAT='F,F' READCOL,'output_0007_flux.dat', lambda4_flux, flux4_flux, FORMAT='F,F' READCOL,'output_0007_rayleigh.dat', lambda4_rayleigh, flux4_rayleigh, FORMAT='F,F' READCOL,'output_0007_normalized.dat', lambda4_norm, flux4_norm, FORMAT='F,F' ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************** SECTION 3 - Combine Data ************************************************** ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;----------------------------------------------------------------------------------------------------------- ; Combine a straight average of the data, no weighting ;----------------------------------------------------------------------------------------------------------- equator_flux = (flux1_flux + flux2_flux)/2.0 equator_rayleigh = (flux1_rayleigh + flux2_rayleigh)/2.0 equator_norm = (flux1_norm + flux2_norm)/2.0 south_flux = (flux3_flux + flux4_flux)/2.0 south_rayleigh = (flux3_rayleigh + flux4_rayleigh)/2.0 south_norm = (flux3_norm + flux4_norm)/2.0 ;----------------------------------------------------------------------------------------------------------- ; Find the max value for the terrestrial and the venusian peaks in ; flux, rayleigh, and normalized units ;----------------------------------------------------------------------------------------------------------- equ_terr_flux_max = -10 equ_terr_rayleigh_max = -10 equ_terr_norm_max = -10 south_terr_flux_max = -10 south_terr_rayleigh_max = -10 south_terr_norm_max = -10 FOR n=13, 19 DO BEGIN if equator_flux[n] GT equ_terr_flux_max THEN equ_terr_flux_max = equator_flux[n] if equator_rayleigh[n] GT equ_terr_rayleigh_max THEN equ_terr_rayleigh_max = equator_rayleigh[n] if equator_norm[n] GT equ_terr_norm_max THEN equ_terr_norm_max = equator_norm[n] if south_flux[n] GT south_terr_flux_max THEN south_terr_flux_max = south_flux[n] if south_rayleigh[n] GT south_terr_rayleigh_max THEN south_terr_rayleigh_max = south_rayleigh[n] if south_norm[n] GT south_terr_norm_max THEN south_terr_norm_max = south_norm[n] ENDFOR equ_venus_flux_max = -10 equ_venus_rayleigh_max = -10 equ_venus_norm_max = -10 south_venus_flux_max = -10 south_venus_rayleigh_max = -10 south_venus_norm_max = -10 FOR n=21, 26 DO BEGIN if equator_flux[n] GT equ_venus_flux_max THEN equ_venus_flux_max = equator_flux[n] if equator_rayleigh[n] GT equ_venus_rayleigh_max THEN equ_venus_rayleigh_max = equator_rayleigh[n] if equator_norm[n] GT equ_venus_norm_max THEN equ_venus_norm_max = equator_norm[n] if south_flux[n] GT south_venus_flux_max THEN south_venus_flux_max = south_flux[n] if south_rayleigh[n] GT south_venus_rayleigh_max THEN south_venus_rayleigh_max = south_rayleigh[n] if south_norm[n] GT south_venus_norm_max THEN south_venus_norm_max = south_norm[n] ENDFOR ;----------------------------------------------------------------------------------------------------------- ; Set error levels (for now, this is just a fraction of the max value ; so there's something to plot) ;----------------------------------------------------------------------------------------------------------- equator_flux_max = MAX(equator_flux) equator_rayleigh_max = MAX(equator_rayleigh) south_flux_max = MAX(south_flux) south_rayleigh_max = MAX(south_rayleigh) equator_flux_err = equator_flux_max*0.1 equator_rayleigh_err = equator_rayleigh_max*0.1 equator_norm_err = 0.1 south_flux_err = south_flux_max*0.1 south_rayleigh_err = south_rayleigh_max*0.1 south_norm_err = 0.1 lambda = lambda1_flux airmass_equator = (airmass[0]+airmass[1])/2.0 ;average airmass of combined exposures airmass_south = (airmass[2]+airmass[3])/2.0 ;will need for accounting for terrestrial Rayleigh scaled to 1 airmass ;------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------ ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************************** FLUX UNITS ************************************************* ; ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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' PRINT,'' PRINT, '' ; Hartley_Oct1 = [equator_flux] ;------------------------------------------------------------------------------------------ ; 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-10] ; Venusian Flux parinfo[1].limits = [5577.4, 5577.58] parinfo[2].limits = [0.,0.2] parinfo[3].limits = [1e-20,1e-10] ; Terrestrial Flux parinfo[4].limits = [5577.2, 5577.4] start = [DOUBLE(equ_venus_flux_max), 5577.57D, 0.079D, DOUBLE(equ_terr_flux_max*1.05), 5577.33D] ;INITIAL GUESSES, Flux 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(equator_flux),1) ;this used to have 2048 instead of N_ELEMENTS fluxfittell = fltarr(N_ELEMENTS(equator_flux),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 ---------------------------------------- ;------------------------------------------------------------------------------------------------------------ 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,.88,.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='16July2012_greenline_equator_combined_flux.ps',/COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES DEVICE,FILENAME=strcompress(string(day[0])+month[0]+'2012_greenline_'+position[0]+'_combined_flux.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_flux[*], equator_flux_err[*], start, PARINFO=parinfo, $ NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1) plot, lambda1_flux, equator_flux[*], xrange=[5576.5, 5578.2], /xstyle, $ yrange=[1.4* MIN(equator_flux[*]), 1.2*MAX(equator_flux[*])], $ 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,lambda,hartley_oct1[*,i],hartley_oct1_err[*,i] ;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=1,COLOR=250,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 - 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_'+exp_numb[k]+'_flux.dat', /GET_LUN ;**** UPDATE VARIABLE **** ; ; for i=0, N_ELEMENTS(lambda)-1 do begin ; printf, lun, lambda1_flux, equator_flux ; endfor ; FREE_LUN, lun ;------------------------------------------------------------------------------------------ DEVICE,/CLOSE ;------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------ ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ************************************** RAYLEIGH UNITS EQUATOR**************************************** ; ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PRINT,'' PRINT, '' print, 'The following calculations will be done in cgs flux units' 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, venus_max*1.2 ] ; Venusian Relative parinfo[1].limits = [5577.4, 5577.55] parinfo[2].limits = [0.,0.08] parinfo[3].limits = [0., terr_max*1.2] ; Terrestrial Relative parinfo[4].limits = [5577.2, 5577.33] start = [DOUBLE(2.0e3), 5577.57D, 0.079D, DOUBLE(1.05*equ_terr_rayleigh_max), 5577.33D] 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 RAYLEIGH UNITS ---------------------------------------- ;------------------------------------------------------------------------------------------------------------ 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,.88,.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]+'2012_greenline_'+position[0]+'_combined_rayleigh.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[*], start, PARINFO=parinfo, $ NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1) plot, lambda, equator_rayleigh[*], xrange=[5576.5, 5578.2], /xstyle, $ yrange=[1.4*MIN(equator_rayleigh[*]), 1.2*MAX(equator_rayleigh[*])], $ 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=[1.2*MIN(equator_rayleigh)/airmass_equator, 1.2*MAX(equator_rayleigh)/airmass_equator], $ ystyle=1, charsize=1.5 xloc = !X.Window[1] + 0.09 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,hartley_oct1[*,i],hartley_oct1_err[*,i] ;plots the error. wavelength = findgen(3000)/1000 + 5576 ;--------------------------------------------------------------------------------------- ; 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. ;----------------------------------------------------------------------------------------- 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=1,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 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_'+exp_numb[k]+'_flux.dat', /GET_LUN ;**** UPDATE VARIABLE **** ; ; for i=0, N_ELEMENTS(lambda)-1 do begin ; printf, lun, lambda1_flux, equator_flux ; 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 cgs flux units' 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, 2] ; Venusian Relative parinfo[1].limits = [5577.4, 5577.55] parinfo[2].limits = [0.,0.079] parinfo[3].limits = [0., 1.0] ; Terrestrial Relative parinfo[4].limits = [5577.2, 5577.33] start = [DOUBLE(equ_venus_norm_max), 5577.56D, 0.079D, DOUBLE(equ_terr_norm_max), 5577.33D] ;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 ---------------------------------------- ;------------------------------------------------------------------------------------------------------------ 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,.88,.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]+'2012_greenline_'+position[0]+'_combined_norm.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[*], start, PARINFO=parinfo, $ NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1) plot, lambda, equator_norm[*], xrange=[5576.5, 5578.2], /xstyle, $ yrange=[1.1*MIN(equator_norm[*]), 1.2*MAX(equator_norm[*])], $ 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 ;--------------------------------------------------------------------------------------- ; 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,hartley_oct1[*,i],hartley_oct1_err[*,i] ;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=1,COLOR=250,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 - 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_'+exp_numb[k]+'_flux.dat', /GET_LUN ;**** UPDATE VARIABLE **** ; ; for i=0, N_ELEMENTS(lambda)-1 do begin ; printf, lun, lambda1_flux, equator_flux ; endfor ; FREE_LUN, lun ;------------------------------------------------------------------------------------------ DEVICE,/CLOSE ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ;------------------------------------------------------------------------------------------------------------------ ; ********************************************* SOUTH LIMB PLOTS *************************************************** ;------------------------------------------------------------------------------------------------------------------ ; ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************************** FLUX UNITS ************************************************* ; ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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' PRINT,'' PRINT, '' parinfo = replicate({fixed:0, limited:[1,1], limits:[1e-20,1e6]},5) parinfo[0].limits = [1e-20,1e-10] ; Venusian Flux parinfo[1].limits = [5577.4, 5577.58] parinfo[2].limits = [0.,0.2] parinfo[3].limits = [1e-20,1e-10] ; Terrestrial Flux parinfo[4].limits = [5577.2, 5577.4] start = [DOUBLE(south_venus_flux_max), 5577.57D, 0.079D, DOUBLE(1.05*south_terr_flux_max), 5577.33D] ;INITIAL GUESSES, Flux 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(equator_flux),1) ;this used to have 2048 instead of N_ELEMENTS fluxfittell = fltarr(N_ELEMENTS(equator_flux),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 ---------------------------------------- ;------------------------------------------------------------------------------------------------------------ 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,.88,.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[2])+month[2]+'2012_greenline_'+position[2]+'_combined_flux.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, south_flux[*], south_flux_err[*], start, PARINFO=parinfo, $ NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1) plot, lambda1_flux, south_flux[*], xrange=[5576.5, 5578.2], /xstyle, $ yrange=[1.4* MIN(south_flux[*]), 1.2*MAX(south_flux[*])], $ 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[2]+' 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,lambda,hartley_oct1[*,i],hartley_oct1_err[*,i] ;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=1,COLOR=250,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 - 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(south_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_'+exp_numb[k]+'_flux.dat', /GET_LUN ;**** UPDATE VARIABLE **** ; ; for i=0, N_ELEMENTS(lambda)-1 do begin ; printf, lun, lambda1_flux, equator_flux ; endfor ; FREE_LUN, lun ;------------------------------------------------------------------------------------------ DEVICE,/CLOSE ;------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------ ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ************************************** RAYLEIGH UNITS SOUTH**************************************** ; ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PRINT,'' PRINT, '' print, 'The following calculations will be done in cgs flux units' 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, venus_max*1.2 ] ; Venusian Relative parinfo[1].limits = [5577.4, 5577.55] parinfo[2].limits = [0.,0.08] parinfo[3].limits = [0., terr_max*1.2] ; Terrestrial Relative parinfo[4].limits = [5577.2, 5577.33] start = [DOUBLE(south_venus_rayleigh_max), 5577.57D, 0.079D, DOUBLE(1.05*south_terr_rayleigh_max), 5577.33D] 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 RAYLEIGH UNITS ---------------------------------------- ;------------------------------------------------------------------------------------------------------------ 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,.88,.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[2])+month[2]+'2012_greenline_'+position[2]+'_combined_rayleigh.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, south_rayleigh[*], south_rayleigh_err[*], start, PARINFO=parinfo, $ NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1) plot, lambda, south_rayleigh[*], xrange=[5576.5, 5578.2], /xstyle, $ yrange=[1.4*MIN(south_rayleigh[*]), 1.2*MAX(south_rayleigh[*])], $ ystyle=9, psym=-4,xtitle='Wavelength ('+angstrom+' )', $ ytitle='Venusian Intesnsity (R/'+angstrom+')', title=strcompress(month[0]+' '+string(day[0])+', 2012 - '+position[2]+' 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=[1.2*MIN(south_rayleigh)/airmass_south, 1.2*MAX(south_rayleigh)/airmass_south], $ ystyle=1, charsize=1.5 xloc = !X.Window[1] + 0.09 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,hartley_oct1[*,i],hartley_oct1_err[*,i] ;plots the error. wavelength = findgen(3000)/1000 + 5576 ;--------------------------------------------------------------------------------------- ; 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. ;----------------------------------------------------------------------------------------- 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=1,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 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(south_rayleigh)) ;!!! chisquared = 0.0 totresidual = 0.0 for j=0,N_ELEMENTS(south_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(south_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_'+exp_numb[k]+'_flux.dat', /GET_LUN ;**** UPDATE VARIABLE **** ; ; for i=0, N_ELEMENTS(lambda)-1 do begin ; printf, lun, lambda1_flux, equator_flux ; endfor ; FREE_LUN, lun ;------------------------------------------------------------------------------------------ DEVICE,/CLOSE ;------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------ ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ************************************** NORMALIZED UNITS SOUTH**************************************** ; ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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' PRINT,'' PRINT, '' ; terr_max = -1 ; FOR n=16, 19 DO BEGIN ; if south_norm[n] GT terr_max THEN terr_max = south_norm[n] ; ENDFOR ; venus_max = -1 ; FOR n=19, 23 DO BEGIN ; if south_norm[n] GT venus_max THEN venus_max = south_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, 2] ; Venusian Relative parinfo[1].limits = [5577.4, 5577.55] parinfo[2].limits = [0.,0.8] parinfo[3].limits = [0., 1.0] ; Terrestrial Relative parinfo[4].limits = [5577.2, 5577.33] start = [DOUBLE(south_venus_norm_max), 5577.56D, 0.079D, DOUBLE(south_terr_norm_max), 5577.33D] ;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(south_norm),1) ;this used to have 2048 instead of N_ELEMENTS fluxfittell = fltarr(N_ELEMENTS(south_norm),1) residual = fltarr(N_ELEMENTS(south_norm),1) totresidual = fltarr(1) chisquared = fltarr(1) redchi = fltarr(1) averesid = fltarr(1) ew5577 = fltarr(1) ewerr5577 = fltarr(1) ;------------------------------------------------------------------------------------------------------------ ; ---------------------------------- PLOTTING IN NORMALIZED UNITS SOUTH-------------------------------------- ;------------------------------------------------------------------------------------------------------------ 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,.88,.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[2])+month[2]+'2012_greenline_'+position[2]+'_combined_norm.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, south_norm[*], south_norm_err[*], start, PARINFO=parinfo, $ NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1) plot, lambda, south_norm[*], xrange=[5576.5, 5578.2], /xstyle, $ yrange=[1.4* MIN(south_norm[*]), 1.2*MAX(south_norm[*])], $ ystyle=1, psym=10,xtitle='Wavelength ('+angstrom+' )', $ ytitle='Normalized Units', title=strcompress(month[0]+' '+string(day[0])+', 2012 - '+position[2]+' 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,lambda,hartley_oct1[*,i],hartley_oct1_err[*,i] ;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=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=1,COLOR=250,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 - 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(south_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(Venus_solsub) ;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_'+exp_numb[k]+'_flux.dat', /GET_LUN ;**** UPDATE VARIABLE **** ; ; for i=0, N_ELEMENTS(lambda)-1 do begin ; printf, lun, lambda1_flux, equator_flux ; endfor ; FREE_LUN, lun ;------------------------------------------------------------------------------------------ DEVICE,/CLOSE STOP end