/*
 * site_writePotential_zphi.cpp
 *
 *  Created on: Feb 17, 2016
 *      Author: Argo
 */

#include "site.h"
#include "math.h"
#include <vector>
#include <complex>
#include "../delphi/delphi_constants.h"


void CSite::writeZetaPhiFile()
{
    size_t MAXWIDTH = 45;
    vector<delphi_real> surfPots;      
    cout << left << setw(MAXWIDTH) << " Surface> Writing Surface Potential values now to "
         << strZetaPhiFile << endl;

    ofstream zetaPhi_file;
    zetaPhi_file.open(strZetaPhiFile);

    int i=0,j=0,k=0,pidx=0,cidx=0;  // ARGO the grid indices (i,j,k)from the zetaSite-file
                                    // written in the space module
    int m=6;                        // num of theta partitions/2
    int n=12;                       // num of phi partitions/2

    //AGO the real-time coords for the zetaSite-file coordinates from space module
    delphi_real px=0.,py=0.,pz=0.,cosine=0.,phi=0.,sine =0.;

    delphi_real simple_avg_SurfPot = 0;

    vector<delphi_real>::iterator coord_itr = surf_grid_coords_v.begin();
    vector<delphi_real>::iterator index_itr = surf_grid_index_v.begin();

    if ( bZetaPhiOut )
    {
        while (index_itr != surf_grid_index_v.end())
        {
            i = *index_itr;
            j = *(index_itr + 1);
            k = *(index_itr + 2);

            index_itr += 3;
            surfPots.push_back(phimap[k][j][i]);
        }

        zetaPhi_file << "#REMARK IMPORTANT INFORMATION AT THE BOTTOM OF THIS FILE" << endl;

        vector<delphi_real>::iterator surfPot_itr = surfPots.begin();

        while ( coord_itr != surf_grid_coords_v.end() )
        {
            px = *coord_itr;
            py = *(coord_itr + 1);
            pz = *(coord_itr + 2);

            zetaPhi_file << setw(10) << setprecision(6) << px
                         << setw(10) << setprecision(6) << py
                         << setw(10) << setprecision(6) << pz
                         << setw(15) << setprecision(6) << *(surfPot_itr)
                         << endl;

            coord_itr += 3;
            surfPot_itr++;
        }

        zetaPhi_file << "# -------------------------------------------------------" << endl;
        // mean RADIAL distance
        // Rh = meanPotential(radial_dist_v);

        // mean surface potential
        simple_avg_SurfPot = meanPotential(surfPots);

        zetaPhi_file << left << setw(MAXWIDTH) << "#REMARK SIMPLE AVERAGE SURFACE POTENTIAL "
                     << " = " << simple_avg_SurfPot << " kT/e" << endl;
        zetaPhi_file << left << setw(MAXWIDTH) << "#REMARK GEOMETRIC CENTER (ANG) " << " = "
                     << fgBoxCenter.nX << "  " << fgBoxCenter.nY << "  " << fgBoxCenter.nZ << endl;
        zetaPhi_file.close();
    }
}

delphi_real CSite::meanPotential ( vector<delphi_real>& pots )
{
    delphi_real mn = 0;
    vector<delphi_real>::iterator itpm = pots.begin();

    while ( itpm != pots.end() )
    {
        mn += *itpm;
        itpm++;
    }

    mn /= pots.size();
    return mn;
}
