/* this is a 1D table, we can splurge for lots of entries and it still 
   won't be too large to store. */

#define TABLESIZE 16384

int initialized=0;
float table[TABLESIZE];


static double RingTexture(double *pos, double spotsize,
			  double inner, double outer, double transition,
			  double varyscale, double ampratio,
			  double low, double high)
{
  double R0, R1, weight;
  int I0, I1;
  
  /* We need a color table premade for us with the accumulated densities. */
  if (!initialized) MakeColorTable(inner, outer, transition,
				   varyscale, ampratio, low, high);
  
  /* OK, lets find our range of radii. The range of radii are centered
     at the point's radius from the origin, and we split the spot radius
     between the two.
     To tweak the antialiasing, we could scale the spotsize up or down
     here before we use it. 
  */
  
  
  R0=sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2])-0.5*spotsize;
  R1=R0+spotsize;
  
  /* R0 may be negative if we  have huge spot sizes, which is obviously
     wrong and awkward. Lets make sure it's at least 0.0. */
  
  R0=MAX(R0, 0.0);
  
  /* OK, we know what range we want to average over. Lets transform those
     into index numbers in our color table. We know that at radius = [outer],
     we want to be at index [TABLESIZE-1]. So we multiply by a
     constant.
     
     Note that we lose a little resolution because we quantize the
     index into an integer. If the table were smaller, we could
     compensate by using a linear interpolation, but with a really big 
     table the effect is negligable, and it definately makes the codxe simpler.
     
  */
  
  I0=(int)(R0*(TABLESIZE-1)/outer);
  I1=(int)(R1*(TABLESIZE-1)/outer);
  /* we make a tiny change to make sure I0 and I1 are not coincident, which
     simplifies tests later. This case will rarely happen anyway. */
  if (I0==I1) I1++;
  
  /* These indexes may be out of range. If so, lets do the right thing. */
  if (I0>=TABLESIZE) return 0.0; /* we're completely outside the rings,
				    no effect */
  if (I1>=TABLESIZE)
    {
      /* our outer range has run "off" of the end of the result. we have to
	 take this into account by figuring the ratio of what's fallen off
	 to what remains, and make sure our final output is properly
	 weighted */
      
      weight=((double)(TABLESIZE-I0))/(I1-I0);
      /* now we change I1 to start within the range. The weight parameter
	 will compensate for the range change. */
      I1=TABLESIZE-1;            
    }
  else weight=1.0; /* we're fully within the rings, we'll use the full
		      weight. */
  
  /* we now want the average value between I0 and I1. This is the
     easy part!  */
  
  return weight*(table[I1]-table[I0])/(I1-I0);    
}


/* routine to build the summed color table itself. This includes the
   computation of the ring density tucked within a simple loop to
   build the sum table as it goes. This is a precomputation that only
   needs to be done once. */

static MakeColorTable(double inner, double outer, double transition,
		      double varyscale, double ampratio,
		      double low, double high)
{
  double R, A, F;
  int i;
  
  /* sweep the radii out to the outer radius. Accumulate samples in
     the summed color table. Point sample the ring density at each sample.
     A radius of 0 is index 0. A radius of [outer] is equal to the table
     size-1.     
  */
  
  table[0]=0.0;  /* start with a 0 table */
  
  for (i=0; i<TABLESIZE-1; i++)
    {
      R=outer*(i+0.5)/(TABLESIZE-1);  /* R varies between 0 and [outer] */
      
      /* compute the simple inner/outer transitions to form an alpha channel
	 of a simple washer-like disk. This will be used to modulate the
	 density of the fine rings and prevent any rings from being too
	 close or too far from the center. */
      
      if (R<=inner) A=0.0;
      else /* in first transition zone? */
	if (R<inner+transition)
	  {
	    A=(R-inner)/transition;  /* linear 0 to 1 ramp */
	    A*=A*(3.0-2.0*A);        /* hermite curve smooth transition */
	  }
	else /* in outer transition zone? */
	  if (R>outer-transition)
	    {
	      A=(outer-R)/transition; /* linear 1 to 0 ramp */
	      A*=A*(3.0-2.0*A);        /* hermite curve smooth transition */
	    }
	  else A=1.0; /* we're in the main body of the ring */
      
      /* now let's compute the ring density. We use a 1D version of
	 fractal noise. We use a 3D noise routine but pass in (R, 0, 0)*/
      
      F=fractal3(R, 0.0, 0.0, varyscale, ampratio);
      
      /* F is now between -1 and 1. But we use the low and high values as
	 a clipping range for the noise. */
      
      if (F<=low) F=0.0; /* F is too low, we set the ring density to 0.0 */
      else if (F>=high) F=1.0; /* full ring density */
      else /* we're in a transition zone */       
	{
	  F=(F-low)/(high-low); /* now a 0 to 1 range */
	  F*=F*(3.0-2.0*F); /* Hermite it to make it smooth */
	}
      

      /* OK, our ring density is F and our shaping alpha value is A. The
	 net density is the PRODUCT of these two. */
      
      table[i+1]=table[i]+R*F;                        	      
   }    
  initialized=1;
  return;
}