diff --git a/qpms/beyn.c b/qpms/beyn.c index 1623666..b33291c 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -56,10 +56,19 @@ static complex double zrandN(double sigma, double mu){ return randN(sigma, mu) + I*randN(sigma, mu); } +static inline double dsq(double a) { return a * a; } + +static _Bool beyn_contour_ellipse_inside_test(struct beyn_contour_t *c, complex double z) { + double rRe = c->z_dz[c->n][0]; + double rIm = c->z_dz[c->n][1]; + complex double zn = z - c->centre; + return dsq(creal(zn)/rRe) + dsq(cimag(zn)/rIm) < 1; +} + beyn_contour_t *beyn_contour_ellipse(complex double centre, double rRe, double rIm, size_t n) { beyn_contour_t *c; - QPMS_CRASHING_MALLOC(c, sizeof(beyn_contour_t) + n*sizeof(c->z_dz[0])); + QPMS_CRASHING_MALLOC(c, sizeof(beyn_contour_t) + (n+1)*sizeof(c->z_dz[0])); c->centre = centre; c->n = n; for(size_t i = 0; i < n; ++i) { @@ -68,6 +77,10 @@ beyn_contour_t *beyn_contour_ellipse(complex double centre, double rRe, double r c->z_dz[i][0] = centre + ct*rRe + I*st*rIm; c->z_dz[i][1] = (I*rRe*st + rIm*ct) / n; } + // We hide the half-axes metadata after the discretisation points. + c->z_dz[n][0] = rRe; + c->z_dz[n][1] = rIm; + c->inside_test = beyn_contour_ellipse_inside_test; return c; } diff --git a/qpms/beyn.h b/qpms/beyn.h index 4cf9213..5eeccf9 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -37,6 +37,8 @@ typedef struct beyn_contour_t { * but it should be "somewhere around" where the contour is. */ complex double centre; + /// Function testing that a point \a z lies inside the contour (optional). + _Bool (*inside_test)(struct beyn_contour_t *, complex double z); complex double z_dz[][2]; ///< Pairs of contour points and derivatives in that points. } beyn_contour_t; diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 7ea32e9..1e7712e 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -549,6 +549,7 @@ cdef extern from "gsl/gsl_matrix.h": cdef extern from "beyn.h": ctypedef struct beyn_contour_t: + bint (*inside_test)(beyn_contour_t *, cdouble z) pass ctypedef struct beyn_result_gsl_t: pass