From b4753e66e29af0b6087f591838069deb9689b78f Mon Sep 17 00:00:00 2001
From: jaseg <git@jaseg.de>
Date: Sun, 26 Mar 2023 00:40:43 +0100
Subject: WIP

---
 svg-flatten/geom_test.py       | 60 ++++++++++++++++++++++++++++++++++++++++++
 svg-flatten/include/geom2d.hpp | 31 +++++++++++++++++++++-
 2 files changed, 90 insertions(+), 1 deletion(-)
 create mode 100644 svg-flatten/geom_test.py

diff --git a/svg-flatten/geom_test.py b/svg-flatten/geom_test.py
new file mode 100644
index 0000000..d49b500
--- /dev/null
+++ b/svg-flatten/geom_test.py
@@ -0,0 +1,60 @@
+#!/usr/bin/env python3
+
+from math import *
+
+def calc(mat):
+    [xx, yx], [xy, yy] = mat
+
+    a = xx**2 + xy**2
+    b = xx*yx + xy*yy
+    c = yy**2 + yx**2
+
+    print(f'{a=:.2f} {c=:.2f} {b=:.2f}')
+    tan_2_alpha = 2*b/(a-c)
+    print(f'atan2={atan2(2*b, a-c)/pi:.2f}*pi')
+    #tan_alpha = tan_2_alpha / (1 + sqrt(1 + tan_2_alpha**2)) # FIXME: bounds?
+    cos_2_alpha = 1/sqrt(1 + tan_2_alpha**2)
+    sin_2_alpha = tan_2_alpha / sqrt(1 + tan_2_alpha**2)
+    print(f'tan(2a)={tan_2_alpha:.2f} cos(2a)={cos_2_alpha:.2f} sin(2a)={sin_2_alpha:.2f}')
+    cos_alpha = sqrt((1 + cos_2_alpha)/2)
+    sin_alpha = sqrt((1 - cos_2_alpha)/2)
+    print(f'cos(a)={cos_alpha:.2f} sin(a)={sin_alpha:.2f}')
+
+    for sgn_cos, sgn_sin in [(-1, -1), (-1, 1), (1, -1), (1, 1)]:
+        p = xx * sgn_cos * cos_alpha + yx * sgn_sin * sin_alpha
+        q = xy * sgn_cos * cos_alpha + yy * sgn_sin * sin_alpha
+        dist = hypot(p, q)
+        yield dist
+
+def gen(sx, sy, m, theta):
+    xx = sx * cos(theta)
+    xy = sx * sin(theta)
+    yy = sy * (cos(theta) + m * sin(theta))
+    yx = sy * (m * cos(theta) - sin(theta))
+
+    mat = [xx, yx], [xy, yy]
+    return mat
+
+
+for sx, sy in [
+        (1, 0.9),
+        (1, 1.0),
+        (1, 1.1),
+        (0.9, 1),
+        (1.0, 1),
+        (1.1, 1)]:
+    for m in [0, 0.1, 1, 10]:
+        for theta in [0, pi/8, pi/4, pi/3, pi/2, pi, 3*pi/4]:
+            print(f'{sx=:.1f} {sy=:.1f} {m=:.1f} theta={theta/pi:.2f}*pi |', end=' ')
+            mat = gen(sx, sy, m, theta)
+
+            try:
+                dists = list(calc(mat))
+                str_dists = ' '.join(f'{x:.2f}' for x in dists)
+                print(f'[{str_dists}] | min={min(dists):.2f} max={max(dists):.2f}')
+            except:
+                print('E')
+            break
+        break
+    break
+
diff --git a/svg-flatten/include/geom2d.hpp b/svg-flatten/include/geom2d.hpp
index 47cf3be..9650575 100644
--- a/svg-flatten/include/geom2d.hpp
+++ b/svg-flatten/include/geom2d.hpp
@@ -113,9 +113,38 @@ namespace gerbolyze {
 
             double doc2phys_skew(double dist_doc) {
                 /* https://math.stackexchange.com/a/3521141 */
+                /* https://stackoverflow.com/a/70381885 */
                 /* xx yx x0
                  * xy yy y0 */
-                s_x = sqrt();
+                double s_x = sqrt(xx*xx + xy*xy);
+
+                if (xx == 0 && xy == 0) {
+                    return std::numeric_limits<double>::infinity;
+                }
+
+                double theta = atan2(xy, xx);
+                double f = (xx*yy - xy*yx);
+
+                if (f == 0) {
+                    return std::numeric_limits<double>::infinity;
+                }
+
+                double m = (xx*yx + yy*xy) / f;
+
+                double f = xx + m*xy;
+                double s_y = 0;
+
+                if (f == 0) {
+                    f = m*xx - xy;
+                    if (f == 0) {
+                        return std::numeric_limits<double>::infinity;
+                    }
+                    s_y = yx*s_x / f;
+                } else {
+                    s_y = yy*s_x / f;
+                }
+
+                return s_x - s_y > 
             }
 
             double doc2phys_min(double dist_doc) {
-- 
cgit