src/hd/rholo3.c

#ifndef lint
static const char       RCSid[] = "$Id: rholo3.c,v 3.44 2018/10/05 19:19:16 greg Exp $";
#endif
/*
 * Routines for tracking beam compuatations
 */

#include "rholo.h"

#ifndef NFRAG2CHUNK
#define NFRAG2CHUNK     4096    /* number of fragments to start chunking */
#endif

#ifndef MAXADISK
#define MAXADISK        10240.  /* maximum holodeck size (Megs) for ambient */
#endif

#ifndef abs
#define abs(x)          ((x) > 0 ? (x) : -(x))
#endif
#ifndef sgn
#define sgn(x)          ((x) > 0 ? 1 : (x) < 0 ? -1 : 0)
#endif

#define rchunk(n)       (((n)+(RPACKSIZ/2))/RPACKSIZ)

int     chunkycmp = 0;          /* clump beams together on disk */

static PACKHEAD *complist=NULL; /* list of beams to compute */
static int      complen=0;      /* length of complist */
static int      listpos=0;      /* current list position for next_packet */
static int      lastin= -1;     /* last ordered position in list */

static void sortcomplist(void);
static void mergeclists(PACKHEAD *cdest, PACKHEAD *cl1, int n1, PACKHEAD *cl2, int n2);
static void view_list(FILE      *fp);
static void ambient_list(void);
static double beamvolume(HOLO   *hp, int        bi);
static void dispbeam(BEAM       *b, HDBEAMI     *hb);


static int
beamcmp(b0, b1)                         /* comparison for compute order */
PACKHEAD        *b0, *b1;
{
        BEAMI   *bip0, *bip1;
        long    c;
                                        /* first check desired quantities */
        if (chunkycmp)
                c = rchunk(b1->nr)*(rchunk(b0->nc)+1L) -
                                rchunk(b0->nr)*(rchunk(b1->nc)+1L);
        else
                c = b1->nr*(b0->nc+1L) - b0->nr*(b1->nc+1L);
        if (c > 0) return(1);
        if (c < 0) return(-1);
                                /* only one file, so skip the following: */
#if 0
                                        /* next, check file descriptors */
        c = hdlist[b0->hd]->fd - hdlist[b1->hd]->fd;
        if (c) return(c);
#endif
                                        /* finally, check file positions */
        bip0 = &hdlist[b0->hd]->bi[b0->bi];
        bip1 = &hdlist[b1->hd]->bi[b1->bi];
                                        /* put diskless beams last */
        if (!bip0->nrd)
                return(bip1->nrd > 0);
        if (!bip1->nrd)
                return(-1);
        c = bip0->fo - bip1->fo;
        return(c < 0 ? -1 : c > 0);
}


int
beamidcmp(b0, b1)                       /* comparison for beam searching */
PACKHEAD        *b0, *b1;
{
        int     c = b0->hd - b1->hd;

        if (c) return(c);
        return(b0->bi - b1->bi);
}


static void
dispbeam(                               /* display a holodeck beam */
        BEAM    *b,
        HDBEAMI *hb
)
{
        static int      n = 0;
        static PACKHEAD *p = NULL;

        if (b == NULL)
                return;
        if (b->nrm > n) {               /* (re)allocate packet holder */
                n = b->nrm;
                if (p == NULL) p = (PACKHEAD *)malloc(packsiz(n));
                else p = (PACKHEAD *)realloc((void *)p, packsiz(n));
                CHECK(p==NULL, SYSTEM, "out of memory in dispbeam");
        }
                                        /* assign packet fields */
        memcpy((void *)packra(p), (void *)hdbray(b), b->nrm*sizeof(RAYVAL));
        p->nr = p->nc = b->nrm;
        for (p->hd = 0; hdlist[p->hd] != hb->h; p->hd++)
                if (hdlist[p->hd] == NULL)
                        error(CONSISTENCY, "unregistered holodeck in dispbeam");
        p->bi = hb->b;
        disp_packet(p);                 /* display it */
        if (n >= 1024) {                /* free ridiculous packets */
                free((void *)p);
                p = NULL; n = 0;
        }
}


void
bundle_set(     /* bundle set operation */
        int     op,
        PACKHEAD        *clist,
        int     nents
)
{
        int     oldnr, n;
        HDBEAMI *hbarr;
        PACKHEAD        *csm;
        int     i;
                                        /* search for common members */
        for (csm = clist+nents; csm-- > clist; )
                csm->nc = -1;
        qsort((void *)clist, nents, sizeof(PACKHEAD), beamidcmp);
        for (i = 0; i < complen; i++) {
                csm = (PACKHEAD *)bsearch((void *)(complist+i), (void *)clist,
                                nents, sizeof(PACKHEAD), beamidcmp);
                if (csm == NULL)
                        continue;
                oldnr = complist[i].nr;
                csm->nc = complist[i].nc;
                switch (op) {
                case BS_ADD:            /* add to count */
                        complist[i].nr += csm->nr;
                        csm->nr = 0;
                        break;
                case BS_MAX:            /* maximum of counts */
                        if (csm->nr > complist[i].nr)
                                complist[i].nr = csm->nr;
                        csm->nr = 0;
                        break;
                case BS_ADJ:            /* reset count */
                        complist[i].nr = csm->nr;
                        csm->nr = 0;
                        break;
                case BS_DEL:            /* delete count */
                        if (csm->nr == 0 || csm->nr >= complist[i].nr)
                                complist[i].nr = 0;
                        else
                                complist[i].nr -= csm->nr;
                        break;
                }
                if (complist[i].nr != oldnr)
                        lastin = -1;    /* flag sort */
        }
                                /* record computed rays for uncommon beams */
        for (csm = clist+nents; csm-- > clist; )
                if (csm->nc < 0)
                        csm->nc = bnrays(hdlist[csm->hd], csm->bi);
                                /* complete list operations */
        switch (op) {
        case BS_NEW:                    /* new computation set */
                listpos = 0; lastin = -1;
                if (complen)            /* free old list */
                        free((void *)complist);
                complist = NULL;
                if (!(complen = nents))
                        return;
                complist = (PACKHEAD *)malloc(nents*sizeof(PACKHEAD));
                if (complist == NULL)
                        goto memerr;
                memcpy((void *)complist, (void *)clist, nents*sizeof(PACKHEAD));
                break;
        case BS_ADD:                    /* add to computation set */
        case BS_MAX:                    /* maximum of quantities */
        case BS_ADJ:                    /* adjust set quantities */
                if (nents <= 0)
                        return;
                sortcomplist();         /* sort updated list & new entries */
                qsort((void *)clist, nents, sizeof(PACKHEAD), beamcmp);
                                        /* what can't we satisfy? */
                for (i = nents, csm = clist; i-- && csm->nr > csm->nc; csm++)
                        ;
                n = csm - clist;
                if (op != BS_ADD) {     /* don't regenerate adjusted beams */
                        for (++i; i-- && csm->nr > 0; csm++)
                                ;
                        nents = csm - clist;
                }
                if (n) {                /* allocate space for merged list */
                        PACKHEAD        *newlist;
                        newlist = (PACKHEAD *)malloc(
                                        (complen+n)*sizeof(PACKHEAD) );
                        if (newlist == NULL)
                                goto memerr;
                                                /* merge lists */
                        mergeclists(newlist, clist, n, complist, complen);
                        if (complen)
                                free((void *)complist);
                        complist = newlist;
                        complen += n;
                }
                listpos = 0;
                lastin = complen-1;     /* list is now sorted */
                break;
        case BS_DEL:                    /* delete from computation set */
                return;                 /* already done */
        default:
                error(CONSISTENCY, "bundle_set called with unknown operation");
        }
        if (outdev == NULL || !nents)   /* nothing to display? */
                return;
                                        /* load and display beams we have */
        hbarr = (HDBEAMI *)malloc(nents*sizeof(HDBEAMI));
        for (i = nents; i--; ) {
                hbarr[i].h = hdlist[clist[i].hd];
                hbarr[i].b = clist[i].bi;
        }
        hdloadbeams(hbarr, nents, dispbeam);
        free((void *)hbarr);
        if (hdfragflags&FF_READ) {
                listpos = 0;
                lastin = -1;            /* need to re-sort list */
        }
        return;
memerr:
        error(SYSTEM, "out of memory in bundle_set");
}


static double
beamvolume(     /* compute approximate volume of a beam */
        HOLO    *hp,
        int     bi
)
{
        GCOORD  gc[2];
        FVECT   cp[4], edgeA, edgeB, cent[2];
        FVECT   crossp[2], diffv;
        double  vol[2];
        int     i;
                                        /* get grid coordinates */
        if (!hdbcoord(gc, hp, bi))
                error(CONSISTENCY, "bad beam index in beamvolume");
        for (i = 0; i < 2; i++) {       /* compute cell area vectors */
                hdcell(cp, hp, gc+i);
                VSUM(edgeA, cp[1], cp[0], -1.0);
                VSUM(edgeB, cp[3], cp[1], -1.0);
                fcross(crossp[i], edgeA, edgeB);
                                        /* compute center */
                cent[i][0] = 0.5*(cp[0][0] + cp[2][0]);
                cent[i][1] = 0.5*(cp[0][1] + cp[2][1]);
                cent[i][2] = 0.5*(cp[0][2] + cp[2][2]);
        }
                                        /* compute difference vector */
        VSUM(diffv, cent[1], cent[0], -1.0);
        for (i = 0; i < 2; i++) {       /* compute volume contributions */
                vol[i] = 0.5*DOT(crossp[i], diffv);
                if (vol[i] < 0.) vol[i] = -vol[i];
        }
        return(vol[0] + vol[1]);        /* return total volume */
}


static void
ambient_list(void)                      /* compute ambient beam list */
{
        unsigned long   wtotal;
        int32   minrt;
        double  frac;
        int     i;
        int     j, k;

        complen = 0;
        for (j = 0; hdlist[j] != NULL; j++)
                complen += nbeams(hdlist[j]);
        complist = (PACKHEAD *)malloc(complen*sizeof(PACKHEAD));
        CHECK(complist==NULL, SYSTEM, "out of memory in ambient_list");
                                        /* compute beam weights */
        k = 0; wtotal = 0;
        for (j = 0; hdlist[j] != NULL; j++) {
                                        /* 512. arbitrary -- adjusted below */
                frac = 512. * VLEN(hdlist[j]->wg[0]) *
                                VLEN(hdlist[j]->wg[1]) *
                                VLEN(hdlist[j]->wg[2]);
                for (i = nbeams(hdlist[j]); i > 0; i--) {
                        complist[k].hd = j;
                        complist[k].bi = i;
                        complist[k].nr = frac*beamvolume(hdlist[j], i) + 0.5;
                        complist[k].nc = bnrays(hdlist[j], i);
                        wtotal += complist[k++].nr;
                }
        }
                                        /* adjust sample weights */
        if (vdef(DISKSPACE))
                frac = 1024.*1024.*vflt(DISKSPACE) / (wtotal*sizeof(RAYVAL));
        else
                frac = 1024.*1024.*MAXADISK / (wtotal*sizeof(RAYVAL));
        minrt = .02*frac*wtotal/complen + 1.1;  /* heuristic mimimum */
        if (minrt > RPACKSIZ)
                minrt = RPACKSIZ;
        for (k = complen; k--; )
                if ((complist[k].nr = frac*complist[k].nr + 0.5) < minrt)
                        complist[k].nr = minrt;
        listpos = 0; lastin = -1;       /* flag initial sort */
}


static void
view_list(                      /* assign beam priority from view list */
        FILE    *fp
)
{
        double  pa = 1.;
        VIEW    curview;
        int     xr, yr;
        char    *err;
        BEAMLIST        blist;

        curview = stdview;
        while (nextview(&curview, fp) != EOF) {
                if ((err = setview(&curview)) != NULL) {
                        error(WARNING, err);
                        continue;
                }
                xr = yr = 1024;
                normaspect(viewaspect(&curview), &pa, &xr, &yr);
                viewbeams(&curview, xr, yr, &blist);
                bundle_set(BS_MAX, blist.bl, blist.nb);
                free((void *)blist.bl);
        }
}


void
init_global(void)                       /* initialize global ray computation */
{
                                        /* free old list and empty queue */
        if (complen > 0) {
                free((void *)complist);
                done_packets(flush_queue());
        }
                                        /* reseed random number generator */
        srandom(time(NULL));
                                        /* allocate beam list */
        if (readinp)
                view_list(stdin);
        else
                ambient_list();
                                        /* no view vicinity */
        myeye.rng = 0;
}


static void
mergeclists(    /* merge two sorted lists */
        PACKHEAD        *cdest,
        PACKHEAD        *cl1,
        int     n1,
        PACKHEAD        *cl2,
        int     n2
)
{
        int     cmp;

        while (n1 | n2) {
                if (!n1) cmp = 1;
                else if (!n2) cmp = -1;
                else cmp = beamcmp(cl1, cl2);
                if (cmp > 0) {
                        *cdest = *cl2;
                        cl2++; n2--;
                } else {
                        *cdest = *cl1;
                        cl1++; n1--;
                }
                cdest++;
        }
}


static void
sortcomplist(void)                      /* fix our list order */
{
        PACKHEAD        *list2;
        int     listlen;
        int     i;

        if (complen <= 0)       /* check to see if there is even a list */
                return;
        if (!chunkycmp)         /* check to see if fragment list is full */
                if (!hdfragOK(hdlist[0]->fd, &listlen, NULL)
#if NFRAG2CHUNK
                                || listlen >= NFRAG2CHUNK
#endif
                                ) {
                        chunkycmp++;    /* use "chunky" comparison */
                        lastin = -1;    /* need to re-sort list */
#ifdef DEBUG
                        error(WARNING, "using chunky comparison mode");
#endif
                }
        if (lastin < 0 || listpos*4 >= complen*3)
                qsort((void *)complist, complen, sizeof(PACKHEAD), beamcmp);
        else if (listpos) {     /* else sort and merge sublist */
                list2 = (PACKHEAD *)malloc(listpos*sizeof(PACKHEAD));
                CHECK(list2==NULL, SYSTEM, "out of memory in sortcomplist");
                memcpy((void *)list2,(void *)complist,listpos*sizeof(PACKHEAD));
                qsort((void *)list2, listpos, sizeof(PACKHEAD), beamcmp);
                mergeclists(complist, list2, listpos,
                                complist+listpos, complen-listpos);
                free((void *)list2);
        }
                                        /* drop satisfied requests */
        for (i = complen; i-- && complist[i].nr <= complist[i].nc; )
                ;
        if (i < 0) {
                free((void *)complist);
                complist = NULL;
                complen = 0;
        } else if (i < complen-1) {
                list2 = (PACKHEAD *)realloc((void *)complist,
                                (i+1)*sizeof(PACKHEAD));
                if (list2 != NULL)
                        complist = list2;
                complen = i+1;
        }
        listpos = 0; lastin = i;
}


/*
 * The following routine works on the assumption that the bundle weights are
 * more or less evenly distributed, such that computing a packet causes
 * a given bundle to move way down in the computation order.  We keep
 * track of where the computed bundle with the highest priority would end
 * up, and if we get further in our compute list than this, we re-sort the
 * list and start again from the beginning.  Since
 * a merge sort is used, the sorting costs are minimal.
 */
int
next_packet(            /* prepare packet for computation */
        PACKET  *p,
        int     n
)
{
        if (listpos > lastin)           /* time to sort the list */
                sortcomplist();
        if (complen <= 0)
                return(0);
        p->hd = complist[listpos].hd;
        p->bi = complist[listpos].bi;
        p->nc = complist[listpos].nc;
        p->nr = complist[listpos].nr - p->nc;
        if (p->nr <= 0)
                return(0);
        DCHECK(n < 1 | n > RPACKSIZ,
                        CONSISTENCY, "next_packet called with bad n value");
        if (p->nr > n)
                p->nr = n;
        complist[listpos].nc += p->nr;  /* find where this one would go */
        if (hdgetbeam(hdlist[p->hd], p->bi) != NULL)
                hdfreefrag(hdlist[p->hd], p->bi);
        while (lastin > listpos && 
                        beamcmp(complist+lastin, complist+listpos) > 0)
                lastin--;
        listpos++;
        return(1);
}
Revision:	3.45
Committed:	Tue Jan 21 22:30:01 2025 UTC (3 months, 2 weeks ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 3.44:	+4 -3 lines
Log Message:	fix(rholo): Minimum # rays/packet was zero in some cases
#	User	Rev	Content
1	gwlarson	3.23	#ifndef lint
2	greg	3.45	static const char RCSid[] = "$Id: rholo3.c,v 3.44 2018/10/05 19:19:16 greg Exp $";
3	gwlarson	3.23	#endif
4	gregl	3.1	/*
5			* Routines for tracking beam compuatations
6			*/
7
8			#include "rholo.h"
9
10	gwlarson	3.27	#ifndef NFRAG2CHUNK
11			#define NFRAG2CHUNK 4096 /* number of fragments to start chunking */
12			#endif
13
14	greg	3.43	#ifndef MAXADISK
15			#define MAXADISK 10240. /* maximum holodeck size (Megs) for ambient */
16			#endif
17
18	gwlarson	3.27	#ifndef abs
19	gregl	3.1	#define abs(x) ((x) > 0 ? (x) : -(x))
20	gwlarson	3.27	#endif
21			#ifndef sgn
22	gregl	3.1	#define sgn(x) ((x) > 0 ? 1 : (x) < 0 ? -1 : 0)
23	gwlarson	3.27	#endif
24	gregl	3.1
25	gwlarson	3.27	#define rchunk(n) (((n)+(RPACKSIZ/2))/RPACKSIZ)
26
27	gwlarson	3.30	int chunkycmp = 0; /* clump beams together on disk */
28
29	gregl	3.4	static PACKHEAD complist=NULL; / list of beams to compute */
30			static int complen=0; /* length of complist */
31			static int listpos=0; /* current list position for next_packet */
32			static int lastin= -1; /* last ordered position in list */
33	gregl	3.1
34	schorsch	3.41	static void sortcomplist(void);
35			static void mergeclists(PACKHEAD cdest, PACKHEAD cl1, int n1, PACKHEAD *cl2, int n2);
36			static void view_list(FILE *fp);
37			static void ambient_list(void);
38			static double beamvolume(HOLO *hp, int bi);
39			static void dispbeam(BEAM b, HDBEAMI hb);
40
41	gregl	3.1
42	schorsch	3.41
43			static int
44	gwlarson	3.23	beamcmp(b0, b1) /* comparison for compute order */
45	greg	3.44	PACKHEAD b0, b1;
46	gregl	3.2	{
47	gwlarson	3.27	BEAMI bip0, bip1;
48	greg	3.44	long c;
49	gwlarson	3.27	/* first check desired quantities */
50	gwlarson	3.30	if (chunkycmp)
51	gwlarson	3.28	c = rchunk(b1->nr)*(rchunk(b0->nc)+1L) -
52			rchunk(b0->nr)*(rchunk(b1->nc)+1L);
53	gwlarson	3.27	else
54	gwlarson	3.28	c = b1->nr(b0->nc+1L) - b0->nr(b1->nc+1L);
55			if (c > 0) return(1);
56			if (c < 0) return(-1);
57	gwlarson	3.27	/* only one file, so skip the following: */
58			#if 0
59			/* next, check file descriptors */
60			c = hdlist[b0->hd]->fd - hdlist[b1->hd]->fd;
61			if (c) return(c);
62			#endif
63			/* finally, check file positions */
64			bip0 = &hdlist[b0->hd]->bi[b0->bi];
65			bip1 = &hdlist[b1->hd]->bi[b1->bi];
66			/* put diskless beams last */
67			if (!bip0->nrd)
68			return(bip1->nrd > 0);
69			if (!bip1->nrd)
70			return(-1);
71			c = bip0->fo - bip1->fo;
72			return(c < 0 ? -1 : c > 0);
73	gregl	3.2	}
74
75
76	gregl	3.14	int
77	gwlarson	3.23	beamidcmp(b0, b1) /* comparison for beam searching */
78	greg	3.44	PACKHEAD b0, b1;
79	gwlarson	3.23	{
80	greg	3.44	int c = b0->hd - b1->hd;
81	gwlarson	3.23
82			if (c) return(c);
83			return(b0->bi - b1->bi);
84			}
85
86
87	schorsch	3.41	static void
88			dispbeam( /* display a holodeck beam */
89	greg	3.44	BEAM *b,
90			HDBEAMI *hb
91	schorsch	3.41	)
92	gregl	3.14	{
93			static int n = 0;
94			static PACKHEAD *p = NULL;
95
96			if (b == NULL)
97			return;
98			if (b->nrm > n) { /* (re)allocate packet holder */
99			n = b->nrm;
100			if (p == NULL) p = (PACKHEAD *)malloc(packsiz(n));
101	greg	3.35	else p = (PACKHEAD )realloc((void )p, packsiz(n));
102	gwlarson	3.31	CHECK(p==NULL, SYSTEM, "out of memory in dispbeam");
103	gregl	3.14	}
104			/* assign packet fields */
105	schorsch	3.38	memcpy((void )packra(p), (void )hdbray(b), b->nrm*sizeof(RAYVAL));
106	gregl	3.14	p->nr = p->nc = b->nrm;
107	gwlarson	3.23	for (p->hd = 0; hdlist[p->hd] != hb->h; p->hd++)
108	gregl	3.14	if (hdlist[p->hd] == NULL)
109			error(CONSISTENCY, "unregistered holodeck in dispbeam");
110	gwlarson	3.23	p->bi = hb->b;
111	gregl	3.14	disp_packet(p); /* display it */
112	gwlarson	3.29	if (n >= 1024) { /* free ridiculous packets */
113	greg	3.34	free((void *)p);
114	gwlarson	3.29	p = NULL; n = 0;
115			}
116	gregl	3.14	}
117
118
119	greg	3.44	void
120	schorsch	3.41	bundle_set( /* bundle set operation */
121			int op,
122			PACKHEAD *clist,
123			int nents
124			)
125	gregl	3.2	{
126	gwlarson	3.23	int oldnr, n;
127			HDBEAMI *hbarr;
128	greg	3.44	PACKHEAD *csm;
129			int i;
130	gwlarson	3.23	/* search for common members */
131			for (csm = clist+nents; csm-- > clist; )
132			csm->nc = -1;
133	greg	3.36	qsort((void *)clist, nents, sizeof(PACKHEAD), beamidcmp);
134	gwlarson	3.23	for (i = 0; i < complen; i++) {
135	greg	3.36	csm = (PACKHEAD )bsearch((void )(complist+i), (void *)clist,
136	gwlarson	3.23	nents, sizeof(PACKHEAD), beamidcmp);
137			if (csm == NULL)
138			continue;
139			oldnr = complist[i].nr;
140			csm->nc = complist[i].nc;
141			switch (op) {
142			case BS_ADD: /* add to count */
143			complist[i].nr += csm->nr;
144			csm->nr = 0;
145			break;
146	gwlarson	3.31	case BS_MAX: /* maximum of counts */
147			if (csm->nr > complist[i].nr)
148			complist[i].nr = csm->nr;
149			csm->nr = 0;
150			break;
151	gwlarson	3.23	case BS_ADJ: /* reset count */
152			complist[i].nr = csm->nr;
153			csm->nr = 0;
154			break;
155			case BS_DEL: /* delete count */
156			if (csm->nr == 0 \|\| csm->nr >= complist[i].nr)
157			complist[i].nr = 0;
158			else
159			complist[i].nr -= csm->nr;
160			break;
161			}
162			if (complist[i].nr != oldnr)
163			lastin = -1; /* flag sort */
164	gregl	3.20	}
165	gwlarson	3.24	/* record computed rays for uncommon beams */
166	gwlarson	3.23	for (csm = clist+nents; csm-- > clist; )
167			if (csm->nc < 0)
168			csm->nc = bnrays(hdlist[csm->hd], csm->bi);
169	gregl	3.20	/* complete list operations */
170	gregl	3.2	switch (op) {
171			case BS_NEW: /* new computation set */
172	gwlarson	3.21	listpos = 0; lastin = -1;
173	gregl	3.20	if (complen) /* free old list */
174	greg	3.34	free((void *)complist);
175	gregl	3.20	complist = NULL;
176			if (!(complen = nents))
177	gregl	3.2	return;
178			complist = (PACKHEAD )malloc(nentssizeof(PACKHEAD));
179			if (complist == NULL)
180			goto memerr;
181	schorsch	3.38	memcpy((void )complist, (void )clist, nents*sizeof(PACKHEAD));
182	gregl	3.2	break;
183			case BS_ADD: /* add to computation set */
184	gwlarson	3.31	case BS_MAX: /* maximum of quantities */
185	gregl	3.11	case BS_ADJ: /* adjust set quantities */
186	gregl	3.2	if (nents <= 0)
187			return;
188	gregl	3.20	sortcomplist(); /* sort updated list & new entries */
189	greg	3.36	qsort((void *)clist, nents, sizeof(PACKHEAD), beamcmp);
190	gregl	3.2	/* what can't we satisfy? */
191	gwlarson	3.23	for (i = nents, csm = clist; i-- && csm->nr > csm->nc; csm++)
192	gregl	3.2	;
193	gwlarson	3.23	n = csm - clist;
194	gwlarson	3.31	if (op != BS_ADD) { /* don't regenerate adjusted beams */
195	gwlarson	3.24	for (++i; i-- && csm->nr > 0; csm++)
196	gwlarson	3.22	;
197	gwlarson	3.24	nents = csm - clist;
198	gwlarson	3.22	}
199	gregl	3.2	if (n) { /* allocate space for merged list */
200			PACKHEAD *newlist;
201			newlist = (PACKHEAD *)malloc(
202			(complen+n)*sizeof(PACKHEAD) );
203			if (newlist == NULL)
204			goto memerr;
205			/* merge lists */
206			mergeclists(newlist, clist, n, complist, complen);
207			if (complen)
208	greg	3.34	free((void *)complist);
209	gregl	3.2	complist = newlist;
210			complen += n;
211			}
212			listpos = 0;
213			lastin = complen-1; /* list is now sorted */
214			break;
215			case BS_DEL: /* delete from computation set */
216	gregl	3.20	return; /* already done */
217	gregl	3.2	default:
218			error(CONSISTENCY, "bundle_set called with unknown operation");
219			}
220	gwlarson	3.27	if (outdev == NULL \|\| !nents) /* nothing to display? */
221	gwlarson	3.22	return;
222			/* load and display beams we have */
223	gwlarson	3.23	hbarr = (HDBEAMI )malloc(nentssizeof(HDBEAMI));
224			for (i = nents; i--; ) {
225			hbarr[i].h = hdlist[clist[i].hd];
226			hbarr[i].b = clist[i].bi;
227	gregl	3.11	}
228	gwlarson	3.23	hdloadbeams(hbarr, nents, dispbeam);
229	greg	3.34	free((void *)hbarr);
230	gwlarson	3.27	if (hdfragflags&FF_READ) {
231			listpos = 0;
232			lastin = -1; /* need to re-sort list */
233			}
234	gregl	3.2	return;
235			memerr:
236			error(SYSTEM, "out of memory in bundle_set");
237			}
238
239
240	schorsch	3.41	static double
241			beamvolume( /* compute approximate volume of a beam */
242			HOLO *hp,
243			int bi
244			)
245	gregl	3.1	{
246	gregl	3.13	GCOORD gc[2];
247			FVECT cp[4], edgeA, edgeB, cent[2];
248	schorsch	3.41	FVECT crossp[2], diffv;
249	gregl	3.13	double vol[2];
250	greg	3.44	int i;
251	gregl	3.13	/* get grid coordinates */
252			if (!hdbcoord(gc, hp, bi))
253			error(CONSISTENCY, "bad beam index in beamvolume");
254			for (i = 0; i < 2; i++) { /* compute cell area vectors */
255			hdcell(cp, hp, gc+i);
256			VSUM(edgeA, cp[1], cp[0], -1.0);
257			VSUM(edgeB, cp[3], cp[1], -1.0);
258			fcross(crossp[i], edgeA, edgeB);
259			/* compute center */
260			cent[i][0] = 0.5*(cp[0][0] + cp[2][0]);
261			cent[i][1] = 0.5*(cp[0][1] + cp[2][1]);
262			cent[i][2] = 0.5*(cp[0][2] + cp[2][2]);
263	gregl	3.1	}
264	gregl	3.13	/* compute difference vector */
265			VSUM(diffv, cent[1], cent[0], -1.0);
266			for (i = 0; i < 2; i++) { /* compute volume contributions */
267	gregl	3.15	vol[i] = 0.5*DOT(crossp[i], diffv);
268	gregl	3.13	if (vol[i] < 0.) vol[i] = -vol[i];
269			}
270			return(vol[0] + vol[1]); /* return total volume */
271	gregl	3.1	}
272
273
274	schorsch	3.41	static void
275			ambient_list(void) /* compute ambient beam list */
276	gregl	3.1	{
277	greg	3.45	unsigned long wtotal;
278			int32 minrt;
279	gregl	3.1	double frac;
280	gregl	3.13	int i;
281	greg	3.44	int j, k;
282	gwlarson	3.31
283	gregl	3.1	complen = 0;
284			for (j = 0; hdlist[j] != NULL; j++)
285			complen += nbeams(hdlist[j]);
286			complist = (PACKHEAD )malloc(complensizeof(PACKHEAD));
287	gwlarson	3.31	CHECK(complist==NULL, SYSTEM, "out of memory in ambient_list");
288	gregl	3.1	/* compute beam weights */
289	gwlarson	3.32	k = 0; wtotal = 0;
290	gregl	3.13	for (j = 0; hdlist[j] != NULL; j++) {
291	greg	3.43	/* 512. arbitrary -- adjusted below */
292	gregl	3.19	frac = 512. * VLEN(hdlist[j]->wg[0]) *
293			VLEN(hdlist[j]->wg[1]) *
294			VLEN(hdlist[j]->wg[2]);
295	gregl	3.1	for (i = nbeams(hdlist[j]); i > 0; i--) {
296			complist[k].hd = j;
297			complist[k].bi = i;
298	gregl	3.13	complist[k].nr = frac*beamvolume(hdlist[j], i) + 0.5;
299	gregl	3.18	complist[k].nc = bnrays(hdlist[j], i);
300	gregl	3.1	wtotal += complist[k++].nr;
301			}
302	gregl	3.13	}
303	gwlarson	3.31	/* adjust sample weights */
304	gregl	3.12	if (vdef(DISKSPACE))
305	gwlarson	3.32	frac = 1024.1024.vflt(DISKSPACE) / (wtotal*sizeof(RAYVAL));
306	gregl	3.12	else
307	greg	3.43	frac = 1024.1024.MAXADISK / (wtotal*sizeof(RAYVAL));
308	greg	3.45	minrt = .02fracwtotal/complen + 1.1; /* heuristic mimimum */
309	gwlarson	3.32	if (minrt > RPACKSIZ)
310			minrt = RPACKSIZ;
311			for (k = complen; k--; )
312			if ((complist[k].nr = frac*complist[k].nr + 0.5) < minrt)
313			complist[k].nr = minrt;
314	gwlarson	3.31	listpos = 0; lastin = -1; /* flag initial sort */
315			}
316
317
318	schorsch	3.41	static void
319			view_list( /* assign beam priority from view list */
320			FILE *fp
321			)
322	gwlarson	3.31	{
323			double pa = 1.;
324			VIEW curview;
325			int xr, yr;
326			char *err;
327			BEAMLIST blist;
328
329	schorsch	3.40	curview = stdview;
330	gwlarson	3.31	while (nextview(&curview, fp) != EOF) {
331			if ((err = setview(&curview)) != NULL) {
332			error(WARNING, err);
333			continue;
334			}
335			xr = yr = 1024;
336			normaspect(viewaspect(&curview), &pa, &xr, &yr);
337			viewbeams(&curview, xr, yr, &blist);
338			bundle_set(BS_MAX, blist.bl, blist.nb);
339	greg	3.34	free((void *)blist.bl);
340	gwlarson	3.31	}
341			}
342
343
344	greg	3.44	void
345	schorsch	3.41	init_global(void) /* initialize global ray computation */
346	gwlarson	3.31	{
347			/* free old list and empty queue */
348			if (complen > 0) {
349	greg	3.34	free((void *)complist);
350	gwlarson	3.31	done_packets(flush_queue());
351			}
352			/* reseed random number generator */
353			srandom(time(NULL));
354			/* allocate beam list */
355			if (readinp)
356			view_list(stdin);
357			else
358			ambient_list();
359	gwlarson	3.25	/* no view vicinity */
360			myeye.rng = 0;
361	gregl	3.1	}
362
363
364	schorsch	3.41	static void
365			mergeclists( /* merge two sorted lists */
366	greg	3.44	PACKHEAD *cdest,
367			PACKHEAD *cl1,
368	schorsch	3.41	int n1,
369	greg	3.44	PACKHEAD *cl2,
370	schorsch	3.41	int n2
371			)
372	gregl	3.1	{
373	greg	3.44	int cmp;
374	gregl	3.1
375			while (n1 \| n2) {
376			if (!n1) cmp = 1;
377			else if (!n2) cmp = -1;
378			else cmp = beamcmp(cl1, cl2);
379			if (cmp > 0) {
380	schorsch	3.40	cdest = cl2;
381	gregl	3.1	cl2++; n2--;
382			} else {
383	schorsch	3.40	cdest = cl1;
384	gregl	3.1	cl1++; n1--;
385			}
386			cdest++;
387			}
388			}
389
390
391	schorsch	3.41	static void
392			sortcomplist(void) /* fix our list order */
393	gregl	3.1	{
394			PACKHEAD *list2;
395	gwlarson	3.27	int listlen;
396	greg	3.44	int i;
397	gregl	3.2
398	gregl	3.3	if (complen <= 0) /* check to see if there is even a list */
399	gregl	3.2	return;
400	gwlarson	3.30	if (!chunkycmp) /* check to see if fragment list is full */
401	gwlarson	3.27	if (!hdfragOK(hdlist[0]->fd, &listlen, NULL)
402			#if NFRAG2CHUNK
403			\|\| listlen >= NFRAG2CHUNK
404			#endif
405			) {
406	gwlarson	3.30	chunkycmp++; /* use "chunky" comparison */
407			lastin = -1; /* need to re-sort list */
408	gwlarson	3.27	#ifdef DEBUG
409			error(WARNING, "using chunky comparison mode");
410			#endif
411			}
412	gregl	3.6	if (lastin < 0 \|\| listpos4 >= complen3)
413	greg	3.36	qsort((void *)complist, complen, sizeof(PACKHEAD), beamcmp);
414	gregl	3.1	else if (listpos) { /* else sort and merge sublist */
415			list2 = (PACKHEAD )malloc(listpossizeof(PACKHEAD));
416	gwlarson	3.31	CHECK(list2==NULL, SYSTEM, "out of memory in sortcomplist");
417	schorsch	3.38	memcpy((void )list2,(void )complist,listpos*sizeof(PACKHEAD));
418	greg	3.36	qsort((void *)list2, listpos, sizeof(PACKHEAD), beamcmp);
419	gregl	3.1	mergeclists(complist, list2, listpos,
420			complist+listpos, complen-listpos);
421	greg	3.34	free((void *)list2);
422	gregl	3.1	}
423	gregl	3.2	/* drop satisfied requests */
424	gregl	3.11	for (i = complen; i-- && complist[i].nr <= complist[i].nc; )
425	gregl	3.2	;
426	gregl	3.4	if (i < 0) {
427	greg	3.34	free((void *)complist);
428	gregl	3.4	complist = NULL;
429			complen = 0;
430			} else if (i < complen-1) {
431	greg	3.35	list2 = (PACKHEAD )realloc((void )complist,
432	gregl	3.2	(i+1)*sizeof(PACKHEAD));
433	gwlarson	3.26	if (list2 != NULL)
434	gregl	3.2	complist = list2;
435	gwlarson	3.26	complen = i+1;
436	gregl	3.2	}
437			listpos = 0; lastin = i;
438	gregl	3.1	}
439
440
441			/*
442			* The following routine works on the assumption that the bundle weights are
443			* more or less evenly distributed, such that computing a packet causes
444			* a given bundle to move way down in the computation order. We keep
445			* track of where the computed bundle with the highest priority would end
446	gwlarson	3.27	* up, and if we get further in our compute list than this, we re-sort the
447	gregl	3.11	* list and start again from the beginning. Since
448			* a merge sort is used, the sorting costs are minimal.
449	gregl	3.1	*/
450	greg	3.44	int
451	schorsch	3.41	next_packet( /* prepare packet for computation */
452	greg	3.44	PACKET *p,
453	schorsch	3.41	int n
454			)
455	gregl	3.1	{
456	gregl	3.10	if (listpos > lastin) /* time to sort the list */
457			sortcomplist();
458	gregl	3.1	if (complen <= 0)
459			return(0);
460			p->hd = complist[listpos].hd;
461			p->bi = complist[listpos].bi;
462	gregl	3.11	p->nc = complist[listpos].nc;
463			p->nr = complist[listpos].nr - p->nc;
464	gregl	3.1	if (p->nr <= 0)
465			return(0);
466	gwlarson	3.28	DCHECK(n < 1 \| n > RPACKSIZ,
467			CONSISTENCY, "next_packet called with bad n value");
468	gregl	3.17	if (p->nr > n)
469			p->nr = n;
470	gregl	3.11	complist[listpos].nc += p->nr; /* find where this one would go */
471	gwlarson	3.27	if (hdgetbeam(hdlist[p->hd], p->bi) != NULL)
472			hdfreefrag(hdlist[p->hd], p->bi);
473	gregl	3.11	while (lastin > listpos &&
474			beamcmp(complist+lastin, complist+listpos) > 0)
475	gregl	3.1	lastin--;
476			listpos++;
477			return(1);
478			}