#include <iostream>
#include <fstream>
#include <iomanip>
#include <stdlib.h>
#include <sys/stat.h>
#include "channelflow/flowfield.h"
#include "channelflow/diffops.h"
#include "channelflow/periodicfunc.h"
#include "channelflow/utilfuncs.h"

using namespace std;
using namespace channelflow;

int main(int argc, char* argv[]) {

  string purpose("save slices of fields for movie-making in Matlab");
  ArgList args(argc, argv, purpose);

  // Assign parameters
  const Real T0        = args.getreal("-T0", "--T0", 0, "start time");
  const Real T1        = args.getreal("-T1", "--T1", 500, "end time");
  const Real dT        = args.getreal("-dT", "--dT", 1.0, "save interval");
  const int  f0        = args.getint("-f0", "--f0", 0, "start frame number");
  const Real ax        = args.getreal("-ax", "-axshift", 0.0, "shift data by ax*Lx in x");
  const Real az        = args.getreal("-az", "-azshift", 0.0, "shift data by az*Lz in z");
  const int  xs        = args.getint ("-xs", "-xstride", 1, "stride for x output");
  const int  ys        = args.getint ("-ys", "-ystride", 1, "stride for x output");
  const int  zs        = args.getint ("-zs", "-zstride", 1, "stride for x output");
  const string idir    = args.getpath("-d",  "--datadir", "data/",  "data dir");
  const string odir    = args.getpath("-o",  "--outdir", "frames/", "output dir");
  const string label   = args.getstr ("-l", "--label", "u", "prefix of field files");

  args.check();
  args.save("./");

  const bool inttime = (abs(dT - int(dT)) < 1e-12) ? true : false;
  const bool shift = (ax != 0.0 || az != 0.0) ? true : false;
  FieldSymmetry tau(1,1,1,ax, az);
  int frame = f0;

  for (Real t=T0; t<=T1; t += dT) {
    string ts = t2s(t, inttime);
    cout << ts << ' ' << flush;

    FlowField u(idir + label + ts);
    if (shift)
      u *= tau;

    plotfield(u, odir, label + i2s(frame),xs,ys,zs);

    ofstream os((odir + label + i2s(frame++) + "_ide.asc").c_str());
    u.makeSpectral();
    u.cmplx(0,0,0,0) += Complex(1.0, 0);
    os << forcing(u) << ' ' << dissipation(u) << ' ' << 0.5*L2Norm2(u) << endl;

    if  (t==T0) {
      Vector x = periodicpoints(u.Nx(), u.Lx());
      Vector y = chebypoints(u.Ny(), u.a(), u.b());
      Vector z = periodicpoints(u.Nz(), u.Lz());
      ofstream osx((odir + "x.asc").c_str());
      ofstream osy((odir + "y.asc").c_str());
      ofstream osz((odir + "z.asc").c_str());
      for (int nx=0; nx<=u.Nx(); nx += xs)
	osx << x(nx) << '\n';
      for (int ny=0; ny<u.Ny(); ny += ys)
	osy << y(ny) << '\n';
      for (int nz=0; nz<=u.Nz(); nz += zs)
	osz << z(nz) << '\n';
    }
  }
  cout << endl;

}