/**************************************************************************
  Copyright (C) 1992 Guy Moreillon
  All rights reserved.

  This software may be freely copied, modified, and redistributed
  provided that this copyright notice is preserved on all copies.

  You may not distribute this software, in whole or in part, as part of
  any commercial product without the express consent of the authors.

  There is no warranty or other guarantee of fitness of this software
  for any purpose.  It is provided solely "as is".
**************************************************************************/
/**************************************************************************

				RAD
		       (Interactive Radiosity)
		    	    Version 1.0
			       1992		    	    
		    	    
		  
		    	  Guy Moreillon
		    	  
**************************************************************************/

/**************************************************************************
  Fichier	: radio.c
  Description	: fonctions de calcul de la radiosite proprement dite
**************************************************************************/

#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <sys/time.h>
#include <malloc.h>
#include "types.h"
#include "macros.h"
#include "rad.h"
#include "ray.h"
#include "misc.h"
#include "radio.h"

#define MIN_B_DIFF	0.05
#define MAX_db		4000
#define MIN_E		0.01
#define MIN_CONTRIB	0.001
#define MIN_START_DIS	0.0001
#define MILLION		1000000

VECTOR	    db[MAX_db];		    /* Tableau des vecteurs db
	   				       utilise pour la subdivision */
VECTOR	    drb[MAX_db];	    /* Tableau des vecteurs 
					       delta received radiosity */
int	    nb_db;		    /* Nombre de vecteurs dans les tableaux */
float 	    Nb_Samples_Sq = (float )(NB_SAMPLES*NB_SAMPLES);
int	    nb_complete = 0;


void B_at_p(TRIANGLE	*shoot,
	    SAMPLE	*sample,
	    VECTOR	B)
/**************************************************************************
  But	: Calcule la radiosite au point sample de la patch shoot
  Entree: shoot :   la patch
	  sample:   l'echantillon
	  B     :   la radiosite (rendue)
  Sortie: neant
**************************************************************************/
{
  if (sample->el == NULL)
    {
      register	VERTEXP	vertex = shoot->tri.poly.vertex;

      V_interp(B, sample->c, vertex[0].dB,
    	          sample->a, vertex[1].dB);
      V_add_mul(B, B, sample->b, vertex[2].dB);
    }
  else
    {
      register	VERTEXP	    array   = shoot->tri.poly.vertex;
      register	short int   *vertex = sample->el->vertex;

      V_interp(B, sample->c, array[vertex[0]].dB,
    	          sample->a, array[vertex[1]].dB);
      V_add_mul(B, B, sample->b, array[vertex[2]].dB);
    };
}

void null_dB(TRIANGLE	*tri)
/**************************************************************************
  But	: Anulle la valeur des dB de la patch (apres la mise a feu de
	  celle-ci)
  Entree: tri		    : la patch
  Sortie: neant
**************************************************************************/
{
  register int	    i;
  register VERTEXP  vertex = tri->tri.poly.vertex;

  for( i = 0; i < tri->tri.poly.nb_v; i++)
    V_null(vertex[i].dB);
}

void comp_surf_el(SUBDIV    *subdiv)
/**************************************************************************
  But	: Calcule la surface de chaque elements d'une subdivision
  Entree: subdiv :  la subdivision
  Sortie: neant
**************************************************************************/
{
  int	    i;
  VECTOR    vab, vac, vh;
  VERTEXP   array = subdiv->dad->tri.poly.vertex;

  for(i = 0; i < 4; i++)
    {
      V_sub(vab, array[subdiv->elements[i].vertex[0]].point,
		 array[subdiv->elements[i].vertex[1]].point);
      V_sub(vac, array[subdiv->elements[i].vertex[0]].point,
		 array[subdiv->elements[i].vertex[2]].point);
      V_cross(vh, vac, vab);
      subdiv->elements[i].surface = V_norm(vh)*0.5;
    };
}

BOOLEAN point_in_el(TRIANGLEP	 tri,
		    SAMPLEP	 sample,
		    ELEMENTP	 el,
		    float 	*a,
		    float 	*b,
		    float 	*c)
/**************************************************************************
  But	: Determine si le sample est dans la subdivision el, et si oui,
	  rend a, b et c, les coordonnees barycentriques du point dans le
	  triangle
  Entree: tri :	    le polygone pere
	  sample :  le point d'echantillon
	  el :	    l'element
	  a, b, c : les coordonnees barycentriques
  Sortie: TRUE si sample est dans el, FALSE sinon
**************************************************************************/
{
  float     naxb, dum;
  VECTOR    ea, eb, axb, bxa, sa, sb, sc, p;
  VERTEXP   array = tri->tri.poly.vertex;

  register  short int	*vertex = el->vertex;

  V_sub(ea, array[vertex[1]].point,
	    array[vertex[0]].point);
  V_sub(eb, array[vertex[2]].point,
	    array[vertex[0]].point);
  V_sub(p, sample->point, array[vertex[0]].point);
  
  V_cross(axb, ea, eb);
  naxb = V_norm(axb);
  dum = 1.0/(naxb * naxb);

  V_cross(sa, eb, axb);
  V_scale(sa, dum, sa);

  *a = V_dot(p, sa);
  if ((*a < 0.0) || (*a > 1.0))
    return FALSE;

  V_cross(bxa, eb, ea);
  V_cross(sb, ea, bxa);
  V_scale(sb, dum, sb);

  *b = V_dot(p, sb);
  if ((*b < 0.0) || (*b > 1.0))
    return FALSE;

  V_add(sc, sa, sb);

  *c = V_dot(p, sc);
  if ((*c < 0.0) || (*c > 1.0))
    return FALSE;

  return TRUE;
}

void find_subdiv(TRIANGLEP  tri,
		 SAMPLEP    sample)
/**************************************************************************
  But	: Determine dans quelle subdivision de tri le point sample se
	  trouve. S'il n'est pas dans la subdivision, les donnees ne sont
	  pas changees.
  Entree: sample :  le point d'echantillon
	  tri :	    le triangle
  Sortie: neant
**************************************************************************/
{
  float  a, b, c;

  register  int	 i;

  for(i = 0; i < 4; i++)
    if (point_in_el(tri->tri.subdiv.dad,
		    sample,
		    &(tri->tri.subdiv.elements[i]),
		    &a, &b, &c))
      if (tri->tri.subdiv.elements[i].sons != NULL)
	find_subdiv(tri->tri.subdiv.elements[i].sons, sample);
      else
        {
	  sample->el = &(tri->tri.subdiv.elements[i]);
	  sample->a  = a;
	  sample->b  = b;
	  sample->c  = c;
	  break;
        };
}

short int create_vertex(TRIANGLEP   tri,
			short int   v1,
			short int   v2,
			float 	    id)
/**************************************************************************
  But	: Cree et rend l'indice du point compris entre v1 et v2.
  Entree: tri :	    le polygone pere
	  v1, v2 :  les deux sommets
	  id :	    l'identificateur du sommet a construire
  Sortie: l'indice du sommet
**************************************************************************/
{
  VECTOR  dum;
  register  int	    top = tri->tri.poly.nb_v;
  register  VERTEXP array;

  tri->tri.poly.nb_v++;
 
  tri->tri.poly.vertex =
    (VERTEXP)arr_realloc(tri->tri.poly.vertex, tri->tri.poly.nb_v,
			 &(tri->tri.poly.nb_tot), sizeof(VERTEX),
			 VERT_BLOCK_SIZE); 
  assert(tri->tri.poly.vertex);
  
  array = tri->tri.poly.vertex;

  array[top].id = id;
  V_add(dum, array[v1].point, array[v2].point);
  V_scale(array[top].point, 0.5, dum);
  V_add(dum, array[v1].B, array[v2].B);
  V_scale(array[top].B, 0.5, dum);
  V_add(dum, array[v1].dB, array[v2].dB);
  V_scale(array[top].dB, 0.5, dum);
  V_add(dum, array[v1].rB, array[v2].rB);
  V_scale(array[top].rB, 0.5, dum);

  array[top].update = FALSE;

  return top;
}

short int get_vertex(TRIANGLE	*tri,
		     short int	v1,
		     short int	v2)
/**************************************************************************
  But	: Rend l'indice du point compris entre v1 et v2. Si le point
	  n'existe pas encore dans la table, il le construit
  Entree: tri :	    le polygone pere
	  v1, v2 :  les deux sommets
  Sortie: l'indice du sommet
**************************************************************************/
{
  float     id, id1, id2, n;

  register  int	    i;
  register  VERTEXP vertex = tri->tri.poly.vertex;

  if (vertex[v1].id < vertex[v2].id)
    id1 = vertex[v1].id, id2 = vertex[v2].id - 1.0;
  else
    id1 = vertex[v2].id, id2 = vertex[v1].id - 1.0;

  n = (float)(id1 + id2);
  id = n*(n+1.0)*0.5 - id2;

  for(i = tri->tri.poly.nb_v - 1; i >= 0; i--)
    if (id == vertex[i].id)
      break;

  if (i >= 0)
    return i;
  else
    return (create_vertex(tri, v1, v2, id));
}

TRIANGLE *create_subdiv(TRIANGLEP   tri,
			short int   v[3])
/**************************************************************************
  But	: Cree une subdivision au triangle tri, en fonction des 3 vertices
	  specifiees
  Entree: tri :	    le triangle
	  v :	    le tableau des indices des sommets
  Sortie: le pointeur sur la subdivision
**************************************************************************/
{
  TRIANGLEP loc;
  short int vertex[3];

  register  int	i, j;
  register  ELEMENT *elements;

  loc = (TRIANGLE *)malloc(sizeof(TRIANGLE));
  assert(loc);

  nb_tri += 3;

  loc->type = subdiv_type;
  loc->dad  = tri->dad;
  if (tri->type == poly_type)
    loc->tri.subdiv.dad = tri;
  else
    loc->tri.subdiv.dad = tri->tri.subdiv.dad;
  elements = loc->tri.subdiv.elements;

  vertex[0] = get_vertex(loc->tri.subdiv.dad, v[0], v[1]);
  vertex[1] = get_vertex(loc->tri.subdiv.dad, v[0], v[2]);
  vertex[2] = get_vertex(loc->tri.subdiv.dad, v[1], v[2]);

  elements[0].vertex[0] = v[0];
  elements[0].vertex[1] = vertex[0];
  elements[0].vertex[2] = vertex[1];

  elements[1].vertex[0] = v[1];
  elements[1].vertex[1] = vertex[2];
  elements[1].vertex[2] = vertex[0];

  elements[2].vertex[0] = v[2];
  elements[2].vertex[1] = vertex[1];
  elements[2].vertex[2] = vertex[2];

  elements[3].vertex[0] = vertex[2];
  elements[3].vertex[1] = vertex[1];
  elements[3].vertex[2] = vertex[0];

  comp_surf_el(&(loc->tri.subdiv));

  elements[0].sons = NULL;
  elements[1].sons = NULL;
  elements[2].sons = NULL;
  elements[3].sons = NULL;

  for(i = 0; i < NB_SAMPLES; i++)
    for(j = 0; j < NB_SAMPLES; j++)
      find_subdiv(loc,
		  &(loc->tri.subdiv.dad->tri.poly.samples->array[i][j]));

  return loc;
}

void sum_dB(TRIANGLE	*tri,
	    VECTOR	s_dB)
/**************************************************************************
  But	: Calcule la somme des dB d'un polygone
  Entree: tri   : le polygone
	  s_dB	: la somme des radiosites
  Sortie: neant
**************************************************************************/
{
  register  int	    i;
  register  VERTEXP vertex = tri->tri.poly.vertex;

  V_null(s_dB);
  for(i = 0; i < tri->tri.poly.nb_v; i++)
    {
      VECTOR	dum;

      V_abs(dum, vertex[i].dB);	
      V_add_self(s_dB, dum);
    }
}

float comp_newE(TRIANGLEP   tri)
/**************************************************************************
  But	: Calcule la valeur d'emmission totale pour un polygone
  Entree: tri	: le polygone
  Sortie: la valeur d'emmission totale
**************************************************************************/
{
  float   inv;
  VECTOR  dum;

  sum_dB(tri, dum);
  inv = 1.0/(float)tri->tri.poly.nb_v;
  V_scale(dum, inv, dum);

  return (SumV(dum) * tri->tri.poly.surface);
}

void radio(SCENE_OP	scene,
	   TRIANGLEP	shoot,
	   VECTOR	point,
	   VECTOR	normal,
	   float	ro,
	   VECTOR	color,
	   VECTOR	B,
	   VECTOR	rB)
/**************************************************************************
  But	: Calcule les effets de la shooting patch sur le point passe
	  en parametre.
  Entree: scene	    : la scene
	  shoot	    : la shooting patch
	  point	    : le point
	  normal    : la normale a la surface en ce point
	  ro	    : le facteur de reflexion de la surface
	  color	    : la couleur de la surface
	  B	    : la radiosite resultante
	  rB	    : la radiosite recue par le point
  Sortie: neant
**************************************************************************/
{
  VECTOR    r, BB, dd, p, dum;
  float     lr, ilr, cos1, cos2, fact;
  float     A2_d_n = shoot->tri.poly.surface/Nb_Samples_Sq;

  register  int	i, j;

  V_null(dd);

  for(i = 0; i < NB_SAMPLES; i++)
    for(j = 0; j < NB_SAMPLES; j++)
      {
	V_copy(p, shoot->tri.poly.samples->array[i][j].point);
        V_sub(r, p, point);
        lr = V_norm(r);
        ilr = 1.0/lr;
        V_scale(r, ilr, r);
        cos1 = V_dot(r, normal);
        if (cos1 <= 0.0)
	  continue;
        cos2 = -V_dot(r, shoot->tri.poly.normal);
        if (cos2 <= 0.0)
	  continue;
        fact = cos1 * cos2 / (M_PI*lr*lr + A2_d_n);
        B_at_p(shoot, &(shoot->tri.poly.samples->array[i][j]), BB);
        V_scale(dum, fact, BB);
	if ((_ABS(dum[0]) < /*MIN_CONTRIB*/0.0) &&
	    (_ABS(dum[1]) < /*MIN_CONTRIB*/0.0) &&
	    (_ABS(dum[2]) < /*MIN_CONTRIB*/0.0))
	  nb_int_cuts++;
	else
	  {
	    RAY	ray;

	    V_copy(ray.dir, r);
	    V_add_mul(ray.org, point, MIN_START_DIS, ray.dir);
	    if (!shadow(&ray, scene, lr - 2 * MIN_START_DIS))
	      V_add_self(dd, dum);
	  }
      };

  V_scale(rB, A2_d_n, dd);  
  V_scale(B, ro, rB);
  V_mul(B, B, color);
}

#define get_indexes(I, E, V) \
  switch (I) {\
    case 0 : E = 0; V = 1; break;\
    case 1 : E = 0; V = 2; break;\
    case 2 : E = 1; V = 1; break;\
    default: Bye_bye("The compiler went nuts !!!");\
  }

void complete_rad(TRIANGLE  *tri,
		  short int v[3])
/**************************************************************************
  But	: Complete les valeurs de radiosite dans la subdivision passee
  Entree: tri   : la subdivision
	  v	: les indices des sommets du triangle contenant la
		  subdivision
  Sortie: neant
**************************************************************************/
{
  VECTOR    B, rB;
  VERTEXP   array;

  register  int	    i, ee, vv;
  register  ELEMENT *elements;
 
  if (tri == NULL)
    return;

  nb_complete++;

  array    = tri->tri.subdiv.dad->tri.poly.vertex;
  elements = tri->tri.subdiv.elements;

  for(i = 0; i < 3; i++)
    {
      get_indexes(i, ee, vv);
      if (array[elements[ee].vertex[vv]].update)
	continue;

      switch (i) {
	case 0 : V_add(B,  db[array[v[0]].db],
			   db[array[v[1]].db]);
		 V_add(rB, drb[array[v[0]].db],
			   drb[array[v[1]].db]);
		 break;
	case 1 : V_add(B,  db[array[v[0]].db],
			   db[array[v[2]].db]);
		 V_add(rB, drb[array[v[0]].db],
			   drb[array[v[2]].db]);
		 break;
	case 2 : V_add(B,  db[array[v[1]].db],
			   db[array[v[2]].db]);
		 V_add(rB, drb[array[v[1]].db],
			   drb[array[v[2]].db]);
		 break;
	default: Bye_bye("The compiler went nuts !!!");
      };

      V_scale(db[nb_db], 0.5, B);
      V_scale(drb[nb_db], 0.5, rB);
      array[elements[ee].vertex[vv]].db = nb_db;
      nb_db++;
      assert(nb_db < MAX_db);

      array[elements[ee].vertex[vv]].update = TRUE;
    };

  for(i = 0; i < 4; i++)
    complete_rad(elements[i].sons, elements[i].vertex);
}

BOOLEAN poly_in_front(TRIANGLEP tri, TRIANGLEP shoot)
/**************************************************************************
  But	: Determine si un polygone se trouve devant la shooting patch
  Entree: tri	    : le polygone
	  shoot	    : la shooting patch
  Sortie: TRUE si tri est devant shoot
**************************************************************************/
{
  int	i;

  for(i = 0; i < 3; i++)
    if (V_dot(shoot->tri.poly.normal, tri->tri.poly.vertex[i].point) >
	-shoot->tri.poly.d)
      return TRUE;

  nb_nif_tests++;
  return FALSE;
}

BOOLEAN object_in_front(OBJECT_OP object, TRIANGLEP shoot)
/**************************************************************************
  But	: Determine si la bounding box d'un objet se trouve devant la
	  shooting patch
  Entree: object    : l'objet
	  shoot	    : la shooting patch
  Sortie: TRUE si object est devant shoot
**************************************************************************/
{
  VECTOR p;
  int	 i;

  for(i = 0; i < 8; i++)
    {
      V_init(p, object->bbox[((i>>2)&1)][0],
		object->bbox[((i>>1)&1)][1],
		object->bbox[(i&1)][2]);

      if (V_dot(shoot->tri.poly.normal, p) > -shoot->tri.poly.d)
        return TRUE;
    }

  nb_nif_tests++;
  return FALSE;
}

void fire_on_subdiv(SCENE_OP	scene,
		    TRIANGLEP	shoot,
		    TRIANGLEP	tri,
		    BOOLEAN	cond[3],
		    short int	v[3],
		    float	min_size_s)
/**************************************************************************
  But	: Calcule les effets de la shooting patch sur la subdivision passee
	  en parametre. Ne recalcule la radiosite d'un nouveau point que si
	  sa condition est vraie (cond[i])
  Entree: scene	    : la scene
	  shoot	    : la shooting patch
	  tri	    : la subdivision (pointeur)
	  cond	    : les conditions des points
	  v	    : les 3 coins du triangle pere
	  min_size_s: la taille min que peut atteindre une subdivision,
		      au carre
  Sortie: neant
**************************************************************************/
{
  VECTOR    B, rB, BB, rBB, d[3];
  BOOLEAN   c[3];
  VERTEXP   array = tri->tri.subdiv.dad->tri.poly.vertex;

  register  int	    i, ee, vv;
  register  ELEMENT *elements = tri->tri.subdiv.elements;

  for(i = 0; i < 3; i++)
    {
      get_indexes(i, ee, vv);
      if (array[elements[ee].vertex[vv]].update)
	continue;

      if (cond[i])
        radio(scene,
	      shoot,
	      array[elements[ee].vertex[vv]].point,
	      tri->tri.subdiv.dad->tri.poly.normal,
	      tri->dad->object.pure.ro, tri->dad->object.pure.color, B, rB);
      else
        {
	  switch (i) {
	    case 0 : V_add(BB,  db[array[v[0]].db],
			        db[array[v[1]].db]);
		     V_add(rBB, drb[array[v[0]].db],
				drb[array[v[1]].db]);
		     break;
	    case 1 : V_add(BB,  db[array[v[0]].db],
			        db[array[v[2]].db]);
		     V_add(rBB, drb[array[v[0]].db],
				drb[array[v[2]].db]);
		     break;
	    case 2 : V_add(BB,  db[array[v[1]].db],
			        db[array[v[2]].db]);
		     V_add(rBB, drb[array[v[1]].db],
				drb[array[v[2]].db]);
		     break;
	    default: Bye_bye("The compiler went nuts !!!");
	  };
	  V_scale(B, 0.5, BB);
	  V_scale(rB, 0.5, rBB);
        };

      V_copy(db[nb_db], B);
      V_copy(drb[nb_db], rB);
      array[elements[ee].vertex[vv]].db = nb_db;
      nb_db++;
      assert(nb_db < MAX_db);

      array[elements[ee].vertex[vv]].update = TRUE;
    };
      
  for(i = 0; i < 4; i++)
    {
      V_sub(d[0], db[array[elements[i].vertex[0]].db],
		  db[array[elements[i].vertex[1]].db]);
      V_sub(d[1], db[array[elements[i].vertex[0]].db],
		  db[array[elements[i].vertex[2]].db]);
      V_sub(d[2], db[array[elements[i].vertex[1]].db],
		  db[array[elements[i].vertex[2]].db]);
      V_abs(d[0], d[0]);
      V_abs(d[1], d[1]);
      V_abs(d[2], d[2]);

      c[0] = ((d[0][0] > MIN_B_DIFF) ||
	      (d[0][1] > MIN_B_DIFF) ||
	      (d[0][2] > MIN_B_DIFF));
      c[1] = ((d[1][0] > MIN_B_DIFF) ||
	      (d[1][1] > MIN_B_DIFF) ||
	      (d[1][2] > MIN_B_DIFF));
      c[2] = ((d[2][0] > MIN_B_DIFF) ||
	      (d[2][1] > MIN_B_DIFF) ||
	      (d[2][2] > MIN_B_DIFF));

      if (c[0] || c[1] || c[2])
	{
	  V_sub(d[0], array[elements[0].vertex[1]].point,
		      array[elements[0].vertex[2]].point);
	  V_sub(d[1], array[elements[0].vertex[1]].point,
		      array[elements[1].vertex[1]].point);
	  V_sub(d[2], array[elements[0].vertex[2]].point,
		      array[elements[1].vertex[1]].point);
	  if ((V_dot(d[0], d[0]) >= min_size_s) &&
	      (V_dot(d[1], d[1]) >= min_size_s) &&
	      (V_dot(d[2], d[2]) >= min_size_s))
	    {
	      if (elements[i].sons == NULL)
	        elements[i].sons = create_subdiv(tri, elements[i].vertex);
	      fire_on_subdiv(scene,
			     shoot,
			     elements[i].sons,
			     c,
			     elements[i].vertex,
			     min_size_s);
	      array = tri->tri.subdiv.dad->tri.poly.vertex;
	    };
	}
      else
	complete_rad(elements[i].sons, elements[i].vertex);
    };
}

void fire_on_patch(SCENE_OP	scene,
		   TRIANGLEP	tri,
		   float	*Em)
/**************************************************************************
  But	: Calcule les effets de la shooting patch sur la patch passee en
	  parametre, et rend la nouvelle valeur d'emission de la patch
	  (multipliee par sa surface -> energie effectivement renvoyee)
	  Cette derniere valeur n'est utilisee que pour determiner si
	  la patch doit etre la nouvelle shooting patch
  Entree: scene : la scene
	  tri   : la patch (pointeur)
	  Em    : la nouvelle emission (pointeur)
  Sortie: neant
**************************************************************************/
{
  TRIANGLEP shoot = scene->state.shooting_p.patch;
  VECTOR    B[3], rB[3], d[3];
  BOOLEAN   c[3];
  short int t_v[3];

  register  int	    i;
  register  VERTEXP vertex = tri->tri.poly.vertex;

  nb_db = 0;
  t_v[0] = 0;
  t_v[1] = 1;
  t_v[2] = 2;

  for(i = 0; i < 3; i++)
    {
      radio(scene,
	    shoot,
	    vertex[i].point, tri->tri.poly.normal,
	    tri->dad->object.pure.ro, tri->dad->object.pure.color, B[i], rB[i]);

      V_copy(db[nb_db], B[i]);
      V_copy(drb[nb_db], rB[i]);
      vertex[i].db = nb_db;
      nb_db++;
      assert(nb_db < MAX_db);

      vertex[i].update = TRUE;
    };

  V_sub(d[0], B[0], B[1]);
  V_sub(d[1], B[0], B[2]);
  V_sub(d[2], B[1], B[2]);
  V_abs(d[0], d[0]);
  V_abs(d[1], d[1]);
  V_abs(d[2], d[2]);

  c[0] = ((d[0][0] > MIN_B_DIFF) ||
	  (d[0][1] > MIN_B_DIFF) ||
	  (d[0][2] > MIN_B_DIFF));
  c[1] = ((d[1][0] > MIN_B_DIFF) ||
	  (d[1][1] > MIN_B_DIFF) ||
	  (d[1][2] > MIN_B_DIFF));
  c[2] = ((d[2][0] > MIN_B_DIFF) ||
	  (d[2][1] > MIN_B_DIFF) ||
	  (d[2][2] > MIN_B_DIFF));

  if (c[0] || c[1] || c[2])
    {
      if (tri->tri.poly.sons == NULL)
	tri->tri.poly.sons = create_subdiv(tri, t_v);
      fire_on_subdiv(scene,
		     shoot,
		     tri->tri.poly.sons,
		     c,
		     t_v,
		     scene->state.min_size_s);
      vertex = tri->tri.poly.vertex;
    }
  else
    complete_rad(tri->tri.poly.sons, t_v);

/*  On additionne les valeurs calculees ici pour que les calculs soient
    corrects lors de la creation de nouveaux sommets (voir create_vertex) */

  for(i = 0; i < tri->tri.poly.nb_v; i++)
    {
      V_add_self(vertex[i].B,  db[vertex[i].db]);
      V_add_self(vertex[i].dB, db[vertex[i].db]);
      V_add_self(vertex[i].rB, drb[vertex[i].db]);
      vertex[i].update = FALSE;
    };

  *Em = comp_newE(tri);
}

OBJECT_OP find_next_object(OBJECT_OP first, SCENE_OP scene)
/**************************************************************************
  But	: Trouve le prochain objet de la liste commencant a first qui soit
	  devant la shooting patch. Met a jour la new_shoot
  Entree: first	    : le premier objet de la liste
	  scene	    : la scene
  Sortie: l'objet suivant a traiter
**************************************************************************/
{
  SHOOT_P   shoot;
  OBJECT_OP obj;

  shoot.E	= 0.0;
  shoot.patch	= NULL;
  shoot.o	= -1;
  shoot.p	= -1;

  for(obj = first; obj != NULL; obj = obj->next)
    if (!object_in_front(obj, scene->state.shooting_p.patch))
      {
        find_shooting_p(&shoot, scene, obj);
	if (_ABS(shoot.E) > _ABS(scene->state.new_shoot.E))
	  scene->state.new_shoot = shoot;
      }
    else
      break;

  return obj;
}

BOOLEAN fire_on_object(SCENE_OP		scene,
		       OBJECT_OP	object,
		       int		depth,
		       struct timeval	*to)
/**************************************************************************
  But	: Calcule les effets de la shooting patch sur l'objet passe en
	  parametre, et rend l'etat du calcul apres un lapse de temps
	  specifie par l'heure de timeout
  Entree: scene	    : la scene
	  object    : l'objet
	  depth	    : la profondeur d'imbriquement s'il sagit d'un assemb.
	  to	    : l'heure de timeout
  Sortie: TRUE si le calcul n'est PAS termine, FALSE si on s'arrete a cause
	  du quota de temps
**************************************************************************/
{
  switch(object->type) {
    case assembly_type:
      {
	BOOLEAN timeout = FALSE;

	if (scene->state.complex[depth] == NULL)
	  scene->state.complex[depth] =
	    find_next_object(object->object.assembly.sons, scene);

	while((scene->state.complex[depth] != NULL) && !timeout)
	  {
	    timeout = fire_on_object(scene,
				     scene->state.complex[depth],
				     (depth+1),
				     to);

	    if (!timeout)
	      scene->state.complex[depth] =
		find_next_object(scene->state.complex[depth]->next, scene);
	  };

	return (timeout);
      }
      break;
    case composite_type:
      {
	BOOLEAN	timeout = FALSE;

	if (scene->state.pure == NULL)
	  scene->state.pure = 
	    find_next_object(object->object.composite.sons, scene);

	while((scene->state.pure != NULL) && !timeout)
	  {
	    timeout = fire_on_object(scene,
				     scene->state.pure,
				     0,
				     to);

	    if (!timeout)
	      scene->state.pure =
		find_next_object(scene->state.pure->next, scene);
	  }

	return (timeout);
      }
      break;
    case pure_type:
      {
	BOOLEAN	timeout = FALSE;

	TRIANGLEP   polys = object->object.pure.polys;

	while((scene->state.p < object->object.pure.nb_polys) && !timeout)
	  {
	    if (&(polys[scene->state.p]) != scene->state.shooting_p.patch)
	      {
	        float	E;

	        if (poly_in_front(&(polys[scene->state.p]), scene->state.shooting_p.patch))
	          {
	            struct timeval	tv_cur;
	            struct timezone	tz_cur;

	            fire_on_patch(scene,
			          &(polys[scene->state.p]),
			          &E);

	            if (_ABS(E) > _ABS(scene->state.new_shoot.E))
		      {
		        scene->state.new_shoot.E     = E;
		        scene->state.new_shoot.patch = &(polys[scene->state.p]);
		        scene->state.new_shoot.o     = object->object.pure.id;
		        scene->state.new_shoot.p     = scene->state.p;
		      }

	            gettimeofday(&tv_cur, &tz_cur);
	            timeout = timercmp(&tv_cur, to, >);	      
	          }
	        else
		  {
		    E = comp_newE(&(polys[scene->state.p]));
	            if (_ABS(E) > _ABS(scene->state.new_shoot.E))
		      {
		        scene->state.new_shoot.E     = E;
		        scene->state.new_shoot.patch = &(polys[scene->state.p]);
		        scene->state.new_shoot.o     = object->object.pure.id;
		        scene->state.new_shoot.p     = scene->state.p;
		      }
		  }
	      }

	    scene->state.p++;
	  }

	if (!timeout)
	  scene->state.p = 0;

	return (timeout);
      }
      break;
    default:
      break;
  }
}

BOOLEAN comp_rad(SCENE_OP scene, int time_quota)
/**************************************************************************
  But	: Evalue une nouvelle solution jusqu'a l'epuisement du quota de
	  temps ou jusqu'a la convergence
  Entree: time_quota : quota de temps imparti
  Sortie: TRUE si le calcul n'est pas termine, FALSE sinon
**************************************************************************/
{
  BOOLEAN	    not_done, ok = TRUE;
  struct timeval    tv_start, timeout;
  struct timezone   tz_dum;
  int		    tq = time_quota * 1000;

  gettimeofday(&tv_start, &tz_dum);
  timeout.tv_sec  = tv_start.tv_sec + (tv_start.tv_usec + tq)/MILLION;
  timeout.tv_usec = (tv_start.tv_usec + tq)%MILLION;

  while((not_done = (_ABS(scene->state.shooting_p.E) > MIN_E)) && ok)
    {
      ok = !fire_on_object(scene, scene->scene, 0, &timeout);

      if (ok)
	{
	  scene->state.o = 1;
          null_dB(scene->state.shooting_p.patch);
          scene->state.shooting_p = scene->state.new_shoot;
          scene->state.new_shoot.E = 0.0;
	}
    }

  if (ok)
    {
      scene->state.shooting_p.E	    = 0.0;
      scene->state.shooting_p.patch = NULL;
      scene->state.shooting_p.o	    = 0;
      scene->state.shooting_p.p	    = 0;
    }

  return (not_done);
}
