src/rt/srcsupp.c

#ifndef lint
static const char       RCSid[] = "$Id: srcsupp.c,v 2.21 2014/12/04 05:26:28 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;

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