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