;------------------------------------------------------------------------------------------ ; Compile: .r Candace_Venus_6300.pro ; execute: Candace_Venus_6300 ;------------------------------------------------------------------------------------------ ; Start with 'function function_name'. The function will use the ; variables followed after the function name defined (lambda, A, F, ; and pder for the gfunct function) ;------------------------------------------------------------------------------------------ ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; **************************** SECTION 1 - GAUSSIAN FUNCITONS *********************************************** ; ; Gaussian functions defined and called on in the code for fitting. ; Two functions are calculated, one for a blended Gaussian ; consisting of two Gaussians, and another of a single Gaussian ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;------------------------------------------------------------------------------------------ ; function gfunct ; define Gaussian function for deblend case. ; This will be called later on in the code and the variable F will be ; returned. ; A[] array is defined near the bottom of the code in: ; parinfo[*].value = [2.31e-13, 5577.01, 0.079, 2.11e-12, 5577.232]; **** UPDATE VALUES ***** ; ; A[0] = Amplitude of Venusian green line ; A[1] = central wavelength of of Venusian green line ; A[2] = sigma (spread in the data, defined by the dispersion of the instrument, same for both lines) ; A[3] = Amplitude of the Terrestrial green line ; A[4] = central wavelength of the Terrestrial green line ; ;------------------------------------------------------------------------------------------ function gfunct,lambda, A, F, pder ;define Gaussian function for deblend case F=A[0]*exp(-(lambda-A[1])^2.0/(2.0*(A[2]^2.0)))+A[3]*exp(-(lambda-A[4])^2.0/(2.0*(A[2]^2.0))) IF N_PARAMS() GE 4 THEN $ pder = [[exp(-(lambda-A[1])^2/(2*(A[2]^2)))], [(lambda-A[1])*A[0]/A[2]^2*exp(-(lambda-A[1])^2/(2*(A[2]^2)))],[(lambda-A[1])*A[0]/A[2]^3*exp(-(lambda-A[1])^2/(2*(A[2]^2)))]] ; 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[5]^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 *********************************************** ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;------------------------------------------------------------------------------------------ ; Open text files with wavelength and flux of standard star. ; This is the TRUE flux, not counts. You can get a data file of ; wavelength vs true flux off the APO website. If you didn't use one ; of the stars on their list, you can just download a star of the same ; type (AOV) and scale the counts of your observed flux standard so it ; matches that of the true flux. ; Magnitude of downloaded star (NOT THE ONE YOU OBSERVED) and scaled ; to the magnitude of the one you observed. ; ex, Vega has mag 0 and star observered was 5, divide by 100. ; Scale goes as 2.512^(magnitude_difference) ;------------------------------------------------------------------------------------------ A0V_mag = 0.0 observed_A0V_mag = 6.32 ;**** UPDATE VARIABLE **** diff = A0V_mag - observed_A0V_mag scale = 2.512^diff openr, lun, 'flux_stand_calib.dat', /GET_LUN ;read in flux calibration file **** UPDATE VARIABLE **** fluxfile=fltarr(2,13) readf, lun, fluxfile FREE_LUN, lun fluxlambda=double(reform(fluxfile[0,*])) flux=double(reform(fluxfile[1,*])) flux = flux*scale ;------------------------------------------------------------------------------------------ ; Text file of airmass and observing time of object ;------------------------------------------------------------------------------------------ openr, lun, './airmass.dat', /GET_LUN; read in airmasses **** UPDATE VARIABLE **** airmass=fltarr(1) readf, lun, airmass FREE_LUN, lun openr, lun, './time.dat', /GET_LUN; read in time **** UPDATE VARIABLE **** time=fltarr(1) readf, lun, time FREE_LUN, lun ;------------------------------------------------------------------------------------------ ; read in Venus night data ; Original code has it read a text file only containing the 47th order ;(the one containing the green line). ; Reads in data file as a lun "logical unit number" ; Venus is a 2D float array with 2 coloumns of length 1596 rows (the ; wavelength and flux). The values of flux and wavelength (lun) are ; input into the Venus array. Each column is then written to a single ; 1D array, lambda and Venus_night. ;------------------------------------------------------------------------------------------ openr, lun, 'Venusnight_54.0008.dat', /GET_LUN ;**** UPDATE VARIABLE **** Venus = fltarr(2,2048) 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 ;------------------------------------------------------------------------------------------ ; read in Venus day side data ; Original code has it read a text file only containing the 47th order ; (the one containing the green line). ;------------------------------------------------------------------------------------------ openr, lun, 'Venusday_54_set2.dat', /GET_LUN ;**** UPDATE VARIABLE **** Venus = fltarr(2,2048) readf, lun, Venus FREE_LUN, lun lambda_day = double(reform(Venus[0,*])) ; wavelength value from Venus array (flux and wavelength) Venus_day = double(reform(Venus[1,*])+1.0) ; Flux values from Venus array ;------------------------------------------------------------------------------------------ ;read in solar analog ref data ;------------------------------------------------------------------------------------------ openr, lun, 'solar_analog.0098.dat', /GET_LUN ; **** UPDATE VARIABLE **** solar = fltarr(2,2048) readf, lun, solar FREE_LUN, lun solarflux = double(reform(solar[1,*])+1.0) ;------------------------------------------------------------------------------------------ ;read in fluxcal data (flux standard star) ;------------------------------------------------------------------------------------------ openr, lun, 'A0V_54.0003.dat', /GET_LUN ; **** UPDATE VARIABLE **** fluxcal=fltarr(2,2048) readf, lun, fluxcal FREE_LUN, lun fluxcalflux=double(reform(fluxcal[1,*])+1.0) ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************** SECTION 3 - ERROR CALCULATION ********************************************* ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;------------------------------------------------------------------------------------------ ; calculate errors ; Interpolate to add more points to the flux standard file so it has ; the number of points as the standard star observed. ; The value used in calculating the error (3.8, 5, and 3.39 are: ; "3.8 is the gain of the CCD, 3.39 is the read noise in electrons, and 5 is ; an estimate of how many pixels are combined to form each data point in the ; spectrum (i.e. what the spatial extent of the slit is in pixels). These ; numbers should not be changed unless specs for the instrument ; change." -Adam ;------------------------------------------------------------------------------------------ ;***********************NOTE***************************************************** ; The subtraction of the day side from the night side will cause ; negative values. When the error of the nightside is caluculated ; this will produce NAN's since you take a square root. Therefore, the ; minimum value is found from the spectrum and added to the data plus ; a little more just to make sure there are no zero's. ; This was my adjustment and not in the orignal code. I think this ; needs to be subtracted off afterward... ;******************************************************************************* Venus_night_min = MIN(Venus_night) Venus_day_min = MIN(Venus_day) IF Venus_night_min LT 0 THEN BEGIN Venus_night = Venus_night + abs(Venus_night_min) + abs(Venus_night_min*0.01) ENDIF IF Venus_day_min LT 0 THEN BEGIN Venus_day = Venus_day + abs(Venus_day_min) + abs(Venus_day_min*0.01) ENDIF Venus_night_err = sqrt(Venus_night/3.8 + 5.0 * 3.39) ;See description above for description of numbers Venus_day_err = sqrt(Venus_day/3.8 + 5.0 * 3.39) solarflux_err = sqrt(solarflux/3.8 + 5.0 * 3.39) fluxcalflux_err = sqrt(fluxcalflux/3.8 + 5.0 * 3.39) fluxcalinterp = interpol(flux, fluxlambda, lambda, /spline) ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************** SECTION 4 - FLUX CALIBRATION ********************************************* ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;------------------------------------------------------------------------------------------ ;perform flux calibration and propogate error for nightside and day ;side of Venus. ; ; ****************************** NOTE ******************************* ; The component of exp(-airmass[0])/exp(-1.10)) is the ; airmass of the object and the flux standard star. As my stars are ; near the horizon to match Venus, they will be at a much higher ; airmass than 1.1. MAKE SURE TO CHANGE THIS TO: ; ; exp(-airmass[0])/exp(-star_airmass)) ; Venus_fluxcal = Venus_night * fluxcalinterp * 0.71 * 1e-16 / (3.0 * fluxcal[0,*] * exp(-airmass[0])/exp(-1.10)) ; 0.71 = seeing, less than 1" can set to 1 (amount of flux in slit), ; 1e-16 = flux conversion factor ; 3.0 = scales for exposure time (object_exposure/star_exposure) ; 1.1 = airmass of star ; ; *********************NOTE!!!****************************** ; Due to the subtraction of the day side from the night side, there ; will be negative values of flux. This will cause the calculation of ; Venus_flux_cal_err to crash. In order to get around this, the flux ; is shifted up by an amount equal to the largest negative value. This ; can then be shifted back down later. This was done when Venus_night ; and Venus_night_err were first defined. ; ;------------------------------------------------------------------------------------------ Venus_fluxcal_night = Venus_night * fluxcalinterp * 0.6 * 1e-16 / (2.0 * fluxcal[1,*] * exp(-airmass[0])/exp(-2.5)) ;******* UPDATE VALUES: seeing and exposure time ***** Venus_fluxcal_day = Venus_day * fluxcalinterp * 0.6 * 1e-16 / (0.0017 * fluxcal[1,*] * exp(-airmass[0])/exp(-2.5)) Venus_fluxcal_night_err = abs(Venus_fluxcal_night) * sqrt((Venus_night_err/Venus_night)^2.0 + (fluxcalflux_err/fluxcal[1,*])^2.0) Venus_fluxcal_day_err = abs(Venus_fluxcal_day) * sqrt((Venus_day_err/Venus_day)^2.0 + (fluxcalflux_err/fluxcal[1,*])^2.0) ;******************************************************************************************************** ; Changed (fluxcal_err/fluxcal)^2.0) to (fluxcalflux_err/fluxcal)^2.0) since fluxcal_err is not defined ;******************************************************************************************************* ;------------------------------------------------------------------------------------------ ; Doppler shift telluric corrected solar standard ; change to Venus wavelength scale (i.e. lambda[i] is the same for ; both). ; Original equations are: ; ; solarlambda = lambda + 63.0/3e5 * lambda ; Doppler shift telluric corrected solar standard. ; Sets the wavelength for the solar analog ; solardopcor = interpol(solartelldiv, solarlambda, lambda, /spline) ; change to comet wavelength scale (i.e. lambda[i] is the same for both) ; solardopcorerr = interpol(solartelldiv_err, solarlambda, lambda, /spline) ; ########################## IMPORTANT ############################ ; Make sure to change the Doppler shift. It's the ; Doppler shift difference of the star and object. Make sure you ; check the signs!! I should rewrite this so the fraction is a ; variable read in at the beginning of the code. ; ################################################################# ; Doppler shift telluric corrected solar standard. Lambda is the ; Venusian wavelength, 63 is the velocity of Venus at the time ; of observations and 3e5 is the speed of light. ; However, the telluric features can be taken out by subtracting off ; the dayside here rather than reading in data where the day and night ; have already been subtracted. Since the day and night exposures ; have the same wavelength and Doppler shift, there's no need to ; interpolate the wavelengths or Doppler shift them. ; ; ******************** NOTE **************************************** ; If you don't have a solar standard or the one you have doesn't ; match up well with the Venus spectra, then the fronhaugher lines ; can be taken off by subtracted off the dayside spectra. ; ****************************************************************** ;------------------------------------------------------------------------------------------ ; solarlambda = lambda + 63.0/3e5 * lambda ; Doppler shift telluric corrected solar standard. ; ; Sets the wavelength for the solar analog ; solardopcor = interpol(solartelldiv, solarlambda, lambda, /spline) ;change to comet wavelength scale (i.e. lambda[i] is the same for both) ;solardopcorerr = interpol(solartelldiv_err, solarlambda, lambda, /spline) ;------------------------------------------------------------------------------------------ ; Perform solar subtraction (don't need to flux calibrate since it is ; scaled to comet flux before subtraction anyway) ; ; ***************************** NOTE ************************************* ; Initial equations were: ; Venus_solsub = Venus_fluxcal - solardopcor * median(Venus_fluxcal/solardopcor) ; Venus_solsub_err = sqrt((Venus_fluxcal_err)^2.0 + (solardopcorerr * median(Venus_fluxcal/solardopcor))^2.0) ; ; However, these were changed since the solar lines were taken out by ; subtracting the Venus dayside rather than a solar analog. ; Solardopcor (doppler shifted solar spectra) is now the Venus day ; side. No need to Doppler shift to match the night side. ; ****************************************************************************** ;------------------------------------------------------------------------------------------ ;----------------------------------------------------------------------------------------- ; Using the flux calibrated dayside to subtract ; Need to subtract the dayside from the nightside, but to do that you ; need to scale down the dayside flux. This is done by finding the ; average night/day ratio (find ratio at each wavelength). However, ; a median of this value will make the ratio too large and too much ; dayside flux is subtracted. To correct for this, I added a loop ; that find the ratio and calculates the standard devation of that ; ratio. Ratio values greater than 6 sigma of the average are ; ignored and a new average is calculated. Thus, you have a more ; acurate flux ratio average. The clipping range should be looked ; over in order to find the best fit. Increasing the clipping ; (ie. from 5*stand_dev to 6*stand_dev) will lower the flux of the data. ;----------------------------------------------------------------------------------------- ratio = Venus_fluxcal_night/Venus_fluxcal_day ratio_avg = avg(double(Venus_fluxcal_night/Venus_fluxcal_day),/NAN) stand_deviation = stddev(double(ratio),/NAN) print, "Begin sigma clipping. If continuum doesn't fit on 0, then change clipping range" FOR i=0, 2047 DO BEGIN IF ratio(i) GT (ratio_avg + 12.0*stand_deviation) then begin ;************ UPDATE VALUES *************** ratio(i) = 0 ENDIF IF ratio(i) LT (ratio_avg - 12.0*stand_deviation) then begin ;************ UPDATE VALUES *************** ratio(i) = 0 ENDIF ENDFOR ratio_avg_clip = avg(double(ratio), /NAN) ;ratio_avg_clip = median(double(ratio)) Venus_solsub = Venus_fluxcal_night - Venus_fluxcal_day * ratio_avg_clip Venus_solsub_err = sqrt((Venus_fluxcal_night_err)^2.0 + (Venus_fluxcal_day_err * ratio_avg_clip)^2.0) ; Venus_solsub_err = sqrt((Venus_fluxcal_night_err)^2.0 + (Venus_fluxcal_night_err * ratio_avg_clip)^2.0) ; Venus_solsub = Venus_fluxcal_night - Venus_fluxcal_day * avg(double(Venus_fluxcal_night/Venus_fluxcal_day),/NAN) ; Original line. Takes Median of night/day. Going to try average(mean) and mode(most frequent value) to see if the continum fit is better ; Venus_solsub = Venus_fluxcal_night - Venus_fluxcal_day * median(Venus_fluxcal_night/Venus_fluxcal_day) ; Venus_solsub_err = sqrt((Venus_fluxcal_night_err)^2.0 + (Venus_fluxcal_night_err * median(Venus_fluxcal_night/Venus_fluxcal_day))^2.0) ;---------------------------------------------------------------------------------------------------------- ; Writing out flux calibrated, day side subtracted data and wavelenght ; out to a 2-D array called Hartley_Oct1 (THIS NEEDS TO BE RENANED TO ; SOMETHING NOT A COMET!) for range around green line. ; ; However, before it is written out, the NAN's, -NAN's, and INFINITIES ; are removed (they cause problems with the error estimates when the ; fit is done. This was not in Adam's original code and it really ; helps beat down the error. He must not have had this problem. ; ; The first set of equations (with WHERE) writes out the wavelength ; (lambda), error in flux calibrated, day subtracted data ; (Venus_solsub_err), and flux calibrated, day subtracted data to the ; same name but only writes out values that are finite. ; ; Next, we only write out values that are withing the wavelength range ; we want (defined to be betweeen 5576 - 5578, but you can change ; that). You can also just do this write out and not do the first ; three equations where you write out only the finite values. ; ; These values are then written out into the 2-d array (Hartly_Oct1 ; and Hartley_Oct1_err). BUT!!!!!!!!!! now the array will only be a ; few points in size (about 21 points.), so the loop following this ; section had to be altered to accout for that. ;---------------------------------------------------------------------------------------------------------- ; Hartley_Oct1 = [Venus_solsub] ; Hartley_Oct1_err = [Venus_solsub_err] ;remove all nans ; lambda = lambda[WHERE(FINITE(Venus_solsub))] ; Venus_solsub_err = Venus_solsub_err[WHERE(FINITE(Venus_solsub))] ; Venus_solsub = Venus_solsub[WHERE(FINITE(Venus_solsub))] ;choose lambda around the red line 6298 A to 6302A linerange = WHERE((lambda GE 6298)AND(lambda LE 6302)) ;********* UPDATE WAVELENGTH RANGE FOR ANALYSIS***** lambda = lambda[linerange] Venus_solsub = Venus_solsub[linerange] Venus_solsub_err = Venus_solsub_err[linerange] Hartley_Oct1 = [Venus_solsub] Hartley_Oct1_err = [Venus_solsub_err] ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************** SECTION 5 - GAUSSIAN FITTING TO DATA*************************************** ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;------------------------------------------------------------------------------------------ ; 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.0709, 2.11e-12, 5577.232]; **** UPDATE VALUES ***** ; parinfo[*].value = [4.5e-13, 5577.01, 0.0709, 2.3e-12, 5577.23, 0.0709]; **** UPDATE VALUES ***** ; parinfo[*].value = [2e-14, 5577.13, 0.0709, 9e-14, 5577.31, 0.0709]; **** UPDATE VALUES ***** ;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-5] parinfo[1].limits = [6300.45, 6300.8] parinfo[2].limits = [0.,100.] parinfo[3].limits = [1e-20,1e-5] parinfo[4].limits = [6300.0, 6300.4] start = [DOUBLE(2.31e-13), 6300.55D, 0.0803D, DOUBLE(2.11e-12), 6300.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(Venus_solsub),1) ;this used to have 2048 instead of N_ELEMENTS fluxfittell = fltarr(N_ELEMENTS(Venus_solsub),1) residual = fltarr(N_ELEMENTS(Venus_solsub),1) totresidual = fltarr(1) chisquared = fltarr(1) redchi = fltarr(1) averesid = fltarr(1) ew5577 = fltarr(1) ewerr5577 = fltarr(1) ;--------------------------------------- ;original array sizes: ;fluxfit=fltarr(1651,1) ;fluxfittell=fltarr(1651,1) ;residual=fltarr(1651,1) ;totresidual=fltarr(1) ;chisquared=fltarr(1) ;redchi=fltarr(1) ;averesid=fltarr(1) ;ew5577=fltarr(1) ;ewerr5577=fltarr(1) ;--------------------------------------- ;------------------------------------------------------------------------------------------ ;fit 5577 nm line ;------------------------------------------------------------------------------------------ print, 'begin loop' LOADCT,39,/SILENT !P.CHARTHICK=2 !P.POSITION=[.13,.1,.95,.95] !X.THICK=2 !Y.THICK=2 SET_PLOT,'PS' DEVICE,FILENAME='redlinefit_pole.0011_dayset2.ps',/COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES 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 ; weight = FLTARR(n_elements(Hartley_Oct1[*,i])) ; weight[WHERE(FINITE(weight))] = 1 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) print,'guesses:',start plot, lambda, Hartley_Oct1[*,i], xrange=[6298.0, 6302.0], psym=-4,xtitle='wavelength',ytitle='flux',THICK=4 ;plots the data *********** UPDATE WAVELENGTH PLOT RANGE TO BE SAME AS ANALYSIS RANGE ****** ;--------------------------------------------------------------------------------------- ; 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. ;oplot, lambda, Venus_solsub+Venus_night_err ;oplot, lambda, Venus_solsub-Venus_night_err wavelength = findgen(3000)/1000 + 6299 ;--------------------------------------------------------------------------------------- ; 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. ;-------------------------------------------------------------------------------------- x_shift = 0.0 ;if wavelength calibration isn't perfect, this is the shift so the green line and the fit match y_shift1 = 0.0 oplot, wavelength - x_shift, ((P[0,i] * exp(-(wavelength-P[1,i])^2.0 / (2.0 * (P[2,i]^2.0))) + P[3,i] * exp(-(wavelength-P[4,i])^2.0 / (2.0 * (P[2,i]^2.0))))) + y_shift1, linestyle=2, COLOR=90,THICK=4 ;oplot, wavelength - x_shift, ((P[0,i] * exp(-(wavelength-P[1,i])^2.0 / (2.0 * (P[2,i]^2.0))) + P[3,i] * exp(-(wavelength-P[4,i])^2.0 / (2.0 * (P[5,i]^2.0))))) + y_shift1, linestyle=2, COLOR=90,THICK=4 ; THIS IS WHAT IT SHOULD BE;;;;;;; ; Terrestrial plus green line gaussians ; oplot, wavelength,(P[0,i] * exp(-(wavelength-P[1,i])^2.0 / (2.0 * (P[2,i]^2.0))) + P[3,i] * exp(-(wavelength-P[4,i])^2.0 / (2.0 * (P[5,i]^2.0)))), linestyle=5 ;--------------------------------------------------------------------------------------- ; 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. ;--------------------------------------------------------------------------------------- ; y_shift2 = -0.4e-10 y_shift2 = 0 oplot, wavelength - x_shift, (P[3]*exp(-(wavelength-P[4])^2.0/(2.0*(P[2]^2.0))))+ y_shift2, linestyle=1,COLOR=250,THICK=4 ;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[*,i] = P[0,i] * exp(-(lambda - P[1,i])^2.0 / (2.0 * (P[2,i]^2.0))) fluxfittell[*,i] = P[3,i] * exp(-(lambda - P[4,i])^2.0 / (2.0 * (P[2,i]^2.0))) residual[*,i] = (fluxfit[*,i] + fluxfittell[*,i] - (Hartley_Oct1[*,i])) / (Hartley_Oct1_err[*,i]) ;fluxfit[*,i]=P[0,i]*exp(-(lambda-P[1,i])^2.0/(2.0*(P[2,i]^2.0))) ;fluxfittell[*,i]=P[3,i]*exp(-(lambda-P[4,i])^2.0/(2.0*(P[2,i]^2.0))) ;residual[*,i]=(fluxfit[*,i]+fluxfittell[*,i]-Hartley_Oct1[*,i])/(Hartley_Oct1_err[*,i]) ; fiterr=replicate(0.0,13) fiterr=replicate(0.0,N_ELEMENTS(Venus_solsub)) ;!!! chisquared[i] = 0.0 totresidual[i] = 0.0 ; for j=1019, 1031 do begin ;range of point in lambda to match greenline region for j=0,N_ELEMENTS(Venus_solsub)-1 do begin ;range of point in lambda to match greenline region totresidual[i] = totresidual[i] + abs(residual[j,i]) chisquared[i] = chisquared[i] + residual[j,i]^2.0 endfor ; print, residual[*,i] redchi[i] = (chisquared[i])/4.0 ew5577[i] = 2.50663 * P[0,i] * P[2,i] if sqrt(redchi[i]) gt 1.0 then begin ewerr5577[i] = 2.50663 * perror[0] * P[2,i] * sqrt(redchi[i]) endif else begin ewerr5577[i] = 2.50663 * perror[0] * P[2,i] endelse ; averesid=totresidual/13.0 ;divided by the number of points you're sampling over (j=1022, 1035) averesid=totresidual/N_ELEMENTS(Venus_solsub) ;divided by the number of points you're sampling over (j=1022, 1035) print, 'average residual:', averesid[i] print, 'reduced chi-squared:',redchi[i] print, 'integrated flux:', ew5577[i] print, 'error in flux:', ewerr5577[i] endfor DEVICE,/CLOSE STOP end