--- ray/src/rt/ambcomp.c 1991/06/13 09:41:17 1.5 +++ ray/src/rt/ambcomp.c 1991/08/13 12:16:41 1.11 @@ -60,7 +60,8 @@ AMBHEMI *h; RAY *r; { RAY ar; - int hlist[4]; + int hlist[3]; + double spt[2]; double xd, yd, zd; double b2; double phi; @@ -71,10 +72,9 @@ RAY *r; hlist[0] = r->rno; hlist[1] = dp->t; hlist[2] = dp->p; - hlist[3] = 0; - zd = sqrt((dp->t+urand(ilhash(hlist,4)+dp->n))/h->nt); - hlist[3] = 1; - phi = 2.0*PI * (dp->p+urand(ilhash(hlist,4)+dp->n))/h->np; + peano(spt, 2, urand(ilhash(hlist,3)+dp->n), .01); + zd = sqrt((dp->t + spt[0])/h->nt); + phi = 2.0*PI * (dp->p + spt[1])/h->np; xd = cos(phi) * zd; yd = sin(phi) * zd; zd = sqrt(1.0 - zd*zd); @@ -86,7 +86,7 @@ RAY *r; rayvalue(&ar); ndims--; addcolor(dp->v, ar.rcol); - if (ar.rt < FHUGE) + if (ar.rt > FTINY && ar.rt < FHUGE) dp->r += 1.0/ar.rt; /* (re)initialize error */ if (dp->n++) { @@ -180,20 +180,26 @@ FVECT pg, dg; addcolor(acol, dp->v); } b = bright(acol); - if (b > FTINY) + if (b > FTINY) { b = ndivs/b; - else - b = 0.0; - if (pg != NULL) { - posgradient(pg, div, &hemi); - for (i = 0; i < 3; i++) - pg[i] *= b; + if (pg != NULL) { + posgradient(pg, div, &hemi); + for (i = 0; i < 3; i++) + pg[i] *= b; + } + if (dg != NULL) { + dirgradient(dg, div, &hemi); + for (i = 0; i < 3; i++) + dg[i] *= b; + } + } else { + if (pg != NULL) + for (i = 0; i < 3; i++) + pg[i] = 0.0; + if (dg != NULL) + for (i = 0; i < 3; i++) + dg[i] = 0.0; } - if (dg != NULL) { - dirgradient(dg, div, &hemi); - for (i = 0; i < 3; i++) - dg[i] *= b; - } free((char *)div); } b = 1.0/ndivs; @@ -207,6 +213,14 @@ FVECT pg, dg; else if (arad < minarad) arad = minarad; arad /= sqrt(r->rweight); + if (pg != NULL) { /* clip pos. gradient if too large */ + d = 4.0*DOT(pg,pg)*arad*arad; + if (d > 1.0) { + d = 1.0/sqrt(d); + for (i = 0; i < 3; i++) + pg[i] *= d; + } + } return(arad); oopsy: if (div != NULL) @@ -221,8 +235,8 @@ RAY *r; { register int i; /* set number of divisions */ - hp->nt = sqrt(ambdiv * r->rweight * 0.5) + 0.5; - hp->np = 2 * hp->nt; + hp->nt = sqrt(ambdiv * r->rweight / PI) + 0.5; + hp->np = PI * hp->nt; /* make axes */ VCOPY(hp->uz, r->ron); hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0; @@ -249,6 +263,11 @@ register AMBHEMI *hp; dp = da; for (i = 0; i < hp->nt; i++) for (j = 0; j < hp->np; j++) { +#ifdef DEBUG + if (dp->t != i || dp->p != j) + error(CONSISTENCY, + "division order in comperrs"); +#endif b = bright(dp[0].v); if (i > 0) { /* from above */ b2 = bright(dp[-hp->np].v) - b;