#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create:
#	ebert_shar
#	musgrave
#	peachey_shar
#	perlin
# This archive created: Mon Sep  7 15:30:00 1992
# By:	David S. Ebert (The Ohio State University Dept of Computer & Information Science)
export PATH; PATH=/bin:/usr/bin:$PATH
echo shar: "extracting 'ebert_shar'" '(40562 characters)'
if test -f 'ebert_shar'
then
	echo shar: "will not over-write existing file 'ebert_shar'"
else
sed 's/^DE//' << \SHAR_EOF > 'ebert_shar'
DE#! /bin/sh
DE# This is a shell archive, meaning:
DE# 1. Remove everything above the #! /bin/sh line.
DE# 2. Save the resulting text in a file.
DE# 3. Execute the file with /bin/sh (not csh) to create:
DE#	ebert
DE# This archive created: Tue Sep  1 11:15:57 1992
DE# By:	David S. Ebert (The Ohio State University Dept of Computer & Information Science)
DEexport PATH; PATH=/bin:/usr/bin:$PATH
DEif test ! -d 'ebert'
DEthen
DE	echo shar: "creating directory 'ebert'"
DE	mkdir 'ebert'
DEfi
DEecho shar: "entering directory 'ebert'"
DEcd 'ebert'
DEecho shar: "extracting 'ease.c'" '(3751 characters)'
DEif test -f 'ease.c'
DEthen
DE	echo shar: "will not over-write existing file 'ease.c'"
DEelse
DEcat << \SHAR_EOF > 'ease.c'
DE
DE
DE/* EASE.C  ease-in/ease-out procedure for parameter                        */
DE/*                                                                         */
DE/* Author: Rick Parent, parent@cis.ohio-state.edu                          */
DE/*                                                                         */
DE/* ease-in is done by scaling the -PI/2 to 0 part of the sin curve to      */
DE/*         into the 0-t1 range.  e1=t1/(PI/2) is the scale factor, so the  */
DE/*         output of the sin function will be 0->2*t1/PI instead of 0->1   */
DE/* ease-out is handled similarly, scaling by e2 = (1-t2)/(PI/2).           */
DE/*         the range between t1 and t2 is mapped into itself with after    */
DE/*         translating by e1.                                              */
DE/* as the input value t goes from 0->1, the output value rt would go       */
DE/*         from 0->(e1 + (t2-t1) + e2).  Since we want the final result    */
DE/*         to still be in the 0->1 range, there is a final scale up to     */
DE/*         get it in this range  					   */
DE
DE/* given an input parametric value, t, it applies the ease function and    */
DE/* returns a new eased parameter, rt, to be used in linear interpolation   */
DE
DE#include  "math.h"
DE
DE/* ===================================================================== */
DE/* -------------------------------------------------------------------- */
DE/* ease-in/ease-out */
DE/* uses sin interpolation to get smooth in and smooth out function      */
DE/* abandoned for now. see ease() below */
DE/* -------------------------------------------------------------------- */
DEdouble easeA(t,t1,t2)
DEdouble  t,t1,t2;
DE{
DE  double  ans,s,e1,e2,nt,rt;
DE
DE  e1 = t1*2/M_PI;     /* scale sin curve: distance after t1 sin acceleration */
DE  e2 = (1-t2)*2/M_PI; /* scale sin curve: distance after t2 sin deceleration */
DE
DE  if (t<t1) {
DE    nt = t/t1;                         /* nt: linearly 0->1 */
DE    s = sin(- M_PI/2 + nt*M_PI/2) + 1; /* s accelerates from 0 to 1 */
DE    rt = s*e1;                         /* distance after s */
DE  }
DE  else if (t>t2) {
DE    nt = (t-t2)/(1.0-t2);                   /* t: linearly 0->1 */
DE    s = sin(nt*M_PI/2);                     /* s decelerates from 0 to 1 */
DE    rt = e1 + (t2-t1) + s*e2;               /* distance after s */
DE  }
DE  else {
DE    rt = e1 + t - t1;
DE  }
DE
DE  rt = rt/(e1+(t2-t1)+e2);   /* new total distance: scale back to 0-1 range */
DE  if (rt>1.0) rt = 1.0;
DE  return(rt);
DE}
DE
DE/* ===================================================================== */
DE/* -------------------------------------------------------------------- */
DE/* ease-in/ease-out */
DE/* using parabolic blending at the end points */
DE/* first leg has constant acceleration from 0 to v during time 0 to t1*/
DE/* second leg has constant velocity of v during time t1 to t2 */
DE/* third leg has constant deceleration from v to 0 during time t2 to 1 */
DE/* these are integrated to get the 'distance' traveled at any time */
DE/* -------------------------------------------------------------------- */
DEdouble ease(t,t1,t2)
DEdouble  t,t1,t2;
DE{
DE  double  ans,s,a,b,c,nt,rt;
DE  double  v,a1,a2;
DE
DE  v = 2/(1+t2-t1);  /* constant velocity attained */
DE  a1 = v/t1;        /* acceleration of first leg */
DE  a2 = -v/(1-t2);   /* deceleration of last leg */
DE
DE  if (t<t1) {
DE    rt = 0.5*a1*t*t;       /* pos = 1/2 * acc * t*t */
DE  }
DE  else if (t<t2) {
DE    a = 0.5*a1*t1*t1;      /* distance from first leg */
DE    b = v*(t-t1);            /* distance = vel * time  of second leg */
DE    rt = a + b;
DE  }
DE  else {
DE    a = 0.5*a1*t1*t1;      /* distance from first leg */
DE    b = v*(t2-t1);           /* distance from second leg */
DE    c = ((v + v + (t-t2)*a2)/2) * (t-t2);  /* distance = ave vel. * time */
DE    rt = a + b + c;
DE  }
DE  return(rt);
DE}
DE
DE
DE
DESHAR_EOF
DEfi
DEecho shar: "extracting 'functions.c'" '(24127 characters)'
DEif test -f 'functions.c'
DEthen
DE	echo shar: "will not over-write existing file 'functions.c'"
DEelse
DEcat << \SHAR_EOF > 'functions.c'
DE/*
DE *****************************************************************************
DE * This file contains functions explained in the SIGGRAPH'92 Course 23 notes *
DE * This C source file, as well as the .h files should compile and work       *
DE * correctly. If you have any problems, please contact the author.           *
DE *****************************************************************************
DE * Author: David S. Ebert, ebert@cis.ohio-state.edu                          *
DE * Date:   July, 1992                                                        *
DE *****************************************************************************
DE * This software was developed as part of on going research by the author in *
DE * the CIS department at The Ohio State University.                          *
DE * Any derivative of this code must credit the original author, The Ohio     *
DE * State University and must leave this paragraph unaltered.                 *
DE * All material is copyright 1992 by David S. Ebert.  Permission for         *
DE * non-profit use of the material is granted provided that all copyright     *
DE * notices are retained and the source of the material is appropriately      *
DE * acknowledged.                                                             *
DE *****************************************************************************
DE */
DE#include <math.h>
DE#include "graphics.h"
DE#include "vol.h"
DE#include "macro.h"
DE
DE#define COS_ERP(x)     ((cos((x)*M_PI)+1.0)*0.5)
DE#define POW_TABLE_SIZE 10000
DE#define RAMP_SIZE      200
DE#define OFFSET_SIZE    1024
DE#define DOWN_AMOUNT    0.011904762
DE#define SWIRL_AMOUNT   0.04986655            /* 2*PI/126  */
DE#define SWIRL_FRAMES   126                   /* should actually be 126 */
DE
DE
DEfloat offset[OFFSET_SIZE];
DE
DErgb_td marble();
DErgb_td marble_color();
DEdouble ease();
DEvoid   transform_XYZ();
DE
DE/*
DE *********************************************************************
DE *                            Basic_gas                              *
DE *********************************************************************
DE *      A simple gas function based on Perlin's turbulent flow       *
DE * Created for SIGGRAPH 92 Course 23.                                *
DE *********************************************************************
DE */
DE
DEbasic_gas(pnt,density,parms)
DE     xyz_td  pnt; 
DE     float   *density,*parms;
DE{
DE  float turb;
DE  int   i;
DE  static float pow_table[POW_TABLE_SIZE];
DE  static int calcd=1;
DE  
DE  if(calcd)
DE    { calcd=0;
DE      for(i=POW_TABLE_SIZE-1; i>=0; i--)
DE	pow_table[i] = (float)pow(((double)(i))/(POW_TABLE_SIZE-1)*
DE				  parms[1]*2.0,(double)parms[2]);
DE    }
DE  turb =fast_turbulence(pnt);
DE  *density = pow_table[(int)(turb*(.5*(POW_TABLE_SIZE-1)))];
DE}
DE
DE/*
DE *********************************************************************
DE *                            Basic_gas_sin                          *
DE *********************************************************************
DE *      A simple variation of the above basic_gas, using the sin     *
DE * function. Created for SIGGRAPH 92 Course 23.                      *
DE *********************************************************************
DE */
DE
DEbasic_gas_sin(pnt,density,parms)
DE     xyz_td  pnt; 
DE     float   *density,*parms;
DE{
DE  float turb;
DE  int   i;
DE  static float pow_table[POW_TABLE_SIZE];
DE  static int calcd=1;
DE  
DE  if(calcd)
DE    { calcd=0;
DE      for(i=POW_TABLE_SIZE-1; i>=0; i--)
DE	pow_table[i] = (float)pow(((double)(i))/(POW_TABLE_SIZE-1)*
DE				  parms[1]*2.0,(double)parms[2]);
DE    }
DE  turb =(1.0 +sin(fast_turbulence(pnt)*M_PI*5))*.5;
DE  *density = pow_table[(int)(turb*(.5*(POW_TABLE_SIZE-1)))];
DE}
DE
DE
DE/*
DE *********************************************************************
DE *                            Steam_slab1                            *
DE *********************************************************************
DE *      Based on basic gas,but has spherical ramp off of the density *
DE * for SIGGRAPH 92 Course  23.                                       *
DE *********************************************************************
DE * 04/08/92                                                          *
DE *********************************************************************
DE */
DE
DE
DEsteam_slab1(pnt, pnt_world, density,parms, vol)
DE     xyz_td  pnt, pnt_world;
DE     float   *density,*parms;
DE     vol_td  vol;
DE{
DE  float        turb, dist_sq, density_max;
DE  int          i, indx;
DE  xyz_td       distance;
DE  static float pow_table[POW_TABLE_SIZE], ramp[RAMP_SIZE];
DE  static int   calcd=1;
DE  
DE  if(calcd)
DE    { calcd=0;
DE      for(i=POW_TABLE_SIZE-1; i>=0; i--)
DE	pow_table[i] = (float)pow(((double)(i))/(POW_TABLE_SIZE-1)*
DE				  parms[1]*2.0,(double)parms[2]);
DE      make_ramp_table(ramp);
DE    }
DE  turb =fast_turbulence(pnt);
DE  *density = pow_table[(int)(turb*0.5*(POW_TABLE_SIZE-1))];
DE  
DE  /* determine distance from center of the slab ^2. */
DE  XYZ_SUB(distance,vol.shape.center, pnt_world);
DE  dist_sq = DOT_XYZ(distance,distance);
DE  density_max = dist_sq*vol.shape.inv_rad_sq.y; 
DE  indx = (int)((pnt.x+pnt.y+pnt.z)*100) & (OFFSET_SIZE -1);
DE  density_max += parms[3]*offset[indx];
DE  
DE  if(density_max >= .25) /* ramp off if > 25% from center */
DE    { i = (density_max -.25)*4/3*RAMP_SIZE; /* get table index 0:RAMP_SIZE-1 */
DE      i=MIN(i,RAMP_SIZE-1);
DE      density_max = ramp[i];
DE      *density *=density_max;
DE    }
DE}
DE
DE
DE/*
DE *********************************************************************
DE *                            Steam_slab_final                       *
DE *********************************************************************
DE *      Based on basic gas,but has spherical ramp off of the density.*
DE * Finally, it has vertical ramp off of the density.                 *
DE * Created for SIGGRAPH 92 Course # 23.                              *
DE *********************************************************************
DE */
DE
DE
DEsteam_slab_final(pnt, pnt_world, density,parms, vol)
DE     xyz_td  pnt, pnt_world;
DE     float   *density,*parms;
DE     vol_td  vol;
DE{
DE  float        turb, dist_sq, density_max, dist, offset2;
DE  int          i, indx;
DE  xyz_td       distance;
DE  static float pow_table[POW_TABLE_SIZE], ramp[RAMP_SIZE];
DE  static int   calcd=1;
DE  
DE  if(calcd)
DE    { calcd=0;
DE      for(i=POW_TABLE_SIZE-1; i>=0; i--)
DE	pow_table[i] = (float)pow(((double)(i))/(POW_TABLE_SIZE-1)*
DE				  parms[1]*2.0,(double)parms[2]);
DE      make_ramp_table(ramp);
DE    }
DE  turb =fast_turbulence(pnt);
DE  *density = pow_table[(int)(turb*0.5*(POW_TABLE_SIZE-1))];
DE  
DE  /* determine distance from center of the slab ^2. */
DE  XYZ_SUB(distance,vol.shape.center, pnt_world);
DE  dist_sq = DOT_XYZ(distance,distance);
DE  density_max = dist_sq*vol.shape.inv_rad_sq.y; 
DE  indx = (int)((pnt.x+pnt.y+pnt.z)*100) & (OFFSET_SIZE -1);
DE  density_max += parms[3]*offset[indx];
DE  
DE  if(density_max >= .25) /* ramp off if > 25% from center */
DE    { i = (density_max -.25)*4/3*RAMP_SIZE; /* get table index 0:RAMP_SIZE-1 */
DE      i=MIN(i,RAMP_SIZE-1);
DE      density_max = ramp[i];
DE      *density *=density_max;
DE    }
DE
DE  /* Vertical ramp off of the density
DE   */
DE  dist = pnt_world.y - vol.shape.center.y;
DE  if(dist > 0.0)
DE    { dist = (dist +offset[indx]*.1)*vol.shape.inv_rad.y;
DE      if(dist > .05)
DE	{ offset2 = (dist -.05)*1.111111;
DE	  offset2 = 1 - (exp(offset2)-1.0)/1.718282;
DE	  offset2 *=parms[1];
DE	  *density *= offset2;
DE	}
DE    }
DE
DE}
DE
DE/*
DE **********************************************************************
DE *                           MARBLE                                   *
DE **********************************************************************
DE * These marble functions were used to create images in the           *
DE * SIGGRAPH'92 course 23 notes.                                       *
DE **********************************************************************
DE*/
DErgb_td marble(pnt)
DE     xyz_td  pnt;
DE{
DE  float y;
DE  y = pnt.y + 3.0*fast_turbulence(pnt);
DE  y = sin(y*M_PI); 
DE  return (marble_color(y));
DE}
DE
DE
DErgb_td marble_color(x)
DE     float x;
DE{
DE  rgb_td  clr;
DE  
DE  x = sqrt(x+1.0)*.7071; 
DE  clr.g = .30 + .8*x;
DE  x=sqrt(x);
DE  clr.r = .30 + .6*x;
DE  clr.b = .60 + .4*x; 
DE  return (clr);
DE}
DE
DE
DE
DE/*
DE **********************************************************************
DE *                           MARBLE_FORMING                           *
DE **********************************************************************
DE * These marble functions were used to create images in the           *
DE * SIGGRAPH'92 course 23 notes. They simulate the formation of marble *
DE **********************************************************************
DE*/
DErgb_td marble_forming(pnt, frame_num, start_frame, end_frame)
DE     xyz_td  pnt;
DE     int     frame_num, start_frame, end_frame;
DE{
DE  float x, turb_percent, displacement;
DE  
DE  if(frame_num < start_frame)
DE    { turb_percent=0;
DE      displacement=0;
DE    }
DE  else if (frame_num >= end_frame)
DE    { turb_percent=1;
DE      displacement= 3;
DE    }
DE  else
DE    { turb_percent= ((float)(frame_num-start_frame))/ (end_frame-start_frame);
DE      displacement = 3*turb_percent;
DE    }
DE  
DE  x = pnt.x + turb_percent*3.0*fast_turbulence(pnt) - displacement;
DE  x = sin(x*M_PI); 
DE  return (marble_color(x));
DE}
DE
DE
DErgb_td marble_forming2(pnt, frame_num, start_frame, end_frame, heat_length)
DE     xyz_td  pnt;
DE     int     frame_num, start_frame, end_frame, heat_length;
DE{
DE  float       x, turb_percent, displacement, glow_percent;
DE  rgb_td      m_color;
DE  
DE  if(frame_num < (start_frame-heat_length/2) || 
DE     frame_num > end_frame+heat_length/2)
DE    glow_percent=0;
DE  else if (frame_num < start_frame + heat_length/2)
DE    glow_percent= 1.0 - 
DE      ease( ((start_frame+heat_length/2-frame_num)/ heat_length),0.4,0.6); 
DE  else if (frame_num > end_frame-heat_length/2)
DE    glow_percent =  
DE      ease( ((frame_num-(end_frame-heat_length/2))/ heat_length),0.4,0.6);
DE  else 
DE    glow_percent=1.0;
DE  
DE  if(frame_num < start_frame)
DE    { turb_percent=0;
DE      displacement=0;
DE    }
DE  else if (frame_num >= end_frame)
DE    { turb_percent=1;
DE      displacement= 3;
DE    }
DE  else
DE    { turb_percent= ((float)(frame_num-start_frame))/(end_frame-start_frame);
DE      turb_percent=ease(turb_percent, 0.3, 0.7);
DE      displacement = 3*turb_percent;
DE    }
DE  
DE  x = pnt.y + turb_percent*3.0*fast_turbulence(pnt) - displacement;
DE  x = sin(x*M_PI); 
DE  m_color=marble_color(x);
DE  glow_percent= .5* glow_percent;
DE  m_color.r= glow_percent*(1.0)+ (1-glow_percent)*m_color.r;	
DE  m_color.g= glow_percent*(0.4)+ (1-glow_percent)*m_color.g;	
DE  m_color.b= glow_percent*(0.8)+ (1-glow_percent)*m_color.b;	
DE  return(m_color);
DE}
DE
DE
DE/*
DE **********************************************************************
DE *                           MOVING_MARBLE                            *
DE **********************************************************************
DE * This marble functions uses a helical path for motion of the point  *
DE * before evaluation of the marble solid texture function. This       *
DE * function was developed for SIGGRAPH'92 course 23 notes.            *
DE **********************************************************************
DE*/
DE
DE
DE#define RAD1 .1
DE#define RAD2 .14  
DE
DErgb_td moving_marble(pnt, frame_num)
DE     xyz_td  pnt;
DE     int     frame_num;
DE{
DE  float        x, tmp, tmp2;
DE  static float down, theta, sin_theta, cos_theta;
DE  xyz_td       hel_path, direction;
DE  static int   calcd=1;
DE
DE  if(calcd)
DE    { theta =(frame_num%SWIRL_FRAMES)*SWIRL_AMOUNT; /* swirling effect */
DE      cos_theta = RAD1 * cos(theta) + 0.5;
DE      sin_theta = RAD2 * sin(theta) - 2.0;
DE      down = (float)frame_num*DOWN_AMOUNT+2.0;
DE      calcd=0;
DE    }
DE  tmp = fast_noise(pnt); /* add some randomness */
DE  tmp2 = tmp*1.75;
DE  
DE  /* calculate the helical path */
DE  hel_path.y = cos_theta  + tmp;
DE  hel_path.x = (- down)   + tmp2;
DE  hel_path.z = sin_theta  - tmp2;
DE  XYZ_ADD(direction, pnt, hel_path);
DE  
DE  x = pnt.y + 3.0*fast_turbulence(direction);
DE  x = sin(x*M_PI); 
DE  return (marble_color(x));
DE}
DE
DE
DE
DE/*
DE **********************************************************************
DE *                                     FOG                            *
DE **********************************************************************
DE * This function can be used for solid textured transparency based    *
DE * fog. This is from the SIGGRAPH92 Course 23 notes.                  *
DE **********************************************************************
DE */
DE
DEvoid fog(pnt, transp, frame_num)
DE     xyz_td  pnt;
DE     float   *transp;
DE     int      frame_num;
DE{
DE  float tmp;
DE  xyz_td direction,cyl;
DE  double theta;
DE  
DE  pnt.x += 2.0 + fast_turbulence(pnt);
DE  tmp = fast_noise(pnt);
DE  pnt.y +=  4+tmp;
DE  pnt.z += -2 - tmp; 
DE  
DE  theta =(frame_num%SWIRL_FRAMES)*SWIRL_AMOUNT; 
DE  cyl.x =RAD1 * cos(theta); 
DE  cyl.z =RAD2 * sin(theta);
DE  
DE  direction.x = pnt.x + cyl.x;
DE  direction.y = pnt.y - frame_num*DOWN_AMOUNT;
DE  direction.z = pnt.z + cyl.z;
DE  
DE  *transp = fast_turbulence(direction);
DE  *transp = (1.0 -(*transp)*(*transp)*.275);
DE  *transp =(*transp)*(*transp)*(*transp);
DE}
DE
DE
DE/*
DE **********************************************************************
DE *                           STEAM MOVING                             *
DE **********************************************************************
DE * This function creates steam rising from teacup or other circularly *
DE * cylindrical object. This is from the SIGGRAPH92 Course 23 notes.   *
DE **********************************************************************
DE */
DE
DE
DEsteam_moving(pnt, pnt_world, density,parms, vol, frame_num)
DE     xyz_td  pnt, pnt_world;
DE     float   *density,*parms;
DE     vol_td  vol;
DE     int     frame_num;
DE{
DE  
DE  float tmp,turb, dist_sq, density_max, offset2, theta, dist;
DE  static float ramp[RAMP_SIZE];
DE  xyz_td direction, diff;
DE  int i, indx;
DE  static float pow_table[POW_TABLE_SIZE];
DE  static int calcd=1;
DE  static float down, cos_theta2, sin_theta2;
DE  
DE  if(calcd)
DE    { calcd=0;
DE      /* determine how to move the point through the space (helical path) */
DE      theta =(frame_num%SWIRL_FRAMES)*SWIRL_AMOUNT;       
DE      down = (float)frame_num*DOWN_AMOUNT*3.0 +4.0;
DE      cos_theta2 = RAD1*cos(theta) +2.0;
DE      sin_theta2 = RAD2*sin(theta) -2.0;
DE      
DE      for(i=POW_TABLE_SIZE-1; i>=0; i--)
DE	pow_table[i] = (float)pow(((double)(i))/(POW_TABLE_SIZE-1)*
DE				  parms[1]*2.0,(double)parms[2]);
DE      make_ramp_table(ramp);
DE    }
DE  
DE  tmp = fast_noise(pnt);
DE  direction.x = pnt.x + cos_theta2 +tmp;
DE  direction.y = pnt.y - down + tmp;
DE  direction.z = pnt.z +sin_theta2 +tmp;
DE  
DE  turb =fast_turbulence(direction);
DE  *density = pow_table[(int)(turb*0.5*(POW_TABLE_SIZE-1))];
DE  
DE  /* determine distance from center of the slab ^2. */
DE  XYZ_SUB(diff,vol.shape.center, pnt_world);
DE  dist_sq = DOT_XYZ(diff,diff);
DE  density_max = dist_sq*vol.shape.inv_rad_sq.y; 
DE  indx = (int)((pnt.x+pnt.y+pnt.z)*100) & (OFFSET_SIZE -1);
DE  density_max += parms[3]*offset[indx];
DE  
DE  if(density_max >= .25) /* ramp off if > 25% from center */
DE    { i = (density_max -.25)*4/3*RAMP_SIZE; /* get table index 0:RAMP_SIZE-1 */
DE      i=MIN(i,RAMP_SIZE-1);
DE      density_max = ramp[i];
DE      *density *=density_max;
DE    }
DE  
DE  /* ramp it off vertically */
DE  dist = pnt_world.y - vol.shape.center.y;
DE  if(dist > 0.0)
DE    { dist = (dist +offset[indx]*.1)*vol.shape.inv_rad.y;
DE      if(dist > .05)
DE	{ offset2 = (dist -.05)*1.111111;
DE	  offset2 = 1 - (exp(offset2)-1.0)/1.718282;
DE	  offset2*=parms[1];
DE	  *density *= offset2;
DE	}
DE    }
DE}
DE
DE
DE
DE/*
DE **********************************************************************
DE *                    SPHERICAL_ATTRACTOR                             *
DE **********************************************************************
DE * This function is a simple spherical attractor. It is an example of *
DE * a functional flow field primitive. This is from the SIGGRAPH92     *
DE * Course 23 notes.                                                   *
DE **********************************************************************
DE */
DE
DE
DEspherical_attractor(point, FF, direction, density_scaling, velocity, 
DE		    percent_to_use)
DE     xyz_td       point, *direction;
DE     flow_func_td FF;
DE     float        *density_scaling, *velocity, *percent_to_use;
DE{
DE  float        dist, d2;
DE  
DE  /*calculate distance and direction from center of attractor  */
DE  XYZ_SUB(*direction, point, FF.center);
DE  dist=sqrt(DOT_XYZ(*direction,*direction));
DE  
DE  /* set the density scaling and the velocity to 1 */
DE  *density_scaling=1.0;
DE  *velocity=1.0;
DE  
DE  /* calculate the falloff factor (cosine) */
DE  if(dist > FF.distance)
DE    *percent_to_use=0;
DE  else if (dist < FF.falloff_start)
DE    *percent_to_use=1.0;
DE  else
DE    { d2 = (dist - FF.falloff_start)/(FF.distance - FF.falloff_start);
DE      *percent_to_use = (cos(d2*M_PI)+1.0)*.5;
DE    }
DE}
DE
DE
DE/*
DE **********************************************************************
DE *                         CALC VORTEX                                *
DE **********************************************************************
DE * This function is a visual simultion of vortex like motion.         *
DE * It is an example of a functional flow field primitive.             *
DE * This is from the SIGGRAPH92 23 notes.                              *
DE **********************************************************************
DE */
DE
DEcalc_vortex(pt, ff, path,percent_to_use, frame_num, velocity)
DE     xyz_td       *pt, *path;
DE     flow_func_td *ff;
DE     float        *percent_to_use, *velocity;
DE     int          frame_num;
DE{
DE  static tran_mat_td  mat={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
DE  xyz_td              dir, pt2, diff;
DE  float               theta, dist, d2, dist2;
DE  float               ratio_mult;
DE  
DE  /*calculate distance from center of vortex */
DE  XYZ_SUB(diff,(*pt), ff->center);
DE  dist=sqrt(DOT_XYZ(diff,diff));
DE  dist2 = dist/ff->distance;
DE  /* calculate angle of rotation about the axis */
DE  theta = (ff->parms[0] *(1+.001*(frame_num)))/(pow((.1+dist2*.9),ff->parms[1])); 
DE  
DE  /* calculate the matrix for rotating about the cylinder's axis */
DE  calc_rot_mat(theta, ff->axis, mat, *ff);
DE  transform_XYZ((long)1,mat,pt,&pt2);
DE  XYZ_SUB(dir,pt2,(*pt));
DE  path->x = dir.x;
DE  path->y = dir.y;
DE  path->z = dir.z;
DE  
DE  /* Have the maximum strength increase from frame parms[4] to 
DE   * parms[5]to a maximum of parms[2]   */
DE  if(frame_num < ff->parms[4])
DE    ratio_mult=0;
DE  else if (frame_num <= ff->parms[5])
DE    ratio_mult = (frame_num - ff->parms[4])/(ff->parms[5] - ff->parms[4])* ff->parms[2];
DE  else
DE    ratio_mult = ff->parms[2];
DE  
DE  /* calculate the falloff factor */
DE  if(dist > ff->distance)
DE    { *percent_to_use=0;
DE      *velocity=1;
DE    }
DE  else if (dist < ff->falloff_start)
DE    { *percent_to_use=1.0 *ratio_mult;
DE      /*calc velocity  */
DE      *velocity= 1.0+(1.0 - (dist/ff->falloff_start));
DE    }
DE  else
DE    { d2 = (dist - ff->falloff_start)/(ff->distance - ff->falloff_start);
DE      *percent_to_use = (cos(d2*M_PI)+1.0)*.5*ratio_mult;
DE    }
DE}
DE
DE
DE
DE/*
DE **********************************************************************
DE *                         STEAM WIND                                 *
DE **********************************************************************
DE * This function simulates a breeze blowing from the left of the      *
DE * screen. It is a flow filed function and needs to be called to move *
DE * the point before the evaluation of the turbulence function in the  *
DE * volume density functions (steam moving).                           *
DE * It is an example of a functional flow field primitive.             *
DE * This is from the SIGGRAPH92 23 notes.                              *
DE **********************************************************************
DE */
DE 
DE      
DE  
DE
DEfloat
DEsteam_wind(point, FF, direction, density_scaling, velocity, 
DE		     percent_to_use, frame_num)
DE     xyz_td       point, *direction;
DE     flow_func_td FF;
DE     float        *density_scaling, *velocity, *percent_to_use;
DE     int          frame_num;
DE{
DE  float        dist, d2, ratio_mult, x_disp;
DE  
DE
DE  /*calculate distance and direction from center of attractor  */
DE  XYZ_SUB(*direction, point, FF.center);
DE  dist=sqrt(DOT_XYZ(*direction,*direction));
DE  
DE  /* set the density scaling and the velocity to 1 */
DE  *density_scaling=1.0;
DE  *velocity=1.0;
DE  
DE  /* calculate the falloff factor (cosine) */
DE  if(dist > FF.distance)
DE    *percent_to_use=0;
DE  else if (dist < FF.falloff_start)
DE    *percent_to_use=1.0;
DE  else
DE    { d2 = (dist - FF.falloff_start)/(FF.distance - FF.falloff_start);
DE      *percent_to_use = (cos(d2*M_PI)+1.0)*.5;
DE    }
DE
DE  /*
DE   ***************************************************************************
DE   * Move the Volume of the Steam 
DE   ***************************************************************************
DE   * the shifting is based on the height above the cup (parms[13]->parms[14])
DE   * and the frame range for increasing the strength of the attractor.
DE   *  This is gotten from ratio_mult that is calculated below. 
DE   ***************************************************************************
DE   */
DE  
DE  /* Have the maximum strength increase from frame parms[4] to 
DE   * parms[5]to a maximum of parms[2]. The use of the ease function has been
DE   * added to the function given in the course notes.
DE   */
DE  if(frame_num < FF.parms[4])
DE    ratio_mult=0;
DE  else if (frame_num <= FF.parms[5])
DE    ratio_mult = ease((double)((frame_num - FF.parms[4])/
DE			       (FF.parms[5]-FF.parms[4])), 0.3, 0.7)* 
DE				 FF.parms[2];
DE  else
DE    ratio_mult=1.0;
DE
DE  /* calculate the displacement of the steam volume */
DE  if(point.y < FF.parms[6])
DE    x_disp=0;
DE  else
DE    { if(point.y <= FF.parms[7])
DE	d2 =COS_ERP((point.y - FF.parms[6])/(FF.parms[7] -FF.parms[6]));
DE    else 
DE      d2=0;
DE      x_disp = (1-d2)*ratio_mult*FF.parms[8]+fast_noise(point)*FF.parms[9];
DE    }
DE  return(x_disp);
DE
DE}
DE
DE
DE
DE
DE
DE
DEmake_ramp_table(ramp)
DE     float *ramp;
DE{
DE  int   i;
DE  float dist;
DE  for(i = 0; i < RAMP_SIZE; i++)
DE    { dist =i/(RAMP_SIZE -1.0);
DE      ramp[i]=(cos(dist*M_PI) +1.0)/2.0;
DE    }
DE}
DE
DE
DE
DE
DE
DE
DE
DE/* transform_xyz: transform an array of xyz_td and output cooresponding
DE * xyz_td pnts which are transformed by the matrix input. Note: this routine
DE * divides through by w just in case w is not 1. It should actually only be 
DE * be used for scale, rotate, and translate.
DE */
DE 
DEvoid
DEtransform_XYZ(num_pts, mat, inpts, outpts)
DE
DE    long        num_pts;
DE    tran_mat_td mat;
DE    xyz_td      *inpts;
DE    xyz_td     *outpts;
DE    
DE{
DE    long i;
DE    xyz_td *tmp, tmp2;
DE
DE    if(num_pts==1)
DE      {
DE	tmp2.x = mat[0][0]*inpts[0].x + mat[1][0]*inpts[0].y
DE	         + mat[2][0]*inpts[0].z + mat[3][0];
DE	tmp2.y = mat[0][1]*inpts[0].x + mat[1][1]*inpts[0].y
DE                 + mat[2][1]*inpts[0].z + mat[3][1];
DE	tmp2.z = mat[0][2]*inpts[0].x + mat[1][2]*inpts[0].y
DE                 + mat[2][2]*inpts[0].z + mat[3][2];
DE	*outpts=tmp2;
DE	return;
DE      }
DE    /* else */
DE    tmp =(xyz_td *)malloc(sizeof(xyz_td)*num_pts);
DE
DE    for (i = 0; i < num_pts; i++)
DE      {
DE	tmp[i].x = mat[0][0]*inpts[i].x + mat[1][0]*inpts[i].y
DE	            + mat[2][0]*inpts[i].z + mat[3][0];
DE	tmp[i].y = mat[0][1]*inpts[i].x + mat[1][1]*inpts[i].y
DE	            + mat[2][1]*inpts[i].z + mat[3][1];
DE	tmp[i].z = mat[0][2]*inpts[i].x + mat[1][2]*inpts[i].y
DE                    + mat[2][2]*inpts[i].z + mat[3][2];
DE	outpts[i].x = tmp[i].x;
DE	outpts[i].y = tmp[i].y;
DE	outpts[i].z = tmp[i].z;
DE      }
DE    
DE}
DE
DE
DE
DE
DE
DEcalc_rot_mat(theta, axis, mat, ff)
DE     float        theta;
DE     xyz_td       axis;
DE     tran_mat_td  mat;
DE     flow_func_td ff;
DE{
DE  float cos_theta, sin_theta, compl_cos;
DE
DE  /*
DE   * calculate the matrix for rotating about the cylinder's axis
DE   * by theta. From Roger's and Adam's Mathematical Elements for 
DE   * Computer graphics pp 55-56 & 207,208
DE   */
DE  cos_theta = cos(theta);
DE  sin_theta = sin(theta);
DE  compl_cos = 1-cos_theta;
DE
DE  mat[0][0]= ff.dcx_2  + (ff.compl_dcx_2)*cos_theta;
DE  mat[1][1]= ff.dcy_2  + (ff.compl_dcy_2)*cos_theta;
DE  mat[2][2]= ff.dcz_2  + (ff.compl_dcz_2)*cos_theta;
DE  mat[0][1]= ff.dcx_dcy *(compl_cos) + ff.dcz*sin_theta;
DE  mat[0][2]= ff.dcx_dcz *(compl_cos) - ff.dcy*sin_theta; 
DE  mat[1][0]= ff.dcx_dcy *(compl_cos) - ff.dcz*sin_theta;
DE  mat[1][2]= ff.dcy_dcz *(compl_cos) + ff.dcx*sin_theta;
DE  mat[2][0]= ff.dcx_dcz *(compl_cos) + ff.dcy*sin_theta;
DE  mat[2][1]= ff.dcy_dcz *(compl_cos) - ff.dcx*sin_theta;
DE  mat[3][3]=1;
DE}
DE
DESHAR_EOF
DEfi
DEecho shar: "extracting 'README'" '(1472 characters)'
DEif test -f 'README'
DEthen
DE	echo shar: "will not over-write existing file 'README'"
DEelse
DEcat << \SHAR_EOF > 'README'
DE      SIGGRAPH '92 COURSE 23: PROCEDURAL MODELING AND RENDERING TECHNIQUES
DE
DE                 Modeling and Animating Gases and Fluids 
DE
DE                        by David S. Ebert.
DE
DEThis directory contains the C functions described in my section of the above
DEcourse.  All files necessary to compile the routines are included in
DEthis directory. The only files not included are my implementations of
DEnoise and turbulence. You can use the implementations provided by Ken
DEPerlin or Darwyn Peachey. My implementations are faster, not as
DEelegant, and probably produce more grid artifacts.
DE
DEIf you have any problems with the software, please contact me. If you
DEhave any enhancements, or extensions, I would be interested in hearing
DEabout them.
DE
DEContact Information:
DE
DE--  Dr. David S. Ebert,   Department of Computer and Information Science  -----
DE---  The Ohio State University; 2036 Neil Ave. Columbus OH USA 43210-1277  ----
DE-- ebert@cis.ohio-state.edu or ..!{att,pyramid,killer}!cis.ohio-state.edu!ebert
DE
DE
DEThis software was developed as part of on going research by the author
DEin the CIS department at The Ohio State University.  Any derivative of
DEthis code must credit the original author, The Ohio State University
DEand must leave this paragraph unaltered.  All material is copyright
DE1992 by David S. Ebert.  Permission for non-profit use of the material
DEis granted provided that all copyright notices are retained and the
DEsource of the material is appropriately acknowledged.
DE
DESHAR_EOF
DEfi
DEecho shar: "extracting 'graphics.h'" '(1404 characters)'
DEif test -f 'graphics.h'
DEthen
DE	echo shar: "will not over-write existing file 'graphics.h'"
DEelse
DEcat << \SHAR_EOF > 'graphics.h'
DE/* This is an extremely edited version of the standard graphics.h file from 
DE * the EDGE software library at the Department of Computer and Infomation 
DE * Science,The Ohio State University.
DE */
DE
DE/* EDGE standard header file for graphics programs */
DE#ifndef GRAPHICS_H
DE#define GRAPHICS_H
DE
DE/* define unsigned types for convenience */
DE/* can only define once or cc complains, so need a flag */
DE
DEtypedef unsigned char	uchar;
DEtypedef unsigned long	ulong;
DE
DE#define NL	'\n'
DE#define PI	3.1415926536
DE#define DtoR	0.0174532925
DE#define RtoD	57.295779515
DE
DEtypedef struct point	{ float x, y, z; }	point;
DEtypedef struct color	{ float r,g,b; }	color;
DE
DEdouble fabs(), sqrt(), sin(), cos(), tan();
DE
DE#  define	TRUE	1
DE#  define	YES	1
DE#  define	FALSE	0
DE#  define	NO	0
DE
DE#  undef	MAX
DE#  define	MAX(a, b)	((a) > (b) ? (a) : (b))
DE#  undef	MIN
DE#  define	MIN(a, b)	((a) < (b) ? (a) : (b))
DE
DE#  define	ABS(a)		((a) > 0 ? (a) : -(a))
DE
DE
DE
DE
DE  typedef int bool;
DE
DE  typedef struct xyz_td
DE          {
DE              float    x,
DE                       y,
DE                       z;
DE	  } xyz_td;
DE
DE
DE
DE  typedef struct xyzw_td
DE          {
DE              float    x,
DE                       y,
DE                       z,
DE                       w;
DE	  } xyzw_td;
DE
DE
DE
DE  typedef struct rgb_td
DE          {
DE              float    r,
DE                       g,
DE                       b;
DE	  } rgb_td;
DE
DE
DE
DE  typedef float tran_mat_td[4][4];
DE
DE
DE
DE#endif GRAPHICS_H
DESHAR_EOF
DEfi
DEecho shar: "extracting 'vol.h'" '(4826 characters)'
DEif test -f 'vol.h'
DEthen
DE	echo shar: "will not over-write existing file 'vol.h'"
DEelse
DEcat << \SHAR_EOF > 'vol.h'
DE
DE/*
DE * This is an edited version of my volume structure definitions. The editing 
DE *  was done to aid in understanding the functions.
DE *  There is probably still some non relevant material in these header files.
DE *  Please ignore these. Good Luck.
DE *  David S. Ebert 07/92.
DE */
DE
DE#define FNAMELEN        256
DE#define MAXLIGHTS       100
DE/*
DE ***********************************
DE * volume geometry types 
DE ***********************************
DE */
DE#define SPHERE          0
DE#define BOX             1
DE#define TORUS           2
DE#define FUNCT           3  
DE#define CYLINDER        4
DE
DE/* another illumination type for gases */
DE#define GASEOUS 3
DE
DE
DE#ifdef cray
DE#  define shad_val_type  float
DE#else
DE#  define shad_val_type  unsigned short
DE#endif
DE
DEtypedef struct sphere_td
DE    {
DE      xyz_td   center;
DE      float    radius;
DE      float    rad_sq, inv_rad_sq;
DE    } sphere_td;
DE
DEtypedef struct cylinder_td
DE    {
DE      xyz_td   top, bot;
DE      float    radius;
DE      float    rad_sq, inv_rad_sq;
DE    } cylinder_td;
DE
DEtypedef struct slab_td
DE    {
DE      float    min[3];
DE      float    max[3];
DE      xyz_td   center, rad, rad_sq, inv_rad, inv_rad_sq;
DE    } slab_td;
DE
DEtypedef struct vol_shape_td
DE    {
DE      sphere_td      sphere;
DE      slab_td        b_box;        /* the Volumes defining parallelepiped   */
DE      xyz_td         scale;        /*ow to scale obj to get into -1:1 space */
DE      xyz_td         trans;        /*how to trans obj to get into -1:1 space*/
DE      xyz_td         center;
DE      xyz_td         inv_rad;
DE      xyz_td         inv_rad_sq;
DE      int            type;
DE      char           name[FNAMELEN];
DE      float          delta_vox, delta_dir,
DE                     pow_c, shadow_delta;
DE      float          parm[10];
DE    } vol_shape_td;
DE
DEtypedef struct shadow_td
DE    {
DE      xyz_td        step_size, inv_step, start;
DE      shad_val_type *shadow_table[MAXLIGHTS];
DE      int            x_tab_2, y_tab_2, z_tab_2; /* log2 of table size for bit
DE						   shifting */
DE      int            x_table_size, y_table_size, z_table_size;
DE    } shadow_td;
DE
DE/*
DE *****************************************************************************
DE *                  FUNCTION FIELD DEFINITIONS
DE *****************************************************************************
DE */
DE
DEtypedef unsigned char flow_func_type;
DE
DEtypedef struct flow_func_td       
DE    {
DE      short         func_type;
DE      xyz_td        center, axis;
DE      float         distance;
DE      float         falloff_start;
DE      short         falloff_type;
DE      float         parms[20];
DE      char          filename[FNAMELEN];
DE      float         angle;
DE      long          angle_falloff;
DE      double        cos_angle;
DE      xyz_td        angle_axis;
DE      /* parameters precomputed to rotate about the axis dcx = direction 
DE       * cosine for rotation about the x axis, etc */
DE      float         dcx,dcy,dcz,dcx_2,dcy_2,dcz_2,compl_dcx_2,compl_dcy_2,
DE                    compl_dcz_2,dcx_dcy,dcx_dcz,dcy_dcz;
DE      tran_mat_td   mat;
DE    } flow_func_td;
DE
DE
DE
DE/*
DE *****************************************************************************
DE *               FLOW FIELD DEFINITIONS
DE *****************************************************************************
DE */
DE#ifdef cray
DE   typedef float v_field_type[5];
DE#  define FLOW_MAX  1.0
DE#  define INV_FLOW_MAX   1.0
DE#else
DE   typedef   short v_field_type[5];
DE#  define FLOW_MAX  32767
DE#  define INV_FLOW_MAX   .000030517578125
DE#endif 
DE
DE
DE/*
DE *****************************************************************************
DE *               VOLUME STRUCTURE DEFINITIONS
DE *****************************************************************************
DE */
DE
DEtypedef struct vol_td
DE    {
DE      short          obj_type;
DE      tran_mat_td    tm,           /* transformation mat for the object  */
DE                     inv_tm;
DE      float          y_obj_min,    /* min & max y per volume             */
DE	             y_obj_max;
DE      float          x_obj_min,    /* min & max x per volume             */
DE	             x_obj_max;
DE      int            shape_type;   /* 0=sphere, 1=box;                   */
DE      int            self_shadow;  /* 0=no, 1=yes,2=table based          */
DE      xyz_td         scale_obj;   /*how to scale obj to get into -1:1 space*/
DE      xyz_td         trans_obj;   /*how to trans obj to get into -1:1 space*/
DE      vol_shape_td   shape;
DE      int           (*funct_name)();
DE      rgb_td         color;
DE      float          amb_coef;
DE      float          diff_coef;
DE      float          spec_coef;
DE      float          spec_power;
DE      int            illum_type;  /* illumination - phong, blinn, cooke-t */
DE      int            color_type;  /* constant, solid*/
DE      float          indx_refrac;
DE
DE      /************* the 3D table information ***************/
DE      shadow_td      shadow;      
DE
DE    } vol_td;
DE	 
DE
DE
DE
DE
DE
DESHAR_EOF
DEfi
DEecho shar: "extracting 'macro.h'" '(3276 characters)'
DEif test -f 'macro.h'
DEthen
DE	echo shar: "will not over-write existing file 'macro.h'"
DEelse
DEcat << \SHAR_EOF > 'macro.h'
DE#define XYZ_ADD(c,a,b)   { (c).x = (a).x + (b).x; (c).y = (a).y + (b).y;\
DE                           (c).z = (a).z + (b).z;}
DE#define XYZ_SUB(c,a,b)   { (c).x = (a).x - (b).x; (c).y = (a).y - (b).y; \
DE                           (c).z = (a).z - (b).z;}
DE#define XYZ_MULT(c,a,b)  { (c).x = (a).x * (b).x; (c).y = (a).y * (b).y; \
DE                           (c).z = (a).z * (b).z;}
DE#define XYZ_DIV(c,a,b)   { (c).x = (a).x / (b).x; (c).y = (a).y / (b).y; \
DE                           (c).z = (a).z / (b).z;}
DE
DE#define XYZ_INC(c,a)     { (c).x += (a).x; (c).y += (a).y;\
DE                           (c).z += (a).z;}
DE#define XYZ_DEC(c,a)     { (c).x -= (a).x; (c).y -= (a).y;\
DE                           (c).z -= (a).z;}
DE#define XYZ_ADDC(c,a,b)  { (c).x = (a).x + (b); (c).y = (a).y + (b);\
DE                           (c).z = (a).z + (b);}
DE#define XYZ_SUBC(c,a,b)  { (c).x = (a).x - (b); (c).y = (a).y - (b); \
DE                           (c).z = (a).z - (b);}
DE#define XYZ_MULTC(c,a,b) { (c).x = (a).x * (b); (c).y = (a).y * (b);\
DE                           (c).z = (a).z * (b);}
DE#define XYZ_DIVC(c,a,b)  { (c).x = (a).x / (b); (c).y = (a).y / (b);\
DE                           (c).z = (a).z / (b);}
DE#define XYZ_COPY(b,a)    { (b).x = (a).x; (b).y = (a).y; (b).z = (a).z; }
DE#define XYZ_COPYC(b,a)   { (b).x = (a); (b).y = (a); (b).z = (a); }
DE#define CROSS_XYZ(c,a,b) { (c).x = (a).y * (b).z - (a).z * (b).y; \
DE			   (c).y = (a).z * (b).x - (a).x * (b).z; \
DE			   (c).z = (a).x * (b).y - (a).y * (b).x; }
DE#define CROSS_3(c,a,b)   { (c)[0] = (a)[1] * (b)[2] - (a)[2] * (b)[1]; \
DE			   (c)[1] = (a)[2] * (b)[0] - (a)[0] * (b)[2]; \
DE			   (c)[2] = (a)[0] * (b)[1] - (a)[1] * (b)[0]; }
DE
DE#define  DOT_XYZ(p1, p2)  ((p1).x * (p2).x + (p1).y * (p2).y + (p1).z * (p2).z)
DE
DE
DE#define NORMALIZE_XYZ(v)  { float __tmpnormval; \
DE			    __tmpnormval = (double)sqrt((v).x*(v).x +  \
DE					    (v).y*(v).y + (v).z*(v).z); \
DE			      (v).x /= __tmpnormval; \
DE			      (v).y /= __tmpnormval;  \
DE			      (v).z /= __tmpnormval;}
DE#define R_NORMALIZE_XYZ(v)  { float __tmpnormval; \
DE			    (((__tmpnormval = (double)sqrt((v).x*(v).x +  \
DE					    (v).y*(v).y + (v).z*(v).z)) \
DE			      == 0.0) ? FALSE : ((v).x /= __tmpnormval,  \
DE						 (v).y /= __tmpnormval,  \
DE						 (v).z /= __tmpnormval,TRUE));}
DE
DE#define RGB_ADD(c,a,e)   { (c).r = (a).r + (e).r; (c).g = (a).g + (e).g; \
DE                           (c).b = (a).b + (e).b; }
DE#define RGB_SUB(c,a,e)   { (c).r = (a).r - (e).r; (c).g = (a).g - (e).g; \
DE                           (c).b = (a).b - (e).b; }
DE#define RGB_MULT(c,a,e)  { (c).r = (a).r * (e).r; (c).g = (a).g * (e).g; \
DE                           (c).b = (a).b * (e).b; }
DE#define RGB_ADDC(c,a,e)  { (c).r = (a).r + (e); (c).g = (a).g + (e); \
DE                           (c).b = (a).b + (e); }
DE#define RGB_SUBC(c,a,e)  { (c).r = (a).r - (e); (c).g = (a).g - (e); \
DE                           (c).b = (a).b - (e); }
DE#define RGB_MULTC(c,a,e) { (c).r = (a).r * (e); (c).g = (a).g * (e); \
DE                           (c).b = (a).b * (e); }
DE#define RGB_DIVC(c,a,e)  { (c).r = (a).r / (e); (c).g = (a).g / (e); \
DE                           (c).b = (a).b / (e); }
DE#define RGB_COPY(c,a)    { (c).r = (a).r; (c).g = (a).g; (c).b = (a).b; }
DE#define RGB_COPYC(c,a)   { (c).r = (a); (c).g = (a); (c).b = (a); }
DE
DESHAR_EOF
DEfi
DEecho shar: "done with directory 'ebert'"
DEcd ..
DEexit 0
DE#	End of shell archive
SHAR_EOF
fi
echo shar: "extracting 'musgrave'" '(17201 characters)'
if test -f 'musgrave'
then
	echo shar: "will not over-write existing file 'musgrave'"
else
sed 's/^DE//' << \SHAR_EOF > 'musgrave'
DE/*
DE * The first section of this file of C code contains routines which did
DE * not appear in the course notes for the SIGGRAPH '92 Procedural Modelling
DE * and rendering course, but which were implemented subsequently to the 
DE * printing of those notes.  They represent improvements to the code which 
DE * appears in those notes; they should supercede that code.
DE * 
DE * F. Kenton Musgrave 8/92
DE */
DE
DE
DE/*
DE * "Earth" planet texture, as re-developed for the sequence of slides 
DE * illustrating the process of development of a complex procedural texture.
DE *
DE * F. Kenton Musgrave, 7/92
DE * All rights reserved.
DE */
DEvoid
DEEarth( texture, intersect, normal, colour, cmap, spectral_exp, lacunarity,
DE	octaves, bump_scale, multifractal, dist_scale, offset, sea_level,
DE	mtn_scale, lat_scale, nonlinear, purt_scale, map_exp, ice_caps,
DE	depth_scale, depth_max, mottle_limit, mottle_scale, mottle_dim,
DE	mottle_mag )
DE
DEVector	texture, intersect, *normal;
DEColour	*colour, *cmap;
DEdouble	spectral_exp, lacunarity, octaves, bump_scale;
DEint	multifractal;
DEdouble	dist_scale, offset, sea_level, mtn_scale, lat_scale, nonlinear,
DE	purt_scale, map_exp, ice_caps, depth_scale, depth_max,
DE	mottle_limit, mottle_scale, mottle_dim, mottle_mag;
DE
DE{
DE	Vector	bump, distort, point;
DE	double	chaos, latitude, purt;
DE	Colour	color;
DE	int	index;
DE
DE/*
DE * Choose between fractal bump functions: "ordinary" or multifractal fBm.
DE */
DE	if ( !multifractal )	/* use a "standard" fBm bump function */
DE		bump = VfBm(texture, spectral_exp, lacunarity, octaves );
DE	else {			/* use a "multifractal" fBm bump function */
DE			/* get "distortion" vector, as used with clouds */
DE		distort = DNoise( texture );
DE			/* scale distortion vector */
DE		SMULT( dist_scale, distort );
DE			/* insert distortion vector */
DE		texture = VADD( distort, texture );
DE			/* compute bump vector using displaced point */
DE		bump = MVfBm(texture, spectral_exp, lacunarity, octaves );
DE	}
DE
DE		/* get the "height" of the bump, displacing by offset */
DE	chaos = -DOT(bump, *normal) + offset;
DE		/* set bump for land masses (i.e., areas above "sea level") */
DE	if ( chaos > sea_level ) {
DE		chaos *= mtn_scale;
DE		sea_level *= mtn_scale;
DE		normal->x += bump_scale * bump.x;
DE		normal->y += bump_scale * bump.y;
DE		normal->z += bump_scale * bump.z;
DE		Unit(normal);
DE	}
DE
DE		/* if there's a colormap associated with the texture, use it */
DE	if( cmap ) {
DE			/* use a scaled "z" value for offset in color map */
DE		latitude = intersect.z;
DE			/* make climate symmetric about equator */
DE		if ( latitude < 0. ) latitude = -latitude;
DE			/* fractally purturb color map offset using "chaos" */
DE			/*  "nonlinear" scales purturbation-by-z */
DE			/*  "purt_scale" scales overall purturbation */
DE		latitude += chaos*(nonlinear*(1.-latitude) + purt_scale);
DE		if ( map_exp )
DE			latitude = lat_scale*pow(latitude,map_exp);
DE		else	latitude *= lat_scale;
DE			/* do oceans, adjusted for polar caps */
DE		if ( chaos > sea_level )
DE			index = 1 + (int)(latitude);
DE		else	index = 0;
DE		if ( ice_caps && latitude>ice_caps && index==0 )
DE			index = 255;
DE			/* clamp color map offset to bounds */
DE		if ( index > 255 ) index = 255;
DE		if ( index < 0 ) index = 0;
DE			/* set surface color to calculated color map entry */
DE		color = cmap[index];
DE
DE			/* darken the "deep waters" */
DE		if ( index == 0 )	{
DE			chaos -= sea_level;
DE			chaos *= depth_scale;
DE			chaos = MAX( chaos, -depth_max );
DE			color.red += chaos * color.red;
DE			color.green += chaos * color.green;
DE			color.blue += chaos * color.blue;
DE		}
DE			/* else mottle the color some */
DE		else if ( index < mottle_limit ) {
DE			point = texture;
DE			SMULT(mottle_scale, point);
DE			purt = fBm( point, mottle_dim, 2., 8. );
DE			color.red += color.red * 0.5*mottle_mag*purt;
DE			color.green += color.green * 0.175*mottle_mag*purt;
DE			color.blue += color.green * 0.5*mottle_mag*purt;
DE			/* too colorful:
DE			bump = VfBm( point, mottle_dim, 2., 8. );
DE			color.red += color.red * 0.5*mottle_mag*bump.x;
DE			color.green += color.green * 0.175*mottle_mag*bump.y;
DE			color.blue += color.green * 0.5*mottle_mag*bump.z;
DE			*/
DE			if ( color.red < 0. )	color.red = 0.;
DE			if ( color.green < 0. )	color.green = 0.;
DE			if ( color.blue < 0. )	color.blue = 0.;
DE			if ( color.red > 1. )	color.red = 1.;
DE			if ( color.green > 1. )	color.green = 1.;
DE			if ( color.blue > 1. )	color.blue = 1.;
DE		}
DE		*colour = color;
DE	}
DE} /* Earth() */
DE
DE
DE/* fBm with proper spectral construction - FKM 7/92 */
DEdouble 
DEfBm( point, spectral_exp, lacunarity, octaves )
DE	Vector	point;
DE	double	spectral_exp, lacunarity, octaves;
DE{
DE	register double	i, result, amplitude, frequency=1.0;
DE
DE	result = Noise3( point );
DE	for( i=octaves-1; i>0.0; i-- ) {
DE		point.x *= lacunarity;
DE		point.y *= lacunarity;
DE		point.z *= lacunarity;
DE		frequency *= lacunarity;
DE                amplitude = pow( frequency, -spectral_exp );
DE		if ( i < 1.0 ) amplitude *= i;	/* for octaves remainder */
DE		result += amplitude * Noise3( point );
DE		if (amplitude < VERY_SMALL) break;
DE	}
DE	return( result );
DE} /* fBm() */
DE
DE
DE/* "multifractal" fBm, with proper spectral construction - 8/92 */
DEdouble 
DEMfBm( point, spectral_exp, lacunarity, octaves, offset, threshold )
DE	Vector	point;
DE	double	spectral_exp, lacunarity, octaves, offset, threshold;
DE{
DE	register double	octave, result, amplitude, frequency=1.0;
DE	register double signal, weight;
DE
DE	signal = result = 0.5*(offset + Noise3( point ));
DE	for( octave=octaves-1; octave>0.0; octave-- ) {
DE		point.x *= lacunarity;
DE		point.y *= lacunarity;
DE		point.z *= lacunarity;
DE		frequency *= lacunarity;
DE                amplitude = pow( frequency, -spectral_exp );
DE		/* take care of "octaves" remainder */
DE		if ( octave < 1.0 ) amplitude *= octave;
DE
DE		/* weight successive contributions by previous signal */
DE		weight = signal / threshold;
DE		if ( weight > 1.0 ) weight = 1.0;
DE		if ( weight < 0.0 ) weight = 0.0;
DE		signal = weight * 0.5*(offset + Noise3( point ));
DE		result += amplitude * signal;
DE		if (weight<VERY_SMALL || amplitude<VERY_SMALL) break;
DE	}
DE	return( result );
DE} /* MfBm() */
DE
DE
DE/****************************************************************************/
DE
DE/* 
DE * This section contains the exact code which appeared in the course notes.
DE * Please be advised that, with the exception of the color map array,
DE * you'd be probably be better off using the preceding, updated & revamped
DE * code.  -FKM 8/92
DE */ 
DE
DE/*
DE * Below is an example of an elaborate procedural texture which, 
DE * when applied to the surface of a sphere using the right 20-something
DE * magic parameter values and a just-so color map, can do a nice imitation
DE * of an Earth-like planet, literally creating (part of) a "virtual world".
DE * Note that this function is surface-coloring and bump-map effect,
DE * thus the sphere remains geometrically smooth.
DE * 
DE * There are 22 arguments to this procedure, stored in the textArg[] array.
DE * The following vector of parameter values (plus a color map) created the
DE * planet seen in the Color Plate 1 of these course notes:
DE *
DE * .25 2 -.48 10 0 .45 300 2.6 .7 20 .75 220 170 20 .6 4 1.125 1.2 -.085 1 .3 0
DE *
DE * The code is in the form of one of the (many) cases in an enormous "switch"
DE * statement in our ray tracer (John Amanatide's and Andrew Woo's 
DE * "Optik" from the University of Toronto's Dynamic Graphics Project).
DE */
DE
DEcase PLANET:
DE/* 
DE * Perform intitialization of fractal (spectral exponent) parameter:
DE *	textArg[1] is lacunarity, or the gap between frequencies
DE *	textArg[2] is the fractal codimension parameter
DE *	textArg[4] serves as static storage for the computed exponent
DE *		   (note that this is an idiosyncrasy of our implementation)
DE *	firstPlanetCall is a static flag variable, initilaized to TRUE
DE */
DE	if(firstPlanetCall) {
DE		textArg[4] = pow(textArg[1],(-0.5-textArg[2]));
DE		firstPlanetCall = FALSE;
DE	}
DE
DE/*
DE * Choose between fractal bump functions, based on the value in textArg[19]
DE * 	(note that the high parameter number -- 19 -- indicates that this was
DE * 	 added quite late in the development of the texture).
DE */
DE	if ( !textArg[19] )	/* use a "standard" fBm bump function */
DE		bump = VfBm(texture, textArg[4], textArg[1], textArg[3]);
DE	else {			/* use a "multifractal" fBm bump function */
DE			/* get "distortion" vector, as used with clouds */
DE		distort = DNoise( texture );
DE			/* scale distortion vector */
DE		SMULT( textArg[20], distort );
DE			/* insert distortion vector */
DE		texture = VADD( distort, texture );
DE			/* compute bump vector using displaced point */
DE		bump = MVfBm(texture, textArg[4], textArg[1], textArg[3]);
DE	}
DE
DE		/* get the "height" of the bump, displacing by textArg[18] */
DE	chaos = -DOT(bump, hit->normal) + textArg[18];
DE		/* set bump for land masses (i.e., areas above "sea level") */
DE	if( chaos > 0. ) {
DE		chaos *= textArg[5];
DE		hit->normal.x += textArg[0] * bump.x;
DE		hit->normal.y += textArg[0] * bump.y;
DE		hit->normal.z += textArg[0] * bump.z;
DE		Normalize(&hit->normal)
DE	}
DE
DE		/* if there's a colormap associated with the texture, use it */
DE	if( cmap ) {
DE			/* use a scaled "z" value for offset in color map */
DE		temp = ABS(hit->intersect.z)*textArg[16];
DE			/* fractally purturb color map offset using "chaos" */
DE			/*  textArg[7] scales purturbation-by-z */
DE			/*  textArg[17] scales overall purturbation */
DE		temp = chaos*(textArg[7]*(1.-temp) + textArg[17]) + temp;
DE		if ( temp > 0. )	/* if above "sea level" */
DE				/* "mountains" appear too "chunky", */
DE				/* so exponentiate the color map offset */
DE			offset = (int)(textArg[6]*pow(temp,textArg[15]));
DE		else	offset = 0;	/* (don't mess with oceans) */
DE
DE			/* now do oceans; textArg[11] sets polar ice caps */
DE		if ((offset < 0.) || ((chaos <= 0.) && (offset < textArg[11])))
DE			offset = 0;
DE			
DE			/* clamp color map offset to upper bound */
DE		if ( offset > 255 ) offset = 255;
DE			/* set surface color to calculated color map entry */
DE		color = cmap[offset];
DE
DE			/* darken the "deep waters" */
DE			/* (note that "chaos" is less than 0 here) */
DE		if ( offset == 0 ) {
DE				/* scale the effect */
DE			chaos *= textArg[9];
DE				/* make the effect nonlinear, according to */
DE				/* the whim encoded in in textArg[21] */
DE			if ( textArg[21] )
DE				chaos *= 1.-texture.z*texture.z;
DE
DE				/* limit how dark deepest waters get */
DE			if ( chaos < -textArg[10] ) chaos = -textArg[10];
DE				/* now darken the color of the deeper waters */
DE			color.red += chaos * color.red;
DE			color.green += chaos * color.green;
DE			color.blue += chaos * color.blue;
DE		}
DE			/* else we are in the landmass areas, where the color */
DE			/* is ... boring, so we'll mottle it with a vector-valued */
DE			/* fBm, interpreted as an RGB value */
DE		else if ( offset < textArg[12] ) {	/* don't mottle snow! */
DE				/* scale size of color-bumps */
DE			SMULT(textArg[13], texture);
DE				/* get the vector-valued fBm */
DE				/* (note that we've hard-coded some constants */
DE				/*  in a feeble effort to fight parameter */
DE				/*  proliferation.) */
DE			bump = VfBm( texture, textArg[14], 2., 8. );
DE				/* using only bump.x is a "feature", */
DE				/* not a bug (don't ask me why!); */
DE				/* more hard-coded constants used, */
DE				/* to lessen parameter proliferation */
DE			color.red += color.red * 0.5*textArg[8]*bump.x;
DE			color.green += color.green * 0.175*textArg[8]*bump.x;
DE			color.blue += color.green * 0.5*textArg[8]*bump.x;
DE
DE				/* now clamp errant color values */
DE			if ( color.red < 0. ) color.red = 0.;
DE			if ( color.green < 0. ) color.green = 0.;
DE			if ( color.blue < 0. ) color.blue = 0.;
DE			if ( color.red > 1. ) color.red = 1.;
DE			if ( color.green > 1. ) color.green = 1.;
DE			if ( color.blue > 1. ) color.blue = 1.;
DE		}
DE	} else {		/* no color map, so just use mottled texture */
DE		color.red *= chaos;
DE		color.green *= chaos;
DE		color.blue *= chaos;
DE	}
DE	break;	/* whew! */
DE
DE
DE/*
DE * And now for some of the auxilliary functions and data referred to above...
DE */
DE
DE/* fBm constructed with VLNoise() */
DEdouble VLfBm( point, omega, lambda, octaves )
DE	Vector	point;
DE	double	omega, lambda, octaves;
DE{
DE	register double	l, a, o;
DE	register int	i;
DE	register Vector	tp;
DE	double VLNoise();
DE
DE	l = lambda;  o = omega;
DE	a = VLNoise( point, 1.0 );
DE	for( i=1; i<octaves; i++ ) {
DE		tp.x = l*point.x;
DE		tp.y = l*point.y;
DE		tp.z = l*point.z;
DE		a += o * VLNoise( tp, 1.0 );
DE		l *= lambda;
DE		o *= omega;
DE		if (o < VERY_SMALL) break;
DE	}
DE	return( a );
DE} /* VLfBm() */
DE
DE/* vector-valued fBm */
DEVector VfBm( point, omega, lambda, octaves )
DE	Vector	point;
DE	double	omega, lambda, octaves;
DE{
DE	register double	l, o;
DE	register int	i;
DE	register Vector	tp, n, a;
DE
DE	l = lambda;  o = omega;
DE	a = DNoise( point );
DE	for( i=1; i<octaves; i++ ) {
DE		tp.x = l*point.x;
DE		tp.y = l*point.y;
DE		tp.z = l*point.z;
DE		n = DNoise( tp );
DE		a.x += o * n.x;
DE		a.y += o * n.y;
DE		a.z += o * n.z;
DE		l *= lambda;
DE		o *= omega;
DE		if (o < VERY_SMALL) break;
DE	}
DE	return( a );
DE} /* VfBm() */
DE
DE/* vector-valued multifractal fBm routine */
DEVector MVfBm( point, omega, lambda, octaves )
DE	Vector	point;
DE	double	omega, lambda, octaves;
DE{
DE	register double	tmp, lacunarity, o, weight;
DE	register int	i;
DE	register Vector	tp, tv, result;
DE	double		Noise(), VLNoise();
DE	Vector		DNoise();
DE
DE	lacunarity = lambda;  o = omega;  tp = point;
DE	result.x = 0.0;  result.y = 0.0;  result.z = 0.0;
DE		/* get intitial value */
DE	weight = VLNoise( tp, 1.5 );
DE	if ( weight < 0. ) weight = -weight;
DE	tv = DNoise(tp);
DE	result = SMULT( weight, tv );
DE	for( i=1; i<octaves; i++ ) {
DE		tp = SMULT( lacunarity, tp );
DE			/* get subsequent values, weighted by previous value */
DE		weight *= o * ( N_OFFSET + Noise(tp) ) ;
DE		if ( weight < 0. ) weight = -weight;
DE		if ( weight > 1.0 ) weight = 1.0;
DE		if ( (weight<VERY_SMALL) && (weight>-VERY_SMALL) ) break;
DE		tv = DNoise(tp);
DE		tmp = MIN( weight, omega );
DE		tv = SMULT( tmp, tv );
DE		result = VADD( tv, result );
DE		o *= omega;
DE	} /* for */
DE
DE	return( result );
DE} /* MVfBm() */
DE
DE/*
DE * And now for the maniacs (I know you're out there!) who'd type this in,
DE * here's the color map used to create the planet seen in Color Plate 1:
DE */
DEchar planet_map[256][3] =
DE	{{1,14,81}, {176,134,80}, {170,123,72}, {164,113,64}, {158,103,56},
DE	{153,93,48}, {153,92,46}, {153,90,44}, {154,88,42}, {154,86,40},
DE	{154,84,38}, {154,83,36}, {155,81,34}, {155,79,32}, {155,77,30},
DE	{156,76,28}, {155,74,26}, {154,72,24}, {153,70,22}, {153,68,21},
DE	{152,66,19}, {151,64,17}, {150,62,15}, {150,60,14}, {149,58,12},
DE	{148,56,10}, {147,54,8}, {147,52,7}, {146,50,5}, {145,48,3},
DE	{144,46,1}, {144,45,0}, {142,45,1}, {141,46,2}, {139,47,3},
DE	{138,47,4}, {137,48,5}, {135,49,6}, {134,49,7}, {133,50,8},
DE	{131,51,9}, {130,51,10}, {128,52,11}, {127,53,12}, {126,53,13},
DE	{124,54,14}, {123,55,15}, {122,56,16}, {121,57,17}, {120,57,17},
DE	{119,58,18}, {118,59,19}, {117,59,20}, {116,60,20}, {114,61,21},
DE	{113,61,22}, {112,62,22}, {111,63,23}, {110,63,24}, {109,64,25},
DE	{108,65,25}, {107,65,26}, {106,66,27}, {105,67,28}, {103,68,28},
DE	{101,68,26}, {99,69,24}, {97,69,23}, {95,70,21}, {93,70,19},
DE	{92,71,18}, {90,72,16}, {88,72,14}, {86,73,13}, {84,73,11},
DE	{83,74,9}, {81,75,8}, {79,75,6}, {77,76,4}, {76,77,3},
DE	{74,76,3}, {73,75,4}, {71,74,5}, {70,74,5}, {68,73,6},
DE	{67,72,7}, {65,71,7}, {64,71,8}, {62,70,9}, {61,69,9},
DE	{59,68,10}, {58,68,11}, {56,67,11}, {55,66,12}, {53,65,13},
DE	{52,65,14}, {51,64,14}, {49,63,15}, {48,62,16}, {46,62,16},
DE	{45,61,17}, {43,60,18}, {42,59,18}, {40,59,19}, {39,58,20},
DE	{37,57,20}, {36,56,21}, {34,56,22}, {33,55,22}, {31,54,23},
DE	{30,53,24}, {29,53,25}, {28,52,25}, {27,51,25}, {25,50,24},
DE	{24,50,24}, {22,49,24}, {21,48,23}, {19,47,23}, {19,47,23},
DE	{19,47,23}, {19,47,23}, {19,47,24}, {19,47,24}, {19,47,24},
DE	{20,47,25}, {20,47,25}, {20,47,25}, {20,47,26}, {20,47,26},
DE	{20,47,26}, {21,47,27}, {21,47,27}, {21,47,27}, {21,47,28},
DE	{21,47,28}, {21,47,28}, {21,47,28}, {21,47,28}, {21,48,28},
DE	{21,48,28}, {21,48,28}, {21,48,28}, {21,49,28}, {21,49,28},
DE	{21,49,28}, {21,49,28}, {21,50,28}, {21,50,28}, {21,50,28},
DE	{21,50,28}, {21,51,28}, {23,52,30}, {26,53,32}, {28,54,34},
DE	{31,55,37}, {34,56,39}, {36,57,41}, {39,58,43}, {41,59,46},
DE	{44,60,48}, {47,61,50}, {49,62,53}, {52,63,55}, {54,64,57},
DE	{57,65,59}, {60,66,62}, {62,67,64}, {65,68,66}, {68,69,69},
DE	{69,69,69}, {69,69,69}, {69,70,70}, {70,70,70}, {70,70,70},
DE	{70,70,70}, {71,71,71}, {71,70,70}, {72,70,69}, {72,69,69},
DE	{73,69,68}, {74,68,68}, {74,68,67}, {75,67,67}, {76,67,66},
DE	{76,66,66}, {77,66,65}, {78,65,65}, {78,65,64}, {79,64,64},
DE	{80,64,63}, {81,64,63}, {82,63,62}, {84,62,61}, {94,74,73},
DE	{104,87,86}, {115,100,99}, {125,113,112}, {135,125,125}, {146,138,138},
DE	{156,151,151}, {166,164,164}, {177,177,177}, {196,196,196}, {216,216,216},
DE	{235,235,235}, {255,255,255}, {255,255,255}, {255,255,255}, {255,255,255},
DE	{255,255,255}, {255,255,255}, {255,255,255}, {255,255,255}, {255,255,255},
DE	{255,255,255}, {255,255,255}, {255,255,255}, {255,255,255}, {255,255,255},
DE	{255,255,255}, {255,255,255}, {255,255,255}, {255,255,255}, {255,255,255},
DE	{255,255,255}, {255,255,255}, {255,255,255}, {255,255,255}, {255,255,255},
DE	{255,255,255}, {255,255,255}, {255,255,255}, {255,255,255}, {255,255,255},
DE	{255,255,255}, {255,255,255}, {255,255,255}, {255,255,255}, {255,255,255},
DE	{255,255,255}, {255,255,255}, {255,255,255}, {255,255,255}, {255,255,255},
DE	{255,255,255}, {255,255,255}, {255,255,255}, {255,255,255}, {255,255,255},
DE	{255,255,255}, {255,255,255}, {255,255,255}, {255,255,255}, {255,255,255},
DE	{255,255,255}};
DE
SHAR_EOF
fi
echo shar: "extracting 'peachey_shar'" '(23711 characters)'
if test -f 'peachey_shar'
then
	echo shar: "will not over-write existing file 'peachey_shar'"
else
sed 's/^DE//' << \SHAR_EOF > 'peachey_shar'
DE#! /bin/sh
DE# This is a shell archive.  Remove anything before this line, then unpack
DE# it by saving it into a file and typing "sh file".  To overwrite existing
DE# files, type "sh file -c".  You can also feed this as standard input via
DE# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
DE# will see the following message at the end:
DE#		"End of shell archive."
DE# Contents:  README blue1.sl blue2.sl blue3.sl brickA.sl brickB.sl
DE#   perturb.sl scnoise.c star.sl starpaper.sl wall.sl
DE# Wrapped by peachey@daffy on Wed Aug 19 17:07:14 1992
DEPATH=/bin:/usr/bin:/usr/ucb ; export PATH
DEif test -f 'README' -a "${1}" != "-c" ; then 
DE  echo shar: Will not clobber existing file \"'README'\"
DEelse
DEecho shar: Extracting \"'README'\" \(615 characters\)
DEsed "s/^X//" >'README' <<'END_OF_FILE'
DEXThis directory contains C programs and RenderMan shading language programs
DEXwhich were used as examples in my presentation on Procedural Textures in
DEXCourse 23 at Siggraph 92 in Chicago.  Some of this might not make much sense
DEXwithout the course notes!
DEX
DEXAll material is copyright 1992 by Pixar.  Permission for non-profit use of
DEXthe material is granted provided that all copyright notices are retained and
DEXthe source of the material is appropriately acknowledged.
DEX
DEXPlease direct questions, comments, suggestions, etc. to peachey@pixar.com.
DEX
DEXDarwyn Peachey, Pixar, 1001 W. Cutting, Richmond, CA 94804  (510) 215-3495
DEX
DEEND_OF_FILE
DEif test 615 -ne `wc -c <'README'`; then
DE    echo shar: \"'README'\" unpacked with wrong size!
DEfi
DE# end of 'README'
DEfi
DEif test -f 'blue1.sl' -a "${1}" != "-c" ; then 
DE  echo shar: Will not clobber existing file \"'blue1.sl'\"
DEelse
DEecho shar: Extracting \"'blue1.sl'\" \(1212 characters\)
DEsed "s/^X//" >'blue1.sl' <<'END_OF_FILE'
DEX/*
DEX * Darwyn Peachey, Pixar, July 1992.
DEX * Copyright 1992 by Pixar.  Permission for non-profit use of this material
DEX * is granted provided that all copyright notices are retained and the
DEX * source of the material is appropriately acknowledged.
DEX *
DEX * Direct questions, comments, suggestions to peachey@pixar.com.
DEX */
DEX
DEX
DEX#define PALE_BLUE		color (0.25, 0.25, 0.35)
DEX#define MEDIUM_BLUE		color (0.10, 0.10, 0.30)
DEX#define DARK_BLUE		color (0.05, 0.05, 0.26)
DEX#define DARKER_BLUE		color (0.03, 0.03, 0.20)
DEX#define NNOISE	4
DEX
DEX/* Simplest of three implementations of a blue marble. */
DEX
DEXcolor
DEXmarble_color(float m)
DEX{
DEX	return m * PALE_BLUE;
DEX}
DEX
DEXsurface
DEXblue1(
DEX	uniform float Ka = 1;
DEX	uniform float Kd = 0.8;
DEX	uniform float Ks = 0.2;
DEX	uniform float texturescale = 1;
DEX	uniform float roughness = 0.1;
DEX     )
DEX{
DEX	color Ct;
DEX	point NN;
DEX	point PP;
DEX	float i, f, marble;
DEX
DEX	NN = normalize(faceforward(N, I));
DEX	PP = transform("shader", P) * 2.5 - (0.5, 0.5, 0.5);
DEX
DEX	marble = 0; f = 1;
DEX	for (i = 0; i < NNOISE; i += 1) {
DEX		marble += noise(PP * f) * 1/f;
DEX		f *= 2;
DEX	}
DEX
DEX	marble = clamp(4*marble - 3, 0, 1);
DEX	Ct = marble_color(marble);
DEX
DEX	Ci = Os * (Ct * (Ka * ambient() + Kd * diffuse(NN))
DEX		+ Ks * specular(NN, normalize(-I), roughness));
DEX}
DEEND_OF_FILE
DEif test 1212 -ne `wc -c <'blue1.sl'`; then
DE    echo shar: \"'blue1.sl'\" unpacked with wrong size!
DEfi
DE# end of 'blue1.sl'
DEfi
DEif test -f 'blue2.sl' -a "${1}" != "-c" ; then 
DE  echo shar: Will not clobber existing file \"'blue2.sl'\"
DEelse
DEecho shar: Extracting \"'blue2.sl'\" \(1226 characters\)
DEsed "s/^X//" >'blue2.sl' <<'END_OF_FILE'
DEX/*
DEX * Darwyn Peachey, Pixar, July 1992.
DEX * Copyright 1992 by Pixar.  Permission for non-profit use of this material
DEX * is granted provided that all copyright notices are retained and the
DEX * source of the material is appropriately acknowledged.
DEX *
DEX * Direct questions, comments, suggestions to peachey@pixar.com.
DEX */
DEX
DEX#define PALE_BLUE		color (0.25, 0.25, 0.35)
DEX#define MEDIUM_BLUE		color (0.10, 0.10, 0.30)
DEX#define DARK_BLUE		color (0.05, 0.05, 0.26)
DEX#define DARKER_BLUE		color (0.03, 0.03, 0.20)
DEX#define NNOISE	4
DEX
DEX/* Second of three implementations of a blue marble. */
DEX
DEXcolor
DEXmarble_color(float m)
DEX{
DEX	return mix(PALE_BLUE, DARKER_BLUE, m);
DEX}
DEX
DEXsurface
DEXblue2(
DEX	uniform float Ka = 1;
DEX	uniform float Kd = 0.8;
DEX	uniform float Ks = 0.2;
DEX	uniform float texturescale = 1;
DEX	uniform float roughness = 0.1;
DEX     )
DEX{
DEX	color Ct;
DEX	point NN;
DEX	point PP;
DEX	float i, f, marble;
DEX
DEX	NN = normalize(faceforward(N, I));
DEX	PP = transform("shader", P) * 2.5 - (0.5, 0.5, 0.5);
DEX
DEX	marble = 0; f = 1;
DEX	for (i = 0; i < NNOISE; i += 1) {
DEX		marble += noise(PP * f) * 1/f;
DEX		f *= 2;
DEX	}
DEX
DEX	marble = clamp(4*marble - 3, 0, 1);
DEX	Ct = marble_color(marble);
DEX
DEX	Ci = Os * (Ct * (Ka * ambient() + Kd * diffuse(NN))
DEX		+ Ks * specular(NN, normalize(-I), roughness));
DEX}
DEEND_OF_FILE
DEif test 1226 -ne `wc -c <'blue2.sl'`; then
DE    echo shar: \"'blue2.sl'\" unpacked with wrong size!
DEfi
DE# end of 'blue2.sl'
DEfi
DEif test -f 'blue3.sl' -a "${1}" != "-c" ; then 
DE  echo shar: Will not clobber existing file \"'blue3.sl'\"
DEelse
DEecho shar: Extracting \"'blue3.sl'\" \(1376 characters\)
DEsed "s/^X//" >'blue3.sl' <<'END_OF_FILE'
DEX/*
DEX * Darwyn Peachey, Pixar, July 1992.
DEX * Copyright 1992 by Pixar.  Permission for non-profit use of this material
DEX * is granted provided that all copyright notices are retained and the
DEX * source of the material is appropriately acknowledged.
DEX *
DEX * Direct questions, comments, suggestions to peachey@pixar.com.
DEX */
DEX
DEX#define PALE_BLUE		color (0.25, 0.25, 0.35)
DEX#define MEDIUM_BLUE		color (0.10, 0.10, 0.30)
DEX#define DARK_BLUE		color (0.05, 0.05, 0.26)
DEX#define DARKER_BLUE		color (0.03, 0.03, 0.20)
DEX#define NNOISE	4
DEX
DEX/* Best of three implementations of a blue marble. */
DEX
DEXcolor
DEXmarble_color(float m)
DEX{
DEX	return color spline(m,
DEX		PALE_BLUE, PALE_BLUE,
DEX		MEDIUM_BLUE, MEDIUM_BLUE, MEDIUM_BLUE,
DEX		PALE_BLUE, PALE_BLUE,
DEX		DARK_BLUE, DARK_BLUE,
DEX		DARKER_BLUE, DARKER_BLUE,
DEX		PALE_BLUE, DARKER_BLUE);
DEX}
DEX
DEXsurface
DEXblue3(
DEX	uniform float Ka = 1;
DEX	uniform float Kd = 0.8;
DEX	uniform float Ks = 0.2;
DEX	uniform float texturescale = 1;
DEX	uniform float roughness = 0.1;
DEX     )
DEX{
DEX	color Ct;
DEX	point NN;
DEX	point PP;
DEX	float i, f, marble;
DEX
DEX	NN = normalize(faceforward(N, I));
DEX	PP = transform("shader", P) * 2.5 - (0.5, 0.5, 0.5);
DEX
DEX	marble = 0; f = 1;
DEX	for (i = 0; i < NNOISE; i += 1) {
DEX		marble += noise(PP * f) * 1/f;
DEX		f *= 2;
DEX	}
DEX
DEX	marble = clamp(4*marble - 3, 0, 1);
DEX	Ct = marble_color(marble);
DEX
DEX	Ci = Os * (Ct * (Ka * ambient() + Kd * diffuse(NN))
DEX		+ Ks * specular(NN, normalize(-I), roughness));
DEX}
DEEND_OF_FILE
DEif test 1376 -ne `wc -c <'blue3.sl'`; then
DE    echo shar: \"'blue3.sl'\" unpacked with wrong size!
DEfi
DE# end of 'blue3.sl'
DEfi
DEif test -f 'brickA.sl' -a "${1}" != "-c" ; then 
DE  echo shar: Will not clobber existing file \"'brickA.sl'\"
DEelse
DEecho shar: Extracting \"'brickA.sl'\" \(1357 characters\)
DEsed "s/^X//" >'brickA.sl' <<'END_OF_FILE'
DEX/*
DEX * Darwyn Peachey, Pixar, July 1992.
DEX * Copyright 1992 by Pixar.  Permission for non-profit use of this material
DEX * is granted provided that all copyright notices are retained and the
DEX * source of the material is appropriately acknowledged.
DEX *
DEX * Direct questions, comments, suggestions to peachey@pixar.com.
DEX */
DEX
DEX#define BRICKWIDTH	0.25
DEX#define BRICKHEIGHT	0.08
DEX#define MORTARTHICKNESS	0.01
DEX
DEX#define BMWIDTH		(BRICKWIDTH+MORTARTHICKNESS)
DEX#define BMHEIGHT	(BRICKHEIGHT+MORTARTHICKNESS)
DEX#define MWF		(MORTARTHICKNESS*0.5/BMWIDTH)
DEX#define MHF		(MORTARTHICKNESS*0.5/BMHEIGHT)
DEX
DEX/* Regular brick: no noise added, no antialiasing. */
DEX
DEXsurface
DEXbrickA(
DEX	uniform float Ka = 1;
DEX	uniform float Kd = 1;
DEX	uniform color Cbrick = color (0.5, 0.15, 0.14);
DEX	uniform color Cmortar = color (0.5, 0.5, 0.5);
DEX     )
DEX{
DEX	color Ct;
DEX	point NN;
DEX	float ss, tt, sbrick, tbrick, w, h;
DEX	float scoord = s;
DEX	float tcoord = t;
DEX
DEX	NN = normalize(faceforward(N, I));
DEX
DEX	ss = scoord / BMWIDTH;
DEX	tt = tcoord / BMHEIGHT;
DEX
DEX	if (mod(tt*0.5,1) > 0.5)
DEX		ss += 0.5;	/* shift alternate rows */
DEX
DEX	tbrick = floor(tt);	/* which brick? */
DEX	sbrick = floor(ss);	/* which brick? */
DEX	ss -= sbrick;
DEX	tt -= tbrick;
DEX
DEX	w = step(MWF,ss) - step(1-MWF,ss);
DEX	h = step(MHF,tt) - step(1-MHF,tt);
DEX
DEX	Ct = mix(Cmortar, Cbrick, w*h);
DEX
DEX	/* "matte" reflection model */
DEX	Ci = Os * Ct * (Ka * ambient() + Kd * diffuse(NN));
DEX}
DEEND_OF_FILE
DEif test 1357 -ne `wc -c <'brickA.sl'`; then
DE    echo shar: \"'brickA.sl'\" unpacked with wrong size!
DEfi
DE# end of 'brickA.sl'
DEfi
DEif test -f 'brickB.sl' -a "${1}" != "-c" ; then 
DE  echo shar: Will not clobber existing file \"'brickB.sl'\"
DEelse
DEecho shar: Extracting \"'brickB.sl'\" \(1900 characters\)
DEsed "s/^X//" >'brickB.sl' <<'END_OF_FILE'
DEX/*
DEX * Darwyn Peachey, Pixar, July 1992.
DEX * Copyright 1992 by Pixar.  Permission for non-profit use of this material
DEX * is granted provided that all copyright notices are retained and the
DEX * source of the material is appropriately acknowledged.
DEX *
DEX * Direct questions, comments, suggestions to peachey@pixar.com.
DEX */
DEX
DEX#define boxstep(a,b,x) clamp(((x)-(a))/((b)-(a)),0,1)
DEX
DEX#define BRICKWIDTH	0.25
DEX#define BRICKHEIGHT	0.08
DEX#define MORTARTHICKNESS	0.01
DEX
DEX#define BMWIDTH		(BRICKWIDTH+MORTARTHICKNESS)
DEX#define BMHEIGHT	(BRICKHEIGHT+MORTARTHICKNESS)
DEX#define MWF		(MORTARTHICKNESS*0.5/BMWIDTH)
DEX#define MHF		(MORTARTHICKNESS*0.5/BMHEIGHT)
DEX
DEX/* Antialiased version of brickA.  Noise has also been added to row offset. */
DEX
DEXsurface
DEXbrickB(
DEX	uniform float Ka = 1;
DEX	uniform float Kd = 1;
DEX	uniform color Cbrick = color (0.5, 0.15, 0.14);
DEX	uniform color Cmortar = color (0.5, 0.5, 0.5);
DEX     )
DEX{
DEX	color Ct;
DEX	point NN;
DEX	float ss, tt, sbrick, tbrick, w, h;
DEX	float scoord = s;
DEX	float tcoord = t;
DEX	float swidth, twidth;
DEX
DEX	NN = normalize(faceforward(N, I));
DEX
DEX	ss = scoord / BMWIDTH;
DEX	tt = tcoord / BMHEIGHT;
DEX
DEX	if (mod(tt*0.5,1) > 0.5)
DEX		ss += 0.5;	/* shift alternate rows */
DEX
DEX	tbrick = floor(tt);	/* which brick? */
DEX	/* Add noise to perturb row position slightly. */
DEX	ss += 0.2 * (noise(tbrick+0.5) - 0.5);
DEX	sbrick = floor(ss);	/* which brick? */
DEX	swidth = abs(Du(ss)*du) + abs(Dv(ss)*dv);
DEX	twidth = abs(Du(tt)*du) + abs(Dv(tt)*dv);
DEX	ss -= sbrick;
DEX	tt -= tbrick;
DEX
DEX	if (swidth >= 1)
DEX		w = 1 - 2 * MWF;
DEX	else
DEX		w = clamp(boxstep(MWF-swidth,MWF,ss),max(1-MWF/swidth,0),1)
DEX		    - clamp(boxstep(1-MWF-swidth,1-MWF,ss),0,2*MWF/swidth);
DEX
DEX	if (twidth >= 1)
DEX		h = 1 - 2 * MHF;
DEX	else
DEX		h = clamp(boxstep(MHF-twidth,MHF,tt),max(1-MHF/twidth,0),1)
DEX		    - clamp(boxstep(1-MHF-twidth,1-MHF,tt),0,2*MHF/twidth);
DEX
DEX	Ct = mix(Cmortar, Cbrick, w*h);
DEX
DEX	/* "matte" reflection model */
DEX	Ci = Os * Ct * (Ka * ambient() + Kd * diffuse(NN));
DEX}
DEEND_OF_FILE
DEif test 1900 -ne `wc -c <'brickB.sl'`; then
DE    echo shar: \"'brickB.sl'\" unpacked with wrong size!
DEfi
DE# end of 'brickB.sl'
DEfi
DEif test -f 'perturb.sl' -a "${1}" != "-c" ; then 
DE  echo shar: Will not clobber existing file \"'perturb.sl'\"
DEelse
DEecho shar: Extracting \"'perturb.sl'\" \(796 characters\)
DEsed "s/^X//" >'perturb.sl' <<'END_OF_FILE'
DEX/*
DEX * Darwyn Peachey, Pixar, July 1992.
DEX * Copyright 1992 by Pixar.  Permission for non-profit use of this material
DEX * is granted provided that all copyright notices are retained and the
DEX * source of the material is appropriately acknowledged.
DEX *
DEX * Direct questions, comments, suggestions to peachey@pixar.com.
DEX */
DEX
DEX/* Perturbed texture access: noise is used to perturb texture coordinates. */
DEX
DEXsurface
DEXperturb(
DEX	uniform float Ka = 1;
DEX	uniform float Kd = 1;
DEX     )
DEX{
DEX	color Ct;
DEX	point NN;
DEX	point Psh;
DEX	float ss, tt;
DEX
DEX	NN = normalize(faceforward(N, I));
DEX
DEX	Psh = transform("shader", P);
DEX
DEX	ss = s + 0.2 * noise(Psh) - 0.1;
DEX	tt = t + 0.2 * noise(Psh+(1.5,6.7,3.4)) - 0.1;
DEX
DEX	Ct = texture("mandrill.tx", ss, tt);
DEX
DEX	/* "matte" reflection model */
DEX	Ci = Os * Ct * (Ka * ambient() + Kd * diffuse(NN));
DEX}
DEEND_OF_FILE
DEif test 796 -ne `wc -c <'perturb.sl'`; then
DE    echo shar: \"'perturb.sl'\" unpacked with wrong size!
DEfi
DE# end of 'perturb.sl'
DEfi
DEif test -f 'scnoise.c' -a "${1}" != "-c" ; then 
DE  echo shar: Will not clobber existing file \"'scnoise.c'\"
DEelse
DEecho shar: Extracting \"'scnoise.c'\" \(6003 characters\)
DEsed "s/^X//" >'scnoise.c' <<'END_OF_FILE'
DEX
DEX/*
DEX * scnoise.c 1.3 (Pixar Animation) 8/19/92
DEX * An implementation of John Lewis' sparse convolution noise algorithm from
DEX * Proc. Siggraph 89, 263-270.  Written by Darwyn Peachey, 28 April 1992.
DEX * Copyright 1992 Pixar.  Permission is granted for nonprofit use of this
DEX * code provided that this notice is retained.
DEX *
DEX * Please direct comments, questions, improvements to peachey@pixar.com.
DEX *
DEX * This implementation is substantially slower than my implementation of
DEX * Ken Perlin's gradient interpolation noise.  On the other hand, I haven't
DEX * spent much time trying to optimize this code.
DEX *
DEX * Spatial coordinates index a voxel which has a number NIMP_PER_VOXEL of
DEX * pseudo-randomly distributed white noise impulses in it.  A filter is
DEX * discretely convolved with all noise impulses inside the filter radius to
DEX * determine the noise value.
DEX *
DEX * Simplification: assume that the filter is isotropic (radially symmetric).
DEX */
DEX
DEX#include <stdio.h>
DEX#include <stdlib.h>
DEX#include <math.h>
DEX
DEX#ifndef PI
DEX#define PI      3.141592653589793238
DEX#endif
DEX
DEX#define FLOOR(x)	(((x)<0 && (x)!=(int)(x))?(int)x-1:(int)x)
DEX#define LERP(t,x0,x1)	((x0) + (t)*((x1)-(x0)))
DEX#define IMOD(x,n,r)	{(r)=(x)%(n); if ((r) < 0) (r) += (n);}
DEX
DEX/*
DEX * Permutation table: a one-to-one mapping of 0:255 -> 0:255.
DEX * TABSIZE must be a power of 2 for the TABMASK to make sense.  The
DEX * predefined permutation in "permute" must be generated anew if TABSIZE
DEX * is changed.  Moreover, since "permute" is an unsigned char array
DEX * TABSIZE can't be larger than 256.  Note also that smaller TABSIZE
DEX * reduces the period of noise along the axes.
DEX * (I'm trying to convince you not to change it lightly.)
DEX */
DEX
DEX#define TABSIZE		(1<<8)		/* must be a power of 2 */
DEX#define TABMASK		(TABSIZE-1)	
DEX#define	PERMUTE(x)	permute[(x) & TABMASK]
DEX
DEXunsigned char permute[TABSIZE] = {
DEX	225,155,210,108,175,199,221,144,203,116, 70,213, 69,158, 33,252,
DEX	  5, 82,173,133,222,139,174, 27,  9, 71, 90,246, 75,130, 91,191,
DEX	169,138,  2,151,194,235, 81,  7, 25,113,228,159,205,253,134,142,
DEX	248, 65,224,217, 22,121,229, 63, 89,103, 96,104,156, 17,201,129,
DEX	 36,  8,165,110,237,117,231, 56,132,211,152, 20,181,111,239,218,
DEX	170,163, 51,172,157, 47, 80,212,176,250, 87, 49, 99,242,136,189,
DEX	162,115, 44, 43,124, 94,150, 16,141,247, 32, 10,198,223,255, 72,
DEX	 53,131, 84, 57,220,197, 58, 50,208, 11,241, 28,  3,192, 62,202,
DEX	 18,215,153, 24, 76, 41, 15,179, 39, 46, 55,  6,128,167, 23,188,
DEX	106, 34,187,140,164, 73,112,182,244,195,227, 13, 35, 77,196,185,
DEX	 26,200,226,119, 31,123,168,125,249, 68,183,230,177,135,160,180,
DEX	 12,  1,243,148,102,166, 38,238,251, 37,240,126, 64, 74,161, 40,
DEX	184,149,171,178,101, 66, 29, 59,146, 61,254,107, 42, 86,154,  4,
DEX	236,232,120, 21,233,209, 45, 98,193,114, 78, 19,206, 14,118,127,
DEX	 48, 79,147, 85, 30,207,219, 54, 88,234,190,122, 95, 67,143,109,
DEX	137,214,145, 93, 92,100,245,  0,216,186, 60, 83,105, 97,204, 52
DEX};
DEX
DEXtypedef struct {
DEX	float fx, fy, fz;
DEX	float value;
DEX} Impulse;
DEX
DEX#define NIMPULSE	1001
DEX#define NEXT_IMPULSE(h)	((h)>=NIMPULSE?0:(h)+1)
DEX#define VOXEL_START(ix,iy,iz)	PERMUTE(ix + PERMUTE(iy + PERMUTE(iz)))
DEX
DEXstatic void make_random_impulse_table(Impulse *table, int nimpulses);
DEX
DEXstatic float
DEXscnoise(float x, float y, float z, int nimpulses, float (*filter)(float))
DEX{
DEX	static Impulse table[NIMPULSE];
DEX	static int initialized;
DEX	Impulse *imp;
DEX	int i, j, k, h, n;
DEX	int ix, iy, iz;
DEX	float sum = 0;
DEX	float fx, fy, fz, dx, dy, dz, distsq;
DEX
DEX	/* Initialize the random impulse table if necessary. */
DEX	if (!initialized) {
DEX		make_random_impulse_table(table, NIMPULSE);
DEX		initialized = 1;
DEX	}
DEX
DEX	/* Find the indices (ix,iy,iz) of the unit voxel containing (x,y,z). */
DEX	ix = FLOOR(x); fx = x - ix;
DEX	iy = FLOOR(y); fy = y - iy;
DEX	iz = FLOOR(z); fz = z - iz;
DEX	
DEX	/* Perform the sparse convolution. */
DEX	for (i = -1; i <= 1; i++) {
DEX	  for (j = -1; j <= 1; j++) {
DEX	    for (k = -1; k <= 1; k++) {
DEX		/* Compute voxel hash code. */
DEX		h = VOXEL_START(ix+i,iy+j,iz+k);
DEX
DEX		for (n = nimpulses; n > 0; n--, h = NEXT_IMPULSE(h)) {
DEX			/* Convolve filter and impulse. */
DEX			imp = &table[h];
DEX			dx = fx - (i + imp->fx);
DEX			dy = fy - (j + imp->fy);
DEX			dz = fz - (k + imp->fz);
DEX			distsq = dx*dx + dy*dy + dz*dz;
DEX			sum += (*filter)(distsq) * imp->value;
DEX		}
DEX	    }
DEX	  }
DEX	}
DEX
DEX	return sum / nimpulses;
DEX}
DEX
DEX#define RANDMASK	0x7FFFFFFF
DEX#define FRANDOM		((random() & RANDMASK) / (double) RANDMASK)
DEX
DEXstatic void
DEXmake_random_impulse_table(Impulse *table, int nimpulses)
DEX{
DEX	int i;
DEX
DEX	srandom(555); /* Set random number generator seed. */
DEX	for (i = 0; i < nimpulses; i++, table++) {
DEX		/*
DEX		 * Generate a new pseudo-random impulse.
DEX		 * Note that the value of each impulse (imp->value) is
DEX		 * in the [0,1] interval and so the noise is also in
DEX		 * that interval (RenderMan noise definition).  If you
DEX		 * want a range of [-1,1], define SIGNED_NOISE.
DEX		 */
DEX		table->fx = FRANDOM;
DEX		table->fy = FRANDOM;
DEX		table->fz = FRANDOM;
DEX#ifdef SIGNED_NOISE
DEX		table->value = 2 * FRANDOM - 1;
DEX#else
DEX		table->value = FRANDOM;
DEX#endif
DEX	}
DEX}
DEX
DEX
DEX/* A very simple filter ... */
DEXstatic float
DEXcone(float x)
DEX{
DEX	if (x >= 1)
DEX		return 0;
DEX	return 1-x;
DEX}
DEX
DEX/* Cosine filter: return 0.5 * (1 + cos(PI * sqrt(x))) */
DEXstatic float
DEXcosfilter(float x)
DEX{
DEX#define NENTRIES 100
DEX	static int initialized = 0;
DEX	static float ftab[NENTRIES+1];
DEX	float f;
DEX	int i;
DEX
DEX	if (!initialized) {
DEX		initialized = 1;
DEX		for (i = 0; i < NENTRIES; i++) {
DEX			ftab[i] = 0.5 * (1+cos(PI * sqrt((double)i/NENTRIES)));
DEX		}
DEX		ftab[NENTRIES] = 0;
DEX	}
DEX	if (x < 0 || x >= 1)
DEX		return 0;
DEX	f = x * NENTRIES;
DEX	i = (int)f;
DEX	f -= i;
DEX	return LERP(f, ftab[i], ftab[i+1]);
DEX#undef  NENTRIES
DEX}
DEX
DEX
DEX/*
DEX * Simple driver to exercise the scnoise routine.
DEX * Evaluate noise at random points and print the result.
DEX */
DEX
DEXmain(int argc, char **argv)
DEX{
DEX	float x, y, z;
DEX	float f;
DEX	int i;
DEX
DEX#define RANDOM_COORD	(20 * ((random() & RANDMASK) / (double) RANDMASK) - 10)
DEX	for (i = 0; i < 20; i++) {
DEX		x = RANDOM_COORD;
DEX		y = RANDOM_COORD;
DEX		z = RANDOM_COORD;
DEX		f = scnoise(x, y, z, 4, cosfilter);
DEX		printf("noise(%g, %g, %g) = %g.\n", x, y, z, f);
DEX	}
DEX
DEX	exit(0);
DEX}
DEEND_OF_FILE
DEif test 6003 -ne `wc -c <'scnoise.c'`; then
DE    echo shar: \"'scnoise.c'\" unpacked with wrong size!
DEfi
DE# end of 'scnoise.c'
DEfi
DEif test -f 'star.sl' -a "${1}" != "-c" ; then 
DE  echo shar: Will not clobber existing file \"'star.sl'\"
DEelse
DEecho shar: Extracting \"'star.sl'\" \(1178 characters\)
DEsed "s/^X//" >'star.sl' <<'END_OF_FILE'
DEX/*
DEX * Darwyn Peachey, Pixar, July 1992.
DEX * Copyright 1992 by Pixar.  Permission for non-profit use of this material
DEX * is granted provided that all copyright notices are retained and the
DEX * source of the material is appropriately acknowledged.
DEX *
DEX * Direct questions, comments, suggestions to peachey@pixar.com.
DEX */
DEX
DEXsurface
DEXstar(
DEX	uniform float Ka = 1;
DEX	uniform float Kd = 1;
DEX	uniform color starcolor = color (1.0000,0.5161,0.0000);
DEX	uniform float npoints = 5;
DEX	uniform float sctr = 0.5;
DEX	uniform float tctr = 0.5;
DEX     )
DEX{
DEX	color Ct;
DEX	point NN;
DEX	float ss, tt, angle, r;
DEX	point d0, d1;
DEX	uniform float rmin = 0.07, rmax = 0.2;
DEX	uniform float starangle = 2*PI/npoints;
DEX	uniform point p0 = rmax*(cos(0),sin(0),0);
DEX	uniform point p1 = rmin*(cos(starangle/2),sin(starangle/2),0);
DEX
DEX	NN = normalize(faceforward(N, I));
DEX
DEX	ss = s - sctr; tt = t - tctr;
DEX	angle = atan(ss, tt) + PI;
DEX	r = sqrt(ss*ss + tt*tt);
DEX	angle /= starangle;
DEX	angle -= floor(angle);
DEX
DEX	if (angle >= 0.5)
DEX		angle = 1 - angle;
DEX	d0 = p1 - p0;
DEX	d1 = r*(cos(angle), sin(angle),0) - p0;
DEX	Ct = mix(Cs, starcolor, step(0, zcomp(d0^d1)));
DEX
DEX	/* "matte" reflection model */
DEX	Ci = Os * Ct * (Ka * ambient() + Kd * diffuse(NN));
DEX}
DEEND_OF_FILE
DEif test 1178 -ne `wc -c <'star.sl'`; then
DE    echo shar: \"'star.sl'\" unpacked with wrong size!
DEfi
DE# end of 'star.sl'
DEfi
DEif test -f 'starpaper.sl' -a "${1}" != "-c" ; then 
DE  echo shar: Will not clobber existing file \"'starpaper.sl'\"
DEelse
DEecho shar: Extracting \"'starpaper.sl'\" \(1494 characters\)
DEsed "s/^X//" >'starpaper.sl' <<'END_OF_FILE'
DEX/*
DEX * Darwyn Peachey, Pixar, July 1992.
DEX * Copyright 1992 by Pixar.  Permission for non-profit use of this material
DEX * is granted provided that all copyright notices are retained and the
DEX * source of the material is appropriately acknowledged.
DEX *
DEX * Direct questions, comments, suggestions to peachey@pixar.com.
DEX */
DEX
DEX#define NCELLS 10
DEX#define CELLSIZE (1/NCELLS)
DEX#define snoise(s,t)	(2*noise((s),(t))-1)
DEX
DEXsurface
DEXstarpaper(
DEX	uniform float Ka = 1;
DEX	uniform float Kd = 1;
DEX	uniform color starcolor = color (1.0000,0.5161,0.0000);
DEX	uniform float npoints = 5;
DEX     )
DEX{
DEX	color Ct;
DEX	point NN;
DEX	float ss, tt, angle, r;
DEX	float sctr, tctr, scell, tcell;
DEX	point d0, d1;
DEX	uniform float rmin = 0.007, rmax = 0.02;
DEX	uniform float starangle = 2*PI/npoints;
DEX	uniform point p0 = rmax*(cos(0),sin(0),0);
DEX	uniform point p1 = rmin*(cos(starangle/2),sin(starangle/2),0);
DEX
DEX	NN = normalize(faceforward(N, I));
DEX
DEX	ss = s * NCELLS; tt = t * NCELLS;
DEX	scell = floor(ss); tcell = floor(tt);
DEX	sctr = CELLSIZE * (scell + 0.5 + 0.6 * snoise(scell+0.5, tcell+0.5));
DEX	tctr = CELLSIZE * (tcell + 0.5 + 0.6 * snoise(scell+3.5, tcell+8.5));
DEX	ss = ss * CELLSIZE - sctr;
DEX	tt = tt * CELLSIZE - tctr;
DEX
DEX	angle = atan(ss, tt) + PI;
DEX	r = sqrt(ss*ss + tt*tt);
DEX	angle /= starangle;
DEX	angle -= floor(angle);
DEX
DEX	if (angle >= 0.5)
DEX		angle = 1 - angle;
DEX	d0 = p1 - p0;
DEX	d1 = r*(cos(angle), sin(angle),0) - p0;
DEX	Ct = mix(Cs, starcolor, step(0, zcomp(d0^d1)));
DEX
DEX	/* "matte" reflection model */
DEX	Ci = Os * Ct * (Ka * ambient() + Kd * diffuse(NN));
DEX}
DEEND_OF_FILE
DEif test 1494 -ne `wc -c <'starpaper.sl'`; then
DE    echo shar: \"'starpaper.sl'\" unpacked with wrong size!
DEfi
DE# end of 'starpaper.sl'
DEfi
DEif test -f 'wall.sl' -a "${1}" != "-c" ; then 
DE  echo shar: Will not clobber existing file \"'wall.sl'\"
DEelse
DEecho shar: Extracting \"'wall.sl'\" \(1376 characters\)
DEsed "s/^X//" >'wall.sl' <<'END_OF_FILE'
DEX/*
DEX * Darwyn Peachey, Pixar, July 1992.
DEX * Copyright 1992 by Pixar.  Permission for non-profit use of this material
DEX * is granted provided that all copyright notices are retained and the
DEX * source of the material is appropriately acknowledged.
DEX *
DEX * Direct questions, comments, suggestions to peachey@pixar.com.
DEX */
DEX
DEX/* A rough stucco wall implemented using bump-mapping. */
DEX
DEXsurface
DEXwall(float Ka = 0.4, Kd = 0.6;
DEX	float scale = 0.4, depth = 0.2)
DEX{
DEX	point PS, PQ, NN, PP;
DEX	float f, bh;
DEX
DEX	PS = transform("shader", P) * scale;
DEX	setxcomp(PQ, xcomp(PS) + ycomp(PS));
DEX	setycomp(PQ, xcomp(PS) - ycomp(PS));
DEX	setzcomp(PQ, zcomp(PS));
DEX	f = noise(PQ) + noise(PS);
DEX
DEX	/* add some higher frequencies to roughen edges */
DEX	f += 0.3 * (noise(PS * 3.713)-0.5);
DEX	f += 0.2 * (noise(PS * 7.436)-0.5);
DEX
DEX	/* threshold the function to get "mesas" in 0:1 range */
DEX	bh = smoothstep(1.1, 1.2, f);
DEX
DEX	/* add higher frequency bumps, only where the low freq ones aren't */
DEX	bh += (1 - bh) * smoothstep(0.7, 0.8, noise(PQ * 2.5));
DEX	bh += (1 - bh) * smoothstep(0.7, 0.8, noise(PS * 2.9));
DEX	bh += (1 - bh) * smoothstep(0.7, 0.8, noise(PS * 4.3));
DEX	bh += (1 - bh) * smoothstep(0.7, 0.8, noise(PQ * 5.7));
DEX
DEX	/* bh is between 0 and 1 */
DEX	NN = normalize(faceforward(N, I));
DEX	PP = P - NN * depth * (1 - bh);
DEX	NN = normalize(calculatenormal(PP));
DEX
DEX	Ci = Os * Cs * (Ka * ambient() + Kd * diffuse(NN));
DEX	Oi = Os;
DEX}
DEEND_OF_FILE
DEif test 1376 -ne `wc -c <'wall.sl'`; then
DE    echo shar: \"'wall.sl'\" unpacked with wrong size!
DEfi
DE# end of 'wall.sl'
DEfi
DEecho shar: End of shell archive.
DEexit 0
DE
SHAR_EOF
fi
echo shar: "extracting 'perlin'" '(6475 characters)'
if test -f 'perlin'
then
	echo shar: "will not over-write existing file 'perlin'"
else
sed 's/^DE//' << \SHAR_EOF > 'perlin'
DE                  EXCERPTED FROM SIGGRAPH 92, COURSE 23
DE                           PROCEDURAL MODELING
DE
DE                               Ken Perlin
DE                           New York University
DE
DE        3.6 TURBULENCE AND NOISE
DE
DE        3.6.1 The turbulence function
DE
DE        The turbulence function, which you use  to  make  marble,
DE        clouds,  explosions,  etc., is just a simple fractal gen-
DE        erating loop built on top of the noise function.   It  is
DE        not a real turbulence model at all.  The key trick is the
DE        use of the fabs() function, which makes the function have
DE        gradient discontinuity "fault lines" at all scales.  This
DE        fools the eye into thinking it is seeing the  results  of
DE        turbulent flow.  The turbulence() function gives the best
DE        results when used as a phase shift, as  in  the  familiar
DE        marble trick:
DE
DE
DE
DE            sin(point + turbulence(point) * point.x);
DE
DE
DE
DE        Note the second argument below, lofreq,  which  sets  the
DE        lowest  desired  frequency  component  of the turbulence.
DE        The third, hifreq, argument is used by  the  function  to
DE        ensure  that  the  turbulence  effect reaches down to the
DE        single pixel level, but no further.  I usually  set  this
DE        argument equal to the image resolution.
DE
DE
DE
DEfloat turbulence(point, lofreq, hifreq)
DEfloat point[3], freq, resolution;
DE{
DE        float noise3(), freq, t, p[3];
DE
DE        p[0] = point[0] + 123.456;
DE        p[1] = point[1];
DE        p[2] = point[2];
DE
DE        t = 0;
DE        for (freq = lofreq ; freq < hifreq ; freq *= 2.) {
DE                t += fabs(noise3(p)) / freq;
DE                p[0] *= 2.;
DE                p[1] *= 2.;
DE                p[2] *= 2.;
DE        }
DE        return t - 0.3; /* readjust to make mean value = 0.0 */
DE}
DE
DE
DE
DE        3.6.2 The noise function
DE
DE        noise3 is a rough approximation to "pink"  (band-limited)
DE        noise,  implemented  by  a  pseudorandom tricubic spline.
DE        Given a vector in  3-space,  it returns a  value  between
DE        -1.0  and 1.0.  There are two principal tricks to make it
DE        run fast:
DE
DE        - Precompute an array of pseudo-random unit  length  gra-
DE        dients g[n].
DE
DE        - Precompute a permutation  array  p[]  of  the  first  n
DE        integers.
DE
DE        Given the above two arrays,  any  integer  lattice  point
DE        (i,j,k)  can be quickly mapped to a pseudorandom gradient
DE        vector by:
DE
DE
DE            g[ (p[ (p[i] + j) % n ] + k) % n]
DE
DE
DE        By extending the g[] and p[] arrays, so that  g[n+i]=g[i]
DE        and  p[n+i]=p[i], the above lookup can be replaced by the
DE        (somewhat faster):
DE
DE
DE            g[ p[ p[i] + j ] + k ]
DE
DE
DE        Now for any point in 3-space,  we just  have  to  do  the
DE        following two steps:
DE
DE        (1) Get the  gradient  for  each  of  its  surrounding  8
DE        integer lattice points as above.
DE
DE        (2) Do a tricubic hermite  spline  interpolation,  giving
DE        each lattice point the value 0.0.
DE
DE        The second step above is just an evaluation of  the  her-
DE        mite  derivative basis function 3 * t * t - 2 * t * t * t
DE        in each by a dot product of the gradient at the lattice.
DE
DE        Here is my implementation in C  of  the  noise  function.
DE        Feel  free  to use it, as long as you reference where you
DE        got it. :-)
DE
DE
DE
DE/* noise function over R3 - implemented by a pseudorandom tricubic spline */
DE
DE#include <stdio.h>
DE#include <math.h>
DE
DE#define DOT(a,b) (a[0] * b[0] + a[1] * b[1] + a[2] * b[2])
DE
DE#define B 256
DE
DEstatic p[B + B + 2];
DEstatic float g[B + B + 2][3];
DEstatic start = 1;
DE
DE#define setup(i,b0,b1,r0,r1) \
DE        t = vec[i] + 10000.; \
DE        b0 = ((int)t) & (B-1); \
DE        b1 = (b0+1) & (B-1); \
DE        r0 = t - (int)t; \
DE        r1 = r0 - 1.;
DE
DEfloat noise3(vec)
DEfloat vec[3];
DE{
DE        int bx0, bx1, by0, by1, bz0, bz1, b00, b10, b01, b11;
DE        float rx0, rx1, ry0, ry1, rz0, rz1, *q, sy, sz, a, b, c, d, t, u, v;
DE        register i, j;
DE
DE        if (start) {
DE                start = 0;
DE                init();
DE        }
DE
DE        setup(0, bx0,bx1, rx0,rx1);
DE        setup(1, by0,by1, ry0,ry1);
DE        setup(2, bz0,bz1, rz0,rz1);
DE
DE        i = p[ bx0 ];
DE        j = p[ bx1 ];
DE
DE        b00 = p[ i + by0 ];
DE        b10 = p[ j + by0 ];
DE        b01 = p[ i + by1 ];
DE        b11 = p[ j + by1 ];
DE
DE#define at(rx,ry,rz) ( rx * q[0] + ry * q[1] + rz * q[2] )
DE
DE#define surve(t) ( t * t * (3. - 2. * t) )
DE
DE#define lerp(t, a, b) ( a + t * (b - a) )
DE
DE        sx = surve(rx0);
DE        sy = surve(ry0);
DE        sz = surve(rz0);
DE
DE
DE        q = g[ b00 + bz0 ] ; u = at(rx0,ry0,rz0);
DE        q = g[ b10 + bz0 ] ; v = at(rx1,ry0,rz0);
DE        a = lerp(sx, u, v);
DE
DE        q = g[ b01 + bz0 ] ; u = at(rx0,ry1,rz0);
DE        q = g[ b11 + bz0 ] ; v = at(rx1,ry1,rz0);
DE        b = lerp(sx, u, v);
DE
DE        c = lerp(sy, a, b);          /* interpolate in y at lo x */
DE
DE        q = g[ b00 + bz1 ] ; u = at(rx0,ry0,rz1);
DE        q = g[ b10 + bz1 ] ; v = at(rx1,ry0,rz1);
DE        a = lerp(sx, u, v);
DE
DE        q = g[ b01 + bz1 ] ; u = at(rx0,ry1,rz1);
DE        q = g[ b11 + bz1 ] ; v = at(rx1,ry1,rz1);
DE        b = lerp(sx, u, v);
DE
DE        d = lerp(sy, a, b);          /* interpolate in y at hi x */
DE
DE        return 1.5 * lerp(sz, c, d); /* interpolate in z */
DE}
DE
DEstatic init()
DE{
DE        long random();
DE        int i, j, k;
DE        float v[3], s;
DE
DE/* Create an array of random gradient vectors uniformly on the unit sphere */
DE
DE        srandom(1);
DE        for (i = 0 ; i < B ; i++) {
DE                do {                            /* Choose uniformly in a cube */
DE                        for (j=0 ; j<3 ; j++)
DE                                v[j] = (float)((random() % (B + B)) - B) / B;
DE                        s = DOT(v,v);
DE                } while (s > 1.0);              /* If not in sphere try again */
DE                s = sqrt(s);
DE                for (j = 0 ; j < 3 ; j++)       /* Else normalize */
DE                        g[i][j] = v[j] / s;
DE        }
DE
DE/* Create a pseudorandom permutation of [1..B] */
DE
DE        for (i = 0 ; i < B ; i++)
DE                p[i] = i;
DE        for (i = B ; i > 0 ; i -= 2) {
DE                k = p[i];
DE                p[i] = p[j = random() % B];
DE                p[j] = k;
DE        }
DE
DE/* Extend g and p arrays to allow for faster indexing */
DE
DE        for (i = 0 ; i < B + 2 ; i++) {
DE                p[B + i] = p[i];
DE                for (j = 0 ; j < 3 ; j++)
DE                        g[B + i][j] = g[i][j];
DE        }
DE}
DE
DE
SHAR_EOF
fi
exit 0
#	End of shell archive