From 5b133a9ac839ae9a0f57f2e56c289b79ba52fc0a Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Tue, 10 Jun 2014 14:18:14 +0200
Subject: [PATCH] Increased SPH number. Changed defaults in python samples

---
 python_sample/build_nbody_skymap.py      | 54 ++++++++++++++++++++----
 python_sample/icgen/test_ic_from_borg.py |  2 +-
 sample/simple3DFilter.cpp                |  2 +-
 3 files changed, 48 insertions(+), 10 deletions(-)

diff --git a/python_sample/build_nbody_skymap.py b/python_sample/build_nbody_skymap.py
index 1aa948a..2569a05 100644
--- a/python_sample/build_nbody_skymap.py
+++ b/python_sample/build_nbody_skymap.py
@@ -1,19 +1,57 @@
+import healpy as hp
+import numpy as np
 import cosmotool as ct
 import h5py as h5
+from matplotlib import pyplot as plt
 
 L=600.
 Nside=128
 
-with h5.File("fields.h5", mode="r") as f:
-    density = f["density"][:]
+INDATA="/nethome/lavaux/Copy/PlusSimulation/BORG/Input_Data/2m++.npy"
+tmpp = np.load(INDATA)
 
-N = density.shape[0]
-ix = (np.arange(N)-0.5)*L/N - 0.5 * L
+def build_sky_proj(density, dmax=60.,dmin=0):
+
+  N = density.shape[0]
+  ix = (np.arange(N)-0.5)*L/N - 0.5 * L
 
 
-dist2 = (ix[:,None,None]**2 + ix[None,:,None]**2 + ix[None,None,:]**2)
+  dist2 = (ix[:,None,None]**2 + ix[None,:,None]**2 + ix[None,None,:]**2)
 
-flux = density / dist2
-projsky1 = ct.spherical_projection(Nside, flux, 0, 52, integrator_id=1)
-projsky0 = ct.spherical_projection(Nside, flux, 0, 52, integrator_id=0)
+  flux = density.transpose().astype(ct.DTYPE) # / dist2
+  dmax=N*dmax/L
+  dmin=N*dmin/L
+  projsky1 = ct.spherical_projection(Nside, flux, dmin, dmax, integrator_id=1)
+#  projsky0 = ct.spherical_projection(Nside, flux, 0, 52, integrator_id=0)
 
+  return projsky1*L/N#,projsky0
+
+l,b = tmpp['gal_long'],tmpp['gal_lat']
+
+l = np.radians(l)
+b = np.pi/2 - np.radians(b)
+
+dcmb = tmpp['velcmb']/100.
+
+idx = np.where((dcmb>10)*(dcmb<60))
+
+plt.figure(1)
+plt.clf()
+if False:
+  with h5.File("fields.h5", mode="r") as f:
+    d = f["density"][:]
+    d /= np.average(np.average(np.average(d,axis=0),axis=0),axis=0)
+    proj = build_sky_proj(f["density"][:], dmin=10,dmax=60.)
+else:
+  d = np.load("icgen/dcic0.npy")
+  proj0 = build_sky_proj(1+d, dmin=10,dmax=60.)
+  d = np.load("icgen/dcic1.npy")
+  proj1 = build_sky_proj(1+d, dmin=10,dmax=60.)
+
+hp.mollview(proj0, fig=1, coord='CG', max=60, cmap=plt.cm.coolwarm)
+hp.projscatter(b[idx], l[idx], lw=0, color='g', s=5.0, alpha=0.8)
+
+plt.figure(2)
+plt.clf()
+hp.mollview(proj1, fig=2, coord='CG', max=60, cmap=plt.cm.coolwarm)
+hp.projscatter(b[idx], l[idx], lw=0, color='g', s=5.0, alpha=0.8)
diff --git a/python_sample/icgen/test_ic_from_borg.py b/python_sample/icgen/test_ic_from_borg.py
index 7f069a8..e04c26c 100644
--- a/python_sample/icgen/test_ic_from_borg.py
+++ b/python_sample/icgen/test_ic_from_borg.py
@@ -50,5 +50,5 @@ Pref, bref = bic.compute_ref_power(L, N0, cosmo, range=(0,1.), bins=150)
 Pcic /= D1_0**2
 Pdens /= D1_0**2
 
-borg_evolved = ct.read_borg_vol("final_density_1380.dat")
+borg_evolved = ct.read_borg_vol("final_density_380.dat")
 
diff --git a/sample/simple3DFilter.cpp b/sample/simple3DFilter.cpp
index c333f2e..b3a74d2 100644
--- a/sample/simple3DFilter.cpp
+++ b/sample/simple3DFilter.cpp
@@ -11,7 +11,7 @@
 using namespace std;
 using namespace CosmoTool;
 
-#define N_SPH 16
+#define N_SPH 32
 
 struct VCoord{
   float v[3];