ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/data.c
(Generate patch)

Comparing ray/src/rt/data.c (file contents):
Revision 2.37 by greg, Tue Mar 12 16:54:51 2024 UTC vs.
Revision 2.40 by greg, Wed Mar 13 07:24:53 2024 UTC

# Line 373 | Line 373 | freedata(                      /* release data array reference */
373          if (dta == NULL) {                      /* free all if NULL */
374                  hval = 0; nents = TABSIZ;
375          } else {
376 +                if (dta->next == dta) {
377 +                        free(dta);              /* unlisted temp array */
378 +                        return;
379 +                }
380                  hval = hash(dta->name); nents = 1;
381                  if (!*dta->name) {              /* not a data file? */
382                          dta->next = dtab[hval];
# Line 405 | Line 409 | data_interp(DATARRAY *dp, double *pt, double coef, DAT
409          DATARRAY        sd;
410          int             stride, i;
411          double          x, c0, c1, y0, y1;
412 +                                        /* unlikely, but may as well check */
413 +        if ((-FTINY <= coef) & (coef <= FTINY))
414 +                return(0.);
415                                          /* set up dimensions for recursion */
416          if (dp->nd > 1) {
417                  sd.name = dp->name;
# Line 471 | Line 478 | data_interp(DATARRAY *dp, double *pt, double coef, DAT
478                          for (i = sd.dim[0].ne; i--; )
479                                  rvec[i] += c0*sd.arr.d[i]
480                                          + c1*sd.arr.d[i+stride];
481 <                        return(0.);
475 <                }
476 <                if (dp->type == SPECTY) {
481 >                } else if (dp->type == SPECTY) {
482                          double  f;
483                          sd.arr.s = dp->arr.s + i*stride;
484 <                        f = ldexp(1.0, (int)sd.arr.s[sd.dim[0].ne]
485 <                                        - (COLXS+8));
486 <                        for (i = sd.dim[0].ne; i--; )
487 <                                rvec[i] += c0*f*(sd.arr.s[i] + 0.5);
484 >                        if ((sd.arr.s[sd.dim[0].ne] > 0) & ((-FTINY>c0)|(c0>FTINY))) {
485 >                                f = c0*ldexp(1., (int)sd.arr.s[sd.dim[0].ne]-(COLXS+8));
486 >                                for (i = sd.dim[0].ne; i--; )
487 >                                        rvec[i] += f*(sd.arr.s[i] + 0.5);
488 >                        }
489                          sd.arr.s += stride;
490 <                        f = ldexp(1.0, (int)sd.arr.s[sd.dim[0].ne]
491 <                                        - (COLXS+8));
490 >                        if ((sd.arr.s[sd.dim[0].ne] > 0) & ((-FTINY>c1)|(c1>FTINY))) {
491 >                                f = c1*ldexp(1., (int)sd.arr.s[sd.dim[0].ne]-(COLXS+8));
492 >                                for (i = sd.dim[0].ne; i--; )
493 >                                        rvec[i] += f*(sd.arr.s[i] + 0.5);
494 >                        }
495 >                } else {
496 >                        sd.arr.c = dp->arr.c + i*stride;
497                          for (i = sd.dim[0].ne; i--; )
498 <                                rvec[i] += c1*f*(sd.arr.s[i] + 0.5);
499 <                        return(0.);
498 >                                rvec[i] += c0*colrval(sd.arr.c[i],sd.type)
499 >                                        + c1*colrval(sd.arr.c[i+stride],sd.type);
500                  }
501 <                sd.arr.c = dp->arr.c + i*stride;
491 <                for (i = sd.dim[0].ne; i--; )
492 <                        rvec[i] += c0*colrval(sd.arr.c[i],sd.type)
493 <                                + c1*colrval(sd.arr.c[i+stride],sd.type);
494 <                return(0.);
501 >                return(0.);             /* return value ignored */
502          }
503                                          /* get dependent variable */
504          if (dp->nd > 1) {
# Line 515 | Line 522 | data_interp(DATARRAY *dp, double *pt, double coef, DAT
522                          y1 = dp->arr.d[i+1];
523                  } else if (dp->type == SPECTY) {
524                          if (dp->arr.s[dp->dim[0].ne]) {
525 <                                double  f = ldexp(1.0, -(COLXS+8) +
526 <                                                (int)dp->arr.s[dp->dim[0].ne]);
527 <                                y0 = (dp->arr.s[i] + 0.5)*f;
528 <                                y1 = (dp->arr.s[i+1] + 0.5)*f;
525 >                                double  f = dp->arr.s[dp->dim[0].ne]
526 >                                        ? ldexp(1., -(COLXS+8) +
527 >                                                (int)dp->arr.s[dp->dim[0].ne])
528 >                                        : 0.;
529 >                                y0 = f*(dp->arr.s[i] + 0.5);
530 >                                y1 = f*(dp->arr.s[i+1] + 0.5);
531                          } else
532                                  y0 = y1 = 0.0;
533                  } else {
# Line 556 | Line 565 | datavector(DATARRAY *dp, double *pt)
565                                  sizeof(DATATYPE)*dp->dim[dp->nd-1].ne);
566          if (newdp == NULL)
567                  error(SYSTEM, "out of memory in datavector");
568 <        newdp->next = NULL;
568 >        newdp->next = newdp;            /* flags us as temp vector */
569          newdp->name = dp->name;
570          newdp->type = DATATY;
571          newdp->nd = 1;                  /* vector data goes here */
# Line 569 | Line 578 | datavector(DATARRAY *dp, double *pt)
578          return(newdp);                  /* will be free'd using free() */
579   }
580  
572
573 #if 0
574 double
575 datavalue(              /* interpolate data value at a point */
576        DATARRAY  *dp,
577        double  *pt
578 )
579 {
580        DATARRAY  sd;
581        int  asize;
582        int  lower, upper;
583        int  i;
584        double  x, y0, y1;
585                                        /* set up dimensions for recursion */
586        if (dp->nd > 1) {
587                sd.name = dp->name;
588                sd.type = dp->type;
589                sd.nd = dp->nd - 1;
590                asize = 1;
591                for (i = 0; i < sd.nd; i++) {
592                        sd.dim[i].org = dp->dim[i+1].org;
593                        sd.dim[i].siz = dp->dim[i+1].siz;
594                        sd.dim[i].p = dp->dim[i+1].p;
595                        asize *= (sd.dim[i].ne = dp->dim[i+1].ne) +
596                                ((sd.type==SPECTY) & (i==sd.nd-1));
597                }
598        }
599                                        /* get independent variable */
600        if (dp->dim[0].p == NULL) {             /* evenly spaced points */
601                x = (pt[0] - dp->dim[0].org)/dp->dim[0].siz;
602                x *= (double)(dp->dim[0].ne - 1);
603                i = x;
604                if (i < 0)
605                        i = 0;
606                else if (i > dp->dim[0].ne - 2)
607                        i = dp->dim[0].ne - 2;
608        } else {                                /* unevenly spaced points */
609                if (dp->dim[0].siz > 0.0) {
610                        lower = 0;
611                        upper = dp->dim[0].ne;
612                } else {
613                        lower = dp->dim[0].ne;
614                        upper = 0;
615                }
616                do {
617                        i = (lower + upper) >> 1;
618                        if (pt[0] >= dp->dim[0].p[i])
619                                lower = i;
620                        else
621                                upper = i;
622                } while (i != (lower + upper) >> 1);
623
624                if (i > dp->dim[0].ne - 2)
625                        i = dp->dim[0].ne - 2;
626
627                x = i + (pt[0] - dp->dim[0].p[i]) /
628                                (dp->dim[0].p[i+1] - dp->dim[0].p[i]);
629        }
630                                        /* get dependent variable */
631        if (dp->nd > 1) {
632                if (dp->type == DATATY) {
633                        sd.arr.d = dp->arr.d + i*asize;
634                        y0 = datavalue(&sd, pt+1);
635                        sd.arr.d += asize;
636                        y1 = datavalue(&sd, pt+1);
637                } else if (dp->type == SPECTY) {
638                        sd.arr.s = dp->arr.s + i*asize;
639                        y0 = datavalue(&sd, pt+1);
640                        sd.arr.s += asize;
641                        y1 = datavalue(&sd, pt+1);
642                } else {
643                        sd.arr.c = dp->arr.c + i*asize;
644                        y0 = datavalue(&sd, pt+1);
645                        sd.arr.c += asize;
646                        y1 = datavalue(&sd, pt+1);
647                }
648        } else {
649                if (dp->type == DATATY) {
650                        y0 = dp->arr.d[i];
651                        y1 = dp->arr.d[i+1];
652                } else if (dp->type == SPECTY) {
653                        if (dp->arr.s[dp->dim[0].ne]) {
654                                double  f = ldexp(1.0, -(COLXS+8) +
655                                                (int)dp->arr.s[dp->dim[0].ne]);
656                                y0 = (dp->arr.s[i] + 0.5)*f;
657                                y1 = (dp->arr.s[i+1] + 0.5)*f;
658                        } else
659                                y0 = y1 = 0.0;
660                } else {
661                        y0 = colrval(dp->arr.c[i],dp->type);
662                        y1 = colrval(dp->arr.c[i+1],dp->type);
663                }
664        }
665        /*
666         * Extrapolate as far as one division, then
667         * taper off harmonically to zero.
668         */
669        if (x > i+2)
670                return( (2*y1-y0)/(x-(i-1)) );
671
672        if (x < i-1)
673                return( (2*y0-y1)/(i-x) );
674
675        return( y0*((i+1)-x) + y1*(x-i) );
676 }
677 #endif

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines