#include <iostream>
#include <fstream>
#include <iomanip>
#include "channelflow/flowfield.h"
#include "channelflow/dns.h"
#include "channelflow/utilfuncs.h"

using namespace std;
using namespace channelflow;

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

  string purpose("interpolate a given flowfield onto a different grid");
  ArgList args(argc, argv, purpose);  // Assign parameters

  const bool dealias  = args.getflag("-d", "--dealias", "set padding modes to zero");
  const bool fixdiv   = args.getflag("-dv", "--divergence", "fix divergence, if 3d and Ny gets reduced");
  const int Nxarg     = args.getint("-Nx", "--Nx", 0, "new # x gridpoints");
  const int Nyarg     = args.getint("-Ny", "--Ny", 0, "new # y gridpoints");
  const int Nzarg     = args.getint("-Nz", "--Nz", 0, "new # z gridpoints");
  const Real alpha    = args.getreal("-a", "--alpha", 0.0, "new alpha == 2pi/Lx");
  const Real gamma    = args.getreal("-g", "--gamma", 0.0, "new gamma == 2pi/Lz");
  const Real lx       = (alpha == 0) ? args.getreal("-lx", "--lx", 0.0,   "new Lx = 2 pi lx") : 1.0/alpha;
  const Real lz       = (gamma == 0) ? args.getreal("-lz", "--lz", 0.0,   "new Lz = 2 pi lz") : 1.0/gamma;

  const string u0name = args.getstr(2, "<infield>",  "input flow field");
  const string u1name = args.getstr(1, "<outfield>", "output flow field");

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

  FlowField u0(u0name);

  const int Nx = (Nxarg == 0 ? u0.Nx() : Nxarg);
  const int Ny = (Nyarg == 0 ? u0.Ny() : Nyarg);
  const int Nz = (Nzarg == 0 ? u0.Nz() : Nzarg);

  FlowField u1(Nx,Ny,Nz,u0.Nd(),u0.Lx(),u0.Lz(),u0.a(),u0.b());
  u1.interpolate(u0);

  if (dealias || (Nx >= (3*u0.Nx())/2 && Nx >= (3*u0.Nz())/2))
    u1.zeroAliasedModes();

  Real Lx = (lx == 0) ? u0.Lx() : 2*pi*lx;
  Real Lz = (lz == 0) ? u0.Lz() : 2*pi*lz;

  if (Lx != u0.Lx() || Lz != u0.Lz())
    u1.rescale(Lx,Lz);

  cout << setprecision(16);
  cout << "L2Norm(u0)  == " << L2Norm(u0) << endl;
  cout << "L2Norm(u1)  == " << L2Norm(u1) << endl;
  cout << "bcNorm(u0)  == " << bcNorm(u0) << endl;
  cout << "bcNorm(u1)  == " << bcNorm(u1) << endl;

  if (u0.Nd() == 3) {
    cout << "divNorm(u0) == " << divNorm(u0) << endl;
    cout << "divNorm(u1) == " << divNorm(u1) << endl;

    if (Ny < u0.Ny() && fixdiv) {
      cout <<"fixing diveregence..." << endl;
      Vector v;
      field2vector(u1,v);
      vector2field(v,u1);
      cout << "divNorm(u1) == " << divNorm(u1) << endl;
      cout << "L2Norm(u1)  == " << L2Norm(u1) << endl;
      cout << "bcNorm(u1)  == " << bcNorm(u1) << endl;
    }
  }
  u1.binarySave(u1name);
}