Use long instead of int for particle identifiers

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})
+

src/cosmopower.cpp

@@ -191,15 +191,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);
@@ -210,8 +210,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);
@@ -226,6 +233,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)
@@ -282,3 +299,8 @@ void CosmoPower::setFunction(CosmoFunction f)
       abort();
     }
 }
+
+void CosmoPower::setNormalization(double A_K)
+{
+  normPower = A_K/power(0.002);
+}

src/cosmopower.hpp

@@ -17,6 +17,8 @@ namespace CosmoTool {
     double SIGMA8;
     double OMEGA_B;
     double OMEGA_C;
+    double omega_B;
+    double omega_C;
     double Theta_27;
     
     // DERIVED VARIABLES    
@@ -44,7 +46,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

@@ -1,6 +1,7 @@
 #ifndef __DETAILS_EUCLIDIAN_SPECTRUM_1D
 #define __DETAILS_EUCLIDIAN_SPECTRUM_1D
 
+#include  <iostream>
 #include <boost/function.hpp>
 
 

src/h5_readFlash.cpp

@@ -185,7 +185,7 @@ void h5_read_flash3_particles (H5File* file,
                         float *vel1,
                         float *vel2,
                         float *vel3,
-                        int    *id)
+                        long    *id)
 {
 
   herr_t   status;
@@ -341,7 +341,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

@@ -29,7 +29,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

@@ -71,7 +71,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

@@ -195,7 +195,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

@@ -384,7 +384,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

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