Merge branch 'master' of bitbucket.org:glavaux/cosmotool

sample/CMakeLists.txt

@@ -55,3 +55,7 @@ if (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
   target_link_libraries(test_healpix_calls ${tolink} ${SHARP_LIBRARIES})
   set_target_properties(test_healpix_calls PROPERTIES COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS})
 endif (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
+
+add_executable(test_cosmopower test_cosmopower.cpp)
+target_link_libraries(test_cosmopower ${tolink})
+

sample/test_cosmopower.cpp

@@ -0,0 +1,29 @@
+#include <cmath>
+#include <iostream>
+#include "cosmopower.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+
+int main()
+{
+   CosmoPower cp;
+   CosmoPower::CosmoFunction f[] = { CosmoPower::POWER_EFSTATHIOU, CosmoPower::HU_WIGGLES, CosmoPower::HU_BARYON, CosmoPower::POWER_SUGIYAMA };
+   int num_F = sizeof(f)/sizeof(f[0]);
+
+   cp.setFunction(f[0]);
+   cp.normalize();
+   for (int ik = 0; ik < 100; ik++)
+    {
+      double k = pow(10.0, 4*ik/100.-3);
+      
+      cout << k << " ";
+      for (int q = 0; q < num_F; q++)
+        {
+          cp.setFunction(f[q]);
+          cout << cp.power(k) << " ";
+        } 
+      cout << endl;
+    }
+   return 0;
+}

src/cosmopower.cpp

@@ -225,15 +225,15 @@ static double gslPowSpecNorm(double k, void *params)
   return c->integrandNormalize(k);
 }
 
-double CosmoPower::integrandNormalize(double k)
+double CosmoPower::integrandNormalize(double x)
 {
+  double k = (1-x)/x;
   double f = tophatFilter(k*8.0/h);
-  return power(k)*k*k*f*f;
+  return power(k)*k*k*f*f/(x*x);
 }
 
 void CosmoPower::normalize()
 {
-  int Nsteps = 30000;
   double normVal = 0;
   double abserr;
   gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION);
@@ -244,8 +244,15 @@ void CosmoPower::normalize()
 
   normPower = 1;
 
-  gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &normVal, &abserr);
+  ofstream ff("PP_k.txt");
+  for (int i = 0; i < 100; i++)
+    {
+       double k = pow(10.0, 4.0*i/100.-2);
+       ff << k << " " << power(k) << endl;
+    }
 
+ // gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &normVal, &abserr);
+  gsl_integration_qag(&f, 0, 1, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &normVal, &abserr);
   gsl_integration_workspace_free(w);
 
   normVal /= (2*M_PI*M_PI);
@@ -260,6 +267,16 @@ void CosmoPower::updateCosmology()
   beta = pow(OMEGA_0, 5./9);
   OmegaEff = OMEGA_0 * h * h;
   Gamma0 = OMEGA_0 * h * h;
+  omega_B = OMEGA_B * h * h;
+  omega_C = OMEGA_C * h * h;
+}
+
+void CosmoPower::updatePhysicalCosmology()
+{
+  OMEGA_B = omega_B / (h*h);
+  OMEGA_C = omega_C / (h*h);
+  OMEGA_0 = Gamma0 / (h*h);
+  beta = pow(OMEGA_0, 5./9);
 }
 
 double CosmoPower::eval_theta_theta(double k)
@@ -316,3 +333,8 @@ void CosmoPower::setFunction(CosmoFunction f)
       abort();
     }
 }
+
+void CosmoPower::setNormalization(double A_K)
+{
+  normPower = A_K/power(0.002);
+}

src/cosmopower.hpp

@@ -51,6 +51,8 @@ namespace CosmoTool {
     double SIGMA8;
     double OMEGA_B;
     double OMEGA_C;
+    double omega_B;
+    double omega_C;
     double Theta_27;
     
     // DERIVED VARIABLES    
@@ -78,7 +80,9 @@ namespace CosmoTool {
     void setFunction(CosmoFunction f);
 
     void updateCosmology();
+    void updatePhysicalCosmology();
     void normalize();
+    void setNormalization(double A_K);
 
  
     double eval_theta_theta(double k);

src/fourier/details/euclidian_spectrum_1d.hpp

@@ -35,6 +35,7 @@ knowledge of the CeCILL license and that you accept its terms.
 #ifndef __DETAILS_EUCLIDIAN_SPECTRUM_1D
 #define __DETAILS_EUCLIDIAN_SPECTRUM_1D
 
+#include  <iostream>
 #include <boost/function.hpp>
 
 

src/h5_readFlash.cpp

@@ -219,7 +219,7 @@ void h5_read_flash3_particles (H5File* file,
                         float *vel1,
                         float *vel2,
                         float *vel3,
-                        int    *id)
+                        long    *id)
 {
 
   herr_t   status;
@@ -375,7 +375,7 @@ void h5_read_flash3_particles (H5File* file,
 
     if (id) {
     for(p=0; p < (pcount); p++) {
-      id[p+poffset] = (int)  *(partBuffer+iptag-1+p*numProps);
+      id[p+poffset] = (long)  *(partBuffer+iptag-1+p*numProps);
     } }
 
     if (pos1 && pos2 && pos3) {

src/h5_readFlash.hpp

@@ -63,7 +63,7 @@ void h5_read_flash3_particles (H5File* file,
                         float *vel1,
                         float *vel2,
                         float *vel3,
-                        int    *id);
+                        long    *id);
 
 void h5_read_flash3_header_info(H5File* file,
 				double* time,                   /* simulation time */

src/loadFlash.cpp

@@ -105,7 +105,7 @@ SimuData *CosmoTool::loadFlashMulti(const char *fname, int id, int loadflags)
     } }
  
     if (loadflags & NEED_GADGET_ID) {
-    data->Id = new int[data->NumPart];
+    data->Id = new long[data->NumPart];
     if (data->Id == 0) {
       delete data;
       return 0;

src/loadGadget.cpp

@@ -229,7 +229,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
     try
       {
 	f->beginCheckpoint();
-	data->Id = new int[data->NumPart];
+	data->Id = new long[data->NumPart];
 	if (data->Id == 0)
 	  {
 	    delete f;

src/loadRamses.cpp

@@ -418,7 +418,7 @@ CosmoTool::SimuData *CosmoTool::loadRamsesSimu(const char *basename, int outputI
 	}
       if (flags & NEED_GADGET_ID)
 	{
-	  data->Id = new int[nPar];
+	  data->Id = new long[nPar];
 	}
 
       for (int k = 0; k < 3; k++)

src/loadSimu.hpp

@@ -72,7 +72,7 @@ namespace CosmoTool
 
     long NumPart;
     long TotalNumPart;
-    int *Id;
+    long *Id;
     float *Pos[3];
     float *Vel[3];
     int *type;