From 581edbb148b5144a9a20acbd36b512cdef764ce1 Mon Sep 17 00:00:00 2001 From: Jaad Tannous Date: Wed, 12 Feb 2025 16:16:06 -0500 Subject: [PATCH] adjusted xsd directory --- lvlspy/calculate/evolve.py | 96 ------------- lvlspy/calculate/isomer.py | 260 ---------------------------------- lvlspy/calculate/weisskopf.py | 256 --------------------------------- 3 files changed, 612 deletions(-) delete mode 100644 lvlspy/calculate/evolve.py delete mode 100644 lvlspy/calculate/isomer.py delete mode 100644 lvlspy/calculate/weisskopf.py diff --git a/lvlspy/calculate/evolve.py b/lvlspy/calculate/evolve.py deleted file mode 100644 index 5abfaed..0000000 --- a/lvlspy/calculate/evolve.py +++ /dev/null @@ -1,96 +0,0 @@ -""" -Module to handle the evolution of the level system based on temperature -""" - -import numpy as np -from scipy.sparse import csc_matrix -from scipy.sparse.linalg import expm_multiply - - -def newton_raphson(sp, temp, y0, time, tol=1e-6): - """ - Evolves a system using the Newton-Raphson method - - Args: - ``sp`` (:obj:`lvlspy.species.Species`): The species containing the levels to be evolved - - ``temp`` (:obj:`float`): The temperature in K to evolve the system at. - - ``y0`` (:obj:`numpy.array`): 1D array containing the initial distribution. - - ``time`` (:obj:`numpy.array`,optional): An array containing the time steps. - ``tol`` (:obj:`float`): The convergence condition for the method. Defaults to 1e-6. - - Returns: - ``y`` (:obj:`numpy.array`): 2D array of size n_levels*n_time containing the evolved system. - - ``fug`` (:obj:`numpy.array`) 2D array containing the fugacities as a function of time - """ - - y = np.empty((len(y0), len(time))) - fug = np.empty((len(y0), len(time))) - y[:, 0] = y0 - - rm = sp.compute_rate_matrix( - temp - ) # calculate the rate matrix of the species - eq_prob = sp.compute_equilibrium_probabilities(temp) - y_dt = y[:, 0] - fug[:, 0] = y_dt / eq_prob - for i in range(1, len(time)): - dt = time[i] - time[i - 1] - matrix = np.identity(len(y_dt)) - dt * rm - delta = np.ones(len(y_dt)) - while max(delta) > tol: - delta = np.linalg.solve( - matrix, -_f_vector(y_dt, y[:, i - 1], matrix) - ) - y_dt = y_dt + delta - y[:, i] = y_dt - fug[:, i] = y[:, i] / eq_prob - - return y, fug - - -def _f_vector(y_dt, y_i, rm): - return np.matmul(rm, y_dt) - y_i - - -def csc(sp, temp, y0, time): - """Evolves a system using sparse solver - - Args: - ``sp`` (:obj:`lvlspy.species.Species`): The species containing the levels to be evolved - - ``temp`` (:obj:`float`): The temperature in K to evolve the system at. - - ``y0`` (:obj:`numpy.array`): Array containing the initial condition. - - ``time`` (:obj:`numpy.array`): An array containing the time stamps to evolve the system - - Returns: - ``sol_expm_solver`` (:obj:`numpy.array`): A 2D array containing the evolved system - - ``fug`` (:obj:`numpy.array`) 2D array containing the fugacities as a function of time - """ - - rm = sp.compute_rate_matrix(temp) - eq_prob = sp.compute_equilibrium_probabilities(temp) - rm_csc = csc_matrix(rm) - sol_expm_solver = np.empty([rm.shape[0], time.shape[0]]) - sol_expm_solver[:, 0] = y0 - fug = np.empty((len(y0), len(time))) - fug[:, 0] = y0 / eq_prob - for i in range(len(time) - 1): - y = expm_multiply( - rm_csc, - y0, - start=time[i], - stop=time[i + 1], - num=2, - endpoint=True, - )[0, :] - sol_expm_solver[:, i + 1] = y - fug[:, i + 1] = y / eq_prob - - return sol_expm_solver, fug diff --git a/lvlspy/calculate/isomer.py b/lvlspy/calculate/isomer.py deleted file mode 100644 index a40a3e0..0000000 --- a/lvlspy/calculate/isomer.py +++ /dev/null @@ -1,260 +0,0 @@ -""" -Module to handle isomer related calculations. Functions are built based on the method in -`Gupta and Meyer (2001) `_ -""" - -import numpy as np - - -def transfer_properties(rate_matrix, level_low, level_high): - """Method that calculatest the transfer properties based on the rate matrix - - Args: - ``rate_matrix`` (:obj:`numpy.array`) A 2D array containing the rate matrix - of a species at a given temperature - - ``level_low`` (:obj:`int`) Integer indicating the lower level the transition - is moving to - - ``level_high`` (:obj:`int`) Integer indicating the higher level the transition - is moving from - - Returns: - On successful return, an array containing the following variable is passed: - - ``tpm`` (:obj:`numpy.array`) A 2D array containing the Transition Probability Matrix - - ``f_low_in`` (:obj:`numpy.array`) An array containing the branching ratio from all the - upper levels into the low level - - ``f_low_out`` (:obj:`numpy.array`) An array containing the branching ratio from the lower - levels into all the other levels - - ``f_high_in`` (:obj:`numpy.array`) An array containing the branching ratio from all the - levels into the high level - - ``f_high_out`` (:obj:`numpy.array`) An containing the branching ration out of the high - level into all the other levels - - ``lambda_sum`` (:obj:`numpy.array`) An array containing the diagonal elements of - the rate matrix - - """ - - # setting in the rates going in to the levels - lambda_low_in = rate_matrix[level_low, :] - lambda_high_in = rate_matrix[level_high, :] - - # remove entries corresponding the rows and columns to be removed - lambda_low_in = np.delete(lambda_low_in, [level_low, level_high]) - lambda_high_in = np.delete(lambda_high_in, [level_low, level_high]) - - # setting the rates going out of the levels - lambda_low_out = rate_matrix[:, level_low] - lambda_high_out = rate_matrix[:, level_high] - - # remove entries corresponding to the rows and columns to be removed - lambda_low_out = np.delete(lambda_low_out, [level_low, level_high]) - lambda_high_out = np.delete(lambda_high_out, [level_low, level_high]) - - # extract the diagonal elements from the rate matrix as they are the - # sum of all the rates into the level - lambda_sum = np.diag(rate_matrix) - # this array is the reduced array above without the removed levels - lambda_red = np.delete(lambda_sum, [level_low, level_high]) - - f_low_out = lambda_low_out / lambda_sum[level_low] - f_high_out = lambda_high_out / lambda_sum[level_high] - - f_low_in = lambda_low_in / lambda_red - f_high_in = lambda_high_in / lambda_red - - # setting up the transfer matrix - tpm = rate_matrix - tpm = tpm.T - # remove the columns - tpm = np.delete(tpm, [level_low, level_high], axis=1) - # remove the rows - tpm = np.delete(tpm, [level_low, level_high], axis=0) - - # Divide the row by the diagonal term - tpm = ( - tpm / lambda_red[:, None] - ) # This only works if the arrays are numpy arrays - - # set the diagonal to 0 - np.fill_diagonal(tpm, 0.0) - - return [tpm, f_low_in, f_low_out, f_high_in, f_high_out, lambda_sum] - - -def effective_rate(t, sp, level_low=0, level_high=1): - """ - Method to calculate the effective transition rates between the isomeric and ground states - - Args: - ``t`` (:obj:`float`) The temperature in K. - - ``sp`` (:obj:`lvlspy.species.Species`) The species of which the level system belongs to. - - ``level_low`` (:obj:`int`, optional) The lower level the effective transition rates are - calculated to. Defaults to 0; the ground state. - - ``level_high`` (:obj:`int`, optional) The higher level the effective transtion rates are - calculated to. Defaults to 1; the first excited state - - Returns: - Upon successful return, the method returns the effective transition rates between the higher - and lower level at temperature T - - ``l_low_high`` (:obj:`float`) The effective transition rate from the lower level to the - higher level - - ``l_high_low`` (:obj:`float`) The effective transition rate from the higher level to the - lower level - """ - - rate_matrix = np.abs(sp.compute_rate_matrix(t)) - trans_props = transfer_properties(rate_matrix, level_low, level_high) - f_n = _partial_sum(trans_props[0]) - - # Lambda_21_eff - l_high_low = trans_props[5][level_high] * np.matmul( - trans_props[4].T, np.matmul(f_n, trans_props[1]) - ) - # Lambda_12_eff - l_low_high = trans_props[5][level_low] * ( - np.matmul( - trans_props[2].T, - np.matmul(f_n, trans_props[3]), - ) - ) - - return ( - l_low_high + rate_matrix[level_high, level_low], - l_high_low + rate_matrix[level_low, level_high], - ) - - -def _partial_sum(tpm): - n_terms = 50000 - f_n = np.identity(tpm.shape[0]) + tpm - f_p = tpm - i = 2 - while i < n_terms: - f_p = np.matmul(f_p, tpm) - f_n += f_p - i += 1 - if np.linalg.norm(f_n, np.inf) < 1e-6: - break - - return f_n - - -def cascade_probabilities(t, sp, level_low=0, level_high=1): - """ - Method to calculate the cascace probability vectores (gammas) - - Args: - ``t`` (:obj:`float`) The temperature in K - - ``sp`` (:obj:`lvlspy.species.Species`) The species of which the - probability vectors are to be calculated for - - ``level_low`` (:obj:`int`, optional) The lower level the effective transition rates are - calculated to. Defaults to 0; the ground state. - - ``level_high`` (:obj:`int`, optional) The higher level the effective transtion rates are - calculated to. Defaults to 1; the first excited state - - Returns: - Upon successful return, the cascade probability vectors will be returned as an array - - ``g1_out`` (:obj:`numpy.array`) cascade vector out of lower level - ``g2_out`` (:obj:`numpy.array`) cascade vector out of higher level - ``g1_in`` (:obj:`numpy.array`) cascade vector into lower level - ``g2_in`` (:obj:`numpy.array`) cascade vector into higher level - """ - - rate_matrix = np.abs(sp.compute_rate_matrix(t)) - trans_props = transfer_properties(rate_matrix, level_low, level_high) - - f_n = _partial_sum(trans_props[0]) - - g1_in = np.matmul(f_n, trans_props[1]) - g2_in = np.matmul(f_n, trans_props[3]) - - g1_out = np.matmul(f_n.T, trans_props[2]) - g2_out = np.matmul(f_n.T, trans_props[4]) - - return [g1_in, g2_in, g1_out, g2_out] - - -def ensemble_weights(t, sp, level_low=0, level_high=1): - """ - Method to calculate the ensemble weights - - Args: - ``t`` (:obj:`float`) The temperature in K - - ``sp`` (:obj:`lvlspy.species.Species`) The species of which the ensemble weights - are to be calculated for - - ``level_low`` (:obj:`int`, optional) The lower level the effective transition rates are - calculated to. Defaults to 0; the ground state. - - ``level_high`` (:obj:`int`, optional) The higher level the effective transtion rates are - calculated to. Defaults to 1; the first excited state - - Returns: - Upon successful return, the ensemble weights and their properties will be returned - as an array - - ``w_low`` (:obj:`numpy.array`) weight factor relative to the low level - ``w_high`` (:obj:`numpy.array`) weight factor relative to the high level - ``W_low`` (:obj:`numpy.float`) Enhancement of ensemble abundance over low level - ``W_high`` (:obj:`numpy.float`) Enhancement of ensemble abundance over high level - ``R_lowk`` (:obj:`numpy.array`) The reverse ratio relative to the low level - ``R_highk`` (:obj:`numpy.array`) The reverse ratio relative to the high level - ``G_low`` (:obj:`numpy.float`) Partition function associated with the low level - ``G_high`` (:obj:`numpy.float`) Patition function associated with the high level - - """ - # calculate the equilibrium probabilities - eq_prob = sp.compute_equilibrium_probabilities(t) - - n = len(eq_prob) - # initialize arrays - w_low, w_high, r_lowk, r_highk = ( - np.empty(n), - np.empty(n), - np.empty(n - 2), - np.empty(n - 1), - ) - - # get the cascade probabilities - # gammas structure = [g1_in, g2_in, g1_out, g2_out] - gammas = cascade_probabilities(t, sp, level_low, level_high) - - for i in range(n - 2): - r_lowk[i] = eq_prob[i + 2] / eq_prob[level_low] - r_highk[i] = eq_prob[i + 2] / eq_prob[level_high] - - for i in range(n): - if i == level_low: - w_low[i] = 1.0 - w_high[i] = 0.0 - elif i == level_high: - w_low[i] = 0.0 - w_high[i] = 1.0 - else: - w_low[i] = gammas[0][i - 2] * r_lowk[i - 2] - w_high[i] = gammas[1][i - 2] * r_highk[i - 2] - - w_low, w_high = np.sum(w_low), np.sum(w_high) - # Calculate the partition functions - levels = sp.get_levels() - g_low = levels[level_low].get_multiplicity() * w_low - g_high = levels[level_high].get_multiplicity() * w_high - - return [w_low, w_high, w_low, w_high, r_lowk, r_highk, g_low, g_high] diff --git a/lvlspy/calculate/weisskopf.py b/lvlspy/calculate/weisskopf.py deleted file mode 100644 index 7324156..0000000 --- a/lvlspy/calculate/weisskopf.py +++ /dev/null @@ -1,256 +0,0 @@ -""" -Module to handle weisskopf calculations -""" - -import numpy as np -import scipy.special as spc -from gslconsts.consts import GSL_CONST_NUM_ZETTA - - -class Weisskopf: - """ - A class for handling Weisskopf related calculations - """ - - def rate_mag(self, e_i, e_f, j, a): - """ - Calculates the transition rate between two levels where - the transition is a magnetic multipole - - Args: - ``e_i`` (:obj:`float`) Energy of the initial state (in keV) - - ``e_f`` (:obj:`float`) Energy of the final state (in keV) - - ``j`` (:obj:`int`) Angular momentum of the gamma ray - - ``a`` (:obj:`int`) Mass number - - - Returns: - The magnetic contribution to the Weisskopf estimate between the two states - """ - - de = e_i - e_f - - s = ( - 1.9 * (j + 1) / (j * np.power(spc.factorial2(2 * j + 1), 2)) - ) * np.power(3 / (j + 3), 2) - - return ( - s - * np.power(de / 197000.0, 2 * j + 1) - * np.power(1.2 * np.power(a, 1.0 / 3.0), 2 * j - 2) - * GSL_CONST_NUM_ZETTA - ) - - def rate_elec(self, e_i, e_f, j, a): - """Calculates the transition rate between two levels where - the transition is an electric multipole. - - Args: - ``e_i`` (:obj:`float`) Energy of the initial state (in keV) - - ``e_f`` (:obj:`float`) Energy of the final state (in keV) - - ``j`` (:obj:`int`) Angular momentum of the gamma ray - - ``a`` (:obj:`int`) Mass number - - - Returns: - The electric contribution to the Weisskopf estimate between the two states - """ - - de = e_i - e_f # energy difference - - s = ( - 4.4 * (j + 1) / (j * np.power(spc.factorial2(2 * j + 1), 2)) - ) * np.power(3 / (j + 3), 2) - - return ( - s - * np.power(de / 197000.0, 2 * j + 1) - * np.power(1.2 * np.power(a, 1.0 / 3.0), 2 * j) - * GSL_CONST_NUM_ZETTA - ) - - def estimate(self, e, j, p, a): - """Calculates the Weisskopf estimate for a transition between two states. - - Args: - ``e`` (:obj:`list`) An array containing the energies of the two levels - - ``j`` (:obj:`list`) An array containing the angular momenta of the two levels - - ``p`` (:obj:`list`) An array containing the parity of both levels - - ``a`` (:obj:`int`) The mass number of the species - - Returns: - ``ein_a`` (:obj:`float`) The Einstein A coefficiet of the downwards transition - """ - ein_a = 0.0 - sm = j[0] + j[1] - df = j[0] - j[1] - j_range = range( - max(1, abs(df)), sm + 1 - ) # range of gamma angular momenta - - for jj in j_range: - ein_a += self._get_rate(jj, p, e, a) - - return ein_a - - def estimate_from_ensdf(self, lvs, tran, a): - """ - Calculates the Weisskopf estimate for a transition between two states based on the - properties available from the ENSDF file. - - Args: - ``lvs`` (:obj:`lvlspy.level.Level`) The levels of the species - - ``tran`` (:obj:`list`) An array containing all the data from - ENSDF regarding a single transition - - Returns: - ``ein_a`` (:obj:`float`) The total estimate for the transition rate - (in per second) using Weisskopf single partice estimate - """ - - e = [lvs[tran[0]].get_energy(), lvs[tran[1]].get_energy()] - j = [ - (lvs[tran[0]].get_multiplicity() - 1) // 2, - (lvs[tran[1]].get_multiplicity() - 1) // 2, - ] - if a % 2 != 0: - j = [ - (lvs[tran[0]].get_multiplicity() - 1) / 2, - (lvs[tran[1]].get_multiplicity() - 1) / 2, - ] - - p = [ - lvs[tran[0]].get_properties()["parity"], - lvs[tran[1]].get_properties()["parity"], - ] - ein_a = 0.0 - j_range = range( - int(max(1, abs(j[0] - j[1]))), int(j[0] + j[1] + 1) - ) # range of gamma angular momenta - - m_r = 0 - - if tran[7] != "": - m_r = float(tran[7]) # mixing ratio - - b = self._get_reduced_trans_prob(tran[16]) - i_tran = [ - np.where(np.strings.find(b, "E") == 1)[0], - np.where(np.strings.find(b, "M") == 1)[0], - ] # first is electric, - # second is magnetic - for c, i in enumerate(i_tran): - - if len(i) == 0: - i = np.array([-1]) - i_tran[c] = i - - if tran[16] == "": - for jj in j_range: - ein_a += self._get_rate(jj, p, e, a) - - else: - for jj in j_range: - - ein_a += self._get_adjusted_rate( - jj, p, e, [a, b, i_tran, m_r, tran] - ) - - return ein_a - - def _get_adjusted_rate(self, jj, p, e, arr): - a = arr[0] - b = arr[1] - i_tran = arr[2] - m_r = arr[3] - tran = arr[4] - b_1 = 1.0 - - if np.power(-1, jj) * p[0] == p[1]: - if b[int(i_tran[0])][1] == "E" and int(b[int(i_tran[0])][2]) == jj: - b_1 = float(b[int(i_tran[0])][5 : len(b[int(i_tran[0])])]) - - if b[int(i_tran[0])][1:3] == "E2" and tran[7] != "": - - b_1 = b_1 * np.power(m_r, 2) / (1.0 + np.power(m_r, 2)) - if b_1 == 1.0: - return self.rate_elec(e[0], e[1], jj, a) / 10 - - return ( - self.rate_elec(e[0], e[1], jj, a) * b_1 / self._b_sp_el(a, jj) - ) - - if b[int(i_tran[1])][1] == "M" and int(b[int(i_tran[1])][2]) == jj: - b_1 = float(b[int(i_tran[1])][5 : len(b[int(i_tran[1])])]) - - if b[int(i_tran[1])][1:3] == "M1" and tran[7] != "": - - b_1 = b_1 / (1.0 + np.power(m_r, 2)) - - if b_1 == 1.0: - return self.rate_mag(e[0], e[1], jj, a) / 10 - - return self.rate_mag(e[0], e[1], jj, a) * b_1 / self._b_sp_ml(jj) - - def _get_rate(self, jj, p, e, a): - - if np.power(-1, jj) * p[0] == p[1]: - return ( - self.rate_elec(e[0], e[1], jj, a) / 10 - ) # Weisskopf estimates in generally over-estimate by a factor of 10 - - return ( - self.rate_mag(e[0], e[1], jj, a) / 10 - ) # Weisskopf estimates in generally over-estimate by a factor of 10 - - def _get_reduced_trans_prob(self, mod_b): - reduced_prob = [] - mods = mod_b.split("$") - if mod_b == "": - reduced_prob.append(mod_b) - return reduced_prob - - if len(mods) == 1: - sp_mods = mods[0].split() - - if len(sp_mods[2]) > 4: - reduced_prob.append(sp_mods[2]) - else: - reduced_prob.append(sp_mods[2] + "=" + sp_mods[4]) - else: - for i, m in enumerate(mods): - sp_mods = m.split() - if i == 0: - if len(sp_mods[2]) > 4: - reduced_prob.append(sp_mods[2]) - else: - reduced_prob.append(sp_mods[2] + "=" + sp_mods[4]) - else: - if len(sp_mods[0]) > 4: - reduced_prob.append(sp_mods[0]) - else: - reduced_prob.append(sp_mods[0] + "=" + sp_mods[2]) - - return reduced_prob - - def _b_sp_ml(self, j): - return ( - 10.0 * np.power(3.0 / (j + 3.0), 2) * np.power(1.2, j - 1) / np.pi - ) - - def _b_sp_el(self, a, j): - return ( - np.power(3.0 / (3.0 + j), 2) - * np.power(1.2 * np.power(a, 1 / 3), 2 * j) - / (4 * np.pi) - )