thinfilm_multilayer_oslpy.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_oslpy.h"

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

float cgf(Film film, float phase) {
    return (phase*phase * film.v * .5);
}

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

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

void thin_film(vector lambda, float phase_a, float phase_b, float alpha_a, float alpha_b, output vector R)
{
    vector l = (1.0/lambda);

    vector _phase_a = phase_a*l;
    vector _phase_b = phase_b*l;

    vector _alpha_a = alpha_a*l*l;
    vector _alpha_b = alpha_b*l*l;

    vector phi = exp(-2.0*(_alpha_a + _alpha_b));
    vector psi = exp(-1.0*(_alpha_a + _alpha_b));

    vector x = (exp(-_alpha_a) + exp(-_alpha_b))/(1.0 + psi);
    vector y = (exp(-_alpha_a) - exp(-_alpha_b))/(1.0 - psi);

    vector _M = (_phase_a + _phase_b) * .5;
    vector _N = (_phase_a - _phase_b) * .5;

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

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

    vector 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 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 Thickness (nm)",
        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 (nm)",
        float min = 0,
        float max = 5000
    ]],

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

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

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

/* -------------------------------------------------------------------------- */
/*                                   Outputs                                  */
/* -------------------------------------------------------------------------- */
    
    output color R = 0
    [[
		string help = ""
	]],
    output color T = 0
    [[
		string help = ""
	]]
)
{   
    Film film_a = Film(film_0_d, pow(film_0_sd, 2.0) / pow(film_0_d, 2.0), film_0_n);
    Film film_b = Film(film_1_d, pow(film_1_sd, 2.0) / pow(film_1_d, 2.0), film_1_n);

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

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

    float Kr, Kt;
    vector vR, vT;
    fresnel(I, Normal, 1.0/eta, Kr, Kt, vR, vT);

    vector lambda_samples[ARRAY_SIZE];
    gen_lambda_samples(P, lambda_samples);

    vector spec_R[ARRAY_SIZE];

    float phase_a = phase(film_a, cos_theta);
    float phase_b = phase(film_b, cos_theta);

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

    thin_film(lambda_samples[0], phase_a, phase_b, alpha_a, alpha_b, spec_R[0]);
    thin_film(lambda_samples[1], phase_a, phase_b, alpha_a, alpha_b, spec_R[1]);
    thin_film(lambda_samples[2], phase_a, phase_b, alpha_a, alpha_b, spec_R[2]);
    thin_film(lambda_samples[3], phase_a, phase_b, alpha_a, alpha_b, spec_R[3]);

    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 = Kr + ( (1.0 - Kr) * SPEC_to_sRGB(spec_R, lambda_samples) ) * lamella_falloff;
    color rgb_T = 1.0 - rgb_R;

    R = rgb_R;
    T = rgb_T;
}