ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/sphere.c
Revision: 1.1
Committed: Thu Feb 2 10:41:42 1989 UTC (35 years, 3 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Initial revision

File Contents

# User Rev Content
1 greg 1.1 /* Copyright (c) 1986 Regents of the University of California */
2    
3     #ifndef lint
4     static char SCCSid[] = "$SunId$ LBL";
5     #endif
6    
7     /*
8     * sphere.c - compute ray intersection with spheres.
9     *
10     * 8/19/85
11     */
12    
13     #include "ray.h"
14    
15     #include "otypes.h"
16    
17    
18     o_sphere(so, r) /* compute intersection with sphere */
19     OBJREC *so;
20     register RAY *r;
21     {
22     double a, b, c; /* coefficients for quadratic equation */
23     double root[2]; /* quadratic roots */
24     int nroots;
25     double t;
26     register double *ap;
27     register int i;
28    
29     if (so->oargs.nfargs != 4 || so->oargs.farg[3] <= FTINY)
30     objerror(so, USER, "bad arguments");
31    
32     ap = so->oargs.farg;
33    
34     /*
35     * We compute the intersection by substituting into
36     * the surface equation for the sphere. The resulting
37     * quadratic equation in t is then solved for the
38     * smallest positive root, which is our point of
39     * intersection.
40     * Because the ray direction is normalized, a is always 1.
41     */
42    
43     a = 1.0; /* compute quadratic coefficients */
44     b = c = 0.0;
45     for (i = 0; i < 3; i++) {
46     t = r->rorg[i] - ap[i];
47     b += 2.0*r->rdir[i]*t;
48     c += t*t;
49     }
50     c -= ap[3] * ap[3];
51    
52     nroots = quadratic(root, a, b, c); /* solve quadratic */
53    
54     for (i = 0; i < nroots; i++) /* get smallest positive */
55     if ((t = root[i]) > FTINY)
56     break;
57     if (i >= nroots)
58     return(0); /* no positive root */
59    
60     if (t < r->rot) { /* found closer intersection */
61     r->ro = so;
62     r->rot = t;
63     /* compute normal */
64     a = ap[3];
65     if (so->otype == OBJ_BUBBLE)
66     a = -a; /* reverse */
67     for (i = 0; i < 3; i++) {
68     r->rop[i] = r->rorg[i] + r->rdir[i]*t;
69     r->ron[i] = (r->rop[i] - ap[i]) / a;
70     }
71     r->rod = -DOT(r->rdir, r->ron);
72     }
73     return(1);
74     }