From db1d63b04f1f48ffbd9148416825e22dd12b1bc1 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <guilhem.lavaux@gmail.com>
Date: Tue, 28 Apr 2015 22:08:39 +0200
Subject: [PATCH] Removed global variables

---
 src/CMakeLists.txt           |  1 -
 src/cosmopower.cpp           | 32 ++++++++++++--------------
 src/cosmopower.hpp           |  3 +++
 src/{tf_fit.c => tf_fit.hpp} | 44 ++++++++----------------------------
 4 files changed, 27 insertions(+), 53 deletions(-)
 rename src/{tf_fit.c => tf_fit.hpp} (90%)

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 21f2691..659d4cc 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -10,7 +10,6 @@ SET(CosmoTool_SRCS
 	growthFactor.cpp
 	cosmopower.cpp
   cic.cpp
-  tf_fit.c
 )
 
 IF (ENABLE_OPENMP)
diff --git a/src/cosmopower.cpp b/src/cosmopower.cpp
index eb36b2c..87baf6d 100644
--- a/src/cosmopower.cpp
+++ b/src/cosmopower.cpp
@@ -41,6 +41,7 @@ knowledge of the CeCILL license and that you accept its terms.
 #include <fstream>
 #include <gsl/gsl_integration.h>
 #include "cosmopower.hpp"
+#include "tf_fit.hpp"
 
 using namespace std;
 using namespace CosmoTool;
@@ -69,6 +70,8 @@ CosmoPower::CosmoPower()
 
   Theta_27 = 2.728/2.7;
 
+  ehu_params = 0;
+  
   updateCosmology();
 }
 
@@ -108,23 +111,6 @@ static double powG(double y)
 
   return y * (-6 * a + (2 + 3 * y) *log((a + 1)/(a - 1)));
 }
-
-
-extern "C" {
-  void TFset_parameters(float omega0hh, float f_baryon, float Tcmb);
-  float TFfit_onek(float k, float *tf_baryon, float *tf_cdm); 
-
-  void TFfit_hmpc(float omega0, float f_baryon, float hubble, float Tcmb,
-	  int numk, float *k, float *tf_full, float *tf_baryon, float *tf_cdm);
-
-  float TFsound_horizon_fit(float omega0, float f_baryon, float hubble);
-  float TFk_peak(float omega0, float f_baryon, float hubble);
-  float TFnowiggles(float omega0, float f_baryon, float hubble, 
-		  float Tcmb, float k_hmpc);
-  float TFzerobaryon(float omega0, float hubble, float Tcmb, float k_hmpc);
-}
-
-
  
 double CosmoPower::powerEfstathiou(double k)
 {
@@ -141,7 +127,10 @@ double CosmoPower::powerEfstathiou(double k)
 
 void CosmoPower::updateHuWigglesOriginal()
 {
-  TFset_parameters( (OMEGA_C+OMEGA_B)*h*h,
+  if (ehu_params == 0)
+    ehu_params = new TF_Transfer();
+  
+  ehu_params->TFset_parameters( (OMEGA_C+OMEGA_B)*h*h,
                      OMEGA_B/(OMEGA_C+OMEGA_B), Theta_27*2.7);
 }
 
@@ -368,6 +357,13 @@ double CosmoPower::power(double k)
 
 double CosmoPower::powerHuWigglesOriginal(double k)
 {
+  float tfb, tfc;
+  
+  ehu_params->TFfit_onek(k, &tfb, &tfc);
+
+  double T_k = OMEGA_B/OMEGA_0 * tfb + OMEGA_C/OMEGA_0 * tfc;  
+
+  return normPower * pow(k,n) * T_k * T_k;
 }
 
 void CosmoPower::setFunction(CosmoFunction f)
diff --git a/src/cosmopower.hpp b/src/cosmopower.hpp
index 21ef45a..3a395b4 100644
--- a/src/cosmopower.hpp
+++ b/src/cosmopower.hpp
@@ -38,6 +38,8 @@ knowledge of the CeCILL license and that you accept its terms.
 
 namespace CosmoTool {
 
+  struct TF_Transfer;
+    
   class CosmoPower
   {
   public:
@@ -74,6 +76,7 @@ namespace CosmoTool {
     };
     
     EHuParams ehu;
+    TF_Transfer *ehu_params;
 
     enum CosmoFunction {
       POWER_EFSTATHIOU,
diff --git a/src/tf_fit.c b/src/tf_fit.hpp
similarity index 90%
rename from src/tf_fit.c
rename to src/tf_fit.hpp
index 5a59649..affcb1d 100644
--- a/src/tf_fit.c
+++ b/src/tf_fit.hpp
@@ -14,17 +14,6 @@ calculate other quantities given in Section 4 of the paper. */
 
 #include <math.h>
 #include <stdio.h>
-void TFset_parameters(float omega0hh, float f_baryon, float Tcmb);
-float TFfit_onek(float k, float *tf_baryon, float *tf_cdm); 
-
-void TFfit_hmpc(float omega0, float f_baryon, float hubble, float Tcmb,
-	int numk, float *k, float *tf_full, float *tf_baryon, float *tf_cdm);
-
-float TFsound_horizon_fit(float omega0, float f_baryon, float hubble);
-float TFk_peak(float omega0, float f_baryon, float hubble);
-float TFnowiggles(float omega0, float f_baryon, float hubble, 
-		float Tcmb, float k_hmpc);
-float TFzerobaryon(float omega0, float hubble, float Tcmb, float k_hmpc);
 
 /* ------------------------ DRIVER ROUTINE --------------------------- */
 /* The following is an example of a driver routine you might use. */
@@ -65,26 +54,6 @@ you would be better off not simply #include'ing this file in your programs,
 but rather compiling it separately, calling only the driver, and using
 extern declarations to access the intermediate quantities. */
 
-/* ----------------------------- DRIVER ------------------------------- */
-
-void TFfit_hmpc(float omega0, float f_baryon, float hubble, float Tcmb,
-	int numk, float *k, float *tf_full, float *tf_baryon, float *tf_cdm)
-/* Remember: k[0..numk-1] is in units of h Mpc^-1. */
-{
-    int j;
-    float tf_thisk, baryon_piece, cdm_piece;
-
-    TFset_parameters(omega0*hubble*hubble, f_baryon, Tcmb);
-
-    for (j=0;j<numk;j++) {
-    	tf_thisk = TFfit_onek(k[j]*hubble, &baryon_piece, &cdm_piece); 
-	if (tf_full!=NULL) tf_full[j] = tf_thisk;
-	if (tf_baryon!=NULL) tf_baryon[j] = baryon_piece;
-	if (tf_cdm!=NULL) tf_cdm[j] = cdm_piece;
-    }
-    return;
-}
-
 /* ------------------------ FITTING FORMULAE ROUTINES ----------------- */
 
 /* There are two routines here.  TFset_parameters() sets all the scalar
@@ -97,6 +66,9 @@ global variables in case you wish to access them, e.g. by declaring
 them as extern variables in your main program. */
 /* Note that all internal scales are in Mpc, without any Hubble constants! */
 
+namespace CosmoTool {
+struct TF_Transfer {
+
 float	omhh,		/* Omega_matter*h^2 */
 	obhh,		/* Omega_baryon*h^2 */
 	theta_cmb,	/* Tcmb in units of 2.7 K */
@@ -117,11 +89,11 @@ float	omhh,		/* Omega_matter*h^2 */
 	alpha_gamma;	/* Gamma suppression in approximate TF */
 
 /* Convenience from Numerical Recipes in C, 2nd edition */
-static float sqrarg;
+ float sqrarg;
 #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
-static float cubearg;
+ float cubearg;
 #define CUBE(a) ((cubearg=(a)) == 0.0 ? 0.0 : cubearg*cubearg*cubearg)
-static float pow4arg;
+ float pow4arg;
 #define POW4(a) ((pow4arg=(a)) == 0.0 ? 0.0 : pow4arg*pow4arg*pow4arg*pow4arg)
 	/* Yes, I know the last one isn't optimal; it doesn't appear much */
 
@@ -335,3 +307,7 @@ and omega0 -> omega0*hubble^2. */
     T_0_C0 = 14.2 + 731.0/(1+62.5*q);
     return T_0_L0/(T_0_L0+T_0_C0*q*q);
 }
+
+};
+
+};