//Thinfilm shader cobbled together from:
//https://docs.chaosgroup.com/display/OSLShaders/Thin+Film+Shader
//and:
//https://www.shadertoy.com/view/ls2Bz1
//copied together by your's truly, Entagma - Dec. 2020
//https://entagma.com/holiday-giveaway-mantra-thinfilm-shader/
//sequentially copied and modified to use with 3delight by Faitel - Dec. 2022

#include "stdosl.h"

surface thinfilm
(
  color col = color(0)
  [[  string help = "Constant color",
      string label = "Color",
      string widget = "number"]],
  float fintensity = 0
  [[  string help = "Fresnel reduction multiplier",
      string label = "Intensity",
      string page = "Fresnel reduction",
      string widget = "number",
      float min = 0
  ]],
    float fior = 1.4
  [[  string help = "Index of reflection for fresnel reduction",
      string label = "Ior",
      string page = "Fresnel reduction",
      string widget = "number",
      float min = 0
  ]],
  int sbiter = 8
  [[  string help = "Number of spectral bands. 8 seems fine, 32 is excessive.",
      string label = "Iterations",
      string page = "Spectral Bands",
      string widget = "number",
  ]],
  int spfunc = 0
  [[
      string label = "Function",
      string page = "Spectral Bands",
      string widget = "mapper",
      string options =
        "Gems:0|"
        "Zucconi:1|"
        "Bruton:2|"
  ]],
  float thickness = 1
  [[  string help = "Thin Film Thickness",
      string label = "Thickness",
      string page = "Thin Film",
      string widget = "number",
      float min = 0,
      float max = 1
  ]],
  float thicknessMin = 100
  [[  string help = "Minimal Thin Film Thickness",
      string label = "Thickness Min",
      string page = "Thin Film",
      string widget = "number",
      float min = 0,
      float max = 1000
  ]],
  float thicknessMax = 900
  [[  string help = "Maximal Thin Film Thickness",
      string label = "Thickness Max",
      string page = "Thin Film",
      string widget = "number",
      float min = 0,
      float max = 1000
  ]],
  float nmedium = 1
  [[  string help = "approximate refractive index of air",
      string label = "nmedium",
      string page = "Thin Film",
      string widget = "number",
      float min = 0,
      float max = 10
  ]],
  float nfilm = 1.4
  [[  string help = "approximate refractive index of water",
      string label = "nfilm",
      string page = "Thin Film",
      string widget = "number",
      float min = 0,
      float max = 10
  ]],
  float ninternal = 1
  [[  string help = "approximate refractive index of the lower material", 
      string label = "ninternal",
      string page = "Thin Film",
      string widget = "number",
      float min = 0,
      float max = 10
  ]],

  output color outColor = 0
)
{

float saturate(float x)
{
    return min(1.0, max(0.0,x));
}
vector saturate (vector x)
{
    return min(color(1.0,1.0,1.0), max(color(0.0,0.0,0.0),x));
}
// --- Spectral Zucconi --------------------------------------------
// By Alan Zucconi
// Based on GPU Gems: https://developer.nvidia.com/sites/all/modules/custom/gpugems/books/GPUGems/gpugems_ch08.html
// But with values optimised to match as close as possible the visible spectrum
// Fits this: https://commons.wikimedia.org/wiki/File:Linear_visible_spectrum.svg
// With weighter MSE (RGB weights: 0.3, 0.59, 0.11)
vector bump3y(vector x, vector yoffset)
{
	vector y = color(1.0,1.0,1.0) - x * x;
	y = saturate(y-yoffset);
	return y;
}
vector spectral_zucconi(float w)
{
    // w: [400, 700]
	// x: [0,   1]
	float x = saturate((w - 400.0)/ 300.0);

	vector cs = color(3.54541723, 2.86670055, 2.29421995);
	vector xs = color(0.69548916, 0.49416934, 0.28269708);
	vector ys = color(0.02320775, 0.15936245, 0.53520021);

	return bump3y (	cs * (x - xs), ys);
}

// --- GPU Gems -------------------------------------------------------
// https://developer.nvidia.com/sites/all/modules/custom/gpugems/books/GPUGems/gpugems_ch08.html
vector bump3(vector x) {
        vector y = color(1.0,1.0,1.0) - x * x;
        y = max(y, {0.,0.,0.});
        return y;
}

vector spectral_gems (float w) {
    // w: [400, 700]
    // x: [0,   1]
    float x = clamp((w - 400.0)/ 300.0, 0., 1.);
    float r = 4. * (x - 0.75);
    float g = 4. * (x - 0.5);
    float b = 4. * (x - 0.25);
    vector rgb = {r,g,b};
    return bump3(rgb);
}

// --- Approximate RGB values for Visible Wavelengths ------------------
// by Dan Bruton
// http://www.physics.sfasu.edu/astro/color/spectra.html
// https://stackoverflow.com/questions/3407942/rgb-values-of-visible-spectrum
vector spectral_bruton(float w)
{
	vector c;

	if (w >= 380. && w < 440.)
		c = color
		(
			-(w - 440.) / (440. - 380.),
			0.0,
			1.0
		);
	else if (w >= 440. && w < 490.)
		c = color
		(
			0.0,
			(w - 440.) / (490. - 440.),
			1.0
		);
	else if (w >= 490. && w < 510.)
		c = color
		(	0.0,
			1.0,
			-(w - 510.) / (510. - 490.)
		);
	else if (w >= 510. && w < 580.)
		c = color
		(
			(w - 510.) / (580. - 510.),
			1.0,
			0.0
		);
	else if (w >= 580. && w < 645.)
		c = color
		(
			1.0,
			-(w - 645.) / (645. - 580.),
			0.0
		);
	else if (w >= 645. && w <= 780.)
		c = color
		(	1.0,
			0.0,
			0.0
		);
	else
		c = color
		(	0.0,
			0.0,
			0.0
		);

	return saturate(c);
}

//--- OSL Thinfilm Shader ---
/* Amplitude reflection coefficient (s-polarized) */
float rs(float n1, float n2, float cosI, float cosT) {
    return (n1 * cosI - n2 * cosT) / (n1 * cosI + n2 * cosT);
}

/* Amplitude reflection coefficient (p-polarized) */
float rp(float n1, float n2, float cosI, float cosT) {
    return (n2 * cosI - n1 * cosT) / (n1 * cosT + n2 * cosI);
}

/* Amplitude transmission coefficient (s-polarized) */
float ts(float n1, float n2, float cosI, float cosT) {
    return 2 * n1 * cosI / (n1 * cosI + n2 * cosT);
}

/* Amplitude transmission coefficient (p-polarized) */
float tp(float n1, float n2, float cosI, float cosT) {
    return 2 * n1 * cosI / (n1 * cosT + n2 * cosI);
}

// cosI is the cosine of the incident angle, that is, cos0 = dot(view angle, normal)
// lambda is the wavelength of the incident light (e.g. lambda = 510 for green)
// From http://www.gamedev.net/page/resources/_/technical/graphics-programming-and-theory/thin-film-interference-for-computer-graphics-r2962
float thinFilmReflectance(float cos0, float lambda, float thickness, float n0, float n1, float n2) {

    float PI=M_PI; //Pi osl constant

    // compute the phase change term (constant)
    float d10 = (n1 > n0) ? 0 : PI;
    float d12 = (n1 > n2) ? 0 : PI;
    float delta = d10 + d12;

    // now, compute cos1, the cosine of the reflected angle
    float sin1 = pow(n0 / n1, 2) * (1 - pow(cos0, 2));
    if (sin1 > 1) return 1.0; // total internal reflection
    float cos1 = sqrt(1 - sin1);

    // compute cos2, the cosine of the final transmitted angle, i.e. cos(theta_2)
    // we need this angle for the Fresnel terms at the bottom interface
    float sin2 = pow(n0 / n2, 2) * (1 - pow(cos0, 2));
    if (sin2 > 1) return 1.0; // total internal reflection
    float cos2 = sqrt(1 - sin2);

    // get the reflection transmission amplitude Fresnel coefficients
    float alpha_s = rs(n1, n0, cos1, cos0) * rs(n1, n2, cos1, cos2); // rho_10 * rho_12 (s-polarized)
    float alpha_p = rp(n1, n0, cos1, cos0) * rp(n1, n2, cos1, cos2); // rho_10 * rho_12 (p-polarized)

    float beta_s = ts(n0, n1, cos0, cos1) * ts(n1, n2, cos1, cos2); // tau_01 * tau_12 (s-polarized)
    float beta_p = tp(n0, n1, cos0, cos1) * tp(n1, n2, cos1, cos2); // tau_01 * tau_12 (p-polarized)

    // compute the phase term (phi)
    float phi = (2 * PI / lambda) * (2 * n1 * thickness * cos1) + delta;

    // finally, evaluate the transmitted intensity for the two possible polarizations
    float ts = pow(beta_s, 2) / (pow(alpha_s, 2) - 2 * alpha_s * cos(phi) + 1);
    float tp = pow(beta_p, 2) / (pow(alpha_p, 2) - 2 * alpha_p * cos(phi) + 1);

    // we need to take into account conservation of energy for transmission
    float beamRatio = (n2 * cos2) / (n0 * cos0);

    // calculate the average transmitted intensity (if you know the polarization distribution of your
    // light source, you should specify it here. if you don't, a 50%/50% average is generally used)
    float myt = beamRatio * (ts + tp) / 2;

    // and finally, derive the reflected intensity
    return 1 - myt;
}

float cos0 = abs(dot(I , N));

float myt = clamp(thickness, 0., 1.);
float thick=thicknessMin*(1.0-myt)+thicknessMax*myt;

color outcol = col; //get color;
float fint = fintensity; //get fresnel multiplier
float fnfilm = fior; //get fresnel ior

int maxiter = sbiter; //number of spectral bands. 8 seems fine, 32 is excessive.
for(int i = 0; i < maxiter; i++)
{
    float lambda = 400. + (300. * (float(i)/float(maxiter))); //Spectral bands ranging from 400nm to 700nm
    float ref = thinFilmReflectance(cos0, lambda, thick, nmedium, nfilm, ninternal);

    //we don't need no fresnel, so we need to divide it out again
    //I know - being cheap, but atm just not in the mood for algebra
    float kr, kt;
    fresnel(-I, N, fnfilm, kr, kt);
    kr = kr;
    if (fint == 0 ) ref = ref;
    if (fint > 0 ) ref = ref / (kr * fint);

    vector rgb;
    if (spfunc == 0) rgb = spectral_gems(lambda);
    if (spfunc == 1) rgb = spectral_zucconi(lambda);
    if (spfunc == 2) rgb = spectral_bruton(lambda);

    rgb = rgb * ref;
    outcol = outcol + rgb;
}

outColor = outcol / float(maxiter);
}