Skip to content

Commit

Permalink
remove shortlist calculation - different pr and not ready yet
Browse files Browse the repository at this point in the history
  • Loading branch information
@jvshields
jvshields committed Dec 7, 2023
1 parent 48dfa6f commit d32976a
Show file tree
Hide file tree
Showing 2 changed files with 15 additions and 309 deletions.
7 changes: 4 additions & 3 deletions stardis/base.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,9 +60,10 @@ def run_stardis(config_fname, tracing_lambdas_or_nus):
from stardis.io.model.mesa import read_mesa_model

raw_mesa_model = read_mesa_model(config.model.fname)
stellar_model = raw_mesa_model.to_stellar_model(
adata, truncate_to_shell_number=config.model.truncate_to_shell
)
if config.model.truncate_to_shell is not -99:
raw_mesa_model.truncate_model(config.model.truncate_to_shell)
else:
stellar_model = raw_mesa_model.to_stellar_model(adata)

else:
raise ValueError("Model type not recognized. Must be either 'marcs' or 'mesa'")
Expand Down
317 changes: 11 additions & 306 deletions stardis/plasma/base.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ class HMinusDensity(ProcessingPlasmaProperty):
Attributes
----------
h_minus_density : Pandas DataFrame, dtype float
Density of H-, indexed by depth point.
Density of H-, indexed by shell.
"""

outputs = ("h_minus_density",)
Expand All @@ -62,11 +62,11 @@ def calculate(self, ion_number_density, t_rad, electron_densities):
class H2Density(ProcessingPlasmaProperty):
"""
Used Kittel and Kroemer "Thermal Physics".
Attributes
----------
h2_density : Pandas DataFrame, dtype float
Density of H2, indexed by depth point.
Density of H2, indexed by shell.
"""

outputs = ("h2_density",)
Expand Down Expand Up @@ -129,302 +129,9 @@ def calculate(
return df


class AlphaLineVald(ProcessingPlasmaProperty):
"""
Attributes
----------
alpha_line_from_linelist : DataFrame
A pandas DataFrame with dtype float. This represents the alpha calculation
for each line from Vald at each depth point. Refer to Rybicki and Lightman
equation 1.80. Voigt profiles are calculated later, and B_12 is substituted
appropriately out for f_lu. This assumes LTE for lower level population.
"""

outputs = ("alpha_line_from_linelist", "lines_from_linelist")
latex_name = (r"\alpha_{\textrm{line, vald}}",)
latex_formula = (
r"\dfrac{\pi e^{2} n_{lower} f_{lu}}{m_{e} c}\
\Big(1-exp(-h \nu / k T) \phi(\nu)\Big)",
)

def calculate(
self,
atomic_data,
ion_number_density,
t_electrons,
g,
ionization_data,
):
# solve n_lower : n * g_i / g_0 * e ^ (-E_i/kT)
# get f_lu : loggf -> use g = 2j+1
# emission_correction = (1-e^(-h*nu / kT))
# alphas = ALPHA_COEFFICIENT * n_lower * f_lu * emission_correction

points = len(t_electrons)

# Need degeneracy of ground state of the ion to calculate n_lower
# So set up the initial dataframe with all of the necessary information, and then merge it with the matching g_0
linelist = atomic_data.linelist.rename(columns={"ion_charge": "ion_number"})[
[
"atomic_number",
"ion_number",
"wavelength",
"log_gf",
"e_low",
"e_up",
"j_lo",
"j_up",
"rad",
"stark",
"waals",
]
].merge(
g.loc[:, :, 0].rename("g_0"),
how="left",
on=["atomic_number", "ion_number"],
)

# Truncate to final atomic number
linelist = linelist[
linelist.atomic_number <= (atomic_data.selected_atomic_numbers.max())
]

# Calculate degeneracies
linelist["g_lo"] = linelist.j_lo * 2 + 1
linelist["g_up"] = linelist.j_up * 2 + 1

exponent_by_point = np.exp(
np.outer(
-linelist.e_low.values * u.eV, 1 / (t_electrons * u.K * const.k_B)
).to(1)
)

# grab densities for n_lower - need to use linelist as the index
linelist_with_densities = linelist.merge(
ion_number_density,
how="left",
on=["atomic_number", "ion_number"],
)

n_lower = (
(
exponent_by_point * linelist_with_densities[np.arange(points)]
).values.T # arange mask of the dataframe returns the set of densities of the appropriate ion for the line at each point
* linelist_with_densities.g_lo.values
/ linelist_with_densities.g_0.values
)

linelist["f_lu"] = (
10**linelist.log_gf / linelist.g_lo
) # vald log gf is "oscillator strength f times the statistical weight g of the parent level" see 1995A&AS..112..525P, section 2. Structure of VALD

line_nus = (linelist.wavelength.values * u.AA).to(
u.Hz, equivalencies=u.spectral()
)

emission_correction = 1 - np.exp(
(
-const.h
/ const.k_B
* np.outer(
line_nus,
1 / (t_electrons * u.K),
)
).to(1)
)

alphas = pd.DataFrame(
(
ALPHA_COEFFICIENT
* n_lower
* linelist.f_lu.values
* emission_correction.T
).T
)

if np.any(np.isnan(alphas)) or np.any(np.isinf(np.abs(alphas))):
raise ValueError(
"Some alpha_line from vald are nan, inf, -inf " " Something went wrong!"
)

alphas["nu"] = line_nus.value
linelist["nu"] = line_nus.value

# Linelist preparation below is taken from opacities_solvers/base/calc_alpha_line_at_nu
# Necessary for correct handling of ion numbers using charge instead of astronomy convention (i.e., 0 is neutral, 1 is singly ionized, etc.)
ionization_energies = ionization_data.reset_index()
ionization_energies["ion_number"] -= 1
linelist = pd.merge(
linelist,
ionization_energies,
how="left",
on=["atomic_number", "ion_number"],
)

linelist["level_energy_lower"] = ((linelist["e_low"].values * u.eV).cgs).value
linelist["level_energy_upper"] = ((linelist["e_up"].values * u.eV).cgs).value

# Radiation broadening parameter is approximated as the einstein A coefficient. Vald parameters are in log scale.
linelist["A_ul"] = 10 ** (
linelist["rad"]
) # see 1995A&AS..112..525P for appropriate units - may be off by a factor of 4pi

# Need to remove autoionization lines - can't handle with current broadening treatment because can't calculate effective principal quantum number
valid_indices = linelist.level_energy_upper < linelist.ionization_energy

return alphas[valid_indices], linelist[valid_indices]


# Properties that haven't been used in creating stellar plasma yet,
# might be useful in future ----------------------------------------------------

class AlphaLineShortlistVald(ProcessingPlasmaProperty):
"""
Attributes
----------
alpha_line_from_linelist : DataFrame
A pandas DataFrame with dtype float. This represents the alpha calculation
for each line from Vald at each depth point. Refer to Rybicki and Lightman
equation 1.80. Voigt profiles are calculated later, and B_12 is substituted
appropriately out for f_lu. This assumes LTE for lower level population.
This is the same as AlphaLineVald, but does not recalculate the lower level population because degeneracies are not available for the shortlist.
"""

outputs = ("alpha_line_from_linelist", "lines_from_linelist")
latex_name = (r"\alpha_{\textrm{line, vald}}",)
latex_formula = (
r"\dfrac{\pi e^{2} n_{lower} f_{lu}}{m_{e} c}\
\Big(1-exp(-h \nu / k T) \phi(\nu)\Big)",
)

def calculate(
self,
atomic_data,
ion_number_density,
t_electrons,
g,
ionization_data,
):
# solve n_lower : n * g_i / g_0 * e ^ (-E_i/kT)
# get f_lu : loggf -> use g = 2j+1
# emission_correction = (1-e^(-h*nu / kT))
# alphas = ALPHA_COEFFICIENT * n_lower * f_lu * emission_correction

points = len(t_electrons)

# Need degeneracy of ground state of the ion to calculate n_lower
# So set up the initial dataframe with all of the necessary information, and then merge it with the matching g_0
linelist = atomic_data.linelist.rename(columns={"ion_charge": "ion_number"})[
[
"atomic_number",
"ion_number",
"wavelength",
"log_gf",
"e_low",
"e_up",
"j_lo",
"j_up",
"rad",
"stark",
"waals",
]
].merge(
g.loc[:, :, 0].rename("g_0"),
how="left",
on=["atomic_number", "ion_number"],
)

# Truncate to final atomic number
linelist = linelist[
linelist.atomic_number <= (atomic_data.selected_atomic_numbers.max())
]

# Calculate degeneracies
linelist["g_lo"] = linelist.j_lo * 2 + 1
linelist["g_up"] = linelist.j_up * 2 + 1

exponent_by_point = np.exp(
np.outer(
-linelist.e_low.values * u.eV, 1 / (t_electrons * u.K * const.k_B)
).to(1)
)

# grab densities for n_lower - need to use linelist as the index
linelist_with_densities = linelist.merge(
ion_number_density,
how="left",
on=["atomic_number", "ion_number"],
)

n_lower = (
(
exponent_by_point * linelist_with_densities[np.arange(points)]
).values.T # arange mask of the dataframe returns the set of densities of the appropriate ion for the line at each point
* linelist_with_densities.g_lo.values
/ linelist_with_densities.g_0.values
)

linelist["f_lu"] = (
10**linelist.log_gf / linelist.g_lo
) # vald log gf is "oscillator strength f times the statistical weight g of the parent level" see 1995A&AS..112..525P, section 2. Structure of VALD

line_nus = (linelist.wavelength.values * u.AA).to(
u.Hz, equivalencies=u.spectral()
)

emission_correction = 1 - np.exp(
(
-const.h
/ const.k_B
* np.outer(
line_nus,
1 / (t_electrons * u.K),
)
).to(1)
)

alphas = pd.DataFrame(
(
ALPHA_COEFFICIENT
* n_lower
* linelist.f_lu.values
* emission_correction.T
).T
)

if np.any(np.isnan(alphas)) or np.any(np.isinf(np.abs(alphas))):
raise ValueError(
"Some alpha_line from vald are nan, inf, -inf " " Something went wrong!"
)

alphas["nu"] = line_nus.value
linelist["nu"] = line_nus.value

# Linelist preparation below is taken from opacities_solvers/base/calc_alpha_line_at_nu
# Necessary for correct handling of ion numbers using charge instead of astronomy convention (i.e., 0 is neutral, 1 is singly ionized, etc.)
ionization_energies = ionization_data.reset_index()
ionization_energies["ion_number"] -= 1
linelist = pd.merge(
linelist,
ionization_energies,
how="left",
on=["atomic_number", "ion_number"],
)

linelist["level_energy_lower"] = ((linelist["e_low"].values * u.eV).cgs).value
linelist["level_energy_upper"] = ((linelist["e_up"].values * u.eV).cgs).value

# Radiation broadening parameter is approximated as the einstein A coefficient. Vald parameters are in log scale.
linelist["A_ul"] = 10 ** (
linelist["rad"]
) # see 1995A&AS..112..525P for appropriate units - may be off by a factor of 4pi

# Need to remove autoionization lines - can't handle with current broadening treatment because can't calculate effective principal quantum number
valid_indices = linelist.level_energy_upper < linelist.ionization_energy

return alphas[valid_indices], linelist[valid_indices]


class InputNumberDensity(DataFrameInput):
"""
Expand Down Expand Up @@ -455,15 +162,14 @@ def calculate(self, number_density):
# Creating stellar plasma ------------------------------------------------------


def create_stellar_plasma(stellar_model, atom_data, config):
def create_stellar_plasma(stellar_model, atom_data):
"""
Creates stellar plasma.
Parameters
----------
stellar_model : stardis.model.base.StellarModel
atom_data : tardis.io.atom_data.base.AtomData
config : stardis.config_reader.Configuration
Returns
-------
Expand All @@ -490,24 +196,23 @@ def create_stellar_plasma(stellar_model, atom_data, config):
)
plasma_modules += helium_lte_properties

plasma_modules.append(AlphaLine)

plasma_modules.append(HMinusDensity)
plasma_modules.append(H2Density)

if config.opacity.line.use_vald_linelist:
plasma_modules.append(AlphaLineVald)
else:
plasma_modules.append(AlphaLine)

# plasma_modules.remove(tardis.plasma.properties.radiative_properties.StimulatedEmissionFactor)
# plasma_modules.remove(tardis.plasma.properties.general.SelectedAtoms)
# plasma_modules.remove(tardis.plasma.properties.plasma_input.Density)

fv_geometry = stellar_model.fv_geometry

return BasePlasma(
plasma_properties=plasma_modules,
t_rad=stellar_model.temperatures.value,
abundance=stellar_model.composition.atomic_mass_fraction,
t_rad=fv_geometry.t.values,
abundance=stellar_model.abundances,
atomic_data=atom_data,
density=stellar_model.composition.density.value,
density=fv_geometry.density.values,
link_t_rad_t_electron=1.0,
nlte_ionization_species=[],
nlte_excitation_species=[],
Expand Down

0 comments on commit d32976a

Please sign in to comment.