# 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

int OPTION=0, ODDEVEN_X, ODDEVEN_Y;

struct point_c {
  double x1;
  double x2;
};
typedef struct point_c POINT;
struct w_a_data {
  double x1;
  double x2;
  double w;
  double a;
};
typedef struct w_a_data SDATA;

double wfunc(double *, double *, POINT *, double *, double *, double *);
double wseries(double *, double *, double *, double *, int *, double *, int *, double *, double *, double *);
double aseries(double *, double *, double *, int *, double *, int *, double *, double *);
double afunc(double *, double *, POINT *, double *);
void rotate(double *, POINT *);
double distance(POINT *, POINT *);
double distance2(POINT *, POINT *);
void setpoint(POINT *, double, double, double *);
  
int main(int argc, char **argv)
{
  /* Wave length , Phase (NUTEE=1 -> Angle=2PI) */
  double LAMBDA=1.0;
  double NUTEE;                       /* phase of time */
  /* Radius of Scatterer*/
  double SR=0.30;
  /* Matrix of Scatterers */
  int nax, nay;                       /* number of scatterers */
  double sax, say;                    /* space in horizontal and vertical */
  double phi;                         /* Rotation angle */
  /* Output */
  POINT rangex, rangey;               /* Range of X and Y */
  int tickx, ticky;                   /* Number of ticks X and Y */
  double spacex, spacey;
  SDATA *waves;
#ifdef _OPENMP
  int numproc;
#endif
  int n, m;
  FILE *OUTPUTFILE_W, *OUTPUTFILE_A;
  char line[256], buf1[256], buf2[256];
  
  --argc;
  ++argv;

  if((argc!=12)&&(argc!=13)) {
    printf("Number of arguments is wrong! Check input parameters.\n");
    exit(1);
      }
    sscanf(*argv, "%lf", &NUTEE); argv++;
    sscanf(*argv, "%d", &nax); argv++;
    sscanf(*argv, "%lf", &sax); argv++;
    sscanf(*argv, "%d", &nay); argv++;
    sscanf(*argv, "%lf", &say); argv++;
    sscanf(*argv, "%lf", &phi); argv++;
    sscanf(*argv, "%lf", &rangex.x1); argv++;
    sscanf(*argv, "%lf", &rangex.x2); argv++;
    sscanf(*argv, "%lf", &rangey.x1); argv++;
    sscanf(*argv, "%lf", &rangey.x2); argv++;
    sscanf(*argv, "%d", &tickx); argv++;
    sscanf(*argv, "%d", &ticky);
    if(argc==13){OPTION=1; argv++;
      sscanf(*argv, "%s", line);
    }

    /* Conversion to radian */
  phi *= TORAD;

  ODDEVEN_X = nax % 2;
  ODDEVEN_Y = nay % 2;
  if(ODDEVEN_X) {nax = (nax / 2) + 1;} else {nax = (nax / 2);};
  if(ODDEVEN_Y) {nay = (nay / 2) + 1;} else {nay = (nay / 2);};
  spacex = (rangex.x2 - rangex.x1)/tickx;
  spacey = (rangey.x2 - rangey.x1)/ticky;

  waves = (SDATA *)malloc(sizeof(SDATA) * (tickx + 1) * (ticky + 1));
  if(waves == NULL) {printf("Memory allocation error!\n");exit(1);};


  if(OPTION==1) {
    strcpy(buf1, line);
    strcat(buf1, "_SctrWave.dat");
    OUTPUTFILE_W = fopen(buf1,"w");

    strcpy(buf2, line);
    strcat(buf2, "_Atom.dat");
    OUTPUTFILE_A = fopen(buf2,"w");
  } else {
    OUTPUTFILE_W = fopen("out_SctrWave.dat","w");
    OUTPUTFILE_A = fopen("out_Atom.dat","w");
  }

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

  #ifdef _OPENMP
  #pragma omp parallel for schedule(static)
  #endif
  for(m=0; m <= ticky; m++){
  #ifdef _OPENMP
  #pragma omp parallel for schedule(static)
  #endif
      for(n=0; n <= tickx; n++){
        waves[m*tickx + n].x1 = spacex*n + rangex.x1;
        waves[m*tickx + n].x2 = spacey*m + rangey.x1;
        waves[m*tickx + n].w = wseries(&waves[m*tickx + n].x1, &waves[m*tickx + n].x2, &NUTEE, &LAMBDA, &nax, &sax, &nay, &say, &phi, &SR);
        waves[m*tickx + n].a = aseries(&waves[m*tickx + n].x1, &waves[m*tickx + n].x2,  &SR, &nax, &sax, &nay, &say, &phi);
      }
  }

  for(m=0; m <= ticky; m++){
      for(n=0; n <= tickx; n++){
	fprintf(OUTPUTFILE_W, "%9.5f %9.5f %9.5f\n", waves[m*tickx + n].x1, waves[m*tickx + n].x2, waves[m*tickx + n].w);
	fprintf(OUTPUTFILE_A, "%9.5f %9.5f %9.5f\n", waves[m*tickx + n].x1, waves[m*tickx + n].x2, waves[m*tickx + n].a);
      }
  }

  fclose(OUTPUTFILE_W);
  fclose(OUTPUTFILE_A);
}

double wfunc(double *x_ptr, double *y_ptr, POINT *P_ptr, double *NUTEE_ptr, double *LAMBDA_ptr, double *SR_ptr) {
  POINT Q;
  double value;
  double d;

  Q.x1 = (*x_ptr);
  Q.x2 = (*y_ptr);
  d = distance(&Q, P_ptr);

  value = ( d < *SR_ptr ? 0.0 : (cos( 2 * PI * ( (*NUTEE_ptr) - ( P_ptr->x1 + d)/(*LAMBDA_ptr) - 1.0/2.0) )/d) );
  return(value);
}

double wseries(double *x_ptr, double *y_ptr, double *NUTEE_ptr, double *LAMBDA_ptr, int *nax_ptr, double *sax_ptr, int *nay_ptr, double *say_ptr, double *phi_ptr, double *SR_ptr) {
  double waveseries=0.0, coordx, coordy, coefx, coefy, counterx_1, countery_1;
  POINT P, Q, R, S;
  int counterx, countery;

  for(counterx=0; counterx < (*nax_ptr); counterx++){
    for(countery=0; countery < (*nay_ptr); countery++){
      if(ODDEVEN_X) {counterx_1 = counterx + 0.0;} else {counterx_1 = counterx + 0.5;};
      if(ODDEVEN_Y) {countery_1 = countery + 0.0;} else {countery_1 = countery + 0.5;};
      coordx = counterx_1 * (*sax_ptr);
      coordy = countery_1 * (*say_ptr);

      setpoint(&P, coordx, coordy, phi_ptr);
      setpoint(&Q, -coordx, -coordy, phi_ptr);
      setpoint(&R, -coordx, coordy, phi_ptr);
      setpoint(&S, coordx, -coordy, phi_ptr);

      if((ODDEVEN_X)&&(counterx_1==0)) {coefx = 0.5;} else {coefx = 1.0;}
      if((ODDEVEN_Y)&&(countery_1==0)) {coefy = 0.5;} else {coefy = 1.0;}
	waveseries = waveseries
	+ coefx * coefy * (wfunc(x_ptr, y_ptr, &P, NUTEE_ptr, LAMBDA_ptr, SR_ptr) + wfunc(x_ptr, y_ptr, &Q, NUTEE_ptr, LAMBDA_ptr, SR_ptr)
	   + wfunc(x_ptr, y_ptr, &R, NUTEE_ptr, LAMBDA_ptr, SR_ptr) + wfunc(x_ptr, y_ptr, &S, NUTEE_ptr, LAMBDA_ptr, SR_ptr));
      }
  }
  waveseries += 1.0;
  return(waveseries);
}


double aseries(double *x_ptr, double *y_ptr, double *SR_ptr, int *nax_ptr, double *sax_ptr, int *nay_ptr, double *say_ptr, double *phi_ptr) {
  double atomseries=0.0, coordx, coordy, coefx, coefy, counterx_1, countery_1;
  POINT P, Q, R, S;
  int counterx, countery;

  for(counterx=0; counterx < (*nax_ptr); counterx++){
    for(countery=0; countery < (*nay_ptr); countery++){
      if(ODDEVEN_X) {counterx_1 = counterx + 0.0;} else {counterx_1 = counterx + 0.5;};
      if(ODDEVEN_Y) {countery_1 = countery + 0.0;} else {countery_1 = countery + 0.5;};
      coordx = counterx_1 * (*sax_ptr);
      coordy = countery_1 * (*say_ptr);

      setpoint(&P, coordx, coordy, phi_ptr);
      setpoint(&Q, -coordx, -coordy, phi_ptr);
      setpoint(&R, -coordx, coordy, phi_ptr);
      setpoint(&S, coordx, -coordy, phi_ptr);

      if((ODDEVEN_X)&&(counterx_1==0)) {coefx = 0.5;} else {coefx = 1.0;}
      if((ODDEVEN_Y)&&(countery_1==0)) {coefy = 0.5;} else {coefy = 1.0;}
      atomseries = atomseries
	+ coefx * coefy * (afunc(x_ptr, y_ptr, &P, SR_ptr) + afunc(x_ptr, y_ptr, &Q, SR_ptr)
	   + afunc(x_ptr, y_ptr, &R, SR_ptr) + afunc(x_ptr, y_ptr, &S, SR_ptr) );
      }
  }
  atomseries += -1.1;
  return(atomseries);
}

double afunc(double *x_ptr, double *y_ptr, POINT *P_ptr, double *SR_ptr) {
  double atomf=0.0;
  POINT Q;

  Q.x1 = (*x_ptr);
  Q.x2 = (*y_ptr);
  atomf = ( distance(&Q, P_ptr)<(*SR_ptr) ? (2.6 - 2*distance2(&Q, P_ptr) ) : 0.0 );
  return(atomf);
}

void rotate(double *phi_ptr, POINT *point_ptr){
  double rx, ry;
  rx = cos(*phi_ptr) * (point_ptr->x1) - sin(*phi_ptr) * (point_ptr->x2);
  ry = sin(*phi_ptr) * (point_ptr->x1) + cos(*phi_ptr) * (point_ptr->x2);
  point_ptr->x1 = rx;
  point_ptr->x2 = ry;
}

double distance2(POINT *S_ptr, POINT *T_ptr)
{
  double dx, dy, r;
  dx = T_ptr->x1 - S_ptr->x1;
  dy = T_ptr->x2 - S_ptr->x2;
  r = dx * dx + dy * dy;
  return(r);
}

double distance(POINT *S_ptr, POINT *T_ptr)
{
  double r;
  r = sqrt( distance2(S_ptr, T_ptr) );
  return(r);
}

void setpoint(POINT *P_ptr, double coordx, double coordy, double *phi_ptr)
{
  P_ptr->x1 = coordx;
  P_ptr->x2 = coordy;
  rotate(phi_ptr, P_ptr);
}