{ THIS CODE WAS WRITTEN BY CLAUDE } { -*- delphi -*- }

// Starting prompt was something vaguely equivalent to "If I have a
// uniform granular mixture of materials, each material consisting of
// an sRGB color, bond albedo, and volume per particle, how would one
// compute the mean color of the pile?"

// Requires: uses Math; (for Power, Min, Max)

function TOreProfile.ComputeColor(Encyclopedia: TMaterialEncyclopedia): TMaterialColor;
const
   MaxReflectance = 0.999; // keeps the Hapke inversion finite
   LumaR = 0.2126;
   LumaG = 0.7152;
   LumaB = 0.0722;

   // the places below where we ignore range errors are because of implied Double/Extended conversions

   function SRGBToLinear(Value: Byte): Double; inline;
   var
      C: Double;
   begin
      C := Value / 255.0;
      if (C <= 0.04045) then
         Result := C / 12.92 // $R-
      else
         Result := Power((C + 0.055) / 1.055, 2.4); // $R-
   end;

   function LinearToSRGB(Value: Double): Byte; inline;
   var
      C: Double;
   begin
      if (Value <= 0.0) then
         exit(0);
      if (Value >= 1.0) then
         exit(255);
      if (Value <= 0.0031308) then
         C := Value * 12.92 // $R-
      else
         C := 1.055 * Power(Value, 1.0 / 2.4) - 0.055; // $R-
      Result := Round(C * 255.0); // $R-
   end;

   // Hapke, isotropic scatterers: R = (1 - g) / (1 + g), g = sqrt(1 - w)
   function ReflectanceToSingleScatteringAlbedo(R: Double): Double; inline;
   var
      Gamma: Double;
   begin
      Gamma := (1.0 - R) / (1.0 + R);
      Result := 1.0 - Gamma * Gamma;
   end;

   function SingleScatteringAlbedoToReflectance(W: Double): Double; inline;
   var
      Gamma: Double;
   begin
      Gamma := Sqrt(Max(0.0, 1.0 - W));
      Result := (1.0 - Gamma) / (1.0 + Gamma);
   end;

var
   Ore: TOres;
   Material: TMaterial;
   VolumeFraction, Weight, TotalWeight: Double;
   Luminance, Scale, ScaleSum, ScaleWeight: Double;
   Reflectance: array[0..2] of Double;
   SSASum: array[0..2] of Double;
   Output: array[0..2] of Byte;
   Channel: Integer;
begin
   Result := $000000;
   if (not HasOre) then
      exit;

   TotalWeight := 0.0;
   ScaleSum := 0.0;
   ScaleWeight := 0.0;
   for Channel := 0 to 2 do
      SSASum[Channel] := 0.0;

   for Ore := Low(TOres) to High(TOres) do
   begin
      VolumeFraction := OreFraction[Ore].ToDouble();
      if (VolumeFraction <= 0.0) then
         continue;
      Material := Encyclopedia.Materials[Ore];

      // Visible cross-section per unit volume scales as (volume per particle)^(-1/3).
      Assert(Material.VolumePerUnit.AsDouble > 0.0);
      Weight := VolumeFraction / Power(Material.VolumePerUnit.AsDouble, 1.0 / 3.0); // $R-

      // Color gives chromaticity; bond albedo gives brightness.
      Reflectance[0] := SRGBToLinear((Material.Color shr 16) and $FF); // $R-
      Reflectance[1] := SRGBToLinear((Material.Color shr 8) and $FF); // $R-
      Reflectance[2] := SRGBToLinear(Material.Color and $FF); // $R-
      Luminance := LumaR * Reflectance[0] + LumaG * Reflectance[1] + LumaB * Reflectance[2]; // $R-
      if (Luminance > 0.0) then
      begin
         Scale := Material.BondAlbedo / Luminance;
         for Channel := 0 to 2 do
            Reflectance[Channel] := Min(Reflectance[Channel] * Scale, MaxReflectance); // $R-
         // Track the color->reflectance scale so we can map back to color space at the end.
         ScaleSum := ScaleSum + Weight * Scale;
         ScaleWeight := ScaleWeight + Weight;
      end;
      // (A pure-black color contributes zero reflectance and no scale.)

      for Channel := 0 to 2 do
         SSASum[Channel] := SSASum[Channel] + Weight * ReflectanceToSingleScatteringAlbedo(Reflectance[Channel]);
      TotalWeight := TotalWeight + Weight;
   end;

   if ((TotalWeight <= 0.0) or (ScaleWeight <= 0.0) or (ScaleSum <= 0.0)) then
      exit;
   Scale := ScaleSum / ScaleWeight;

   for Channel := 0 to 2 do
      Output[Channel] := LinearToSRGB(SingleScatteringAlbedoToReflectance(SSASum[Channel] / TotalWeight) / Scale);

   Result := (Output[0] shl 16) or (Output[1] shl 8) or Output[2]; // $R-
end;