Merge branch 'master' of /home/guilhem/Dropbox/gitRoot/CosmoToolbox

sample/CMakeLists.txt

@@ -1,5 +1,5 @@
 SET(tolink CosmoTool ${GSL_LIBRARY} ${GSLCBLAS_LIBRARY})
-include_directories(${CMAKE_SOURCE_DIR}/src)
+include_directories(${CMAKE_SOURCE_DIR}/src ${NETCDF_INCLUDE_PATH} ${GSL_INCLUDE_PATH})
 
 add_executable(testBQueue testBQueue.cpp)
 target_link_libraries(testBQueue ${tolink})

sample/scale.cpp

@@ -1,20 +0,0 @@
-#include "yorick.hpp"
-
-using namespace CosmoTool;
-
-int main(int argc, char **argv)
-{
-  uint32_t *dimList;
-  uint32_t dimRank;
-  ProgressiveInput<float> input = ProgressiveInput<float>::loadArrayProgressive("displacement.nc", dimList, dimRank);
-  ProgressiveOutput<float> output = ProgressiveOutput<float>::saveArrayProgressive("scaledDispl.nc", dimList, dimRank);
-
-  uint32_t N = 3 * dimList[1] * dimList[2] * dimList[3];
-
-  for (uint32_t i = 0; i < N; i++)
-   {
-     output.put(input.read() * 500./.65);
-   } 
-
-  return 0;
-}

src/CMakeLists.txt

@@ -8,6 +8,7 @@ SET(CosmoTool_SRCS
 	powerSpectrum.cpp
 	yorick.cpp
 	miniargs.cpp
+	growthFactor.cpp
 )
 
 SET(CosmoTool_SRCS ${CosmoTool_SRCS}
@@ -31,8 +32,10 @@ SET(CosmoTool_SRCS ${CosmoTool_SRCS}
 	sparseGrid.hpp
 	sphSmooth.hpp
 	yorick.hpp
+	growthFactor.hpp
 )
 
+include_directories(${GSL_INCLUDE_PATH} ${NETCDF_INCLUDE_PATH})
 add_library(CosmoTool ${CosmoTool_SRCS})
 target_link_libraries(CosmoTool ${NETCDF_LIBRARY} ${NETCDFCPP_LIBRARY} ${GSL_LIBRARY} ${GSLCBLAS_LIBRARY})
 

src/growthFactor.cpp

@@ -0,0 +1,111 @@
+#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;
+  
+  a_input[0] = 0;
+  D_output[0] = 0;
+
+  for (int i = 1; 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

src/interpolate.cpp

@@ -125,7 +125,7 @@ Interpolate CosmoTool::buildInterpolateFromFile(const char *fname)
     {
       MyPair m;
 
-      if ((f >> m.a >> m.b).eof())
+      if (!(f >> m.a >> m.b))
 	break;
 
       allData.push_back(m);

src/octTree.hpp

@@ -54,9 +54,7 @@ namespace CosmoTool
     template<typename FunT>
     void walkTree(FunT f)
     {
-      OctCoords rootCenter = { octCoordCenter, octCoordCenter, octCoordCenter };
-
-      walkTreeElements(f, unconditioned, 0, rootCenter, octCoordCenter);
+      walkTree(f, unconditioned);
     }
 
     template<typename FunT, typename CondT>
@@ -80,8 +78,14 @@ namespace CosmoTool
     float xMin[3];
 
 
-    template<typename FunT,typename CondT>
-    void walkTreeElements(FunT f, CondT condition, octPtr node, 
+    static bool unconditioned()
+    {
+      return true;
+    }
+
+    template<typename FunT, typename CondT>
+    void walkTreeElements(FunT f, CondT condition,
+			  octPtr node, 
 			  const OctCoords& icoord,
 			  octCoordType halfNodeLength)
     {
@@ -127,6 +131,7 @@ namespace CosmoTool
 		false, true);
 	      continue;
 	    }
+
 	  walkTreeElements(f, condition, cells[node].children[i], newCoord, halfNodeLength/2);
 	}