#ifndef SIMULATION2D_HPP
#define SIMULATION2D_HPP
#include "Rod2d.hpp"
#include "LiveAgg.hpp"
#include <fstream>
#include <iostream>
/* This is used to calculate the convergence order of the particles
*/

inline void printProgress(double progress)
{
  const int barWidth = 70;
  std::cout << "[";
  const int pos = static_cast<int>(barWidth * progress);
  for (int i = 0; i < barWidth; ++i) {
    if (i < pos) {
      std::cout << "=";
    } else if (i == pos) {
      std::cout << ">";
    } else {
      std::cout << " ";
    }
  }
  std::cout << "] " << int(progress * 100.0) << " %   \r";
  std::cout.flush();
}

template<void (*T)(Rod2d &)>
void simulate(double L, int seed, const std::string &filename)
{
  // we now simulate the msd after 10 steps


  Rod2d rod(L, 1e-5, seed);
  const int sampleWidth = 1000;
  std::vector<LiveAgg> msd(sampleWidth, LiveAgg());
  std::vector<LiveAgg> oaf(sampleWidth, LiveAgg());
  std::vector<LiveAgg> empxx(sampleWidth, LiveAgg());
  std::vector<LiveAgg> empyy(sampleWidth, LiveAgg());
  std::vector<double> xPos(sampleWidth, 0.0);
  std::vector<double> yPos(sampleWidth, 0.0);
  std::vector<double> e1(sampleWidth, 0.0);
  std::vector<double> e2(sampleWidth, 0.0);
  std::vector<double> targetMSD(sampleWidth, 0.0);
  std::vector<double> targetOAF(sampleWidth, 0.0);
  std::vector<double> targetEMPXX(sampleWidth, 0.0);
  std::vector<double> targetEMPYY(sampleWidth, 0.0);
  std::vector<size_t> samplePoints(sampleWidth, 1);
  size_t temp = 1;
  const double Dmean = 0.5 * (rod.Dpara + rod.Dortho);
  const double deltaD = rod.Dortho - rod.Dpara;
  for (size_t i = 0; i < samplePoints.size(); ++i) {
    samplePoints[i] = temp;
    temp += 100;
    const double time = static_cast<double>(temp) * rod.deltaT;
    targetMSD[i] = 4.0 * 0.5 * (rod.Dpara + rod.Dortho) * time;
    targetOAF[i] = std::exp(-rod.Drot * time);
    const double u = 4.0;
    targetEMPXX[i] = Dmean - deltaD / 2.0 * (1 - exp(-u * rod.Drot * time)) / u / rod.Drot / time;
    targetEMPYY[i] = Dmean + deltaD / 2.0 * (1 - exp(-u * rod.Drot * time)) / u / rod.Drot / time;
  }
  unsigned long long int iteration = 0;
  //const unsigned long long int maxItr = ((long int)10000)*1000000;
  const unsigned long long int maxItr = 10000000;
  const unsigned long long int updateInter = 10000;
  while (iteration < maxItr) {
    T(rod);
    iteration++;
    for (size_t i = 0; i < samplePoints.size(); ++i) {
      const size_t sample = samplePoints[i];
      if (iteration % sample == 0) {
        msd[i].feed((rod.x1 - xPos[i]) * (rod.x1 - xPos[i]) + (rod.x2 - yPos[i]) * (rod.x2 - yPos[i]));
        oaf[i].feed(rod.e1 * e1[i] + rod.e2 * e2[i]);
        empxx[i].feed(std::pow((rod.x1 - xPos[i]) * e1[i] + (rod.x2 - yPos[i]) * e2[i], 2));
        empyy[i].feed(std::pow((rod.x1 - xPos[i]) * e2[i] - (rod.x2 - yPos[i]) * e1[i], 2));
        xPos[i] = rod.x1;
        yPos[i] = rod.x2;
        e1[i] = rod.e1;
        e2[i] = rod.e2;
      }
    }
    if (iteration % updateInter == 0) {
      printProgress(static_cast<double>(iteration)/ maxItr);
    }
  }

  std::cout << "\n";

  std::ofstream fileMSD(filename + "_" + std::to_string(L) + "_MSD.out");
  std::ofstream fileOAF(filename + "_" + std::to_string(L) + "_OAF.out");
  std::ofstream fileXX(filename + "_" + std::to_string(L) + "_XX.out");
  std::ofstream fileYY(filename + "_" + std::to_string(L) + "_YY.out");

  for (size_t i = 0; i < samplePoints.size(); ++i) {
    fileMSD << static_cast<double>(samplePoints[i]) * rod.deltaT << "\t" << targetMSD[i] << "\t" << msd[i].getMean() << "\t" << msd[i].getSEM() << "\t" << msd[i].getSD() << "\t" << msd[i].getNumPoints() << "\n";
    fileOAF << static_cast<double>(samplePoints[i]) * rod.deltaT << "\t" << targetOAF[i] << "\t" << oaf[i].getMean() << "\t" << oaf[i].getSEM() << "\t" << oaf[i].getSD() << "\t" << oaf[i].getNumPoints() << "\n";
    fileXX << static_cast<double>(samplePoints[i]) * rod.deltaT << "\t" << targetEMPXX[i] << "\t" << empxx[i].getMean() << "\t" << empxx[i].getSEM() << "\t" << empxx[i].getSD() << "\t" << empxx[i].getNumPoints() << "\n";
    fileYY << static_cast<double>(samplePoints[i]) * rod.deltaT << "\t" << targetEMPYY[i] << "\t" << empyy[i].getMean() << "\t" << empyy[i].getSEM() << "\t" << empyy[i].getSD() << "\t" << empyy[i].getNumPoints() << "\n";
  }
  fileMSD.close();
  fileOAF.close();
  fileXX.close();
  fileYY.close();
}


template<void (*T)(Rod2d &)>
void konvergenz(double L, int seed, double deltaT, const std::string &filename)
{
  // we now simulate the msd after 10 steps
  Rod2d rod(L, deltaT, seed);
  const double time = 2;
  const int steps = static_cast<int>(time / deltaT);
  const int minSteps = 1000;

  //const double targetMSD = 4.0*0.5*(rod.Dpara+rod.Dortho)*time;
  const double targetOAF = std::exp(-rod.Drot * time);
  /*const double u = 4.0;
    const double Dmean = 0.5*(rod.Dpara+rod.Dortho);
    const double deltaD = rod.Dortho - rod.Dpara;
    const double targetEMPXX = Dmean-deltaD/2.0*(1 - exp(-u*rod.Drot*time))/u/rod.Drot / time;
    const double targetEMPYY = Dmean+deltaD/2.0*(1 - exp(-u*rod.Drot*time))/u/rod.Drot / time;*/
  //LiveAgg msd;
  LiveAgg oaf;
  //LiveAgg empxx;
  //LiveAgg empyy;
  int meanCount = 0;
  while (meanCount < minSteps or std::abs(targetOAF - oaf.getMean()) < oaf.getSEM() * 10) {
    //reset the position
    rod.reset();
    for (int i = 0; i < steps; ++i) {
      T(rod);
    }
    //msd.feed(rod.x1*rod.x1+rod.x2*rod.x2);
    oaf.feed(rod.e1);
    //.feed(rod.x1*rod.x1/2/time);
    //empyy.feed(rod.x2*rod.x2/2/time);
    if (meanCount % 100000 == 0) {
      //std::cout<<msd.getSEM()<<" "<<std::abs(targetMSD-msd.getMean());
      std::cout << "\t" << oaf.getSEM() << " " << std::abs(targetOAF - oaf.getMean()) << std::endl;
      //std::cout<<"\t"<<empxx.getSEM()<<" "<<std::abs(targetEMPXX-empxx.getMean());
      //std::cout<<"\t"<<empyy.getSEM()<<" "<<std::abs(targetEMPYY-empyy.getMean())<<std::endl;;
    }
    if (meanCount % 1000000 == 0) {
      std::ofstream fileOAF(filename + "_length=" + std::to_string(L) + "_dt=" + std::to_string(rod.deltaT) + "_OAF.out");
      fileOAF << rod.deltaT << "\t" << std::abs(targetOAF - oaf.getMean()) << "\t" << oaf.getSEM() << "\n";
      fileOAF.close();
    }
    meanCount++;
  }
  std::ofstream fileOAF(filename + "_length=" + std::to_string(L) + "_dt=" + std::to_string(rod.deltaT) + "_OAF.out");
  fileOAF << rod.deltaT << "\t" << std::abs(targetOAF - oaf.getMean()) << "\t" << oaf.getSEM() << "\n";
  fileOAF.close();
}


#endif// SIMULATION2D_HPP