vec3 physicalSun(in vec3 sceneDir) {
    float cosTheta = dot(sceneDir, sunVector);
    float centerToEdge = clamp(acos(cosTheta) / sunAngularRadius, 0.0, 1.0);
    float cosSunRadius = cos(sunAngularRadius);
    if(cosTheta < cosSunRadius) return vec3(0.0);

    vec3 a = vec3(0.397, 0.503, 0.652);

    float  mu = sqrt(1.0 - centerToEdge*centerToEdge);
    vec3 factor = 1.0 - (1.0 - pow(vec3(mu), a));
    vec3 finalLuminance = sunLuminance * cube(factor);

    return finalLuminance * rcp(pi);
}

vec3 physicalMoon(in vec3 sceneDir) {
	vec3 viewPos = vec3(0.0, planetRadius, 0.0);
    vec2 dists = RSI(-moonVector, sceneDir, moonAngularRadius);

    vec3 moonNormal = fNormalize(sceneDir*dists.x - moonVector);

    vec3 diffuse = moonColor * dot(sunVector, moonNormal) / pi;

	return dists.y > 0.0 ? (diffuse * sunIlluminance) : vec3(0.0);
}