#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <algorithm>
#include "ChemPlugin.h"

int exit_client(int status)
{
   std::cin.get();
   return status;
}

void write_line(std::ofstream& f, ChemPlugin *cp, int nx,
                double gap, double& then)
{
   double now = cp[0].Report1("Time", "years");
   if ((then - now) < gap / 1e4) {
      f << now;
      for (int i=0; i<nx; i++)
         f << "\t" << cp[i].Report1("concentration Na+", "mmol/kg");
      f << std::endl;
      then += gap;
   }
}

int main(int argc, char** argv) {
   std::cout << "Model diffusion in one dimension"
             << std::endl << std::endl;

   // Simulation parameters.
   int nx = 100; // number of instances along x
   double length = 100; // cm
   double deltax = length / nx; // cm
   double deltay = 1.0, deltaz = 1.0; // cm
   double porosity = 0.25;  // volume fraction

   double diffcoef = 1e-6;  // cm2/s
   double trans = deltay * deltaz * porosity * diffcoef / deltax;  // cm3/s

   double xstable = 0.9;

   double time_end = 15.0;  // years
   double delta_years = time_end / 3; // years
   double next_output = 0.0;  // years

   // Open output file and write instance positions on the first line.
   std::ofstream f;
   f.open("Diffusion.txt");
   if (f.is_open()) {
      f << "years";
      for (int i=0; i<nx; i++)
         f << "\t" << (i+0.5) * deltax;
      f << std::endl;
   }
   else {
      std::cout << "Failed to open output file" << std::endl;
      return exit_client(-1);
   }
   
   // Configure and initialize the instances.
   std::string cmd0 = "volume = " + std::to_string(deltax * deltay * deltaz) +
                      " cm3; porosity = " + std::to_string(porosity) +
                      "; time end = " + std::to_string(time_end) +
                      " years; Xstable = " + std::to_string(xstable);
   std::string cmd1 = "Na+ = 1 mmol/kg; Cl- = 1 mmol/kg";
   std::string cmd2 = "Na+ = 0.001 mmol/kg; Cl- = 0.001 mmol/kg";

   ChemPlugin *cp = new ChemPlugin[nx];
   cp[0].Console("stdout");
   cp[0].Config("pluses = banner");

   for (int i=0; i<nx; i++) {
      cp[i].Config(cmd0);
      if (i < nx/2)
         cp[i].Config(cmd1);
      else
         cp[i].Config(cmd2);

      if (cp[i].Initialize()) {
         std::cout << "Instance " << i << " failed to initialize" << std::endl;
         return exit_client(-1);
      }
   }
   write_line(f, cp, nx, delta_years, next_output);

   // Link the instances.
   for (int i=1; i<nx; i++) {
      CpiLink link = cp[i].Link(cp[i-1]);
      link.Transmissivity(trans, "cm3/s");
   }

   // Time marching loop.
   while (true) {
      double deltat = 1e99; 
      for (int i=0; i<nx; i++)
         deltat = std::min(deltat, cp[i].ReportTimeStep());
      for (int i=0; i<nx; i++)
         if (cp[i].AdvanceTimeStep(deltat)) return exit_client(0);
      for (int i=0; i<nx; i++)
         if (cp[i].AdvanceTransport()) return exit_client(-1);
      for (int i=0; i<nx; i++)
         if (cp[i].AdvanceChemical()) return exit_client(-1);

      write_line(f, cp, nx, delta_years, next_output);
   }

   // Never gets here.
   return 0;
}