mx_blackbody.osl

void mx_blackbody(float temperature, output color color_value)
{
    float xc, yc;
    float t, t2, t3, xc2, xc3;

    t = 1000.0 / temperature;
    t2 = t * t;
    t3 = t * t * t;

    // if value outside valid range of approximation clamp to accepted temperature range
    clamp(temperature, 1667.0, 25000.0);


    // Cubic spline approximation for Kelvin temperature to sRGB conversion
    // (https://en.wikipedia.org/wiki/Planckian_locus#Approximation)
    if (temperature < 4000.0) {  // 1667K <= temperature < 4000K
      xc = -0.2661239 * t3 - 0.2343580 * t2 + 0.8776956 * t + 0.179910;
    }
    else {  // 4000K <= temperature <= 25000K
      xc = -3.0258469 * t3 + 2.1070379 * t2 + 0.2226347 * t + 0.240390;
    }
    xc2 = xc * xc;
    xc3 = xc * xc * xc;

    if (temperature < 2222.0) {  // 1667K <= temperature < 2222K
      yc = -1.1063814 * xc3 - 1.34811020 * xc2 + 2.18555832 * xc - 0.20219683;
    }
    else if (temperature < 4000.0) {  // 2222K <= temperature < 4000K
      yc = -0.9549476 * xc3 - 1.37418593 * xc2 + 2.09137015 * xc - 0.16748867;
    }
    else {  // 4000K <= temperature <= 25000K
      yc = 3.0817580 * xc3 - 5.87338670 * xc2 + 3.75112997 * xc - 0.37001483;
    }

    if (yc <= 0.0) {  // avoid division by zero
      color_value = color(1.0);
      return;
    }

    vector XYZ = vector(xc / yc, 1.0, (1 - xc - yc) / yc);

    /// XYZ to Rec.709 RGB colorspace conversion
    matrix XYZ_to_RGB = matrix( 3.2406, -0.9689, 0.0557,
                               -1.5372, 1.8758, -0.2040,
                               -0.4986, 0.0415, 1.0570);

    color_value = transform(XYZ_to_RGB, XYZ);
    color_value = max(color_value, vector(0.0));
}