;------------------------------------------------------------------------------------------ ; Compile: .r Candace_Venus_5577_part17_publish.pro ; execute: Candace_Venus_5577_part17_publish ;------------------------------------------------------------------------------------------ ; Start with 'function function_name'. The function will use the ; variables followed after the function name defined (lambda, A, F, ; and pder for the gfunct function) ;------------------------------------------------------------------------------------------ ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; **************************** SECTION 1 - GAUSSIAN FUNCITONS *********************************************** ; ; Gaussian functions defined and called on in the code for fitting. ; Two functions are calculated, one for a blended Gaussian ; consisting of two Gaussians, and another of a single Gaussian ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;------------------------------------------------------------------------------------------ ; function gfunct ; define Gaussian function for deblend case. ; This will be called later on in the code and the variable F will be ; returned. ; A[] array is defined near the bottom of the code in: ; parinfo[*].value = [2.31e-13, 5577.01, 0.079, 2.11e-12, 5577.232]; **** UPDATE VALUES ***** ; ; A[0] = Amplitude of Venusian green line ; A[1] = central wavelength of of Venusian green line ; A[2] = sigma (spread in the data, defined by the dispersion of the instrument, same for both lines) ; A[3] = Amplitude of the Terrestrial green line ; A[4] = central wavelength of the Terrestrial green line ; ;------------------------------------------------------------------------------------------ function gfunct,lambda, A, F, pder ;define Gaussian function for deblend case F=A[0]*exp(-(lambda-A[1])^2.0/(2.0*(A[2]^2.0)))+A[3]*exp(-(lambda-A[4])^2.0/(2.0*(A[2]^2.0))) IF N_PARAMS() GE 4 THEN $ pder = [[exp(-(lambda-A[1])^2/(2*(A[2]^2)))], [(lambda-A[1])*A[0]/A[2]^2*exp(-(lambda-A[1])^2/(2*(A[2]^2)))], $ [(lambda-A[1])*A[0]/A[2]^3*exp(-(lambda-A[1])^2/(2*(A[2]^2)))]] return, F END ;------------------------------------------------------------------------------------------ ; function gfunct2 ; define Gaussian function for individual line ;------------------------------------------------------------------------------------------ function gfunct2,lambda, A, F, pder ;define Gaussian function for individual line F=A[0]*exp(-(lambda-A[1])^2.0/(2.0*(A[2]^2.0))) IF N_PARAMS() GE 4 THEN $ pder = [[exp(-(lambda-A[1])^2/(2*(A[2]^2)))], [(lambda-A[1])*A[0]/A[2]^2*exp(-(lambda-A[1])^2/(2*(A[2]^2)))],$ [(lambda-A[1])*A[0]/A[2]^3*exp(-(lambda-A[1])^2/(2*(A[2]^2)))]] return, F END ;------------------------------------------------------------------------------------------ ; function OIflux ; A[0]=CO2 flux contribution, A[1]=co2red, A[2]=h2ored, A[3]=co2green, A[4]=h2ogreen ;------------------------------------------------------------------------------------------ function OIflux, lineratio, A, F; A[0]=CO2 flux contribution, A[1]=co2red, A[2]=h2ored, A[3]=co2green, A[4]=h2ogreen B = A[1] / (A[2]+A[1]) C = (lineratio*A[2] - A[4]) / (A[3] - lineratio*A[1]) F = A[0]/(B*C) return, F END ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************** SECTION 2 - READING IN DATA *********************************************** ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;------------------------------------------------------------------------------------------ ; Read in info on nightside exposures, dayside exposures, the ; synthetic solar and and observed standard stars, and initial guesses ; on the Gaussian fits. ;------------------------------------------------------------------------------------------ READCOL, 'nightside_input_pub.dat', date, radial_vel, venus_night_name, exp_numb, position, venus_night_exp,$ venus_night_airmass, Venus_night_angle, $ FORMAT = 'A,F,A,A,A,F,F,F' READCOL, 'dayside_input_pub.dat', venus_day_name, venus_day_exp, venus_day_airmass, doppler_wind, $ stand_dev_scale_up, stand_dev_scale_down, day_scale_extra, $ FORMAT = 'A,F,F,F,F,F,F' READCOL, 'standard_star_input.dat', rad_vel_shift, synth_solar_stand_dev, solar_extra_scale, $ stand_star, mag, star_exp, star_airmass, trans, trans_err, $ FORMAT = 'F,F,F,A,F,F,F,F,F' READCOL, 'gaussian_input.dat', gaussian_flux_terr, gaussian_flux_venus, gaussian_ray_terr, gaussian_ray_venus, $ gaussian_norm_terr,gaussian_norm_venus, $ FORMAT='F,F,F,F,F,F' READCOL, 'wavelength_input_pub.dat', feature, feature_true_lambda, feature_lambda, order, dispersion, solar_contin_up, $ solar_contin_down, fraunhofer, $ FORMAT = 'A,F,F,I,F,F,F,F' READCOL, 'slit_input.dat', slit_center_distance, disk_radius, $ FORMAT = 'F, F' ;------------------------------------------------------------------------------------------ ; Set up counters for file lengths so that it doesn't have to be specified ; when the file is read in because it really slows the code down. ;------------------------------------------------------------------------------------------ FOR k=0, 1 DO BEGIN ; Number of exposures to be analyzed READCOL, 'flux_stand_calib.dat', counter_flux, FORMAT = 'F' READCOL, venus_night_name[0], counter_venus, FORMAT = 'F' ;read in READCOL, stand_star[0], counter_stand_star, FORMAT = 'F' ;read in ENDFOR ;------------------------------------------------------------------------------------------ plot_upper_bound = 1.0 ; Angstroms added to feature_wavelength to define upper bound of wavelength range when plotting plot_lower_bound = 0.7 ; Angstroms subtracted from feature_wavelength to define lower bound of wavelength range when plotting feature_upper_bound = 0.07 ; Upper wavelength bracket around emission feature; cant be too large or it will grab terrestrail feature_lower_bound = 0.07 ; Lower wavelegnth braket around emission Venus_radius = 6052.0 ;------------------------------------------------------------------------------------------ ; Open file that will contain info on the peak 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 printf, lun_1, '#-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------' printf, lun_1, '# File Position Flux ---------------------- Line Range (A)-------------------- --------- Continuum -------- ----- Green Line peak ----- Venus # Sigma' printf, lun_1, '# Name Type --------- Lower Bounds ------ ------ Upper Bounds----- Mean Stand Dev Terrestrial Venus above noise' printf, lun_1, '#-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------' ;------------------------------------------------------------------------------------------ ; Open file that will contain the integrated flux of the Venusian and ; Terrestrial green line calculated in Rayleighs. ;------------------------------------------------------------------------------------------ openw, lun_2, 'output_integrated_rayleigh.dat', /GET_LUN, width = 250 printf, lun_2, '#---------------------------------------------------------------------------------------------------' printf, lun_2, '# File Position Flux ----------- Rayleigh (Integrated) ------------ ' printf, lun_2, '# Name Type Terrestrial Error Venus Error ' printf, lun_2, '#---------------------------------------------------------------------------------------------------' ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ************************ START MAIN LOOP OF CODE OVER ALL EXPOSURES ************************************** ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% doppler_shift = fltarr(N_elements(date)) Venus_zenith_angle = fltarr(N_elements(date)) sec_theta = fltarr(N_elements(date)) alpha = fltarr(N_elements(date)) slit_dist = fltarr(N_elements(date)) ;FOR k=0, 1 DO BEGIN ; Number of exposures to be analyzed FOR k=0, N_ELEMENTS(venus_night_exp)-1 DO BEGIN ; Number of exposures to be analyzed ;------------------------------------------------------------------------------------------ ; Set wavelength bounds for contiuum area near feature, (this is now defined ; in the input file wavelength_input.dat) used for ; dividing the solar spectrum and normalizing the dayside to the ; nightside. ; ; Second, define the position of a fraunhofer line, so as to scale ; the solar spectrum to the nightside. However, the Fraunhofer lines ; are Doppler shifted so that is corrected. ; ; For the green line, there is a large Franhofer line centered at 5576.095 ; and we define a width of 0.4 of the feature, and then doppler shift it ; so it will match Venus. ;------------------------------------------------------------------------------------------ fraunhofer_lower = (fraunhofer[k] - 0.2) + (radial_vel[k] * fraunhofer[k] / 3.0e5) fraunhofer_upper = (fraunhofer[k] + 0.2) + (radial_vel[k] * fraunhofer[k] / 3.0e5) ;------------------------------------------------------------------------------------------ ; 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) ; ; Here, I used HR7001 as the standard star, which was found from ; ftp://ftp.eso.org/pub/stecf/standards/hststan/fhr7001.dat ; It has a magnitude of 0.0, so if you change which A0V flux calibrated ; star you use, you need to change the magnitude. ; The columns read in are: ; wavelength Flux Flux ; (A) (ergs/cm/cm/s/A * 10**16 ) (Mili-Jy) ; ; But we are not using the mili-Jy flux, obviously. ;------------------------------------------------------------------------------------------ A0V_mag = 0.0 ; Magnitude of flux calibrated standard star HR7001 observed_A0V_mag = mag[k] ;**** 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 for HR7001 ; fluxfile=fltarr(3,794) ; Columns: wavel ; readf, lun, fluxfile ; FREE_LUN, lun openr, lun, 'flux_stand_calib.dat', /GET_LUN ; read in flux calibration file for HR7001 fluxfile = fltarr(2, N_ELEMENTS(counter_flux)) readf, lun, fluxfile FREE_LUN, lun fluxlambda = double(reform(fluxfile[0,*])) flux = double(reform(fluxfile[1,*])) flux = flux*scale ;scales the flux of a synthetic A0V star to the magnitude of the A0V you observed. ;------------------------------------------------------------------------------------------ ; read in Venus night data by doing: ; - READCOL, venus_night_name[k], counter, FORMAT = 'F' ; This readcols the data file and sets the first column to the ; variable counter. This variable won't actually be used except ; to be used to define the length of the data file in the next ; step. ; - openr, lun, venus_night_name[k], /GET_LUN ; 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. ;------------------------------------------------------------------------------------------ ; READCOL, venus_night_name[k], counter, FORMAT = 'F' ;read in openr, lun, venus_night_name[k], /GET_LUN ;**** UPDATE VARIABLE **** Venus = fltarr(2,N_ELEMENTS(counter_venus)) 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). ; There will be a doppler shift between the night and day side ; spectra due to the wind velocity of the clouds towards and away on ; either side of the planet. To correct for this shift, a shift ; constant (shift_day) is added to the spectra of the Venusian ; dayside. An interpolation is done to shift the flux values for the ; wavlength shift. If the wavelenth shift is +, you are moving ; backward in the array, and the shift value will be -. ; ; Venus_day_new = interpol(Venus_day, shift_day, lambda_day, /spline) ; ; This interpolates the flux values at the shifted wavelength ; positions. ; Venus_day_new = shifted flux ; Venus_day = unshifted flux ; shift_day = shifted wavelength ; lambda = nightside wavelength your matching to ;------------------------------------------------------------------------------------------ ; READCOL, venus_day_name[k], counter, FORMAT = 'F' ;read in openr, lun, venus_day_name[k], /GET_LUN ;**** UPDATE VARIABLE **** Venus = fltarr(2,N_ELEMENTS(counter_venus)) 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, 1.0 is added so that the zero's at the ; end of the file won't be zero, saving headacheds for error analysis shift_day = lambda_day + doppler_wind[k] ;sign will be opposite of direction you want to shift Venus_day = interpol(Venus_day, shift_day, lambda, /spline) ;see description above ;------------------------------------------------------------------------------------------ ; read in observed standard star data ;------------------------------------------------------------------------------------------ ; READCOL, stand_star[k], counter, FORMAT = 'F' ;read in openr, lun, stand_star[k], /GET_LUN ; **** UPDATE VARIABLE **** fluxcal=fltarr(2,N_ELEMENTS(counter_stand_star)) readf, lun, fluxcal FREE_LUN, lun fluxcalflux=double(reform(fluxcal[1,*])+1.0) ;------------------------------------------------------------------------------------------ ;read in solar spectrum from BASS 2000 data base ; This is NOT the full solar spectrum. A file from 3200 - 10 would be too large ; for the system to hold. Instead, this file is the solar spectrum from the ; BASS 2000 website of only 6 features +/- 15 Angstroms on either side. These ; features are: ; - 3464.4 ; - 3914.4 ; - 5200.7 ; - 5577.3 ; - 5754.8 ; - 6300.23 ; - 6363.7 ; - 7600 ; synth_solarflux_res = interpol(synth_solarflux, synth_solarlambda, lambda, /spline) ; ; Above is an interpolation to reduce the resolution of the synthetic ; solar spectrum from the BASS 2000 website. By reducing the ; resolution, the width of the line will widen, and the depth will ; decrease (due to conserving full width-half max). This needs to be ; shifted to the nightside of Venus, which will be the Doppler shift ; of Venus itself (not the clouds doppler shift) ;------------------------------------------------------------------------------------------ READCOL, 'solar_spectrum_all.dat', counter, FORMAT = 'F' ;read in openr, lun, 'solar_spectrum_all.dat', /GET_LUN ; **** UPDATE VARIABLE **** synth_solar=fltarr(2,N_ELEMENTS(counter)) readf, lun, synth_solar FREE_LUN, lun synth_solarlambda = double(reform(synth_solar[0,*])) synth_solarflux = double(reform(synth_solar[1,*])+1.0) ;-------------------------------------------------------------------------------- ; Shift the synthetic solar spectrum in wavelength so ; it matches the Doppler shifte Venus. There will still ; be a small additional shift due to the Doppler motion from ; the winds. ; Flux is interpolated so the values are doppler shifted to their ; proper position to match Venus' lines. ;-------------------------------------------------------------------------------- doppler_shift[k] = radial_vel[k] * feature_lambda[k]/ 3e5 synth_solarlambda_shift = synth_solarlambda + doppler_shift[k] + rad_vel_shift[k] ;sign will be opposite the way you want to shift synth_solarflux_res = interpol(synth_solarflux, synth_solarlambda_shift, lambda, /spline) ;---------------------------------------------------------------------------------------------------------- ; Calculate the position of desired Venusian feature given the relative velocity of Venus. ; Based on the direction of the speed, the position will be red or blue shifted of the ; Venusian green line. The upper and lower bounds of the green line wavelength are ; considered. ;---------------------------------------------------------------------------------------------------------- Venus_lambda_guess = feature_lambda[k]+ doppler_shift[k] ; Set Veusian greenline wavelength guess Venus_lambda_upper = feature_lambda[k]+ doppler_shift[k] + feature_upper_bound ; Upper wavelength bracket around emission, ; cant be too large or it will grab terrestrail Venus_lambda_lower = feature_lambda[k]+ doppler_shift[k] - feature_lower_bound ; Lower wavelegnth braket around emission ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ;************************* SECTION 3 - ERROR CALCULATION FROM POISSON STATISTICS *************************** ; ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;------------------------------------------------------------------------------------------------- ; Before doing an error calculation, I'm going to clip the data so we only consider ; a specified range around around the feature in question. This will keep us from ; considering areas around the edges of the order where things get weird. Just so long as ; you look at the spectra and see there are no large spikes positive or negative in the ; range you specify, you should be fine. ; However, it should be noted, the range you chose will be the range where statistics are ; calculated. So, if the range is large, you could have large error bars caused from distant ; features (i.e. the 1 sigma error bars that look larger than they should around the feature. ; This was initally at +/- 12 angstroms, but I changed that to 5, so we are looking a 10 A range ;-------------------------------------------------------------------------------------------------- linerange = WHERE((lambda GE (double(feature_lambda[k] - 6.0)))AND(lambda LE (feature_lambda[k] + 6.0))) ;Specify wavelength range we will analyze lambda = lambda[linerange] ; redefining lambda range, so lambda array is now only the range just specificed Venus_night = Venus_night[linerange] Venus_day = Venus_day[linerange] synth_solarflux = synth_solarflux[linerange] fluxcalflux = fluxcalflux[linerange] synth_solarflux_res = synth_solarflux_res[linerange] ;------------------------------------------------------------------------------------------ ; 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 ;------------------------------------------------------------------------------------------ Venus_night_err = sqrt(abs(Venus_night/3.8 + 5.0 * 3.39)) ;See description above for description of numbers Venus_day_err = sqrt(abs(Venus_day/3.8 + 5.0 * 3.39)) ; It is the error from Poisson statistics = sqrt(counts*other) fluxcalflux_err = sqrt(fluxcalflux/3.8 + 5.0 * 3.39) fluxcalinterp = interpol(flux, fluxlambda, lambda, /spline) ; flux = # of flux points in synthetic A0V star, will only be a few data points ; fluxlambda = # of wavelength data points in synthetic A0V spectra, ; (will only be a few points) ; lambda = # wavelength points in Venus spectra, will be thousands of points ; interpolate to make the number of syntheic A0V spectra have the same number ; of data points as those in the observed data, will will have thousands of points ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ;************************* SECTION 4 - SUBTRACTING SCALED DAY SIDE FROM NIGHT ******************************** ; ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;----------------------------------------------------------------------------------------- ; 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. ;----------------------------------------------------------------------------------------- linerange_upper = feature_lambda[k] + 6.0 linerange_lower = feature_lambda[k] - 6.0 linerange = WHERE((lambda GE linerange_lower) AND (lambda LE linerange_upper)) ;Specify wavelength range we will analyze ratio_linerange = Venus_night[linerange]/Venus_day[linerange] ratio_avg_linerange = avg(double(ratio_linerange)) ; When it comes to adjusting scale by the standard ; deviation, the avg value is much better to use than ; the median value ; ratio_med_linerange = median(double(ratio_linerange)) ratio_stand_dev = stddev(double(ratio_linerange),/NAN) ; standar deviation of the ratio between the wavelengths defined by linerange FOR i=0, N_ELEMENTS(ratio_linerange)-1 DO BEGIN ; done in wavelength range between 5575 - 5579 IF ratio_linerange(i) GT (ratio_avg_linerange + stand_dev_scale_down[k] * ratio_stand_dev) then begin print, ratio_linerange(i) ratio_linerange(i) = 0 ; anything outside the upper bounds set by the standard deviation is set to 0 and not considered ENDIF IF ratio_linerange(i) LT (ratio_avg_linerange - stand_dev_scale_up[k] * ratio_stand_dev) then begin print, ratio_linerange(i) ratio_linerange(i) = 0 ; anything outside the lower bounds set by the standard deviation is set to 0 and not considered ENDIF ENDFOR ; ratio_med_clip = median(double(ratio_linerange)) ratio_avg_clip = avg(double(ratio_linerange)) Venus_day_scaled = Venus_day*(ratio_avg_clip + day_scale_extra[k]) ;This may change for each feature and may need to be included ; in the wavelength_input file instead. ;---------------------------------------------------------------------------------- ; quick test in DIVIDING the scaled dayside from the nightside, then multiplying ; the contiuum of the dayside so as to get the flux units back. ; Nope, that looks like shit, we aren't doing that. It also does a poor job ; of matching the shape of the continuum. So, we are sticking with subtracting. ;---------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------- ; Original calculation of subtracting the scaled dayside. ;---------------------------------------------------------------------------------- Venus_diff = Venus_night - Venus_day*(ratio_avg_clip + day_scale_extra[k]) ; true night side spectrum, after scaled dayside is subtracted Venus_diff_err = sqrt(abs((Venus_night_err)^2.0 + (Venus_day_err * (ratio_avg_clip+day_scale_extra[k]))^2.0)) ; error in residuals, using err = sqrt[(val1)^2 + (val2)^2] ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************** SECTION 5 - FLUX CALIBRATION ********************************************* ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;----------------------------------------------------------------------------------------- ;Calculate flux calibration and error on already subtracted night - day spectra ;----------------------------------------------------------------------------------------- ;------------------------------------------------------------------------------------------ ;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 (transmission) = trans[k]), ; 1e-16 = flux conversion factor ; 3.0 = scales for exposure time (object_exposure/star_exposure) ; 1.1 = airmass of star ; airmass[0] = airmass of Venus ; fluxcal[0,*] = flux of standar star observed ; ; *********************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_diff = Venus_diff * fluxcalinterp * trans[k] * 1e-16 / $ ((venus_night_exp[k]/star_exp[k]) * fluxcalflux * exp(-venus_night_airmass[k])/exp(-star_airmass[k])) Venus_fluxcal_night_diff_err = abs(Venus_fluxcal_night_diff) * $ sqrt((Venus_diff_err/Venus_diff)^2.0 + (fluxcalflux_err/fluxcalflux)^2.0 + (trans_err[k]/trans[k])^2.0) LOADCT,39 SET_PLOT,'PS' DEVICE, FILENAME = date[k]+'_'+feature[k]+'_'+exp_numb[k]+'_'+position[k]+'_day_night_diff.ps',/COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES plot, lambda, Venus_fluxcal_night_diff, xrange=[(feature_lambda[k] - 3), (feature_lambda[k] + 3)], /xstyle, yrange=[MIN(Venus_fluxcal_night_diff)*0.6, MAX(Venus_fluxcal_night_diff)*1.1], /ystyle, psym=0 oplot, [feature_lambda[k]-6,feature_lambda[k]+6], [0, 0], linestyle = 0, color =0 ;------------------------------------------------------------------------------------------------ ; Legend. ; oplot for the lines ; xyouts for the words ;------------------------------------------------------------------------------------------------ oplot, [feature_lambda[k]+0.4, feature_lambda[k]+0.6], [MAX(Venus_fluxcal_night_diff)*1.04, MAX(Venus_fluxcal_night_diff)*1.04], linestyle = 0, color = 0 xyouts, feature_lambda[k]+0.7, MAX(Venus_fluxcal_night_diff)*1.04, 'Night - Scaled day', color = 0 DEVICE,/CLOSE FREE_LUN, lun ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************** SECTION 6 - RESIDUAL SOLAR SUBTRACTION************************************** ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;---------------------------------------------------------------------------------------------------- ; Not all of the solar continuum is subtracted from the subtraction of the dayside, the rest must ; be removed by considering a synthetic solar spectrum, and even then, it won't remove everything. ; In order to scale the synthetic solar spectrum and subtract it from the (night - scaled_day) ; we need to look at the fraunhofer lines. One strong Franhofer line near the feature in question ; is chosen and defined at the begining of the code, along with the feature, which is doppler ; shifted. ; ; The solar spectrum is an absorption spectrum, so the maximum values will be at the conitnuum. ; Therefore, the wavelength range of the synthetic spectrum is defined at the top part of the code ; and the maximum value in that line range is defined as the continuum level. This max value ; is subtracted from the synthetic solar spectrum to make the contniuum be zero. ; ; Once the continuum is zero, we can multiply the the continuum by a constant, so as to set ; the depth of the Franhofer lines equal to the peak in the )night - scaled_day) spectrum. ; To do this, the max value in the Fraunhofer line range for the (night - scaled_day) ; and the min value for the solar spectrum is found. The ratio of the value is found and ; used to scale the solar spectrum, which is then subtracted (by adding the negative ; of the spectrum) from the (night - scaled_day) spectrum. ; ; A plot of the night-day without solar subtraction, with solar subtration, and the ; scaled solar spectrum are made for viewing. ;------------------------------------------------------------------------------------------------ synth_solarflux_zero = synth_solarflux_res - MAX(synth_solarflux_res) ; sets synthetic solar continuum to 0 so you can ; multiply by a constant and change depth of franhofer ; lines without changing continuum linerange = WHERE((lambda GE fraunhofer_lower) AND (lambda LE fraunhofer_upper)) Venus_fluxcal_night_diff_fran = MAX(Venus_fluxcal_night_diff(linerange)) synth_solarflux_fran = MIN(synth_solarflux_zero(linerange)) solar_scale = abs(Venus_fluxcal_night_diff_fran/synth_solarflux_fran) synth_solarflux_scaled = synth_solarflux_zero * solar_scale * solar_extra_scale[k] ; scales Franhofer lines so they match ; the magnitude of the nightside franhofer Venus_solsub = Venus_fluxcal_night_diff + synth_solarflux_scaled Venus_solsub_err = Venus_fluxcal_night_diff_err ;--------------------------------------------------------------------------------- ; plot of the Venus nightside with the sythnetic solar spectrum divided out. ; Use this plot as a guide to makde adjustmetns to the synthetic solar scaling. ;--------------------------------------------------------------------------------- linerange = WHERE((lambda GE (feature_lambda[k]- 6.0))AND(lambda LE (feature_lambda[k]+ 6.0))) ;Specify wavelength range we will analyze flux_max = MAX(Venus_solsub(linerange), index) flux_max_err = Venus_solsub_err(linerange(index)) lambda_max = lambda[linerange (index)] flux_min = MIN(Venus_solsub(linerange), index) flux_min_err = Venus_solsub_err(linerange(index)) lambda_min = lambda[linerange (index)] Venus_day_max = MAX(Venus_day(linerange)) Venus_day_min = MIN(Venus_day(linerange)) Venus_night_max = MAX(Venus_night(linerange)) Venus_night_min = MIN(Venus_night(linerange)) LOADCT,39 SET_PLOT,'PS' DEVICE, FILENAME = date[k]+'_'+feature[k]+'_'+exp_numb[k]+'_'+position[k]+'_solar_scale.ps',/COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES plot, lambda, Venus_solsub, xrange=[(feature_lambda[k] - 3), (feature_lambda[k] + 3)], /xstyle, yrange=[MIN(Venus_fluxcal_night_diff)*1.1, MAX(Venus_fluxcal_night_diff)*1.1], /ystyle, psym=0 oplot, lambda, Venus_fluxcal_night_diff, linestyle = 2, color = 50 oplot, lambda, synth_solarflux_scaled, linestyle = 4, color =250 oplot, [feature_lambda[k]-6,feature_lambda[k]+6], [0, 0], linestyle = 0, color =0 ;--------------------------------------------------------------------------------------- ; Plot vertical lines on where green line and venusian green line should be ; ;--------------------------------------------------------------------------------------- oplot, fltarr(2) + feature_lambda[k], !Y.CRANGE, linestyle=2,COLOR=30,THICK=2 oplot, fltarr(2) + doppler_shift[k] + feature_lambda[k], !y.crange, linestyle=2,COLOR=30,THICK=2 ; doppler_shift[k] = radial_vel[k] ;------------------------------------------------------------------------------------------------ ; Legend. ; oplot for the lines ; xyouts for the words ;------------------------------------------------------------------------------------------------ oplot, [feature_lambda[k]+0.4, feature_lambda[k]+0.6], [MAX(Venus_solsub)*1.04, MAX(Venus_solsub)*1.04], linestyle = 0, color = 0 oplot, [feature_lambda[k]+0.4, feature_lambda[k]+0.6], [MAX(Venus_solsub)*1.00, MAX(Venus_solsub)*1.00], linestyle = 2, color = 50 oplot, [feature_lambda[k]+0.4, feature_lambda[k]+0.6], [MAX(Venus_solsub)*0.96, MAX(Venus_solsub)*0.96], linestyle = 4, color = 250 xyouts, feature_lambda[k]+0.7, MAX(Venus_solsub)*1.04, 'Solar Subtracted', color = 0 xyouts, feature_lambda[k]+0.7, MAX(Venus_solsub)*1.0, 'No Solar Subt', color = 50 xyouts, feature_lambda[k]+0.7, MAX(Venus_solsub)*0.96, 'Scaled Solar', color = 250 DEVICE,/CLOSE FREE_LUN, lun ;------------------------------------------------------------------------------------- ; plot to adjust slight doppler shift between the nightside and dayside. ; This difference is due to the motion of the clouds since it's moving toward you on ; one side and away from you on the other side, since limb measurements are take on ; both. This is done manually using the dayside input file. This should be ; automated at some point, but I'm not sure how. ;------------------------------------------------------------------------------------- LOADCT,39 SET_PLOT,'PS' DEVICE, FILENAME = date[k]+'_'+feature[k]+'_'+exp_numb[k]+'_'+position[k]+'_day_night.ps',/COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES plot, lambda, Venus_night, xrange=[(feature_lambda[k] - 6.0), (feature_lambda[k] + 6.0)], /xstyle, yrange=[MIN(Venus_night)*1.1, MAX(Venus_night)*1.1], /ystyle, psym=0 oplot, lambda, Venus_day, psym = 0, color = 215 ; plot dayside in orange oplot, lambda, Venus_day_scaled, linestyle = 2, color = 250 ;plot scaled dayside in dashed red ;--------------------------------------------------------------------------------------- ; Plot vertical lines on where green line and venusian green line should be ; ;--------------------------------------------------------------------------------------- oplot, fltarr(2) + feature_lambda[k], !Y.CRANGE, linestyle=2,COLOR=30,THICK=2 oplot, fltarr(2) + doppler_shift[k] + feature_lambda[k], !y.crange, linestyle=2,COLOR=30,THICK=2 ; doppler_shift[k] = radial_vel[k] ;---------------------------------------------------------------------------------- ; make legend ;---------------------------------------------------------------------------------- oplot, [feature_lambda[k]+0.4, feature_lambda[k]+0.6], [Venus_night_max*1.075, Venus_night_max*1.075], linestyle = 0, color = 0 oplot, [feature_lambda[k]+0.4, feature_lambda[k]+0.6], [Venus_night_max*1.055, Venus_night_max*1.055], linestyle = 0, color = 215 oplot, [feature_lambda[k]+0.4, feature_lambda[k]+0.6], [Venus_night_max*1.035, Venus_night_max*1.035], linestyle = 2, color = 250 xyouts, feature_lambda[k]+0.7, Venus_night_max*1.07, 'Venus Night', color = 0 xyouts, feature_lambda[k]+0.7, Venus_night_max*1.05, 'Venus Day', color = 215 xyouts, feature_lambda[k]+0.7, Venus_night_max*1.03, 'Venus Day Scaled', color = 250 DEVICE,/CLOSE FREE_LUN, lun ;---------------------------------------------------------------------------------------------------------- ; 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. ;---------------------------------------------------------------------------------------------------------- linerange = WHERE((lambda GE (feature_lambda[k] - 2.0))AND(lambda LE feature_lambda[k] + 2.0)) ;Specify wavelength range we will analyze flux_max = MAX(Venus_solsub(linerange), index) flux_max_err = Venus_solsub_err(linerange(index)) lambda_max = lambda[linerange (index)] flux_min = MIN(Venus_solsub(linerange), index) flux_min_err = Venus_solsub_err(linerange(index)) lambda_min = lambda[linerange (index)] lambda = lambda[linerange] ; redefining lambda range, so lambda array is now only the range just specificed Venus_solsub = Venus_solsub[linerange] ; matching flux at new range to that point, so flux at 5577 really is at 5577 Venus_solsub_err = Venus_solsub_err[linerange] synth_solarflux_res = synth_solarflux_res[linerange] ; matching solar flux to it's proper wavelength (otherwise you'll ; plot 5577 and the first flux point, which will not be at 5577 ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************** SECTION 7 - GAUSSIAN FITTING TO DATA*************************************** ; ; This is done in units of Flux, Rayleighs, or Relative depending on what the user specifies. ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************************** 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, '' Venus_night_final = [Venus_solsub] ;make into an array with wavelength Venus_night_final_err = [Venus_solsub_err] ;------------------------------------------------------------------------------------------ ; Find error due to scatter in the point around 0. Do not consider the terrestrial or Venusian ; green line peaks, only the noise. Not sure how to add this to the true error, so for now ; I will just write it out to file. ;------------------------------------------------------------------------------------------ linerange = WHERE((lambda GE (feature_lambda[k] - 0.1))AND(lambda LE (feature_lambda[k] + 0.1))) ;set line range around terr peak terr_flux_max = MAX(Venus_night_final(linerange), index) ; find peak of terrestrial emssion terr_lambda_max = lambda[linerange(index)] ; match wavelength to peak flux linerange_venus = WHERE((lambda GE (Venus_lambda_lower)) AND (lambda LE Venus_lambda_upper)) Venus_flux_max = MAX(Venus_night_final(linerange_venus), index) ; Find the peak Venusian flux and wavelength Venus_lambda_max = lambda[linerange_venus (index)] IF Venus_lambda_guess LT feature_lambda[k] THEN BEGIN ; Set lower wavelength bounds to be Venus if lambda_lower = Venus_lambda_guess ; Venus is blue shifted of terrestrial lambda_upper = terr_lambda_max ENDIF IF Venus_lambda_guess GT feature_lambda[k] THEN BEGIN ; Set lower wavelength bounds to be terr if lambda_lower = terr_lambda_max ; Venus is red shifted of terrestrial lambda_upper = Venus_lambda_guess ENDIF lambda_continuum_lower = lambda_lower - 0.15 ; Define the area around the Venusian and terrestrial features lambda_continuum_upper = lambda_upper + 0.15 ; The continuum will be analyzed before the lower limit and above the ; upper limit, which are centered on the terrestrial/Venusian line, ; so 0.15 A is added so as to not consider the width of the feature linerange_continuum = WHERE((lambda GE (double(feature_lambda[k] - 1.0)))AND(lambda LE (feature_lambda[k] + 0.9))) ;Specify wavelength range we will analyze lambda_continuum = lambda[linerange_continuum] ; redefining lambda range, so lambda array is now only the range just specificed feature_continuum = WHERE ( (lambda_continuum LE (lambda_continuum_lower)) OR (lambda_continuum GE (lambda_continuum_upper + 0.15))) Venus_night_final_continuum = Venus_night_final[linerange_continuum] mean_continuum = MEAN((Venus_night_final_continuum(feature_continuum))) sigma_continuum = stddev(Venus_night_final_continuum(feature_continuum)) Venus_greenline_confidence = (Venus_flux_max - mean_continuum)/sigma_continuum print, '-------------------------------------------------------------------' print, feature ; print, lambda(feature_continuum) ; print, Venus_night_final(feature_continuum) print, 'Mean Continuum', mean_continuum print, 'Sigma Continuum', sigma_continuum print, '-------------------------------------------------------------------' ;------------------------------------------------------------------------------------------ ; Print information to file ;------------------------------------------------------------------------------------------ printf, lun_1, venus_night_name[k], position[k], feature[k], 'Flux', lambda[0], lambda_continuum_lower, $ lambda_continuum_upper, lambda[N_elements(lambda)-1], mean_continuum, $ sigma_continuum, terr_flux_max, Venus_flux_max, Venus_greenline_confidence, $ FORMAT='(A,5X, A, 5X, A, 5X, 15A, 5X, F,F, X, e5.2, e5.2, 10X, e5.2, X, e5.2, X, f5.2)' ;------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------ ; 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[2].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 ; This can be done by inserting the actual value, or by a variable. ; Here, we are going to calculate the shift using the radial vel ; of Venus. ;------------------------------------------------------------------------------------------ parinfo = replicate({fixed:0, limited:[1,1], limits:[1e-20,1e6]},5) parinfo[0].limits = [1e-20,1e-8] ; Venusian Flux parinfo[1].limits = [DOUBLE(Venus_lambda_lower),DOUBLE(Venus_lambda_upper)] ; Venusian wavelength parinfo[2].LIMITS = [dispersion[k]-0.01, dispersion[k]+0.01] ; Dispersion of echelle parinfo[3].limits = [1e-20,1e-8] ; Terrestrial Flux ; parinfo[3].limits = [DOUBLE(terr_flux_max*0.8), DOUBLE(terr_flux_max*1.2)] ; Terrestrial Flux parinfo[4].limits = [DOUBLE(feature_lambda[k] - 0.15), DOUBLE(feature_lambda[k] + 0.15)] ; Terrestrial wavelength ; parinfo[0].value = DOUBLE(gaussian_flux_venus[k]) ; Venusian Flux ; parinfo[1].value = 5577.55D ; Venusian wavelength ; parinfo[2].value = 0.79 ; Dispersion of echelle ; parinfo[3].value = 1e-14 ; Terrestrial Flux ; parinfo[4].value = 5577.3D ; Terrestrial wavelength start = [DOUBLE(gaussian_flux_venus[k]), DOUBLE(Venus_lambda_guess), DOUBLE(dispersion[k]), DOUBLE(gaussian_flux_terr[k]), DOUBLE(feature_lambda[k])] ;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(Venus_night_final),1) ;this used to have 2048 instead of N_ELEMENTS fluxfittell = fltarr(N_ELEMENTS(Venus_night_final),1) residual = fltarr(N_ELEMENTS(Venus_night_final),1) totresidual = fltarr(1) chisquared = fltarr(1) redchi = fltarr(1) averesid = fltarr(1) ew5577 = fltarr(1) ew5577_terr = 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=3 !P.POSITION=[.19,.12,.85,.95] ;set lower left and upper right corners of plotting box !X.CHARSIZE = 1.45 !Y.CHARSIZE = 1.45 !X.THICK=4 !Y.THICK=4 SET_PLOT,'PS' DEVICE, FILENAME = date[k]+'_'+feature[k]+'_'+exp_numb[k]+'_'+position[k]+'_flux_FINAL.ps',/COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES ;------------------------------------------------------------------------------------------------------------- ;--------------------------------------- BEGIN GAUSSIAN FIT FOR FLUX ----------------------------------------- ;------------------------------------------------------------------------------------------------------------- ;------------------------------------------------------------------------------------------------------------- ; The fittin program, MPFITFUN, will fit a gaussian to the data given the wavelength, flux, and error to the ; data. It does consider the error in the data when attempting to make the best fit. Thus, if you have poor ; quatlity data, and the majority of the error is in the transmission, then thepoints with the most flux will ; have the largest error. Therefore, the terrestrial and Venusian green line flux peaks will have the largest ; error compared to the surround flux points. Given MPFITFUN will consider this error, it will likely fit a ; gaussian function to the wrong data point, choosing one in the continuum because of it's small error bar. To ; avoid this problem, the error for all points is set to 1 and then scaled down to a very small value, 1e-18. ; This makes all points equally weighted and the green line peaks can then be found with the function. ; The origianl fitting function read as: ; P[*,i] = MPFITFUN('gfunct', lambda, Venus_night_final[*,i], Venus_night_final_err[*,i], start, PARINFO=parinfo, ; NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1) ; Venus_night_final_err[*,i] has been modified to Venus_night_final_err[*,i]/ Venus_night_final_err[*,i] * 1e-18 ; to make all errors equal and small. This is NOT be the error plotted. The true error will be plotted, this ; adjustment is made only for Gaussian fitting purposes. ;------------------------------------------------------------------------------------------------------------- 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, Venus_night_final[*,i], Venus_night_final_err[*,i]/ Venus_night_final_err[*,i] * 1e-18, $ start, PARINFO=parinfo, NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN, quiet=1) ; P[*,i] = MPFITFUN('gfunct', lambda, Venus_night_final[*,i], Venus_night_final_err[*,i], start, PARINFO=parinfo, NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN) ;, quiet=1) plot, lambda, Venus_night_final[*,i], xrange=[feature_lambda[k] - plot_lower_bound, feature_lambda[k] + plot_upper_bound], /xstyle, $ yrange=[(flux_min - flux_min_err*1.1),(flux_max + flux_max_err*1.1)], /ystyle, psym=4, $ xtitle='Wavelength ('+angstrom+')', ytitle='Flux (ergs cm!u-2 !n s!u-1 !n '+angstrom+'!u-1!n)', $ title = date[k]+' - '+position[k]+' '+exp_numb[k] +' ' + feature[k], THICK=8 ; charsize=1.5 ;--------------------------------------------------------------------------------------- ; Plot error on each point. ; If the errors are off the scale, multiple by a small number ; just to compare relative sizes of errors on each point. ;-------------------------------------------------------------------------------------- oploterror, lambda, Venus_night_final[*,i], Venus_night_final_err[*,i], psym=4, thick = 4 ;plots the error. Thick for connecting lines ;--------------------------------------------------------------------------------------- ; Plot the synthetic solar spectrum. Comment this out for the final plot. ;--------------------------------------------------------------------------------------- ; oplot,lambda, synth_solarflux_res_scaled, psym = 0 ;is not plotting correctly ;--------------------------------------------------------------------------------------- ; Plots Guassian of terrestrail green line: ; P[0,i] * exp(-(wavelength-P[1,i])^2.0 / (2.0 * (P[2,i]^2.0)) ; ; and the Doppler shifted green line: ; P[3,i] * exp(-(wavelength-P[4,i])^2.0 / (2.0 * (P[5,i]^2.0))) ; ; added together. It's drawn as a dashed line. If it's ; not matched in flux to the observed emission, something went ; wrong somewhere. ;----------------------------------------------------------------------------------------- wavelength = findgen(3000)/1000 + feature_lambda[k] - 1.5 ; Where the x range will start 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[2,i]^2.0))))), $ linestyle=2, COLOR=82,THICK=8 ;--------------------------------------------------------------------------------------- ; Guassian plot of just the terrestral green line, so you can check ; how it compares to the observed line. If this line fits better ; than the blended line, then you don't have green line emission. ;--------------------------------------------------------------------------------------- oplot, wavelength, (P[3]*exp(-(wavelength-P[4])^2.0/(2.0*(P[2]^2.0)))), $ linestyle=4,COLOR=250,THICK=8 ;Terrestrail Gaussian only oplot, wavelength, (P[0]*exp(-(wavelength-P[1])^2.0/(2.0*(P[2]^2.0)))), $ linestyle=4,COLOR=150,THICK=8 ;Terrestrail Gaussian only print, 'parameters:', P print, 'parameter errors:', perror ;--------------------------------------------------------------------------------------- ; Plot the upper and lower bounds to a 1 sigma stand deviation around the mean ; of the scatter of the continuum ;--------------------------------------------------------------------------------------- sigma_continuum_upper = mean_continuum + sigma_continuum sigma_continuum_lower = mean_continuum - sigma_continuum oplot, [5575, 5579], [sigma_continuum_upper, sigma_continuum_upper], THICK = 4, LINESTYLE = 1, COLOR = 10 oplot, [5575, 5579], [sigma_continuum_lower, sigma_continuum_lower], THICK = 4, LINESTYLE = 1, COLOR = 10 oplot, [5575, 5579], [mean_continuum, mean_continuum], 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[*,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] - (Venus_night_final[*,i])) / (Venus_night_final_err[*,i]) fiterr=replicate(0.0,N_ELEMENTS(Venus_night_final)) ;!!! chisquared[i] = 0.0 totresidual[i] = 0.0 for j=0,N_ELEMENTS(Venus_night_final)-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 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/N_ELEMENTS(Venus_night_final) ;divided by the number of points you're sampling over (N_elements..) ; print, 'average residual:', averesid[i] ; print, 'reduced chi-squared:',redchi[i] ; print, 'integrated flux:', ew5577[i] ; print, 'error in flux:', ewerr5577[i] endfor ;------------------------------------------------------------------------------------------ ; Print data to file ;------------------------------------------------------------------------------------------ openw, lun, 'output_'+string(exp_numb[k])+'_'+string(feature[k])+'_flux.dat', /GET_LUN ;**** UPDATE VARIABLE **** printf, lun, '#---------------------------------------------------' printf, lun, '# Wavelength Flux Error' printf, lun, '# (A) (cgs)' printf, lun, '#---------------------------------------------------' for n=0, N_ELEMENTS(lambda)-1 do begin printf, lun, lambda[n], Venus_night_final[n,0], Venus_night_final_err[n,0] endfor FREE_LUN, lun ;------------------------------------------------------------------------------------------ ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************************** RAYLEIGH UNITS ********************************************* ; ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; Endif Else If(tmpread EQ '2')Then Begin PRINT,'' PRINT, '' print, 'The following calculations will be done in Rayleighs units' PRINT,'' PRINT, '' ;--------------------------------------------------------------------- ; There are two conversions to Rayleighs that must be made: one for ; Earth, where we are inside the shell of the emitting layer, and one ; for Venus, where we are observing outside the shell. I think the best ; way to do this is to first convert to Rayleighs and then to the two ; zenith angles. The order shouldn't matter because they are both ; values multipled to the observed flux. ; ; When plotting the data, we have to plot one of the rayleigh values. ; I think it would be easiest if we plotted the Venusian flux in ; Rayleighs, and then convert that to the Terrestrial by undoing ; the sec(@) correction and doing the terrestrial correction. ;--------------------------------------------------------------------- ;--------------------------------------------------------------------- ; Convert flux to Rayleighs on Venus ;--------------------------------------------------------------------- ; ; In order to convert the observed Venusian flux to Rayleighs, the ; local zenith angle on Venus must be considered. As Rayleighs are ; defined at a zenith angle of 0 (the sub-Earth point, or nadir), the ; emission must be corrected for the local Venus zenith angle that we ; are looking at. ; ; We will assume a plane-parallel, infinitally long emitting ; layer (see Crisp et al, JGR 1996 for a full description). This ; assumption does not work for local zenith angles greater than 85 ; (i.e., the limb). However, our slit is so large that the angle ; between 85 - 90 will only be a small portion of the slit, less than ; half. ; ; Here, we calculate what the local Venus zenith angle is at the ; center of the slit and will use that as our emission angle to be ; converted to emission at angle 0. ; The position of the slit is found by looking at an ecam image of ; different contrasts, to see the edge of Venus is one and the ; slit in the other. A circle is fit to Venus and the center of ; the disk marked, then the position of the slit from the disk ; center is measured. You have to do all this outside the code. ; The disk radius is normalized to 1 and the distance to the slit found ; by dist/radius. The Venus zenith angle at the center point of the ; slit is found by taking the arc sin (sin-1) of the distance. This ; will be given in radians. The Venus flux is then corrected to the ; Venus zenith by flux*cos(zenith_angle) ;---------------------------------------------------------------------- slit_dist[k] = slit_center_distance[k]/disk_radius[k] ;normalizing Venus_zenith_angle[k] = asin(slit_dist[k]) Venus_zenith_flux = Venus_solsub * cos(Venus_zenith_angle[k]) ; F(zenith) = F(obs)/sec(@) = F(obs)*cos(a) ; where we have to convert @ back to radians for idl ; We never had to put it in degress, I just like degrees. ;------------------------------------------------------------------------- ; Convert Venus zenith flux to Rayleighs ;------------------------------------------------------------------------- conv_const = 3.71546e-14 ; h and some other stuff peak = feature_lambda[k] ; wavelength slit_area = 5.12 ; 1.6*3.2 slit size Rayleigh = Venus_zenith_flux * peak / (slit_area*conv_const) Venus_solsub = Rayleigh Venus_solsub_err = Venus_solsub_err * (peak / (slit_area * conv_const)) * cos(Venus_zenith_angle[k]) ;------------------------------------------------------------------------- ; Convert flux to Rayleighs on Earth ;------------------------------------------------------------------------- ; ; As Rayleighs are defined at zenith, and nightglow emission increases ; intensity with airmass, you need to convert the intensity seen at a ; given angle to what it would be at zenith. The correction is: ; ; F(obs) = F(zenith)*sec(theta) ; ; Where F(obs) is the observed Flux, F(zenith) is the flux you would ; see at zenith, and theta is the angle mesasured from the observing ; layer to the observer. In order to find the angle theta, the ; van Rhijn method is used (van Rhijn 1921a, see Chamberlain's Physics ; of Aurora and nightglow book). The van Rhijn method is: ; ; The Rhihn method give the relation between theta and the observed ; zenith angle, alpha, where 'a' is the radius of Earth and z is the ; altitdue of the emitting layer. ; ; sec(theta) = 1 / [ 1 - (a/(a+z))^2 * sin^2(alpha) ]^(1/2) ; ; The y-axis on the right side is then adjusted by this amount ; ; First, the flux must be scaled to the flux at zenith using the ; Rhijn method. Next, the zenith flux is converted to Raylieghs. ; ;------------------------------------------------------------------------- ; conv_const = 3.71546e-14 ; peak = feature_lambda[k] ; airmass = venus_night_airmass[k] ; slit_area = 5.12 ; Rayleigh = Venus_solsub*peak/(slit_area*conv_const) ; Venus_solsub = Rayleigh ; Venus_solsub_err = Venus_solsub_err * peak / (slit_area * conv_const) ;------------------------------------------------------------------------ ; Find max and min Rayleigh value and errors ;------------------------------------------------------------------------ rayleigh_max = MAX(Venus_solsub,index) rayleigh_max_err = Venus_solsub_err(index) lambda_max = lambda[index] rayleigh_min = MIN(Venus_solsub, index) rayleigh_min_err = Venus_solsub_err(index) lambda_min = lambda[index] ;------------------------------------------------------------------------ ; Make flux and error values into an array ;------------------------------------------------------------------------ Venus_night_final = [Venus_solsub] ;make into an array with wavelength Venus_night_final_err = [Venus_solsub_err] ;------------------------------------------------------------------------ ; Calculate error of mean around continuum and print to file ;------------------------------------------------------------------------- linerange = WHERE((lambda GE (feature_lambda[k] - 0.1))AND(lambda LE (feature_lambda[k] + 0.1))) ; set line range around terr peak terr_flux_max = MAX(Venus_night_final(linerange), index) ; find peak of terrestrial emssion terr_lambda_max = lambda[linerange(index)] ; match wavelength to peak flux linerange_venus = WHERE((lambda GE (Venus_lambda_lower)) AND (lambda LE Venus_lambda_upper)) Venus_flux_max = MAX(Venus_night_final(linerange_venus), index) ; Find the peak Venusian flux and wavelength Venus_lambda_max = lambda[linerange_venus (index)] IF Venus_lambda_guess LT feature_lambda[k] THEN BEGIN ; Set lower wavelength bounds to be Venus if lambda_lower = Venus_lambda_guess ; Venus is blue shifted of terrestrial lambda_upper = terr_lambda_max ENDIF IF Venus_lambda_guess GT feature_lambda[k] THEN BEGIN ; Set lower wavelength bounds to be terr if lambda_lower = 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 GE (double(feature_lambda[k] - 1.0)))AND(lambda LE (feature_lambda[k] + 0.9))) ;Specify wavelength range we will analyze lambda_continuum = lambda[linerange_continuum] ; redefining lambda range, so lambda array is now only the range just specificed feature_continuum = WHERE ( (lambda_continuum LE (lambda_continuum_lower)) OR (lambda_continuum GE (lambda_continuum_upper + 0.15))) Venus_night_final_continuum = Venus_night_final[linerange_continuum] mean_continuum = MEAN((Venus_night_final_continuum(feature_continuum))) sigma_continuum = stddev(Venus_night_final_continuum(feature_continuum)) Venus_greenline_confidence = (Venus_flux_max - mean_continuum)/sigma_continuum print, '-------------------------------------------------------------------' ; print, lambda(feature_continuum) ; print, Venus_night_final(feature_continuum) print, 'Mean Continuum', mean_continuum print, 'Sigma Continuum', sigma_continuum print, '-------------------------------------------------------------------' ; linerange_continuum = WHERE ( (lambda LE (lambda_continuum_lower)) OR (lambda GE (lambda_continuum_upper + 0.15))) ; mean_continuum = MEAN((Venus_night_final(linerange_continuum))) ; sigma_continuum = stddev(Venus_night_final(linerange_continuum)) ; Venus_greenline_confidence = (Venus_flux_max - mean_continuum)/sigma_continuum ;------------------------------------------------------------------------------------------ ; Print information to file ;------------------------------------------------------------------------------------------ printf, lun_1, venus_night_name[k], position[k], feature[k], 'Flux', lambda[0], lambda_continuum_lower, $ lambda_continuum_upper, lambda[N_elements(lambda)-1], mean_continuum, $ sigma_continuum, terr_flux_max, Venus_flux_max, Venus_greenline_confidence, $ FORMAT='(A,5X, A, 5X, A, 15A, 5X, F,F, X, e5.2, e5.2, 10X, e5.2, X, e5.2, X, f5.2)' ;------------------------------------------------------------------------------------------ ; 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 = [1.0, 5.0e4] ; Venusian Rayleigh parinfo[1].limits = [DOUBLE(Venus_lambda_lower), DOUBLE(Venus_lambda_upper)] ; Venusian wavelength parinfo[2].limits = [dispersion[k]-0.01, dispersion[k]+0.01] ; Dispersion of echelle parinfo[3].limits = [1.0, 5.0e4] ; Terrestrial Rayleigh parinfo[4].limits = [DOUBLE(feature_lambda[k] - 0.1), DOUBLE(feature_lambda[k] + 0.1)] ; Terrestrial wavelength start = [DOUBLE(gaussian_ray_venus[k]), DOUBLE(Venus_lambda_guess), DOUBLE(dispersion[k]), DOUBLE(gaussian_ray_terr[k]), DOUBLE(feature_lambda[k])] ;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(Venus_night_final),1) ;this used to have 2048 instead of N_ELEMENTS fluxfittell = fltarr(N_ELEMENTS(Venus_night_final),1) residual = fltarr(N_ELEMENTS(Venus_night_final),1) totresidual = fltarr(1) chisquared = fltarr(1) redchi = fltarr(1) averesid = fltarr(1) ew5577 = fltarr(1) ew5577_terr = fltarr(1) ewerr5577 = fltarr(1) ewerr5577_terr = fltarr(1) ;------------------------------------------------------------------------------------------------------------ ; --------------------------------------- PLOTTING IN RAYLEIGH UNITS ---------------------------------------- ;------------------------------------------------------------------------------------------------------------ angstrom = '!6!sA!r!u!9 %!6!n' ;define the angstrom character print, 'begin loop' LOADCT,39,/SILENT !P.CHARTHICK=3 !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 = date[k]+'_'+feature[k]+'_'+exp_numb[k]+'_'+position[k]+'_rayleigh_FINAL.ps',/COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES ;------------------------------------------------------------------------------------------------------------- ;--------------------------------------- BEGIN GAUSSIAN FIT FOR RAYLEIGHS ----------------------------------------- ;------------------------------------------------------------------------------------------------------------- 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, Venus_night_final[*,i], Venus_night_final_err[*,i]/ Venus_night_final_err[*,i] * 1e-18, start, PARINFO=parinfo, $ NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN, quiet=1) plot, lambda, Venus_night_final[*,i], linestyle=0, xrange=[(feature_lambda[k] - plot_lower_bound), (feature_lambda[k] + plot_upper_bound)], /xstyle,$ yrange=[(rayleigh_min - rayleigh_min_err*1.1),(rayleigh_max + rayleigh_max_err*1.1)], ystyle=9, psym=4, $ xtitle='Wavelength ('+angstrom+' )', ytitle='Venusian Intensity (R/'+angstrom+')', $ title = date[k]+' - '+position[k]+' '+exp_numb[k]+' '+feature[k]+'', charsize=1.5, THICK=8 ;----------------------------------------------------------------------------------------------------- ; Set up different range on right axis for the Terrestrial intensity. ; First, the value plotted, the Venus flux converted to the zenith value ; and then to Rayleighs, must be converted back to normal flux and rayleigh, ; so undo the correction to Venus zenith. ; ; F(obs) = F(zenith_venus) * sec(@) ; F(zenith_venus) = F(obs) / sec(@) ; = F(obs) * cos(@) ; ; Venus_rayleigh = F(zenith_venus) * Rayleigh conversion ; ; The Rayleigh conversion is already done, so undo the cos(@) and then multiply by ; cos(theta), where (theta) is the angle of the emitting layer on Earth to the observer. ; The angle & is find by using the van Rhijn method: which gives the relation between ; the solar zenith angle on earth and the emitting layer angle. ; ; Terr_rayleigh = (Venus_rayleigh / cos(@)) / sec(theta) ; ; = Venus_rayleigh / (cos(@) * sec(theta)) ; sec(theta) = 1 / [ 1 - (a/(a+z))^2 * sin^2(alpha) ]^(1/2) ; ; a = radius of Earth in km ; z = emitting altitude in km ; alpha = observed angle measured from zenith ; ; Divide by sec(theta) to convert observerd intensity to zenith ; ; To get the axis of the plot to go down and not up, you need to do a XYOutS command and change ; the orientation by 90 degrees the opposite way. ;----------------------------------------------------------------------------------------------------- a = 6371.0 ; radius of earth in KM z = 95.0 ; Earth's nightglow emitting layer in km (is this true for N2??) alpha[k] = Venus_night_angle[k] * 2.0*3.15159 / 360.0 ; Venus_night_angle is the Earth's solar zenith angle in degrees, ; which is changed to radians since idl only uses radians sec_theta[k] = 1.0 / [ 1.0 - (a/(a+z))^2.0 * (sin(alpha[k]))^2.0 ]^(0.5) axis, yaxis=1, yrange=[(rayleigh_min - rayleigh_min_err*1.1)/ (sec_theta[k] * cos(Venus_zenith_angle[k])), $ (rayleigh_max + rayleigh_max_err*1.1)/ (sec_theta[k] * cos(Venus_zenith_angle[k]))], $ ;y axis values ystyle=1, charsize=1.5 xloc = !X.Window[1] + 0.102 ;define x position for label yloc = (0.9 - 0.15) / 2.0 + 0.15 ;define y position for label XYOutS, xloc, yloc, 'Terrestrial Intensity (R/' +angstrom+')!C', 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. ; NOTE!!!! When you plot the error, it replots the lines ; connecting the points, so you have to chnage the line ; thickness here, not in the initial plot sections ;-------------------------------------------------------------------------------------- oploterror, lambda, Venus_night_final[*,i], Venus_night_final_err[*,i], psym = 4, thick=4 ;plots the error. Thick for connecting lines ;--------------------------------------------------------------------------------------- ; Plot the synthetic solar spectrum. Comment this out for the final plot. ;--------------------------------------------------------------------------------------- ; oplot,lambda, synth_solarflux_res_scaled, psym = 0 ;--------------------------------------------------------------------------------------- ; Plots Guassian of terrestrail green line: ; P[0,i] * exp(-(wavelength-P[1,i])^2.0 / (2.0 * (P[2,i]^2.0)) ; ; and the Doppler shifted green line: ; P[3,i] * exp(-(wavelength-P[4,i])^2.0 / (2.0 * (P[5,i]^2.0))) ; ; added together. It's drawn as a dashed line. If it's ; not matched in flux to the observed emission, something went ; wrong somewhere. ;----------------------------------------------------------------------------------------- wavelength = findgen(3000)/1000 + feature_lambda[k] - 1.5 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[2,i]^2.0))))), $ linestyle=2, COLOR=82,THICK=8 ;--------------------------------------------------------------------------------------- ; Guassian plot of just the terrestral green line, so you can check ; how it compares to the observed line. If this line fits better ; than the blended line, then you don't have green line emission. ;--------------------------------------------------------------------------------------- oplot, wavelength, (P[3]*exp(-(wavelength-P[4])^2.0/(2.0*(P[2]^2.0)))), $ linestyle=4,COLOR=250,THICK=8 ;Terrestrail Gaussian only 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 ;--------------------------------------------------------------------------------------- sigma_continuum_upper = mean_continuum + sigma_continuum sigma_continuum_lower = mean_continuum - sigma_continuum oplot, [5575, 5579], [sigma_continuum_upper, sigma_continuum_upper], THICK = 4, LINESTYLE = 1, COLOR = 10 oplot, [5575, 5579], [sigma_continuum_lower, sigma_continuum_lower], THICK = 4, LINESTYLE = 1, COLOR = 10 oplot, [5575, 5579], [mean_continuum, mean_continuum], 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[*,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] - Venus_night_final[*,i]) / (Venus_night_final_err[*,i]) fiterr=replicate(0.0,N_ELEMENTS(Venus_night_final)) ;!!! chisquared[i] = 0.0 totresidual[i] = 0.0 for j=0,N_ELEMENTS(Venus_night_final)-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 redchi[i] = (chisquared[i])/4.0 ew5577[i] = 2.50663 * P[0,i] * P[2,i] ew5577_terr[i] = 2.50663 * P[3,i] * P[2,i] / (sec_theta[k]* cos(Venus_zenith_angle[k])) 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 if sqrt(redchi[i]) gt 1.0 then begin ewerr5577_terr[i] = 2.50663 * perror[0] * P[2,i] * sqrt(redchi[i]) endif else begin ewerr5577_terr[i] = 2.50663 * perror[0] * P[2,i] endelse averesid=totresidual/N_ELEMENTS(Venus_night_final) ;divided by the number of points you're sampling over (N_elements..) ; print, 'average residual:', averesid[i] ; print, 'reduced chi-squared:', redchi[i] ; print, 'integrated flux Venus:', ew5577[i] ; print, 'integrated flux Earth:', ew5577_terr[i] ; print, 'error in flux:', ewerr5577[i] ;------------------------------------------------------------------------------------------ ; Print integrated intensity to file ;------------------------------------------------------------------------------------------ printf, lun_2, venus_night_name[k], position[k], feature[k], 'Rayleigh', ew5577_terr[i], ewerr5577_terr[i], ew5577[i], ewerr5577[i],$ FORMAT='(A,5X, A, 5X, A, 5x, 15A, 5X, F,F, X, F,F)' ;------------------------------------------------------------------------------------------ DEVICE,/CLOSE endfor ;------------------------------------------------------------------------------------------ ; Print data to file ;------------------------------------------------------------------------------------------ openw, lun, 'output_'+string(exp_numb[k])+'_'+feature[k]+'_rayleigh.dat', /GET_LUN ;**** UPDATE VARIABLE **** printf, lun, '#---------------------------------------------------' printf, lun, '# Wavelength Rayleighs Error' printf, lun, '# (A) ' printf, lun, '#---------------------------------------------------' for n=0, N_ELEMENTS(lambda)-1 do begin printf, lun, lambda[n], Venus_night_final[n,0], Venus_night_final_err[n,0] endfor FREE_LUN, lun DEVICE,/CLOSE ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; ; ****************************************** NORMALIZED UNITS ******************************************* ; ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ; Endif Else If(tmpread EQ '3')Then Begin PRINT,'' PRINT, '' print, 'The following calculations will be done in Relative units, scaled to the terrestrial green line' PRINT,'' PRINT, '' ;---------------------------------------------------------------------------------------------------- ; For Normalized units, divide flux by the Max value relative units ; But if the max value is the venusian peak, then this will scale to Venus ; and not Earth's line. To fix, find the max in the range around the terrestrial ;---------------------------------------------------------------------------------------------------- linerange = WHERE((lambda GE (Venus_lambda_lower))AND(lambda LE Venus_lambda_upper)) Venus_flux_max = MAX(Venus_solsub(linerange), index) Venus_lambda_max = lambda[linerange (index)] linerange = WHERE((lambda GE (feature_lambda[k] - 0.1))AND(lambda LE (feature_lambda[k] + 0.1))) terr_flux_max = MAX(Venus_solsub(linerange), index) terr_lambda_max = lambda[linerange (index)] Venus_solsub_scale = Venus_solsub/terr_flux_max Venus_night_final = [Venus_solsub_scale] ;-------------------------------------------------------- ; want min and max value of the whole set of points to set ; the y range for plotting ;------------------------------------------------------- flux_max = MAX(Venus_solsub_scale,index) flux_max_err = (Venus_solsub_err(index)) * Venus_solsub_scale(index) lambda_max = lambda[index] flux_min = MIN(Venus_solsub_scale, index) flux_min_err = (Venus_solsub_err(index)) * Venus_solsub_scale(index) lambda_min = lambda[index] ;------------------------------------------------------------ ; Calculate error of mean around continuum and print to file ;------------------------------------------------------------ linerange = WHERE((lambda GE (feature_lambda[k] - 0.1))AND(lambda LE (feature_lambda[k] + 0.1))) ;set line range around terr peak terr_flux_max = MAX(Venus_night_final(linerange), index) ; find peak of terrestrial emssion terr_lambda_max = lambda[linerange(index)] ; match wavelength to peak flux linerange_venus = WHERE((lambda GE (Venus_lambda_lower)) AND (lambda LE Venus_lambda_upper)) Venus_flux_max = MAX(Venus_night_final(linerange_venus), index) ; Find the peak Venusian flux and wavelength Venus_lambda_max = lambda[linerange_venus (index)] IF Venus_lambda_guess LT feature_lambda[k] THEN BEGIN ; Set lower wavelength bounds to be Venus if lambda_lower = Venus_lambda_guess ; Venus is blue shifted of terrestrial lambda_upper = terr_lambda_max ENDIF IF Venus_lambda_guess GT feature_lambda[k] THEN BEGIN ; Set lower wavelength bounds to be terr if lambda_lower = 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 GE (double(feature_lambda[k] - 1.0)))AND(lambda LE (feature_lambda[k] + 0.9))) ;Specify wavelength range we will analyze lambda_continuum = lambda[linerange_continuum] ; redefining lambda range, so lambda array is now only the range just specificed feature_continuum = WHERE ( (lambda_continuum LE (lambda_continuum_lower)) OR (lambda_continuum GE (lambda_continuum_upper + 0.15))) Venus_night_final_continuum = Venus_night_final[linerange_continuum] mean_continuum = MEAN((Venus_night_final_continuum(feature_continuum))) sigma_continuum = stddev(Venus_night_final_continuum(feature_continuum)) Venus_greenline_confidence = (Venus_flux_max - mean_continuum)/sigma_continuum print, '-------------------------------------------------------------------' ; print, lambda(feature_continuum) ; print, Venus_night_final(feature_continuum) print, 'Mean Continuum', mean_continuum print, 'Sigma Continuum', sigma_continuum print, '-------------------------------------------------------------------' ;------------------------------------------------------------------------------------------ ; Print information to file ;------------------------------------------------------------------------------------------ printf, lun_1, venus_night_name[k], position[k], feature[k], 'Flux', lambda[0], lambda_continuum_lower, $ lambda_continuum_upper, lambda[N_elements(lambda)-1], mean_continuum, $ sigma_continuum, terr_flux_max, Venus_flux_max, Venus_greenline_confidence, $ FORMAT='(A,5X, A, 5X, A, 5x, 15A, 5X, F,F, X, e5.2, e5.2, 10X, e5.2, X, e5.2, X, f5.2)' ; 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))) ; mean_continuum = MEAN((Venus_night_final(linerange_continuum))) ; sigma_continuum = stddev(Venus_night_final(linerange_continuum)) ; Venus_greenline_confidence = (Venus_flux_max/terr_flux_max - mean_continuum)/sigma_continuum ; printf, lun_1, venus_night_name[k], position[k], 'Norm', lambda[0], lambda_continuum_lower, $ ; lambda_continuum_upper, lambda[N_elements(lambda)-1], mean_continuum, $ ; sigma_continuum, terr_flux_max, Venus_flux_max, Venus_greenline_confidence, $ ; FORMAT='(A,5X,A, 5X, 15A, 5X, F,F, X, F,F, 10X, F, X, e2.2, X, F)' ;------------------------------------------------------------ ; 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, 4] ; Venusian Relative parinfo[1].limits = [Venus_lambda_lower, Venus_lambda_upper] parinfo[2].LIMITS = [dispersion[k]-0.01, dispersion[k]+0.01] parinfo[3].limits = [0, 1] ; Terrestrial Relative parinfo[4].limits = [DOUBLE(feature_lambda[k] - 0.1), DOUBLE(feature_lambda[k] + 0.1)] start = [DOUBLE(gaussian_norm_venus[k]), DOUBLE(Venus_lambda_guess), double(dispersion[k]), DOUBLE(gaussian_norm_terr[k]), DOUBLE(feature_lambda[k])] 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(Venus_night_final),1) ;this used to have 2048 instead of N_ELEMENTS fluxfittell = fltarr(N_ELEMENTS(Venus_night_final),1) residual = fltarr(N_ELEMENTS(Venus_night_final),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 LOADCT,39,/SILENT !P.CHARTHICK=3 !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 = date[k]+'_'+feature[k]+'_'+exp_numb[k]+'_'+position[k]+'_normalized_FINAL_publish.ps',/COLOR,/LANDSCAPE,YSIZE=7,XSIZE=9.5,/COURIER,/INCHES ;------------------------------------------------------------------------------------------------------------- ;--------------------------------------- BEGIN GAUSSIAN FIT FOR NORMALIZED------------------------------------ ;------------------------------------------------------------------------------------------------------------- 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, Venus_night_final[*,i], Venus_night_final_err[*,i]/ Venus_night_final_err[*,i] * 1e-18, start, PARINFO=parinfo, $ NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, /NAN, quiet=1) plot, lambda, Venus_night_final[*,i], xrange=[5576.6,5578.2], /xstyle, $ yrange=[-0.4,1.2], $ /ystyle, psym=10, xtitle='Wavelength ('+angstrom+')', $ ytitle='Normalized Units', title=date[k]+' - '+position[k]+' '+exp_numb[k]+''+feature[k]+'', $ charsize=1.5, THICK=8 ; xrange=[feature_lambda[k] - plot_lower_bound, feature_lambda[k] + plot_upper_bound] ;--------------------------------------------------------------------------------------- ; Plot the synthetic solar spectrum. Comment this out for the final plot. ;--------------------------------------------------------------------------------------- ; oplot,lambda, synth_solarflux_res_scaled, psym = 0 ;--------------------------------------------------------------------------------------- ; 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,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))))), $ linestyle=2, COLOR=82, THICK=8 ;--------------------------------------------------------------------------------------- ; Guassian plot of just the terrestral green line, so you can check ; how it compares to the observed line. If this line fits better ; than the blended line, then you don't have green line emission. ;--------------------------------------------------------------------------------------- oplot, wavelength, (P[3]*exp(-(wavelength-P[4])^2.0/(2.0*(P[2]^2.0)))), $ linestyle=4,COLOR=250,THICK=8 ;Terrestrail Gaussian only oplot, wavelength, (P[0]*exp(-(wavelength-P[1])^2.0/(2.0*(P[2]^2.0)))), $ linestyle=4,COLOR=206,THICK=8 ;Venusain Gaussian only ;--------------------------------------------------------------------------------------- ; Plot vertical lines on where green line and venusian green line should be ; ;--------------------------------------------------------------------------------------- oplot, fltarr(2) + feature_lambda[k], !Y.CRANGE, linestyle=2,COLOR=30,THICK=2 oplot, fltarr(2) + doppler_shift[k] + feature_lambda[k], !y.crange, linestyle=2,COLOR=30,THICK=2 ; doppler_shift[k] = radial_vel[k] 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 ;--------------------------------------------------------------------------------------- sigma_continuum_upper = mean_continuum + sigma_continuum sigma_continuum_lower = mean_continuum - sigma_continuum oplot, [5575, 5579], [sigma_continuum_upper, sigma_continuum_upper], THICK = 4, LINESTYLE = 1, COLOR = 10 oplot, [5575, 5579], [sigma_continuum_lower, sigma_continuum_lower], THICK = 4, LINESTYLE = 1, COLOR = 10 oplot, [5575, 5579], [mean_continuum, mean_continuum], THICK = 4, LINESTYLE = 0, COLOR = 10 print, '-------------------------------------------------------------------' print, feature_lambda[k] ; print, lambda(feature_continuum) print, Venus_night_final(feature_continuum) print, 'Mean Continuum', mean_continuum print, 'Sigma Continuum', sigma_continuum print, '-------------------------------------------------------------------' ;------------------------------------------------------------------------------------------ ; 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] - (Venus_night_final[*,i])) / (Venus_night_final_err[*,i]) fiterr=replicate(0.0,N_ELEMENTS(Venus_night_final)) ;!!! chisquared[i] = 0.0 totresidual[i] = 0.0 for j=0,N_ELEMENTS(Venus_night_final)-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 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/N_ELEMENTS(Venus_night_final) ;divided by the number of points you're sampling over (N_elements..) ; print, 'average residual:', averesid[i] ; print, 'reduced chi-squared:',redchi[i] ; print, 'integrated flux:', ew5577[i] ; print, 'error in flux:', ewerr5577[i] endfor ;------------------------------------------------------------------------------------------ ; Print data to file ;------------------------------------------------------------------------------------------ openw, lun, 'output_'+string(exp_numb[k])+'_'+feature[k]+'_normalized.dat', /GET_LUN ;**** UPDATE VARIABLE **** printf, lun, '#---------------------------------------------------' printf, lun, '# Wavelength Normalized Error' printf, lun, '# (A) Units ' printf, lun, '#---------------------------------------------------' for n=0, N_ELEMENTS(lambda)-1 do begin printf, lun, lambda[n], Venus_night_final[n,0], Venus_night_final_err[n,0] endfor FREE_LUN, lun ;------------------------------------------------------------------------------------------ DEVICE,/CLOSE ; Endif Else Print,'' ; Print,'WRONG. End program' Print,'' Print,'' ENDFOR FREE_LUN, lun_1 FREE_LUN, lun_2 STOP end