src/rt/srcsupp.c

#ifndef lint
static const char       RCSid[] = "$Id: srcsupp.c,v 2.19 2014/06/17 18:57:41 greg Exp $";
#endif
/*
 *  Support routines for source objects and materials
 *
 *  External symbols declared in source.h
 */

#include "copyright.h"

#include  "ray.h"

#include  "otypes.h"

#include  "source.h"

#include  "cone.h"

#include  "face.h"

#define SRCINC          32              /* realloc increment for array */

SRCREC  *source = NULL;                 /* our list of sources */
int  nsources = 0;                      /* the number of sources */

SRCFUNC  sfun[NUMOTYPE];                /* source dispatch table */


void
initstypes(void)                        /* initialize source dispatch table */
{
        extern VSMATERIAL  mirror_vs, direct1_vs, direct2_vs;
        static SOBJECT  fsobj = {fsetsrc, flatpart, fgetplaneq, fgetmaxdisk};
        static SOBJECT  ssobj = {ssetsrc, nopart};
        static SOBJECT  sphsobj = {sphsetsrc, nopart};
        static SOBJECT  cylsobj = {cylsetsrc, cylpart};
        static SOBJECT  rsobj = {rsetsrc, flatpart, rgetplaneq, rgetmaxdisk};

        sfun[MAT_MIRROR].mf = &mirror_vs;
        sfun[MAT_DIRECT1].mf = &direct1_vs;
        sfun[MAT_DIRECT2].mf = &direct2_vs;
        sfun[OBJ_FACE].of = &fsobj;
        sfun[OBJ_SOURCE].of = &ssobj;
        sfun[OBJ_SPHERE].of = &sphsobj;
        sfun[OBJ_CYLINDER].of = &cylsobj;
        sfun[OBJ_RING].of = &rsobj;
}


int
newsource(void)                 /* allocate new source in our array */
{
        if (nsources == 0)
                source = (SRCREC *)malloc(SRCINC*sizeof(SRCREC));
        else if (nsources%SRCINC == 0)
                source = (SRCREC *)realloc((void *)source,
                                (unsigned)(nsources+SRCINC)*sizeof(SRCREC));
        if (source == NULL)
                return(-1);
        source[nsources].sflags = 0;
        source[nsources].nhits = 1;
        source[nsources].ntests = 2;    /* initial hit probability = 50% */
#if SHADCACHE
        source[nsources].obscache = NULL;
#endif
        return(nsources++);
}


void
setflatss(                              /* set sampling for a flat source */
        SRCREC  *src
)
{
        double  mult;
        int  i;

        src->ss[SV][0] = src->ss[SV][1] = src->ss[SV][2] = 0.0;
        for (i = 0; i < 3; i++)
                if (src->snorm[i] < 0.6 && src->snorm[i] > -0.6)
                        break;
        src->ss[SV][i] = 1.0;
        fcross(src->ss[SU], src->ss[SV], src->snorm);
        mult = .5 * sqrt( src->ss2 / DOT(src->ss[SU],src->ss[SU]) );
        for (i = 0; i < 3; i++)
                src->ss[SU][i] *= mult;
        fcross(src->ss[SV], src->snorm, src->ss[SU]);
}


void
fsetsrc(                        /* set a face as a source */
        SRCREC  *src,
        OBJREC  *so
)
{
        FACE  *f;
        int  i, j;
        double  d;
        
        src->sa.success = 2*AIMREQT-1;          /* bitch on second failure */
        src->so = so;
                                                /* get the face */
        f = getface(so);
        if (f->area == 0.0)
                objerror(so, USER, "zero source area");
                                                /* find the center */
        for (j = 0; j < 3; j++) {
                src->sloc[j] = 0.0;
                for (i = 0; i < f->nv; i++)
                        src->sloc[j] += VERTEX(f,i)[j];
                src->sloc[j] /= (double)f->nv;
        }
        if (!inface(src->sloc, f))
                objerror(so, USER, "cannot hit source center");
        src->sflags |= SFLAT;
        VCOPY(src->snorm, f->norm);
        src->ss2 = f->area;
                                                /* find maximum radius */
        src->srad = 0.;
        for (i = 0; i < f->nv; i++) {
                d = dist2(VERTEX(f,i), src->sloc);
                if (d > src->srad)
                        src->srad = d;
        }
        src->srad = sqrt(src->srad);
                                                /* compute size vectors */
        if (f->nv == 4)                         /* parallelogram case */
                for (j = 0; j < 3; j++) {
                        src->ss[SU][j] = .5*(VERTEX(f,1)[j]-VERTEX(f,0)[j]);
                        src->ss[SV][j] = .5*(VERTEX(f,3)[j]-VERTEX(f,0)[j]);
                }
        else
                setflatss(src);
}


void
ssetsrc(                        /* set a source as a source */
        SRCREC  *src,
        OBJREC  *so
)
{
        double  theta;
        
        src->sa.success = 2*AIMREQT-1;          /* bitch on second failure */
        src->so = so;
        if (so->oargs.nfargs != 4)
                objerror(so, USER, "bad arguments");
        src->sflags |= (SDISTANT|SCIR);
        VCOPY(src->sloc, so->oargs.farg);
        if (normalize(src->sloc) == 0.0)
                objerror(so, USER, "zero direction");
        theta = PI/180.0/2.0 * so->oargs.farg[3];
        if (theta <= FTINY)
                objerror(so, USER, "zero size");
        src->ss2 = 2.0*PI * (1.0 - cos(theta));
                                        /* the following is approximate */
        src->srad = sqrt(src->ss2/PI);
        VCOPY(src->snorm, src->sloc);
        setflatss(src);                 /* hey, whatever works */
        src->ss[SW][0] = src->ss[SW][1] = src->ss[SW][2] = 0.0;
}


void
sphsetsrc(                      /* set a sphere as a source */
        SRCREC  *src,
        OBJREC  *so
)
{
        int  i;

        src->sa.success = 2*AIMREQT-1;          /* bitch on second failure */
        src->so = so;
        if (so->oargs.nfargs != 4)
                objerror(so, USER, "bad # arguments");
        if (so->oargs.farg[3] <= FTINY)
                objerror(so, USER, "illegal source radius");
        src->sflags |= SCIR;
        VCOPY(src->sloc, so->oargs.farg);
        src->srad = so->oargs.farg[3];
        src->ss2 = PI * src->srad * src->srad;
        for (i = 0; i < 3; i++)
                src->ss[SU][i] = src->ss[SV][i] = src->ss[SW][i] = 0.0;
        for (i = 0; i < 3; i++)
                src->ss[i][i] = 0.7236 * so->oargs.farg[3];
}


void
rsetsrc(                        /* set a ring (disk) as a source */
        SRCREC  *src,
        OBJREC  *so
)
{
        CONE  *co;
        
        src->sa.success = 2*AIMREQT-1;          /* bitch on second failure */
        src->so = so;
                                                /* get the ring */
        co = getcone(so, 0);
        if (CO_R1(co) <= FTINY)
                objerror(so, USER, "illegal source radius");
        VCOPY(src->sloc, CO_P0(co));
        if (CO_R0(co) > 0.0)
                objerror(so, USER, "cannot hit source center");
        src->sflags |= (SFLAT|SCIR);
        VCOPY(src->snorm, co->ad);
        src->srad = CO_R1(co);
        src->ss2 = PI * src->srad * src->srad;
        setflatss(src);
}


void
cylsetsrc(                      /* set a cylinder as a source */
        SRCREC  *src,
        OBJREC  *so
)
{
        CONE  *co;
        int  i;
        
        src->sa.success = 4*AIMREQT-1;          /* bitch on fourth failure */
        src->so = so;
                                                /* get the cylinder */
        co = getcone(so, 0);
        if (CO_R0(co) <= FTINY)
                objerror(so, USER, "illegal source radius");
        if (CO_R0(co) > .2*co->al)              /* heuristic constraint */
                objerror(so, WARNING, "source aspect too small");
        src->sflags |= SCYL;
        for (i = 0; i < 3; i++)
                src->sloc[i] = .5 * (CO_P1(co)[i] + CO_P0(co)[i]);
        src->srad = .5*co->al;
        src->ss2 = 2.*CO_R0(co)*co->al;
                                                /* set sampling vectors */
        for (i = 0; i < 3; i++)
                src->ss[SU][i] = .5 * co->al * co->ad[i];
        src->ss[SV][0] = src->ss[SV][1] = src->ss[SV][2] = 0.0;
        for (i = 0; i < 3; i++)
                if (co->ad[i] < 0.6 && co->ad[i] > -0.6)
                        break;
        src->ss[SV][i] = 1.0;
        fcross(src->ss[SW], src->ss[SV], co->ad);
        normalize(src->ss[SW]);
        for (i = 0; i < 3; i++)
                src->ss[SW][i] *= .8559 * CO_R0(co);
        fcross(src->ss[SV], src->ss[SW], co->ad);
}


SPOT *
makespot(                       /* make a spotlight */
        OBJREC  *m
)
{
        SPOT  *ns;

        if ((ns = (SPOT *)m->os) != NULL)
                return(ns);
        if ((ns = (SPOT *)malloc(sizeof(SPOT))) == NULL)
                return(NULL);
        if (m->oargs.farg[3] <= FTINY)
                objerror(m, USER, "zero angle");
        ns->siz = 2.0*PI * (1.0 - cos(PI/180.0/2.0 * m->oargs.farg[3]));
        VCOPY(ns->aim, m->oargs.farg+4);
        if ((ns->flen = normalize(ns->aim)) == 0.0)
                objerror(m, USER, "zero focus vector");
        m->os = (char *)ns;
        return(ns);
}


int
spotout(                        /* check if we're outside spot region */
        RAY  *r,
        SPOT  *s
)
{
        double  d;
        FVECT  vd;
        
        if (s == NULL)
                return(0);
        if (s->flen < -FTINY) {         /* distant source */
                vd[0] = s->aim[0] - r->rorg[0];
                vd[1] = s->aim[1] - r->rorg[1];
                vd[2] = s->aim[2] - r->rorg[2];
                d = DOT(r->rdir,vd);
                /*                      wrong side?
                if (d <= FTINY)
                        return(1);      */
                d = DOT(vd,vd) - d*d;
                if (PI*d > s->siz)
                        return(1);      /* out */
                return(0);      /* OK */
        }
                                        /* local source */
        if (s->siz < 2.0*PI * (1.0 + DOT(s->aim,r->rdir)))
                return(1);      /* out */
        return(0);      /* OK */
}


double
fgetmaxdisk(            /* get center and squared radius of face */
        FVECT  ocent,
        OBJREC  *op
)
{
        double  maxrad2;
        double  d;
        int  i, j;
        FACE  *f;
        
        f = getface(op);
        if (f->area == 0.)
                return(0.);
        for (i = 0; i < 3; i++) {
                ocent[i] = 0.;
                for (j = 0; j < f->nv; j++)
                        ocent[i] += VERTEX(f,j)[i];
                ocent[i] /= (double)f->nv;
        }
        d = DOT(ocent,f->norm);
        for (i = 0; i < 3; i++)
                ocent[i] += (f->offset - d)*f->norm[i];
        maxrad2 = 0.;
        for (j = 0; j < f->nv; j++) {
                d = dist2(VERTEX(f,j), ocent);
                if (d > maxrad2)
                        maxrad2 = d;
        }
        return(maxrad2);
}


double
rgetmaxdisk(            /* get center and squared radius of ring */
        FVECT  ocent,
        OBJREC  *op
)
{
        CONE  *co;
        
        co = getcone(op, 0);
        VCOPY(ocent, CO_P0(co));
        return(CO_R1(co)*CO_R1(co));
}


double
fgetplaneq(                     /* get plane equation for face */
        FVECT  nvec,
        OBJREC  *op
)
{
        FACE  *fo;

        fo = getface(op);
        VCOPY(nvec, fo->norm);
        return(fo->offset);
}


double
rgetplaneq(                     /* get plane equation for ring */
        FVECT  nvec,
        OBJREC  *op
)
{
        CONE  *co;

        co = getcone(op, 0);
        VCOPY(nvec, co->ad);
        return(DOT(nvec, CO_P0(co)));
}


int
commonspot(             /* set sp1 to intersection of sp1 and sp2 */
        SPOT  *sp1,
        SPOT  *sp2,
        FVECT  org
)
{
        FVECT  cent;
        double  rad2, cos1, cos2;

        cos1 = 1. - sp1->siz/(2.*PI);
        cos2 = 1. - sp2->siz/(2.*PI);
        if (sp2->siz >= 2.*PI-FTINY)            /* BIG, just check overlap */
                return(DOT(sp1->aim,sp2->aim) >= cos1*cos2 -
                                        sqrt((1.-cos1*cos1)*(1.-cos2*cos2)));
                                /* compute and check disks */
        rad2 = intercircle(cent, sp1->aim, sp2->aim,
                        1./(cos1*cos1) - 1.,  1./(cos2*cos2) - 1.);
        if (rad2 <= FTINY || normalize(cent) == 0.)
                return(0);
        VCOPY(sp1->aim, cent);
        sp1->siz = 2.*PI*(1. - 1./sqrt(1.+rad2));
        return(1);
}


int
commonbeam(             /* set sp1 to intersection of sp1 and sp2 */
        SPOT  *sp1,
        SPOT  *sp2,
        FVECT  dir
)
{
        FVECT  cent, c1, c2;
        double  rad2, d;
                                        /* move centers to common plane */
        d = DOT(sp1->aim, dir);
        VSUM(c1, sp1->aim, dir, -d);
        d = DOT(sp2->aim, dir);
        VSUM(c2, sp2->aim, dir, -d);
                                        /* compute overlap */
        rad2 = intercircle(cent, c1, c2, sp1->siz/PI, sp2->siz/PI);
        if (rad2 <= FTINY)
                return(0);
        VCOPY(sp1->aim, cent);
        sp1->siz = PI*rad2;
        return(1);
}


int
checkspot(                      /* check spotlight for behind source */
        SPOT  *sp,      /* spotlight */
        FVECT  nrm      /* source surface normal */
)
{
        double  d, d1;

        d = DOT(sp->aim, nrm);
        if (d > FTINY)                  /* center in front? */
                return(1);
                                        /* else check horizon */
        d1 = 1. - sp->siz/(2.*PI);
        return(1.-FTINY-d*d < d1*d1);
}


double
spotdisk(               /* intersect spot with object op */
        FVECT  oc,
        OBJREC  *op,
        SPOT  *sp,
        FVECT  pos
)
{
        FVECT  onorm;
        double  offs, d, dist;

        offs = getplaneq(onorm, op);
        d = -DOT(onorm, sp->aim);
        if (d >= -FTINY && d <= FTINY)
                return(0.);
        dist = (DOT(pos, onorm) - offs)/d;
        if (dist < 0.)
                return(0.);
        VSUM(oc, pos, sp->aim, dist);
        return(sp->siz*dist*dist/PI/(d*d));
}


double
beamdisk(               /* intersect beam with object op */
        FVECT  oc,
        OBJREC  *op,
        SPOT  *sp,
        FVECT  dir
)
{
        FVECT  onorm;
        double  offs, d, dist;

        offs = getplaneq(onorm, op);
        d = -DOT(onorm, dir);
        if (d >= -FTINY && d <= FTINY)
                return(0.);
        dist = (DOT(sp->aim, onorm) - offs)/d;
        VSUM(oc, sp->aim, dir, dist);
        return(sp->siz/PI/(d*d));
}


double
intercircle(                            /* intersect two circles */
        FVECT  cc,              /* midpoint (return value) */
        FVECT  c1,              /* circle centers */
        FVECT  c2,
        double  r1s,            /* radii squared */
        double  r2s
)
{
        double  a2, d2, l;
        FVECT  disp;

        VSUB(disp, c2, c1);
        d2 = DOT(disp,disp);
                                        /* circle within overlap? */
        if (r1s < r2s) {
                if (r2s >= r1s + d2) {
                        VCOPY(cc, c1);
                        return(r1s);
                }
        } else {
                if (r1s >= r2s + d2) {
                        VCOPY(cc, c2);
                        return(r2s);
                }
        }
        a2 = .25*(2.*(r1s+r2s) - d2 - (r2s-r1s)*(r2s-r1s)/d2);
                                        /* no overlap? */
        if (a2 <= 0.)
                return(0.);
                                        /* overlap, compute center */
        l = sqrt((r1s - a2)/d2);
        VSUM(cc, c1, disp, l);
        return(a2);
}
Revision:	2.20
Committed:	Sun Jun 22 18:05:28 2014 UTC (9 years, 10 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad4R2P2, rad4R2, rad4R2P1
Changes since 2.19:	+95 -73 lines
Log Message:	ANSIfication
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.20	static const char RCSid[] = "$Id: srcsupp.c,v 2.19 2014/06/17 18:57:41 greg Exp $";
3	greg	1.1	#endif
4			/*
5			* Support routines for source objects and materials
6	greg	2.9	*
7			* External symbols declared in source.h
8			*/
9
10	greg	2.10	#include "copyright.h"
11	greg	1.1
12			#include "ray.h"
13
14			#include "otypes.h"
15
16			#include "source.h"
17
18			#include "cone.h"
19
20			#include "face.h"
21
22	greg	2.19	#define SRCINC 32 /* realloc increment for array */
23	greg	1.1
24			SRCREC source = NULL; / our list of sources */
25			int nsources = 0; /* the number of sources */
26
27			SRCFUNC sfun[NUMOTYPE]; /* source dispatch table */
28
29
30	greg	2.9	void
31	greg	2.20	initstypes(void) /* initialize source dispatch table */
32	greg	1.1	{
33	greg	1.9	extern VSMATERIAL mirror_vs, direct1_vs, direct2_vs;
34	greg	1.14	static SOBJECT fsobj = {fsetsrc, flatpart, fgetplaneq, fgetmaxdisk};
35			static SOBJECT ssobj = {ssetsrc, nopart};
36			static SOBJECT sphsobj = {sphsetsrc, nopart};
37			static SOBJECT cylsobj = {cylsetsrc, cylpart};
38			static SOBJECT rsobj = {rsetsrc, flatpart, rgetplaneq, rgetmaxdisk};
39	greg	1.1
40			sfun[MAT_MIRROR].mf = &mirror_vs;
41	greg	1.9	sfun[MAT_DIRECT1].mf = &direct1_vs;
42			sfun[MAT_DIRECT2].mf = &direct2_vs;
43	greg	1.1	sfun[OBJ_FACE].of = &fsobj;
44			sfun[OBJ_SOURCE].of = &ssobj;
45			sfun[OBJ_SPHERE].of = &sphsobj;
46	greg	1.14	sfun[OBJ_CYLINDER].of = &cylsobj;
47	greg	1.1	sfun[OBJ_RING].of = &rsobj;
48			}
49
50
51	greg	1.2	int
52	greg	2.20	newsource(void) /* allocate new source in our array */
53	greg	1.1	{
54			if (nsources == 0)
55	greg	1.13	source = (SRCREC )malloc(SRCINCsizeof(SRCREC));
56			else if (nsources%SRCINC == 0)
57	greg	2.11	source = (SRCREC )realloc((void )source,
58	greg	1.13	(unsigned)(nsources+SRCINC)*sizeof(SRCREC));
59	greg	1.1	if (source == NULL)
60	greg	1.2	return(-1);
61	greg	1.1	source[nsources].sflags = 0;
62	greg	2.14	source[nsources].nhits = 1;
63			source[nsources].ntests = 2; /* initial hit probability = 50% */
64	greg	2.12	#if SHADCACHE
65			source[nsources].obscache = NULL;
66			#endif
67	greg	1.2	return(nsources++);
68	greg	1.1	}
69
70
71	greg	2.9	void
72	greg	2.20	setflatss( /* set sampling for a flat source */
73			SRCREC *src
74			)
75	greg	1.14	{
76			double mult;
77	greg	2.20	int i;
78	greg	1.14
79			src->ss[SV][0] = src->ss[SV][1] = src->ss[SV][2] = 0.0;
80			for (i = 0; i < 3; i++)
81			if (src->snorm[i] < 0.6 && src->snorm[i] > -0.6)
82			break;
83			src->ss[SV][i] = 1.0;
84			fcross(src->ss[SU], src->ss[SV], src->snorm);
85			mult = .5 * sqrt( src->ss2 / DOT(src->ss[SU],src->ss[SU]) );
86			for (i = 0; i < 3; i++)
87			src->ss[SU][i] *= mult;
88			fcross(src->ss[SV], src->snorm, src->ss[SU]);
89			}
90
91
92	greg	2.9	void
93	greg	2.20	fsetsrc( /* set a face as a source */
94			SRCREC *src,
95			OBJREC *so
96			)
97	greg	1.1	{
98	greg	2.20	FACE *f;
99			int i, j;
100	greg	1.14	double d;
101	greg	1.1
102			src->sa.success = 2AIMREQT-1; / bitch on second failure */
103			src->so = so;
104			/* get the face */
105			f = getface(so);
106	greg	2.16	if (f->area == 0.0)
107			objerror(so, USER, "zero source area");
108	greg	1.1	/* find the center */
109			for (j = 0; j < 3; j++) {
110			src->sloc[j] = 0.0;
111			for (i = 0; i < f->nv; i++)
112			src->sloc[j] += VERTEX(f,i)[j];
113			src->sloc[j] /= (double)f->nv;
114			}
115			if (!inface(src->sloc, f))
116	greg	2.16	objerror(so, USER, "cannot hit source center");
117	greg	1.1	src->sflags \|= SFLAT;
118			VCOPY(src->snorm, f->norm);
119			src->ss2 = f->area;
120	greg	1.14	/* find maximum radius */
121			src->srad = 0.;
122			for (i = 0; i < f->nv; i++) {
123			d = dist2(VERTEX(f,i), src->sloc);
124			if (d > src->srad)
125			src->srad = d;
126			}
127			src->srad = sqrt(src->srad);
128			/* compute size vectors */
129	greg	2.7	if (f->nv == 4) /* parallelogram case */
130	greg	1.14	for (j = 0; j < 3; j++) {
131			src->ss[SU][j] = .5*(VERTEX(f,1)[j]-VERTEX(f,0)[j]);
132			src->ss[SV][j] = .5*(VERTEX(f,3)[j]-VERTEX(f,0)[j]);
133			}
134			else
135			setflatss(src);
136	greg	1.1	}
137
138
139	greg	2.9	void
140	greg	2.20	ssetsrc( /* set a source as a source */
141			SRCREC *src,
142			OBJREC *so
143			)
144	greg	1.1	{
145			double theta;
146
147			src->sa.success = 2AIMREQT-1; / bitch on second failure */
148			src->so = so;
149			if (so->oargs.nfargs != 4)
150			objerror(so, USER, "bad arguments");
151	greg	2.17	src->sflags \|= (SDISTANT\|SCIR);
152	greg	1.1	VCOPY(src->sloc, so->oargs.farg);
153			if (normalize(src->sloc) == 0.0)
154			objerror(so, USER, "zero direction");
155			theta = PI/180.0/2.0 * so->oargs.farg[3];
156			if (theta <= FTINY)
157			objerror(so, USER, "zero size");
158			src->ss2 = 2.0PI (1.0 - cos(theta));
159	greg	1.14	/* the following is approximate */
160			src->srad = sqrt(src->ss2/PI);
161			VCOPY(src->snorm, src->sloc);
162			setflatss(src); /* hey, whatever works */
163			src->ss[SW][0] = src->ss[SW][1] = src->ss[SW][2] = 0.0;
164	greg	1.1	}
165
166
167	greg	2.9	void
168	greg	2.20	sphsetsrc( /* set a sphere as a source */
169			SRCREC *src,
170			OBJREC *so
171			)
172	greg	1.1	{
173	greg	2.20	int i;
174	greg	1.14
175	greg	1.1	src->sa.success = 2AIMREQT-1; / bitch on second failure */
176			src->so = so;
177			if (so->oargs.nfargs != 4)
178			objerror(so, USER, "bad # arguments");
179			if (so->oargs.farg[3] <= FTINY)
180	greg	2.16	objerror(so, USER, "illegal source radius");
181	greg	2.17	src->sflags \|= SCIR;
182	greg	1.1	VCOPY(src->sloc, so->oargs.farg);
183	greg	1.14	src->srad = so->oargs.farg[3];
184			src->ss2 = PI * src->srad * src->srad;
185			for (i = 0; i < 3; i++)
186			src->ss[SU][i] = src->ss[SV][i] = src->ss[SW][i] = 0.0;
187			for (i = 0; i < 3; i++)
188	greg	2.17	src->ss[i][i] = 0.7236 * so->oargs.farg[3];
189	greg	1.1	}
190
191
192	greg	2.9	void
193	greg	2.20	rsetsrc( /* set a ring (disk) as a source */
194			SRCREC *src,
195			OBJREC *so
196			)
197	greg	1.1	{
198	greg	2.20	CONE *co;
199	greg	1.1
200			src->sa.success = 2AIMREQT-1; / bitch on second failure */
201			src->so = so;
202			/* get the ring */
203			co = getcone(so, 0);
204	greg	2.16	if (CO_R1(co) <= FTINY)
205			objerror(so, USER, "illegal source radius");
206	greg	1.1	VCOPY(src->sloc, CO_P0(co));
207			if (CO_R0(co) > 0.0)
208	greg	2.16	objerror(so, USER, "cannot hit source center");
209	greg	2.17	src->sflags \|= (SFLAT\|SCIR);
210	greg	1.1	VCOPY(src->snorm, co->ad);
211	greg	1.14	src->srad = CO_R1(co);
212			src->ss2 = PI * src->srad * src->srad;
213			setflatss(src);
214			}
215
216
217	greg	2.9	void
218	greg	2.20	cylsetsrc( /* set a cylinder as a source */
219			SRCREC *src,
220			OBJREC *so
221			)
222	greg	1.14	{
223	greg	2.20	CONE *co;
224			int i;
225	greg	1.14
226			src->sa.success = 4AIMREQT-1; / bitch on fourth failure */
227			src->so = so;
228			/* get the cylinder */
229			co = getcone(so, 0);
230	greg	2.16	if (CO_R0(co) <= FTINY)
231			objerror(so, USER, "illegal source radius");
232	greg	1.14	if (CO_R0(co) > .2co->al) / heuristic constraint */
233			objerror(so, WARNING, "source aspect too small");
234	greg	1.15	src->sflags \|= SCYL;
235	greg	1.14	for (i = 0; i < 3; i++)
236			src->sloc[i] = .5 * (CO_P1(co)[i] + CO_P0(co)[i]);
237	greg	1.15	src->srad = .5*co->al;
238	greg	1.14	src->ss2 = 2.CO_R0(co)co->al;
239			/* set sampling vectors */
240			for (i = 0; i < 3; i++)
241			src->ss[SU][i] = .5 * co->al * co->ad[i];
242			src->ss[SV][0] = src->ss[SV][1] = src->ss[SV][2] = 0.0;
243			for (i = 0; i < 3; i++)
244			if (co->ad[i] < 0.6 && co->ad[i] > -0.6)
245			break;
246			src->ss[SV][i] = 1.0;
247			fcross(src->ss[SW], src->ss[SV], co->ad);
248			normalize(src->ss[SW]);
249			for (i = 0; i < 3; i++)
250	greg	1.15	src->ss[SW][i] = .8559 CO_R0(co);
251	greg	1.14	fcross(src->ss[SV], src->ss[SW], co->ad);
252	greg	1.1	}
253
254
255			SPOT *
256	greg	2.20	makespot( /* make a spotlight */
257			OBJREC *m
258			)
259	greg	1.1	{
260	greg	2.20	SPOT *ns;
261	greg	1.1
262	greg	2.5	if ((ns = (SPOT *)m->os) != NULL)
263			return(ns);
264	greg	1.1	if ((ns = (SPOT *)malloc(sizeof(SPOT))) == NULL)
265			return(NULL);
266	greg	2.16	if (m->oargs.farg[3] <= FTINY)
267			objerror(m, USER, "zero angle");
268	greg	1.1	ns->siz = 2.0PI (1.0 - cos(PI/180.0/2.0 * m->oargs.farg[3]));
269			VCOPY(ns->aim, m->oargs.farg+4);
270			if ((ns->flen = normalize(ns->aim)) == 0.0)
271			objerror(m, USER, "zero focus vector");
272	greg	2.5	m->os = (char *)ns;
273	greg	1.1	return(ns);
274			}
275
276
277	greg	2.9	int
278	greg	2.20	spotout( /* check if we're outside spot region */
279			RAY *r,
280			SPOT *s
281			)
282	greg	2.5	{
283			double d;
284			FVECT vd;
285
286			if (s == NULL)
287			return(0);
288	greg	2.8	if (s->flen < -FTINY) { /* distant source */
289	greg	2.5	vd[0] = s->aim[0] - r->rorg[0];
290			vd[1] = s->aim[1] - r->rorg[1];
291			vd[2] = s->aim[2] - r->rorg[2];
292			d = DOT(r->rdir,vd);
293			/* wrong side?
294			if (d <= FTINY)
295			return(1); */
296			d = DOT(vd,vd) - d*d;
297			if (PI*d > s->siz)
298			return(1); /* out */
299			return(0); /* OK */
300			}
301			/* local source */
302			if (s->siz < 2.0PI (1.0 + DOT(s->aim,r->rdir)))
303			return(1); /* out */
304			return(0); /* OK */
305			}
306
307
308	greg	1.1	double
309	greg	2.20	fgetmaxdisk( /* get center and squared radius of face */
310			FVECT ocent,
311			OBJREC *op
312			)
313	greg	1.1	{
314			double maxrad2;
315	greg	1.5	double d;
316	greg	2.20	int i, j;
317			FACE *f;
318	greg	1.1
319			f = getface(op);
320	greg	1.5	if (f->area == 0.)
321			return(0.);
322	greg	1.1	for (i = 0; i < 3; i++) {
323			ocent[i] = 0.;
324			for (j = 0; j < f->nv; j++)
325			ocent[i] += VERTEX(f,j)[i];
326			ocent[i] /= (double)f->nv;
327			}
328	greg	1.5	d = DOT(ocent,f->norm);
329			for (i = 0; i < 3; i++)
330			ocent[i] += (f->offset - d)*f->norm[i];
331	greg	1.1	maxrad2 = 0.;
332			for (j = 0; j < f->nv; j++) {
333	greg	1.5	d = dist2(VERTEX(f,j), ocent);
334			if (d > maxrad2)
335			maxrad2 = d;
336	greg	1.1	}
337			return(maxrad2);
338			}
339
340
341			double
342	greg	2.20	rgetmaxdisk( /* get center and squared radius of ring */
343			FVECT ocent,
344			OBJREC *op
345			)
346	greg	1.1	{
347	greg	2.20	CONE *co;
348	greg	1.1
349			co = getcone(op, 0);
350			VCOPY(ocent, CO_P0(co));
351			return(CO_R1(co)*CO_R1(co));
352			}
353
354
355			double
356	greg	2.20	fgetplaneq( /* get plane equation for face */
357			FVECT nvec,
358			OBJREC *op
359			)
360	greg	1.1	{
361	greg	2.20	FACE *fo;
362	greg	1.1
363			fo = getface(op);
364			VCOPY(nvec, fo->norm);
365			return(fo->offset);
366			}
367
368
369			double
370	greg	2.20	rgetplaneq( /* get plane equation for ring */
371			FVECT nvec,
372			OBJREC *op
373			)
374	greg	1.1	{
375	greg	2.20	CONE *co;
376	greg	1.1
377			co = getcone(op, 0);
378			VCOPY(nvec, co->ad);
379			return(DOT(nvec, CO_P0(co)));
380	greg	1.4	}
381
382
383	greg	2.9	int
384	greg	2.20	commonspot( /* set sp1 to intersection of sp1 and sp2 */
385			SPOT *sp1,
386			SPOT *sp2,
387			FVECT org
388			)
389	greg	1.4	{
390			FVECT cent;
391			double rad2, cos1, cos2;
392
393			cos1 = 1. - sp1->siz/(2.*PI);
394			cos2 = 1. - sp2->siz/(2.*PI);
395			if (sp2->siz >= 2.PI-FTINY) / BIG, just check overlap */
396			return(DOT(sp1->aim,sp2->aim) >= cos1*cos2 -
397			sqrt((1.-cos1cos1)(1.-cos2*cos2)));
398			/* compute and check disks */
399			rad2 = intercircle(cent, sp1->aim, sp2->aim,
400			1./(cos1cos1) - 1., 1./(cos2cos2) - 1.);
401			if (rad2 <= FTINY \|\| normalize(cent) == 0.)
402			return(0);
403			VCOPY(sp1->aim, cent);
404			sp1->siz = 2.PI(1. - 1./sqrt(1.+rad2));
405			return(1);
406			}
407
408
409	greg	2.9	int
410	greg	2.20	commonbeam( /* set sp1 to intersection of sp1 and sp2 */
411			SPOT *sp1,
412			SPOT *sp2,
413			FVECT dir
414			)
415	greg	1.4	{
416			FVECT cent, c1, c2;
417			double rad2, d;
418			/* move centers to common plane */
419			d = DOT(sp1->aim, dir);
420	greg	2.18	VSUM(c1, sp1->aim, dir, -d);
421	greg	1.4	d = DOT(sp2->aim, dir);
422	greg	2.18	VSUM(c2, sp2->aim, dir, -d);
423	greg	1.4	/* compute overlap */
424			rad2 = intercircle(cent, c1, c2, sp1->siz/PI, sp2->siz/PI);
425			if (rad2 <= FTINY)
426			return(0);
427			VCOPY(sp1->aim, cent);
428			sp1->siz = PI*rad2;
429			return(1);
430			}
431
432
433	greg	2.9	int
434	greg	2.20	checkspot( /* check spotlight for behind source */
435			SPOT sp, / spotlight */
436			FVECT nrm /* source surface normal */
437			)
438	greg	1.4	{
439			double d, d1;
440
441			d = DOT(sp->aim, nrm);
442			if (d > FTINY) /* center in front? */
443	greg	1.8	return(1);
444	greg	1.4	/* else check horizon */
445			d1 = 1. - sp->siz/(2.*PI);
446	greg	1.8	return(1.-FTINY-dd < d1d1);
447	greg	1.4	}
448
449
450			double
451	greg	2.20	spotdisk( /* intersect spot with object op */
452			FVECT oc,
453			OBJREC *op,
454			SPOT *sp,
455			FVECT pos
456			)
457	greg	1.6	{
458			FVECT onorm;
459			double offs, d, dist;
460
461			offs = getplaneq(onorm, op);
462			d = -DOT(onorm, sp->aim);
463			if (d >= -FTINY && d <= FTINY)
464			return(0.);
465			dist = (DOT(pos, onorm) - offs)/d;
466			if (dist < 0.)
467			return(0.);
468	greg	2.18	VSUM(oc, pos, sp->aim, dist);
469	greg	1.6	return(sp->sizdistdist/PI/(d*d));
470			}
471
472
473			double
474	greg	2.20	beamdisk( /* intersect beam with object op */
475			FVECT oc,
476			OBJREC *op,
477			SPOT *sp,
478			FVECT dir
479			)
480	greg	1.6	{
481			FVECT onorm;
482			double offs, d, dist;
483
484			offs = getplaneq(onorm, op);
485			d = -DOT(onorm, dir);
486			if (d >= -FTINY && d <= FTINY)
487			return(0.);
488			dist = (DOT(sp->aim, onorm) - offs)/d;
489	greg	2.18	VSUM(oc, sp->aim, dir, dist);
490	greg	1.6	return(sp->siz/PI/(d*d));
491			}
492
493
494			double
495	greg	2.20	intercircle( /* intersect two circles */
496			FVECT cc, /* midpoint (return value) */
497			FVECT c1, /* circle centers */
498			FVECT c2,
499			double r1s, /* radii squared */
500			double r2s
501			)
502	greg	1.4	{
503			double a2, d2, l;
504			FVECT disp;
505
506	greg	2.18	VSUB(disp, c2, c1);
507	greg	1.4	d2 = DOT(disp,disp);
508			/* circle within overlap? */
509			if (r1s < r2s) {
510			if (r2s >= r1s + d2) {
511			VCOPY(cc, c1);
512			return(r1s);
513			}
514			} else {
515			if (r1s >= r2s + d2) {
516			VCOPY(cc, c2);
517			return(r2s);
518			}
519			}
520			a2 = .25(2.(r1s+r2s) - d2 - (r2s-r1s)*(r2s-r1s)/d2);
521			/* no overlap? */
522			if (a2 <= 0.)
523			return(0.);
524			/* overlap, compute center */
525			l = sqrt((r1s - a2)/d2);
526	greg	2.18	VSUM(cc, c1, disp, l);
527	greg	1.4	return(a2);
528	greg	1.1	}