const vec3 BetaR = vec3(5.802e-6, 13.558e-6, 33.1e-6);
const float BetaM = 3.996e-6;
const float BetaMA = 4.40e-6;

const float Anisotropy = 0.8;

// Atmosphere height
const float Hr = 7994;
const float Hm = 1200;

const float EarthRad = 6360e3;
const float AtmRad = 6420e3;

const float TEXEL_SIZE_INV = 1.0 / 64.0; // Used to prevent interpolation between rayleigh and mie

const float ISOTROPIC_PHASE = 1 / (4 * PI);

float rayleigh_phase(float Mu) {
    return 8.0/10.0 * (7.0/5.0 + Mu/2);
}

float cs_phase(float Mu, float g) {
    float g2 = g*g;
    float A = 3 * (1 - g2) * (1 + Mu * Mu);
    float B = 8 * PI * (2 + g2) * pow(1 + g2 - 2 * g * Mu, 1.5);
    return A / B;
}

float cs_phase(float Mu) {
    return cs_phase(Mu, Anisotropy);
}

float get_delta(vec3 LightPosN) {
    return dot(gbufferModelView[1].xyz, LightPosN);
}

float get_phi(vec3 ViewPosN) {
    return dot(gbufferModelView[1].xyz, ViewPosN);
}

// takes atm_height, VdotUp, LdotUp
vec3 to_uvw(float h, float p, float d) {
    const float EarthRad2 = EarthRad*EarthRad;
    float u = sqrt((h*h - EarthRad2) / (AtmRad*AtmRad - EarthRad2));
    float v = (1 + p) / 2;
    float w = (1-exp(-2.8 * d - 0.8)) / (1 - exp(-3.6));
    return vec3(u, v, w);
}


vec3 sample_scattering(vec3 UVW, const bool IsMoon) {
    vec3 SkyColor = texture(scatteringTexture, UVW).rgb; // Divide by 8 to reduce color banding. There's probably a better solution for this
    if(!IsMoon)
        return SkyColor * DayAmbientColor;
    else
        return SkyColor * NightAmbientColor;
}

vec3 get_direct_color(const bool IsMoon) {
    const int STEP_COUNT = 8;
    const vec3 ORIGIN = vec3(0, 300 + EarthRad, 0);
    vec3 LightPosN = IsMoon ? -sunPosN : sunPosN;

    float t0Light, t1Light;
    raySphereIntersect(ORIGIN, view_player(LightPosN), AtmRad, t0Light, t1Light);
    float Step = t1Light/STEP_COUNT;
    float tCurrent = 0.0, RlhT = 0.0, MieT = 0.0;
    for(int i = 0; i < STEP_COUNT; i++) {
        vec3 P = ORIGIN + (tCurrent + Step * 0.5) * view_player(LightPosN);
        float Len = length(P) - EarthRad;
        RlhT += exp(-Len / Hr) * Step;
        MieT += exp(-Len / Hm) * Step;
        tCurrent += Step;
    }

    vec3 SkyColor = exp(-BetaR*RlhT - (BetaM+BetaMA)*MieT);

    if(!IsMoon)
        return SkyColor*SunColor;
    else
        return SkyColor*MoonColor;
}

vec3 get_shadowlight_color() {
    vec3 LightColorDirect;
    if(sunAngle < 0.5) 
        LightColorDirect = dataBuf.SunColor;
    else
        LightColorDirect = dataBuf.MoonColor;

	float LHeight = sin(sunAngle * TAU); // to_player_pos(sunPosN).y;
    LightColorDirect *= smoothstep(0, 0.05, abs(LHeight));
    return LightColorDirect;
}

vec3 get_ambient_color() {
    const float STEP_COUNT = 6.0;
    vec3 UVW = to_uvw(cameraPosition.y * 5 + EarthRad, 0.0, get_delta(sunPosN));
    UVW.z = clamp(UVW.z, TEXEL_SIZE_INV, 1 - TEXEL_SIZE_INV);
    UVW.z *= 0.5;
    vec3 Ambient = vec3(0);
    for(float i = 0.0; i < STEP_COUNT; i++) {
        UVW.y = i / STEP_COUNT / 2 + 0.5;
        Ambient += sample_scattering(UVW, false);
        Ambient += sample_scattering(vec3(UVW.xy, UVW.z + 0.5), false);
        float WMoon = clamp(to_uvw(0, 0, get_delta(-sunPosN)).z*0.5, TEXEL_SIZE_INV, 1 - TEXEL_SIZE_INV);
        Ambient += sample_scattering(vec3(UVW.xy, WMoon), true);
        Ambient += sample_scattering(vec3(UVW.xy, WMoon + 0.5), true);
    }
    return Ambient / STEP_COUNT;
}
