/* Particule de systeme de particules           */
/*                                              */
/* Auteur: Nicolas JANEY                        */
/* nicolas.janey@univ-fcomte.fr                 */
/* Octobre 2010                                 */

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

#include "Position3D.h"
#include "Direction3D.h"
#include "Particule.h"

static float random(void) {
  return(((float) rand()/RAND_MAX));
}

Particule::Particule(void) {
}

Particule::~Particule(void) {
}

double Particule::getNextEvent(double y) {
  double nextEvent = 1000000.0;
  double delta = by*by+-4.0*ay*(cy-y);
  double r = (-by-sqrt(delta))/2.0/ay;
  if ( r > 0.0 )
    nextEvent = r+start;
    else {
    r = (-ay+sqrt(delta))/2.0/ay;
    if ( r > 0.0 )
      nextEvent = r+start; }
  return(nextEvent);
}

void Particule::init(double time) {
  rebond = 0;
  start = startIni = time;
  ax = 2.0*random()+1.0;
  bx = random();
  ay = -5.0;
  by = 0.0;
  cy = 2.0*random();
  az = 2.0*random()+1.0;
  bz = random();
  nextEvent = getNextEvent(-(20.0+2*random()));
  update(time);
}

void Particule::update(double time) {
  if ( time > nextEvent ) {
    rebond++;
    double t = nextEvent-start;
    Position3D pr;
    pr.c[0] = ax*t+bx;
    pr.c[1] = ay*t*t+cy;
    pr.c[2] = az*t+bz;
    Direction3D dr;
    dr.c[0] = ax;
    dr.c[1] = 2.0*ay;
    dr.c[2] = az;
    ax = dr.c[0];
    ay = -5.0;
    az = dr.c[2];
    t = dr.c[1]/((3.0+3.0*random())*ay);
    by = 0.0;
    bx = pr.c[0]+ax*t;
    cy = pr.c[1]-ay*t*t;
    bz = pr.c[2]+az*t;
    start = nextEvent+t;
    nextEvent = getNextEvent(-20.0-rebond*(10.0+random())); }
  double t = time-startIni;
  if ( t > 10.0 )
    init(time);
  t = time-start;
  p.c[0] = ax*t+bx;
  p.c[1] = ay*t*t+cy;
  p.c[2] = az*t+bz;
  d.c[0] = ax;
  d.c[1] = 2.0*ay;
  d.c[2] = az;
  v = sqrt(d.c[0]*d.c[0]+d.c[1]*d.c[1]+d.c[2]*d.c[2]);
  d.c[0] /= v;
  d.c[1] /= v;
  d.c[2] /= v;
}