Skip to content

Commit

Permalink
Address merge conflicts for PR #23
Browse files Browse the repository at this point in the history
  • Loading branch information
@astrochun
astrochun committed May 30, 2020
2 parents 25622d0 + c2e443f commit 677472c
Show file tree
Hide file tree
Showing 2 changed files with 132 additions and 77 deletions.
38 changes: 30 additions & 8 deletions Metallicity_Stack_Commons/analysis/attenuation.py
View file
Original file line number Diff line number Diff line change
@@ -1,36 +1,58 @@
from astropy.io import ascii as asc
from astropy.table import Table
import numpy as np
from os.path import join

from .. import k_dict

HgHb_CaseB = 0.468 # Hg/Hb ratio for zero reddening

HgHb_CaseB = 0.468 # Hg/Hb ratio for zero reddening
HaHb_CaseB = 2.86 # Ha/Hb ratio for zero reddening

k_HBETA = k_dict['HBETA']
k_HGAMMA = k_dict['HGAMMA']


def compute_EBV(fitspath, combine_asc):
"""
Purpose:
Determines E(B-V) from Hg/Hb flux ratio
:param fitspath: str containing root path
:param combine_asc: Astropy table containing emission-line flux
"""

ID = combine_asc['ID']
HBETA = combine_asc['HBETA_Flux_Observed'].data
HGAMMA = combine_asc['HGAMMA_Flux_Observed'].data

EBV = -2.5 * np.log10((HGAMMA / HBETA) / (HgHb_CaseB)) / (k_HGAMMA - k_HBETA)
EBV = -2.5 * np.log10((HGAMMA / HBETA) / HgHb_CaseB) / (k_HGAMMA - k_HBETA)

out_ascii = fitspath + '/dust_attenuation_values.tbl'
out_ascii = join(fitspath, 'dust_attenuation_values.tbl')

tab1 = Table([ID, EBV], names=('ID', 'E(B-V)'))
asc.write(tab1, out_ascii, format='fixed_width_two_line')


def compute_A(EBV):
"""
Purpose:
Compute A(Lambda) for all possible emission lines
:param EBV: float value of E(B-V)
Has not been configured to handle a large array. Some array handling would be needed
:return A_dict: dict containing A(lambda) with keys identical to k_dict
"""

k_arr = np.array(list(k_dict.values()))

A_arr = k_arr * EBV
A_dict = dict(zip(list(k_dict.keys()),A_arr))
A_dict = dict(zip(list(k_dict.keys()), A_arr))

return A_dict


def line_ratio_atten(ratio, EBV, wave_top, wave_bottom):

k_top = k_dict[wave_top]
Expand All @@ -40,9 +62,9 @@ def line_ratio_atten(ratio, EBV, wave_top, wave_bottom):

return ratio_atten

def Hb_SFR(log_LHb, EBV):

'''
def Hb_SFR(log_LHb, EBV):
"""
Purpose:
Determine dust-corrected SFR using the H-beta luminosity and a
measurement for nebular attenuation
Expand All @@ -56,8 +78,8 @@ def Hb_SFR(log_LHb, EBV):
:return logSFR: numpy array or float containing the SFR in
logarithmic units of M_sun/yr
'''
"""

logSFR = np.log10(4.4e-42 * 2.86) + 0.4*EBV*k_dict['HBETA'] + log_LHb
logSFR = np.log10(4.4e-42 * HaHb_CaseB) + 0.4*EBV*k_HBETA + log_LHb

return logSFR
171 changes: 102 additions & 69 deletions Metallicity_Stack_Commons/plotting/balmer.py
View file
Original file line number Diff line number Diff line change
@@ -1,56 +1,59 @@
'''
"""
balmer
------
Generates plots illustrating Balmer recombination lines
This code was created from:
https://github.com/astrochun/Zcalbase_gal/blob/master/Analysis/DEEP2_R23_O32/balmer_plots.py
'''

"""

import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.io import ascii as asc
from matplotlib.backends.backend_pdf import PdfPages

# This needs to be updated with new definitions and assignments
from ..analysis.fitting import gauss, double_gauss

# This needs to be updated with new definitions and assignments
from .. import scalefact, wavelength_dict

nrows = 3
ncols = 3
n_rows = 3
n_cols = 3


def extract_fit(astropy_table, line_name, balmer=False):
'''
"""
Purpose:
Extract best fit from table and returns a dictionary
Extract best fit from table, return a list of fitting parameters
:param astropy_table: Astropy table
:param astropy_table: Astropy table containing fitting result
:param line_name: line to extract fit results
:param balmer: boolean to indicate whether line is a Balmer line
:return:
param_list: list containing fit to pass in
'''
"""

try:
xbar = combine_asc[line_name + '_X_bar']
astropy_table[line_name + '_Center'].data
except KeyError:
print("Line not present in table")
print("Exiting!!!")
return

xbar = combine_asc[line_name + '_X_bar']
sp = combine_asc[line_name + '_Pos_Sig']
ap = combine_asc[line_name + '_Pos_Amp']
con = combine_asc[line_name + '_Const']
xbar = astropy_table[line_name + '_Center'].data
sp = astropy_table[line_name + '_Sigma'].data
ap = astropy_table[line_name + '_Norm'].data
con = astropy_table[line_name + '_Median'].data

param_list = [xbar, sp, ap, con]

if balmer:
sn = combine_asc[line_name + '_Neg_Sig']
an = combine_asc[line_name + '_Neg_Amp']
sn = astropy_table[line_name + '_Abs_Sigma'].data
an = astropy_table[line_name + '_Abs_Norm'].data

param_list += [sn, an]

Expand All @@ -60,12 +63,27 @@ def extract_fit(astropy_table, line_name, balmer=False):
else:
return param_list


def fitting_result(wave, y_norm, lambda_cen, balmer_fit, balmer_fit_neg):
"""
Purpose:
Returns fitting results based on inputs of best fit
:param wave: numpy array of rest-frame wavelength
:param y_norm: Normalize 1-D spectra in units of 10^-17 erg/s/cm2/AA
:param lambda_cen: Central wavelength in Angstroms
:param balmer_fit: list containing Balmer emission fits
:param balmer_fit_neg: list containing the absorption ("stellar") Balmer fit
:return:
"""

dx = wave[2] - wave[1]

x_sigsnip = np.where(np.abs((wave - lambda_cen)) / balmer_fit[1] <= 2.5)[0]
gauss0 = double_gauss(wave, *balmer_fit)
neg0 = gauss(wave, *balmer_fit_neg)
gauss0_diff = gauss0 - neg0
y_norm_diff = y_norm[x_sigsnip] - Bneg0[x_sigsnip]
y_norm_diff = y_norm[x_sigsnip] - neg0[x_sigsnip]

# Residuals
x_sigsnip_2 = np.where(np.abs((wave - lambda_cen)) / balmer_fit[1] <= 3.0)[0]
Expand All @@ -77,96 +95,111 @@ def fitting_result(wave, y_norm, lambda_cen, balmer_fit, balmer_fit_neg):

return gauss0, resid, x_sigsnip_2, flux_g, flux_s

def HbHgHd_fits(fitspath, nrow, ncol,Stack_name,combine_flux_tab, out_pdf):

stack2D, header = fits.getdata(Stack_name,header=True)
# noinspection PyUnboundLocalVariable
def HbHgHd_fits(stack_name, astropy_table_file, out_pdf):
"""
Purpose:
Generate PDF plots that illustrate H-delta, H-gamma, and H-beta line
profiles and best fit
:param stack_name: filename of the stack spectra (str)
:param astropy_table_file: astropy table containing ID
:param out_pdf: Name of output file (str)
"""

stack2D, header = fits.getdata(stack_name, header=True)
wave = header['CRVAL1'] + header['CDELT1']*np.arange(header['NAXIS1'])
dispersion = header['CDELT1']

combine_asc = asc.read(combine_flux_tab)
astropy_table = asc.read(astropy_table_file)

ID = combine_asc['ID']
ID = astropy_table['ID'].data

pdfpages = PdfPages(out_pdf)
pdf_pages = PdfPages(out_pdf)

for ii in range(len(ID)):

if ii % nrows ==0:
fig, ax_arr = plt.subplots(nrows=nrows, ncols=ncols, squeeze = False)
if ii % n_rows == 0:
fig, ax_arr = plt.subplots(nrows=n_rows, ncols=n_cols, squeeze=False)

y0 = stack2D[ii]
y_norm = y0/scalefact
dx = wave[2]-wave[1]

Hb_fit, Hb_fit_neg = extract_fit(combine_asc[ii], 'HBETA', balmer=True)
Hg_fit, Hg_fit_neg = extract_fit(combine_asc[ii], 'Hgamma', balmer=True)
Hd_fit, Hd_fit_neg = extract_fit(combine_asc[ii], 'HDELTA', balmer=True)
Hb_fit, Hb_fit_neg = extract_fit(astropy_table[ii], 'HBETA', balmer=True)
Hg_fit, Hg_fit_neg = extract_fit(astropy_table[ii], 'HGAMMA', balmer=True)
Hd_fit, Hd_fit_neg = extract_fit(astropy_table[ii], 'HDELTA', balmer=True)

# This will need to be updated
wave_beta = wavelength_dict['HBETA']
wave_gamma = wavelength_dict['HGAMMA']
wave_delta = wavelength_dict['HDELTA']

##Beta
fit_result_in = [wave, y_norm, wavelength_dict['HBETA'], Hb_fit, Hb_fit_neg]
# Beta
fit_result_in = [wave, y_norm, wave_beta, Hb_fit, Hb_fit_neg]
Bgauss0, Bresid, Bx_sigsnip_2, Bflux_g, Bflux_s = fitting_result(fit_result_in)

##Gamma
fit_result_in = [wave, y_norm, wavelength_dict['HGAMMA'], Hg_fit, Hg_fit_neg]
# Gamma
fit_result_in = [wave, y_norm, wave_gamma, Hg_fit, Hg_fit_neg]
Ggauss0, Gresid, Gx_sigsnip_2, Gflux_g, Gflux_s = fitting_result(fit_result_in)

##Delta
fit_result_in = [wave, y_norm, wavelength_dict['HDELTA'], Hd_fit, Hd_fit_neg]
# Delta
fit_result_in = [wave, y_norm, wave_delta, Hd_fit, Hd_fit_neg]
Dgauss0, Dresid, Dx_sigsnip_2, Dflux_g, Dflux_s = fitting_result(fit_result_in)

row = ii % nrows
row = ii % n_rows

txt0 = r'ID: %i' % (ID[ii]) +'\n'
txt0 += r'+$\sigma$: %.3f, -$\sigma$: %.3f '% (Hb_fit[1], Hb_fit_neg[1]) + '\n'
txt0 += 'F_G: %.3f F_S: %.3f' %(Bflux_g, Bflux_s)
# The below code could be refactored or simplified
txt0 = r'ID: %i' % (ID[ii]) + '\n'
txt0 += r'+$\sigma$: %.3f, -$\sigma$: %.3f ' % (Hb_fit[1], Hb_fit_neg[1]) + '\n'
txt0 += 'F_G: %.3f F_S: %.3f' % (Bflux_g, Bflux_s)

ax_arr[row][2].plot(wave, y_norm,'k', linewidth=0.3, label= 'Emission')
ax_arr[row][2].plot(wave,Bgauss0, 'm', linewidth= 0.25, label= 'Beta Fit')
ax_arr[row][2].plot(wave, y_norm, 'k', linewidth=0.3, label='Emission')
ax_arr[row][2].plot(wave, Bgauss0, 'm', linewidth=0.25, label='Beta Fit')
ax_arr[row][2].set_xlim(wave_beta-50, wave_beta+50)

ax_arr[row][2].annotate(txt0, [0.95,0.95], xycoords='axes fraction',
va='top', ha='right', fontsize= '5')
ax_arr[row][2].plot(wave[Bx_sigsnip_2],Bresid, 'r', linestyle='dashed',
linewidth = 0.2, label= 'Residuals')
ax_arr[row][2].annotate(txt0, [0.95, 0.95], xycoords='axes fraction',
va='top', ha='right', fontsize='5')
ax_arr[row][2].plot(wave[Bx_sigsnip_2], Bresid, 'r', linestyle='dashed',
linewidth=0.2, label='Residuals')

txt1 = r'ID: %i' % (ID[ii]) +'\n'
txt1 += r'+$\sigma$: %.3f, -$\sigma$: %.3f '% (Hg_fit[1], Hg_fit_neg[1]) + '\n'
txt1 += 'F_G: %.3f F_S: %.3f' %(Gflux_g, Gflux_s)
txt1 = r'ID: %i' % (ID[ii]) + '\n'
txt1 += r'+$\sigma$: %.3f, -$\sigma$: %.3f ' % (Hg_fit[1], Hg_fit_neg[1]) + '\n'
txt1 += 'F_G: %.3f F_S: %.3f' % (Gflux_g, Gflux_s)

ax_arr[row][1].plot(wave, y_norm,'k', linewidth=0.3, label= 'Emission')
ax_arr[row][1].plot(wave,Ggauss0, 'm', linewidth= 0.25, label= 'Gamma Fit')
ax_arr[row][1].plot(wave, y_norm, 'k', linewidth=0.3, label='Emission')
ax_arr[row][1].plot(wave, Ggauss0, 'm', linewidth=0.25, label='Gamma Fit')
ax_arr[row][1].set_xlim(wave_gamma-50, wave_gamma+50)

ax_arr[row][1].annotate(txt1, [0.95,0.95], xycoords='axes fraction',
va='top', ha='right', fontsize= '5')
ax_arr[row][1].plot(wave[Gx_sigsnip_2],Gresid, 'r', linestyle='dashed',
linewidth = 0.2, label= 'Residuals')
ax_arr[row][1].annotate(txt1, [0.95, 0.95], xycoords='axes fraction',
va='top', ha='right', fontsize='5')
ax_arr[row][1].plot(wave[Gx_sigsnip_2], Gresid, 'r', linestyle='dashed',
linewidth=0.2, label='Residuals')

txt2 = r'ID: %i' % (ID[ii]) +'\n'
txt2 += r'+$\sigma$: %.3f, -$\sigma$: %.3f '% (Hd_fit[1], Hd_fit_neg[1]) + '\n'
txt2 += 'F_G: %.3f F_S: %.3f' %(Dflux_g, Dflux_s)
txt2 = r'ID: %i' % (ID[ii]) + '\n'
txt2 += r'+$\sigma$: %.3f, -$\sigma$: %.3f ' % (Hd_fit[1], Hd_fit_neg[1]) + '\n'
txt2 += 'F_G: %.3f F_S: %.3f' % (Dflux_g, Dflux_s)

ax_arr[row][0].plot(wave, y_norm,'k', linewidth=0.3, label= 'Emission')
ax_arr[row][0].plot(wave,Dgauss0, 'm', linewidth= 0.25, label= 'Detla Fit')
ax_arr[row][0].plot(wave, y_norm, 'k', linewidth=0.3, label='Emission')
ax_arr[row][0].plot(wave, Dgauss0, 'm', linewidth=0.25, label='Delta Fit')
ax_arr[row][0].set_xlim(wave_delta-50, wave_delta+50)

ax_arr[row][0].set_ylim(0,1.5)
ax_arr[row][0].set_ylim(0, 1.5)

ax_arr[row][0].annotate(txt0, [0.95,0.95], xycoords='axes fraction',
va='top', ha='right', fontsize= '5')
ax_arr[row][0].plot(wave[Dx_sigsnip_2],Dresid, 'r', linestyle='dashed',
linewidth = 0.2, label= 'Residuals')

ax_arr[row][0].annotate(txt0, [0.95, 0.95], xycoords='axes fraction',
va='top', ha='right', fontsize='5')
ax_arr[row][0].plot(wave[Dx_sigsnip_2], Dresid, 'r', linestyle='dashed',
linewidth=0.2, label='Residuals')

ax_arr[row][0].set_yticklabels([0,0.5,1,1.5])
ax_arr[row][0].set_yticklabels([0, 0.5, 1, 1.5])
ax_arr[row][1].set_yticklabels([])
ax_arr[row][2].set_yticklabels([])

if row != nrows-1 and ii != stack2D.shape[0]-1:
if row != n_rows-1 and ii != stack2D.shape[0]-1:
ax_arr[row][0].set_xticklabels([])
ax_arr[row][1].set_xticklabels([])
ax_arr[row][2].set_xticklabels([])

if ii % nrows == nrows-1: fig.savefig(pdfpages, format='pdf')
if ii % n_rows == n_rows-1:
fig.savefig(pdf_pages, format='pdf')

pdfpages.close()
pdf_pages.close()

0 comments on commit 677472c

Please sign in to comment.