Fix eclipse artifacts

This reflects the fix in CalcMySky commit
afe131c57b14f984e3cdcd090fa905695b02ffce.

Fixes #3729

atmosphere/default/shaders/double-scattering-eclipsed/precomputation/0/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/double-scattering-eclipsed/precomputation/1/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/double-scattering-eclipsed/precomputation/2/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/double-scattering-eclipsed/precomputation/3/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/double-scattering-eclipsed/precomputed/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/eclipsed-zero-order-scattering/0/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/eclipsed-zero-order-scattering/1/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/eclipsed-zero-order-scattering/2/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/eclipsed-zero-order-scattering/3/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/light-pollution/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/multiple-scattering/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/single-scattering-eclipsed/on-the-fly/0/aerosols/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/single-scattering-eclipsed/on-the-fly/0/molecules/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/single-scattering-eclipsed/on-the-fly/1/aerosols/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/single-scattering-eclipsed/on-the-fly/1/molecules/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/single-scattering-eclipsed/on-the-fly/2/aerosols/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/single-scattering-eclipsed/on-the-fly/2/molecules/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/single-scattering-eclipsed/on-the-fly/3/aerosols/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/single-scattering-eclipsed/on-the-fly/3/molecules/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/single-scattering-eclipsed/precomputation/0/aerosols/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/single-scattering-eclipsed/precomputation/0/molecules/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/single-scattering-eclipsed/precomputation/1/aerosols/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/single-scattering-eclipsed/precomputation/1/molecules/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/single-scattering-eclipsed/precomputation/2/aerosols/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/single-scattering-eclipsed/precomputation/2/molecules/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }
 

atmosphere/default/shaders/single-scattering-eclipsed/precomputation/3/aerosols/common-functions.frag

@@ -186,8 +186,12 @@ float circlesIntersectionArea(float R1, float R2, float d)
     if(d>=R1+R2) return 0.;
 
     // Return area of the lens with radii R1 and R2 and offset d
-    return sqr(R1)*acos(clamp( (sqr(d)+sqr(R1)-sqr(R2))/(2*d*R1) ,-1.,1.)) +
-           sqr(R2)*acos(clamp( (sqr(d)+sqr(R2)-sqr(R1))/(2*d*R2) ,-1.,1.)) -
+    // Note: R1^2-R2^2 must be computed before adding d^2 to the result,
+    // otherwise catastrophic cancellation will occur when R1≈R2 and d≈0,
+    // which will result in visual artifacts.
+    CONST float dRR = sqr(R1)-sqr(R2);
+    return sqr(R1)*acos(clampCosine( (sqr(d)+dRR)/(2*d*R1) )) +
+           sqr(R2)*acos(clampCosine( (sqr(d)-dRR)/(2*d*R2) )) -
            0.5*sqrt(max( (-d+R1+R2)*(d+R1-R2)*(d-R1+R2)*(d+R1+R2) ,0.));
 }