From 827fb28967d50ea1fbaab7b7a716502cbec083f9 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Sat, 8 Sep 2018 21:47:39 +0200
Subject: [PATCH] Corrected angle computationw

---
 methods/entropic_eq_srt.py                 |  3 ++
 phasefield/contact_angle_circle_fitting.py | 38 +++++++++++++---------
 2 files changed, 25 insertions(+), 16 deletions(-)

diff --git a/methods/entropic_eq_srt.py b/methods/entropic_eq_srt.py
index 448be382..7a7b039c 100644
--- a/methods/entropic_eq_srt.py
+++ b/methods/entropic_eq_srt.py
@@ -6,6 +6,9 @@ from lbmpy.methods.conservedquantitycomputation import DensityVelocityComputatio
 
 
 class EntropicEquilibriumSRT(AbstractLbMethod):
+    """Equilibrium from 'Minimal entropic kinetic models for hydrodynamics'
+    Ansumali, S. ; Karlin, I. V;  Öttinger, H. C, (2003)
+    """
     def __init__(self, stencil, relaxation_rate, force_model, conserved_quantity_calculation):
         super(EntropicEquilibriumSRT, self).__init__(stencil)
 
diff --git a/phasefield/contact_angle_circle_fitting.py b/phasefield/contact_angle_circle_fitting.py
index d5bae464..52681f6a 100644
--- a/phasefield/contact_angle_circle_fitting.py
+++ b/phasefield/contact_angle_circle_fitting.py
@@ -122,19 +122,25 @@ def liquid_lens_neumann_angles(concentration, drop_phase_idx=2, enclosing_phase1
 
     y1, y2 = (i[1] for i in intersections)
     assert np.abs(y1 - y2) < 0.1, "Liquid lens is not aligned in y direction (rotated?)"
-    y_horizontal = (y1 + y2) / 2
-
-    xs = np.sqrt(circles[0].radius ** 2 - (y_horizontal - circles[0].center[1]) ** 2)
-    angle = np.rad2deg(np.arctan(-xs / np.sqrt(circles[0].radius ** 2 - xs ** 2))) % 180
-    angle = 180 - angle
-    neumann3_upper = angle
-    neumann1 = 180 - neumann3_upper
-
-    xs = np.sqrt(circles[1].radius ** 2 - (y_horizontal - circles[1].center[1]) ** 2)
-    angle = np.rad2deg(np.arctan(-xs / np.sqrt(circles[1].radius ** 2 - xs ** 2))) % 180
-    angle = 180 - angle
-    neumann3_lower = angle
-    neumann2 = 180 - neumann3_lower
-
-    neumann3 = neumann3_upper + neumann3_lower
-    return neumann1, neumann2, neumann3
+
+    left_intersection = sorted(intersections, key=lambda e: e[0])[0]
+    lower_circle, upper_circle = sorted(circles, key=lambda c: c.center[1])
+
+    def rotate90_ccw(vector):
+        return np.array([-vector[1], vector[0]])
+
+    def rotate90_cw(vector):
+        return np.array([vector[1], -vector[0]])
+
+    def normalize(vector):
+        return vector / np.linalg.norm(vector)
+
+    vector_20 = normalize(rotate90_ccw(lower_circle.center - left_intersection))
+    vector_01 = np.array([-1, 0])
+    vector_12 = normalize(rotate90_cw(upper_circle.center - left_intersection))
+
+    angles = [np.rad2deg(np.arccos(np.dot(v1, v2))) for v1, v2 in [(vector_20, vector_01),
+                                                                   (vector_01, vector_12),
+                                                                   (vector_12, vector_20)]
+              ]
+    return angles
-- 
GitLab