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

/*
    [Credits]:
        Wavelength-dependent 5th order polynomial limb darkening equation taken from NASA's SolarSoft code (https://hesperia.gsfc.nasa.gov/ssw/gen/idl/solar/darklimb_correct.pro)
        http://www.physics.hmc.edu/faculty/esin/a101/limbdarkening.pdf
*/

float physicalSun(vec3 viewDirection, vec3 sunDirection, float wavelength) {
    if (dot(viewDirection, sunDirection) < cos(sunAngularRadius)) return 0.0;
    float cosTheta   = dot(viewDirection, sunDirection);
    float distToEdge = sin(acos(cosTheta)) / sin(sunAngularRadius);

    wavelength *= 10.0; // To Angströms

    const float[6] polynomialCoeffsPS = float[](-8.9829751, 0.0069093916, -1.8144591e-6, 2.2540875e-10, -1.3389747e-14, 3.0453572e-19);
    const float[6] polynomialCoeffsNL = float[](9.2891180, -0.0062212632, 1.5788029e-6, -1.9359644e-10, 1.1444469e-14, -2.599494e-19);

    float PS, NL;
    for (int k = 0; k < 6; k++) {
        float power = pow(wavelength, k);
        PS += polynomialCoeffsPS[k] * power;
        NL += polynomialCoeffsNL[k] * power;
    }
    return 1.0 - PS - NL + PS * cos(asin(distToEdge)) + NL * pow2(cos(asin(distToEdge)));
}

#include "/include/fragment/hapke.glsl"

uniform sampler2D moon;

float getMoonAlbedo(vec3 normal, float wavelength) {
    vec2 coords = projectSphereToUV(normal);
    vec3 albedo = texture(moon, coords).rgb;
            albedo = srgbToLinear(albedo);

    return srgbToSpectrum(albedo, wavelength);
}

float physicalMoon(vec3 viewDirection, vec3 sunDirection, vec3 moonDirection, float solarIrradiance, float wavelength) {
    vec2 sphere = intersectSphere(-moonDirection, viewDirection, moonAngularRadius);

    if (sphere.y >= 0.0) {
        vec3 normal = normalize(viewDirection * sphere.x - moonDirection);
        float NdotL = saturate(dot(normal, sunDirection));

        float albedo = getMoonAlbedo(normal, wavelength);

        float brdf = Hapke_2012(normal, -viewDirection, sunDirection, albedo) * NdotL;

        if (isnan(brdf)) return 0.0;
        if (isinf(brdf)) return 1.0;

        return brdf * solarIrradiance;
    }
    return 0.0;
}
