thinfilm_multilayer.osl

/* ------------------------------------------------------------------------- *
*
*    Copyright (C) 2024 Jake Kurtz
*
*    This program is free software: you can redistribute it and/or modify
*    it under the terms of the GNU General Public License as published by
*    the Free Software Foundation, either version 3 of the License, or
*    (at your option) any later version.
*
*    This program is distributed in the hope that it will be useful,
*    but WITHOUT ANY WARRANTY; without even the implied warranty of
*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
*    GNU General Public License for more details.
*
*    You should have received a copy of the GNU General Public License
*    along with this program. If not, see <https://www.gnu.org/licenses/>.
*
* ------------------------------------------------------------------------- */

#include "spectral.h"

struct Film
{
    float d; // Film thickness (nm)
    float v; // Film thickness variance (nm)
    float n; // IOR
};

float cgf(Film film, float phase) {
    return (pow2(phase) * film.v / (pow2(film.d)*2.0));
}

float phase(Film film, float lambda, float cos_theta) 
{
    float _sin_theta = (1.0 - (pow2(cos_theta))) / pow(film.n, 2.0);
    float _cos_theta = sqrt(1.0 - _sin_theta);

    return (M_4PI * film.n * film.d * _cos_theta) / lambda;
}

void thin_film(float lambda, Film film_a, Film film_b, float cos_theta, output float R)
{
    float phase_a = phase(film_a, lambda, cos_theta);
    float phase_b = phase(film_b, lambda, cos_theta);

    float alpha_a = cgf(film_a, phase_a);
    float alpha_b = cgf(film_b, phase_b);

    float phi = exp(-2.0*(alpha_a + alpha_b));
    float psi = exp(-1.0*(alpha_a + alpha_b));

    float x = (exp(-alpha_a) + exp(-alpha_b))/(1.0 + psi);
    float y = (exp(-alpha_a) - exp(-alpha_b))/(1.0 - psi);

    float _M = (phase_a + phase_b) * .5;
    float _N = (phase_a - phase_b) * .5;

    float numer = 1 - cos(_M)*cos(_N)*x + sin(_M)*sin(_N)*y;
    float denom = 1 + phi - 2.0 * cos(2.0*_M) * psi;

    float intensity = (1.0 - phi) * numer / denom;

    float peak_max = 2.0 * (1+psi) / (1-psi);

    R = (intensity / peak_max);
}

shader ML_ThinFilm
(
/* -------------------------------------------------------------------------- */
/*                                   Inputs                                   */
/* -------------------------------------------------------------------------- */

    float film_0_n = 1.56
    [[ 
        string help = "",
        string label = "Film 0 IOR",
        float min = 0,
        float max = 25
    ]],
    
    float film_0_sd = 16.0
    [[ 
        string help = "",
        string label = "Film 0 Sigma (nm)",
        float min = 0,
        float max = 1000
    ]],
    
    float film_0_d = 120.0
    [[ 
        string help = "",
        string label = "Film 0 Thickness",
        float min = 0,
        float max = 5000
    ]],

    float film_1_n = 1.52
    [[ 
        string help = "",
        string label = "Film 1 IOR",
        float min = 0,
        float max = 25
    ]],
    
    float film_1_sd = 16.0
    [[ 
        string help = "",
        string label = "Film 1 Sigma (nm)",
        float min = 0,
        float max = 1000
    ]],
    
    float film_1_d = 90.0
    [[ 
        string help = "",
        string label = "Film 1 Thickness",
        float min = 0,
        float max = 5000
    ]],

    float fall_off = 15.0
    [[ 
        string help = "",
        string label = "Fall Off Speed"
    ]],

    normal Normal = N
    [[ 
        string label = "Normal"
    ]],

    normal LNormal = N
    [[ 
        string label = "Lamella Normal"
    ]],

/* -------------------------------------------------------------------------- */
/*                                   Outputs                                  */
/* -------------------------------------------------------------------------- */

    output color R = 0
    [[
		string help = ""
	]],
    output color T = 0
    [[
		string help = ""
	]]
)
{
    Film film_a = Film(film_0_d, pow2(film_0_sd), film_0_n);
    Film film_b = Film(film_1_d, pow2(film_1_sd), film_1_n);

    float cos_theta = dot(N, I);
    float Lcos_theta = dot(LNormal, I);

    float eta = backfacing() ? 1.0/film_a.n : film_a.n;

    float Fr = fresnel(cos_theta, eta);

    float lambda_samples[LAMBDA_SAMPLES];
    gen_lambda_samples(I+P, lambda_samples);

    float spec_R[LAMBDA_SAMPLES];

    for (int i = 0; i < LAMBDA_SAMPLES; i++) {
        thin_film(lambda_samples[i], film_a, film_b, Lcos_theta, spec_R[i]);
    }

    float lamella_falloff = 1.0;
    if (Lcos_theta < 0.0) {
        lamella_falloff = mix(0.0, 1.0, pow(1.0 + Lcos_theta, fall_off));
    }

    color rgb_R = Fr + ( (1.0 - Fr) * SPEC_to_sRGB(spec_R, lambda_samples) ) * lamella_falloff;
    color rgb_T = 1.0 - rgb_R;

    R = rgb_R;
    T = rgb_T;
}