float vc_getLD(in vec3 rpos, const int steps, in vec3 lvec) {
    float density   = 1.0;

    float stepMod   = 1.0-(linStep(rpos.y, vc_lowEdge, vc_highEdge)*0.97)*(timeNoon);
    vec3 dir        = normalize(mat3(gbufferModelViewInverse)*lvec);
    float stepSize  = (s_vcThickness/steps)*stepMod;
    vec3 rstep      = dir*stepSize;
        rpos        += rstep;

    float transmittance = 0.0;

    for (int i = 0; i<steps; i++, rpos += rstep) {
        if (rpos.y>vc_highEdge || rpos.y<vc_lowEdge) break;
        float coverage = vc_getCoverage(rpos);
        if (coverage <= 0.0) continue;

        #ifdef s_vcFastLighting
            float oD    = vc_getShapeLQ(rpos, coverage)*stepSize;
        #else
            float oD    = vc_getShape(rpos, coverage)*stepSize;
        #endif

        transmittance += oD;
    }
    return transmittance*density;
}
float vc_getLD(in vec3 rpos, const int steps) {
    return vc_getLD(rpos, steps, mix(vec.light, vec.up, timeLightTransition));
}

float vc_getScatter(in float lD, in float powder, in float vdotl) {
    float transmittance = exp2(-lD);
    float phase         = c_miePhase(vdotl)*0.93+0.07;

    return max(powder*phase*transmittance*invLog2, 0.0);
}

#ifndef refPass
float vc_getLDexp(in vec3 rpos, const int steps, in vec3 lvec) {
    float density   = 1.0;

    float stepMod   = 1.0-(linStep(rpos.y, vc_lowEdge, vc_highEdge)*0.9)*(timeNoon);
    vec3 dir        = normalize(mat3(gbufferModelViewInverse)*lvec);    //this can be moved outside of the raymarching loop, because matrix expensive
    float stepSize  = (33.3/steps)*stepMod;
    vec3 rstep      = dir;

    float transmittance = 0.0;

    for (int i = 0; i<steps; i++, rpos += dir * stepSize) {
        if (rpos.y>vc_highEdge || rpos.y<vc_lowEdge) break;
        float coverage = vc_getCoverage(rpos);
        if (coverage <= 0.0) continue;

        #ifdef s_vcFastLighting
            float oD    = vc_getShapeLQ(rpos, coverage)*stepSize;
        #else
            float oD    = vc_getShape(rpos, coverage)*stepSize;
        #endif

        transmittance += oD;
		stepSize *= 2.2;
    }
    return transmittance*density;
}
float vc_getLDexp(in vec3 rpos, const int steps) {
    return vc_getLDexp(rpos, steps, mix(vec.light, vec.up, timeLightTransition));   //move lvec mix outside of loop here
}
float vc_getScatter(in float lD, in float powder, in float vdotl, in float ldM, in float phaseM) {
    float transmittance = exp2(-lD*ldM);
    float phase         = c_miePhase(vdotl, phaseM)*0.93+0.07;

    return max(powder*phase*transmittance*invLog2, 0.0);
}

float vc_getSkylight(in vec3 rpos) {
    #ifdef s_vcSkylightHQ
        float lD    = vc_getLD(rpos, s_vcSkylightSamples, vec.up);
        return exp2(-lD*invLog2*0.0114);
    #else
        float light = sstep(rpos.y, vc_lowEdge, vc_highEdge)*0.85+0.15;
            light   = sqrt(light);
            //light  *= 1.0-saturate(oD*0.2)*0.2;

        return light;
    #endif
}

void vc_scatter0(inout vec2 scatter, in float oD, in vec3 rpos, in float scattercoeff, in float vdotl, in float transmittance, in float stept, in float mie, in float pmie) {
    float lD    = vc_getLD(rpos, 3);
    float integral = scatterIntegral(stept, 1.0);
    float powder = 1.0-exp(-oD);
        powder  = mix(powder, 1.0, pmie);

    float s  = vc_getScatter(lD*0.1, powder, vdotl)*integral*scattercoeff;

    scatter.x += s*transmittance*0.9;
    scatter.y += vc_getSkylight(rpos)*integral*transmittance;
}

void vc_scatter(inout vec2 scatter, in float oD, in vec3 rpos, in float scattercoeff, in float vdotl, in float transmittance, in float stept, in float mie, in float pmie) {
    float lD    = vc_getLD(rpos, s_vcLightSamples);
    float integral = scatterIntegral(stept, 1.0);
    float powder = 1.0-exp(-oD);
        powder  = mix(powder, 1.0, pmie);

    float s  = vc_getScatter(lD*0.1, powder, vdotl)*integral*scattercoeff;

    scatter.x += s*transmittance;
    scatter.y += vc_getSkylight(rpos)*integral*transmittance;
}

void vc_multiscatter0(inout vec2 scatter, in float oD, in vec3 rpos, in float scattercoeff, in float vdotl, in float transmittance, in float stept, in float mie, in float pmie) {
    float lD    = vc_getLDexp(rpos, s_vcLightSamples);
    float integral = scatterIntegral(stept, 1.0);
    float powder = exp(-oD -lD);
        powder  = mix(1.0-powder, 1.0+powder*0.5, pmie);
    float s     = 0.0;
    float n     = 0.0;

    for (int i = 0; i<3; i++) {
        float scatterM  = pow(0.5, float(i));
        float ldM       = pow(0.15, float(i));
        float phaseM    = pow(0.85, float(i));

        scattercoeff *= scatterM;

        s  += vc_getScatter(lD, powder, vdotl, ldM, phaseM)*integral*scattercoeff;
        n  += scatterM;
    }

    scatter.x += s*transmittance*(1.0/n);
    scatter.y += vc_getSkylight(rpos)*integral*transmittance;
}

void vc_multiscatter1(inout vec2 scatter, in float oD, in vec3 rpos, in float scattercoeff, in float vdotl, in float transmittance, in float stept, in float mie, in float pmie) {
    float lD    = vc_getLDexp(rpos, s_vcLightSamples)*1.11;
    float integral = scatterIntegral(stept, 1.11);
    float s     = 0.0;
    float n     = 0.0;

    for (int i = 0; i<3; i++) {
        float scatterM  = pow(0.5, float(i));
        float ldM       = pow(0.15, float(i));
        float phaseM    = pow(0.6, float(i));
        float powder    = exp(-oD -lD*ldM);
            powder      = mix(1.0-powder, 1.0+powder*0.5, pmie);

        scattercoeff *= scatterM;

        s  += vc_getScatter(lD, powder, vdotl, ldM, phaseM)*integral*scattercoeff;
        n  += scatterM;
    }

    scatter.x += s*transmittance*(1.0/n);
    scatter.y += vc_getSkylight(rpos)*integral*transmittance;
}
#endif

void vc_scatterRef(inout vec2 scatter, in float oD, in vec3 rpos, in float scattercoeff, in float vdotl, in float transmittance, in float stept, in float mie) {
    float lD    = vc_getLD(rpos, 3);
    float integral = scatterIntegral(stept, 1.0);
    float powder = 1.0-exp(-oD);
        powder  = mix(powder, 1.0, vdotl*0.5+0.5);

    float s  = vc_getScatter(lD*0.1, powder, vdotl)*integral*scattercoeff;

    float light = sstep(rpos.y, vc_lowEdge, vc_highEdge)*0.85+0.15;
        light   = sqrt(light);

    scatter.x += s*transmittance;
    scatter.y += light*integral*transmittance;
}