--- ray/src/rt/fprism.c	1993/09/29 10:40:31	2.1
+++ ray/src/rt/fprism.c	2003/02/22 02:07:28	2.4
@@ -1,9 +1,6 @@
-/* Copyright (c) 1993 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: fprism.c,v 2.4 2003/02/22 02:07:28 greg Exp $";
 #endif
-
 /* Ce programme calcule les directions et les energies des rayons lumineux
    resultant du passage d'un rayon au travers d'un vitrage prismatique
  
@@ -24,7 +21,6 @@ double x;
 		return(0.);
 	return(sqrt(x));
 }
-#define sqrt(x)	Sqrt(x)
 
 /* definitions de macros utiles */
 
@@ -103,7 +99,7 @@ extern long eclock;
 
 /* Definition des routines */
 
-#define term(a,b) a/sqrt(a*a+b*b)
+#define term(a,b) a/Sqrt(a*a+b*b)
 static
 prepare_matrices()
 {
@@ -209,7 +205,7 @@ FVECT v,v_out;
  FVECT v_temp;
  double det;
 
- det = sqrt( (YY(v)*YY(v)+ZZ(v)*ZZ(v))*(YY(v)*YY(v)+ZZ(v)*ZZ(v))+
+ det = Sqrt( (YY(v)*YY(v)+ZZ(v)*ZZ(v))*(YY(v)*YY(v)+ZZ(v)*ZZ(v))+
          (XX(v)*XX(v)*YY(v)*YY(v))+(XX(v)*XX(v)*ZZ(v)*ZZ(v)) );
  XX(v_temp) = (YY(v)*YY(v)+ZZ(v)*ZZ(v))/det;
  YY(v_temp) = -( XX(v)*YY(v) )/det;
@@ -227,7 +223,7 @@ FVECT v,v_out;
  FVECT v_temp;
  double det;
 
- det = sqrt( (ZZ(v)*ZZ(v)+YY(v)*YY(v)) );
+ det = Sqrt( (ZZ(v)*ZZ(v)+YY(v)*YY(v)) );
  XX(v_temp) = 0.;
  YY(v_temp) = -ZZ(v)/det;
  ZZ(v_temp) = YY(v)/det;
@@ -367,7 +363,23 @@ TRAYON r_initial;
  
 
 
+static int
+compare(r1,r2,marge)
+TRAYON r1, r2;
+double marge;
 
+{
+ double arctg1, arctg2;
+
+ arctg1 = atan2(Y(r1),X(r1));
+ arctg2 = atan2(Y(r2),X(r2));
+ if ((arctg1 - marge <= arctg2) && (arctg1 + marge >= arctg2)) return 1;
+ else return 0;
+}
+
+
+
+
 static
 sortie(r)
 TRAYON r;
@@ -411,9 +423,9 @@ TRAYON r;
 {
  double det;
  
- det = sqrt(X(r)*X(r)+Y(r)*Y(r)+Z(r)*Z(r));
- sinus = sqrt(Y(r)*Y(r)+Z(r)*Z(r))/det;
- cosinus = sqrt(X(r)*X(r))/det;
+ det = Sqrt(X(r)*X(r)+Y(r)*Y(r)+Z(r)*Z(r));
+ sinus = Sqrt(Y(r)*Y(r)+Z(r)*Z(r))/det;
+ cosinus = Sqrt(X(r)*X(r))/det;
  if (r.n == 1.) rapport = prism.np * prism.np;
  else rapport = 1./(prism.np * prism.np);
  return;
@@ -432,7 +444,7 @@ TRAYON r_incident;
  X(r_reflechi) = -X(r_incident);
  Y(r_reflechi) = Y(r_incident);
  Z(r_reflechi) = Z(r_incident);
- if(sinus > sqrt(rapport) || r_incident.dest == tot_ref)
+ if(sinus > Sqrt(rapport) || r_incident.dest == tot_ref)
 	{
 	 r_reflechi.ppar1 = r_incident.ppar1;
 	 r_reflechi.pper1 = r_incident.pper1;
@@ -442,14 +454,14 @@ TRAYON r_incident;
 	}
  else
 	{ 
-	 r_reflechi.ppar1 = r_incident.ppar1*(rapport*cosinus-sqrt(rapport-
-		 (sinus*sinus)))/(rapport*cosinus+sqrt(rapport-(sinus*sinus)));
- 	 r_reflechi.pper1 = r_incident.pper1*(cosinus-sqrt
-	        (rapport-(sinus*sinus)))/(cosinus+sqrt(rapport-(sinus*sinus)));
- 	 r_reflechi.ppar2 = r_incident.ppar2*(rapport*cosinus-sqrt(rapport-
-		 (sinus*sinus)))/(rapport*cosinus+sqrt(rapport-(sinus*sinus)));
- 	 r_reflechi.pper2 = r_incident.pper2*(cosinus-sqrt
-	        (rapport-(sinus*sinus)))/(cosinus+sqrt(rapport-(sinus*sinus)));
+	 r_reflechi.ppar1 = r_incident.ppar1*(rapport*cosinus-Sqrt(rapport-
+		 (sinus*sinus)))/(rapport*cosinus+Sqrt(rapport-(sinus*sinus)));
+ 	 r_reflechi.pper1 = r_incident.pper1*(cosinus-Sqrt
+	        (rapport-(sinus*sinus)))/(cosinus+Sqrt(rapport-(sinus*sinus)));
+ 	 r_reflechi.ppar2 = r_incident.ppar2*(rapport*cosinus-Sqrt(rapport-
+		 (sinus*sinus)))/(rapport*cosinus+Sqrt(rapport-(sinus*sinus)));
+ 	 r_reflechi.pper2 = r_incident.pper2*(cosinus-Sqrt
+	        (rapport-(sinus*sinus)))/(cosinus+Sqrt(rapport-(sinus*sinus)));
 	 r_reflechi.e = r_incident.e *(((r_reflechi.ppar1*r_reflechi.ppar1+
          r_reflechi.pper1*r_reflechi.pper1)/(r_incident.ppar1*r_incident.ppar1+
 	 r_incident.pper1*r_incident.pper1))+((r_reflechi.ppar2*r_reflechi.ppar2
@@ -471,22 +483,22 @@ TRAYON r_incident;
 
  r_transmis = r_incident;
  trigo(r_incident); 
- if (sinus <= sqrt(rapport) && r_incident.dest != tot_ref) 
+ if (sinus <= Sqrt(rapport) && r_incident.dest != tot_ref) 
  	{
  	 X(r_transmis) = (X(r_incident)/(fabs(X(r_incident))))*
-			 (sqrt(1.-(Y(r_incident)*Y(r_incident)+Z(r_incident)*
+			 (Sqrt(1.-(Y(r_incident)*Y(r_incident)+Z(r_incident)*
 	         	           Z(r_incident))/rapport));
- 	 Y(r_transmis) = Y(r_incident)/sqrt(rapport);
-	 Z(r_transmis) = Z(r_incident)/sqrt(rapport);
-	 r_transmis.ppar1 = r_incident.ppar1*2.*sqrt(rapport)*cosinus/
-		 	   (sqrt(rapport-sinus*sinus)+rapport*cosinus);
-	 r_transmis.pper1 = r_incident.pper1*2.*cosinus/(cosinus+sqrt(rapport
+ 	 Y(r_transmis) = Y(r_incident)/Sqrt(rapport);
+	 Z(r_transmis) = Z(r_incident)/Sqrt(rapport);
+	 r_transmis.ppar1 = r_incident.ppar1*2.*Sqrt(rapport)*cosinus/
+		 	   (Sqrt(rapport-sinus*sinus)+rapport*cosinus);
+	 r_transmis.pper1 = r_incident.pper1*2.*cosinus/(cosinus+Sqrt(rapport
 			   - sinus*sinus));
-	 r_transmis.ppar2 = r_incident.ppar2*2.*sqrt(rapport)*cosinus/
-		 	   (sqrt(rapport-sinus*sinus)+rapport*cosinus);
-	 r_transmis.pper2 = r_incident.pper2*2.*cosinus/(cosinus+sqrt(rapport
+	 r_transmis.ppar2 = r_incident.ppar2*2.*Sqrt(rapport)*cosinus/
+		 	   (Sqrt(rapport-sinus*sinus)+rapport*cosinus);
+	 r_transmis.pper2 = r_incident.pper2*2.*cosinus/(cosinus+Sqrt(rapport
 			   - sinus*sinus));
-	 r_transmis.e = (r_incident.e/2)*(sqrt(rapport-sinus*sinus)/cosinus)
+	 r_transmis.e = (r_incident.e/2)*(Sqrt(rapport-sinus*sinus)/cosinus)
 	    *(((r_transmis.ppar1*r_transmis.ppar1+r_transmis.pper1*
 	        r_transmis.pper1)
 	    /(r_incident.ppar1*r_incident.ppar1+r_incident.pper1*
@@ -506,22 +518,6 @@ TRAYON r_incident;
 
 
 
-static int
-compare(r1,r2,marge)
-TRAYON r1, r2;
-double marge;
-
-{
- double arctg1, arctg2;
-
- arctg1 = atan2(Y(r1),X(r1));
- arctg2 = atan2(Y(r2),X(r2));
- if ((arctg1 - marge <= arctg2) && (arctg1 + marge >= arctg2)) return 1;
- else return 0;
-}
-
-
-
 #define ensuite(rayon,prob_passage,destination) r_suite = rayon; \
 				 r_suite.e = prob_passage(rayon)*rayon.e; \
 				 r_suite.dest = destination; \
@@ -633,63 +629,8 @@ TRAYON *r1,*r2;
 
 
 
-static double
-l_get_val()
-
-{
- int val, dir, i, trouve, curseur;
- int nb;
- double valeur;
- TRAYON *rayt, raynull;
- 
- if (prismclock < 0 || prismclock < eclock) setprism();
- if (bidon == BADVAL) {
-	errno = EDOM;
-	return(0.0);
- }
- val = (int)(argument(1) + .5);
- dir = (int)(argument(2) + .5);
- nb = (int)(argument(3) + .5);
- X(raynull) = bidon;
- Y(raynull) = Z(raynull) = 0.;
- raynull.e = 0.;
- trouve = curseur = 0;
- if ( !nosource && nb==2 ) nb=1; /* on est en train de tracer la source
-				     a partir de sa seconde source virtuelle */ 
-#ifdef DEBUG
- fprintf(stderr, " On considere le rayon no: %d\n", nb);
-#endif
- for(i=0; i < nbrayons &&!trouve; i++)
-  {
-   if(ray[i].v[0] * dir * sens >= 0.) curseur ++;
-   if(curseur == nb)
-   {
-    rayt = &ray[i];
-    trouve = 1;
-   }
-  }
- if(!trouve) rayt = &raynull; 
- switch(val) {
-	case 0 : valeur = rayt->v[0];
-	         break;
-	case 1 : valeur = rayt->v[1];
-		 break;
-	case 2 : valeur = rayt->v[2];
-		 break;
-	case 3 : valeur = rayt->e;
-		 break;
-	default : errno = EDOM; return(0.0);
-    } 
-#ifdef DEBUG
-  fprintf(stderr, "get_val( %i, %i, %i) = %lf\n",val,dir,nb,valeur);
-#endif
-  return valeur;
-}
-
-
 static
 setprism()
-
 {
  double d;
  TRAYON r_initial,rsource;
@@ -784,6 +725,61 @@ if ( X(r_initial) != 0.)
  bidon = BADVAL;
  return;
 }
+
+
+static double
+l_get_val()
+
+{
+ int val, dir, i, trouve, curseur;
+ int nb;
+ double valeur;
+ TRAYON *rayt, raynull;
+ 
+ if (prismclock < 0 || prismclock < eclock) setprism();
+ if (bidon == BADVAL) {
+	errno = EDOM;
+	return(0.0);
+ }
+ val = (int)(argument(1) + .5);
+ dir = (int)(argument(2) + .5);
+ nb = (int)(argument(3) + .5);
+ X(raynull) = bidon;
+ Y(raynull) = Z(raynull) = 0.;
+ raynull.e = 0.;
+ trouve = curseur = 0;
+ if ( !nosource && nb==2 ) nb=1; /* on est en train de tracer la source
+				     a partir de sa seconde source virtuelle */ 
+#ifdef DEBUG
+ fprintf(stderr, " On considere le rayon no: %d\n", nb);
+#endif
+ for(i=0; i < nbrayons &&!trouve; i++)
+  {
+   if(ray[i].v[0] * dir * sens >= 0.) curseur ++;
+   if(curseur == nb)
+   {
+    rayt = &ray[i];
+    trouve = 1;
+   }
+  }
+ if(!trouve) rayt = &raynull; 
+ switch(val) {
+	case 0 : valeur = rayt->v[0];
+	         break;
+	case 1 : valeur = rayt->v[1];
+		 break;
+	case 2 : valeur = rayt->v[2];
+		 break;
+	case 3 : valeur = rayt->e;
+		 break;
+	default : errno = EDOM; return(0.0);
+    } 
+#ifdef DEBUG
+  fprintf(stderr, "get_val( %i, %i, %i) = %lf\n",val,dir,nb,valeur);
+#endif
+  return valeur;
+}
+
 
 setprismfuncs()
 {