diff --git a/python/_cosmo_power.pyx b/python/_cosmo_power.pyx
index bf4bd4b..78305f8 100644
--- a/python/_cosmo_power.pyx
+++ b/python/_cosmo_power.pyx
@@ -39,10 +39,14 @@ cdef extern from "cosmopower.hpp" namespace "CosmoTool":
     double OmegaEff
     double Gamma0
     double normPower
+    double A_BAO
+    double r_BAO
+    double k_D_BAO
     
         
     CosmoPower()
     void setFunction(CosmoFunction)
+    void setFunction_BAO(CosmoFunction,double,double,double)
     void updateCosmology()
     void updatePhysicalCosmology()
     void normalize(double,double)
@@ -101,8 +105,7 @@ cdef class CosmologyPower:
     """
     self.power.SIGMA8 = s8
     self.power.normalize(k_min, k_max)
-    
-    
+   
   def setFunction(self,funcname):
     """setFunction(self, funcname)
     
diff --git a/src/cosmopower.cpp b/src/cosmopower.cpp
index 5e58c2c..9c50588 100644
--- a/src/cosmopower.cpp
+++ b/src/cosmopower.cpp
@@ -249,6 +249,19 @@ double CosmoPower::matterTransferFunctionHu(double k)
   return T_k;
 }
 
+double CosmoPower::noWiggleTk(double k)
+{
+  double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
+  double alpha_Gamma = 1 - 0.328 * log(431 * OmegaEff) * OMEGA_B / OMEGA_0 + 0.38 * log(22.3 * OmegaEff) * pow(OMEGA_B / OMEGA_0, 2);
+  double GammaEff = OMEGA_0 * h * (alpha_Gamma + (1 - alpha_Gamma)/(1 + pow(0.43 * k * s, 4)));
+  double q = k/(h*GammaEff) * pow(Theta_27, 2);
+  double L_0 = log(2 * M_E + 1.8 * q);
+  double C_0 = 14.2 + 731 / (1 + 62.5 * q);
+  double T0 = L_0 / (L_0 + C_0 * q * q);
+
+  return T0;
+}
+
 double CosmoPower::powerHuBaryons(double k)
 {
   double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
@@ -447,6 +460,8 @@ void CosmoPower::setFunction(CosmoFunction f)
     case POWER_BDM:
       eval = &CosmoPower::powerBDM;
       break;
+    case NOWIGGLE_TK:
+      eval = &CosmoPower::noWiggleTk;
     case POWER_TEST:
       eval = &CosmoPower::powerTest;
       break;
diff --git a/src/cosmopower.hpp b/src/cosmopower.hpp
index ba0bffc..768992c 100644
--- a/src/cosmopower.hpp
+++ b/src/cosmopower.hpp
@@ -66,6 +66,10 @@ namespace CosmoTool {
     double Gamma0;
     double normPower;
 
+    double A_BAO;
+    double r_BAO;
+    double k_D_BAO;
+
     struct EHuParams {
       double k_silk;
       double s;
@@ -89,14 +93,14 @@ namespace CosmoTool {
       POWER_SUGIYAMA,
       POWER_BDM,
       POWER_TEST,
-      HU_WIGGLES_ORIGINAL
+      HU_WIGGLES_ORIGINAL,
+      NOWIGGLE_TK
     };
 
     CosmoPower();
     ~CosmoPower();
 
     void setFunction(CosmoFunction f);
-
     void updateCosmology();
     void updatePhysicalCosmology();
     void normalize(double k_min = -1, double k_max = -1);
@@ -108,6 +112,7 @@ namespace CosmoTool {
     double power(double k);
 
     double integrandNormalize(double k);
+
   private:
     double (CosmoPower::*eval)(double);
 
@@ -122,6 +127,8 @@ namespace CosmoTool {
     double powerBDM(double k);
     double powerTest(double k);
     double powerHuWigglesOriginal(double k);
+    double noWiggleTk(double k);
+
   };
 
 };
diff --git a/src/powerSpectrum.cpp b/src/powerSpectrum.cpp
index 777654c..7786080 100644
--- a/src/powerSpectrum.cpp
+++ b/src/powerSpectrum.cpp
@@ -58,6 +58,7 @@ using namespace std;
 #define POWER_TEST 8
 
 #define POWER_SPECTRUM POWER_EFSTATHIOU
+#define SAMPLE_WIGGLES 9
 
 namespace Cosmology {
 
@@ -206,6 +207,31 @@ double powG(double y)
  */
 double powerSpectrum(double k, double normPower)
 {
+#if POWER_SPECTRUM == SAMPLE_WIGGLES
+  // BAO wiggle parameterization for reconstruction
+  // Babic et al. 2022, https://arxiv.org/abs/2203.06177
+
+  // No-wiggle transfer function
+
+  double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
+  double alpha_Gamma = 1 - 0.328 * log(431 * OmegaEff) * OMEGA_B / OMEGA_0 + 0.38 * log(22.3 * OmegaEff) * pow(OMEGA_B / OMEGA_0, 2);
+  double GammaEff = OMEGA_0 * h * (alpha_Gamma + (1 - alpha_Gamma)/(1 + pow(0.43 * k * s, 4)));
+  double q = k/(h*GammaEff) * pow(Theta_27, 2);
+  double L_0 = log(2 * M_E + 1.8 * q);
+  double C_0 = 14.2 + 731 / (1 + 62.5 * q);
+  double T0 = L_0 / (L_0 + C_0 * q * q);
+
+  // Wiggle parameterization
+
+  double A = 0;
+  double r_s = 10;
+  double k_D = 2 * M_PI / 100;
+
+  double param = 1 + A * sin(k * r_s) * exp(- k / k_D);
+
+  return normPower * pow(k, n) * T0 * T0 * param;
+#endif
+
 #if POWER_SPECTRUM == POWER_EFSTATHIOU
   double a = 6.4/Gamma0;
   double b = 3/Gamma0;