From eaf6dbd74d4470d99831fc974a458a84e93d0f0e Mon Sep 17 00:00:00 2001
From: Chichi Lalescu <chichilalescu@gmail.com>
Date: Tue, 29 May 2018 09:46:22 +0200
Subject: [PATCH] add spectrum plot

---
 bfps/test/test_Parseval.py | 17 ++++++++++++++++-
 1 file changed, 16 insertions(+), 1 deletion(-)

diff --git a/bfps/test/test_Parseval.py b/bfps/test/test_Parseval.py
index 6bc19574..52552628 100644
--- a/bfps/test/test_Parseval.py
+++ b/bfps/test/test_Parseval.py
@@ -8,6 +8,8 @@ import sys
 import bfps
 from bfps import DNS
 
+import matplotlib.pyplot as plt
+
 def main():
     niterations = 10
     c = DNS()
@@ -25,7 +27,20 @@ def main():
              '--wd', './'] +
              sys.argv[1:])
     c.compute_statistics()
-    print((c.statistics['energy(t)'] - c.statistics['renergy(t)']) / c.statistics['renergy(t)'])
+    print((c.statistics['energy(t)']*3 - c.statistics['renergy(t)']) / c.statistics['renergy(t)'])
+    print(list(c.statistics.keys()))
+    energyk = np.mean(c.statistics['energy(t, k)'], axis = 0)
+    nshell = c.get_data_file()['kspace/nshell'].value
+    renergy = np.mean(c.statistics['renergy(t)'])
+    print(renergy, np.trapz(energyk[:-2], c.statistics['kshell'][:-2]))
+
+    f = plt.figure()
+    a = f.add_subplot(111)
+    a.plot(c.statistics['kshell'], energyk)
+    a.plot(c.statistics['kshell'], (energyk / nshell)*(4*np.pi*c.statistics['kshell']**2))
+    a.set_yscale('log')
+    a.set_xscale('log')
+    f.savefig('spectrum.pdf')
     return None
 
 if __name__ == '__main__':
-- 
GitLab