QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Test and fix the rounded half-ellipse contour.. 

Former-commit-id: 448c30d375b3f9d0abab9442ae00bfcd123e5cb9
This commit is contained in:
Marek Nečada 2019-09-30 08:53:57 +03:00
parent d6084f3264
commit f940e62e52
4 changed files with 39 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
qpms/beyn.c
View File

@ -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;

6
qpms/beyn.h
View File

@ -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?).

4
tests/CMakeLists.txt
View File

@ -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 .. )

19
tests/kidneycontour.c Normal file
View File

@ -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;
}