Upgraded gadget load to support variable masses. Adjusted simple3DFilter for that too.

python_sample/icgen/borgicgen.py

@@ -156,7 +156,7 @@ def whitify(density, L, cosmo, supergenerate=1, func='HU_WIGGLES'):
 
     print np.where(np.isnan(density_hat_super))[0].size
 
-  cube = CubeFT(Ns, L)
+  cube = CubeFT(L, Ns)
   cube.dhat = density_hat_super
   return np.fft.irfftn(density_hat_super)*Ns**1.5
 

python_sample/icgen/gen_ic_from_borg.py

@@ -9,10 +9,10 @@ cosmo['omega_B_0']=0.049
 cosmo['SIGMA8']=0.8344
 cosmo['ns']=0.9624
 
-supergen=2
+supergen=1
 zstart=50
 astart=1/(1.+zstart)
 halfPixelShift=False
 
 if __name__=="__main__":
-    bic.write_icfiles(*bic.run_generation("initial_condition_borg.dat", 0.001, astart, cosmo, supersample=2, shiftPixel=halfPixelShift, do_lpt2=False), supergenerate=supergen)
+    bic.write_icfiles(*bic.run_generation("initial_density_988.dat", 0.001, astart, cosmo, supersample=1, shiftPixel=halfPixelShift, do_lpt2=False), supergenerate=supergen)

sample/gadgetToArray.cpp

@@ -50,9 +50,11 @@ int main(int argc, char **argv)
   uint32_t res;
   char *fname;
   int id;
+  double MPC;
 
   MiniArgDesc desc[] = {
     { "SNAPSHOT", &fname, MINIARG_STRING },
+    { "MPC", &MPC, MINIARG_DOUBLE },
     { 0, 0, MINIARG_NULL }
   };
 
@@ -63,7 +65,7 @@ int main(int argc, char **argv)
   
 
   SimuData *p = loadGadgetMulti(fname, 0, 0);
-  double L0 = p->BoxSize/1000;
+  double L0 = p->BoxSize/MPC;
   array_type parts(boost::extents[p->TotalNumPart][7]);
   uint64_t q = 0;
   
@@ -72,13 +74,12 @@ int main(int argc, char **argv)
       p = loadGadgetMulti(fname, cpuid, NEED_POSITION|NEED_VELOCITY|NEED_MASS);
       for (uint32_t i = 0; i < p->NumPart; i++)
         {
-          parts[q][0] = p->Pos[0][i]/1000;
-          parts[q][1] = p->Pos[1][i]/1000;
-          parts[q][2] = p->Pos[2][i]/1000;
           for (int j = 0; j < 3; j++)
             {
+              parts[q][j] = p->Pos[j][i]/MPC;
               while (parts[q][j] < 0) parts[q][j] += L0;
               while (parts[q][j] >= L0) parts[q][j] -= L0;
+              parts[q][j] -= L0/2;
             }
           parts[q][3] = p->Vel[0][i];
           parts[q][4] = p->Vel[1][i];

sample/simple3DFilter.cpp

@@ -114,7 +114,7 @@ int main(int argc, char **argv)
   rLimit2 = rLimit*rLimit;
 
   hdf5_read_array(in_f, "particles", v1_data);
-  assert(v1_data.shape()[1] == 6);
+  assert(v1_data.shape()[1] == 7);
 
   N1_points = v1_data.shape()[0];
   
@@ -132,6 +132,7 @@ int main(int argc, char **argv)
         allCells_1[i].coord[j] = v1_data[i][j];
       for (int k = 0; k < 3; k++)
         allCells_1[i].val.pValue.v[k] = v1_data[i][3+k];
+      allCells_1[i].val.pValue.mass = v1_data[i][6];
       allCells_1[i].active = true;
       allCells_1[i].val.weight = 0.0;
 

src/loadGadget.cpp

@@ -262,9 +262,17 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
   }
 
   if (loadflags & NEED_MASS) {
+    bool do_load = false;
+
+    for (int k = 0; k < 6; k++)
+      {
+         do_load = do_load || ((h.mass[k] == 0)&&(h.npart[k]>0));
+      }
+
     try
       {
         long l = 0;
+	if (do_load)
           f->beginCheckpoint();
         data->Mass = new float[NumPart];
         for (int k = 0; k < 6; k++)
@@ -281,6 +289,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
                 }
             }
           }
+        if (do_load)
           f->endCheckpoint();
       }
     catch (const InvalidUnformattedAccess& e)
@@ -290,6 +299,11 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
 	delete data;
 	return 0;
       }
+    catch (const EndOfFileException& e)
+      {
+        for (int k = 0; k < 6; k++)
+          cerr << "mass[" << k << "] = " << h.mass[k] << endl;
+      }
   } else {
     f->skip(2*4);
     for (int k = 0; k < 6; k++)