Back to index

radiance  4R0+20100331
sphere.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char RCSid[] = "$Id: sphere.c,v 2.7 2004/03/30 16:13:01 schorsch Exp $";
00003 #endif
00004 /*
00005  *  sphere.c - compute ray intersection with spheres.
00006  */
00007 
00008 #include "copyright.h"
00009 
00010 #include  "ray.h"
00011 #include  "otypes.h"
00012 #include  "rtotypes.h"
00013 
00014 
00015 extern int
00016 o_sphere(                   /* compute intersection with sphere */
00017        OBJREC  *so,
00018        register RAY  *r
00019 )
00020 {
00021        double  a, b, c;     /* coefficients for quadratic equation */
00022        double  root[2];     /* quadratic roots */
00023        int  nroots;
00024        double  t;
00025        register RREAL  *ap;
00026        register int  i;
00027 
00028        if (so->oargs.nfargs != 4)
00029               objerror(so, USER, "bad # arguments");
00030        ap = so->oargs.farg;
00031        if (ap[3] < -FTINY) {
00032               objerror(so, WARNING, "negative radius");
00033               so->otype = so->otype == OBJ_SPHERE ?
00034                             OBJ_BUBBLE : OBJ_SPHERE;
00035               ap[3] = -ap[3];
00036        } else if (ap[3] <= FTINY)
00037               objerror(so, USER, "zero radius");
00038 
00039        /*
00040         *     We compute the intersection by substituting into
00041         *  the surface equation for the sphere.  The resulting
00042         *  quadratic equation in t is then solved for the
00043         *  smallest positive root, which is our point of
00044         *  intersection.
00045         *     Since the ray is normalized, a should always be
00046         *  one.  We compute it here to prevent instability in the
00047         *  intersection calculation.
00048         */
00049                             /* compute quadratic coefficients */
00050        a = b = c = 0.0;
00051        for (i = 0; i < 3; i++) {
00052               a += r->rdir[i]*r->rdir[i];
00053               t = r->rorg[i] - ap[i];
00054               b += 2.0*r->rdir[i]*t;
00055               c += t*t;
00056        }
00057        c -= ap[3] * ap[3];
00058 
00059        nroots = quadratic(root, a, b, c); /* solve quadratic */
00060        
00061        for (i = 0; i < nroots; i++)              /* get smallest positive */
00062               if ((t = root[i]) > FTINY)
00063                      break;
00064        if (i >= nroots)
00065               return(0);                  /* no positive root */
00066 
00067        if (t >= r->rot)
00068               return(0);                  /* other is closer */
00069 
00070        r->ro = so;
00071        r->rot = t;
00072                                    /* compute normal */
00073        a = ap[3];
00074        if (so->otype == OBJ_BUBBLE)
00075               a = -a;                     /* reverse */
00076        for (i = 0; i < 3; i++) {
00077               r->rop[i] = r->rorg[i] + r->rdir[i]*t;
00078               r->ron[i] = (r->rop[i] - ap[i]) / a;
00079        }
00080        r->rod = -DOT(r->rdir, r->ron);
00081        r->rox = NULL;
00082        r->pert[0] = r->pert[1] = r->pert[2] = 0.0;
00083        r->uv[0] = r->uv[1] = 0.0;
00084 
00085        return(1);                  /* hit */
00086 }