From 2e7b925be2fc53ae7debee2f6a6fd825ca549389 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Thu, 26 Sep 2019 11:36:41 +0300
Subject: [PATCH] WIP trying to fix the half-ellipse contour.

Former-commit-id: cbba6fffacbb38322a590478899f38c7d7dafe8e
---
 qpms/beyn.c          |  4 ++--
 tests/CMakeLists.txt |  5 +++++
 tests/tbeyn3.c       | 11 +++++++++--
 3 files changed, 16 insertions(+), 4 deletions(-)

diff --git a/qpms/beyn.c b/qpms/beyn.c
index bbc9f2b..9bb0974 100644
--- a/qpms/beyn.c
+++ b/qpms/beyn.c
@@ -125,7 +125,7 @@ beyn_contour_t *beyn_contour_halfellipse(complex double centre, double rRe,
   }
   for(size_t i = 0; i < nline; ++i) {
     double t = 0.5 * (1 - (double) nline) + i;
-    complex double z = centre + faktor * t * l;
+    complex double z = centre + faktor * t * 2 * l / nline;
     switch(or) { // ensure correct zero signs; CHECKME!!!
       case BEYN_CONTOUR_HALFELLIPSE_RE_PLUS:
         if(creal(z) == 0 && signbit(creal(z)))
@@ -146,7 +146,7 @@ beyn_contour_t *beyn_contour_halfellipse(complex double centre, double rRe,
       default: QPMS_WTF;
     }
     c->z_dz[narc + i][0] = z;
-    c->z_dz[narc + i][1] = faktor * l / nline;
+    c->z_dz[narc + i][1] = faktor * 2 * l / nline;
   }
   // We hide the half-axes metadata after the discretisation points.
   c->z_dz[n][0] = rRe;
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 8185d79..5d39526 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -28,6 +28,11 @@ target_link_libraries( tbeyn3a qpms gsl lapacke amos m )
 target_include_directories( tbeyn3a PRIVATE .. )
 target_compile_definitions( tbeyn3a PRIVATE VARIANTA )
 
+add_executable( tbeyn3a_implus tbeyn3.c )
+target_link_libraries( tbeyn3a_implus qpms gsl lapacke amos m )
+target_include_directories( tbeyn3a_implus PRIVATE .. )
+target_compile_definitions( tbeyn3a_implus PRIVATE VARIANTA IMPLUS )
+
 add_executable( tbeyn3b tbeyn3.c )
 target_link_libraries( tbeyn3b qpms gsl lapacke amos m )
 target_include_directories( tbeyn3b PRIVATE .. )
diff --git a/tests/tbeyn3.c b/tests/tbeyn3.c
index 0c983b6..fd0260e 100644
--- a/tests/tbeyn3.c
+++ b/tests/tbeyn3.c
@@ -1,8 +1,10 @@
+#define _GNU_SOURCE
 #include <qpms/beyn.h>
 #include <stdio.h>
 #include <string.h>
 #include <stdlib.h>
 #include <math.h>
+#include <fenv.h>
 
 static double randU(double a, double b) {return a + (b-a) * random() * (1. / RAND_MAX); }
 // Normal distribution via Box-Muller transform
@@ -45,13 +47,14 @@ int M_function(complex double *target, const size_t m, const complex double z, v
 }
 
 int main(int argc, char **argv) {
-  complex double z0 = 0;
+  feenableexcept(FE_INVALID | FE_OVERFLOW);
+  complex double z0 = 0+3e-1*I;
 #ifdef RXSMALL
   double Rx = .1;
 #else 
   double Rx = .3; // Variant B will fail in this case due to large number of eigenvalues (>30)
 #endif
-  double Ry = .3;
+  double Ry = .25;
 #ifdef VARIANTF
   int L = 10, N = 150, dim = 10;
 #else
@@ -59,7 +62,11 @@ int main(int argc, char **argv) {
 #endif
   if (argc > 1) N = atoi(argv[1]);
   if (argc > 2) L = atoi(argv[2]);
+#ifdef IMPLUS
+  beyn_contour_t *contour = beyn_contour_halfellipse(z0, Rx, Ry, N, BEYN_CONTOUR_HALFELLIPSE_IM_PLUS);
+#else
   beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N);
+#endif
   struct param p;
   p.T0 = malloc(dim*dim*sizeof(double));
   p.T1 = malloc(dim*dim*sizeof(double));