Implemented standard linear growth factor computation

src/CMakeLists.txt

@@ -8,6 +8,7 @@ SET(CosmoTool_SRCS
 	powerSpectrum.cpp
 	yorick.cpp
 	miniargs.cpp
+	growthFactor.cpp
 )
 
 SET(CosmoTool_SRCS ${CosmoTool_SRCS}
@@ -31,6 +32,7 @@ SET(CosmoTool_SRCS ${CosmoTool_SRCS}
 	sparseGrid.hpp
 	sphSmooth.hpp
 	yorick.hpp
+	growthFactor.hpp
 )
 
 add_library(CosmoTool ${CosmoTool_SRCS})

src/growthFactor.cpp

@@ -0,0 +1,108 @@
+#include <cmath>
+#include <gsl/gsl_integration.h>
+#include "interpolate.hpp"
+#include "growthFactor.hpp"
+
+using namespace CosmoTool;
+
+#define AMIN 1e-5
+#define AMAX 1.0
+#define NUM_WORK 5000
+#define TOLERANCE 1e-6
+
+typedef struct {
+  double OmegaLambda;
+  double OmegaMatter;
+  double Hubble;
+} Cosmology;
+
+static double computeOmegaMatter(Cosmology *cosmo, double a)
+{
+  return cosmo->OmegaMatter / (cosmo->OmegaMatter + a*a*a * cosmo->OmegaLambda);
+}
+
+static double computeHdotH(Cosmology *cosmo, double a)
+{
+  return -1.5 * cosmo->OmegaMatter / (a * (cosmo->OmegaMatter + a*a*a*cosmo->OmegaLambda));
+}
+
+static double computeE(double OmegaMatter, double OmegaLambda, double a)
+{
+  double H2;
+  double OmegaK = (1 - OmegaMatter - OmegaLambda);
+  
+  H2 = OmegaMatter/(a*a*a) + OmegaLambda + OmegaK/(a*a);
+
+  return sqrt(H2);
+}
+
+static double computeEprime(Cosmology *cosmo, double a)
+{
+  double H2;
+  double OmegaK = (1 - cosmo->OmegaMatter - cosmo->OmegaLambda);
+  
+  H2 = -3*cosmo->OmegaMatter/(a*a*a*a) - 2*OmegaK/(a*a*a);
+
+  return 0.5*H2/computeE(cosmo->OmegaMatter, cosmo->OmegaLambda, a);  
+}
+
+static inline double cube(double x)
+{
+  return x*x*x;
+}
+
+static double integrandGrowthFactor(double a, void *params)
+{
+  Cosmology *cosmo = (Cosmology *)params;
+
+  return 1/cube(computeE(cosmo->OmegaMatter, cosmo->OmegaLambda, a)*a);
+}
+
+Interpolate CosmoTool::buildLinearGrowth(double OmegaLambda, double OmegaMatter, double Hubble, int numPoints)
+{
+  Cosmology cosmology;
+  gsl_integration_workspace *work = gsl_integration_workspace_alloc(NUM_WORK);
+  gsl_function f;
+  double *a_input, *D_output;
+
+  cosmology.OmegaLambda = OmegaLambda;
+  cosmology.OmegaMatter = OmegaMatter;
+  cosmology.Hubble = Hubble;
+
+  a_input = new double[numPoints];
+  D_output = new double[numPoints];
+
+  f.params = &cosmology;
+  f.function = integrandGrowthFactor;
+  
+  for (int i = 0; i < numPoints; i++)
+    {
+      double a_dest = 0 + 1.0*i/(numPoints-1);
+      double result, abserr;
+      double E = computeE(cosmology.OmegaMatter, cosmology.OmegaLambda, a_dest);
+      double Eprime = computeEprime(&cosmology, a_dest);
+
+      gsl_integration_qag(&f, 0, a_dest, 0, TOLERANCE, NUM_WORK, 
+			  GSL_INTEG_GAUSS61, work, &result, &abserr);
+
+
+      result *= 2.5 * computeE(cosmology.OmegaMatter, cosmology.OmegaLambda, a_dest) * OmegaMatter;
+
+      D_output[i] = result;
+      a_input[i] = a_dest;
+
+    }
+  gsl_integration_workspace_free(work);
+
+  for (int i = 0; i < numPoints; i++)
+    {
+      D_output[i] /= D_output[numPoints-1];
+    }
+
+  Interpolate p(a_input, D_output, numPoints, true, false, true);
+
+  delete[] a_input;
+  delete[] D_output;
+
+  return p;
+}

src/growthFactor.hpp

@@ -0,0 +1,12 @@
+#ifndef COSMO_GROWTH_FACTOR_HPP
+#define COSMO_GROWTH_FACTOR_HPP
+
+#include "interpolate.hpp"
+
+namespace CosmoTool
+{
+  Interpolate buildLinearGrowth(double OmegaLambda, double OmegaMatter, double Hubble, int numPoints = 10000);
+
+};
+
+#endif