/* Mathematiques de l'informatique graphique    */
/* Sphere                                       */
/*                                              */
/* Auteur: Nicolas JANEY                        */
/* nicolas.janey@univ-fcomte.fr                 */
/* Novembre 2011                                */

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

#include "Sphere.h"
#include "Position3D.h"
#include "Direction3D.h"
#include "Rayon.h"

Sphere::Sphere(void) {
  centre = new Position3D();
  rayon = 1.0;
}

Sphere::Sphere(double cx,double cy,double cz,double rayon) {
  centre = new Position3D(cx,cy,cz);
  this->rayon = rayon;
}

Sphere::Sphere(Position3D *c,double rayon) {
  centre = new Position3D(c);
  this->rayon = rayon;
}

Sphere::~Sphere(void) {
  delete(centre);
}

bool Sphere::intersection(Rayon *r) {
  double dpx = centre->c[0]-r->pos->c[0];
  double dpy = centre->c[1]-r->pos->c[1];
  double dpz = centre->c[2]-r->pos->c[2];
  double a = r->dir->c[0]*r->dir->c[0]+
             r->dir->c[1]*r->dir->c[1]+
             r->dir->c[2]*r->dir->c[2];
  double b = -2.0*(r->dir->c[0]*dpx+
                   r->dir->c[1]*dpy+
                   r->dir->c[2]*dpz);
  double c = dpx*dpx+dpy*dpy+dpz*dpz-rayon*rayon;
  double delta = b*b-4.0*a*c;
  if ( delta < 0.0 )
    return(false);
  if ( delta == 0.0 ) {
    double t = -b/(2.0*a);
    return(t > 0.0); }
  double t1 = (-b-sqrt(delta))/(2.0*a);
  if ( t1 > 0.0 )
    return(true);
  double t2 = (-b+sqrt(delta))/(2.0*a);
  return(t2 > 0.0);
}

bool Sphere::intersection(Rayon *r,double *tt) {
  double dpx = centre->c[0]-r->pos->c[0];
  double dpy = centre->c[1]-r->pos->c[1];
  double dpz = centre->c[2]-r->pos->c[2];
  double a = r->dir->c[0]*r->dir->c[0]+
             r->dir->c[1]*r->dir->c[1]+
             r->dir->c[2]*r->dir->c[2];
  double b = -2.0*(r->dir->c[0]*dpx+
                   r->dir->c[1]*dpy+
                   r->dir->c[2]*dpz);
  double c = dpx*dpx+dpy*dpy+dpz*dpz-rayon*rayon;
  double delta = b*b-4.0*a*c;
  if ( delta < 0.0 )
    return(false);
  if ( delta == 0.0 ) {
    double t = -b/(2.0*a);
    if ( t > 0.0 ) {
      *tt = t;
      return(true); }
      else
      return(false); }
  double t1 = (-b-sqrt(delta))/(2.0*a);
  double t2 = (-b+sqrt(delta))/(2.0*a);
  double min = ( t1 < t2 ) ? t1 : t2;
  double max = ( t1 > t2 ) ? t1 : t2;
  if ( ( t1 <= 0.0 ) && ( t2 <= 0.0 ) )
    return(false);
  if ( t1 <= 0.0 )
    *tt = t2;
    else
    *tt = t1;
  return(true);
}