From 6aaae42b09892495615bae80079bbe11ac98bd4b Mon Sep 17 00:00:00 2001 From: JPVentura Date: Wed, 5 Apr 2017 11:55:36 -0400 Subject: [PATCH 1/5] added updated files --- EqWi.py | 54 ++++++------ README.md | 12 +-- ewidth.py | 244 +++++++++++++++++++++--------------------------------- 3 files changed, 125 insertions(+), 185 deletions(-) diff --git a/EqWi.py b/EqWi.py index 663ba12..04d7256 100644 --- a/EqWi.py +++ b/EqWi.py @@ -1,42 +1,40 @@ import matplotlib.pyplot as plt -import numpy as np +import numpy as numpy from astropy.io import fits -import pyspeckit datfile = fits.open('filename') hdr = datfile[1].header -"Kelle was here" -#try: -# data = datfile[1].data -# except IndexError: -# pass -# else: -# good = False - names = datfile[1].dtype.names + +try: + data = datfile[1].data + except IndexError: + pass + else: + good = False + names = data.dtype.names if 'loglam' in names and 'flux' in names: flux = datfile[1].data['flux'] - wlen = 10**np.array(datflike[1].data['loglam']) + wlen = 10**np.array(data.loglam) indx = np.where((wlen > 8250) & (wlen < 8350)) avgval = flux[indx].mean() nflux = flux/avgval - - # good = True + + good = True #elif 'wavelength' in names and 'flux' in names: - # wa = data.wavelength - # fl = data.flux - # good = True + wa = data.wavelength + fl = data.flux + good = True #if good: - # er = np.ones_like(flux) - # try: - # err = datfile[1].data['ivar'] - # except AttributeError: - # pass - # co = np.empty_like(fl) * np.nan - # try: - # co = data.co - # except AttributeError: - # pass - # return Spectrum(wa=wa, fl=fl, er=er, co=co, filename=filename) -# + er = np.ones_like(fl) + try: + er = data.er + except AttributeError: + pass + co = np.empty_like(fl) * np.nan + try: + co = data.co + except AttributeError: + pass + return Spectrum(wa=wa, fl=fl, er=er, co=co, filename=filename) diff --git a/README.md b/README.md index b8d0369..690a80c 100644 --- a/README.md +++ b/README.md @@ -1,10 +1,10 @@ -# Equivalent_widths +# EqWi (Equivalent Widths): -Python code to calculate the equivalent widths given spectrum of various file extensions using [PySpecKit](https://github.com/pyspeckit/pyspeckit). By: Jean-Paul Ventura +Python script to calculate the equivalent widths of spectral line emission & absorption features for various spectral data file types using the pyspeckit module [PyspecKit](https://github.com/pyspeckit/pyspeckit). By Jean-Paul Ventura -### Information -Created on February 16, 2017 +### Information ### -### Referencing +### Referencing ### + +### License ### -### License diff --git a/ewidth.py b/ewidth.py index 8203f65..ddd423d 100644 --- a/ewidth.py +++ b/ewidth.py @@ -1,4 +1,4 @@ - + # Jean-paul Ventura # January 12, 2017 # Modified: February 20th, 2016 @@ -6,177 +6,119 @@ from astropy.io import fits import numpy as np + # Calculate equivalent width def eqwidth(filename,wavelength1,wavelength2,wavelength3,wavelength4): -""" -This function measures the equivalent width of spectral emission/absorption features. + ''' -==================== -Function parameters: -==================== + This function measures the equivalent width of spectral emission/absorption features. + + ==================== + Function parameters: + ==================== + + filename - SDSS spectra .fits filename entered as a string. You may specify entire + path as a string, or if file is located in the current working directory, + simply the filename (e.g. 'SDSS-12345.fits') + + wavelength1 - wavelength value where continuum region to the left of the spectral feature begins. + + wavelength2 - wavelength value where continuum region to the left of the spectral feature ends. + + wavelength3 - wavelength value where continuum region to the right of the spectral feature begins. + + wavelength4 - wavelength value where continuum region to the right of the spectral feature ends. + + ================= + Function returns: + ================= + + The measured equivalent width of the spectral emission/absorption feature -filename - SDSS spectra .fits filename entered as a string. You may specify entire - path as a string, or if file is located in the current working directory, - simply the filename (e.g. 'SDSS-12345.fits') + ''' -wavelength1 - wavelength value where continuum region to the left of the spectral feature begins. + # Open .fits file + hdu = fits.open(filename) -wavelength2 - wavelength value where continuum region to the left of the spectral feature ends. + # Define wavelength and flux arrays by accessing .fits sloan data + lmbda = 10**np.array(hdu[1].data['loglam']) + flux = hdu[1].data['flux'] -wavelength3 - wavelength value where continuum region to the right of the spectral feature begins. -wavelength4 - wavelength value where continuum region to the right of the spectral feature ends. + # Normalize flux by dividing flux array by mean flux value over relatively quiet + # continuum region + indx = np.where((lmbda > 8250) & (lmbda < 8350)) + avgval = flux[indx].mean() + nflux = flux/avgval -================= -Function returns: -================= -The measured equivalent width of the spectral emission/absorption feature -""" + # Index for continuum region where line feature is found + fuindx = (lmbda >= wavelength1) & (lmbda <= wavelength4) - hdu = fits.open(filename) + # Index for the line feature itself + coindx = ((lmbda >= wavelength1) & (lmbda <= wavelength2)) | ((lmbda >= wavelength3) & (lmbda <= wavelength4)) - lambda = 10**np.array(hdu[1].data['loglam']) - flux = hdu[1].data['flux'] - - #designate fundamental and continuum indices - fuind = (lambda >= wavelength1) & (lambda <= wavelength4) - coind = ((lambda >= wavelength1) & (lambda <= wavelength2)) | ((lambda >= wavelength3) & (lambda <= wavelength4)) - - x1 = lambda[coind] - y1 = flux[coind] + # Isolate the line feature wavelength region + x1 = lmbda[coindx] + + #Isolate the normalized flux over the line feature wavelength region + y1 = nflux[coindx] + + #fit 1st degree polynomial to the line feature m,c = np.polyfit(x1,y1,1) - # #print m,c - # rms = sqrt( mean( (m*x1 + c - y1)**2 ) ) - indx = np.where((lambda > 8250) & (lambda < 8350)) - avgval = flux[indx].mean() - nflux = flux/avgval - lambda = lambda[fuind] - flux = flux[fuind] + #print m,c + #rms = sqrt( mean( (m*x1 + c - y1)**2 ) ) + + lmbda = lmbda[fuindx] + flux = nflux[fuindx] - dlambda= lambda[2]-lambda[1] - fitl = m*lambda + c + + dlmbda= lmbda[2]-lmbda[1] + fitl = m*lmbda + c - ewind = np.where((lambda >= w2) & (lambda <= w3)) - eqwi = sum(1-nflux[ewind])*dlambda + ewindx = np.where((lmbda >= wavelength2) & (lmbda <= wavelength3)) + eqwi = sum(1-nflux[ewindx])*dlmbda # For error #fx1, err = mcerr(wav,dat,w2,w3,m,c,rms) return eqwi -#========================= - -# def mcerr(wav,dat,w1,w2,m,c,rms): -# siz = 3000 -# guess = [] -# dlam = wav[2]-wav[1] -# -# for i in range(siz): -# fitl = m*wav + c + uniform(-1,1)*rms -# ndat = dat / fitl -# ewind = (wav >= w1) & (wav <= w2) -# fx = sum(1-ndat[ewind]) *dlam -# guess.append(fx) -# -# return mean(guess),std(guess) - -# ========================= - - - - -##hdu.close() -# -## Halpha -#w1 = 6548 -#w2 = 6556. -#w3 = 6570 -#w4 = 6580 -#st = 2 # spread to consider -#cc = 0.002 - -#value = eqwidth(wave,data,w1,w2,w3,w4) -#print value -#val1, err1 = mcew(wave,data,w1,w2,w3,w4,st,cc) -#val1 = ewidth(wave,data,w1,w2,w3,w4) -#print val1 - -# -#fuind = (wave >= w1) & (wave <= w4) -#coind = ((wave >= w1) & (wave <= w2)) | ((wave >= w3) & (wave <= w4)) -# -#x1 = wave[coind] -#y1 = data[coind] -#m,c = polyfit(x1,y1,1) -# -#wav1 = wave[fuind] -#dat1 = data[fuind] -#fitl1 = m*wav1 + c -#ndat1 = dat1 / fitl1 -# -# -#fuind = (wave >= 6540) & (wave <= 6750) -#data = data[fuind] -#wave = wave[fuind] -# -#plot(wave,data,"k") -#plot(wav1,dat1,"b") -#plot(wav1,fitl1,"r") -# -#axvspan(w2,w3,facecolor="0.8") - - - -# Li -#w1 = 6700.5 -#w2 = 6706. -#w3 = 6711.5 -#w4 = 6713. -#st = 1. -#cc = 0.000 -# -##val2, err2 = mcew(wave,data,w1,w2,w3,w4,st,cc) -#val2, err2 = ewidth(wave,data,w1,w2,w3,w4) -#print val2, err2 -# -# -#fuind = (wave >= w1) & (wave <= w4) -#coind = ((wave >= w1) & (wave <= w2)) | ((wave >= w3) & (wave <= w4)) -# -#x1 = wave[coind] -#y1 = data[coind] -#m,c = polyfit(x1,y1,1) -# -#wav2 = wave[fuind] -#dat2 = data[fuind] -#fitl2 = m*wav2 + c -#ndat2 = dat2 / fitl2 -# -# -#plot(wav2,dat2,"b") -#plot(wav2,fitl2,"r") -#axvspan(w2,w3,facecolor="0.8") -# -#tout = '%s: %2.3f+/-%0.3f, %1.3f+/-%0.3f' % (name[0], val1, err1, val2, err2) -# -#xlabel('Wavelength (Ang)') -#ylabel('Flux') -#title(tout) -#ranges = axis() -#x = [ranges[0],ranges[1]] -#x = [6540,6750] -#y = [ranges[2],ranges[3]] -##y = [-70,-20] -##x.reverse() -#ranges = concatenate((x,y)) -#axis(ranges) -#grid() -# -#outname = name[0]+'_'+str(num)+'.png' -##savefig(outname) - -show() + +#Calculate line height of spectral emission feature +def line_height(filename,wavelength1,wavelength2,wavelength3,wavelength4): + + # Open sloan .fits file + hdu = fits.open(filename) + + #Construct wavelength (lmbda) and flux arrays from data + lmbda = 10**np.array(hdu[1].data['loglam']) + flux = hdu[1].data['flux'] + + # Normalize flux by dividing flux array by mean flux value over relatively quiet + # continuum region + indx = np.where((lmbda > 8250) & (lmbda < 8350)) + avgval = flux[indx].mean() + nflux = flux/avgval + + # designate fundamental and continuum indices + fuind = (lmbda >= wavelength1) & (lmbda <= wavelength4) + coind = ((lmbda >= wavelength1) & (lmbda <= wavelength2)) | ((lmbda >= wavelength3) & (lmbda <= wavelength4)) + + lmbda = lmbda[fuind] + flux = nflux[fuind] + x1 = lmbda[coind] + y1 = nflux[coind] + m,c = np.polyfit(x1,y1,1) + #print m,c + rms = sqrt( mean( (m*x1 + c - y1)**2 ) ) + + + lnheight = max(y1)-min(y1) + + return lnheight + + From ee2c1194f7c6c405e2188312d580e3b09cbbc3f9 Mon Sep 17 00:00:00 2001 From: JPVentura Date: Tue, 23 May 2017 22:59:19 -0400 Subject: [PATCH 2/5] added docstrings and revised comments for all functions --- ewidth.py | 167 +++++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 128 insertions(+), 39 deletions(-) diff --git a/ewidth.py b/ewidth.py index ddd423d..d5deb04 100644 --- a/ewidth.py +++ b/ewidth.py @@ -7,9 +7,8 @@ import numpy as np -# Calculate equivalent width def eqwidth(filename,wavelength1,wavelength2,wavelength3,wavelength4): - ''' + """ This function measures the equivalent width of spectral emission/absorption features. @@ -33,9 +32,9 @@ def eqwidth(filename,wavelength1,wavelength2,wavelength3,wavelength4): Function returns: ================= - The measured equivalent width of the spectral emission/absorption feature + eqwi - The measured equivalent width of the spectral emission/absorption feature - ''' + """ # Open .fits file hdu = fits.open(filename) @@ -45,51 +44,128 @@ def eqwidth(filename,wavelength1,wavelength2,wavelength3,wavelength4): flux = hdu[1].data['flux'] - # Normalize flux by dividing flux array by mean flux value over relatively quiet - # continuum region - indx = np.where((lmbda > 8250) & (lmbda < 8350)) - avgval = flux[indx].mean() - nflux = flux/avgval - - # Index for continuum region where line feature is found + # Index for continuum region where entire line feature is found (wings + emission line) fuindx = (lmbda >= wavelength1) & (lmbda <= wavelength4) - # Index for the line feature itself + # Index for both of the wings of the line feature coindx = ((lmbda >= wavelength1) & (lmbda <= wavelength2)) | ((lmbda >= wavelength3) & (lmbda <= wavelength4)) - # Isolate the line feature wavelength region + # Isolate the wavelength subarrays pertaining to both of the wings of the line feature x1 = lmbda[coindx] - #Isolate the normalized flux over the line feature wavelength region + #Isolate the normalized flux subarrays of the line features wings. y1 = nflux[coindx] - #fit 1st degree polynomial to the line feature + # Generate first degree polynomial coefficients that correspond to the line wings m,c = np.polyfit(x1,y1,1) - #print m,c - #rms = sqrt( mean( (m*x1 + c - y1)**2 ) ) - lmbda = lmbda[fuindx] - flux = nflux[fuindx] + # calculate the RMSE for the baseline for later use to generate error array. + rmse = sqrt( mean( (m*x1 + c - y1)**2 ) ) + # Assign the wavelength and normalized flux arrays over the entire feature (wings + line) to new variables + lmbda = lmbda[fuindx] + flux = flux[fuindx] - dlmbda= lmbda[2]-lmbda[1] - fitl = m*lmbda + c + # using the coefficients of the deg(1) polynomial generated earlier from the wings, fit a line to them + # in order to establish the baseline for the line feature. + fitl = np.array(m*lmbda + c) + # Normalize the flux array by dividing it by the fitline. + nflux = flux/fitline + # assign the differential wavelength element to a variable. + dlmbda= lmbda[2]-lmbda[1] + + # Create index for the wavelength of the line feature itself (line without the wings) ewindx = np.where((lmbda >= wavelength2) & (lmbda <= wavelength3)) + + # Sum continuosly over the normalized flux array and assign result to a variable. eqwi = sum(1-nflux[ewindx])*dlmbda + + # Calculate the uncertainty in the equivalent with measurement + eqwi1, err = mcerror(lmbda,flux,wavelength2,wavelength3,m,c,rmse) + + # Return the result of the equivalent width measurement and associated error + return eqwi, err + + +def mcerror(lmbda,flux,wavelength2,wavelength3,m,c,rmse): + """ + This function calculates the error of an equivalent width measurement. The error is characterized + by the value of the RMSE above or below the baseline generated for the feature. The error is calculated + via the mean of multiple error values generated via 3000 iterations of the monte carlo method for the + single emission feature. - # For error - #fx1, err = mcerr(wav,dat,w2,w3,m,c,rms) + ========== + ARGUMENTS: + ========== - return eqwi + lmbda - Wavelength array of entire feature (wings + line). + + flux - Normalized flux array of the entire feature (wings + line) + + wavelength2 - wavelength value where continuum region to the left of the spectral feature ends. + + wavelength3 - wavelength value where continuum region to the right of the spectral feature begins. + + m,c - Coefficients of a first degree polynomial fit to the wings of the entire line feature. + + rmse - RMSE error calculated over the first degree polynomial line fit to the wings(baseline). + + ======== + RETURNS: + ======== + + error and standard deviation of the equivalent width measurement taken over 3000 monte carlo iterations + + """ + + size = 3000 + guess= [] + dlmbda = lmbda[2]-lmbda[1] + + for i in range(size): + fitl = m*lmbda + c + np.random.uniform(-1,1)*rms + nflux = flux / fitl + ewindx = (lmbda >= wavelength2) & (lmbda <= wavelength3) + fx = sum(1-nflux[ewindx]) *dlmbda + guess.append(fx) + + return np.mean(guess),np.std(guess) #Calculate line height of spectral emission feature def line_height(filename,wavelength1,wavelength2,wavelength3,wavelength4): + """ + Tnis function calculates the emission line feature height as the vertical-distance between the + maximum flux value and minimum flux value over the feature interval. + + ========== + ARGUMENTS: + ========== + + filename - SDSS spectra .fits filename entered as a string. You may specify entire + path as a string, or if file is located in the current working directory, + simply the filename (e.g. 'SDSS-12345.fits') + + wavelength1 - wavelength value where continuum region to the left of the spectral feature begins. + + wavelength2 - wavelength value where continuum region to the left of the spectral feature ends. + + wavelength3 - wavelength value where continuum region to the right of the spectral feature begins. + + wavelength4 - wavelength value where continuum region to the right of the spectral feature ends. + + ======== + RETURNS: + ======== + + lnheight - line height of spectral emission feature. + + """ # Open sloan .fits file hdu = fits.open(filename) @@ -98,27 +174,40 @@ def line_height(filename,wavelength1,wavelength2,wavelength3,wavelength4): lmbda = 10**np.array(hdu[1].data['loglam']) flux = hdu[1].data['flux'] - # Normalize flux by dividing flux array by mean flux value over relatively quiet - # continuum region - indx = np.where((lmbda > 8250) & (lmbda < 8350)) - avgval = flux[indx].mean() - nflux = flux/avgval - - # designate fundamental and continuum indices - fuind = (lmbda >= wavelength1) & (lmbda <= wavelength4) - coind = ((lmbda >= wavelength1) & (lmbda <= wavelength2)) | ((lmbda >= wavelength3) & (lmbda <= wavelength4)) + # Index for continuum region where entire line feature is found (wings + emission line) + fuindx = (lmbda >= wavelength1) & (lmbda <= wavelength4) + + # Index for both of the wings of the line feature + coindx = ((lmbda >= wavelength1) & (lmbda <= wavelength2)) | ((lmbda >= wavelength3) & (lmbda <= wavelength4)) - lmbda = lmbda[fuind] - flux = nflux[fuind] - x1 = lmbda[coind] - y1 = nflux[coind] + # Isolate the wavelength subarrays pertaining to both of the wings of the line feature + x1 = lmbda[coindx] + + #Isolate the normalized flux subarrays of the line features wings. + y1 = nflux[coindx] + + # Generate first degree polynomial coefficients that correspond to the line wings m,c = np.polyfit(x1,y1,1) #print m,c - rms = sqrt( mean( (m*x1 + c - y1)**2 ) ) + # calculate the RMSE for the baseline for later use to generate error array. + rmse = sqrt( mean( (m*x1 + c - y1)**2 ) ) + + # Assign the wavelength and normalized flux arrays over the entire feature (wings + line) to new variables + lmbda = lmbda[fuindx] + flux = flux[fuindx] + + # using the coefficients of the deg(1) polynomial generated earlier from the wings, fit a line to them + # in order to establish the baseline for the line feature. + fitl = np.array(m*lmbda + c) + + # Normalize the flux array by dividing it by the fitline. + nflux = np.array(flux/fitline) - lnheight = max(y1)-min(y1) + ewindx = np.where((lmbda >= wavelength2) & (lmbda <= wavelength3)) + lnheight = max(nflux[ewindx])-min(nflux[ewindx]) + return lnheight From 8b34cd847e8fb75a9f7d5c3cacf1d76818746a7c Mon Sep 17 00:00:00 2001 From: JPVentura Date: Tue, 23 May 2017 23:12:16 -0400 Subject: [PATCH 3/5] added comments to the mcerror function --- ewidth.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/ewidth.py b/ewidth.py index d5deb04..20cb752 100644 --- a/ewidth.py +++ b/ewidth.py @@ -122,13 +122,19 @@ def mcerror(lmbda,flux,wavelength2,wavelength3,m,c,rmse): error and standard deviation of the equivalent width measurement taken over 3000 monte carlo iterations """ - + # assign the number of monte carlo iterations to variable 'size'. size = 3000 + + # Create an empty list to which eqmeasurements will be appended. guess= [] + + # Define the differential wavelength element and assign to variable 'dlmbda'. dlmbda = lmbda[2]-lmbda[1] + # Iterate error calculations over equivalent width measurement 3000 times by drawing randomly from + # a uniform distribution about the RMSE. Then take the mean and std dev of those values and return values. for i in range(size): - fitl = m*lmbda + c + np.random.uniform(-1,1)*rms + fitl = m*lmbda + c + np.random.uniform(-1,1)*rmse nflux = flux / fitl ewindx = (lmbda >= wavelength2) & (lmbda <= wavelength3) fx = sum(1-nflux[ewindx]) *dlmbda @@ -204,8 +210,10 @@ def line_height(filename,wavelength1,wavelength2,wavelength3,wavelength4): # Normalize the flux array by dividing it by the fitline. nflux = np.array(flux/fitline) + # Create index for the wavelength of the line feature itself (line without the wings) ewindx = np.where((lmbda >= wavelength2) & (lmbda <= wavelength3)) + # Calculate the line height for the featrue and assign to a variable. lnheight = max(nflux[ewindx])-min(nflux[ewindx]) return lnheight From 9dc88e83dab2166d2827177af68f4604a26bd1b2 Mon Sep 17 00:00:00 2001 From: JPVentura Date: Wed, 24 May 2017 15:28:29 -0400 Subject: [PATCH 4/5] made edits to README.md file --- README.md | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index 690a80c..5ae6d11 100644 --- a/README.md +++ b/README.md @@ -1,10 +1,6 @@ -# EqWi (Equivalent Widths): +# EqWiS (Equivalent Widths): -Python script to calculate the equivalent widths of spectral line emission & absorption features for various spectral data file types using the pyspeckit module [PyspecKit](https://github.com/pyspeckit/pyspeckit). By Jean-Paul Ventura +Python module to calculate the equivalent widths of spectral line emission features, emission feature line height, and uncertainties via the Monte Carlo Method of Sloan Digital Sky Survey (SDSS) Spectra. By Jean-Paul Ventura -### Information ### -### Referencing ### - -### License ### From ca97967e6812ad28cd6149326f8fa1fc94365df1 Mon Sep 17 00:00:00 2001 From: Jean-paul Ventura Date: Wed, 24 May 2017 15:32:58 -0400 Subject: [PATCH 5/5] Delete EqWi.py Deleted preliminary test code. --- EqWi.py | 40 ---------------------------------------- 1 file changed, 40 deletions(-) delete mode 100644 EqWi.py diff --git a/EqWi.py b/EqWi.py deleted file mode 100644 index 04d7256..0000000 --- a/EqWi.py +++ /dev/null @@ -1,40 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as numpy -from astropy.io import fits - -datfile = fits.open('filename') -hdr = datfile[1].header - - - -try: - data = datfile[1].data - except IndexError: - pass - else: - good = False - names = data.dtype.names - if 'loglam' in names and 'flux' in names: - flux = datfile[1].data['flux'] - wlen = 10**np.array(data.loglam) - indx = np.where((wlen > 8250) & (wlen < 8350)) - avgval = flux[indx].mean() - nflux = flux/avgval - - good = True - #elif 'wavelength' in names and 'flux' in names: - wa = data.wavelength - fl = data.flux - good = True - #if good: - er = np.ones_like(fl) - try: - er = data.er - except AttributeError: - pass - co = np.empty_like(fl) * np.nan - try: - co = data.co - except AttributeError: - pass - return Spectrum(wa=wa, fl=fl, er=er, co=co, filename=filename)