#include<bibli_fonctions.h>

typedef struct {
  double masse, charge, *pos;
} particule;

particule* part;
int Np = 30;

void init_part(double* q) {
  int i;
  part = (particule*)malloc(Np * sizeof(particule));
  for(i = 0; i < Np; i++) {
    part[i].masse = 4.*drand48() + 1.;
    part[i].charge = 2.*drand48() + 1.;
    part[i].pos = &(q[2*i]);
    part[i].pos[0] = drand48();
    part[i].pos[1] = drand48();
  }
}

void systeme(double* q, double t, double* qp, int n) {
  int i, j, j0;
  double** r = (double**)malloc(Np * sizeof(double*));
  for(i = 0; i < Np; i++) {
    r[i] = (double*)malloc(Np * sizeof(double));
    for(j = 0; j < Np; j++)
      if(j > i)
	r[i][j] = sqrt(pow(part[i].pos[0]-part[j].pos[0], 2) 
		       + pow(part[i].pos[1]-part[j].pos[1], 2)); 
      else if(j < i)
	r[i][j] = r[j][i];
  }
  for(i = 0; i < n/2; i++)
    qp[i] = q[i+n/2];
  for(i = n/2; i < n; i+=2) {
    qp[i] = 0.; qp[i+1] = 0.;
    j0 = (i-n/2)/2;
    for(j = 0; j < Np; j++)
      if(j != j0) {
	qp[i] += part[j0].charge * part[j].charge / part[j0].masse
	         * (part[j0].pos[0]-part[j].pos[0]) / pow(r[j0][j], 3); 
	qp[i+1] += part[j0].charge * part[j].charge / part[j0].masse
	         * (part[j0].pos[1]-part[j].pos[1]) / pow(r[j0][j], 3); 
      }
  }
  for(i = 0; i < Np; i++)
    free(r[i]);
  free(r);
}

int main() {
  int i, j, n = 4*Np, Nt = 200;
  double t = 0, tfin = 0.2, dt = (tfin - t) / (Nt - 1);
  double* q = (double*)malloc(n * sizeof(double));
  fstream fich("particules.res", ios::out);
  srand48(time(NULL));
  init_part(q);
  for(i = n/2; i < n; i++) 
    q[i] = 0.;
  for(i = 0; i < Nt; i++) {
    for(j = 0; j < Np; j++)
      fich << part[j].pos[0] << " " << part[j].pos[1] << " ";
    fich << endl;
    rk4(systeme, q, t, dt, n);
    for(j = 0; j < n/2; j++) {
      if(q[j] < 0.) {
	q[j] = -q[j]; q[j+n/2] = -q[j+n/2];
      }
      if(q[j] > 1.) {
	q[j] = 2.-q[j]; q[j+n/2] = -q[j+n/2];
      }
      if(j%2 == 0 && sqrt(q[j+n/2]*q[j+n/2]+q[j+n/2+1]*q[j+n/2+1]) > 0.1/dt) { 
	q[j+n/2] = 0.; q[j+n/2+1] = 0.;
      }
    }
  }
  free(q); free(part);
  fich.close();
  ostringstream pyth;
  pyth << "A = loadtxt('particules.res')\n"
       << "for i in range(" << Nt-1 << "):\n"
       << "  clf()\n"
       << "  xlim(0,1)\n"
       << "  ylim(0,1)\n"
       << "  for j in range(" << Np << "):\n"
       << "    plot(A[i:i+2,2*j], A[i:i+2,2*j+1], linewidth=4,"
       <<          "color=cm.Spectral(j/" << Np << ".))\n"
       << "  pause(0.001)\n";
  make_plot_py(pyth);
  return 0;
}