/***********************************************/
/*           Copyright (c) 2025 Belmu          */
/*             All Rights Reserved             */
/***********************************************/

/*
    [Credits]:
        Jessie - providing the Klein-Nishina phase and phase sampling functions (https://github.com/Jessie-LC)

    [References]:
        Nishita, T. (1993). Display of the earth taking into account atmospheric scattering. http://nishitalab.org/user/nis/cdrom/sig93_nis.pdf
        Frisvad, R., J. (2011). Importance sampling the Rayleigh phase function. https://backend.orbit.dtu.dk/ws/portalfiles/portal/6314521/3D9A1d01.pdf
        THE HENYEY-GREENSTEIN PHASE FUNCTION. (2011). https://www.astro.umd.edu/~jph/HG_note.pdf

    [Notes]:
        Phase functions represent the angular distribution of scattered radiation.
*/

const float isotropicPhase = 0.25 / PI;

float rayleighPhase(float cosTheta) {
    const float normalization = 3.0 / (16.0 * PI);
    return normalization * (1.0 + pow2(cosTheta));
}

float cornetteShanksPhase(float cosTheta, float g) {
    const float normalization = 3.0 / (4.0 * TAU);
    float gg = g * g;

    float num   = (1.0 - gg) * (1.0 + pow2(cosTheta));
    float denom = (2.0 + gg) * pow(1.0 + gg - 2.0 * g * cosTheta, 1.5);
    return normalization * (num / denom);
}

float henyeyGreensteinPhase(float cosTheta, float g) {
    const float normalization = 0.25 * RCP_PI;
    float gg = g * g;

    return normalization * (1.0 - gg) * rcp(pow(1.0 + gg - 2.0 * g * cosTheta, 1.5));
}

float kleinNishinaPhase(float cosTheta, float e) {
    return e / (TAU * (e * (1.0 - cosTheta) + 1.0) * log(e * 2.0 + 1.0));
}

#if defined STAGE_FRAGMENT || defined STAGE_COMPUTE

    vec3 sampleRayleighPhase() {
        float xi = randF();

        float u        = -pow(2.0 * (xi * 2.0 - 1.0) + pow(4.0 * pow2(xi * 2.0 - 1.0) + 1.0, 0.5), 1.0 / 3.0);
        float cosTheta = u - 1.0 / u;

        return generateUnitVector(vec2(randF(), cosTheta * 0.5 + 0.5));
    }

    vec3 sampleHenyeyGreensteinPhase(float g) {
        float P  = randF();
        float s  = P * 2.0 - 1.0;
        float gg = g * g;

        float cosTheta = (0.5 / g) * (1.0 + gg - pow2((1.0 - gg) / (1.0 + g * s)));

        return generateUnitVector(vec2(randF(), cosTheta * 0.5 + 0.5));
    }

    vec3 sampleKleinNishinaPhase(float e) {
        float cosTheta = (-pow(e * 2.0 + 1.0, 1.0 - randF()) + e + 1.0) / e;
        return generateUnitVector(vec2(randF(), cosTheta * 0.5 + 0.5));
    }

#endif
