/*
====================================================================================================

    Copyright (C) 2020 RRe36

    All Rights Reserved unless otherwise explicitly stated.


    By downloading this you have agreed to the license and terms of use.
    These can be found inside the included license-file
    or here: https://rre36.com/copyright-license

    Violating these terms may be penalized with actions according to the Digital Millennium
    Copyright Act (DMCA), the Information Society Directive and/or similar laws
    depending on your country.

====================================================================================================
*/

#include "const.glsl"
#include "density.glsl"

vec2 rsi(vec3 pos, vec3 dir, float r) {
    float b     = dot(pos, dir);
    float det   = sqr(b) - dot(pos, pos) + sqr(r);

    if (det < 0.0) return vec2(-1.0);

        det     = sqrt(det);

    return vec2(-b) + vec2(-det, det);
}

vec2 airPhaseFunction(float cosTheta) {
    return vec2(rayleighPhase(cosTheta), mieCS(cosTheta, airMieG));
}

vec3 atmosphericScattering(vec3 direction, mat2x3 lightDirection, vec3 illuminance, mat2x3 illuminanceMod) {
    vec3 position   = vec3(0.0, eyeAltitude + planetRad, 0.0);

    vec2 airDist    = rsi(position, direction, atmosRad);
    vec2 planetDist = rsi(position, direction, planetRad - 5e3);

    bool isPlanet = planetDist.y >= 0.0;

    vec2 dist   = vec2(0.0);
        dist.x  = isPlanet && planetDist.x < 0.0 ? planetDist.y : max(airDist.x, 0.0);
        dist.y  = isPlanet && planetDist.x > 0.0 ? planetDist.x : airDist.y;

    float stepSize  = (dist.y - dist.x) / airScatterIterations;

    float rayLocation = dist.x;

    vec3 rayPos     = position + direction * (dist.x + stepSize * 0.5);

    vec3 airmass  = vec3(0.0);

    vec2 cosTheta = vec2(dot(direction, lightDirection[0]), dot(direction, lightDirection[1]));
    vec4 airPhase = vec4(airPhaseFunction(cosTheta.x), airPhaseFunction(cosTheta.y));
    const float phaseIso = 0.25 * pi;

    mat2x3 sunScattering    = mat2x3(0.0);
    mat2x3 moonScattering   = mat2x3(0.0);
    vec3 multiScattering    = vec3(0.0);

    vec3 transmittance      = vec3(1.0);

    vec2 illuminanceIntensity = vec2(0.0);

    for (uint i = 0; i < airScatterIterations; ++i) {
        if (airmass.y > 1e35) break;
        float elevation = length(rayPos) - planetRad;

        vec3 density    = getAirDensity(elevation) * stepSize;

        illuminanceIntensity += exp(-elevation / illuminanceFalloff) * stepSize;

            airmass    += density;

        vec3 stepOpticalDepth = airExtinctMat * density;

        vec3 stepTransmittance = saturate(expf(-stepOpticalDepth));
        vec3 stepTransmittedFraction = saturate((stepTransmittance - 1.0) / -stepOpticalDepth);
        vec3 visScattering  = transmittance * stepTransmittedFraction;

        vec3 sunAirmass = getAirmass(rayPos, lightDirection[0], airmassIterations);
        vec3 moonAirmass = getAirmass(rayPos, lightDirection[1], airmassIterations);

        vec3 sunlightAtten = exp(-airExtinctMat * sunAirmass) * visScattering;
        vec3 moonlightAtten = exp(-airExtinctMat * moonAirmass) * visScattering;

        sunScattering[0]   += sunlightAtten * density.x;
        sunScattering[1]   += sunlightAtten * density.y;

        moonScattering[0]  += moonlightAtten * density.x;
        moonScattering[1]  += moonlightAtten * density.y;

        multiScattering    += (airScatterMat * density.xy) * visScattering;

        rayPos     += direction * stepSize;

        transmittance *= stepTransmittance;
    }

    illuminanceIntensity    = max0(illuminanceIntensity / atmosDepth);

    //return vec3(illuminanceIntensity.xyy);

    sunScattering[0]   *= airPhase.x;
    sunScattering[1]   *= airPhase.y;

    moonScattering[0]  *= airPhase.z;
    moonScattering[1]  *= airPhase.w;

    vec3 sunColor   = airScatterMat[0] * sunScattering[0] + airScatterMat[1] * sunScattering[1];
        sunColor   *= sunIllum;

    vec3 moonColor  = airScatterMat[0] * moonScattering[0] + airScatterMat[1] * moonScattering[1];
        moonColor  *= moonIllum;

    vec2 illumPhase = mix(airPhase.xz, vec2(1.0), 0.5);

    vec3 multiIllum = getLuma(illuminance) * pow(normalize(illuminance), vec3(rpi)) * illuminanceIntensity.x * halfPi;
        multiIllum += illuminance.z * illuminanceIntensity.y / halfPi;
        multiIllum *= illuminanceMod[0] * illumPhase.x + illuminanceMod[1] * illumPhase.y;
        //multiIllum *= direction.x > 0.0 ? 1.0 : 0.0;

    vec3 multiColor = multiScattering * (illuminance * sqrt2 * skyIlluminanceMult + multiIllum * skyMultiscatterMult) * phaseIso;

    return sunColor + moonColor + multiColor;
}