diff --git a/zeus21/constants.py b/zeus21/constants.py
index 6016a68..ae60756 100644
--- a/zeus21/constants.py
+++ b/zeus21/constants.py
@@ -12,6 +12,7 @@
 ######################
 
 
+
 precisionboost = 1.0 #makes integrals take more points for boost in precision. Baseline = 1.0
 FLAG_FORCE_LINEAR_CF = 0 #0 to do standard calculation, 1 to force linearization of correlation function
 
@@ -54,7 +55,7 @@
 EN_ION_HeI = 24.59 #eV
 sigma0norm = 1e-18 #cm^2 normalization of xray cross sections
 
-ZMAX_INTEGRAL = 35.0 #at which z we start the integrals. We take 35. as fiducial since there is not much SFRD. Make sure to test if you have exotic cosmology/astrophysics
+ZMAX_INTEGRAL = 50.0 #35.0 #at which z we start the integrals. We take 35. as fiducial since there is not much SFRD. Make sure to test if you have exotic cosmology/astrophysics
 
 #LyA related constants
 wavelengthLyC = 91.1753 ##lyman continuum in nm
@@ -75,6 +76,15 @@
 Tstar_21 = 0.0682 #T* in K for the 21cm transition
 A10_21 = 2.85e-15 #1/s, Einstein 10 coeff for HI
 
+############################################################################################
+
+Tk_HH = [1, 2, 4, 6, 8, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100, 200, 300, 500, 700, 1000, 2000, 3000, 5000, 7000, 10000]
+k_HH = [1.38e-13, 1.43e-13, 2.71e-13, 6.60e-13, 1.47e-12, 2.88e-12, 9.10e-12, 1.78e-11, 2.73e-11, 3.67e-11, 5.38e-11, 6.86e-11, 8.14e-11, 9.25e-11, 1.02e-10, 1.11e-10, 1.19e-10, 1.75e-10, 2.09e-10, 2.56e-10, 2.91e-10, 3.31e-10, 4.27e-10, 4.97e-10, 6.03e-10, 6.87e-10, 7.87e-10]
+Tk_He = [1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 3000, 5000, 7000, 10000, 15000, 20000]
+k_He = [2.39e-10, 3.37e-10, 5.30e-10, 7.46e-10, 1.05e-9, 1.63e-9, 2.26e-9, 3.11e-9, 4.59e-9, 5.92e-9, 7.15e-9, 7.71e-9, 8.17e-9, 8.32e-9, 8.37e-9, 8.29e-9, 8.11e-9]
+
+############################################################################################
+
 
 #whether to renormalize the C2 oefficients (appendix in 2302.08506)
 C2_RENORMALIZATION_FLAG = 1 - FLAG_FORCE_LINEAR_CF
diff --git a/zeus21/sfrd.py b/zeus21/sfrd.py
index 938f87a..5616fc3 100644
--- a/zeus21/sfrd.py
+++ b/zeus21/sfrd.py
@@ -39,8 +39,11 @@ def __init__(self, Cosmo_Parameters, ClassCosmo, Astro_Parameters, HMF_interpola
         self.SFRD_avg = np.zeros_like(self.zintegral)
         self.SFRDbar2D = np.zeros((self.Nzintegral, Cosmo_Parameters.NRs)) #SFR at z=zprime when averaged over a radius R (so up to a higher z)
         self.gamma_index2D = np.zeros_like(self.SFRDbar2D) #index of SFR ~ exp(\gamma delta)
+        self.gamma2_index2D = np.zeros_like(self.SFRDbar2D) #quadratic fit gamma like ~ exp(\gamma \delta + \gamma_2 \delta^2)
+        
         self.niondot_avg = np.zeros_like(self.zintegral) #\dot nion at each z (int d(SFRD)/dM *fesc(M) dM)/rhobaryon
-        self.gamma_Niondot_index2D = np.zeros_like(self.SFRDbar2D) #index of SFR ~ exp(\gamma delta)
+        self.gamma_Niondot_index2D = np.zeros_like(self.SFRDbar2D) #index of SFR ~ exp(\gamma \delta)
+        self.gamma2_Niondot_index2D = np.zeros_like(self.SFRDbar2D) #index of SFR ~ exp(\gamma \delta + \gamma_2 \delta^2)
 
         self.ztabRsmoo = np.zeros_like(self.SFRDbar2D) #z's that correspond to each Radius R around each zp
 
@@ -64,8 +67,8 @@ def __init__(self, Cosmo_Parameters, ClassCosmo, Astro_Parameters, HMF_interpola
 
 
         #and EPS factors
-        Nsigmad = 0.2 #how many sigmas we explore
-        Nds = 2 #how many deltas
+        Nsigmad = .2 #how many sigmas we explore (.2 originally)
+        Nds = 3 #how many deltas (2 originally)
         deltatab_norm = np.linspace(-Nsigmad,Nsigmad,Nds)
 
 
@@ -197,6 +200,7 @@ def __init__(self, Cosmo_Parameters, ClassCosmo, Astro_Parameters, HMF_interpola
                     deltatab = deltatab_norm * sigmaR
 
                 SFRD_delta = np.zeros_like(deltatab)
+                Niondot_delta = np.zeros_like(deltatab)
 
                 nu0 = Cosmo_Parameters.delta_crit_ST/sigmacurr #the EPS delta = 0 result. Note 1/sigmacurr and not sigmamod
                 nu0[indextoobig]=1.0 #set to 1 to avoid under/overflows, we don't sum over those masses since they're too big
@@ -220,10 +224,38 @@ def __init__(self, Cosmo_Parameters, ClassCosmo, Astro_Parameters, HMF_interpola
 
                     integrand_EPS *= (1.0 + dd) #to convert from Lagrangian to Eulerian
                     SFRD_delta[idelta] = np.trapz(integrand_EPS, HMF_interpolator.logtabMh)
-
-
-
-                self.gamma_index2D[izp,iR] = np.log(SFRD_delta[-1]/SFRD_delta[0])/(deltatab[-1] - deltatab[0]) #notice its der wrt delta, not delta/sigma or growth
+                    Niondot_delta[idelta] = np.trapz(integrand_EPS * fesctab, HMF_interpolator.logtabMh) #ionizing emissivity (Niondot) is basically SFRD but with an escape fraction that varies with halo mass
+
+
+                midpoint = len(SFRD_delta)//2
+
+                self.gamma_index2D[izp,iR] = np.log(SFRD_delta[midpoint+1]/SFRD_delta[midpoint-1])/(deltatab[midpoint+1] - deltatab[midpoint-1]) #notice its der wrt delta, not delta/sigma or growth
+                
+                self.gamma_Niondot_index2D[izp, iR] = np.log(Niondot_delta[midpoint+1]/Niondot_delta[midpoint-1])/(deltatab[midpoint+1] - deltatab[midpoint-1])
+                
+                
+                '''if (izp == 35)*(iR == 16): #z~10 and R~10Mpc 
+                    print(izp)
+                    print(iR)
+                    print(zp)
+                    print(RR)
+                    self.testSFRD = SFRD_delta
+                    self.testdeltatab = deltatab
+                    self.testNiondot = Niondot_delta
+                    self.testfesctab = fesctab'''
+                    
+                
+                
+                der1 =  np.log(SFRD_delta[midpoint]/SFRD_delta[midpoint-1])/(deltatab[midpoint] - deltatab[midpoint-1]) #ln(y2/y1)/(x2-x1)
+                der2 =  np.log(SFRD_delta[midpoint+1]/SFRD_delta[midpoint])/(deltatab[midpoint+1] - deltatab[midpoint]) #ln(y3/y2)/(x3-x2)
+                
+                self.gamma2_index2D[izp,iR] = (der2 - der1)/(deltatab[midpoint+1] - deltatab[midpoint-1]) #second derivative: (der2-der1)/((x3-x1)/2)
+                
+                der1_N = np.log(Niondot_delta[midpoint]/Niondot_delta[midpoint-1])/(deltatab[midpoint] - deltatab[midpoint-1])
+                der2_N =  np.log(Niondot_delta[midpoint+1]/Niondot_delta[midpoint])/(deltatab[midpoint+1] - deltatab[midpoint])
+                
+                self.gamma2_Niondot_index2D[izp, iR] = (der2_N - der1_N)/(deltatab[midpoint+1] - deltatab[midpoint-1])
+                
 
 
             self.coeff2LyAzpRR[izp] = self.Rtabsmoo * self.dlogRR * self.SFRDbar2D[izp,:] * LyAintegral/ constants.yrTos/constants.Mpctocm**2
@@ -305,6 +337,11 @@ def __init__(self, Cosmo_Parameters, ClassCosmo, Astro_Parameters, HMF_interpola
             self._coeff_Ja_xa_0 = 1.66e11/(1+self.zintegral) #They use a fixed (and slightly ~10% off) value.
         else:
             self._coeff_Ja_xa_0 = 8.0*np.pi*(constants.wavelengthLyA/1e7)**2 * constants.widthLyA * constants.Tstar_21/(9.0*constants.A10_21*self.T_CMB) #units of (cm^2 s Hz sr), convert from Ja to xa. should give 1.81e11/(1+self.zintegral) for Tcmb_0=2.725 K
+            
+        #collision coefficient fh(1-x_e)
+        self.xc_HH = Cosmo_Parameters.f_H * (1.0 - self.xe_avg) * cosmology.n_baryon(Cosmo_Parameters, self.zintegral) * np.interp(self.Tk_avg, constants.Tk_HH, constants.k_HH) / constants.A10_21 * constants.Tstar_21 / cosmology.Tcmb(ClassCosmo, self.zintegral)
+        self.xc_He = self.xe_avg * cosmology.n_baryon(Cosmo_Parameters, self.zintegral) * np.interp(self.Tk_avg, constants.Tk_He, constants.k_He) / constants.A10_21 * constants.Tstar_21 / cosmology.Tcmb(ClassCosmo, self.zintegral) #xe
+        self.xc_avg = self.xc_HH + self.xc_He
 
 
         if(constants.FLAG_WF_ITERATIVE==True): #iteratively find Tcolor and Ts. Could initialize one to zero, but this should converge faster
@@ -314,7 +351,7 @@ def __init__(self, Cosmo_Parameters, ClassCosmo, Astro_Parameters, HMF_interpola
             self.invTcol_avg = 1.0 / self.Tk_avg
             self.coeff_Ja_xa = self._coeff_Ja_xa_0 * Salpha_exp(self.zintegral, self.Tk_avg, self.xe_avg)
             self.xa_avg = self.coeff_Ja_xa * self.Jalpha_avg
-            self._invTs_avg = (1.0/self.T_CMB+self.xa_avg*self.invTcol_avg)/(1+self.xa_avg)
+            self._invTs_avg = (1.0/self.T_CMB+self.xa_avg*self.invTcol_avg+self.xc_avg*self.invTk_avg)/(1+self.xa_avg+self.xc_avg)
 
             _invTs_tryfirst = self._invTs_avg #so the loop below does not trigger
 
@@ -332,7 +369,7 @@ def __init__(self, Cosmo_Parameters, ClassCosmo, Astro_Parameters, HMF_interpola
             self.invTcol_avg = 1.0/self.Tk_avg + constants.gcolorfactorHirata * 1.0/self.Tk_avg * (_invTs_tryfirst - 1.0/self.Tk_avg)
 
             #and finally Ts^-1
-            self._invTs_avg = (1.0/self.T_CMB+self.xa_avg * self.invTcol_avg)/(1+self.xa_avg)
+            self._invTs_avg = (1.0/self.T_CMB+self.xa_avg * self.invTcol_avg + self.xc_avg * 1.0/self.Tk_avg)/(1+self.xa_avg+self.xc_avg)
 
 
 
@@ -376,7 +413,7 @@ def __init__(self, Cosmo_Parameters, ClassCosmo, Astro_Parameters, HMF_interpola
 
 
         #and finally, get the signal
-        self.T21avg = cosmology.T021(Cosmo_Parameters,self.zintegral) * self.xa_avg/(1.0 + self.xa_avg) * (1.0 - self.T_CMB * self.invTcol_avg) * self.xHI_avg
+        self.T21avg = cosmology.T021(Cosmo_Parameters,self.zintegral) * (self.xa_avg + self.xc_avg)/(1.0 + self.xa_avg + self.xc_avg) * (1.0 - self.T_CMB * self.invTcol_avg) * self.xHI_avg