

mat2x3 atmosphericScattering(vec3 rayPos, vec3 rayDir, vec3 sunDir, vec3 moonDir, vec3 sunIrradiance, vec3 moonIrradiance, vec3 ambient, const int iSteps) {
    const float rSteps = rcp(float(iSteps));

    vec2 airDist = rsi(rayPos, rayDir, planetRadius + atmosHeight);
    vec2 planetDist = rsi(rayPos, rayDir, atmosLowerLimit);
         airDist.y = planetDist.x > 0.0 ? min(airDist.y, planetDist.x) : airDist.y;
	float dd = (airDist.y - max0(airDist.x)) * rSteps;
    vec3 rayStep = rayDir*dd;
    rayPos = rayPos + rayStep*(0.33);

    vec2 cosTheta = vec2(dot(rayDir, sunDir), dot(rayDir, moonDir));
    vec4 phase = vec4(rayleighPhase(cosTheta.x), miePhase(cosTheta.x, mieg),
                      rayleighPhase(cosTheta.y), miePhase(cosTheta.y, mieg));

    vec3 transmittance = vec3(1.0);
    vec3 sunInScatterRayl = vec3(0.0);
    vec3 sunInScatterMie = vec3(0.0);
    vec3 moonInScatterRayl = vec3(0.0);
    vec3 moonInScatterMie = vec3(0.0);
    vec3 multipleScattering = vec3(0.0);

    float distanceTraveled = 0.0;

    for(int i = 0; i < iSteps; i++) {
        float h = getElevation(rayPos);
        vec3 density = getAtmosDensity(h);
        vec3 stepOpticalDepth = atmosExtinctionMat * (density * dd);
        vec3 stepTransmittance = clamp01(expf(-stepOpticalDepth));

        vec3 stepIntegral = clamp01((1.0 - stepTransmittance)/stepOpticalDepth);
        vec3 visibleScattering = transmittance * stepIntegral;

        vec3 sunOpticalDepth = getAtmosOpticalDepth(dot(normalize(rayPos), sunDir), h);//getAtmosOpticalDepth(rayPos, sunDir, jSteps);
        vec3 moonOpticalDepth = getAtmosOpticalDepth(dot(normalize(rayPos), moonDir), h);//getAtmosOpticalDepth(rayPos, moonDir, jSteps);

        vec3 sunTransmittance = exp(-sunOpticalDepth) * visibleScattering;
        vec3 moonTransmittance = exp(-moonOpticalDepth) * visibleScattering;

        sunInScatterRayl   += sunTransmittance * density.x * dd;
        sunInScatterMie   += sunTransmittance * density.y * dd;

        moonInScatterRayl  += moonTransmittance * density.x * dd;
        moonInScatterMie  += moonTransmittance * density.y * dd;

        vec3 stepScattering = atmosScatterMat * density.xy * dd;
        vec3 stepScatteringAlbedo = stepScattering / stepOpticalDepth;
        vec3 multipleScatteringFactor = 0.84 * stepScatteringAlbedo;
        vec3 multipleScatteringEnergy = multipleScatteringFactor / (1.0 - multipleScatteringFactor);
        multipleScattering += multipleScatteringEnergy * visibleScattering * stepScattering;
        // multipleScattering += (atmosScatterMat*density.xy) * visibleScattering * dd;


        transmittance *= stepTransmittance;

        rayPos += rayStep;
    }

    float k = 2.0;
    sunInScatterRayl   = k*(atmosScatterMat[0] * (phase.x*sunInScatterRayl)) * sunIrradiance;
    sunInScatterMie    = k*(atmosScatterMat[1] * (phase.y*sunInScatterMie)) * sunIrradiance;
    moonInScatterRayl  = k*(atmosScatterMat[0] * (phase.z*moonInScatterRayl)) * moonIrradiance;
    moonInScatterMie   = k*(atmosScatterMat[1] * (phase.w*moonInScatterMie)) * moonIrradiance;
    multipleScattering = multipleScattering * ambient * 0.5;

    

    mat2x3 scatter_trans = mat2x3(max(sunInScatterRayl + sunInScatterMie + moonInScatterRayl + moonInScatterMie + multipleScattering, 0.0), transmittance);
    return scatter_trans;
}

// mat2x3 atmosphericScatteringAerialPerspective(vec3 rayPos, vec3 rayDir, float dist, vec3 sunDir, vec3 sunIrradiance, vec3 ambient, const int iSteps, const int jSteps) {
//     const float rSteps = rcp(float(iSteps));

//     vec2 planetDist = rsi(rayPos, rayDir, planetRadius-1e3);
//          dist = planetDist.x > 0.0 ? min(dist, planetDist.x) : dist;
// 	float dd = (dist) * rSteps;
//     vec3 rayStep = rayDir*dd;
//     rayPos = rayPos + rayDir*(dd*0.5);

//     float cosTheta = dot(rayDir, sunDir);
//     vec2 phase = vec2(rayleighPhase(cosTheta.x), miePhase(cosTheta.x, mieg));

//     vec3 transmittance = vec3(1.0);
//     vec3 sunInScatterRayl = vec3(0.0);
//     vec3 sunInScatterMie = vec3(0.0);
//     vec3 multipleScattering = vec3(0.0);

//     float distanceTraveled = 0.0;

//     for(int i = 0; i < iSteps; i++) {
//         vec3 density = getAtmosDensity(getElevation(rayPos));
//         vec3 stepOpticalDepth = atmosExtinctionMat * (density * dd);
//         vec3 stepTransmittance = clamp01(expf(-stepOpticalDepth));

//         vec3 stepIntegral = clamp01((1.0 - stepTransmittance)/stepOpticalDepth);
//         vec3 visibleScattering = transmittance * stepIntegral;

//         vec3 sunOpticalDepth = getAtmosOpticalDepth(rayPos, sunDir, jSteps);

//         vec3 sunTransmittance = exp(-sunOpticalDepth) * visibleScattering;

//         sunInScatterRayl   += sunTransmittance * density.x * dd;
//         sunInScatterMie   += sunTransmittance * density.y * dd;

//         vec3 stepScattering = atmosScatterMat * density.xy * dd;
//         vec3 stepScatteringAlbedo = stepScattering / stepOpticalDepth;
//         vec3 multipleScatteringFactor = 0.84 * stepScatteringAlbedo;
//         vec3 multipleScatteringEnergy = multipleScatteringFactor / (1.0 - multipleScatteringFactor);
//         multipleScattering += multipleScatteringEnergy * visibleScattering * stepScattering;
//         // multipleScattering += (atmosScatterMat*density.xy) * visibleScattering * dd;


//         transmittance *= stepTransmittance;

//         rayPos += rayStep;
//     }

//     sunInScatterRayl   = (atmosScatterMat[0] * (phase.x*sunInScatterRayl)) * sunIrradiance;
//     sunInScatterMie    = (atmosScatterMat[1] * (phase.y*sunInScatterMie)) * sunIrradiance;
//     multipleScattering = multipleScattering * ambient * 0.5;

    

//     mat2x3 scatter_trans = mat2x3(max(sunInScatterRayl + sunInScatterMie + multipleScattering, 0.0), transmittance);
//     return scatter_trans;
// }
