PBR metal shader with accurate Fresnel based on complex refractive index (script)

UPDATE: merged with thin film interference, please see this thread: https://blenderartists.org/forum/showthread.php?403299-UPDATE-Cycles-PBR-thin-film-interference-and-metals

This is the result of messing around with complex numbers within OSL. Blend file with different metals / transition metals / semiconductors here: https://dl.dropboxusercontent.com/u/5326890/interference.blend

Examples:



/* OSL shader to calculate thin-film interference
   Use at your own risk
   v1.4, 2016 prutser. License for the *code*: CC BY-NC-SA 4.0
   Note: the output of the code may explicitly be used for commercial renders.
   
   
   v1.4 take TIR correctly into account by reverting to complex math, making the shader suitable 
        for glass->vacuum reflections. 
        Promote IORs to complex and color-dependent. 
        Handle physically-correct thin layers on metals
   v1.3: improved correctness by taking air to substrate or substrate to air directions into account
   v1.2: realized I could immensly simplify the implementation, ditching the complex math
   v1.1: fixed bug in complex math functions (pass by reference issues)*/

#include "stdosl.h"

struct complex {
        color r;
        color i;
  };
  

    
shader thin_film(
    float d = 1,
    color IORfilm = 1.6,
    color kfilm = 0,
    color IORsubstrate = 1.45,
    color ksubstrate = 0,
    normal Normal = N,
    output color F = color(0),
    output color F0 = color(0))
{
       
    void cmult(complex x, complex y, output complex z)
    {
      complex tmp; //pass by reference issues
      tmp.r = x.r*y.r-x.i*y.i;
      tmp.i = x.r*y.i + x.i*y.r;
      z.r = tmp.r; z.i = tmp.i;
    }
    
    void cmult(color x, complex y, output complex z)
    {
      z.r = x*y.r;
      z.i = x*y.i;
    }
    
    void cmult(complex x, color y, output complex z)
    {
      z.r = x.r*y;
      z.i = x.i*y;
    }

    void cdiv(complex x, complex y, output complex z)
    {
      complex ystar;
      ystar.r = y.r;
      ystar.i = -y.i;
      
      cmult(x,ystar,z); 
      z.r = z.r/color((ystar.r*ystar.r+ystar.i*ystar.i));
      z.i = z.i/color((ystar.r*ystar.r+ystar.i*ystar.i));
    }
    
    void cadd(complex x, complex y, output complex z)
    {
      z.r = x.r + y.r;
      z.i = x.i + y.i;
    }
    
    void cadd(complex x, color y, output complex z)
    {
      z.r = x.r + y;
      z.i = x.i;
    }
    
    void cadd(color x, complex y, output complex z)
    {
      z.r = x + y.r;
      z.i = y.i;
    }
    
    void csub(complex x, complex y, output complex z)
    {
      z.r = x.r - y.r;
      z.i = x.i - y.i;
    }

    void csub(color x, complex y, output complex z)
    {
      z.r = x - y.r;
      z.i = -y.i;
    }
    
    void csub(complex x, color y, output complex z)
    {
      z.r = x.r - y;
      z.i = x.i;
    }
    
    void csqrt(complex x, output complex z)
    {
      color r = sqrt(sqrt(color(x.r*x.r) + color(x.i*x.i)));
      color phi = atan2(x.i,x.r)/2;
      
      z.r = r*cos(phi);
      z.i = r*sin(phi);
    }
    
    void cexp(complex x, output complex z)
    {
        complex tmp;
        tmp.r = exp(x.r)*cos(x.i);
        tmp.i = exp(x.r)*sin(x.i);
        z.r=tmp.r;z.i=tmp.i;
    }  
    
    void kz(complex epsilon, color k, color kx, output complex kz_)
    {
      complex root;
      cmult(epsilon,k*k,root);          //eps*k^2
      csub(root, kx*kx, root);    //eps*k^2-kx^2
      csqrt(root, kz_);
    }
    
    void kx_k(float cosi, color wavelength, color n, output color kx, output color k)
    {
      k = M_PI*2/wavelength;
      
      kx = sin(acos(cosi))*k*n;
    }
    
    void rp(complex kzi, complex kzk, complex ei, complex ek, output complex _rp)
    {
     //return (-kzi/ei + kzk/ek)/(kzi/ei + kzk/ek);
     complex ki_ei; complex kk_ek;
     cdiv(kzi, ei, ki_ei);
     
     cdiv(kzk, ek, kk_ek);
     
     complex N; complex D;
     csub(kk_ek, ki_ei, N);
     cadd(ki_ei, kk_ek, D);
     
     cdiv(N,D,_rp);
    }
    
    void rs(complex kzi, complex kzk, output complex _rs)
    {
     //return (kzi - kzk)/(kzi + kzk);
     complex N; complex D;
     csub(kzi, kzk, N);
     cadd(kzi, kzk, D);
     
     cdiv(N,D,_rs);
    }
      
              
    void RperpRpar(float cosi, complex eps_0, complex eps_1, complex eps_2, color wavelength, float d, output color Rpar, output color Rperp)
    {    
        complex im;
        im.r = 0; im.i = 1;
        
        color kx; color omega;color k;
        
        kx_k(cosi, wavelength, sqrt(eps_0.r), kx, k);
             
        complex kz0; complex kz1; complex kz2;       
        
        kz(eps_0, k, kx, kz0);
        kz(eps_1, k, kx, kz1);
        kz(eps_2, k, kx, kz2);
   
                
        complex r01; complex r12;
        if (d != 0 && (eps_0.r != eps_1.r || eps_0.i != eps_1.i) && (eps_1.r != eps_2.r || eps_1.i != eps_2.i)) //thin film
        {
      
            rp(kz0,kz1,eps_0,eps_1, r01); 
            rp(kz1,kz2,eps_1, eps_2, r12); 
            // R = (r01+r12.*exp(2*1i.*kz1*d))./(1+r01.*r12.*exp(2*1i.*kz1.*d));
            
            complex exponent; complex ex;
            cmult(2,im,exponent); //2i;
            cmult(exponent, kz1,exponent); //2i*kz1;
            cmult(color(d), exponent, exponent); //2i*kz1*d

            cexp(exponent, exponent); //exp(2i*kz1*d)
                    
            cmult(r12, exponent, exponent); //r12*exp(2i*kz1*d)
            
            complex Num;
            cadd(r01,exponent, Num); // N = r01+r12.*exp(2*1i.*kz1*d)
            
            cmult(r01,exponent, exponent); //exponent = r01*r12.*exp(2*1i.*kz1*d)
            
            
            complex D;
            cadd(1, exponent,D); // D = 1+r01*r12.*exp(2*1i.*kz1*d)
            
            complex Rp;
            cdiv(Num,D, Rp);
            
            Rpar = (Rp.r*Rp.r+Rp.i*Rp.i);
           
            
            // Now S-polarized light
            rs(kz0,kz1, r01);
            rs(kz1,kz2, r12); 

            // R = (r01+r12.*exp(2*1i.*kz1*d))./(1+r01.*r12.*exp(2*1i.*kz1.*d));
            
            cmult(2,im,exponent); //2i;
            cmult(exponent, kz1,exponent); //2i*kz1;
            cmult(color(d), exponent, exponent); //2i*kz1*d
            cexp(exponent, exponent);
            cmult(r12, exponent, exponent);
            
            cadd(r01,exponent, Num); //r01+r12.*exp(2*1i.*kz1*d)
            
            cmult(r01,exponent, exponent); //r01*r12.*exp(2*1i.*kz1*d)
            
            
            cadd(1, exponent,D); //1+r01*r12.*exp(2*1i.*kz1*d)
            
            
            cdiv(Num,D, Rp);
            
            Rperp = color(Rp.r*Rp.r+Rp.i*Rp.i); 
        }
        else //no thin film, speed up calculation
        {
            
            rp(kz0,kz2,eps_0,eps_2, r01);
            Rpar = r01.r*r01.r + r01.i*r01.i;
            rs(kz0,kz2, r01);
            Rperp = r01.r*r01.r + r01.i*r01.i;
        }
        

    }       
    
    float cosi = dot(I,Normal);
    color wavelength = 1e-3*color(635,532,470); //um
    
    color Rperp; color Rpar;
    complex nsub = {IORsubstrate, ksubstrate}; 
    complex nfilm = {IORfilm, kfilm};

    complex eps1;
    cmult(nfilm,nfilm,eps1);
    complex eps2; complex eps0;
        
    if (backfacing())
    {
      cmult(nsub,nsub,eps0);
      eps2.r=1;eps2.i=0;
    }
    else
    {
      cmult(nsub,nsub,eps2);
      eps0.r=1;eps0.i=0;
    }
    
    RperpRpar(cosi, eps0, eps1, eps2, wavelength, d, Rpar, Rperp)  ;
    
    
    F = (Rpar+Rperp)*0.5;
 
    if (isconnected(F0))
    {
      RperpRpar(1, eps0, eps1, eps2, wavelength, d, Rpar, Rperp)  ;
      F0 = (Rpar+Rperp)*0.5;
    }
    
}

Cheers!

thanks vry much,
in the blender there’s no script node setup?

Yes there should be, did you look in the group node?

ahh. got it now, thanks