From f940e62e5272a1df5dcb39a79181146fdcc9b621 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Mon, 30 Sep 2019 08:53:57 +0300
Subject: [PATCH] Test and fix the rounded half-ellipse contour..

Former-commit-id: 448c30d375b3f9d0abab9442ae00bfcd123e5cb9
---
 qpms/beyn.c           | 20 ++++++++++----------
 qpms/beyn.h           |  6 ++++++
 tests/CMakeLists.txt  |  4 ++++
 tests/kidneycontour.c | 19 +++++++++++++++++++
 4 files changed, 39 insertions(+), 10 deletions(-)
 create mode 100644 tests/kidneycontour.c

diff --git a/qpms/beyn.c b/qpms/beyn.c
index 9f6b081..28911be 100644
--- a/qpms/beyn.c
+++ b/qpms/beyn.c
@@ -216,24 +216,24 @@ beyn_contour_t *beyn_contour_kidney(complex double centre, double rRe,
   }
   // First small arc
   t_lo = t_hi; t_hi = t_lo + (M_PI_2 + alpha) * y;
-  for(; t = i * dt, dt < t_hi; ++i) {
+  for(; t = i * dt, t < t_hi; ++i) {
     double phi = (t - t_lo) / y - M_PI_2;
     c->z_dz[i][0] = centre + (ar * (x + y * cos(phi)) + y * (1 + sin(phi)) * I) * faktor;
-    c->z_dz[i][1] = dt * y * (- ar * sin(phi) + cos(phi) * I) * faktor;
+    c->z_dz[i][1] = dt * (- ar * sin(phi) + cos(phi) * I) * faktor;
   }
   // Big arc
-  t_lo = t_hi; t_hi = t_lo + (M_PI + 2 * alpha) * h;
-  for(; t = i * dt, dt < t_hi; ++i) {
-    double phi = (t - t_lo) / h + M_PI_2 + alpha;
-    c->z_dz[i][0] = centre + (ar * (h * cos(phi)) + h * (1 + sin(phi)) * I) * faktor;
-    c->z_dz[i][1] = dt * h * (- ar * sin(phi) + cos(phi) * I) * faktor;
+  t_lo = t_hi; t_hi = t_lo + (M_PI - 2 * alpha) * h;
+  for(; t = i * dt, t < t_hi; ++i) {
+    double phi = (t - t_lo) / h + alpha;
+    c->z_dz[i][0] = centre + (ar * (h * cos(phi)) + h * sin(phi) * I) * faktor;
+    c->z_dz[i][1] = dt * (- ar * sin(phi) + cos(phi) * I) * faktor;
   }
   // Second small arc
   t_lo = t_hi; t_hi = t_lo + (M_PI_2 + alpha) * y;
-  for(; t = i * dt, dt < t_hi; ++i) {
+  for(; t = i * dt, t < t_hi; ++i) {
     double phi = (t - t_lo) / y + M_PI - alpha;
-    c->z_dz[i][0] = centre + (ar * (x + y * cos(phi)) + y * (1 + sin(phi)) * I) * faktor;
-    c->z_dz[i][1] = dt * y * (- ar * sin(phi) + cos(phi) * I) * faktor;
+    c->z_dz[i][0] = centre + (ar * (- x + y * cos(phi)) + y * (1 + sin(phi)) * I) * faktor;
+    c->z_dz[i][1] = dt * (- ar * sin(phi) + cos(phi) * I) * faktor;
   }
   // Straight line, second part
   t_lo = t_hi; t_hi = tmax;
diff --git a/qpms/beyn.h b/qpms/beyn.h
index c822ecd..25bc1c7 100644
--- a/qpms/beyn.h
+++ b/qpms/beyn.h
@@ -60,6 +60,12 @@ typedef enum {
 beyn_contour_t *beyn_contour_halfellipse(complex double centre, double halfax_re, double halfax_im, size_t npoints,
 		beyn_contour_halfellipse_orientation or);
 
+/// Similar to halfellipse but with rounded corners.
+beyn_contour_t *beyn_contour_kidney(complex double centre, double halfax_re, double halfax_im, 
+	double rounding, ///< Must be in interval [0, 0.5)
+	size_t n, beyn_contour_halfellipse_orientation or);
+
+
 /// Beyn algorithm result structure (GSL matrix/vector version).
 typedef struct beyn_result_gsl_t {
 	size_t neig; ///< Number of eigenvalues found (a bit redundant?).
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 5d39526..a8513dc 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -15,6 +15,10 @@ add_executable( test_scalar_ewald32 test_scalar_ewald32.c )
 target_link_libraries( test_scalar_ewald32 qpms gsl lapacke amos m )
 target_include_directories( test_scalar_ewald32 PRIVATE .. )
 
+add_executable( kidneycontour kidneycontour.c )
+target_link_libraries( kidneycontour qpms gsl lapacke amos m )
+target_include_directories( kidneycontour PRIVATE .. )
+
 add_executable( tbeyn tbeyn.c )
 target_link_libraries( tbeyn qpms gsl lapacke amos m )
 target_include_directories( tbeyn PRIVATE .. )
diff --git a/tests/kidneycontour.c b/tests/kidneycontour.c
new file mode 100644
index 0000000..8f4e635
--- /dev/null
+++ b/tests/kidneycontour.c
@@ -0,0 +1,19 @@
+#include <qpms/beyn.h>
+#include <stdio.h>
+
+#define CPAIR(x) creal(x), cimag(x)
+
+int main(int argc, char **argv)
+{
+  double rRe = 2e3, rIm = 1.5e3, rounding = 0.2;
+  complex double centre = 1e3 * I;
+  size_t n = 100;
+
+  beyn_contour_t *c = beyn_contour_kidney(centre, rRe, rIm, rounding, n, BEYN_CONTOUR_HALFELLIPSE_IM_PLUS);
+
+  for(size_t i = 0; i < n; ++i) 
+    printf("%g\t%g\t%g\t%g\n", CPAIR(c->z_dz[i][0]), CPAIR(c->z_dz[i][1]));
+
+  free(c);
+  return 0;
+}