From 3312bc61e2e95863d41935e4a6fad6f455f56a10 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 26 Jan 2018 15:07:29 +0200 Subject: [PATCH] =?UTF-8?q?Pr=C3=A1ce=20na=20vswf=20;=20jdu=20dom=C5=AF?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Former-commit-id: 112c65dcb0de85908b43cd488c63ac14072dc97e --- notes/ewald.lyx | 74 +++++++++++++++++++++++++++++++++++++++++++---- qpms/qpms_types.h | 2 ++ qpms/vswf.c | 18 ++++++------ qpms/vswf.h | 14 ++++----- 4 files changed, 86 insertions(+), 22 deletions(-) diff --git a/notes/ewald.lyx b/notes/ewald.lyx index 1a97fc1..4ee8591 100644 --- a/notes/ewald.lyx +++ b/notes/ewald.lyx @@ -915,7 +915,69 @@ reference "eq:2d long range regularisation problem statement" \end_inset might lead to satisfactory results; see below. - +\end_layout + +\begin_layout Standard +\begin_inset Note Note +status open + +\begin_layout Plain Layout +In natural/dimensionless units; +\begin_inset Formula $x=k_{0}r$ +\end_inset + +, +\begin_inset Formula $\tilde{k}=k/k_{0}$ +\end_inset + +, +\begin_inset Formula $č=c/k_{0}$ +\end_inset + + +\begin_inset Formula +\[ +s_{q}(x)\equiv e^{ix}x^{-q} +\] + +\end_inset + + +\begin_inset Formula +\[ +\tilde{\rho}_{\kappa,č}(x)\equiv\left(1-e^{-čx}\right)^{\text{\kappa}}=\left(1-e^{-\frac{c}{k_{0}}k_{0}r}\right)^{\kappa}=\left(1-e^{-cr}\right)^{\kappa}=\rho_{\kappa,c}(r) +\] + +\end_inset + + +\begin_inset Formula +\[ +s_{q}^{\textup{L}}\left(x\right)\equiv s_{q}(x)\tilde{\rho}_{\kappa,č}(x)=e^{ix}x^{-q}\left(1-e^{-čx}\right)^{\text{\kappa}} +\] + +\end_inset + + +\begin_inset Formula +\begin{eqnarray*} +\pht n{s_{q}^{\textup{L}}}\left(\tilde{k}\right) & = & \int_{0}^{\infty}s_{q}^{\textup{L}}\left(x\right)xJ_{n}\left(\tilde{k}x\right)\ud x=\int_{0}^{\infty}s_{q}\left(x\right)\tilde{\rho}_{\kappa,č}(x)xJ_{n}\left(\tilde{k}x\right)\ud x\\ + & = & \int_{0}^{\infty}s_{k_{0},q}\left(r\right)\rho_{\kappa,c}(r)\left(k_{0}r\right)J_{n}\left(kr\right)\ud\left(k_{0}r\right)\\ + & = & k_{0}^{2}\int_{0}^{\infty}s_{k_{0},q}\left(r\right)\rho_{\kappa,c}(r)rJ_{n}\left(kr\right)\ud r\\ + & = & k_{0}^{2}\pht n{s_{k_{0},q}^{\textup{L}}}\left(k\right) +\end{eqnarray*} + +\end_inset + + +\end_layout + +\end_inset + + +\end_layout + +\begin_layout Standard \begin_inset Note Note status open @@ -3037,7 +3099,7 @@ where the spherical Hankel transform 2) \begin_inset Formula \[ -\bsht lg(k)\equiv\int_{0}^{\infty}\ud r\, r^{2}g(r)j_{l}\left(kr\right). +\bsht lg(k)\equiv\int_{0}^{\infty}\ud r\,r^{2}g(r)j_{l}\left(kr\right). \] \end_inset @@ -3047,7 +3109,7 @@ Using this convention, the inverse spherical Hankel transform is given by 3) \begin_inset Formula \[ -g(r)=\frac{2}{\pi}\int_{0}^{\infty}\ud k\, k^{2}\bsht lg(k)j_{l}(k), +g(r)=\frac{2}{\pi}\int_{0}^{\infty}\ud k\,k^{2}\bsht lg(k)j_{l}(k), \] \end_inset @@ -3060,7 +3122,7 @@ so it is not unitary. An unitary convention would look like this: \begin_inset Formula \begin{equation} -\usht lg(k)\equiv\sqrt{\frac{2}{\pi}}\int_{0}^{\infty}\ud r\, r^{2}g(r)j_{l}\left(kr\right).\label{eq:unitary 3d Hankel tf definition} +\usht lg(k)\equiv\sqrt{\frac{2}{\pi}}\int_{0}^{\infty}\ud r\,r^{2}g(r)j_{l}\left(kr\right).\label{eq:unitary 3d Hankel tf definition} \end{equation} \end_inset @@ -3114,8 +3176,8 @@ where the Hankel transform of order is defined as \begin_inset Formula \begin{eqnarray} -\pht mg\left(k\right) & = & \int_{0}^{\infty}\ud r\, g(r)J_{m}(kr)r\label{eq:unitary 2d Hankel tf definition}\\ - & = & \left(-1\right)^{m}\int_{0}^{\infty}\ud r\, g(r)J_{\left|m\right|}(kr)r +\pht mg\left(k\right) & = & \int_{0}^{\infty}\ud r\,g(r)J_{m}(kr)r\label{eq:unitary 2d Hankel tf definition}\\ + & = & \left(-1\right)^{m}\int_{0}^{\infty}\ud r\,g(r)J_{\left|m\right|}(kr)r \end{eqnarray} \end_inset diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index 22aff7e..b59375f 100644 --- a/qpms/qpms_types.h +++ b/qpms/qpms_types.h @@ -18,7 +18,9 @@ typedef enum { // Normalisations typedef enum { + // As in TODO QPMS_NORMALIZATION_XU = 3, // NI! + // As in http://www.eit.lth.se/fileadmin/eit/courses/eit080f/Literature/book.pdf, power-normalised QPMS_NORMALIZATION_KRISTENSSON = 2, // NI! QPMS_NORMALIZATION_POWER = QPMS_NORMALIZATION_KRISTENSSON, // NI! QPMS_NORMALIZATION_TAYLOR = 1, diff --git a/qpms/vswf.c b/qpms/vswf.c index 5195dee..661cda6 100644 --- a/qpms/vswf.c +++ b/qpms/vswf.c @@ -2,7 +2,7 @@ #include // Legendre functions also for negative m, see DLMF 14.9.3 -qpms_errno_t qpms_legendre_deriv_y_fill(double *target, double *target_deriv, double x, qpms_l_t lMax, +int qpms_errno_t qpms_legendre_deriv_y_fill(double *target, double *target_deriv, double x, qpms_l_t lMax, gsl_sf_legendre_t lnorm, double csphase) { size_t n = gsl_sf_legenre_array_n(lMax); @@ -12,7 +12,7 @@ qpms_errno_t qpms_legendre_deriv_y_fill(double *target, double *target_deriv, do lnorm, (size_t)lMax, x, csphase, legendre_tmp,legendre_tmp_deriv); for (qpms_l_t l = 0; l <= lMax; ++l) for (qpms_m_t m = 0; m <= l; ++m) { - qpms_y_t y = gpms_mn2y(m,l); + qpms_y_t y = qpms_mn2y(m,l); size_t i = gsl_sf_legenre_array_index(l,m); target[y] = legendre_tmp[i]; target_deriv[y] = legendre_deriv_tmp[i]; @@ -21,7 +21,7 @@ qpms_errno_t qpms_legendre_deriv_y_fill(double *target, double *target_deriv, do case GSL_SF_LEGEDRE_NONE: for (qpms_l_t l = 0; l <= lMax; ++l) for (qpms_m_t m = 1; m <= l; ++m) { - qpms_y_t y = gpms_mn2y(-m,l); + qpms_y_t y = qpms_mn2y(-m,l); size_t i = gsl_sf_legenre_array_index(l,m); // viz DLMF 14.9.3, čert ví, jak je to s cs fasí. double factor = exp(lgamma(l-m+1)-lgamma(n+m+1))*((m%2)?-1:1); @@ -34,7 +34,7 @@ qpms_errno_t qpms_legendre_deriv_y_fill(double *target, double *target_deriv, do case GSL_SF_LEGENDRE_FULL: for (qpms_l_t l = 0; l <= lMax; ++l) for (qpms_m_t m = 1; m <= l; ++m) { - qpms_y_t y = gpms_mn2y(-m,l); + qpms_y_t y = qpms_mn2y(-m,l); size_t i = gsl_sf_legenre_array_index(l,m); // viz DLMF 14.9.3, čert ví, jak je to s cs fasí. double factor = ((m%2)?-1:1); // this is the difference from the unnormalised case @@ -51,13 +51,13 @@ qpms_errno_t qpms_legendre_deriv_y_fill(double *target, double *target_deriv, do return QPMS_SUCCESS; } -// FIXME ZAPOMNĚL JSEM NA POLE DERIVACÍ (TÉŽ HLAVIČKOVÝ SOUBOR) -double *qpms_legendre_deriv_y_get(double x, qpms_l_t lMax, gsle_sf_legendre_t lnorm, +int qpms_legendre_deriv_y_get(double **target, double **dtarget, double x, qpms_l_t lMax, gsl_sf_legendre_t lnorm, double csphase) { - double *ar = malloc(sizeof(double)*qpms_lMaxnelem(lMax)); - if (qpms_legendre_deriv_y_fill(ar, x, lMax, lnorm, csphase)) - abort(); + + *target = malloc(sizeof(double)*qpms_lMax2nelem(lMax)); + *dtarget = malloc(sizeof(double)*qpms_lMax2nelem(lMax)); + return qpms_legendre_deriv_y_fill(ar, x, lMax, lnorm, csphase); } diff --git a/qpms/vswf.h b/qpms/vswf.h index 8f989bc..4903a86 100644 --- a/qpms/vswf.h +++ b/qpms/vswf.h @@ -4,10 +4,10 @@ #include // Electric wave N; NI -complex double qpms_vswf_single_el(int m, int n, sph_t kdlj, +csphvec_t qpms_vswf_single_el(int m, int n, sph_t kdlj, qpms_bessel_t btyp, qpms_normalization_t norm); // Magnetic wave M; NI -complex double qpms_vswf_single_mg(int m, int n, sph_t kdlj, +csphvec_t qpms_vswf_single_mg(int m, int n, sph_t kdlj, qpms_bessel_t btyp, qpms_normalization_t norm); // Set of electric and magnetic VSWF in spherical coordinate basis @@ -24,12 +24,12 @@ typedef struct { * N.B. for the norm definitions, see * https://www.gnu.org/software/gsl/manual/html_node/Associated-Legendre-Polynomials-and-Spherical-Harmonics.html * ( gsl/specfunc/legendre_source.c and 7.24.2 of gsl docs - * + */ -double *qpms_legendre_deriv_y_get(double x, qpms_l_t lMax, - gsle_sf_legendre_t lnorm, double csphase); //FIXME what about derivative? -qpms_errno_t qpms_legendre_deriv_y_fill(double *where, double *where_deriv, double x, - qpms_l_t lMax, gsl_sf_legendre_t lnorm, double csphase); //NI +int qpms_legendre_deriv_y_get(double **result, double **result_deriv, double x, qpms_l_t lMax, + gsle_sf_legendre_t lnorm, double csphase); +int qpms_errno_t qpms_legendre_deriv_y_fill(double *where, double *where_deriv, double x, + qpms_l_t lMax, gsl_sf_legendre_t lnorm, double csphase); qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj, qpms_bessel_t btyp, qpms_normalization_t norm);//NI