Phong + Freshnel

// PHONG EXPLICIT DIFF OR SPEC SHLICK FRESNEL
    else if (hit_material->shader == 9) {

        float gamma = hit_material->shininess;
        Vector3 omega_i, omega_r;

        if(hit_material->diffuse.ToColor3f().AsLinear() + hit_material->specular.ToColor3f().AsLinear() > 1) return Color3f{0, 0, 0};

        // BEGIN IMPORTANCE SAMPLING
        bool sample_diffuse;
        float random0 = rng.uniformThread() * (hit_material->diffuse.MaxComponent() + hit_material->specular.MaxComponent());
        if( random0 < hit_material->diffuse.MaxComponent())
            sample_diffuse = true;
        else
            sample_diffuse = false;
        // END IMPORTANCE SAMPLING

        if(sample_diffuse) {
            omega_i = CosineWeightedSampling(normal, pdf);

            RTCRay rayIncident = PrepareRay(hit_position, omega_i);
            Color3f L_i = TraceRay(rayIncident, reflectCount + 1, n1, isFromAir);

            float cosTheta_o = omega_o.DotProduct(normal);
            float cosTheta_i = omega_i.DotProduct(normal);
            float cosTheta_r = omega_i.DotProduct(omega_r);
            float raisedCosTheta_r = pow(cosTheta_r, gamma);

            // SPECULAR REFLECTANCE
            Color3f ps = hit_material->specular.ToColor3f().AsLinear();
            Color3f F = Fresnel(ps, cosTheta_o);

            // DIFFUSE REFLECTANCE
            Color3f pd = hit_material->diffuse.ToColor3f().AsLinear();
            Color3f R_d = pd * (1.0f - F.MaxComponent()) / (1.0f - ps.MaxComponent());
            Color3f kd = R_d * M_1_PI;

            Color3f f_r = kd;
            Color3f L_r_indir = L_i * f_r * cosTheta_i / (pdf * alpha);

            if (isnan(L_r_indir.r) || isnan(L_r_indir.g) || isnan(L_r_indir.b)){
                printf("OPSIE: L_i * f_r * %.3f / (%.3f * %.3f)\n", cosTheta_i, pdf, alpha);
                printf("L_i: (%.3f, %.3f, %.3f)\n", L_i.r, L_i.g, L_i.b);
                printf("f_r: (%.3f, %.3f, %.3f)\n", f_r.r, f_r.g, f_r.b);
                return Color3f{ 0.0, 0.0, 0.0};
            }

            return L_r_indir + L_e * f_r;
        } else {
            omega_r = ReflectVector(omega_o, normal);
            omega_i = CosineLobeSampling(omega_r, pdf, gamma);
            RTCRay rayIncident = PrepareRay(hit_position, omega_i);
            Color3f L_i = TraceRay(rayIncident, reflectCount + 1, n1, isFromAir);

            float cosTheta_o = omega_o.DotProduct(normal);
            float cosTheta_i = omega_i.DotProduct(normal);
            float cosTheta_r = omega_i.DotProduct(omega_r);
            float raisedCosTheta_r = pow(cosTheta_r, gamma);

            // ENERGY CONSERVING
            float I_M_1 = (gamma + 2.0f) * 0.5f * M_1_PI;
            float I_M = (gamma + 2.0f) / 2.0f * M_PI;
            //float I_M_1_calc = 1 / calc_I_M(normal.DotProduct(-omega_o), gamma);

            // SPECULAR REFLECTANCE
            Color3f ps = hit_material->specular.ToColor3f().AsLinear();
            Color3f F = Fresnel(ps, cosTheta_o);
            Color3f ks = F * raisedCosTheta_r * I_M_1 / cosTheta_i;

            Color3f f_r = ks;
            Color3f L_r_indir = L_i * f_r * cosTheta_i / (pdf * alpha);

            if (isnan(L_r_indir.r) || isnan(L_r_indir.g) || isnan(L_r_indir.b)){
                printf("OPSIE: L_i * f_r * %.3f / (%.3f * %.3f)\n", cosTheta_i, pdf, alpha);
                printf("L_i: (%.3f, %.3f, %.3f)\n", L_i.r, L_i.g, L_i.b);
                printf("f_r: (%.3f, %.3f, %.3f)\n", f_r.r, f_r.g, f_r.b);
                return Color3f{ 0.0, 0.0, 0.0};
            }

            return L_r_indir + L_e * f_r;
        }



    }