CollimatedSpecularRadiance: PUBLIC PROCEDURE [surfaceNormal: Vector, lightSource: LightSource, albedo: REAL ← 1] RETURNS [radiance: REAL] = {radiance ← 180};
 -- Collimated light source at lightSource. Specular surface with given albedo oriented by given surfaceNormal.

ElementWiseProduct: PRIVATE PROC [v1, v2: Vector] RETURNS [prod: Vector] = {
 prod[1] ← v1[1] * v2[1];
 prod[2] ← v1[2] * v2[2];
 prod[3] ← v1[3] * v2[3];
 };

DiffuseReflectance: PUBLIC PROCEDURE [surfaceNormal: Vector,
                 surfacePt: Point3d, surfaceColor: Color,
                 lightSources: LightSourceList]
          RETURNS [completeColor: Color] =
{
 OPEN SVVector3d;
 r, g, b, intensity, ambientScalar: REAL;
 scaleFactor: REAL;
 surfaceToLight, surfaceClr, surfaceClrTemp, liteClr: Vector;
 colorSum: Vector ← [0,0,0];

 [r,g,b] ← GraphicsColor.ColorToRGB[surfaceColor];
 surfaceClr ← [r,g,b];
 surfaceNormal ← Normalize[surfaceNormal];  -- all vectors must be kept normalized

 ambientScalar ← DotProduct[surfaceNormal,[0,1,0]];  -- overhead ambient light
 IF ambientScalar > 0 THEN
  colorSum ← Scale[surfaceClr, 0.25 * DotProduct[surfaceNormal,[0,1,0]]];
  -- colorSum ← Add[colorSum, Scale[surfaceClr, 0.25 ... ] is clearer but slower.  

 FOR lite: LightSourceList ← lightSources, lite.rest UNTIL lite = NIL  -- for each light source
  DO
  [r,g,b] ← GraphicsColor.ColorToRGB[lite.first.color];
 liteClr ← [r,g,b];
 surfaceToLight ← Normalize[Difference[lite.first.position, surfacePt]];
  -- vector from surface to light source
 intensity ← DotProduct[surfaceNormal, surfaceToLight];  -- lambertian reflectance
 IF intensity < 0.0 THEN intensity ← 0.0;      -- no negative reflectance allowed
 surfaceClrTemp ← Scale[surfaceClr, intensity];  -- scale surface color by reflectance
 surfaceClrTemp ← ElementWiseProduct[surfaceClrTemp, liteClr];
 
  colorSum ← Add[colorSum, surfaceClrTemp];   -- sum reflectances for all lights
  ENDLOOP;

 scaleFactor ← MAX[colorSum[1], colorSum[2], colorSum[3]]; -- scale sum to max. value to avoid
 IF scaleFactor > 1.0 THEN colorSum ← Scale[colorSum, 1.0/scaleFactor];
  -- desaturation at bright spots
 r ← colorSum[1];
 g ← colorSum[2];
 b ← colorSum[3];
 completeColor ← GraphicsColor.RGBToColor[r,g,b];
};


DiffuseAndSpecularReflectance: PUBLIC PROCEDURE [eyePoint: Point3d, surfaceNormal: Vector,
                 surfacePt: Point3d, surfaceColor: Color,
                 lightSources: LightSourceList]
          RETURNS [completeColor: Color] =
{
 OPEN SVVector3d;
 r, g, b, intensity, hilitValue, ambientScalar: REAL;
 scaleFactor: REAL;
 surfaceToLight, surfaceToEye, surfaceClr, surfaceClrTemp, finalClr, halfAngle, liteClr: Vector; colorSum: Vector ← [0,0,0];

 [r,g,b] ← GraphicsColor.ColorToRGB[surfaceColor];
 surfaceClr ← [r,g,b];
 surfaceNormal ← Normalize[surfaceNormal];           -- all vectors must be kept normalized

 ambientScalar ← DotProduct[surfaceNormal,[0,1,0]];  -- overhead ambient light
 IF ambientScalar > 0 THEN
  colorSum ← Scale[surfaceClr, 0.25 * DotProduct[surfaceNormal,[0,1,0]]];  -- colorSum ← Add[colorSum, Scale[surfaceClr, 0.25 ... ] is clearer but slower.  

 FOR lite: LightSourceList ← lightSources, lite.rest UNTIL lite = NIL       -- for each light source
  DO
  [r,g,b] ← GraphicsColor.ColorToRGB[lite.first.color];
 liteClr ← [r,g,b];
 surfaceToLight ← Normalize[Difference[lite.first.position, surfacePt]]; -- vector from surface to light source
 intensity ← DotProduct[surfaceNormal, surfaceToLight];         -- lambertian reflectance
 IF intensity < 0.0 THEN intensity ← 0.0;             -- no negative reflectance allowed
 surfaceClrTemp ← Scale[surfaceClr, intensity];         -- scale surface color by reflectance
 surfaceClrTemp ← ElementWiseProduct[surfaceClrTemp, liteClr];
 
 surfaceToEye ← Normalize[Difference[eyePoint, surfacePt]];
 halfAngle ← Normalize[Add[surfaceToLight, surfaceToEye]];
 hilitValue ← DotProduct[surfaceNormal, halfAngle];         -- specular reflectance
  IF hilitValue < 0.95 THEN hilitValue ← 0.0
   ELSE hilitValue ← RealFns.Power[hilitValue, 100];          -- all surfaces have same specular power

  IF hilitValue > .01
 THEN                       -- add in light source color for highlight
   finalClr ← Add[surfaceClrTemp, Scale[Difference[liteClr, surfaceClrTemp], hilitValue]]
 ELSE finalClr ← surfaceClrTemp;
  
  colorSum ← Add[colorSum, finalClr];            -- sum reflectances for all lights
  ENDLOOP;

 scaleFactor ← MAX[colorSum[1], colorSum[2], colorSum[3]];        -- scale sum to max. value to avoid
 IF scaleFactor > 1.0 THEN colorSum ← Scale[colorSum, 1.0/scaleFactor];   -- desaturation at bright spots
 r ← colorSum[1];
 g ← colorSum[2];
 b ← colorSum[3];
 completeColor ← GraphicsColor.RGBToColor[r,g,b];
};


Power: PRIVATE PROC [base: REAL, exponent: NAT] RETURNS [baseToExp: REAL] = {
 baseToExp ← 1.0;
 FOR i: NAT IN [1..exponent] DO
  baseToExp ← baseToExp * base;
 ENDLOOP;
 };

SimpleFilmIrradiance: PUBLIC PROCEDURE [radiance: REAL, omegaAperture: REAL] RETURNS [irradiance: REAL] = {
 irradiance ← radiance*omegaAperture;
 };

SolidAngle: PUBLIC PROCEDURE [surfacePoint: Point3d, areaAperture: REAL] RETURNS [omegaAperture: REAL] = {
 -- solid angle is just Acos(theta)/r^2, where r is the distance from surface point to aperture (at [0,0,0], A is the areaAperture, and theta is the angle between the aperture normal [0,0,-1] and the line from surfacePoint to [0,0,0];
 fromApertureToSurface: Vector ← surfacePoint; -- same as surfacePoint because aperture at [0,0,0]
 r: REAL ← SVVector3d.Magnitude[fromApertureToSurface];
 cosTheta: REAL ← -fromApertureToSurface[3]/r;
 omegaAperture ← (areaAperture*cosTheta)/(r*r);
 };

SpecFilmIrradiance: PUBLIC PROCEDURE [surfaceNormal: Vector, lightSource: LightSource, albedo: REAL ← 1, omegaAperture: REAL] RETURNS [irradiance: REAL] = {

 -- A convenience routine which combines CollimatedSpecularRadiance and SimpleFilmIrradiance.

 radiance: REAL ← CollimatedSpecularRadiance[surfaceNormal,lightSource,albedo];
 irradiance ← SimpleFilmIrradiance[radiance,omegaAperture];
 };
 
END.