atomic_orbital_shader.osl

#include "stdosl.h"

struct complex
{
    float real;
    float imag;
};

color phase2rgb(float phase,float offset)
{
    float pi = M_PI;
    return color("hsv",mod(phase/(2*pi)+offset,1.),1.0,1.0);
}

int fac(int arg)
{
    int tmp = arg;
    int fac = 1;
    while (tmp>0)
    {
        fac*=tmp;
        tmp-=1;
    }
    return fac;
}

float laguerrePol(int alpha,int n,float x)
{
    float tmp=0.;
    for (int i=0; i<=n; i++)
    {
        tmp+=fac(n+alpha)*pow(-1,i)/fac(alpha+i)/fac(n-i)*pow(x,i)/fac(i);
    }
    return tmp;
}





complex sphericalHarm_(int m, int l, float sin_theta, float cos_theta, float phi)
{
    float pi = M_PI;
    
    float tmp=0.;
    for (int q=0; q<=(l-abs(m))/2.; q++)
    {
        tmp+=pow(-1,l+q) * fac(l)/(fac(l-q)*fac(q)) * fac(2*(l-q))/((fac(2*(l-q) - (l+abs(m)))*fac(l+abs(m)))) * pow(cos_theta, l-abs(m)-2*q);
    }
    tmp*=pow(2,-l)*pow(sin_theta,abs(m));
    tmp*=sqrt((2*l+1)*0.5 + fac(l-abs(m))/fac(l) + fac(l+abs(m))/fac(l) );
    tmp/=sqrt(2*pi);
    complex z;
    z.real = tmp*cos(m*phi);
    z.imag = tmp*sin(m*phi);
    return z;
}



complex Psi_(int n, int l, int m, int Z,float r, float phi, float sin_theta, float cos_theta, float a)
{
    float pi = M_PI;
    float k_ = Z/(n*a);
    float R = sqrt(pow(2*k_,3)*fac(n-l-1)/(2*n*fac(n+l)))*exp(-r*k_)*pow(2*r*k_,l)*laguerrePol(2*l+1,n-l-1,2*r*k_);

    complex tmp = sphericalHarm_(m,l,sin_theta,cos_theta,phi);
    tmp.real*=R;
    tmp.imag*=R;
    return tmp;
}

float Phi(float X, float Y)
{
    float pi = M_PI;
    return atan(Y/X)*(X > 0) + (atan(Y/X) + pi)*(X < 0);
}

shader generated_volume(
    vector origin = vector(0.0,0.0,0.0),
    int n = 1,
    int l = 0,
    int m = 0,
    int Z = 1,
    float a = 0.02,
    float density = 1.0,
    float offset = 15,

    //output closure color BSDF = diffuse(N),
    output closure color PA = diffuse(N),
    output closure color PD = diffuse(N),

    )
{
    //float pi = M_PI;
    
    float x = P[0]-origin[0];
    float y = P[1]-origin[1];
    float z = P[2]-origin[2];

    
    float r2 = pow(x,2)+pow(y,2)+pow(z,2);
    float r = pow(r2,0.5);
    
    float sin_theta = pow((pow(x,2)+pow(y,2))/r2,0.5);
    float cos_theta = pow(pow(z,2)/r2,0.5);

    float phi = Phi(x,y);
    
    
    complex psi_ = Psi_(n,l,m,Z,r,phi,sin_theta,cos_theta,a);
    color rgb = phase2rgb(Phi(psi_.real, psi_.imag),offset/100);
    color rgb_inverted = color("rgb",1-rgb[0], 1-rgb[1], 1-rgb[2]);
        
    PA = (rgb_inverted*absorption()+rgb*emission())*pow((pow(density*psi_.real,2)+pow(density*psi_.imag,2)),0.5);

    psi_.real = pow((pow(density*psi_.real,2)+pow(density*psi_.imag,2)),1);
    psi_.imag = 0;
    

    color rgb_ = phase2rgb(Phi(psi_.real, psi_.imag),offset/100);
    color rgb_inverted_ = color("rgb",1-rgb_[0], 1-rgb_[1], 1-rgb_[2]);

    
    PD = (rgb_inverted_*absorption()+rgb_*emission())*pow((pow(density*psi_.real,2)+pow(density*psi_.imag,2)),1);

    
    
}