
float planarCloudNoiseTex(sampler2D tex, vec3 pos) {
    float p = floor(pos.z);
    float f = pos.z - p;
    
    vec2 invNoisTexRes = 1.0/textureSize(tex, 0);
    // float zStretch = 10.0 * invNoisTexRes;
    
    vec2 sampleCoord = pos.xy*invNoisTexRes + p * 17.0 * invNoisTexRes;
    
    vec2 noise = texture(tex, sampleCoord).xy;//vec2(texture(tex, sampleCoord)[channel], texture(tex, sampleCoord + 17.0 * invNoisTexRes)[channel]);
    
    return mix(noise.x, noise.y, f);
}



float cloudNoiseTex(sampler2D tex, vec3 pos) {
    vec3 p = floor(pos);
    vec3 f = fract(pos);
    
    vec2 invNoisTexRes = 1.0/textureSize(tex, 0);
    // float zStretch = 10.0 * invNoisTexRes;
    
    vec2 sampleCoord = (p.xy-vec2(97.0)*p.z) + f.xy;//pos.xy*invNoisTexRes + p * 17.0 * invNoisTexRes;
    
    vec2 noise = texture(tex, sampleCoord*invNoisTexRes).xy;//vec2(texture(tex, sampleCoord)[channel], texture(tex, sampleCoord + 17.0 * invNoisTexRes)[channel]);
    
    return cube_smooth(mix(noise.x, noise.y, f.z));
}


// float cloudNoiseTex(sampler2D tex, vec3 pos, int channel) {
//     vec3 p = floor(pos);
//     vec3 f = fract(pos);
    
//     vec2 invNoisTexRes = 1.0/textureSize(tex, 0);
//     // float zStretch = 10.0 * invNoisTexRes;
    
//     vec2 sampleCoord = p.xy-vec2(97.0)*p.z + f.xy;//pos.xy*invNoisTexRes + p * 17.0 * invNoisTexRes;
    
//     vec2 noise = texture(tex, sampleCoord*invNoisTexRes).xy;//vec2(texture(tex, sampleCoord)[channel], texture(tex, sampleCoord + 17.0 * invNoisTexRes)[channel]);
    
//     return cube_smooth(mix(noise.x, noise.y, f.z));
// }

// float cloudShape(float alt) {
//     return billow(alt);
//     // return alt < 0.3 ? 4.0*pow(alt-0.3, 2.0) : pow(alt-0.3, 0.7);
//     // return clamp(0.0, bowl_base(alt)*0 + clamp01(pow(alt-0.3, 0.8))*0.12, 1.0);//mix(bowl_base(alt), abs(pow(alt-0.3, 0.85)) * 0.25, float(alt > 0.3)) ;
// }

float time = frameTimeCounter * CLOUD_SPEED*0.5;
vec3 displacement = vec3(time, 0.0, time*0.4);

float getCloudDensity(sampler2D tex, vec3 pos) {
    vec3 originalPos = pos;
    pos.y = length(pos + atmosCamPos) - atmosCamPos.y;
    if(pos.y > CLOUD_TOP || pos.y < CLOUD_HEIGHT) return 0.0;

    float alt = (pos.y - CLOUD_HEIGHT) / CLOUD_THICKNESS;


    pos *= CLOUD_ZOOM;
    pos += displacement;
    pos *= 1.8;

    float heightShaping = clamp01(map(alt, 0.0, 0.07, 0.0, 1.0))*clamp01(map(alt, 0.2, 1.0, 1.0, 0.0));
    float coverage_offset = 0.47 - wetness*0.24 - CLOUD_COVERAGE - cloudCoverageOffset;
    #ifdef CLOUD_KEBAB_STYLE
        coverage_offset -=  0.070;
    #endif

    #ifndef CLOUD_KEBAB_STYLE
        vec2 coverageZoom = vec2(0.018);
    #else
        vec2 coverageZoom = vec2(0.007, 0.026);
    #endif

    float coverage = texture(colortex5, pos.xz*coverageZoom).b;
        coverage = clamp01((coverage - coverage_offset) / (1.0 - clamp01(coverage_offset)));
        coverage = clamp(coverage * mix(heightShaping, 1.0, wetness*0.3), 0.0, 1.0);// - cloudShape(alt)
    if (coverage <= 0.0) return 0.0;

    float cloud = -coverage;

    cloud += cloudNoiseTex(tex, pos*3.72) * 0.3 /  1.17;
    if(cloud >= 0.0) return 0.0;

    cloud += cloudNoiseTex(tex, pos*14.88) * 0.135 * cloudGranularityFactor * CLOUD_GRANULARITY;
    // cloud += (1.0 - texture3D(colortex5, pos/128.0*14.88*3.0).r) * 0.1375;
    // cloud += (1.0 - texture3D(colortex5, pos/128.0*30.88*3.0).r) * 0.07;
    
    cloud += 0.5*0.037;//cloudNoiseTex(tex, pos*34.2) * 0.037;

    // float fade     = 1.0 - max(0.7*pow(alt-0.2, 1.0), 0.1);

    float densityHeightAlter = alt * mix(clamp01(map(alt, 0.0, 0.15, 0.0, 1.0)) * clamp01(map(alt, 0.9, 1.0, 1.0, 0.0)), 1.0, wetness);
    float distFalloff = smoothstep(100000.0, 70000.0, length(originalPos.xz));
    return max(-cloud*densityHeightAlter*22.2*CLOUD_DENSITY*cloudDensityFactor, 0.0);// - cloudShape(alt2)*0.15
}

// float getPlanarCloudDensity(sampler2D tex, vec3 pos) {
//     if(pos.y > PLANAR_CLOUD_TOP || pos.y < PLANAR_CLOUD_HEIGHT) return 0.0;
//     float time = frameTimeCounter * CLOUD_SPEED;
//     vec3 displacement = vec3(time, 0.0, time*0.4);

//     vec3 cloudPos = pos * CLOUD_ZOOM + displacement;

// 	float noise = planarCloudNoiseTex(tex, cloudPos) * 0.5;
//     	  noise += planarCloudNoiseTex(tex, cloudPos * 2.0 + displacement        ) * 0.25;
//     	  noise += planarCloudNoiseTex(tex, cloudPos * 7.0 - displacement        ) * 0.125;
//     	  noise += planarCloudNoiseTex(tex, cloudPos * 16.0 + displacement * 16.0) * 0.0625;

//     const float top = 0.004;
//     const float bottom = 0.01;

//     float horizonHeight = pos.y - PLANAR_CLOUD_HEIGHT;
//     float density = smoothstep(0.55, 0.6, noise);// * (1.0 - exp2(-bottom * horizonHeight)) * exp2(-top * horizonHeight);

//     return clamp(density, 0.0, 1.0) * CLOUD_DENSITY;
// }

float getPlanarCloudDensity(sampler2D tex, vec3 pos) {
    pos.y = length(pos + atmosCamPos) - atmosCamPos.y;
    if(pos.y > PLANAR_CLOUD_TOP || pos.y < PLANAR_CLOUD_HEIGHT) return 0.0;

    float alt = (pos.y - PLANAR_CLOUD_HEIGHT) / PLANAR_CLOUD_THICKNESS;


    pos *= PLANAR_CLOUD_ZOOM;
    pos += displacement*10.0 + vec3(10000.0, 0,0);

    float heightShaping = clamp01(map(alt, 0.0, 0.07, 0.0, 1.0))*clamp01(map(alt, 0.2, 1.0, 1.0, 0.0));
    float coverage_offset = 0.42;
    const vec3 stretchFactor = vec3(1.0, 1.0, 1.0);
    float coverage = texture(colortex5, pos.xz*0.007*stretchFactor.xz).r;
        coverage = clamp01((coverage - coverage_offset) / (1.0 - clamp01(coverage_offset)));
        coverage = clamp(coverage * heightShaping, 0.0, 1.0);// - cloudShape(alt)
    if (coverage <= 0.0) return 0.0;

    float cloud = -coverage;

    cloud += cloudNoiseTex(tex, pos*vec3(1.22, 1.0, 1.22)*stretchFactor)*0.2;
    if(cloud >= 0.0) return 0.0;

    cloud += cloudNoiseTex(tex, pos*vec3(1.88, 1.0, 1.88)*pow(stretchFactor, vec3(1.5))) * 0.085;
    
    cloud += cloudNoiseTex(tex, pos*vec3(4.88, 1.0, 4.88)*stretchFactor) * 0.057;
    cloud *= smoothstep(0.3, 0.95, texture(colortex5, pos.xz*0.008*stretchFactor.xz).r);//smoothstep(0.4, 0.6, texture(tex, pos.xz*CLOUD_ZOOM*stretchFactor.xz).r);

    // float fade     = 1.0 - max(0.7*pow(alt-0.2, 1.0), 0.1);

    // float densityHeightAlter = alt * mix(clamp01(map(alt, 0.0, 0.15, 0.0, 1.0)) * clamp01(map(alt, 0.9, 1.0, 1.0, 0.0)), 1.0, wetness);

    return smoothstep(0.0, 0.5, -cloud*0.35*cloudDensityFactor);// - cloudShape(alt2)*0.15
}



float powder(float od) {
    return 8.0 * (1.0 - 0.97 * expf(-10.0 * od));
    // return 1.0 - expf(-od * 2.0);
}


float henyeyGreensteinPhase(float cosTheta, float g) {
	const float norm = 1.0/TAU;//0.25/PI;

	float gg = g * g;
	return norm * ((1.0-gg) / pow(1.0+gg-2.0*g*cosTheta, 3.0/2.0));
}

// TODO: ADJUST THIS
float cloudPhase(float cosTheta, in float an, in float od) {
    float frontLobe = henyeyGreensteinPhase(cosTheta*an, pow(0.45, od + 1.0));
    float backLobe  = henyeyGreensteinPhase(cosTheta*an, -pow(0.45, od + 1.0));
    float peakLobe  = henyeyGreensteinPhase(cosTheta*an, pow(0.9, od + 1.0));
    return mix(mix(frontLobe, backLobe, 0.25), peakLobe, 0.15);
    // return mix(henyeyGreensteinPhase(cosTheta, 0.8*an), henyeyGreensteinPhase(cosTheta, -0.6*an), 0.4);
}

float cloudAmbientPhase(float cosTheta, in float an, in float od) {
    float frontLobe = henyeyGreensteinPhase(cosTheta*an, pow(0.35, od + 1.0));
    float backLobe  = henyeyGreensteinPhase(cosTheta*an, -pow(0.35, od + 1.0));
    return mix(frontLobe, backLobe, 0.5);
}

float lightOpticalDepth(in sampler2D tex, vec3 r0, in vec3 rd, in int steps, in float dither) {
    // vec3 endPos = r0 + rd*intersetCloudPlaneTop(r0, rd);
    // const float stepExp = 2.0;
    const float basestep = mix(150.0, 20.0, clamp01(float(steps-1)/4.0));
    const float stepExp = mix(1.6, 2.0, wetness);// mix(2.2, 1.1, clamp01(float(steps-8)/8.0));
    float dd = 22.0;//CLOUD_THICKNESS/float(steps);
    float ddPrev = dd;
    vec3 rayStep = rd;
    vec3 rayPos = r0 + rayStep*dd*dither;

    float od = 0.0;
    for(uint i = 0; i < steps; i++, rayPos += rayStep*dd) {

        if(rayPos.y > CLOUD_TOP || rayPos.y < CLOUD_HEIGHT) continue;

        ddPrev = dd;
        dd *= stepExp;

        float density = getCloudDensity(tex, rayPos);
        if(density <= 0.0) continue;
        od += density*ddPrev;
    }
    return od;
}
float ambientOpticalDepth(in sampler2D tex, vec3 r0, in vec3 rd, in int steps, in float dither) {
    float dd = 0.2*CLOUD_THICKNESS/float(steps);
    vec3 rayStep = rd;
    vec3 rayPos = r0 + rayStep*dd*dither;

    float od = 0.0;
    for(uint i = 0; i < steps; i++, rayPos += rayStep*dd) {

        float density = getCloudDensity(tex, rayPos);
        if(density <= 0.0) continue;
        od += density*dd;
    }
    return od;
}

float lightPlanarOpticalDepth(in sampler2D tex, vec3 r0, in vec3 rd, in int steps, in float dither) {
    // vec3 endPos = r0 + rd*intersetCloudPlaneTop(r0, rd);
    // const float stepExp = 2.0;
    const float basestep = mix(150.0, 20.0, clamp01(float(steps-1)/4.0));
    const float stepExp = mix(1.6, 2.0, wetness);// mix(2.2, 1.1, clamp01(float(steps-8)/8.0));
    float dd = 5.0;//CLOUD_THICKNESS/float(steps);
    float ddPrev = dd;
    vec3 rayStep = rd;
    vec3 rayPos = r0 + rayStep*dd*dither;

    float od = 0.0;
    for(uint i = 0; i < steps; i++, rayPos += rayStep*dd) {

        if(rayPos.y > PLANAR_CLOUD_TOP || rayPos.y < PLANAR_CLOUD_HEIGHT) continue;

        ddPrev = dd;
        dd *= stepExp;

        float density = getPlanarCloudDensity(tex, rayPos);
        if(density <= 0.0) continue;
        od += density*ddPrev;
    }
    return od;
}

float nightFactor = mix(1.0, 0.2, night);

const int M = 8;

vec4 cloudScattering(vec3 r0, in vec3 rd, in vec3 dirToSun, in vec3 skylight, float dither, out float cloudDist) {
    if(rd.y < 0.0) {
        cloudDist = -1.0;
        return vec4(vec3(0.0), 1.0);
    }

    // float startDist = rsi(r0 + atmosCamPos, rd, planetRadius + CLOUD_HEIGHT).y;
    // float endDist = rsi(r0 + atmosCamPos, rd, planetRadius + CLOUD_TOP).y;
    // vec2 planetDist = rsi(r0 + atmosCamPos, rd, planetRadius);

    // if(planetDist.x > 0.0) return vec4(0,0,0,1.0);
    float startDist = intersetCloudPlane(r0, rd);
    float endDist = intersetCloudPlaneTop(r0, rd);
    vec3 startPos = r0 + rd*startDist;
    vec3 endPos = r0 + rd*endDist;
    
    vec3 rayStep = (endPos - startPos)/float(CLOUD_PRIMARY_STEPS);
    float dd = length(rayStep);
    vec3 rayPos = startPos + rayStep*dither;

    float cloudDistWeights = 0.0000000001;
    cloudDist = 0.0;
    float distanceTraveled = distance(startPos, r0);
    // cloudDist = endDist;

    vec2 scatter = vec2(0.0);
    float transmittance = 1.0;
    float cosTheta = dot(rd, dirToSun);
    float cosThetaSky = dot(rd, vec3(0,1,0));
    float isoPhase = 1.0/(4.0*PI);

    //Phase and multiscatter precomp
    float[M] an;
    for(uint n = 0; n < M; n++) {
        an[n] = pow(0.5, float(n));
    }

    for(uint i = 0; i < CLOUD_PRIMARY_STEPS; i++, rayPos += rayStep, distanceTraveled += dd) {
        float density = getCloudDensity(colortex4, rayPos);
        if(density <= 0.0 || transmittance <= 1e-5) continue;
        float stepOpticalDepth = CLOUD_EXTINC * density * dd;
        float stepTransmittance = expf(-stepOpticalDepth);
        float powder = powder(CLOUD_EXTINC * density);
        float powderSun = mix(powder, 1.0, cosTheta * 0.5 + 0.5);
        float powderSky = mix(powder, 1.0, cosThetaSky * 0.5 + 0.5);

        float sunOpticalDepth = lightOpticalDepth(colortex4, rayPos, dirToSun, CLOUD_SECONDARY_STEPS, dither);
        #ifndef CLOUD_FAST_AMBIENT
        float ambientOpticalDepth = ambientOpticalDepth(colortex4, rayPos, vec3(0,1,0), CLOUD_AMBIENT_STEPS, dither);
        #else
        float ambientOpticalDepth = density * dd;
        #endif

        float S_integral = (1.0 - stepTransmittance) / CLOUD_EXTINC;
        float visibleScattering = S_integral * transmittance;

        // vec3 scatteringSun = shadowLightTransmitted*CLOUD_SCATTER*expf(-sunOpticalDepth)*powder;
        // vec3 scatteringAmbient = skylight*CLOUD_SCATTER*powder;
        // float skyVisibility = expf(-ambientOpticalDepth);

        // float scatteringAlbedo = clamp01(CLOUD_SCATTER / CLOUD_EXTINC);

        // float multipleScatteringFactor = scatteringAlbedo * 0.84;

        // vec3 scatterMultiple  = scatteringSun;
        //     scatterMultiple += scatteringAmbient * skyVisibility;
        //     scatterMultiple *= isoPhase;
        //     scatterMultiple *= multipleScatteringFactor / (1.0 - multipleScatteringFactor);

        // vec3 scattering  = scatteringSun * phase;
        //     scattering += scatteringAmbient * isoPhase * skyVisibility;
        //     scattering += scatterMultiple*1.3;
        //     scattering *= S_integral;

        for(uint n = 0; n < M; n++) {
            // float an = pow(0.5, n);

            float phase = cloudPhase(cosTheta, an[n], sunOpticalDepth);//*1.25;
            float phaseSky = cloudAmbientPhase(cosThetaSky, an[n], ambientOpticalDepth)*PI;
            vec2 S = vec2(phase*powderSun*expf(-sunOpticalDepth*CLOUD_EXTINC*an[n]),
                             phaseSky*powderSky*expf(-ambientOpticalDepth*CLOUD_EXTINC*an[n]));
            
            scatter += visibleScattering * S * an[n] * 1.7;
        }

        cloudDist += distanceTraveled * transmittance;
        cloudDistWeights += transmittance;


        transmittance *= stepTransmittance;
    }
    #ifdef PLANAR_CLOUDS_ENABLED
        float planarCloudsDist = intersetPlanarCloudPlane(r0, rd) + PLANAR_CLOUD_THICKNESS*0.5;
        rayPos = r0 + rd*planarCloudsDist;
        float density = getPlanarCloudDensity(colortex4, rayPos);

        float stepOpticalDepth = CLOUD_EXTINC * density * PLANAR_CLOUD_THICKNESS;
        float stepTransmittance = expf(-stepOpticalDepth);
        float powder = powder(CLOUD_EXTINC * density);
        float powderSun = mix(powder, 1.0, cosTheta * 0.5 + 0.5);
        float powderSky = mix(powder, 1.0, cosThetaSky * 0.5 + 0.5);

        float sunOpticalDepth = lightPlanarOpticalDepth(colortex4, rayPos, dirToSun, 3, dither);
        float ambientOpticalDepth = density * PLANAR_CLOUD_THICKNESS;

        float S_integral = (1.0 - stepTransmittance) / CLOUD_EXTINC;
        float visibleScattering = S_integral * transmittance;


        for(uint n = 0; n < M; n++) {
            // float an = pow(0.5, n);

            float phase = cloudPhase(cosTheta, an[n], sunOpticalDepth);
            float phaseSky = cloudAmbientPhase(cosThetaSky, an[n], ambientOpticalDepth)*PI;
            vec2 S = vec2(phase*powderSun*expf(-sunOpticalDepth*CLOUD_EXTINC*an[n]),
                             phaseSky*powderSky*expf(-ambientOpticalDepth*CLOUD_EXTINC*an[n]));
            
            scatter += visibleScattering * S * an[n] * 1.7;
        }
        transmittance *= stepTransmittance;

        // if(density >= 1e-5) {
            cloudDist += (distanceTraveled + planarCloudsDist + PLANAR_CLOUD_THICKNESS) * transmittance;
            cloudDistWeights += transmittance;
        // }
        
    #endif
    
    vec3 finalScatter = scatter.x * nightFactor * CLOUD_SCATTER * shadowLightTransmitted;
            finalScatter += scatter.y * skylight * CLOUD_SCATTER;
    
    cloudDist /= cloudDistWeights;


    return clamp(vec4(finalScatter, transmittance), 0.0, 999999.0);
}