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

/*function outputting the parameter file: all the settings*/
void 
print_parfile(int xres, int yres, int zres,
              double dx, double dy, double dz,
              int comp, int threads, char *filebase,
              int hz, int athickness, double wavelength, double period)
{
    int verbose = 1;
    int bnd = 10;
    int xw = 2;
    int yw = 2;
    int zw = 2;
    char parfile[100];
    char matfile[100];
    char outfile[100];

    /*all the filenames will be generated from the same input string*/
    sprintf(parfile, "%s.par", filebase);
    sprintf(matfile, "%s.mat", filebase);
    sprintf(outfile, "%s_%g.gwy", filebase, wavelength*1e9);

    FILE *fw = fopen(parfile, "w");
    if (!fw) fprintf(stderr, "Error opening file %s\n", parfile);

    /*basic settings*/
    fprintf(fw, "POOL\n");
    fprintf(fw, "%d %d %d %g %g %g\n\n", xres, yres, zres, dx, dy, dz);

    fprintf(fw, "COMP\n");
    fprintf(fw, "%d\n\n", comp);

    fprintf(fw, "VERBOSE\n");
    fprintf(fw, "%d\n\n", verbose);

    fprintf(fw, "THREADS\n");
    fprintf(fw, "%d\n\n", threads);

    /*PLCR/ADE multiplier*/
    fprintf(fw, "DT_MULT\n");
    fprintf(fw, "0.8\n\n");

    /*graphics card use*/
    fprintf(fw, "GPU\n");
    fprintf(fw, "0\n\n");

    fprintf(fw, "UGPU\n");
    fprintf(fw, "0\n\n");

    /*Total/scattered field source with a single infinite layer*/
    fprintf(fw, "SOURCE_LTSF\n");
    fprintf(fw, "%d %d %d %d %d %d 0 0 1.5708 1 %d 14.943 1 3209 0 1 %g 1\n\n", bnd, bnd, bnd, 
                  xres-bnd, yres-bnd, zres-bnd, hz, wavelength);

    fprintf(fw, "SOURCE_WAVELENGTH\n");
    fprintf(fw, "%g %g %g\n\n", wavelength, wavelength, wavelength);
 
    /*remove TSF boundary at the bottom*/
    fprintf(fw, "LTSF_SKIP\n");
    fprintf(fw, "kn\n\n");

    /*periodic boundaries in xy, if requested*/

    fprintf(fw, "MBOUNDARY_X0\n");
    fprintf(fw, "periodic %d\n\n", 4*bnd);

    fprintf(fw, "MBOUNDARY_XN\n");
    fprintf(fw, "periodic %d\n\n", xres - 4*bnd);

    fprintf(fw, "MBOUNDARY_Y0\n");
    fprintf(fw, "periodic %d\n\n", 4*bnd);

    fprintf(fw, "MBOUNDARY_YN\n");
    fprintf(fw, "periodic %d\n\n", yres - 4*bnd);


    /*simple absorbing boundaries everywhere*/
    fprintf(fw, "BOUNDARY_ALL\n");
    fprintf(fw, "liao\n\n");

    /*material file placement*/
    fprintf(fw, "MEDIUM_VECTOR\n");
    fprintf(fw, "%s\n\n", matfile);

    /*output settings for images*/
    fprintf(fw, "OUT_FILE\n");
    fprintf(fw, "%s\n\n", outfile);

    fprintf(fw, "OUT_IMAGE\n");
    fprintf(fw, "Ex 1000 -1 -1 %d ex_xy\n\n", hz-athickness+1);

    fprintf(fw, "OUT_IMAGE\n");
    fprintf(fw, "Ex 1000 %d -1 -1 ex_yz\n\n", xres/2);

    fprintf(fw, "OUT_IMAGE\n");
    fprintf(fw, "Ex 1000 -1 %d -1 ex_xz\n\n", yres/2);

    /*after simulation is converged, output set of ASCII data files with intensity cross-sections for further analysis*/
    fprintf(fw, "OUT_PLANE\n");
    fprintf(fw, "Ex 10 5000 6000 1 -1 -1 %d ex_xy\n\n", hz-athickness+1);

    fprintf(fw, "OUT_PLANE\n");
    fprintf(fw, "Ey 10 5000 6000 1 -1 -1 %d ey_xy\n\n", hz-athickness+1);

    /*volume output for Paraview*/
    fprintf(fw, "OUT_VOLUME\n");
    fprintf(fw, "Material 1000 0 1 2 mat.vtk\n\n");

    /*maximum value output*/
    fprintf(fw, "OUT_SUM\n");
    fprintf(fw, "Max 100 %d %d %d %d %d %d -1 0 0 0 xmax_%g.txt\n\n", 4*bnd, 4*bnd, 0, xres-4*bnd, yres-4*bnd, zres, wavelength*1e9);


    fclose(fw);
}
    
/*function outputting the material file: definition of material placement in the computational domain. This uses already existing list of top sphere and cylinder base positions*/
void  
print_matfile(char *filebase, int xres, int yres, int zres, double dx, double thickness, double period, double radius, double radius2, int hz, double *xpos, double *ypos, double *xtop, double *ytop, int npos)
{
    int i;
    char matfile[100];

    //material properties
    char aldata[] = "5 1.0000 0 2.0598e16 2.2876e14 5.2306 -0.51202 2.2694e15 3.2867e14 5.2704 0.42503 2.4668e15 1.7731e15";
    char audata[] = "5 1.1431 0 1.3202e16 1.0805e14 0.26698 -1.2371 3.8711e15 4.4642e14 3.0834 -1.0968 4.1684e15 2.3555e15";
    char agdata[] = "5 15.833 0 1.3861e16 4.5841e13 1.0171 -0.93935 6.6327e15 1.6666e15 15.797 1.8087 9.2726e17 2.3716e17";
    char sidata[] = "99 Si_111";

    sprintf(matfile, "%s.mat", filebase);


    FILE *fw = fopen(matfile, "w");
    if (!fw) fprintf(stderr, "Error opening file %s\n", matfile);


    for (i=0; i<npos; i++) {
         fprintf(fw, "7 %g %g %g %g %g %g %g %s\n", xtop[i], ytop[i], (double)hz-thickness/dx, xpos[i], ypos[i], (double)hz+5, radius/dx, sidata);
         fprintf(fw, "4 %g %g %g %g %s\n", xtop[i], ytop[i], (double)hz-thickness/dx, radius2/dx, audata);
    }
  
    fprintf(fw, "8 %g %g %g %g %g %g %s\n", 0.0, 0.0, (double)hz, (double)xres, (double)yres, (double)zres, sidata);
    fclose(fw);
}

/*lj potential between two particles for revise steps*/
static double
get_lj_potential_spheres(double ax, double ay, double bx, double by, double asize, double bsize)
{
    double sigma = 0.82*(asize+bsize);
    double dist = ((ax-bx)*(ax-bx)
                    + (ay-by)*(ay-by));

    if ((asize>0 && bsize>0) && dist > asize/100) {
        double s2 = sigma*sigma, s4 = s2*s2, s6 = s4*s2, s12 = s6*s6;
        double d3 = dist*dist*dist, d6 = d3*d3;
        return (asize)*2e-5*(s12/d6 - s6/d3); //corrected for particle size
    }
    //return (asize)*3e-5*(pow(sigma, 12)/pow(dist, 6) - pow(sigma, 6)/pow(dist, 3)); //corrected for particle size
    else return 0;
}

/*run a simple relaxation to create something more realistic, if requested*/
void revise(double *rx, double *ry, int ndata, int xres, int yres, int nsteps, double radius2)
{
   int k, m, i;
   double fx[17], fy[17];
   double vx[17], vy[17];
   double ax[17], ay[17];
   double rdisizes[17];
   double diff = 0.1;
   double timestep = 30, mass = 1; 

   for (k=0; k<ndata; k++) {
      rdisizes[k] = radius2+3;
      vx[k] = vy[k] = 0;
      ax[k] = ay[k] = 0;
   }

   for (i=0; i<nsteps; i++) {
        /*test succesive LJ steps on substrate*/
        for (k = 0; k < ndata; k++) {
            fx[k] = fy[k] = 0;
            /*calculate forces for all particles on substrate*/

            for (m = 0; m < ndata; m++) {
                if (m == k)
                    continue;

                fx[k] -= (get_lj_potential_spheres(rx[m], ry[m], rx[k]+diff, ry[k], rdisizes[k], rdisizes[m])
                              -get_lj_potential_spheres(rx[m], ry[m], rx[k]-diff, ry[k], rdisizes[k], rdisizes[m]))/2/diff;
                fy[k] -= (get_lj_potential_spheres(rx[m], ry[m], rx[k], ry[k]+diff, rdisizes[k], rdisizes[m])
                              -get_lj_potential_spheres(rx[m], ry[m], rx[k], ry[k]-diff, rdisizes[k], rdisizes[m]))/2/diff;

            }
        }

        for (k=0; k<ndata; k++)
        {
            /*move all particles*/
            rx[k] += vx[k]*timestep + 0.5*ax[k]*timestep*timestep;
            vx[k] += 0.5*ax[k]*timestep;
            ax[k] = fx[k]/mass;
            vx[k] += 0.5*ax[k]*timestep;
            vx[k] *= 0.9;
            if (fabs(vx[k])>0.01) vx[k] = 0;

            ry[k] += vy[k]*timestep + 0.5*ay[k]*timestep*timestep;
            vy[k] += 0.5*ay[k]*timestep;
            ay[k] = fy[k]/mass;
            vy[k] += 0.5*ay[k]*timestep;
            vy[k] *= 0.9;
            if (fabs(vy[k])>0.01) vy[k] = 0;

            if (rx[k]<rdisizes[k]) rx[k] = rdisizes[k];
            if (ry[k]<rdisizes[k]) ry[k] = rdisizes[k];
            if (rx[k]>(xres-rdisizes[k])) rx[k] = xres-rdisizes[k];
            if (ry[k]>(yres-rdisizes[k])) ry[k] = yres-rdisizes[k];
        }
  }

}

int
create_linear_model(double *xtop, double *ytop, double *xpos, double *ypos, int *xres, int *yres, int *zres, double period, double randomize, double sthickness, double thickness, double dx)
{
    int i, ax, ay, az, hz, n;
    
    /*calculate the needed size*/
    *xres = 3*period/dx + 80;
    *yres = 2*period/dx + 80;
    *zres = (sthickness+thickness)/dx + 40;

    ax = (*xres)/2;
    ay = (*yres)/2;
    az = (*zres) - sthickness/dx - thickness/2.0/dx;
    hz = (*zres) - sthickness/dx;

    //create the model
    printf("creating the model\n");

    n=0;
    xtop[n] = xpos[n] = ax - period/dx/2.0 + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
    ytop[n] = ypos[n] = ay + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
    n++;

    xtop[n] = xpos[n] = ax + period/dx/2.0 + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
    ytop[n] = ypos[n] = ay + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
    n++;

    return n;
}
int
create_triangular_model(double *xtop, double *ytop, double *xpos, double *ypos, int *xres, int *yres, int *zres, double period, double randomize, double sthickness, double thickness, double dx)
{
    int i, ax, ay, az, hz, n;

    /*calculate the needed size*/
    *xres = 3*(period)/dx + 80;
    *yres = 3*(sqrt(3)/2*period)/dx + 80;
    *zres = (sthickness+thickness)/dx + 40;

    ax = (*xres)/2;
    ay = (*yres)/2;
    az = (*zres) - sthickness/dx - thickness/2.0/dx;
    hz = (*zres) - sthickness/dx;

    //create the model
    printf("creating the model\n");
    n=0;

    xtop[n] = xpos[n] = ax + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
    ytop[n] = ypos[n] = ay + sqrt(3)/4.0*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
    n++;

    for (i=0; i<=1; i++) {
       xtop[n] = xpos[n] = ax + (double)i*period/dx - 0.5*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5) ;
       ytop[n] = ypos[n] = ay - sqrt(3)/4.0*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5) ;
       n++;
    }

    return n;
}

int
create_rectangular_model(double *xtop, double *ytop, double *xpos, double *ypos, int *xres, int *yres, int *zres, double period, double randomize, double sthickness, double thickness, double dx)
{
    int i, ax, ay, az, hz, n;
    
    /*calculate the needed size*/
    *xres = 3*period/dx + 80;
    *yres = 3*period/dx + 80;
    *zres = (sthickness+thickness)/dx + 40;

    ax = (*xres)/2;
    ay = (*yres)/2;
    az = (*zres) - sthickness/dx - thickness/2.0/dx;
    hz = (*zres) - sthickness/dx;

    //create the model
    printf("creating the model\n");

    n=0;
    xtop[n] = xpos[n] = ax - period/dx/2.0 + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
    ytop[n] = ypos[n] = ay - period/dx/2.0 + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
    n++;

    xtop[n] = xpos[n] = ax - period/dx/2.0 + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
    ytop[n] = ypos[n] = ay + period/dx/2.0 + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
    n++;

    xtop[n] = xpos[n] = ax + period/dx/2.0 + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
    ytop[n] = ypos[n] = ay - period/dx/2.0 + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
    n++;

    xtop[n] = xpos[n] = ax + period/dx/2.0 + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
    ytop[n] = ypos[n] = ay + period/dx/2.0 + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
    n++;

    return n;
}

int
create_circular_model(double *xtop, double *ytop, double *xpos, double *ypos, int *xres, int *yres, int *zres, double period, double randomize, double sthickness, double thickness, double dx, int nobjects)
{
    int i, ax, ay, az, hz, n;
    double radius;
    
    /*calculate the needed size*/
    radius = ((double)nobjects*period)/(2*M_PI);
    *xres = (3*radius+period)/dx + 80;
    *yres = (3*radius+period)/dx + 80;
    *zres = (sthickness+thickness)/dx + 40;

    ax = (*xres)/2;
    ay = (*yres)/2;
    az = (*zres) - sthickness/dx - thickness/2.0/dx;
    hz = (*zres) - sthickness/dx;

    //create the model
    printf("creating the model  %g\n", radius/dx);

    for (n=0; n<nobjects; n++) {
         xtop[n] = xpos[n] = ax + radius/dx*sin(2*M_PI*((double)n/((double)nobjects))) + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
         ytop[n] = ypos[n] = ay + radius/dx*cos(2*M_PI*((double)n/((double)nobjects))) + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
    }

    return n;
}

int
create_hexagonal_model_7(double *xtop, double *ytop, double *xpos, double *ypos, int *xres, int *yres, int *zres, double period, double randomize, double sthickness, double thickness, double dx)
{
    int i, ax, ay, az, hz, n;
    
    /*calculate the needed size*/
    *xres = 3*(period)/dx + 80;
    *yres = 4*(sqrt(3)/2*period)/dx + 80;
    *zres = (sthickness+thickness)/dx + 40;

    ax = (*xres)/2;
    ay = (*yres)/2;
    az = (*zres) - sthickness/dx - thickness/2.0/dx;
    hz = (*zres) - sthickness/dx;

    //create the model
    printf("creating the model\n");
    n=0;
    for (i=-1; i<=1; i++) {
       xtop[n] = xpos[n] = ax + (double)i*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
       ytop[n] = ypos[n] = ay + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
       n++;
    }
    for (i=0; i<=1; i++) {
       xtop[n] = xpos[n] = ax + (double)i*period/dx - 0.5*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5) ;
       ytop[n] = ypos[n] = ay - sqrt(3)/2.0*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5) ;
       n++;
    }

    for (i=0; i<=1; i++) {
       xtop[n] = xpos[n] = ax + (double)i*period/dx - 0.5*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5) ;
       ytop[n] = ypos[n] = ay + sqrt(3)/2.0*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5) ;
       n++;
    }

    return n;
}


int
create_hexagonal_model_17(double *xtop, double *ytop, double *xpos, double *ypos, int *xres, int *yres, int *zres, double period, double randomize, double sthickness, double thickness, double dx)
{
    int i, ax, ay, az, hz, n;
    
    /*calculate the needed size*/
    *xres = 5*(period)/dx + 80;
    *yres = 6*(sqrt(3)/2*period)/dx + 80;
    *zres = (sthickness+thickness)/dx + 40;

    ax = (*xres)/2;
    ay = (*yres)/2;
    az = (*zres) - sthickness/dx - thickness/2.0/dx;
    hz = (*zres) - sthickness/dx;

    //create the model
    printf("creating the model\n");
    n=0;
    for (i=-1; i<=1; i++) {
       xtop[n] = xpos[n] = ax + (double)i*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
       ytop[n] = ypos[n] = ay + randomize*(((double)rand())/((double)RAND_MAX) - 0.5);
       n++;
    }
    for (i=-1; i<=2; i++) {
       xtop[n] = xpos[n] = ax + (double)i*period/dx - 0.5*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5) ;
       ytop[n] = ypos[n] = ay - sqrt(3)/2.0*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5) ;
       n++;
    }

    for (i=-1; i<=2; i++) {
       xtop[n] = xpos[n] = ax + (double)i*period/dx - 0.5*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5) ;
       ytop[n] = ypos[n] = ay + sqrt(3)/2.0*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5) ;
       n++;
    }
    for (i=-1; i<=1; i++) {
       xtop[n] = xpos[n] = ax + (double)i*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5) ;
       ytop[n] = ypos[n] = ay - 2.0*sqrt(3)/2.0*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5) ;
       n++;
    }
    for (i=-1; i<=1; i++) {
       xtop[n] = xpos[n] = ax + (double)i*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5) ;
       ytop[n] = ypos[n] = ay + 2.0*sqrt(3)/2.0*period/dx + randomize*(((double)rand())/((double)RAND_MAX) - 0.5) ;
       n++;
    }

    return n;
}


int main(int argc, char *argv[])
{
    int i, n, xres, yres, zres, hz, relaxsteps, model;
    double dx, period, thickness, radius, radius2;    
    int comp, threads, substrate;
    double sthickness;
    double wavelength;
    double xpos[17];
    double ypos[17];
    double xtop[17];
    double ytop[17];
    int npos = 17;
    double randomize; 
    
    wavelength = 633e-9; //source wavelength
    comp = 6000;  //number of computation steps
    threads = 2;  //max number of threads
    sthickness = 140e-9; //substrate thickness


    if (argc<10) {
      fprintf(stderr, "syntax: ./sersgen filename height periodicity radius radius2 dx randomize relax model\n");
      fprintf(stderr, "where the individual parameters have the following meaning:\n");
      fprintf(stderr, "   filename: basis of all the output filenames (without suffix)\n");
      fprintf(stderr, "   height: height of the rods in nm\n");
      fprintf(stderr, "   periodicity: rod-rod closest distance in nm\n");
      fprintf(stderr, "   radius: cylinder radius in nm\n");
      fprintf(stderr, "   radius2: sphere radius in nm\n");
      fprintf(stderr, "   dx: voxel spacing in nm\n");
      fprintf(stderr, "   randomize: randomization factor in voxels\n");
      fprintf(stderr, "   relax: number of relax steps\n");
      fprintf(stderr, "   model: used model: 0: linear, 1: triangular: 2: rectangular, 3: pentagonal, 4: six-rod circle, 5: hexagonal 7, 6: hexagonal 17\n");



      return 0;
    }

    /*load parameters and convert them to meters*/
    dx = atof(argv[6])*1e-9;
    thickness = atof(argv[2])*1e-9;
    period = atof(argv[3])*1e-9;
    radius = atof(argv[4])*1e-9;
    radius2 = atof(argv[5])*1e-9;
    randomize = atof(argv[7]);
    relaxsteps = atoi(argv[8]);
    model = atoi(argv[9]);

 
    /*initialize random numbers generator*/
    srand(1);

    //create positions 
    if (model==0) npos = create_linear_model(xtop, ytop, xpos, ypos, &xres, &yres, &zres, period, randomize, sthickness, thickness, dx);
    else if (model==1) npos = create_triangular_model(xtop, ytop, xpos, ypos, &xres, &yres, &zres, period, randomize, sthickness, thickness, dx);
    else if (model==2) npos = create_rectangular_model(xtop, ytop, xpos, ypos, &xres, &yres, &zres, period, randomize, sthickness, thickness, dx);
    else if (model==3) npos = create_circular_model(xtop, ytop, xpos, ypos, &xres, &yres, &zres, period, randomize, sthickness, thickness, dx, 5);
    else if (model==4) npos = create_circular_model(xtop, ytop, xpos, ypos, &xres, &yres, &zres, period, randomize, sthickness, thickness, dx, 6);
    else if (model==5) npos = create_hexagonal_model_7(xtop, ytop, xpos, ypos, &xres, &yres, &zres, period, randomize, sthickness, thickness, dx);
    else if (model==6) npos = create_hexagonal_model_17(xtop, ytop, xpos, ypos, &xres, &yres, &zres, period, randomize, sthickness, thickness, dx);
    else {
       printf("Error: unsupported model!\n");
       return 1;
    }

    hz = (zres) - sthickness/dx;


    //let the sphere positions relax if needed
    if (relaxsteps) {
       printf("relaxing positions\n");
       revise(xtop, ytop, npos, xres, yres, relaxsteps, radius2/dx); 
    }

    printf("dimensions: %d x %d x %d    hz %d  size %d\n", xres, yres, zres, hz, (int)(thickness/dx));

    printf("output everything\n");

    //output everything
    print_parfile(xres, yres, zres,
              dx, dx, dx,
              comp, threads, argv[1],
              hz, thickness/dx, wavelength, period);

    print_matfile(argv[1], xres, yres, zres, dx, thickness, period, radius, radius2, hz, xpos, ypos, xtop, ytop, npos);

}