src/rt/srcsupp.c

#ifndef lint
static const char       RCSid[] = "$Id: srcsupp.c,v 2.22 2015/05/21 05:54:54 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.)
                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 == NULL)
                objerror(so, USER, "illegal source");
        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 == NULL)
                objerror(so, USER, "illegal source");
        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);
        if (co == NULL)
                return(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);
        if (co == NULL) {
                memset(nvec, 0, sizeof(FVECT));
                return(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.23
Committed:	Thu Apr 21 00:40:35 2016 UTC (8 years ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R2, rad5R1, rad5R3
Changes since 2.22:	+12 -2 lines
Log Message:	Made zero cone and ring radii non-fatal (degenerate -> ignore)
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.23	static const char RCSid[] = "$Id: srcsupp.c,v 2.22 2015/05/21 05:54:54 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.23	if (f->area == 0.)
102	greg	2.16	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.23	if (co == NULL)
200			objerror(so, USER, "illegal source");
201	greg	2.16	if (CO_R1(co) <= FTINY)
202			objerror(so, USER, "illegal source radius");
203	greg	1.1	VCOPY(src->sloc, CO_P0(co));
204			if (CO_R0(co) > 0.0)
205	greg	2.16	objerror(so, USER, "cannot hit source center");
206	greg	2.17	src->sflags \|= (SFLAT\|SCIR);
207	greg	1.1	VCOPY(src->snorm, co->ad);
208	greg	1.14	src->srad = CO_R1(co);
209			src->ss2 = PI * src->srad * src->srad;
210			setflatss(src);
211			}
212
213
214	greg	2.9	void
215	greg	2.20	cylsetsrc( /* set a cylinder as a source */
216			SRCREC *src,
217			OBJREC *so
218			)
219	greg	1.14	{
220	greg	2.20	CONE *co;
221			int i;
222	greg	1.14
223			src->sa.success = 4AIMREQT-1; / bitch on fourth failure */
224			src->so = so;
225			/* get the cylinder */
226			co = getcone(so, 0);
227	greg	2.23	if (co == NULL)
228			objerror(so, USER, "illegal source");
229	greg	2.16	if (CO_R0(co) <= FTINY)
230			objerror(so, USER, "illegal source radius");
231	greg	1.14	if (CO_R0(co) > .2co->al) / heuristic constraint */
232			objerror(so, WARNING, "source aspect too small");
233	greg	1.15	src->sflags \|= SCYL;
234	greg	1.14	for (i = 0; i < 3; i++)
235			src->sloc[i] = .5 * (CO_P1(co)[i] + CO_P0(co)[i]);
236	greg	1.15	src->srad = .5*co->al;
237	greg	1.14	src->ss2 = 2.CO_R0(co)co->al;
238			/* set sampling vectors */
239			for (i = 0; i < 3; i++)
240			src->ss[SU][i] = .5 * co->al * co->ad[i];
241	greg	2.22	getperpendicular(src->ss[SW], co->ad, rand_samp);
242	greg	1.14	for (i = 0; i < 3; i++)
243	greg	1.15	src->ss[SW][i] = .8559 CO_R0(co);
244	greg	1.14	fcross(src->ss[SV], src->ss[SW], co->ad);
245	greg	1.1	}
246
247
248			SPOT *
249	greg	2.20	makespot( /* make a spotlight */
250			OBJREC *m
251			)
252	greg	1.1	{
253	greg	2.20	SPOT *ns;
254	greg	1.1
255	greg	2.5	if ((ns = (SPOT *)m->os) != NULL)
256			return(ns);
257	greg	1.1	if ((ns = (SPOT *)malloc(sizeof(SPOT))) == NULL)
258			return(NULL);
259	greg	2.16	if (m->oargs.farg[3] <= FTINY)
260			objerror(m, USER, "zero angle");
261	greg	1.1	ns->siz = 2.0PI (1.0 - cos(PI/180.0/2.0 * m->oargs.farg[3]));
262			VCOPY(ns->aim, m->oargs.farg+4);
263			if ((ns->flen = normalize(ns->aim)) == 0.0)
264			objerror(m, USER, "zero focus vector");
265	greg	2.5	m->os = (char *)ns;
266	greg	1.1	return(ns);
267			}
268
269
270	greg	2.9	int
271	greg	2.20	spotout( /* check if we're outside spot region */
272			RAY *r,
273			SPOT *s
274			)
275	greg	2.5	{
276			double d;
277			FVECT vd;
278
279			if (s == NULL)
280			return(0);
281	greg	2.8	if (s->flen < -FTINY) { /* distant source */
282	greg	2.5	vd[0] = s->aim[0] - r->rorg[0];
283			vd[1] = s->aim[1] - r->rorg[1];
284			vd[2] = s->aim[2] - r->rorg[2];
285			d = DOT(r->rdir,vd);
286			/* wrong side?
287			if (d <= FTINY)
288			return(1); */
289			d = DOT(vd,vd) - d*d;
290			if (PI*d > s->siz)
291			return(1); /* out */
292			return(0); /* OK */
293			}
294			/* local source */
295			if (s->siz < 2.0PI (1.0 + DOT(s->aim,r->rdir)))
296			return(1); /* out */
297			return(0); /* OK */
298			}
299
300
301	greg	1.1	double
302	greg	2.20	fgetmaxdisk( /* get center and squared radius of face */
303			FVECT ocent,
304			OBJREC *op
305			)
306	greg	1.1	{
307			double maxrad2;
308	greg	1.5	double d;
309	greg	2.20	int i, j;
310			FACE *f;
311	greg	1.1
312			f = getface(op);
313	greg	1.5	if (f->area == 0.)
314			return(0.);
315	greg	1.1	for (i = 0; i < 3; i++) {
316			ocent[i] = 0.;
317			for (j = 0; j < f->nv; j++)
318			ocent[i] += VERTEX(f,j)[i];
319			ocent[i] /= (double)f->nv;
320			}
321	greg	1.5	d = DOT(ocent,f->norm);
322			for (i = 0; i < 3; i++)
323			ocent[i] += (f->offset - d)*f->norm[i];
324	greg	1.1	maxrad2 = 0.;
325			for (j = 0; j < f->nv; j++) {
326	greg	1.5	d = dist2(VERTEX(f,j), ocent);
327			if (d > maxrad2)
328			maxrad2 = d;
329	greg	1.1	}
330			return(maxrad2);
331			}
332
333
334			double
335	greg	2.20	rgetmaxdisk( /* get center and squared radius of ring */
336			FVECT ocent,
337			OBJREC *op
338			)
339	greg	1.1	{
340	greg	2.20	CONE *co;
341	greg	1.1
342			co = getcone(op, 0);
343	greg	2.23	if (co == NULL)
344			return(0.);
345	greg	1.1	VCOPY(ocent, CO_P0(co));
346			return(CO_R1(co)*CO_R1(co));
347			}
348
349
350			double
351	greg	2.20	fgetplaneq( /* get plane equation for face */
352			FVECT nvec,
353			OBJREC *op
354			)
355	greg	1.1	{
356	greg	2.20	FACE *fo;
357	greg	1.1
358			fo = getface(op);
359			VCOPY(nvec, fo->norm);
360			return(fo->offset);
361			}
362
363
364			double
365	greg	2.20	rgetplaneq( /* get plane equation for ring */
366			FVECT nvec,
367			OBJREC *op
368			)
369	greg	1.1	{
370	greg	2.20	CONE *co;
371	greg	1.1
372			co = getcone(op, 0);
373	greg	2.23	if (co == NULL) {
374			memset(nvec, 0, sizeof(FVECT));
375			return(0.);
376			}
377	greg	1.1	VCOPY(nvec, co->ad);
378			return(DOT(nvec, CO_P0(co)));
379	greg	1.4	}
380
381
382	greg	2.9	int
383	greg	2.20	commonspot( /* set sp1 to intersection of sp1 and sp2 */
384			SPOT *sp1,
385			SPOT *sp2,
386			FVECT org
387			)
388	greg	1.4	{
389			FVECT cent;
390			double rad2, cos1, cos2;
391
392			cos1 = 1. - sp1->siz/(2.*PI);
393			cos2 = 1. - sp2->siz/(2.*PI);
394			if (sp2->siz >= 2.PI-FTINY) / BIG, just check overlap */
395			return(DOT(sp1->aim,sp2->aim) >= cos1*cos2 -
396			sqrt((1.-cos1cos1)(1.-cos2*cos2)));
397			/* compute and check disks */
398			rad2 = intercircle(cent, sp1->aim, sp2->aim,
399			1./(cos1cos1) - 1., 1./(cos2cos2) - 1.);
400			if (rad2 <= FTINY \|\| normalize(cent) == 0.)
401			return(0);
402			VCOPY(sp1->aim, cent);
403			sp1->siz = 2.PI(1. - 1./sqrt(1.+rad2));
404			return(1);
405			}
406
407
408	greg	2.9	int
409	greg	2.20	commonbeam( /* set sp1 to intersection of sp1 and sp2 */
410			SPOT *sp1,
411			SPOT *sp2,
412			FVECT dir
413			)
414	greg	1.4	{
415			FVECT cent, c1, c2;
416			double rad2, d;
417			/* move centers to common plane */
418			d = DOT(sp1->aim, dir);
419	greg	2.18	VSUM(c1, sp1->aim, dir, -d);
420	greg	1.4	d = DOT(sp2->aim, dir);
421	greg	2.18	VSUM(c2, sp2->aim, dir, -d);
422	greg	1.4	/* compute overlap */
423			rad2 = intercircle(cent, c1, c2, sp1->siz/PI, sp2->siz/PI);
424			if (rad2 <= FTINY)
425			return(0);
426			VCOPY(sp1->aim, cent);
427			sp1->siz = PI*rad2;
428			return(1);
429			}
430
431
432	greg	2.9	int
433	greg	2.20	checkspot( /* check spotlight for behind source */
434			SPOT sp, / spotlight */
435			FVECT nrm /* source surface normal */
436			)
437	greg	1.4	{
438			double d, d1;
439
440			d = DOT(sp->aim, nrm);
441			if (d > FTINY) /* center in front? */
442	greg	1.8	return(1);
443	greg	1.4	/* else check horizon */
444			d1 = 1. - sp->siz/(2.*PI);
445	greg	1.8	return(1.-FTINY-dd < d1d1);
446	greg	1.4	}
447
448
449			double
450	greg	2.20	spotdisk( /* intersect spot with object op */
451			FVECT oc,
452			OBJREC *op,
453			SPOT *sp,
454			FVECT pos
455			)
456	greg	1.6	{
457			FVECT onorm;
458			double offs, d, dist;
459
460			offs = getplaneq(onorm, op);
461			d = -DOT(onorm, sp->aim);
462			if (d >= -FTINY && d <= FTINY)
463			return(0.);
464			dist = (DOT(pos, onorm) - offs)/d;
465			if (dist < 0.)
466			return(0.);
467	greg	2.18	VSUM(oc, pos, sp->aim, dist);
468	greg	1.6	return(sp->sizdistdist/PI/(d*d));
469			}
470
471
472			double
473	greg	2.20	beamdisk( /* intersect beam with object op */
474			FVECT oc,
475			OBJREC *op,
476			SPOT *sp,
477			FVECT dir
478			)
479	greg	1.6	{
480			FVECT onorm;
481			double offs, d, dist;
482
483			offs = getplaneq(onorm, op);
484			d = -DOT(onorm, dir);
485			if (d >= -FTINY && d <= FTINY)
486			return(0.);
487			dist = (DOT(sp->aim, onorm) - offs)/d;
488	greg	2.18	VSUM(oc, sp->aim, dir, dist);
489	greg	1.6	return(sp->siz/PI/(d*d));
490			}
491
492
493			double
494	greg	2.20	intercircle( /* intersect two circles */
495			FVECT cc, /* midpoint (return value) */
496			FVECT c1, /* circle centers */
497			FVECT c2,
498			double r1s, /* radii squared */
499			double r2s
500			)
501	greg	1.4	{
502			double a2, d2, l;
503			FVECT disp;
504
505	greg	2.18	VSUB(disp, c2, c1);
506	greg	1.4	d2 = DOT(disp,disp);
507			/* circle within overlap? */
508			if (r1s < r2s) {
509			if (r2s >= r1s + d2) {
510			VCOPY(cc, c1);
511			return(r1s);
512			}
513			} else {
514			if (r1s >= r2s + d2) {
515			VCOPY(cc, c2);
516			return(r2s);
517			}
518			}
519			a2 = .25(2.(r1s+r2s) - d2 - (r2s-r1s)*(r2s-r1s)/d2);
520			/* no overlap? */
521			if (a2 <= 0.)
522			return(0.);
523			/* overlap, compute center */
524			l = sqrt((r1s - a2)/d2);
525	greg	2.18	VSUM(cc, c1, disp, l);
526	greg	1.4	return(a2);
527	greg	1.1	}