From 06ed21cae00e1083e0e22cae3f005079799fddec Mon Sep 17 00:00:00 2001
From: Oliver Hahn <oliverjhahn@gmail.com>
Date: Sun, 19 May 2019 12:05:04 +0200
Subject: [PATCH] on the fly output of power spectra

---
 include/grid_fft.hh |  6 ++++--
 src/grid_fft.cc     | 37 +++++++++++++++++++++++++++++++++----
 src/main.cc         |  9 +++++++++
 3 files changed, 46 insertions(+), 6 deletions(-)

diff --git a/include/grid_fft.hh b/include/grid_fft.hh
index 0ad191b..36a6e30 100644
--- a/include/grid_fft.hh
+++ b/include/grid_fft.hh
@@ -430,9 +430,11 @@ public:
 
     void Write_to_HDF5(std::string fname, std::string datasetname);
 
-    void ComputePowerSpectrum(std::vector<double> &bin_k, std::vector<double> &bin_P, std::vector<double> &bin_eP, int nbins);
+    void Write_PowerSpectrum( std::string ofname );
 
-    void Compute_PDF(std::string ofname, int nbins = 1000, double scale = 1.0, double rhomin = 1e-3, double rhomax = 1e3);
+    void Compute_PowerSpectrum(std::vector<double> &bin_k, std::vector<double> &bin_P, std::vector<double> &bin_eP, std::vector<size_t> &bin_count, int nbins);
+
+    void Write_PDF(std::string ofname, int nbins = 1000, double scale = 1.0, double rhomin = 1e-3, double rhomax = 1e3);
 
     void zero_DC_mode(void)
     {
diff --git a/src/grid_fft.cc b/src/grid_fft.cc
index 0b4d621..6c58218 100644
--- a/src/grid_fft.cc
+++ b/src/grid_fft.cc
@@ -501,7 +501,7 @@ void Grid_FFT<data_t>::Write_to_HDF5(std::string fname, std::string datasetname)
 #include <iomanip>
 
 template <typename data_t>
-void Grid_FFT<data_t>::Compute_PDF( std::string ofname, int nbins, double scale, double vmin, double vmax )
+void Grid_FFT<data_t>::Write_PDF( std::string ofname, int nbins, double scale, double vmin, double vmax )
 {
     double logvmin = std::log10(vmin);
     double logvmax = std::log10(vmax);
@@ -551,15 +551,44 @@ void Grid_FFT<data_t>::Compute_PDF( std::string ofname, int nbins, double scale,
 }
 
 template <typename data_t>
-void Grid_FFT<data_t>::ComputePowerSpectrum(std::vector<double> &bin_k, std::vector<double> &bin_P, std::vector<double> &bin_eP, int nbins)
+void Grid_FFT<data_t>::Write_PowerSpectrum( std::string ofname )
+{
+    std::vector<double> bin_k, bin_P, bin_eP;
+    std::vector<size_t> bin_count;
+    int nbins = 4*std::max( nhalf_[0], std::max(nhalf_[1], nhalf_[2]) );
+    this->Compute_PowerSpectrum( bin_k, bin_P, bin_eP, bin_count, nbins );
+#if defined(USE_MPI)
+    if (CONFIG::MPI_task_rank == 0)
+    {
+#endif
+        std::ofstream ofs( ofname.c_str());
+        std::size_t numcells = size(0) * size(1) * size(2);
+
+        //ofs << "# a = " << aexp << std::endl;
+        ofs << "# " << std::setw(14) << "k" << std::setw(16) << "P(k)" << std::setw(16) << "err. P(k)"
+            << std::setw(16) <<"#modes" << "\n";
+
+        for( int ibin=0; ibin<nbins; ++ibin ){
+            if( bin_count[ibin] > 0 )
+            ofs << std::setw(16) << bin_k[ibin]
+                << std::setw(16) << bin_P[ibin]
+                << std::setw(16) << bin_eP[ibin]
+                << std::setw(16) << bin_count[ibin]
+                << std::endl;
+        }
+#if defined(USE_MPI)
+    }
+#endif
+}
+
+template <typename data_t>
+void Grid_FFT<data_t>::Compute_PowerSpectrum(std::vector<double> &bin_k, std::vector<double> &bin_P, std::vector<double> &bin_eP, std::vector<size_t> &bin_count, int nbins)
 {
     real_t kmax = std::max(std::max(kfac_[0] * nhalf_[0], kfac_[1] * nhalf_[1]),
                            kfac_[2] * nhalf_[2]),
            kmin = std::min(std::min(kfac_[0], kfac_[1]), kfac_[2]),
            dklog = log10(kmax / kmin) / nbins;
 
-    std::vector<size_t> bin_count;
-
     bin_count.assign(nbins,0);
     bin_k.assign(nbins, 0);
     bin_P.assign(nbins, 0);
diff --git a/src/main.cc b/src/main.cc
index e0220cf..551f6c0 100644
--- a/src/main.cc
+++ b/src/main.cc
@@ -15,6 +15,7 @@
 #include <output_plugin.hh>
 #include <cosmology_calculator.hh>
 
+// initialise with "default" values
 namespace CONFIG{
 int  MPI_thread_support = -1;
 int  MPI_task_rank = 0;
@@ -94,6 +95,7 @@ int main( int argc, char** argv )
 
     //...
     const std::string fname_hdf5 = the_config.GetValueSafe<std::string>("output", "fname_hdf5", "output.hdf5");
+    const std::string fname_analysis = the_config.GetValueSafe<std::string>("output", "fbase_analysis", "output");
     //////////////////////////////////////////////////////////////////////////////////////////////
 
     std::unique_ptr<CosmologyCalculator>  the_cosmo_calc;
@@ -276,6 +278,13 @@ int main( int argc, char** argv )
                 }
             }
         }
+
+        delta.Write_PowerSpectrum(fname_analysis+"_"+"power_delta1.txt");
+        delta2.Write_PowerSpectrum(fname_analysis+"_"+"power_delta2.txt");
+        delta3a.Write_PowerSpectrum(fname_analysis+"_"+"power_delta3a.txt");
+        delta3b.Write_PowerSpectrum(fname_analysis+"_"+"power_delta3b.txt");
+        delta3.Write_PowerSpectrum(fname_analysis+"_"+"power_delta3.txt");
+
         phi.FourierTransformBackward();
         phi2.FourierTransformBackward();
         phi3a.FourierTransformBackward();