+
+#define testDeltaT (3 * 86400. * 1e-6) // usec
+
+vec3 differential(vec3 p, float deltaT) {
+    if (deltaT == 0)
+        return p;
+
+    float phi = atan(p.x, p.z);
+    float theta = acos(p.y);
+
+    float sin2l = sin(theta);
+    sin2l *= sin2l;
+    float sin4l = sin2l * sin2l;
+    phi -= deltaT * (0.343 * sin2l + 0.474 * sin4l);
+
+    return vec3(sin(theta) * sin(phi), cos(theta), sin(theta) * cos(phi));
+}
 
 float intersectPlane(const vec4 quat, const vec4 vecin, bool hideBack) {
     vec3 altnormal = rotate_vector(quat, vec3(0., 0., 1.));

    bool onDisk = radius2 <= 1;

    float factor;
    vec3 hitPoint = vec3(0.), rotatedHitPoint = vec3(0.);

    if (onDisk) {
        hitPoint = vec3(up1.x, up1.y, sqrt(1. - radius2));
        rotatedHitPoint = rotate_vector_inverse(cameraDifferenceRotationQuat, hitPoint);
        factor = 1.;
        gl_FragDepth = 0.5 - hitPoint.z * CLIP_SCALE_NARROW;
    } else {

    if (rotatedHitPoint.z <= 0.) { // off-limb or back
        hitPoint = vec3(up1.x, up1.y, intersectPlane(cameraDifferenceRotationQuat, up1, onDisk));
        rotatedHitPoint = rotate_vector_inverse(cameraDifferenceRotationQuat, hitPoint);

        if (length(rotatedHitPoint) <= 1) // central disk
            discard;

        if (calculateDepth != 0) // intersecting Eufhoria planes
            gl_FragDepth = 0.5 - hitPoint.z * CLIP_SCALE_WIDE;
