//===========================================================
//===========================================================
//===========================================================
//
//  Molecular dynamics simulation of hardspheres
//
//===========================================================
//===========================================================
//===========================================================

#include <iostream>
#include <math.h>
#include <fstream>
#include <vector>
#include <time.h>
#include <string.h>
#include <stdlib.h>
#include <cstdlib>

#include "box.h"
#include "sphere.h"
#include "event.h"
#include "heap.h"


double *spherespp(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate)
{
	int eventspercycle = 20;              // # events per sphere per cycle
	double initialpf = 0.01;              // initial packing fraction
	double temp = 0.2;                    // initial temp., use 0 for zero velocities
	double growthrate = 0.001;            // growth rate
	double massratio = 1.;               // ratio of sphere masses
	int hardwallBC = 0;                   // 0 for periodic, 1 for hard wall BC

  double d, r;   // initial diameter and radius of spheres

  r = pow(initialpf*pow(SIZE, DIM)/(N*VOLUMESPHERE), 1.0/((double)(DIM)));

  box b(N, r, growthrate, maxpf, bidispersityratio, 
	bidispersityfraction, massratio, hardwallBC);

	b.CreateSpheres(temp);
  
  
  while ((b.collisionrate < maxcollisionrate) && (b.pf < maxpf) && (b.pressure < maxpressure)) 
    {
      b.Process(eventspercycle*N);
      b.Synchronize(true);
    }
  
	double *data = (double *)malloc(DIM * N * sizeof(double));

  b.WriteConfiguration(data);
  
  return data;
}

extern "C" {
	double *spheres(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate) {
		return spherespp(N, bidispersityratio, bidispersityfraction, maxpf, maxpressure, maxcollisionrate);
	}
}