//
// Program: cond_piano.cxx
//
// Description: Evaluate the potential and electric fields for a parallel
//              plate condenser.
//
// Author: A. Garfagnini
//
// Date: 16 October 2013
//

#include <iostream>
#include <fstream>
#include <iomanip>
#include <unistd.h>
#include <libgen.h>
#include <stdlib.h>
#include <cmath>

using namespace std;

bool debug = false;

void usage(const char * pname)
{
    cerr << "Usage: " << pname << " [-d] [-o file_name]" << endl;
    cerr << "\nEvaluate the Potential for a parallel plate condenser";
    cerr << "\nusing the Gauss-Seidel method\n\nOptions:";
    cerr << "\n -d : enable debug printout";
    cerr << "\nOther options:\n";
    cerr << " -o file_name : save output in a user definable file name\n";
}

inline int index(int i, int j, int n) {
    return i + n*j;
}

void print_efield(ostream & out, double * A, int n,
                  double xmin, double ymin, double delta)
{
    double E_x = 0.0;
    double E_y = 0.0;
    double x = xmin;
    double y = ymin;
    for (int i=0; i<n-1; i++) {
        for (int j=0; j<n-1; j++) {
            E_x = (A[index(i+1,j,n)] - A[index(i,j,n)])/delta;
            E_y = (A[index(i,j+1,n)] - A[index(i,j,n)])/delta;
            out << x + i*delta << " "
                << y + j*delta << " "
                << E_x << " " << E_y << endl;
        }
    }
}

void print_potential(ostream & out, double * A, int n,
                     double xmin, double ymin, double delta)
{
    double x = xmin;
    double y = ymin;
    for (int i=0; i<n-1; i++) {
        for (int j=0; j<n-1; j++) {
            out << x + i*delta << " "
                << y + j*delta << " "
                << A[index(i+1,j,n)] << endl;
        }
    }
}

void print_matrix(ostream & out, double * A, int n)
{
    for (int i=0; i<n; i++) {
        if (debug)
            out << "[";
        for (int j=0; j<n; j++) {
            out << setw(5) << A[index(i,j,n)] << " ";
        }
        if (debug)
            out << "]";
        out << endl;
    }
}

void print_map(ostream & out, double * A, int n,
               double xmin, double delta_x, double ymin, double delta_y)
{
    double x = xmin;
    double y = ymin;
    for (int i=0; i<n; i++) {
        for (int j=0; j<n; j++) {
            out << x + delta_x*i << " "
                << y + delta_y*j << " ";
            out << setw(5) << A[index(i,j,n)] << endl;
        }
    }
}

void copy_matrix(double * A_source, double * A_dest, int n)
{
    for (int i=0; i<n; i++) {
        for (int j=0; j<n; j++) {
            A_dest[index(i,j,n)] = A_source[index(i,j,n)];
        }
    }
}

int main(int argc, char * argv[])
{
   int opt;
   ofstream out_file;
   bool time_accuracy = false;
   bool final_matrix = false;
   bool electric_field = false;

    // Handle command line options
    while((opt = getopt(argc, argv, "do:h?ame")) != -1) {

        switch (opt) {

        case 'a':
            time_accuracy = true;
            break;

        case 'e':
            electric_field = true;
            break;

        case 'm':
            final_matrix = true;
            break;

        case 'd':
            debug = true;
            break;

        case 'o':
            out_file.open(optarg);
            break;

        case 'h':
        case '?':
        default:
            usage( basename( argv[0] ) );
            return 1;
        }
    }

    cout << "Grid Points inside the parallel plate condenser: ";
    int N;
    cin >> N;

    double a = -1.0;
    double b = +1.0;
    double dx = (b - a)/(N+1);

    if (out_file.is_open()) {
        out_file << "# Grid Division: " << N << endl;
        out_file << "# Grid Step: " << dx
                 << " from " << a << " to " << b << endl;
    }

    cout << "Delta_V accuracy to stop computation: ";
    double dv_prec;
    cin >> dv_prec;
    if (out_file.is_open()) {
        out_file << "# Delta_V accuracy to stop computation: " << dv_prec
                 << endl;
    }

    ostream & myout = (out_file.is_open()) ? out_file : cout;

    // Create a matrix for the grid
    int M = N+1;
    double * V = new double[M*M];
    double * Vnew = new double[M*M];


    // - Init the system: set the starting conditions
    for (int i=0; i<M; i++) {
        V[index(i,   0, M)] = -1.0;
        V[index(i,   N, M)] = +1.0;
        for (int j=1; j<N; j++) {
            V[index(i,  j, M)] = 0.0;
        }
    }

    // Copy the Potential to the buffer and print it on the screen
    copy_matrix(V, Vnew, M);
    if (debug)
        print_matrix(myout, V, M);

    int iter = 0;
    double dv_global = 0.0;
    int grid_points = 0; 
    do {

        iter++;
        dv_global = 0.0;
        grid_points = 0;

        // - Evaluate the new field
        for (int i=0; i<M; i++) {
            for (int j=1; j<N; j++) {

                int i_left  = i - 1;
                int i_right = i + 1;
                int j_up   = (j - 1 + M) % M;
                int j_down = (j + 1) % M;
                Vnew[index(i,j,M)] = (V[index(i_left,  j, M)] +
                                      V[index(i_right, j, M)] +
                                      V[index(i, j_up,   M)] +
                                      V[index(i, j_down, M)]) / 4.0;

                double delta = fabs(Vnew[index(i,j,M)] - V[index(i,j,M)]);
                if (delta > 0.0) {
                    grid_points ++;
                    dv_global += delta;
                }
            }
        }

        if (debug) {
            cout << endl;
            cout << "Grid points: " << grid_points << endl;
            cout << "DeltaV_global: " << dv_global / grid_points << endl;

        }
        if (time_accuracy) {
            myout << iter << " " << dv_global << " "
                  << grid_points << " " << dv_global / grid_points << " "
                  << dv_prec << endl;
        }

        if (debug)
            print_matrix(myout, Vnew, M);

        copy_matrix(Vnew, V, M);

    } while (dv_global / grid_points > dv_prec);

    if (electric_field) {
        myout << "# Final Electric field\n";
        print_efield(myout, Vnew, M, a, a, dx);
    } else if (final_matrix) {
        myout << "# Final Potential\n";
        print_potential(myout, Vnew, M, a, a, dx);
    }

    if (out_file.is_open()) {
        out_file.close();
    }

    return 0;
}