Skip to content

Commit

Permalink
Rename cumtrapz
Browse files Browse the repository at this point in the history
  • Loading branch information
vatai committed Mar 14, 2024
1 parent 72bbc69 commit f1ec24d
Show file tree
Hide file tree
Showing 2 changed files with 7 additions and 12 deletions.
4 changes: 3 additions & 1 deletion radicalpy/simulation.py
Original file line number Diff line number Diff line change
Expand Up @@ -713,7 +713,9 @@ def product_probability(self, obs: State, rhos: np.ndarray) -> np.ndarray:
@staticmethod
def product_yield(product_probability, time, k):
"""Calculate the product yield and the product yield sum."""
product_yield = k * sp.integrate.cumtrapz(product_probability, time, initial=0)
product_yield = k * sp.integrate.cumulative_trapezoid(
product_probability, time, initial=0
)
product_yield_sum = k * np.trapz(product_probability, dx=time[1])
return product_yield, product_yield_sum

Expand Down
15 changes: 4 additions & 11 deletions tests/radpy.py
Original file line number Diff line number Diff line change
Expand Up @@ -572,14 +572,12 @@ def TimeEvolution(
space="Hilbert",
model="Exponential",
):

# time = np.linspace(t_min, t_max, t_stepsize)
# dt = time[1] - time[0]
time = np.arange(0, t_max, dt)

match space:
case "Hilbert":

HZ = HamiltonianZeeman_RadicalPair(spins, B)
H_total = H + HZ
rho0 = Hilbert_initial(initial, spins, H_total)
Expand All @@ -601,14 +599,15 @@ def TimeEvolution(
K = Kinetics(spins, k, time, model=model)
evol_reaction = evol * K

ProductYield = integrate.cumtrapz(evol_reaction, time, initial=0) * k
ProductYield = (
integrate.cumulative_trapezoid(evol_reaction, time, initial=0) * k
)
ProductYieldSum = np.max(ProductYield)
# print('Product yield: ', '%.2f' % ProductYieldSum)

return [time, evol_reaction, ProductYield, ProductYieldSum, rho]

case "Liouville":

HZ = HamiltonianZeeman_RadicalPair(spins, B)
HZ = Hilbert2Liouville(HZ)
H_total = H + HZ
Expand All @@ -627,7 +626,7 @@ def TimeEvolution(
evol[i] = np.real(np.trace(np.matmul(Pobs.T, rhotL)))
rho.append(rhotL)

ProductYield = integrate.cumtrapz(evol, time, initial=0) * k
ProductYield = integrate.cumulative_trapezoid(evol, time, initial=0) * k
ProductYieldSum = np.max(ProductYield)
# print('Product yield: ', '%.2f' % ProductYieldSum)

Expand Down Expand Up @@ -747,7 +746,6 @@ def BhalfTheoretical(radicalA_effectiveHFC, radicalB_effectiveHFC):
def SingletYieldTwoSite(
B, gamma, theta, phi, r12A, r13A, r23A, r12B, r13B, r23B, kS, kF, tau, spins
):

P12S = projop(spins, "S")
P12SA = np.concatenate((np.ndarray.flatten(P12S), np.zeros(len(P12S) ** 2)))
P12SB = np.concatenate((np.zeros(len(P12S) ** 2), np.ndarray.flatten(P12S)))
Expand Down Expand Up @@ -795,7 +793,6 @@ def SingletYieldTwoSite(
def SingletYieldTwoSiteRelaxation(
B, gamma, theta, phi, r12A, r13A, r23A, r12B, r13B, r23B, kS, kF, tau, spins
):

P12S = projop(spins, "S")
rho0 = P12S
M = 2
Expand Down Expand Up @@ -881,7 +878,6 @@ def SingletYieldTwoSiteApprox(
tstep,
spins,
):

P12S = projop(spins, "S")
HZ = HamiltonianZeeman3D(spins, B, theta, phi, gamma)

Expand Down Expand Up @@ -921,7 +917,6 @@ def SingletYieldTwoSiteTimePropagation(
tstep,
spins,
):

P12S = projop(spins, "S")
P12SA = np.concatenate((np.ndarray.flatten(P12S), np.zeros(len(P12S) ** 2)))
P12SB = np.concatenate((np.zeros(len(P12S) ** 2), np.ndarray.flatten(P12S)))
Expand Down Expand Up @@ -989,7 +984,6 @@ def SingletYieldTwoSiteTimePropagationApprox(
tstep,
spins,
):

P12S = projop(spins, "S")
HZ = HamiltonianZeeman3D(spins, B, theta, phi, gamma)

Expand Down Expand Up @@ -1067,7 +1061,6 @@ def TimeEvolutionPlot2D(
def DensityMatrixPlot2D(
rhot, x_axis_labels, y_axis_labels, space="Hilbert", colourmap="viridis"
):

match space:
case "Hilbert":
fig = plt.figure()
Expand Down

0 comments on commit f1ec24d

Please sign in to comment.