diff --git a/python/_cosmo_power.pyx b/python/_cosmo_power.pyx
index c23da1f..b417bd9 100644
--- a/python/_cosmo_power.pyx
+++ b/python/_cosmo_power.pyx
@@ -17,7 +17,8 @@ cdef extern from "cosmopower.hpp" namespace "CosmoTool":
       POWER_BARDEEN "CosmoTool::CosmoPower::POWER_BARDEEN",
       POWER_SUGIYAMA "CosmoTool::CosmoPower::POWER_SUGIYAMA",
       POWER_BDM,
-      POWER_TEST
+      POWER_TEST,
+      HU_WIGGLES_ORIGINAL "CosmoTool::CosmoPower::HU_WIGGLES_ORIGINAL"
 
   cdef cppclass CosmoPower:
     double n
@@ -113,6 +114,8 @@ cdef class CosmologyPower:
       f = POWER_BARDEEN
     elif funcname=='SUGIYAMA':
       f = POWER_SUGIYAMA
+    elif funcname=='HU_WIGGLES_ORIGINAL':
+      f = HU_WIGGLES_ORIGINAL
 
     self.power.setFunction(f)
     
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 3d59230..a6f241b 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -9,7 +9,6 @@ SET(CosmoTool_SRCS
 	growthFactor.cpp
 	cosmopower.cpp
   cic.cpp
-  tf_fit.c
 )
 
 IF (Boost_FOUND)
diff --git a/src/cosmopower.cpp b/src/cosmopower.cpp
index c6bec8c..9848597 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,9 +70,17 @@ CosmoPower::CosmoPower()
 
   Theta_27 = 2.728/2.7;
 
+  ehu_params = 0;
+  
   updateCosmology();
 }
 
+CosmoPower::~CosmoPower()
+{
+  if (ehu_params)
+    delete ehu_params;
+}
+
 /*
  * This is \hat{tophat}
  */
@@ -108,23 +117,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 +133,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 +363,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..a04442b 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,
@@ -88,6 +91,7 @@ namespace CosmoTool {
     };
 
     CosmoPower();
+    ~CosmoPower();
 
     void setFunction(CosmoFunction f);
 
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);
 }
+
+};
+
+};