From 9e424d562d7b6faa20284cea07cfc21c4416f65a Mon Sep 17 00:00:00 2001
From: Alex Reustle <areustledev@gmail.com>
Date: Thu, 2 Dec 2021 10:10:14 -0500
Subject: [PATCH] Power Spectrum LogEnergy Distribution

---
 src/nuspacesim/config.py                         |  2 +-
 src/nuspacesim/simulation/spectra/local_plots.py |  7 ++++++-
 src/nuspacesim/simulation/spectra/spectra.py     | 15 +++++++--------
 3 files changed, 14 insertions(+), 10 deletions(-)

diff --git a/src/nuspacesim/config.py b/src/nuspacesim/config.py
index bd4174c..d018666 100644
--- a/src/nuspacesim/config.py
+++ b/src/nuspacesim/config.py
@@ -135,7 +135,7 @@ def __call__(self) -> dict:
 
 @dataclass
 class PowerSpectrum:
-    index: float = 0.5
+    index: float = 2.0
     """Power Law Log Energy of the tau neutrinos in GeV."""
 
     lower_bound: float = 6.0
diff --git a/src/nuspacesim/simulation/spectra/local_plots.py b/src/nuspacesim/simulation/spectra/local_plots.py
index 573d9a3..3c95447 100644
--- a/src/nuspacesim/simulation/spectra/local_plots.py
+++ b/src/nuspacesim/simulation/spectra/local_plots.py
@@ -41,9 +41,14 @@ def spectra_histogram(inputs, results, *args, **kwargs):
     log_e_nu = results
 
     color = "g"
-    fig, ax = plt.subplots(1, constrained_layout=True)
+    fig = plt.figure(figsize=(8, 7), constrained_layout=True)
+    ax = fig.add_subplot(211)
     ax.hist(log_e_nu, 100, log=False, facecolor=color)
     ax.set_xlabel(f"log(E_nu) of {N} events")
 
+    ax = fig.add_subplot(212)
+    ax.hist(log_e_nu, 100, log=True, facecolor=color)
+    ax.set_xlabel(f"log(E_nu) of {N} events")
+
     fig.suptitle("Energy Spectra Histogram, Log(E_nu)")
     plt.show()
diff --git a/src/nuspacesim/simulation/spectra/spectra.py b/src/nuspacesim/simulation/spectra/spectra.py
index 4c46b4a..5923428 100644
--- a/src/nuspacesim/simulation/spectra/spectra.py
+++ b/src/nuspacesim/simulation/spectra/spectra.py
@@ -55,10 +55,12 @@ def energy_spectra(
         return np.full(shape=(N), fill_value=spectra.log_nu_tau_energy)
 
     if isinstance(spectra, PowerSpectrum):
-        delta = 10 ** spectra.upper_bound - 10 ** spectra.lower_bound
-        return 10 ** spectra.lower_bound + delta * np.random.power(
-            spectra.index, size=N
-        )
+        p = spectra.index
+        a = 10 ** spectra.lower_bound
+        b = 10 ** spectra.upper_bound
+        mp = 1 - p
+        u = np.random.uniform(0.0, 1.0 + np.finfo(np.float64).eps, size=N)
+        return np.reciprocal(mp) * np.log10(u * (b ** mp - a ** mp) + a ** mp)
 
     if isinstance(spectra, Callable):
         return spectra(*args, size=N, **kwargs)
@@ -74,7 +76,4 @@ def __init__(self, config: NssConfig):
         self.config = config
 
     def __call__(self, N, *args, **kwargs):
-        es = energy_spectra(N, self.config.simulation.spectrum, *args, **kwargs)
-        print(np.min(es), np.max(es))
-        print(np.log10(np.min(es)), np.log10(np.max(es)))
-        return np.log10(es)
+        return energy_spectra(N, self.config.simulation.spectrum, *args, **kwargs)