home *** CD-ROM | disk | FTP | other *** search
/ DP Tool Club 8 / CDASC08.ISO / NEWS / RADIANCE / SRC / RT / SPHERE.C < prev    next >
C/C++ Source or Header  |  1993-10-07  |  2KB  |  86 lines
  1. /* Copyright (c) 1992 Regents of the University of California */
  2.  
  3. #ifndef lint
  4. static char SCCSid[] = "@(#)sphere.c 2.2 10/24/92 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 FLOAT  *ap;
  27.     register int  i;
  28.  
  29.     if (so->oargs.nfargs != 4)
  30.         objerror(so, USER, "bad # arguments");
  31.     ap = so->oargs.farg;
  32.     if (ap[3] < -FTINY) {
  33.         objerror(so, WARNING, "negative radius");
  34.         so->otype = so->otype == OBJ_SPHERE ?
  35.                 OBJ_BUBBLE : OBJ_SPHERE;
  36.         ap[3] = -ap[3];
  37.     } else if (ap[3] <= FTINY)
  38.         objerror(so, USER, "zero radius");
  39.  
  40.     /*
  41.      *    We compute the intersection by substituting into
  42.      *  the surface equation for the sphere.  The resulting
  43.      *  quadratic equation in t is then solved for the
  44.      *  smallest positive root, which is our point of
  45.      *  intersection.
  46.      *    Since the ray is normalized, a should always be
  47.      *  one.  We compute it here to prevent instability in the
  48.      *  intersection calculation.
  49.      */
  50.                 /* compute quadratic coefficients */
  51.     a = b = c = 0.0;
  52.     for (i = 0; i < 3; i++) {
  53.         a += r->rdir[i]*r->rdir[i];
  54.         t = r->rorg[i] - ap[i];
  55.         b += 2.0*r->rdir[i]*t;
  56.         c += t*t;
  57.     }
  58.     c -= ap[3] * ap[3];
  59.  
  60.     nroots = quadratic(root, a, b, c);    /* solve quadratic */
  61.     
  62.     for (i = 0; i < nroots; i++)        /* get smallest positive */
  63.         if ((t = root[i]) > FTINY)
  64.             break;
  65.     if (i >= nroots)
  66.         return(0);            /* no positive root */
  67.  
  68.     if (t >= r->rot)
  69.         return(0);            /* other is closer */
  70.  
  71.     r->ro = so;
  72.     r->rot = t;
  73.                     /* compute normal */
  74.     a = ap[3];
  75.     if (so->otype == OBJ_BUBBLE)
  76.         a = -a;            /* reverse */
  77.     for (i = 0; i < 3; i++) {
  78.         r->rop[i] = r->rorg[i] + r->rdir[i]*t;
  79.         r->ron[i] = (r->rop[i] - ap[i]) / a;
  80.     }
  81.     r->rod = -DOT(r->rdir, r->ron);
  82.     r->rox = NULL;
  83.  
  84.     return(1);            /* hit */
  85. }
  86.