/***************************************************************/
/***************************************************************/
/***************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <complex>

#include "ctl-math.h"
#include "ctlgeom.h"
#include "meep.hpp"
#include "meep/meepgeom.hpp"
#include "meep/vec.hpp"

using namespace meep;
using namespace std;

typedef std::complex<double> cdouble;

vector3 v3(double x, double y = 0.0, double z = 0.0) {
  vector3 v;
  v.x = x;
  v.y = y;
  v.z = z;
  return v;
}

/***************************************************************/
/* dummy material function needed to pass to structure( )      */
/* constructor as a placeholder before we can call             */
/* set_materials_from_geometry                                 */
/***************************************************************/
double dummy_eps(const vec &) { return 1.0; }

// Pretty=true  --> {x,y,z}
// Pretty=false -->  x y z
void fprint_vec(FILE *f, vec v, bool Pretty = false) {
  const char *s = Pretty ? "{" : " ";
  LOOP_OVER_DIRECTIONS(v.dim, d) {
    fprintf(f, "%s%e", s, v.in_direction(d));
    if (Pretty) s = ",";
  }
  fprintf(f, "%s", Pretty ? "} " : " ");
}

// Pretty=true  --> {x,y,z}
// Pretty=false -->  x y z
void fprint_ivec(FILE *f, ivec v, bool Pretty = false) {
  const char *s = Pretty ? "{" : " ";
  LOOP_OVER_DIRECTIONS(v.dim, d) {
    fprintf(f, "%s%i", s, v.in_direction(d));
    if (Pretty) s = ",";
  }
  fprintf(f, "%s", Pretty ? "} " : " ");
}

/***************************************************************/
/***************************************************************/
/***************************************************************/
static void print_chunk_info(fields_chunk *fc, int ichunk, component cgrid, ivec is, ivec ie,
                             vec s0, vec s1, vec e0, vec e1, double dV0, double dV1, ivec shift,
                             complex<double> shift_phase, const symmetry &S, int sn, void *data_) {
  (void)s0;
  (void)s1;
  (void)e0;
  (void)e1;
  (void)dV0;
  (void)dV1;
  (void)shift_phase;
  (void)data_;

  int iproc = fc->n_proc();
  char FileName[100];
  snprintf(FileName, 100, "/tmp/ChunkInfo_%i_%i", iproc, ichunk);
  FILE *f = fopen(FileName, "w");

  // write some info on the chunk
  fprintf(f, "# chunk %i (process %i) \n", ichunk, iproc);
  fprintf(f, "# is={%i,%i,%i}\n", is.in_direction(X), is.in_direction(Y), is.in_direction(Z));
  fprintf(f, "# ie={%i,%i,%i}\n", ie.in_direction(X), ie.in_direction(Y), ie.in_direction(Z));
  fprintf(f, "# Ex=%s\n", component_name(S.transform(Ex, -sn)));
  fprintf(f, "# Ey=%s\n", component_name(S.transform(Ey, -sn)));
  fprintf(f, "# Hx=%s\n", component_name(S.transform(Hx, -sn)));
  fprintf(f, "# Hy=%s\n", component_name(S.transform(Hy, -sn)));
  fprintf(f, "# \n");
  fprintf(f, "# columns below: \n");
  fprintf(f, "# 1,2   process, chunk index\n");
  fprintf(f, "# 3-5   grid-point indices     (physical, i.e. before symmetry)\n");
  fprintf(f, "# 6-8   grid-point coordinates (physical, i.e. before symmetry)\n");
  fprintf(f, "# 9-11  grid-point indices     (logical,  i.e. after symmetry)\n");
  fprintf(f, "# 12-14 grid-point coordinates (logical,  i.e. after symmetry)\n");

  // loop over all grid points
  vec rshift(shift * (0.5 * fc->gv.inva));
  LOOP_OVER_IVECS(fc->gv, is, ie, idx) {
    fprintf(f, "%i %i ", iproc, ichunk);

    IVEC_LOOP_ILOC(fc->gv, iloc);
    fprint_ivec(f, iloc);
    iloc = S.transform(iloc, sn) + shift;
    fprint_ivec(f, iloc);

    IVEC_LOOP_LOC(fc->gv, loc);
    fprint_vec(f, loc);
    loc = S.transform(loc, sn) + rshift;
    fprint_vec(f, loc);

    fprintf(f, "\n");
  }
  fclose(f);
}

/***************************************************************/
/***************************************************************/
/***************************************************************/
structure create_structure(double resolution, int symmetries) {
  double sx = 8.0;   // size of cell in X direction
  double sy = 6.0;   // size of cell in Y direction
  double w = 3.0;    // width of waveguide
  double dpml = 1.0; // PML thickness

  geometry_lattice.size.x = sx;
  geometry_lattice.size.y = sy;
  geometry_lattice.size.z = 0.0;
  grid_volume gv = voltwo(sx, sy, resolution);
  gv.center_origin();
  symmetry S = (symmetries == 1 ? mirror(Y, gv)
                                : symmetries == 2 ? mirror(X, gv) + mirror(Y, gv) : symmetry());
  structure the_structure(gv, dummy_eps, pml(dpml), S);

  vector3 e1 = v3(1.0, 0.0, 0.0);
  vector3 e2 = v3(0.0, 1.0, 0.0);
  vector3 e3 = v3(0.0, 0.0, 1.0);
  meep_geom::material_type dielectric = meep_geom::make_dielectric(12.0);

  vector3 center = v3(0.0, 0.0, 0.0);
  vector3 size = v3(sx, w, 0.0);
  geometric_object objects[1];
  objects[0] = make_block(dielectric, center, e1, e2, e3, size);
  geometric_object_list g = {1, objects};
  meep_geom::set_materials_from_geometry(&the_structure, g);

  return the_structure;
}

/***************************************************************/
/***************************************************************/
/***************************************************************/
int main(int argc, char *argv[]) {
  initialize mpi(argc, argv);

  /***************************************************************/
  /***************************************************************/
  /***************************************************************/
  double resolution = 5.0;
  int symmetries = 0; // {0,1,2} for none, Y-mirror, XY-mirror
  for (int narg = 1; narg < argc; narg++)
    if (!strcasecmp(argv[narg], "--resolution")) sscanf(argv[narg + 1], "%le", &resolution);
  for (int narg = 1; narg < argc - 1; narg++)
    if (!strcasecmp(argv[narg], "--symmetries")) sscanf(argv[narg + 1], "%i", &symmetries);

  /***************************************************************/
  /***************************************************************/
  /***************************************************************/
  structure s = create_structure(resolution, symmetries);
  fields f(&s);
  f.step();
  const char *prefix = "WriteChunkInfo";
  // f.output_hdf5(Dielectric, f.total_volume(), 0, false, true, prefix);

  f.loop_in_chunks(print_chunk_info, 0, f.total_volume());
}