diff --git a/EqWi.py b/EqWi.py deleted file mode 100644 index 663ba12..0000000 --- a/EqWi.py +++ /dev/null @@ -1,42 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np -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 - if 'loglam' in names and 'flux' in names: - flux = datfile[1].data['flux'] - wlen = 10**np.array(datflike[1].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(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) -# diff --git a/README.md b/README.md index b8d0369..5ae6d11 100644 --- a/README.md +++ b/README.md @@ -1,10 +1,6 @@ -# Equivalent_widths +# EqWiS (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 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 -Created on February 16, 2017 -### Referencing -### License diff --git a/ewidth.py b/ewidth.py index 8203f65..20cb752 100644 --- a/ewidth.py +++ b/ewidth.py @@ -1,4 +1,4 @@ - + # Jean-paul Ventura # January 12, 2017 # Modified: February 20th, 2016 @@ -6,177 +6,216 @@ 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: + ================= + + eqwi - 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. -================= -Function returns: -================= + # 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)) + + # 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 + + # 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 = flux/fitline -The measured equivalent width of the spectral emission/absorption feature -""" + # 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. + + ========== + ARGUMENTS: + ========== + + 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 + + """ + # 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)*rmse + 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) - lambda = 10**np.array(hdu[1].data['loglam']) + #Construct wavelength (lmbda) and flux arrays from data + lmbda = 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)) + # 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)) - x1 = lambda[coind] - y1 = flux[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 ) ) - indx = np.where((lambda > 8250) & (lambda < 8350)) - avgval = flux[indx].mean() - nflux = flux/avgval - lambda = lambda[fuind] - flux = flux[fuind] - - - dlambda= lambda[2]-lambda[1] - fitl = m*lambda + c - - - - ewind = np.where((lambda >= w2) & (lambda <= w3)) - eqwi = sum(1-nflux[ewind])*dlambda - - # 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() + #print m,c + + # 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) + + # 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 + +