#include<bibli_fonctions.h>

double** p;
int Np = 100;

void init_charges() {
  int i;
  double r, th;
  p = (double**)malloc(Np * sizeof(double*));
  for(i = 0; i < Np; i++) {
    p[i] = (double*)malloc(2 * sizeof(double));
    r = sqrt(drand48());
    th = 2.*M_PI * drand48();
    p[i][0] = r * cos(th);
    p[i][1] = r * sin(th);
  }
}

void systeme(double* q, double t, double* qp, int n) {
  int i;
  double* r = (double*)malloc(Np * sizeof(double));
  for(i = 0; i < Np; i++)
    r[i] = sqrt((q[0]-p[i][0])*(q[0]-p[i][0]) + (q[1]-p[i][1])*(q[1]-p[i][1])); 
  qp[0] = q[2];
  qp[1] = q[3];
  qp[2] = 0; qp[3] = 0;
  for(i = 0; i < Np; i++) {
    qp[2] += (q[0] - p[i][0])/r[i]/r[i]/r[i];
    qp[3] += (q[1] - p[i][1])/r[i]/r[i]/r[i];
  }
  free(r);
}

int main() {
  int i, n = 4, Nt = 10000;
  double t = 0, tfin = 10, dt = (tfin - t) / (Nt - 1);
  double* q = (double*)malloc(n * sizeof(double));
  fstream fich("solution_charges.res", ios::out);
  srand48(time(NULL));
  init_charges();
  q[0] = 0; q[1] = 0; q[2] = 0; q[3] = 0;
  for (i = 0; i < Np; i++)
    fich << i << " " << p[i][0] << " " << p[i][1] << endl;
  for (i = 0; i < Nt; i++) {
    fich << t << " " << q[0] << " " << q[1] << endl;
    if(q[0]*q[0]+q[1]*q[1] > 4) break;
    rk4(systeme, q, t, dt, n);
    t += dt;
  }
  for(i = 0; i < Np; i++)
    free(p[i]);
  free(p);
  free(q);
  fich.close();
  ostringstream pyth;
  pyth << "A = loadtxt('solution_charges.res')\n"
       << "plot(A[:" << Np << ",1], A[:" << Np << ",2], 'o')\n"
       << "plot(A[" << Np << ":,1], A[" << Np << ":,2])\n"
       << "x = linspace(-1,1,100)\n"
       << "plot(x,sqrt(1-x*x),'r')\n"
       << "plot(x,-sqrt(1-x*x),'r')\n"
       << "xlim(-1,1)\n"
       << "ylim(-1,1)\n";
  make_plot_py(pyth);
  return 0;
}