src/rt/srcsupp.c

#ifndef lint
static const char       RCSid[] = "$Id: srcsupp.c,v 2.24 2022/08/04 22:43:46 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 if (f->nv == 3) {                        /* triangle case */
                int     near0 = 2;
                double  dmin = dist2line(src->sloc, VERTEX(f,2), VERTEX(f,0));
                for (i = 0; i < 2; i++) {
                        double  d2 = dist2line(src->sloc, VERTEX(f,i), VERTEX(f,i+1));
                        if (d2 >= dmin)
                                continue;
                        near0 = i;
                        dmin = d2;                      /* radius = min distance */
                }
                if (dmin < .08*f->area)
                        objerror(so, WARNING, "triangular source with poor aspect");
                i = (near0 + 1) % 3;
                for (j = 0; j < 3; j++)
                        src->ss[SU][j] = VERTEX(f,i)[j] - VERTEX(f,near0)[j];
                normalize(src->ss[SU]);
                dmin = sqrt(dmin);
                for (j = 0; j < 3; j++)
                        src->ss[SU][j] *= dmin;
                fcross(src->ss[SV], f->norm, src->ss[SU]);
        } else
                setflatss(src);                         /* hope for convex! */
}


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