/* Auteur: Nicolas JANEY */
/* nicolas.janey@univ-fcomte.fr */
/* Avril 2001 */
/* Calcul des rayons transmis */
/* et reflechis sur une sphere */
#include <stdio.h>
#include <math.h>
typedef struct coord_3D {
float x ;
float y ;
float z ; } coord_3D ;
typedef struct direction_3D {
float dx ;
float dy ;
float dz ; } direction_3D ;
typedef struct rayon {
coord_3D pos ;
direction_3D dir ; } rayon ;
/* En entree le rayon r teste. En sortie
le nombre d'intersection(s) */
int intersectionSphere(rayon *r) {
float a = r->dir.dx*r->dir.dx +
r->dir.dy*r->dir.dy +
r->dir.dz*r->dir.dz ;
float b = 2*(r->dir.dx*r->pos.x +
r->dir.dy*r->pos.y +
r->dir.dz*r->pos.z) ;
float c = r->pos.x*r->pos.x +
r->pos.y*r->pos.y +
r->pos.z*r->pos.z - 1 ;
float delta = b*b - 4*a*c ;
if ( delta > 0 )
return(2) ;
if ( delta == 0 )
return(1) ;
return(0) ;
}
int posIntersectionSphere(rayon *r,
coord_3D *p) {
float a = r->dir.dx*r->dir.dx +
r->dir.dy*r->dir.dy +
r->dir.dz*r->dir.dz ;
float b = 2*(r->dir.dx*r->pos.x +
r->dir.dy*r->pos.y +
r->dir.dz*r->pos.z) ;
float c = r->pos.x*r->pos.x +
r->pos.y*r->pos.y +
r->pos.z*r->pos.z - 1 ;
float delta = b*b - 4*a*c ;
if ( delta > 0 ) {
float rcn =(float) pow((double) delta,0.5) ;
float d = (-b-rcn)/2/a ;
if ( d > 0 ) {
p->x = r->pos.x + d * r->dir.dx ;
p->y = r->pos.y + d * r->dir.dy ;
p->z = r->pos.z + d * r->dir.dz ;
return(1) ; }
else {
float d = (-b+rcn)/2/a ;
if ( d > 0 ) {
p->x = r->pos.x + d * r->dir.dx ;
p->y = r->pos.y + d * r->dir.dy ;
p->z = r->pos.z + d * r->dir.dz ;
return(1) ; } }
return(0) ; }
if ( delta == 0 ) {
float d = -b/2/a ;
p->x = r->pos.x + d * r->dir.dx ;
p->y = r->pos.y + d * r->dir.dy ;
p->z = r->pos.z + d * r->dir.dz ;
return(1) ; }
return(0) ;
}
float produitScalaire(direction_3D *d1,
direction_3D *d2) {
return(d1->dx*d2->dx +
d1->dy*d2->dy +
d1->dz*d2->dz) ;
}
void calculReflexion(rayon *r,
float d,
rayon *ref) {
ref->pos.x = r->pos.x + d * r->dir.dx ;
ref->pos.y = r->pos.y + d * r->dir.dy ;
ref->pos.z = r->pos.z + d * r->dir.dz ;
direction_3D n = { ref->pos.x,
ref->pos.y,
ref->pos.z } ;
direction_3D i = { -r->dir.dx,
-r->dir.dy,
-r->dir.dz } ;
float scal = produitScalaire(&n,&i) ;
ref->dir.dx = 2 * n.dx * scal - i.dx ;
ref->dir.dy = 2 * n.dy * scal - i.dy ;
ref->dir.dz = 2 * n.dz * scal - i.dz ;
}
/* En entree le rayon r teste. En sortie
un booleen indiquant l'existence d'un
rayon reflechi, le rayon reflechi ref */
int rayonReflechiSphere(rayon *r,
rayon *ref) {
float a = r->dir.dx*r->dir.dx +
r->dir.dy*r->dir.dy +
r->dir.dz*r->dir.dz ;
float b = 2*(r->dir.dx*r->pos.x +
r->dir.dy*r->pos.y +
r->dir.dz*r->pos.z) ;
float c = r->pos.x*r->pos.x +
r->pos.y*r->pos.y +
r->pos.z*r->pos.z - 1 ;
float delta = b*b - 4*a*c ;
a *= 2 ;
if ( delta > 0 ) {
float rcn =(float) pow((double) delta,0.5) ;
float d = (-b-rcn)/a ;
if ( d > 0 ) {
calculReflexion(r,d,ref) ;
return(1) ; }
else {
float d = (-b+rcn)/2 ;
if ( d > 0 ) {
calculReflexion(r,d,ref) ;
return(1) ; } }
return(0) ; }
if ( delta == 0 ) {
float d = -b/a ;
calculReflexion(r,d,ref) ;
return(1) ; }
return(0) ;
}
int main(int argc,char **argv) {
/* rayon r = { {3.0F,-0.2F,-0.3F},
{-0.95F,0.22F,0.22F} } ;*/
/* rayon r = { {3.0F,0.7071067811F,
0.7071067811F},
{-1.0F,0.0F,0.0F} } ;*/
/* rayon r = { {3.0F,0.7071067811F,0.0F},
{-1.0F,0.0F,0.0F} } ;*/
rayon r = { {3.0F,0.0F,0.7071067811F},
{-1.0F,0.0F,0.0F} } ;
printf("Rayon initial : %f %f %f\n",
r.pos.x,r.pos.y,r.pos.z) ;
printf(" %f %f %f\n",
r.dir.dx,r.dir.dy,r.dir.dz) ;
int n = intersectionSphere(&r) ;
printf("Nombre intersections : %d\n",n) ;
coord_3D p ;
int inter = posIntersectionSphere(&r,&p) ;
if ( inter ) {
printf("Pos intersection : %f %f %f\n",
p.x,p.y,p.z) ; }
else {
printf("Pas d'intersection\n") ; }
rayon ref ;
inter = rayonReflechiSphere(&r,&ref) ;
if ( inter ) {
printf("Rayon reflechi : %f %f %f\n",
ref.pos.x,ref.pos.y,ref.pos.z) ;
printf(" %f %f %f\n",
ref.dir.dx,
ref.dir.dy,
ref.dir.dz) ; }
else {
printf("Pas de reflexion\n") ; }
return(0);
}
RETOUR