ThinFilmInterference - 图1

    1. /* OSL shader to calculate thin-film interference
    2. v3: Modified : 03.23.2021 by Saul Espinosa for Redshift 3d GPU Renderer, added metadata, new labels, ACES output.
    3. v1.6, 2016 prutser. License for the *code*: CC BY-NC-SA 4.0
    4. Note: the output of the code may explicitly be used for commercial renders.
    5. v1.6: Allow double layers for more flexibility. Use the #define statements below to make the code only calculate
    6. */
    8. #include "stdosl.h"
    10. #define SINGLELAYER
    11. #define DOUBLELAYER
    12. #define NONVACUUM
    13. #define EXPOSE_WAVELENGTH
    15. #ifdef DOUBLELAYER
    16. #define SINGLELAYER
    17. #endif
    19. #ifndef NONVACUUM
    20. #define IOR_outside 1
    21. #endif
    23. #ifndef EXPOSE_WAVELENGTH
    24. #define wavelength color(0.65,0.532,0.45)
    25. #endif
    29. struct complex {
    30.         color r;
    31.         color i;
    32.   };
    36. shader thin_film
    37. [[  string help = "Advanced Multi-Layered Thin Film Interference",
    38.     string label = "Thin Film Interference" ]]
    39. (
    41. #ifdef SINGLELAYER
    42.     float thickness_T = 1
    43.     [[
    44.         string help = "thickness of layer 1 in micrometers",
    45.         string label = "Top Layer Thickness",
    46.         string page = "1 : Top Layer",
    47.         float min = 0,
    48.         float max = 25
    49.     ]],
    50.     float IOR_top = 1.6
    51.     [[
    52.         string help = "IOR for top layer",
    53.         string label = "Top Layer IOR",
    54.         string page = "1 : Top Layer",
    55.         float min = 0,
    56.         float max = 25
    57.     ]],
    58.     color k_top = 0
    59.     [[
    60.         string help = "Extinction coefficient of top layer",
    61.         string label = "Top Layer Extinction",
    62.         string page = "1 : Top Layer",
    63.         color min = 0,
    64.         color max = 1
    65.     ]],
    66. #endif
    68. #ifdef DOUBLELAYER
    69.     float thickness_M = 0
    70.     [[
    71.         string help = "thickness of layer 2 in micrometers",
    72.         string label = "Middle Layer Thickness",
    73.         string page = "2 : Middle Layer",
    74.         float min = 0,
    75.         float max = 25
    76.     ]],
    77.     float IOR_middle = 1
    78.     [[
    79.         string help = "IOR for middle layer",
    80.         string label = "Mid Layer IOR",
    81.         string page = "2 : Middle Layer",
    82.         float min = 0,
    83.         float max = 25
    84.     ]],
    85.     color k_middle = 0
    86.     [[
    87.         string help = "Extinction coefficient of top layer",
    88.         string label = "Mid Layer Extinction",
    89.         string page = "2 : Middle Layer",
    90.         color min = 0,
    91.         color max = 1
    92.     ]],
    93. #endif
    95. // Substrate Parameters
    97.     float IOR_base_material = 1
    98.     [[
    99.         string help = "IOR for Substrate Base layer",
    100.         string label = "Base Layer IOR",
    101.         string page = "3 : Base Layer",
    102.         float min = 0,
    103.         float max = 25
    104.     ]],
    105.     color k_base_material = 0
    106.     [[
    107.         string help = "Extinction coefficient of substrate layer",
    108.         string label = "Base Layer Extinction",
    109.         string page = "3 : Base Layer",
    110.         color min = 0,
    111.         color max = 1
    112.     ]],
    114. #ifdef EXPOSE_WAVELENGTH
    115.     color wavelength = color(0.65,0.532,0.45)
    116.     [[
    117.         string help = "RGB wavelengths to use for the interference",
    118.         string label = "Wavelength",
    119.         string page = "4 : Extra",
    120.         color min = 0,
    121.         color max = 1
    122.     ]],
    123. #endif
    125. #ifdef NONVACUUM
    126.     color IOR_outside = 1
    127.     [[
    128.         string help = "the index of refraction of the ‘outside world’. Typically 1, but under water would be 1.33, for example",
    129.         string label = "IOR of Environment",
    130.         string page = "4 : Extra",
    131.         color min = 0,
    132.         color max = 3
    133.     ]],
    135.     normal Normal = N
    136.     [[
    137.     string page = "4 : Extra"
    138.     ]],
    140.     int   aces = 0
    141.     [[ string widget = "checkBox",
    142.     string label = "ACES",
    143.     string page = "4 : Extra",
    144.     int connectable = 0 ]],
    146. #endif
    148. // Outputs
    150.     output color outColor = 0
    151.     [[
    152.         string help = "outColor: Reflected percentage of the light, per color (RGB), at the angle between the normal and the ray of light"
    153.     ]],
    154.     output color outColor_N = 0
    155.     [[
    156.         string help = "outColor_N: reflected percentage of the light, per color (RGB), at normal incidence",
    157.     ]],
    158.     output color Transmission = 1
    159.         [[
    160.         string help = "Transmission: transmitted percentage of light, per color. Can be hooked up to a refractive shader"
    161.     ]],
    162.     output color Trans_N = 1
    163.     [[
    164.         string help = "Trans_N: transmitted percentage of light at normal incidence"
    165.     ]]
    167.     )
    169. {
    171.         // ACES sRGB Transform
    172.         matrix aces_tm = matrix(
    173.         0.6131, 0.0701, 0.0206, 0,
    174.         0.3395, 0.9164, 0.1096, 0,
    175.         0.0474, 0.0135, 0.8698, 0,
    176.         0, 0, 0, 1);
    178.     void cmult(complex x, complex y, output complex z)
    179.     {
    180.       complex tmp; //pass by reference issues, if in reality x == z
    181.       tmp.r = x.r*y.r-x.i*y.i;
    182.       tmp.i = x.r*y.i + x.i*y.r;
    183.       z.r = tmp.r; z.i = tmp.i;
    184.     }
    186.     void cmult(color x, complex y, output complex z)
    187.     {
    188.       z.r = x*y.r;
    189.       z.i = x*y.i;
    190.     }
    192.     void cmult(complex x, color y, output complex z)
    193.     {
    194.       z.r = x.r*y;
    195.       z.i = x.i*y;
    196.     }
    198.     void cdiv(complex x, complex y, output complex z)
    199.     {
    200.       complex ystar;
    201.       ystar.r = y.r;
    202.       ystar.i = -y.i;
    204.       cmult(x,ystar,z); 
    205.       z.r = z.r/color((ystar.r*ystar.r+ystar.i*ystar.i));
    206.       z.i = z.i/color((ystar.r*ystar.r+ystar.i*ystar.i));
    207.     }
    209.     void cadd(complex x, complex y, output complex z)
    210.     {
    211.       z.r = x.r + y.r;
    212.       z.i = x.i + y.i;
    213.     }
    215.     void cadd(complex x, color y, output complex z)
    216.     {
    217.       z.r = x.r + y;
    218.       z.i = x.i;
    219.     }
    221.     void cadd(color x, complex y, output complex z)
    222.     {
    223.       z.r = x + y.r;
    224.       z.i = y.i;
    225.     }
    227.     void csub(complex x, complex y, output complex z)
    228.     {
    229.       z.r = x.r - y.r;
    230.       z.i = x.i - y.i;
    231.     }
    233.     void csub(color x, complex y, output complex z)
    234.     {
    235.       z.r = x - y.r;
    236.       z.i = -y.i;
    237.     }
    239.     void csub(complex x, color y, output complex z)
    240.     {
    241.       z.r = x.r - y;
    242.       z.i = x.i;
    243.     }
    245.     void csqrt(complex x, output complex z)
    246.     {
    247.       color r = sqrt(sqrt(color(x.r*x.r) + color(x.i*x.i)));
    248.       color phi = atan2(x.i,x.r)/2;
    250.       z.r = r*cos(phi);
    251.       z.i = r*sin(phi);
    252.     }
    254.     void cexp(complex x, output complex z)
    255.     {
    256.         complex tmp;
    257.         tmp.r = exp(x.r)*cos(x.i);
    258.         tmp.i = exp(x.r)*sin(x.i);
    259.         z.r=tmp.r;z.i=tmp.i;
    260.     }  
    262.     void kz(complex epsilon, color k, color kx, output complex kz_)
    263.     {
    264.       complex root;
    265.       cmult(epsilon,k*k,root);          //eps*k^2
    266.       csub(root, kx*kx, root);    //eps*k^2-kx^2
    267.       csqrt(root, kz_);
    268.     }
    270.     // void kx_k(float cosi, color wavelength, color n, output color kx, output color k)
    271.     // {
    272.       // k = M_PI*2/wavelength;
    274.       // kx = sin(acos(cosi))*k*n;
    275.     // }
    277.     void rp(complex kzi, complex kzk, complex ei, complex ek, output complex _rp)
    278.     {
    279.      //return (ek kzi - ei kzk)/(ek kzi + ei kzk);
    280.      complex ek_kzi; complex ei_kzk;
    281.      cmult(ek, kzi, ek_kzi);
    282.      cmult(ei, kzk, ei_kzk);
    284.      complex N; complex D;
    285.      csub(ek_kzi, ei_kzk, N);
    286.      cadd(ek_kzi, ei_kzk, D);
    288.      cdiv(N,D,_rp);
    289.     }
    291.     void rp2tp(complex rp, complex ei, complex ek, output complex _tp)
    292.     { 
    293.      complex tmp;
    294.      cdiv(ei, ek, tmp); //tmp = ei/ek
    295.      csqrt(tmp,tmp);    //tmp = sqrt(ei/ek)
    297.      cadd(1,rp,_tp);    //tp = 1 + rp
    299.      cmult(_tp,tmp,_tp);    //tp = (1+rp)*sqrt(ei/ek)
    301.     }
    302.     void rs(complex kzi, complex kzk, output complex _rs)
    303.     {
    304.      //return (kzi - kzk)/(kzi + kzk);
    305.      complex N; complex D;
    306.      csub(kzi, kzk, N);
    307.      cadd(kzi, kzk, D);
    309.      cdiv(N,D,_rs);
    310.     }
    312. #ifdef SINGLELAYER
    313.     void FR(complex r01, complex r12, complex kz1, color d, output complex _R, output complex D)
    314.     {
    315.         // R = (r01+r12*exp(2*1i*kz1*d))/(1+r01*r12*exp(2i*kz1*d))
    316.         complex im = {color(0),color(1)};
    318.         complex exponent; 
    319.         cmult(2,im,exponent); //2i;
    320.         cmult(exponent, kz1,exponent); //2i*kz1;
    321.         cmult(d, exponent, exponent); //2i*kz1*d
    323.         cexp(exponent, exponent); //exp(2i*kz1*d)
    325.         cmult(r12, exponent, exponent); //r12*exp(2i*kz1*d)
    327.         complex Num;
    328.         cadd(r01, exponent, Num); // N = r01+r12.*exp(2*1i.*kz1*d)
    330.         cmult(r01,exponent, exponent); //exponent = r01*r12.*exp(2*1i.*kz1*d)
    332.         cadd(1, exponent,D); // D = 1+r01*r12.*exp(2*1i.*kz1*d)
    334.         cdiv(Num,D, _R);
    335.     }
    337.     void FT(complex t01, complex t12, complex kz1, color d, complex D, output complex _T)
    338.     {
    339.         complex im = {color(0),color(1)};
    340.         complex Num;
    341.         cmult(im, kz1,Num); //1i*kz1;
    342.         cmult(d, Num, Num); //1i*kz1*d
    344.         cexp(Num, Num); //exp(1i*kz1*d)
    346.         cmult(t01,Num,Num);
    347.         cmult(t12,Num,Num);
    348.         cdiv(Num,D,_T);
    349.     }
    350. #endif
    352. #if defined(SINGLELAYER) && !defined(DOUBLELAYER)
    353.     void RperpRpar(color kx, color k, complex eps_0, complex eps_1, complex eps_2, float thickness_T, output complex Rp, output complex Rs, output complex Tp, output complex Ts, int calcT)
    354.     {    
    355.         complex kz0; complex kz1; complex kz2;       
    357.         kz(eps_0, k, kx, kz0);
    358.         kz(eps_1, k, kx, kz1);
    359.         kz(eps_2, k, kx, kz2);
    361.         complex r01; complex r12; complex D; // Denominator, reuse for speed
    362.         complex t01; complex t12;
    364.         rp(kz0,kz1,eps_0,eps_1, r01); 
    365.         rp(kz1,kz2,eps_1, eps_2, r12); 
    367.         FR(r01,r12,kz1,thickness_T,Rp, D);    
    369.         if (calcT)
    370.         {
    371.                 //Transmission = t01*t12*(exp(1i*kz1*thickness_T))/(1+r01*r12*(exp(2i*kz1*thickness_T)))
    372.             rp2tp(r01, eps_0,eps_1, t01); 
    373.             rp2tp(r12, eps_1, eps_2, t12); 
    375.             FT(t01,t12,kz1,thickness_T,D, Tp);
    376.         }
    378.         // Now S-polarized light
    379.         rs(kz0,kz1, r01);
    380.         rs(kz1,kz2, r12); 
    381.         FR(r01,r12,kz1,thickness_T,Rs, D);
    383.         if (calcT)
    384.         {
    385.             cadd(1,r01,t01);
    386.             cadd(1,r12,t12);
    387.             FT(t01,t12,kz1,thickness_T,D, Ts);
    388.         } 
    390.     }
    391. #endif
    393. #ifdef DOUBLELAYER
    394. void RperpRpar(color kx, color k, complex eps_0, complex eps_1, complex eps_2, complex eps_3, float thickness_T, float thickness_M, output complex Rp, output complex Rs, output complex Tp, output complex Ts, int calcT)
    395.     {    
    396.         complex kz0; complex kz1; complex kz2; complex kz3;       
    398.         kz(eps_0, k, kx, kz0);
    399.         kz(eps_1, k, kx, kz1);
    400.         kz(eps_2, k, kx, kz2);
    401.         kz(eps_3, k, kx, kz3);
    403.         complex r01; complex r12; complex r23; complex D; // Denominator, reuse for speed
    404.         complex t01; complex t12; complex t23;
    406.         rp(kz0,kz1,eps_0, eps_1, r01); 
    407.         rp(kz1,kz2,eps_1, eps_2, r12); 
    408.         rp(kz2,kz3,eps_2, eps_3, r23);
    410.         // calculate last layer first
    411.         FR(r12,r23,kz2,thickness_M,Rp, D);    
    413.         // because we reuse D, do not continue with thin film calculation before calculating partial Tp        
    414.         if (calcT)
    415.         {
    416.                 //Transmission = t01*t12*(exp(1i*kz1*thickness_T))/(1+r01*r12*(exp(2i*kz1*thickness_T)))
    417.             rp2tp(r12, eps_1,eps_2, t12); 
    418.             rp2tp(r23, eps_2, eps_3, t23); 
    420.             FT(t12,t23,kz2,thickness_M,D, Tp);
    421.         }
    423.         FR(r01,Rp,kz1,thickness_T,Rp,D);
    425.         if (calcT)
    426.         {
    427.             //Transmission = t01*t12*(exp(1i*kz1*thickness_T))/(1+r01*r12*(exp(2i*kz1*thickness_T)))
    428.             rp2tp(r01, eps_0,eps_1, t01); 
    430.             FT(t01,Tp,kz1,thickness_T,D, Tp);
    431.         }
    433.         // Now S-polarized light
    434.         rs(kz0,kz1, r01);
    435.         rs(kz1,kz2, r12); 
    436.         rs(kz2,kz3, r23);
    437.         FR(r12,r23,kz2,thickness_M,Rs, D);
    440.         if (calcT)
    441.         {
    442.             cadd(1,r12,t12);
    443.             cadd(1,r23,t23);
    444.             FT(t12,t23,kz2,thickness_M,D, Ts);
    445.         } 
    447.         FR(r01,Rs,kz1,thickness_T,Rs, D);
    449.         if (calcT)
    450.         {
    451.             cadd(1,r01,t01);
    452.             FT(t01,Ts,kz1,thickness_T,D, Ts);
    453.         }        
    454.     }
    455. #endif
    457. #if !defined(SINGLELAYER)
    459.     void RperpRpar(color kx, color k, complex eps_0, complex eps_2, output complex Rp, output complex Rs, output complex Tp, output complex Ts, int calcT)
    460.     {    
    461.         complex kz0; complex kz2;       
    463.         kz(eps_0, k, kx, kz0);
    464.         kz(eps_2, k, kx, kz2);
    466.         rp(kz0,kz2,eps_0,eps_2, Rp);
    467.         rs(kz0,kz2, Rs);
    469.         if (calcT)
    470.         {
    471.                 rp2tp(Rp,eps_0,eps_2, Tp); 
    472.                 cadd(Rs,1,Ts);        
    473.         }
    474.     }    
    475. #endif
    477.     int calcT = isconnected(Transmission);
    479.     color k = M_PI*2/wavelength; //denominator is wavelength in um
    481.     complex Rp; complex Rs; complex Ts; complex Tp;
    482. #ifdef DOUBLELAYER
    483.     float local_d1 = thickness_T; float local_d2 = thickness_M;
    484. #endif
    485.     complex nsub = {IOR_base_material, k_base_material}; 
    488.     complex eps3; complex eps0;
    490.     if (backfacing())
    491.     {
    492.       cmult(nsub,nsub,eps0);
    493.       eps3.r=IOR_outside; eps3.i=0;
    494.     }
    495.     else
    496.     {
    497.       cmult(nsub,nsub,eps3);
    498.       eps0.r=IOR_outside; eps0.i=0;
    499.     }
    501.     // the following is nonsense for backfacing metallic substrates  
    502.     // but so is a ray of light coming out of a metallic substrate
    503.     color kx = sin(acos(dot(I,Normal)))*k*sqrt(eps0.r);
    505. #ifdef SINGLELAYER
    506.     complex nfilm = {IOR_top, k_top};
    508.     complex eps1;
    509.     cmult(nfilm,nfilm,eps1);
    510. #endif 
    512. #ifdef DOUBLELAYER
    513.     nfilm.r = IOR_middle; nfilm.i =k_middle;
    515.     complex eps2;
    516.     cmult(nfilm,nfilm,eps2);
    517.     if (backfacing())
    518.     {
    519.       complex tmp = eps1;
    520.       eps1 = eps2;
    521.       eps2 = tmp;
    522.       local_d1 = thickness_M; local_d2 = thickness_T;
    523.     }
    524. #endif 
    526. #ifdef DOUBLELAYER
    527.     RperpRpar(kx, k, eps0, eps1, eps2, eps3, local_d1, local_d2, Rp, Rs, Tp, Ts, calcT);
    528. #elif defined(SINGLELAYER)
    529.     RperpRpar(kx, k, eps0, eps1, eps3, thickness_T, Rp, Rs, Tp, Ts, calcT);
    530. #else //no layer
    531.     RperpRpar(kx, k, eps0, eps3, Rp, Rs, Tp, Ts, calcT);
    532. #endif
    534.     color out_film = (Rp.r*Rp.r+Rp.i*Rp.i+Rs.r*Rs.r+Rs.i*Rs.i)*0.5;
    536.     float r = out_film[0], g = out_film[1], b = out_film[2];
    538.     // ACES Output
    539.     if (aces == 0)
    540.     outColor = out_film;
    541.     else
    542.     {
    543.     outColor = transform(aces_tm, vector(r,g,b));
    544.     }    
    546.     if(calcT)
    547.     {
    548. #ifdef DOUBLELAYER
    549.           if (k_top == 0 && k_middle == 0) //no absorption in the film
    550. #endif
    551. #if defined(SINGLELAYER) && !defined(DOUBLELAYER)
    552.           if (k_top == 0) //no absorption in the film
    553. #endif
    555. #ifdef SINGLELAYER
    556.                 Transmission = 1-outColor;      // Quick shortcut
    557.           else
    558.           {                  // I'm going to be pragmatic here. 
    559.                           // Strictly speaking, the angle of refraction depends not only on the IOR, but also on the absorbance of the substrate (k). 
    560.                           // But, since this isn't supported in ray tracers anyway, I'm going to let it slide.
    561.                           // Moreover, this is only really an issue when k ~ 1, and then the light is absorbed within about a wavelength anyway.
    563.                 Transmission = (Tp.r*Tp.r+Tp.i*Tp.i+Ts.r*Ts.r+Ts.i*Ts.i)*0.5;
    564.                 Transmission = backfacing()? Transmission/IOR_base_material:Transmission*IOR_base_material;
    565.             }
    566. #endif
    567. #ifndef SINGLELAYER
    568.         Transmission = 1-outColor;
    569. #endif 
    570.     }
    572.     calcT = isconnected(Trans_N);
    574.     if (isconnected(outColor_N) || calcT)
    575.     {
    577. #ifdef DOUBLELAYER
    578.     RperpRpar(0, k, eps0, eps1, eps2, eps3, thickness_T, thickness_M, Rp, Rs, Tp, Ts, calcT);
    579. #elif defined(SINGLELAYER)
    580.     RperpRpar(0, k, eps0, eps1, eps3, thickness_T, Rp, Rs, Tp, Ts, calcT);
    581. #else //no layer
    582.     RperpRpar(0, k, eps0, eps3, Rp, Rs, Tp, Ts, calcT);
    583. #endif
    585.       color outColor_N = (Rp.r*Rp.r+Rp.i*Rp.i+Rs.r*Rs.r+Rs.i*Rs.i)*0.5;
    587.       if (calcT)
    588.       {
    589. #ifdef DOUBLELAYER
    590.           if (k_top == 0 && k_middle == 0) //no absorption in the film
    591. #elif defined(SINGLELAYER)
    592.           if (k_top == 0) //no absorption in the film
    593. #else
    594.           color Trans_N = 1-outColor_N;
    595. #endif
    597. #ifdef SINGLELAYER  
    598.                 color Trans_N = 1-outColor_N;      // Quick shortcut
    599.           else
    600.           {                  
    601.                 Trans_N = (Tp.r*Tp.r+Tp.i*Tp.i+Ts.r*Ts.r+Ts.i*Ts.i)*0.5;
    602.                 Trans_N = backfacing()? Trans_N/IOR_base_material:Trans_N*IOR_base_material;
    603.           }
    604. #endif
    605.       }
    606.     }
    607. }