#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>

#ifdef _OPENMP
#include <omp.h>
#endif

#define PI 3.14159265358979323846264338328
#define TORAD PI/180.0
#define BUF_SIZE (1<<20)

/* Structure of Arrays_*/
typedef struct
{
  double *x;
  double *y;
  double *w;  /* weight */
  int n;
} AtomsA;

typedef struct
{
  double min, max;
} Range;

typedef struct
{
  double x;
  double y;
  double w;
  double a;
} SDATA;

int build_atoms_centered (AtomsA * atoms, const int naa, const int nab,
			  const double la, const double lb, const double gamma,
			  const double phi);
static inline void calc_point (double x, double y,
			       const AtomsA * restrict atoms, const double NUT,
			       const double inv_lambda, const double two_pi,
			       const double SR2, double *wave, double *atom);
int allocate_atoms (AtomsA * a, int n);
void free_atoms (AtomsA * a);

int
main (int argc, char **argv)
{
  int OPTION = 0;
  /* Wave length 1.0 */
  const double LAMBDA = 1.0, inv_lambda = 1.0 / LAMBDA, two_pi = 2.0 * PI;
  double NUT;			/* phase of time (NUT=1 -> Angle=2PI) */
  /* Radius of Scatterer */
  const double SR = 0.30;
  const double SR2 = SR * SR;
  /* Matrix parameter of Scatterers */
  int naa, nab;			/* number of scatterers */
  double la, lb;		/* length of axis a and axis b */
  double gamma, phi;		/* Angle between a and b, Rotation angle */
  /* Output */
  Range rangex, rangey;		/* Range of X and Y */
  int tickx, ticky;		/* Number of ticks X and Y */
  SDATA *waves;
#ifdef _OPENMP
  int numproc;
#endif
  int n, m;
  FILE *OUTPUTFILE_W, *OUTPUTFILE_A;
  char line[256], buf1[512], buf2[512];
  char bufW[BUF_SIZE];
  char bufA[BUF_SIZE];
  size_t posW = 0;
  size_t posA = 0;


  --argc;
  ++argv;

  if ((argc != 13) && (argc != 14)) {
    printf ("Number of arguments is wrong! Check input parameters.\n");
    exit (1);
  }
  sscanf (*argv, "%lf", &NUT);
  argv++;
  sscanf (*argv, "%d", &naa);
  argv++;
  sscanf (*argv, "%d", &nab);
  argv++;
  sscanf (*argv, "%lf", &la);
  argv++;
  sscanf (*argv, "%lf", &lb);
  argv++;
  sscanf (*argv, "%lf", &gamma);
  argv++;
  sscanf (*argv, "%lf", &phi);
  argv++;
  sscanf (*argv, "%lf", &rangex.min);
  argv++;
  sscanf (*argv, "%lf", &rangex.max);
  argv++;
  sscanf (*argv, "%lf", &rangey.min);
  argv++;
  sscanf (*argv, "%lf", &rangey.max);
  argv++;
  sscanf (*argv, "%d", &tickx);
  argv++;
  sscanf (*argv, "%d", &ticky);
  if (argc == 14) {
    OPTION = 1;
    argv++;
    sscanf (*argv, "%s", line);
  }

  /* Conversion to Radian */
  gamma *= TORAD;
  phi *= TORAD;

  if (tickx <= 0 || ticky <= 0) {
    fprintf (stderr, "Error: tickx and ticky must be positive\n");
    exit (1);
  }
  /* Grid spacings */
  const double g_lx = (rangex.max - rangex.min) / tickx;
  const double g_ly = (rangey.max - rangey.min) / ticky;

  /* Memory allocation */
  size_t total_points = (size_t) (ticky + 1) * (tickx + 1);
  waves = (SDATA *) malloc (sizeof (SDATA) * total_points);
  if (waves == NULL) {
    fprintf (stderr, "Error: Memory allocation error!\n");
    goto cleanup;
  }

  /* Output file name */
  if (OPTION == 1) {
    snprintf (buf1, sizeof (buf1), "%s_SctrWave.dat", line);
    snprintf (buf2, sizeof (buf2), "%s_Atom.dat", line);
  } else {
    strcpy (buf1, "out_SctrWave.dat");
    strcpy (buf2, "out_Atom.dat");
  }

  if (!(OUTPUTFILE_W = fopen (buf1, "w"))
      || !(OUTPUTFILE_A = fopen (buf2, "w"))) {
    fprintf (stderr, "Error: Could not open output files.\n");
  }


  const int n_max = naa * nab;
  AtomsA atoms = { NULL, NULL, NULL, 0 };
  if (allocate_atoms (&atoms, n_max) != 0) {
    fprintf (stderr, "Error: Failed to allocate atom memory.\n");
    goto cleanup;
  }
  /* Creation of Scatterers(matrix) */
  const int n_atoms =
    build_atoms_centered (&atoms, naa, nab, la, lb, gamma, phi);


#ifdef _OPENMP
  numproc = omp_get_num_procs ();
#ifdef DEBUG
  printf ("The number of processors is %d\n", numproc);
#endif
  omp_set_num_threads (numproc);
#endif

  /* Parallelize over observation points */
#ifdef _OPENMP
#pragma omp parallel for collapse(2) schedule(static) private (m, n)
#endif
  for (m = 0; m <= ticky; m++) {
    for (n = 0; n <= tickx; n++) {
      int id = m * (tickx + 1) + n;
      double x = g_lx * n + rangex.min;
      double y = g_ly * m + rangey.min;
      waves[id].x = x;
      waves[id].y = y;
      calc_point (x, y, &atoms, NUT, inv_lambda, two_pi, SR2, &waves[id].w,
		  &waves[id].a);
    }
  }


  setvbuf (OUTPUTFILE_W, NULL, _IOFBF, 1 << 20);
  setvbuf (OUTPUTFILE_A, NULL, _IOFBF, 1 << 20);
  for (size_t id = 0; id < total_points; id++) {
    int lenW = snprintf (bufW + posW, BUF_SIZE - posW, "%9.5f %9.5f %9.5f\n",
			 waves[id].x, waves[id].y, waves[id].w);
    int lenA = snprintf (bufA + posA, BUF_SIZE - posA, "%9.5f %9.5f %9.5f\n",
			 waves[id].x, waves[id].y, waves[id].a);
    if (lenW > 0)
      posW += lenW;
    if (lenA > 0)
      posA += lenA;

    if (posW > BUF_SIZE - 256) {
      fwrite (bufW, 1, posW, OUTPUTFILE_W);
      posW = 0;
    }
    if (posA > BUF_SIZE - 256) {
      fwrite (bufA, 1, posA, OUTPUTFILE_A);
      posA = 0;
    }
  }
  if (posW > 0)
    fwrite (bufW, 1, posW, OUTPUTFILE_W);
  if (posA > 0)
    fwrite (bufA, 1, posA, OUTPUTFILE_A);

cleanup:
  if (OUTPUTFILE_W)
    fclose (OUTPUTFILE_W);
  if (OUTPUTFILE_A)
    fclose (OUTPUTFILE_A);
  if (waves)
    free (waves);
  free_atoms (&atoms);
  return 0;
}



/* メモリ確保 */
int
allocate_atoms (AtomsA * a, int n)
{
  a->n = n;
  a->x = (double *) malloc (sizeof (double) * n);
  a->y = (double *) malloc (sizeof (double) * n);
  a->w = (double *) malloc (sizeof (double) * n);
  if (!a->x || !a->y || !a->w)
    return -1;
  return 0;
}

/* メモリ解放 */
void
free_atoms (AtomsA * a)
{
  if (a) {
    free (a->x);
    free (a->y);
    free (a->w);
    a->x = a->y = a->w = NULL;
  }
}

/* Precalculation of atomic coordinates */
int
build_atoms_centered (AtomsA * atoms, const int naa, const int nab,
		      const double la, const double lb, const double gamma,
			  const double phi)
{
  int i, j;

  const int nx = naa, ny = nab;
  const double cg = cos (gamma);
  const double sg = sin (gamma);
  /* Center of matrix */
  const double cx = 0.5 * ((nx - 1) * la + (ny - 1) * lb * cg);
  const double cy = 0.5 * (ny - 1) * lb * sg;

  const double cp = cos (phi);
  const double sp = sin (phi);

  int count = 0;
  for (i = 0; i < nx; i++) {
    for (j = 0; j < ny; j++) {

      double x = i * la + j * lb * cg - cx;
      double y = j * lb * sg - cy;

      atoms->x[count] = cp * x - sp * y;
      atoms->y[count] = sp * x + cp * y;
      atoms->w[count] = 1.0;
      count++;
    }
  }
  return count;
}


/* Calc. of wave and atom*/
static inline void
calc_point (double x, double y,
	    const AtomsA * restrict atoms,
	    const double NUT, const double inv_lambda,
	    const double two_pi, const double SR2, double *wave, double *atom)
{
  double sum_wave = 1.0;
  double sum_atom = -1.1;
  const double term1 = two_pi * (NUT - 0.5);
  const double inv_L_2pi = inv_lambda * two_pi;
  const int n_atoms = atoms->n;

  for (int i = 0; i < n_atoms; i++) {
    double dx = x - atoms->x[i];
    double dy = y - atoms->y[i];
    double r2 = dx * dx + dy * dy;
    if (r2 >= SR2) {
      double r = sqrt (r2);
      double sqrtr = sqrt (r);
      double phase = term1 - (atoms->x[i] + r) * inv_L_2pi;
      sum_wave += atoms->w[i] * cos (phase) / sqrtr;
    } else {
      sum_atom += atoms->w[i] * (2.6 - 2.0 * r2);
    }
  }
  *wave = sum_wave;
  *atom = sum_atom;
}