;+ ; Project : Course Computer Project ; ; Name : HUBBLE ; ; Purpose : Driver procedure to determine Hubble's constant from spectra. ; ; Explanation : This is the main driving procedure to figure out Hubble's ; constant from the galaxy spectra of McQuade, Calzetti, and ; Kinney 1995, ApJS, 97, 331 using Gaussian line fitting and ; linear least-square analysis of the velocity determined from ; the lines and the given distances. ; ; Use : HUBBLE, GALAXIES, HUBCONST, HERROR ; ; Inputs : None. ; ; Opt. Inputs : None. ; ; Outputs : GALAXIES: A structure array containing the following tags ; for each galaxy that was analyzed: ; ; NAME: The name of the galaxy. ; ; DIST: The distance to the galaxy in Mpc. ; ; D_UNC: The uncertainty in this distance in Mpc. ; ; VEL: The radial velocity of the galaxy in km/s. ; ; V_UNC: The uncertainty in this velocity in km/s. ; ; Z: The redshift of the galaxy. ; ; COMMENTS: Any comments concerning the uncertainties, etc. ; ; HUBCONST: Hubble's constant in units of km/s/Mpc. ; ; HERROR: The error in the Hubble constant. ; ; Opt. Outputs: None. ; ; Keywords : SPECDIR: The directory of the spectral data if it is not ; stored in the same subdirectory as this procedure. ; ; DISTDIR: The directory where the distance data file is stored. ; ; GLINE: A structure array of size (NGAL*MAXNLINE) containing ; information of the lines used for the radial velocity ; calculations. Here NGAL is the number of galaxies in ; the sample and MAXNLINE is the maximum number of lines ; used for any of the galaxies in the sample. The tags ; are: ; ; GALAXY: The name of the galaxy = GALAXIES.NAME. ; ; LINE: The ID of the measured line (i.e., Ca II K). ; ; TYPE: Whether the galaxies spectrum is either an ; emission or an absorption line spectrum. ; ; WAVELAB: The laboratory rest wavelength in Angstroms. ; ; WAVEOBS: The observed line wavelength returned by ; GAUSSFIT. ; ; GAUSSWIDTH: The Gaussian width of the line in ; Angstroms = SQRT(2) * SIGMA, the standard ; deviation. ; ; FWHM: The Full-Width-at-Half-Maximum in Angstroms ; = 2.354 * SIGMA. ; ; FLUX_PEAK: The peak flux of the emission line or ; lowest flux level of an absorption line. ; ; VRAD: The radial velocity of the line in km/s. ; ; SAVEDATA: The name of the file to save results of each galaxy ; as it is reduced. This way, users can retrieve ; previously reduced data if a run crashed prematurely. ; The default file name "hubblesavedata.txt" is used ; if this keyword is not set. Note that data is saved ; to an external file whether or not this keyword is set. ; (This has not yet been implemented but will be shortly.) ; ; ERRMSG: If defined and passed, then any error messages will ; be returned to the user in this parameter rather than ; being printed internally by this procedure. If no ; errors are encountered, then a null string is ; returned. In order to use this feature, the string ; ERRMSG must be defined first, e.g., ; ; ERRMSG = '' ; HUBBLE, GALAXIES, HUBCONST, HERROR, $ ; ERRMSG=ERRMSG ; IF ERRMSG(0) NE '' THEN ... ; ; Calls : DECIDEOPT, GALFLSCALE, RDGALDIST, RDGALSPEC, SELECTGAL, ; SPECTYPE. ; ; Common : None. ; ; Restrictions: None. ; ; Side effects: The galaxy spectra are stored in ASCII text files with ; the file extension ".txt" and are stored in the SPECDIR ; subdirectory. The galaxy distance data are stored in a ; file named "galdist.txt" located in the DISTDIR directory. ; Pre-reduced spectrum data is stored in the SAVEDATA file ; in the current directory. ; ; Category : Class project, cosmology, Hubble's law, galaxies, spectra. ; ; Prev. Hist. : None. ; ; Written : Donald G. Luttermoser, ETSU/Physics, 30 Aug 2002. ; ; Modified : Version 1, Donald G. Luttermoser, ETSU/Physics, 30 Aug 2002 ; Initial program. ; Version 2, Donald G. Luttermoser, ETSU/Physics, 29 Jan 2004 ; Minor modifications made to student comments. ; Version 3, Donald G. Luttermoser, ETSU/Physics, 31 Mar 2007 ; Rewrote my complete procedure since my original one was not ; recovered from torgo (my AlphaStation 200) before it died. ; Got rid of the MESSAGE command and changed the internal ; EMESSAGE string variable to EMESS. ; Version 4, Donald G. Luttermoser, ETSU/Physics, 10 Oct 2009 ; Added a GUI to preselect galaxies to be reduced. Added ; an output file to store previous results so one doesn't ; have to start from scratch if a run crashes (see new ; keyword SAVEDATA. Changed DATDIR keyword to SPECDIR and ; added new keyword DISTDIR. Procedure updated for IDL ; version 6.1+. ; ; Version : Version 4, 10 October 2009. ; ;- PRO HUBBLE, GALAXIES, HUBCONST, HERROR, GLINE=GLINE, SPECDIR=SPECDIR, $ DISTDIR=DISTDIR, SAVEDATA=SAVEDATA, ERRMSG=ERRMSG ; ;ON_ERROR, 2 ; Return to caller if error is encountered. EMESS = '' ; Set to non-null string if error is encountered. GALAXIES = {DATA, NAME: '', DIST: 0., D_UNC: 0., VEL: 0., V_UNC: 0., Z: 0., $ COMMENT: ''} HUBCONST = 0. & HERROR = 0. AFONT = [-1, 0, 0] & PASS = 0 ; ; Check input. ; IF N_PARAMS() NE 3 THEN BEGIN EMESS = 'Syntax: HUBBLE, GALAXIES, HUBCONST, HERROR' GOTO, HANDLE_ERROR ENDIF ; ; Check the machine type and IDL version release number. ; TERM = STRUPCASE(!D.NAME) CASE TERM OF 'X': BEGIN OSTYPE = 'Unix' & DELIM = '/' END 'WIN': BEGIN OSTYPE = 'Windows' & DELIM = '\' END ELSE: BEGIN EMESS = 'Unknown windows and machine type: ' + TERM GOTO, HANDLE_ERROR END ENDCASE DPOS = STRPOS(!VERSION.RELEASE, '.') IDLVERS = FLOAT(STRMID(!VERSION.RELEASE, 0, DPOS+2)) ; ; Is the "save data" file name passed? Does the file exist? ; IF KEYWORD_SET(SAVEDATA) THEN BEGIN SAVEFILE = SAVEDATA[0] ASCK = SIZE(SAVEFILE) IF ASCK[N_ELEMENTS(ASCK)-2] NE 7 THEN BEGIN WERR = DIALOG_MESSAGE(['SAVEDATA name must be a scalar string.', $ 'Resetting name to "hubblesavedata.txt".'], /ERROR) SAVEFILE = 'hubblesavedata.txt' ENDIF ENDIF ELSE SAVEFILE = 'hubblesavedata.txt' ; IF IDLVERS LT 6.1 THEN BEGIN CKSAVE = FINDFILE(SAVEFILE, COUNT=QSAVE) ENDIF ELSE BEGIN QRES = DIALOG_MESSAGE(['You are using IDL release version '+ $ !VERSION.RELEASE+'.', 'Please edit "hubble.pro" and remove the', $ 'comment (";") symbol on the FILE_SEARCH', 'function lines.', $ 'Also do this in the "rdgaldist.pro" file.'], /INFORMATION) ; CKSAVE = FILE_SEARCH(SAVEFILE, COUNT=QSAVE) ENDELSE ; ; Retrieve spectrum file names. ; IF KEYWORD_SET(SPECDIR) THEN BEGIN SFILES = SPECDIR & SLEN = STRLEN(SFILES) IF STRMID(SFILES, SLEN-1, 1) NE DELIM THEN SFILES = SFILES + DELIM SFILES = SFILES + '*.txt' ENDIF ELSE BEGIN SPECDIR = '' & SFILES = SPECDIR + '*.txt' ENDELSE ; ; Check to see if the spectra files exist. ; JSS = 0 GETSPECFILES: ; ; This will prevent infinite loops from occurring if the user is ; having trouble. ; IF JSS GT 3 THEN BEGIN EMESS = 'Unable to ascertain location of galaxy spectra.' GOTO, HANDLE_ERROR ENDIF IF IDLVERS LT 6.1 THEN BEGIN CKSFILE = FINDFILE(SFILES, COUNT=NCT) ENDIF ELSE BEGIN ; CKSFILE = FILE_SEARCH(SFILES, COUNT=NCT) ENDELSE ; ; This assumes there are at least 5 spectrum files present. ; IF NCT LT 5 THEN BEGIN IF SPECDIR EQ '' THEN STXT = 'current' ELSE STXT = SPECDIR WRES = DIALOG_MESSAGE(['Galaxy spectrum files not found in ', $ STXT + ' directory.', $ 'You will be asked to select a new directory.']) WAIT, 1 SPECDIR = DIALOG_PICKFILE(/DIRECTORY, TITLE='Please select the galaxy '+ $ 'spectra directory:') SFILES = SPECDIR + '*.txt' JSS = JSS + 1 GOTO, GETSPECFILES ENDIF ; ; In case other ".txt" are located in the directory where the spectra are ; located, only keep those files whose filenames start with "ic", "mrk", ; "ngc", and "ugc". ; ICGALS = WHERE(STRPOS(CKSFILE, DELIM+'ic') GT 0, NIC) IF NIC GT 0 THEN SPECFILES = CKSFILE[ICGALS] ; MRKGALS = WHERE(STRPOS(CKSFILE, DELIM+'mrk') GT 0, NMRK) IF NMRK GT 0 THEN SPECFILES = [SPECFILES, CKSFILE[MRKGALS]] ; NGCGALS = WHERE(STRPOS(CKSFILE, DELIM+'ngc') GT 0, NNGC) IF NNGC GT 0 THEN SPECFILES = [SPECFILES, CKSFILE[NGCGALS]] ; UGCGALS = WHERE(STRPOS(CKSFILE, DELIM+'ugc') GT 0, NUGC) IF NUGC GT 0 THEN SPECFILES = [SPECFILES, CKSFILE[UGCGALS]] ; NCT = NIC + NMRK + NNGC + NUGC ; Total number of galaxies to reduce. ; ; Note that the galaxy distances data file is searched for and read in ; procedure RDGALDIST. ; IF KEYWORD_SET(DISTDIR) THEN BEGIN DFILE = DISTDIR & DLEN = STRLEN(DFILE) IF STRMID(DFILE, DLEN-1, 1) NE DELIM THEN DFILE = DFILE + DELIM DFILE = DFILE + 'galdist.txt' ENDIF ELSE BEGIN DISTDIR = '' & DFILE = DISTDIR + 'galdist.txt' ENDELSE ; ; Read in the distances (in Mpc) to the galaxies. This data is stored in an ; IDL structure array. GDIST is a structure array, see rdgaldist.pro. ; RDGALDIST, GDIST, NGAL, DISTFILE=DFILE, ERRMSG=EMESS IF EMESS NE '' THEN GOTO, HANDLE_ERROR IF NGAL NE NCT THEN BEGIN EMESS = ['Number of galaxies listed in "distance" file (' + $ STRTRIM(NGAL, 2) + ') is not the same as', 'those listed ' + $ 'in the spectra subdirectory (' + STRTRIM(NCT, 2) + ').'] GOTO, HANDLE_ERROR ENDIF ; ; Sort the data in GDIST such that it is in the same order as the order of the ; files. ; GINDSORT = INTARR(NGAL) FOR I = 0, NGAL-1 DO BEGIN ISRT = WHERE(STRPOS(STRLOWCASE(SPECFILES), GDIST[I].NAME) GT -1, NSRT) IF NSRT EQ 0 THEN BEGIN EMESS = 'The galaxy ' + GDIST[I].NAME + ' does not have ' + $ 'a spectrum!' GOTO, HANDLE_ERROR ENDIF IF NSRT GT 1 THEN BEGIN EMESS = 'The galaxy ' + GDIST[I].NAME + ' seems to have ' + $ 'more than one spectrum!' GOTO, HANDLE_ERROR ENDIF GINDSORT[I] = ISRT[0] ENDFOR SPECFILES = SPECFILES[GINDSORT] ; ; Map the entries in GDIST into the GALAXIES structure at this point. ; GALAXIES = REPLICATE(GALAXIES, NGAL) GALAXIES[*].NAME = GDIST[*].NAME GALAXIES[*].DIST = GDIST[*].D GALAXIES[*].D_UNC = GDIST[*].UNC ; ; Have user select which galaxy spectra will be analyzed. ; (Still to be added.) ; ;SELECTGAL, GALAXIES, SAVEFILE, QSAVE, 0, GSELECT, ERRMSG=EMESS ;IF EMESS NE '' THEN GOTO, HANDLE_ERROR ;STOP ; ; Set up internal parameters. ; STDFWHM = [5.0, 10.0] ; ; Emission lines (note that all students only have to do the first 3): ; ;STDELINES = ['H-beta', '[O III]', '[O III]', 'H-alpha', '[N II]'] STDELINES = ['H-beta', '[O III]', '[O III]'] ; ; Note that there is a weaker [N II] just shortward of H-alpha at 6548.1. ; ;STDEWAVE = [4861.3, 4958.9, 5006.9, 6562.8, 6583.4] STDEWAVE = [4861.3, 4958.9, 5006.9] NELINES = N_ELEMENTS(STDEWAVE) ; ; Absorption lines (note that all students only have to do the Ca II lines): ; ;STDALINES = ['Ca II K', 'Ca II H', 'Na I D2', 'Na I D1'] STDALINES = ['Ca II K', 'Ca II H'] ; ; Note that the sodium D-lines are blended. ; ;STDAWAVE = [3933.7, 3968.5, 5890.0, 5895.9] STDAWAVE = [3933.7, 3968.5] NALINES = N_ELEMENTS(STDAWAVE) NLINES = MAX([NELINES, NALINES]) ; ; The next part of the code will take a while to accomplish. Each spectrum ; file is read in and then analyzed, one at a time. Results are stored in the ; GLINE structure array. The FLUX_PEAK tab is not needed for this project, but ; we'll save the result anyway since GAUSSFIT returns it. ; GLINE = {SPECTRUM, GALAXY: '', LINE: '', TYPE: '', WAVELAB: 0.0D0, WAVEOBS: 0.0D0, $ GAUSSWIDTH: 0., FWHM: 0.0, FLUX_PEAK: 0.0, VRAD: 0.} GLINE = REPLICATE(GLINE, NCT*NLINES) ; ; Setup plot appearance characteristics. ; IF OSTYPE EQ 'Unix' THEN BEGIN !P.CHARSIZE = 2.2 & !X.CHARSIZE = 1.0 & !Y.CHARSIZE = 1.0 ENDIF ELSE BEGIN !P.CHARSIZE = 1.6 & !X.CHARSIZE = 1.0 & !Y.CHARSIZE = 1.0 ENDELSE !P.THICK = 2 & !X.THICK = 2 & !Y.THICK = 2 ; ; Bring up a plot window in the upper-right corner of the screen. ; WINDOW, 0, TITLE='Gaussian Fitting Plot Window' ; ; We are going to require to save one Gauss fit plot for an absorption line ; and one for an emission line. The following variables are set to yes when ; we chose an appropriate plot to save. ; QEPLOT = 'no' & QAPLOT = 'no' ; ; Read in each spectrum and analyze. ; FOR I = 0, NGAL-1 DO BEGIN ; ; First, read the spectrum and scale the flux. ; PRINT, 'Reading spectrum: ' + CKSFILE[I] ; ; Make a call to a spectrum reading procedure here. ; RDGALSPEC, CKSFILE[I], GALNAME, WAVE, FLUX, ERRMSG=EMESS IF EMESS NE '' THEN GOTO, HANDLE_ERROR ; ; Make a call to a flux scale factor procedure here (return FSC and FLAB). ; GALFLSCALE, FLUX, FSC, FLAB, ERRMSG=EMESS IF EMESS NE '' THEN GOTO, HANDLE_ERROR ; ; Plot the full spectrum and ask whether it is an absorption or emission line ; spectrum. Note that the WSHOW commands are placed before each PLOT command to ; insure that the plot window is clearly visible. ; WSHOW PLOT, WAVE, FLUX*FSC, XTITLE='!6Wavelength (A)', YTITLE=FLAB, $ TITLE='!17'+GALNAME, FONT=AFONT[PASS] STYPE = SPECTYPE() ; STARTAGAIN: ; IF STYPE EQ 'absorption' THEN BEGIN NLINES = NALINES & STDLINES = STDALINES & STDWAVE = STDAWAVE ENDIF ELSE BEGIN NLINES = NELINES & STDLINES = STDELINES & STDWAVE = STDEWAVE ENDELSE VRARR = DBLARR(NLINES) ; Temporary array to store radial velocities. ; ; Analyze each of the spectral lines for each of the galaxies. ; In this section, fit each line of interest to a Gaussian using the GAUSSFIT ; procedure. ; WVBLUE = -200. & WVRED = 200. ; FOR J = 0, NLINES-1 DO BEGIN ; REPLOTLINE: ; ; Make a call to a flux scale factor procedure here (return FSC and FLAB). ; We need to do this each time since different parts of the spectrum may have vastly ; different flux levels and we want to keep the y-axis label accurate ; ILEFT = WHERE(WAVE GE STDWAVE[J]+WVBLUE, NLEFT) IF NLEFT EQ 0 THEN ILEFT = 0 ELSE ILEFT = ILEFT[0] IRIGHT = WHERE(WAVE GT STDWAVE[J]+WVRED, NRIGHT) IF NRIGHT EQ 0 THEN IRIGHT = N_ELEMENTS(WAVE)-1 ELSE $ IRIGHT = (IRIGHT[0]-1) > ILEFT+1 GALFLSCALE, FLUX(ILEFT:IRIGHT), FSC, FLAB, ERRMSG=EMESS IF EMESS NE '' THEN GOTO, HANDLE_ERROR ; PLOT, WAVE, FLUX*FSC, XTITLE='!6Wavelength (A)', YTITLE=FLAB, $ XRANGE=[STDWAVE[J]+WVBLUE,STDWAVE[J]+WVRED], TITLE='!17'+GALNAME, $ FONT=AFONT[PASS] ; PRINT, '**********************' PRINT, 'Move cursor to the '+STDLINES[J]+' and click left mouse button,' PRINT, ' or click right mouse button to move plot window redward,' PRINT, ' or click center button to move plot window blueward.' WSHOW CURSOR, XL, YL WAIT, 1.0 CASE !ERR OF 4: BEGIN ; Right button = Move redward. WVBLUE = WVBLUE + 300. WVRED = WVRED + 300. GOTO, REPLOTLINE END 2: BEGIN ; Center button = Move blueward. WVBLUE = WVBLUE - 300. WVRED = WVRED - 300. GOTO, REPLOTLINE END ; ; Gaussian stuff here. Save your Gaussian output in arrays for later use -- ; perhaps in the array structure GLINE. ; ELSE: BEGIN ; RESETLIMITS: ; ; Plot the line in question in a +/-40 Angstrom range from the cursor click. ; PLOT, WAVE, FLUX*FSC, XTITLE='!6Wavelength (A)', YTITLE=FLAB, $ XRANGE=[XL-40.,XL+40.], FONT=AFONT[PASS], $ TITLE='!17'+GALNAME+': '+STDLINES[J]+' Line' ; PRINT, ' ' PRINT, '--------------' PRINT, 'Move cursor to LEFT side of '+STDLINES[J]+' and click mouse.' WSHOW CURSOR, XB, YB OPLOT, [XB, XB], [MAX(FLUX*FSC), MIN(FLUX*FSC)], LINE=2 WAIT, 0.5 ; PRINT, '++++++++++++++' PRINT, 'Move cursor to RIGHT side of '+STDLINES[J]+' and click mouse.' WSHOW CURSOR, XR, YR OPLOT, [XR, XR], [MAX(FLUX*FSC), MIN(FLUX*FSC)], LINE=2 ; ; Isolate the line from the rest of the spectrum. ; IBLUE = WHERE(WAVE GE XB, NBLUE) IF NBLUE EQ 0 THEN IBLUE = 0 ELSE IBLUE = IBLUE[0] IRED = WHERE(WAVE GT XR, NRED) IF NRED EQ 0 THEN IRED = N_ELEMENTS(WAVE)-1 ELSE $ IRED = (IRED[0]-1) > IBLUE+1 ; ; Fit the profile to a Gaussian. ; FGFIT = GAUSSFIT(WAVE[IBLUE:IRED], FLUX[IBLUE:IRED]*FSC, ACOEFF, $ NTERMS=6) OPLOT, WAVE[IBLUE:IRED], FGFIT, PSYM=2 ; ; Ask the user to acknowledge that the fit is good. Calculate the radial velocity ; to show to the user. ; DELLAMB = ACOEFF[1] - STDWAVE[J] VRAD = 2.997925D5 * DELLAMB / STDWAVE[J] QRES = DIALOG_MESSAGE([STDLINES[J]+' has a radial velocity of '+$ STRTRIM(STRING(VRAD, FORMAT='(F10.2)'), 2)+' km/s and a', $ 'Gaussian width of '+STRTRIM(STRING(SQRT(2.)*ACOEFF[2], $ FORMAT='(F10.2)'),2)+' Angstroms.', $ 'Is the Gaussian fit acceptable?'], /QUESTION) QSAVE = 'no' IF QRES EQ 'Yes' THEN BEGIN QSAVE = 'yes' ENDIF ELSE BEGIN QRES = DECIDEOPT(STYPE) CASE STRMID(QRES, 0, 6) OF 'do not': ; Do not save data, continue to next line. 'save r': QSAVE = 'yes' ; Save the results anyway. 'choose': GOTO, RESETLIMITS ; Reset the line profile boundaries. 'change': BEGIN ; Reset spectrum type. IF STYPE EQ 'absorption' THEN STEMP = 'emission' IF STYPE EQ 'emission' THEN STEMP = 'absorption' STYPE = STEMP GOTO, STARTAGAIN END 'wrong ': GOTO, REPLOTLINE ELSE: ; An error was encountered, do not save data. ENDCASE ENDELSE IF QSAVE EQ 'yes' THEN BEGIN GLINE[I*3+J].WAVEOBS = ACOEFF[1] GLINE[I*3+J].WAVELAB = STDWAVE[J] VRARR[J] = VRAD GLINE[I*3+J].VRAD = VRAD GLINE[I*3+J].GAUSSWIDTH = SQRT(2.) * ACOEFF[2] GLINE[I*3+J].FWHM = 2.354 * ACOEFF[2] GLINE[I*3+J].TYPE = STYPE GLINE[I*3+J].LINE = STDLINES[J] GLINE[I*3+J].GALAXY = GALNAME ; ; Save the Gaussfit plot. One for absorption and one for emission. ; ; IF (STYPE EQ 'absorption') AND (QAPLOT EQ 'no') THEN BEGIN ; QPSAV = DIALOG_MESSAGE(['Do you wish ??? ; ENDIF ENDIF END ENDCASE ENDFOR ; ; Calculate average radial velocities for all of the lines in this galaxy spectrum ; and store it along with the uncertainty in the GALAXIES structure array. ; IUSE = WHERE(VRARR NE 0., NUSE) IF NUSE GT 0 THEN BEGIN VRAD = TOTAL(VRARR) / NUSE VUNC = MAX(ABS(VRARR)) - ABS(VRAD) GALAXIES[I].VEL = VRAD GALAXIES[I].V_UNC = VUNC GALAXIES[I].Z = VRAD / 2.997925D5 IF NUSE EQ 1 THEN SLE = ' line' ELSE SLE = ' lines' GALAXIES[I].COMMENT = STRTRIM(STRING(NUSE), 2) + SLE + $ ' used in v(rad) calculation' ENDIF ELSE GALAXIES[I].COMMENT = 'No useful lines found in spectrum for ' + $ 'v(rad) calculation' ; ENDFOR ; WAIT, 1.0 ; Wait a bit so we can see the last spectrum plotted. ; ; At this point, make a plot of the average radial velocity for each galaxy ; as a function of its distance (and include the error bars in both axes) both ; to the terminal and to a hardcopy device. ; FINALPLOT: ; IF PASS EQ 0 THEN WSHOW IUSE = WHERE(GALAXIES.VEL NE 0., NUSE) PLOT, GALAXIES[IUSE].DIST, GALAXIES[IUSE].VEL/1000, XTITLE='!6Distance (Mpc)', $ YTITLE='!6Radial Velocity (1000 km/s)', TITLE="!17Hubble's Law Observations", $ PSYM=2, FONT=AFONT[PASS] ; ; Now draw uncertainty lines for each point. ; FOR I = 0, NUSE-1 DO BEGIN J = IUSE[I] & VUNC = GALAXIES[J].V_UNC & DUNC = GALAXIES[J].D_UNC ; ; First velocity uncertainty, then distance uncertainty. ; OPLOT, [GALAXIES[J].DIST, GALAXIES[J].DIST], (GALAXIES[J].VEL+[-VUNC, VUNC])/1000 OPLOT, GALAXIES[J].DIST+[-DUNC, DUNC], [GALAXIES[J].VEL, GALAXIES[J].VEL]/1000 ENDFOR ; ; Now use LINFIT to use a least squares analysis to determine the slope of the ; line (this will be Hubble's constant) and its uncertainty. Convolve this ; uncertainty with the uncertainties in velocity and distance. ; IF PASS EQ 0 THEN BEGIN IF IDLVERS LT 6.1 THEN BEGIN HRES = LINFIT(GALAXIES[IUSE].DIST, GALAXIES[IUSE].VEL, CHISQ=CHISQ, $ PROB=PROB, SIGMA=HSIGMA) HLAW = HRES[0] + HRES[1] * GALAXIES[IUSE].DIST ENDIF ELSE BEGIN HRES = LINFIT(GALAXIES[IUSE].DIST, GALAXIES[IUSE].VEL, CHISQ=CHISQ, $ PROB=PROB, SIGMA=HSIGMA, YFIT=HLAW) ENDELSE ; HUBCONST = HRES[1] ; This is Hubble's constant that is returned to outside. ; ; Calculate total uncertainty. ; ; !!!!! STUDENT'S WRITTEN SECTION: ; I have left in the calculation of the normalized error for the ; velocity component of the uncertainty in Hubble's constant. ; Students, you are to write a line under VERRNORM for a calculation ; of the normalized error for the uncertainty in the distance called ; DERRNORM, and the normalized error for the slope of the fitted ; straight line to the data called HERRNORM. These 3 terms are then ; used in the equation to calculate the uncertainty in Hubble's ; constant = HERROR. ; VERRNORM = TOTAL((GALAXIES[IUSE].V_UNC/GALAXIES[IUSE].VEL)^2) ; DERRNORM = (put your distance error equation here) ; HERRNORM = (put your slope error equation here) ; ; !!!!! END OF STUDENT WRITTEN SECTION. ; HERROR = SQRT(VERRNORM + DERRNORM + HERRNORM) * HUBCONST ; Total error in H. ; PRINT, ' ' PRINT, ">>>>> Hubble's Constant for fitted line: " + STRTRIM(STRING(HUBCONST, $ FORMAT='(F9.1)'), 2) + ' km/s/Mpc +/- ' + STRTRIM(STRING(HERROR, $ FORMAT='(F9.1)'), 2) + ' km/s/Mpc' PRINT, '>>>>> Chi-Square = ', CHISQ PRINT, '>>>>> Goodness-of-fit Probabilty = ', PROB IF PROB GT 0.1 THEN PRINT, '>>>>> The linear fit to the data is "believable"!' $ ELSE PRINT, '>>>>> The linear fit to this data should NOT be trusted!' ENDIF ; ; Plot the fit to the data. ; OPLOT, GALAXIES[IUSE].DIST, HLAW/1000 CASE PASS OF 1: BEGIN PTSIZE = 1.4 XT1 = 0.50 & YT1 = 0.30 XB1 = 0.60 & YB1 = 0.23 END 2: BEGIN PTSIZE = 1.6 XT1 = 0.50 & YT1 = 0.30 XB1 = 0.60 & YB1 = 0.25 END ELSE: BEGIN PTSIZE = 1.4 XT1 = 0.45 & YT1 = 0.30 XB1 = 0.50 & YB1 = 0.25 END ENDCASE XYOUTS, XT1, YT1, '!6H!Do!N = ' + STRTRIM(STRING(HUBCONST, FORMAT='(F9.1)'), 2) + $ ' km/s/Mpc', /NORM, SIZE=PTSIZE, FONT=AFONT[PASS] XYOUTS, XB1, YB1, '!6+/- ' + STRTRIM(STRING(HERROR, FORMAT='(F9.1)'), 2) + $ ' km/s/Mpc', /NORM, SIZE=PTSIZE, FONT=AFONT[PASS] ; ; Save to hardcopy file. First make an encapsulated PS plot then a standard PS plot. ; CASE PASS OF 0: BEGIN ; Encapsulated postscript file. SET_PLOT, 'PS' DEVICE, /PORTRAIT, /ENCAPSULATE, FILENAME='hubblelaw.eps' DEVICE, /SCHOOLBOOK, /BOLD, FONT_INDEX=17 DEVICE, /SYMBOL, FONT_INDEX=7 DEVICE, /BKMAN, /DEMI, FONT_INDEX=6 DEVICE, /HELVETICA, /BOLD, FONT_INDEX=5 !P.CHARSIZE = 1.2 & !X.CHARSIZE = 1.0 & !Y.CHARSIZE = 1.0 !P.THICK = 5 & !X.THICK = 5 & !Y.THICK = 5 PASS = 1 GOTO, FINALPLOT END 1: BEGIN ; Normal postscript file. DEVICE, /CLOSE_FILE PRINT, ' ' PRINT, '!!! Encapsulated postscript file stored in: "hubblelaw.eps"' PRINT, '!!! Use this file in your LaTeX manuscript.' PRINT, ' ' SET_PLOT, 'PS' DEVICE, /LANDSCAPE, FILENAME='hubblelaw.ps' DEVICE, /SCHOOLBOOK, /BOLD, FONT_INDEX=17 DEVICE, /SYMBOL, FONT_INDEX=7 DEVICE, /BKMAN, /DEMI, FONT_INDEX=6 DEVICE, /HELVETICA, /BOLD, FONT_INDEX=5 !P.CHARSIZE = 2.2 & !X.CHARSIZE = 1.0 & !Y.CHARSIZE = 1.0 !P.THICK = 5 & !X.THICK = 5 & !Y.THICK = 5 PASS = 2 GOTO, FINALPLOT END 2: BEGIN DEVICE, /CLOSE_FILE PRINT, ' ' PRINT, '!!! Normal postscript file stored in: "hubblelaw.ps"' PRINT, '!!! You can send this plot directly to a printer.' PRINT, ' ' SET_PLOT, TERM END ENDCASE ; ; Let the user know that the program is complete. ; PRINT, ' ' PRINT, '*********************************************' PRINT, 'MYHUBBLE procedure has successfully finished!' PRINT, '*********************************************' ; IF EMESS NE '' THEN GOTO, HANDLE_ERROR IF N_ELEMENTS(ERRMSG) GT 0 THEN ERRMSG=EMESS RETURN ; ; Error handling portion of the procedure. ; HANDLE_ERROR: ; IF N_ELEMENTS(ERRMSG) EQ 0 THEN PRINT, EMESS ERRMSG = EMESS RETURN ; END