working implementation for forward/backwards transform + shaders, how…

…ever some issues with display outside of the projection ellipse

src/core/StelCore.cpp

@@ -432,6 +432,9 @@ StelProjectorP StelCore::getProjection(StelProjector::ModelViewTranformP modelVi
 		case ProjectionHammer:
 			prj = StelProjectorP(new StelProjectorHammer(modelViewTransform));
 			break;
+		case ProjectionMollweide:
+			prj = StelProjectorP(new StelProjectorMollweide(modelViewTransform));
+			break;
 		case ProjectionCylinder:
 			prj = StelProjectorP(new StelProjectorCylinder(modelViewTransform));
 			break;

src/core/StelProjectorClasses.cpp

@@ -552,76 +552,86 @@ QString StelProjectorMollweide::getDescriptionI18() const {
     return q_("The Mollweide projection is an equal-area map projection introduced by Karl Mollweide in 1805.");
 }
 
-bool StelProjectorMollweide::forward(Vec3f &v) const {
-    // Convert from 3D sphere to (x,y) in projection space.
+bool StelProjectorMollweide::forward(Vec3f &v) const 
+{
     const float r = std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
-    // Compute latitude (φ) and longitude (λ)
-    const float phi = std::asin(v[1] / r);
-    const float lambda = std::atan2(v[0], -v[2]);
-    // Solve for theta in: 2θ + sin(2θ) = π sinφ
-    float theta = phi; // initial guess
-    const float piSinPhi = static_cast<float>(M_PI) * std::sin(phi);
-    const int maxIter = 10;
-    const float tol = 1e-7f;
-    for (int i = 0; i < maxIter; ++i) {
-        float f = 2.f*theta + std::sin(2.f*theta) - piSinPhi;
-        float fprime = 2.f + 2.f*std::cos(2.f*theta);
-        float dtheta = f / fprime;
-        theta -= dtheta;
-        if (std::fabs(dtheta) < tol)
-            break;
+    const float lambda = std::atan2(v[0], -v[2]); // Longitude
+    const float sinDelta = v[1] / r; // Latitude (delta) sine
+    const float delta = std::asin(sinDelta);
+
+    // Solve 2θ + sin(2θ) = π sinδ using Newton-Raphson
+    float theta = delta;
+    for (int i=0; i<5; ++i) {
+        float f = 2.f*theta + sin(2.f*theta) - M_PIf * sinDelta;
+        float df = 2.f + 2.f*cos(2.f*theta);
+        theta -= f / df;
     }
-    // Compute projection coordinates:
-    // widthStretch is a member (as in the Hammer code)
-    float x = (2.f * static_cast<float>(M_SQRT2) / static_cast<float>(M_PI)) * lambda * std::cos(theta) * widthStretch;
-    float y = static_cast<float>(M_SQRT2) * std::sin(theta);
-    v[0] = x;
-    v[1] = y;
-    v[2] = r; // store original radius if needed
+
+    const float cosTheta = cos(theta);
+    v[0] = (2.0f * M_SQRT2 / M_PIf) * lambda * cosTheta * widthStretch;
+    v[1] = M_SQRT2 * sin(theta);
+    v[2] = r;
     return true;
 }
 
+
 bool StelProjectorMollweide::backward(Vec3d &v) const {
-    // Inverse: from projection (x,y) back to sphere.
-    v[0] /= widthStretch;
-    // Recover theta from y = √2 sinθ
-    double theta = std::asin(v[1] / static_cast<double>(M_SQRT2));
-    double cosTheta = std::cos(theta);
-    // Avoid division by zero near the poles
-    if (std::fabs(cosTheta) < 1e-10)
+	v[0] /= widthStretch;
+    const double x = v[0];
+    const double y = v[1];
+
+    // Check if point lies within Mollweide ellipse bounds
+    const bool ok = (x*x + 4.0*y*y <= 8.0);
+
+    // Compute theta (exit early if invalid)
+    if (std::abs(y) > M_SQRT2) 
         return false;
-    // Recover λ from x = (2√2/π) λ cosθ
-    double lambda = (v[0] * static_cast<double>(M_PI)) / (2.0 * static_cast<double>(M_SQRT2) * cosTheta);
-    // Recover φ from: 2θ + sin(2θ) = π sinφ  =>  φ = arcsin((2θ + sin2θ)/π)
-    double phi = std::asin((2.0*theta + std::sin(2.0*theta)) / static_cast<double>(M_PI));
-    double cosPhi = std::cos(phi);
-    // Convert spherical (φ,λ) to 3D vector (using same convention as in the Hammer code)
-    v[0] = cosPhi * std::sin(lambda);
-    v[1] = std::sin(phi);
-    v[2] = - cosPhi * std::cos(lambda);
-    return true;
+    const double theta = asin(y / M_SQRT2);
+
+    // Reconstruct spherical coordinates (lambda, delta)
+    const double cosTheta = std::cos(theta);
+    const double lambda = (x * M_PI) / (2.0 * M_SQRT2 * cosTheta);
+    const double sinDelta = (2.0 * theta + sin(2.0*theta)) / M_PI;
+    const double delta = std::asin(sinDelta);
+
+    // Convert to 3D Cartesian
+    const double cosDelta = std::cos(delta);
+    v[0] = cosDelta * std::sin(lambda);
+    v[1] = sinDelta;
+    v[2] = -cosDelta * std::cos(lambda);
+    return ok;
 }
 
 QByteArray StelProjectorMollweide::getForwardTransformShader() const
 {
-	return modelViewTransform->getForwardTransformShader() + R"(
+    return modelViewTransform->getForwardTransformShader() + R"(
 #line 1 102
 // Forward transform shader snippet:
 uniform float PROJECTOR_FWD_widthStretch;
-vec3 projectorForwardTransform(vec3 v) {
+vec3 projectorForwardTransform(vec3 v) 
+{
+    const float M_SQRT2 = 1.41421356;
+    const float PI = 3.14159265;
+    float widthStretch = PROJECTOR_FWD_widthStretch;
+    
     float r = length(v);
-    float phi = asin(v.y / r);
-    float lambda = atan(v.x, -v.z);
-    float theta = phi;
-    // Iterative solution for theta (10 iterations)
-    for (int i = 0; i < 10; i++) {
-        float f = 2.0 * theta + sin(2.0 * theta) - 3.14159265 * sin(phi);
-        float fprime = 2.0 + 2.0 * cos(2.0 * theta);
-        theta -= f / fprime;
+    float lambda = atan(v[0], -v[2]);
+    float sinDelta = v[1] / r;
+    float delta = asin(sinDelta);
+    
+    // Newton-Raphson iterations for theta
+    float theta = delta;
+    for (int i=0; i<5; i++) {
+        float f = 2.0*theta + sin(2.0*theta) - PI * sinDelta;
+        float df = 2.0 + 2.0*cos(2.0*theta);
+        theta = theta - f/df;
     }
-    float x = (2.82842712 / 3.14159265) * lambda * cos(theta) * PROJECTOR_FWD_widthStretch; // 2√2 = 2.82842712
-    float y = 1.41421356 * sin(theta); // √2 = 1.41421356
-    return vec3(x, y, r);
+    
+    float cosTheta = cos(theta);
+    v[0] = (2.0 * M_SQRT2 / PI) * lambda * cosTheta * widthStretch;
+    v[1] = M_SQRT2 * sin(theta);
+    v[2] = r;
+    return v;
 }
 #line 1 0
 )";
@@ -633,22 +643,28 @@ QByteArray StelProjectorMollweide::getBackwardTransformShader() const
 #line 1 103
 // Backward transform shader snippet:
 uniform float PROJECTOR_FWD_widthStretch;
-vec3 projectorBackwardTransform(vec3 v, out bool ok) {
-    v.x /= PROJECTOR_FWD_widthStretch;
-    float theta = asin(v.y / 1.41421356);
+vec3 projectorBackwardTransform(vec3 v, out bool ok) 
+{
+    float widthStretch = PROJECTOR_FWD_widthStretch;
+    v[0] /= widthStretch;
+    float x = v[0];
+    float y = v[1];
+    
+    // Validity check: lies within Mollweide ellipse (x^2 +4y^2 <=8)
+    ok = (x*x + 4.0*y*y <= 8.0);
+
+    float theta = asin(y / 1.41421356); // y ~ M_SQRT2*sin(theta)
     float cosTheta = cos(theta);
-    if (abs(cosTheta) < 0.00001) {
-        ok = false;
-        return v;
-    }
-    float lambda = (v.x * 3.14159265) / (2.82842712 * cosTheta);
-    float phi = asin((2.0 * theta + sin(2.0 * theta)) / 3.14159265);
-    float cosPhi = cos(phi);
-    float x = cosPhi * sin(lambda);
-    float y = sin(phi);
-    float z = - cosPhi * cos(lambda);
-    ok = true;
-    return vec3(x, y, z);
+    float lambda = (x * 3.14159265) / (2.0 * 1.41421356 * cosTheta);
+    
+    float sinDelta = (2.0*theta + sin(2.0*theta)) / 3.14159265;
+    float delta = asin(sinDelta);
+    
+    float cosDelta = cos(delta);
+    v[0] = cosDelta * sin(lambda);
+    v[1] = sinDelta;
+    v[2] = -cosDelta * cos(lambda);
+    return v;
 }
 #line 1 0
 )";

src/core/StelProjectorClasses.hpp

@@ -113,19 +113,32 @@ class StelProjectorHammer : public StelProjector
 	}
 };
 
-class StelProjectorMollweide : public StelProjector 
-{
+class StelProjectorMollweide : public StelProjector {
 public:
 	StelProjectorMollweide(ModelViewTranformP func) : StelProjector(func) {}
 	QString getNameI18() const override;
 	QString getDescriptionI18() const override;
-	float getMaxFov() const override { return 180.f; }
+	float getMaxFov() const override { return 185.f; }
 	bool forward(Vec3f &v) const override;
 	bool backward(Vec3d &v) const override;
+	// float fovToViewScalingFactor(float fov) const override;
+	//float viewScalingFactorToFov(float vsf) const override;
+	//float deltaZoom(float fov) const override;
 	QByteArray getForwardTransformShader() const override;
 	QByteArray getBackwardTransformShader() const override;
 protected:
+// TODO: take a look at the discontinuities => might be the reason for weird display glitches at high FOV
 	bool hasDiscontinuity() const override { return true; }
+	bool intersectViewportDiscontinuityInternal(const Vec3d& p1, const Vec3d& p2) const override {
+		return p1[0] * p2[0] < 0 && !(p1[2] < 0 && p2[2] < 0);
+	}
+	bool intersectViewportDiscontinuityInternal(const Vec3d& capN, double capD) const override {
+		static const SphericalCap cap1(1, 0, 0);
+		static const SphericalCap cap2(-1, 0, 0);
+		static const SphericalCap cap3(0, 0, -1);
+		SphericalCap cap(capN, capD);
+		return cap.intersects(cap1) && cap.intersects(cap2) && cap.intersects(cap3);
+	}
 };
 
 class StelProjectorCylinder : public StelProjector