3dgs_sh.,vex

float C0 = 0.28209479177387814;
float C1 = 0.4886025119029199;
float C2[] = {1.0925484305920792, -1.0925484305920792, 0.31539156525252005, -1.0925484305920792, 0.5462742152960396};
float C3[] = {-0.5900435899266435, 2.890611442640554, -0.4570457994644658, 0.3731763325901154, -0.4570457994644658, 1.445305721320277, -0.5900435899266435};

// Function to evaluate spherical harmonics at unit directions
function float[] eval_sh(int deg; float sh[]; vector dir; float C0; float C1; float C2[]; float C3[])
{
    float x, y, z;
    x = dir.x;
    y = dir.y;
    z = dir.z;
    
    float result[];
    float xx, yy, zz, xy, yz, xz;

    // Evaluate spherical harmonics
    for (int i = 0; i < len(sh); ++i)
    {
        append(result, C0 * sh[i]);
    }
    
    if (deg > 0)
    {
        for (int i = 0; i < len(sh); ++i)
        {
            result[i] -= C1 * y * sh[1];
        }

        for (int i = 0; i < len(sh); ++i)
        {
            result[i] += C1 * z * sh[2];
        }

        for (int i = 0; i < len(sh); ++i)
        {
            result[i] -= C1 * x * sh[3];
        }

        if (deg > 1)
        {
            xx = x * x;
            yy = y * y;
            zz = z * z;
            xy = x * y;
            yz = y * z;
            xz = x * z;

            for (int i = 0; i < len(sh); ++i)
            {
                result[i] += C2[0] * xy * sh[4];
            }

            for (int i = 0; i < len(sh); ++i)
            {
                result[i] += C2[1] * yz * sh[5];
            }

            for (int i = 0; i < len(sh); ++i)
            {
                result[i] += C2[2] * (2.0 * zz - xx - yy) * sh[6];
            }

            for (int i = 0; i < len(sh); ++i)
            {
                result[i] += C2[3] * xz * sh[7];
            }

            for (int i = 0; i < len(sh); ++i)
            {
                result[i] += C2[4] * (xx - yy) * sh[8];
            }

            if (deg > 2)
            {
                for (int i = 0; i < len(sh); ++i)
                {
                    result[i] += C3[0] * y * (3 * xx - yy) * sh[9];
                }

                for (int i = 0; i < len(sh); ++i)
                {
                    result[i] += C3[1] * xy * z * sh[10];
                }

                for (int i = 0; i < len(sh); ++i)
                {
                    result[i] += C3[2] * y * (4 * zz - xx - yy) * sh[11];
                }

                for (int i = 0; i < len(sh); ++i)
                {
                    result[i] += C3[3] * z * (2 * zz - 3 * xx - 3 * yy) * sh[12];
                }

                for (int i = 0; i < len(sh); ++i)
                {
                    result[i] += C3[4] * x * (4 * zz - xx - yy) * sh[13];
                }

                for (int i = 0; i < len(sh); ++i)
                {
                    result[i] += C3[5] * z * (xx - yy) * sh[14];
                }

                for (int i = 0; i < len(sh); ++i)
                {
                    result[i] += C3[6] * x * (xx - 3 * yy) * sh[15];
                }
            }
        }
    }

    return result;
}

int active_sh_degree = 3;

// Lookup diction.
vector dir = normalize(v@P - v@P + v@N);
v@dir = dir;

// Evaluate SH and extract color.
float sh_values[] = eval_sh(active_sh_degree, f[]@sh, dir, C0, C1, C2, C3);
float r_component = C1 * sh_values[0]
                  + C1 * sh_values[1]
                  + C1 * sh_values[2]
                  + C2[0] * sh_values[3]
                  + C2[1] * sh_values[4]
                  + C2[2] * sh_values[5]
                  + C2[3] * sh_values[6]
                  + C2[4] * sh_values[7]
                  + C3[0] * sh_values[8]
                  + C3[1] * sh_values[9]
                  + C3[2] * sh_values[10]
                  + C3[3] * sh_values[11]
                  + C3[4] * sh_values[12]
                  + C3[5] * sh_values[13]
                  + C3[6] * sh_values[14];

float g_component = C1 * sh_values[15]
                  + C1 * sh_values[16]
                  + C1 * sh_values[17]
                  + C2[0] * sh_values[18]
                  + C2[1] * sh_values[19]
                  + C2[2] * sh_values[20]
                  + C2[3] * sh_values[21]
                  + C2[4] * sh_values[22]
                  + C3[0] * sh_values[23]
                  + C3[1] * sh_values[24]
                  + C3[2] * sh_values[25]
                  + C3[3] * sh_values[26]
                  + C3[4] * sh_values[27]
                  + C3[5] * sh_values[28]
                  + C3[6] * sh_values[29];

float b_component = C1 * sh_values[30]
                  + C1 * sh_values[31]
                  + C1 * sh_values[32]
                  + C2[0] * sh_values[33]
                  + C2[1] * sh_values[34]
                  + C2[2] * sh_values[35]
                  + C2[3] * sh_values[36]
                  + C2[4] * sh_values[37]
                  + C3[0] * sh_values[38]
                  + C3[1] * sh_values[39]
                  + C3[2] * sh_values[40]
                  + C3[3] * sh_values[41]
                  + C3[4] * sh_values[42]
                  + C3[5] * sh_values[43]
                  + C3[6] * sh_values[44];

// Adjust with base SH coefficient and add director color.                  
vector rgb_result = set(r_component, g_component, b_component);
v@Cd = set(rgb_result.x, rgb_result.y, rgb_result.z);
//v@Cd *= C0;
//v@Cd += .5;