From f2a7806f5f0c2aa6d17dc7c2811a7dd7432c3f90 Mon Sep 17 00:00:00 2001
From: "Timur A. Fatkhullin" <amirych@gmail.com>
Date: Sun, 20 Apr 2025 01:26:45 +0300
Subject: [PATCH] ...

---
 cxx/mount_pz.h | 84 ++++++++++++++++++++++++++++++++++++++++++++------
 1 file changed, 74 insertions(+), 10 deletions(-)

diff --git a/cxx/mount_pz.h b/cxx/mount_pz.h
index 483187b..3dc0ad2 100644
--- a/cxx/mount_pz.h
+++ b/cxx/mount_pz.h
@@ -61,30 +61,94 @@ private:
 };
 
 
-class MccElliipsePZ
+/*
+ * a general planar ellipse equation:
+ *     A*(x-xc)^2 + B*(x-xc)(y-yc) + C*(y-yc)^2 = 1
+ *
+ *     A = cos(t)^2/a^2 + sin(t)^2/b^2
+ *     B = sin(2*t)/a^2 - sin(2*t)/b^2
+ *     C = cos(t)^2/b^2 + sin(t)^2/a^2
+ *
+ *     t - angle between X-axis and big semi-axis (a)
+ *
+ *
+ *    Ellipse on unit sphere (Az, Alt):
+ *        x^2/a^2 + y^2/b^2 = 1,
+ *        where a = tan(alpha), b = tan(beta),
+ *              a and b - big and small spherical semi-axis
+ *        also:
+ *             let delta is a spherical distance between ellipse focuses, then
+ *             d = tan(delta) and b^2 = (a^2 - d^2)/(1 + d^2)
+ *
+ *    tangent coordinates:
+ *    x = tan(Az)
+ *    y = tan(Alt)*sqrt(1+tan(Az)^2)
+ *
+ *
+ *    --------------------------------------------
+ *
+ *    let P - point with (Az_P, Zd_P),
+ *    (Az_C, Zd_C) - center of ellipse, a and b big and small semi-axis,
+ *    vec_a and vec_b - unit vectors along a an b (it lie in tangent surface to point C!!!), then
+ *
+ *    ((vec_P-vec_C)*vec_b)^2/a^2 + ((vec_P-vec_C)*vec_b)^2/b^2 <= 1.0
+ *    {((vec_P-vec_C)*vec_b)^2/tan(a)^2 + ((vec_P-vec_C)*vec_b)^2/tan(b)^2 <= 1.0}
+ *    * - dot-product
+ *
+ *    vec_P = (sin(Zd_P)*cos(Az_P), sin(Zd_P)*sin(Az_P), cos(Zd_P))
+ *    vec_C = (sin(Zd_C)*cos(Az_C), sin(Zd_C)*sin(Az_C), cos(Zd_C))
+ *
+ */
+class MccEllipsePZ
 {
 public:
-    MccElliipsePZ(const MccCoordinate& xc, const MccCoordinate& yc, const MccCoordinate& a, const MccCoordinate& b)
-        : _xc(xc), _yc(yc), _a(a), _b(b), _a2(a * a), _b2(b * b)
+    MccEllipsePZ(const MccCoordinate& xc,
+                 const MccCoordinate& yc,
+                 const MccCoordinate& a,
+                 const MccCoordinate& b,
+                 const MccCoordinate& theta = 0.0)
+        : _xc(xc), _yc(yc), _a(a), _b(b), _theta(theta), _tanXc(std::tan(xc)), _tanYc(std::tan(yc))
     {
+        _tanYc *= std::sqrt(1.0 + _tanXc * _tanXc);
+
+        auto _tan2A = tan(a);
+        auto _tan2B = tan(b);
+
+        _tan2A *= _tan2A;
+        _tan2B *= _tan2B;
+
+        auto ct = cos(theta);
+        auto ct2 = ct * ct;
+        auto st = sin(theta);
+        auto st2 = st * st;
+        auto sin2t = sin(2.0 * theta);
+
+
+        cxx = ct2 / _tan2A + st2 / _tan2B;
+        cyy = st2 / _tan2A + ct2 / _tan2B;
+
+        cxy = sin2t / _tan2A - sin2t / _tan2B;
     }
 
     // circle
-    MccElliipsePZ(const MccCoordinate& xc, const MccCoordinate& yc, const MccCoordinate& a)
-        : mcc::MccElliipsePZ(xc, yc, a, a)
+    MccEllipsePZ(const MccCoordinate& xc, const MccCoordinate& yc, const MccCoordinate& a)
+        : mcc::MccEllipsePZ(xc, yc, a, a)
     {
     }
 
 
 private:
-    double _xc, _yc, _a, _b, _a2, _b2;
+    double _xc, _yc, _a, _b, _theta;
+    double _tanXc, _tanYc, cxx, cxy, cyy;
 
-    bool inZoneImpl(const MccCoordinate& alt, const MccCoordinate& zd)
+    bool inZoneImpl(const MccCoordinate& x, const MccCoordinate& y)
     {
-        auto x2 = alt - _xc;
-        auto y2 = zd - _yc;
+        auto tanX = tan(x);
+        auto tanY = tan(y) * sqrt(1.0 + tanX * tanX);
+        auto xr = tanX - _tanXc;
+        auto yr = tanY - _tanYc;
 
-        return (x2 * x2 / _a2 + y2 * y2 / _b2) <= 1.0;
+        return (cxx * xr * xr + cxy * xr * yr + cyy * yr * yr) <= 1.0;
     }
 };