#include <iostream>
#include <fstream>
#include <complex>
#include <cmath>
#ifdef USE_IMAGINARY
#include "include/Imaginary.hpp"
#endif
#include "boost/array.hpp"
#include "boost/timer.hpp"

/**
 * \file ImaginaryExample.cpp
 * \brief Example code to present the strength of the Imaginary<T> class
 *
 * The timed section of code is ~2.3 times faster with g++ 4.6.1 on
 * an Intel core i7-2820QM when compiled in -O2 or -O3 mode and with -DUSE_IMAGINARY
 * than without without this last compiling option
 * \version 0.99
 * \author Matthieu Schaller
 * \date 13.11.2011
 */

//Potential barrier height
const double V0 = 1.;

//Dirack's constant
const double hbar = 1.;

//Particle mass
const double m = 0.5;

//Domain
const double xl = -500.;
const double xr = 1000.;

//Number of points on the interval
const size_t nbPoints = 2000;

//Time step size
const double deltaT = 0.005;

//Imaginary unit
#ifdef USE_IMAGINARY
const Imaginary<double> ii(1.);
#else
const std::complex<double> ii(0.,1.);
using std::exp;
#endif

//Potential
double V(double x)
{
    return fabs(x-600.) < 100. ? 1. : 0.;
}

int main()
{
    //Usefull constants
    const double deltaX = (xr-xl)/static_cast<double>(nbPoints);
    const double x0 = xl + (xr-xl)*0.35;
    const double sigma0 = (xr-xl)*0.0375;
    const int n = 32;
    const double zk = 2.*M_PI*n/(xr-xl);

    //Spatial Domain
    boost::array<double, nbPoints+4> x;
    for(size_t i(0); i<nbPoints+4; ++i)
        x[i] = -500. + (static_cast<double>(i)-2.)*deltaX;

    //Wavefunction at t-\Delta t, t and t+\Delta t
    boost::array<std::complex<double>, nbPoints+4> Psi_old, Psi_cur, Psi_new;

    //Build an initial wavepacket
    for(size_t i(0); i<nbPoints+4; ++i)
    {
        Psi_old[i] = exp(ii*zk*x[i]) * exp(-(x[i]-x0)*(x[i]-x0)/(2.*sigma0*sigma0));
        Psi_cur[i] = Psi_old[i];
    }

    boost::timer timer;
    double computingTime(0.);

    for(double t(0.); t<2000.; t+=deltaT)
    {
        //Evolve the wavefunction
        timer.restart();
        for(size_t i(2); i<nbPoints+2; ++i)
        {
            Psi_new[i] = Psi_old[i]
                         + ii * hbar * deltaT / (12.*m*deltaX*deltaX) * (-Psi_cur[i-2] + 16.*Psi_cur[i-1] - 30.*Psi_cur[i] + 16.*Psi_cur[i+1] - Psi_cur[i+2])
                         - 2. * ii * V(x[i]) / hbar * deltaT * Psi_cur[i];

        }
        computingTime += timer.elapsed();
        //Periodic boundary condition
        Psi_new[0] = Psi_new[nbPoints];
        Psi_new[1] = Psi_new[nbPoints+1];
        Psi_new[nbPoints+2] = Psi_new[2];
        Psi_new[nbPoints+3] = Psi_new[3];

        //Swap fields
        Psi_cur.swap(Psi_old);
        Psi_new.swap(Psi_cur);
    }

    //Wrtie output to a file
    std::ofstream out("output.dat");
    for(size_t i(2); i<nbPoints+2; ++i)
    {
        out << x[i];
        out << ' ' << real(Psi_cur[i]);
        out << ' ' << imag(Psi_cur[i]);
        out << ' ' << norm(Psi_cur[i]) << '\n';
    }

#ifdef USE_IMAGINARY
    std::cout << "Using Imaginary<double>(1.) as sqrt(-1)" << std::endl;
#else
    std::cout << "Using std::complex<double>(0.,1.) as sqrt(-1)" << std::endl;
#endif
    std::cout << "Time spent evolving the wavefunctions: " << computingTime << std::endl;

    return 0;
}