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 #include #include #include #include +#include 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));