diff --git a/src/nuspacesim/config.py b/src/nuspacesim/config.py
index a4e0611..83af5ee 100644
--- a/src/nuspacesim/config.py
+++ b/src/nuspacesim/config.py
@@ -316,7 +316,10 @@ def serialize_rad(self, x: float) -> str:
     """ Maximum Azimuthal Angle (Radians). """
     angle_from_limb: float = np.radians(7)
     """ Angle From Limb. Default (Radians). """
-    cherenkov_light_engine: Literal["Default"] = "Default"  # "CHASM", "EASCherSim"
+    cherenkov_light_engine: Literal[
+        "Default", "Greisen", "Giaser-Hillas"
+    ] = "Default"  # "CHASM", "EASCherSim"
+
     ionosphere: Optional[Ionosphere] = Ionosphere()
     tau_shower: NuPyPropShower = NuPyPropShower()
     """ Tau Shower Generator. """
@@ -375,16 +378,22 @@ def config_from_fits(filename: str) -> NssConfig:
     def v(key: str):
         fullkey = "Config " + key
         if fullkey not in h:
-            raise KeyError(fullkey)
+            raise KeyError(f"Missing required key '{fullkey}' in FITS header.")
         return h[fullkey]
 
     # header (d)etector config value assocciated with partial key string.
     def d(key: str):
-        return v("detector " + key)
+        try:
+            return v("detector " + key)
+        except KeyError as e:
+            raise KeyError(f"Detector configuration key error: {e}")
 
     # header (s)etector config value assocciated with partial key string.
     def s(key: str):
-        return v("simulation " + key)
+        try:
+            return v("simulation " + key)
+        except KeyError as e:
+            raise KeyError(f"Simulation configuration key error: {e}")
 
     c = {
         "detector": {
diff --git a/src/nuspacesim/simulation/eas_optical/cphotang.py b/src/nuspacesim/simulation/eas_optical/cphotang.py
index 1c31b49..2a365df 100644
--- a/src/nuspacesim/simulation/eas_optical/cphotang.py
+++ b/src/nuspacesim/simulation/eas_optical/cphotang.py
@@ -43,10 +43,11 @@
 
 import dask.bag as db
 import numpy as np
-from dask.diagnostics import ProgressBar
+from dask.diagnostics.progress import ProgressBar
 from numpy.polynomial import Polynomial
 
 from .detector_geometry import distance_to_detector
+from .shower_properties import gaisser_hillas_particle_count, greisen_particle_count
 
 # Wrapped in try-catch block as a hack to enable sphinx documentation to be generated
 # on ReadTheDocs without pre-compiling.
@@ -55,13 +56,18 @@
 except ImportError:
     pass
 
+try:
+    from importlib.resources import as_file, files
+except ImportError:
+    from importlib_resources import as_file, files
+
 __all__ = ["CphotAng"]
 
 
 class CphotAng:
     r"""Cherenkov Photon Angle"""
 
-    def __init__(self, detector_altitude):
+    def __init__(self, detector_altitude, longitudinal_profile_func="Greisen"):
         r"""CphotAng: Cherenkov photon density and angle determination class.
 
         Iterative summation of cherenkov radiation reimplemented in numpy and
@@ -244,8 +250,20 @@ def __init__(self, detector_altitude):
         self.zMaxZ = self.dtype(65.0)
         self.RadE = self.dtype(6378.14)
 
-        Zair = self.dtype(7.4)
-        self.ecrit = self.dtype(0.710 / (Zair + 0.96))
+        # Longitudinal Profile Funciton selection
+        if longitudinal_profile_func == "Greisen":
+            self.particle_count = greisen_particle_count
+        elif longitudinal_profile_func == "Gaisser-Hillas":
+            with as_file(
+                files("nuspacesim.data.CONEX_table")
+                / "dumpGH_conex_pi_E17_95deg_0km_eposlhc_1394249052_211.dat"
+            ) as file:
+                CONEX_table = np.loadtxt(file, usecols=(4, 5, 6, 7, 8, 9))
+                self.particle_count = (
+                    lambda *args, **kwargs: gaisser_hillas_particle_count(
+                        CONEX_table, *args, **kwargs
+                    )
+                )
 
     def theta_view(self, ThetProp):
         """
@@ -512,19 +530,19 @@ def valid_arrays(self, zsave, delgram, gramsum, gramz, ZonZ, ThetPrpA, Eshow):
         t = np.zeros_like(zsave, dtype=self.dtype)
         t[mask] = gramsum[mask] / self.dtype(36.66)
 
-        greisen_beta = self.dtype(np.log(np.float64(Eshow) / np.float64(self.ecrit)))
+        Zair = self.dtype(7.4)
+        ecrit = self.dtype(0.710 / (Zair + 0.96))
+        greisen_beta = self.dtype(np.log(np.float64(Eshow) / np.float64(ecrit)))
         s = np.zeros_like(zsave, dtype=self.dtype)
         s[mask] = self.dtype(3) * t[mask] / (t[mask] + self.dtype(2) * greisen_beta)
 
-        RN = np.zeros_like(zsave, dtype=self.dtype)
-        RN[mask] = (
-            self.dtype(0.31)
-            / np.sqrt(greisen_beta, dtype=self.dtype)
-            * np.exp(
-                t[mask] * (1 - self.dtype(3 / 2) * np.log(s[mask], dtype=self.dtype)),
-                dtype=self.dtype,
-            )
+        # Greisen  or  Gaisser-Hillas  Longitudinal Paramaterization
+        rn, mask = self.particle_count(
+            mask=mask, t=t, s=s, greisen_beta=greisen_beta, gramsum=gramsum, Eshow=Eshow
         )
+
+        RN = np.zeros_like(zsave, dtype=self.dtype)
+        RN[mask] = rn
         RN[RN < 0] = self.dtype(0)
 
         mask &= ~((RN < 1) & (s > 1))
diff --git a/src/nuspacesim/simulation/eas_optical/eas.py b/src/nuspacesim/simulation/eas_optical/eas.py
index 2033fc9..63b113a 100644
--- a/src/nuspacesim/simulation/eas_optical/eas.py
+++ b/src/nuspacesim/simulation/eas_optical/eas.py
@@ -53,7 +53,10 @@ class EAS:
 
     def __init__(self, config: NssConfig):
         self.config = config
-        self.CphotAng = CphotAng(self.config.detector.initial_position.altitude)
+        self.CphotAng = CphotAng(
+            self.config.detector.initial_position.altitude,
+            self.config.simulation.cherenkov_light_engine,
+        )
 
     @decorators.nss_result_store("altDec", "lenDec")
     def altDec(self, beta, tauBeta, tauLorentz, u=None, *args, **kwargs):
@@ -88,7 +91,7 @@ def __call__(
         init_long,
         *args,
         cloudf=None,
-        **kwargs
+        **kwargs,
     ):
         """
         Electromagnetic Air Shower operation.
diff --git a/src/nuspacesim/simulation/eas_optical/shower_properties.py b/src/nuspacesim/simulation/eas_optical/shower_properties.py
index aadbaae..89316f9 100644
--- a/src/nuspacesim/simulation/eas_optical/shower_properties.py
+++ b/src/nuspacesim/simulation/eas_optical/shower_properties.py
@@ -97,7 +97,7 @@ def shower_age(T):
     return 3.0 * T / (T + 41.77325895999150334743982471)
 
 
-def greisen_particle_count(T, s):
+def greisen_particle_count(t, s, greisen_beta, mask, *args, dtype=np.float32, **kwargs):
     r"""Particle count as a function of radiation length from atmospheric depth
 
     Hillas 1461 eqn (6)
@@ -106,12 +106,88 @@ def greisen_particle_count(T, s):
     (0.31 / sqrt (10^8 / (0.710 / 8.36)))
     = 0.0678308895484773316048795658058110209448440898800928880798622962...
     """
+
     # , param_beta=np.log(10 ** 8 / (0.710 / 8.36))
     # N_e = (0.31 / np.sqrt(param_beta)) * np.exp(T * (1.0 - 1.5 * np.log(s)))
     # N_e[N_e < 0] = 0.0
-    N_e = 0.067830889548477331 * np.exp(T * (1.0 - 1.5 * np.log(s)))
-    N_e[N_e < 0] = 0.0
-    return N_e
+    alpha = dtype(0.31) / np.sqrt(greisen_beta, dtype=dtype)
+
+    RN = alpha * np.exp(
+        t * (dtype(1) - dtype(1.5) * np.log(s, dtype=dtype)), dtype=dtype
+    )
+    RN[RN < 0] = 0.0
+    return RN, mask
+
+
+def gaisser_hillas_particle_count(
+    CONEX_table, gramsum, Eshow, mask, *args, dtype=np.float32, **kwargs
+):
+    # Gaisser Hillas Values from CONEX File
+    idx = np.random.randint(low=0, high=CONEX_table.shape[0])
+    Nm, Xm, X0, p1, p2, p3 = CONEX_table[idx]
+    X0 = 0.0 if X0 < 0.0 else X0
+
+    # Masking Gramsum
+    gramsum_mask = gramsum > X0
+    mask &= gramsum_mask
+    Xmask = gramsum[gramsum_mask]
+
+    Nmax = Nm * (Eshow * 1.0e-8)
+    Xmax = Xm + (70.0 * np.log10(Eshow * 1.0e-8))
+    λ = p1 + p2 * Xmask + p3 * Xmask * Xmask
+    λ[λ > 100.0] = 100.0
+    λ[λ < 0.0] = 1.0
+
+    # Parametric Form Parameters
+    x = (Xmask - X0) / λ
+    m = (Xmax - X0) / λ
+    return Nmax * np.exp(m * (np.log(x) - np.log(m)) - (x - m)), mask
+
+
+# def
+# # Implamenting CONEX Data
+# with as_file(
+#     files("nuspacesim.data.CONEX_table")
+#     / "dumpGH_conex_pi_E17_95deg_0km_eposlhc_1394249052_211.dat"
+# ) as file:
+#     self.CONEX_table = np.loadtxt(file, usecols=(4, 5, 6, 7, 8, 9))
+
+# Gaisser Hillas Values from CONEX File
+# idx=np.random.randint(low=0,high=self.CONEX_table.shape[0])
+# Nm, Xm, X0, p1, p2, p3 = self.CONEX_table[idx]
+# X0=0.0 if X0<0.0 else X0
+# Nmax=Nm*(Eshow/1.e8)
+# Xmax = Xm + (70. * np.log10(Eshow / 1.e8))
+# λ=p1 + p2 * Xmask + p3 * Xmask * Xmask # val7+val8*t+val9*t*t
+# λ[λ > 100] = 100
+# λ[λ < 0] = 1
+
+# Parameters for Gaisser Hillas (Nmax and Xmax Scaled)
+# Nmax= 0.045 * (1.+0.0217*(np.log(Eshow/1.e5)))*(Eshow/0.074)
+# Xmax= 36.* np.log(Eshow/0.074)
+# X0= 0.
+# invlam = 1/70.
+
+# p Parameter for scaled equation
+# p=(Xmax/λ)-1
+
+# Masking Gramsum
+# gramsum_mask = gramsum > X0
+# mask &= gramsum_mask
+# Xmask=(gramsum[gramsum_mask])
+
+# Parametric Form Parameters
+# x=(Xmask-X0)/λ
+# m=(Xmax-X0)/λ
+
+# RN = np.zeros_like(zsave, dtype=self.dtype)
+# RN[mask] = (
+
+# Gaisser Hillas Parametric Form Equation (to get rid of overflow errors)
+# Nmax * np.exp(m*(np.log(x)-np.log(m))-(x-m))
+
+# Gaisser Hillas Equation From Shower Properties
+# gaisser_hillas_particle_count(gramsum, Nmax, X0, Xmax, invlam)
 
 
 def shower_age_of_greisen_particle_count(target_count, x0=2):
@@ -130,13 +206,6 @@ def rns(s):
     return newton(rns, x0)
 
 
-def gaisser_hillas_particle_count(X, Nmax, X0, Xmax, invlam):
-    # return ((X - X0) / (Xmax - X0)) ** xmax * np.exp((Xmax - X) * invlam)
-    xmax = (Xmax - X0) * invlam
-    x = (X - X0) * invlam
-    return Nmax * (x / xmax) ** xmax * np.exp(xmax - x)
-
-
 def slant_depth_trig_approx(z_lo, z_hi, theta_tr, z_max=100.0):
     rho = us_std_atm_density
     r0 = rho(0)