Permalink
Browse files

replaced some old NR code by some of OOURA.

implemented a first attempt to calculate the Hankel transform of I(q) to obtain G(z)
  • Loading branch information...
Kohlbrecher committed Jun 24, 2016
1 parent 62843a4 commit b1d5c67a4566f5c08930f6c690aaeeb66c30e32a
View
@@ -31,3 +31,6 @@ SASfit is provided "as is" and with no warranty.
For license information see COPYING.txt.
A paper about SASfit has been published in
J. Appl. Cryst. (2015). 48, 1587-1598
doi:10.1107/S1600576715016544
View
@@ -4300,4 +4300,26 @@ @Article{Miltz1972
timestamp = {2016.05.11},
}
@Article{Wieder1999,
author = {Wieder, Thomas},
title = {Algorithm 794: Numerical Hankel Transform by the Fortran Program HANKEL},
year = {1999},
journal = {ACM Trans. Math. Softw.},
volume = {25},
number = {2},
month = jun,
pages = {240--250},
issn = {0098-3500},
doi = {10.1145/317275.317284},
url = {http://doi.acm.org/10.1145/317275.317284},
acmid = {317284},
address = {New York, NY, USA},
issue_date = {June 1999},
keywords = {Hankel transform, numerical analysis},
numpages = {11},
owner = {kohlbrecher},
publisher = {ACM},
timestamp = {2016.06.22},
}
@Comment{jabref-meta: databaseType:biblatex;}
@@ -36,5 +36,4 @@ if(NOT WIN32)
set(CPACK_GENERATOR "TBZ2")
else(NOT WIN32)
set(CPACK_GENERATOR "ZIP")
endif(NOT WIN32)
endif(NOT WIN32)
@@ -46,6 +46,7 @@ set(SOURCE_sasfit_common
sasfit_utils.c
sasfit_trapzd.c
sasfit_integrate.c
sasfit_intDE2.c
sasfit_mem.c
sasfit_message.c
sasfit_function.c
@@ -378,6 +378,7 @@ src/sasfit_common/sasfit_common.c
src/sasfit_common/sasfit_function.c
src/sasfit_common/sasfit_gamma.c
src/sasfit_common/sasfit_integrate.c
src/sasfit_common/sasfit_intDE2.c
src/sasfit_common/sasfit_mem.c
src/sasfit_common/sasfit_message.c
src/sasfit_common/sasfit_timer.c
@@ -172,6 +172,24 @@ typedef struct
double (*sasfit_pfq) (double *p_r, double *p_i, double *q_r, double *q_i, int ip, int iq,
double z_r, double z_i, int ln_pFq, int ix,
double *pFq_r, double *pFq_i, int nsigfig, sasfit_param * param);/* 95 */
void *reserved96; /* 96 */
void *reserved97; /* 97 */
void *reserved98; /* 98 */
void *reserved99; /* 99 */
void *reserved100; /* 100 */
void *reserved101; /* 101 */
void *reserved102; /* 102 */
void *reserved103; /* 103 */
void *reserved104; /* 104 */
void *reserved105; /* 105 */
void *reserved106; /* 106 */
void *reserved107; /* 107 */
void *reserved108; /* 108 */
void *reserved109; /* 109 */
void (* sasfit_intdeiini) (int lenaw, double tiny, double eps, double *aw); /* 110 */
void (* sasfit_intdei) (double (*f)(double, void *), double a, double *aw, double *i, double *err, void *fparams); /* 111 */
void (* sasfit_intdeoini) (int lenaw, double tiny, double eps, double *aw); /* 112 */
void (* sasfit_intdeo) (double (*f)(double, void *), double a, double omega, double *aw, double *i, double *err, void *fparams); /* 113 */
} sasfit_common_stubs_t;
#if defined(MAKE_SASFIT_PLUGIN)
@@ -516,6 +534,22 @@ typedef struct
#define sasfit_pfq \
(sasfit_common_stubs_ptr->sasfit_pfq) /* 94 */
#endif
#ifndef sasfit_intdeiini
#define sasfit_intdeiini \
(sasfit_common_stubs_ptr->sasfit_intdeiini) /* 110 */
#endif
#ifndef sasfit_intdei
#define sasfit_intdei \
(sasfit_common_stubs_ptr->sasfit_intdei) /* 111 */
#endif
#ifndef sasfit_intdeoini
#define sasfit_intdeoini \
(sasfit_common_stubs_ptr->sasfit_intdeoini) /* 112 */
#endif
#ifndef sasfit_intdeo
#define sasfit_intdeo \
(sasfit_common_stubs_ptr->sasfit_intdeo) /* 113 */
#endif
#endif /* defined(MAKE_SASFIT_PLUGIN) */
@@ -45,6 +45,7 @@
* \ingroup sasfit_eps
* Contains various numerical constants.
*/
typedef struct
{
scalar aniso;
@@ -115,6 +115,10 @@ sasfit_common_DLLEXP scalar sasfit_rod_fc(scalar Q, scalar R);
sasfit_common_DLLEXP scalar sasfit_sphere_fc(scalar Q, scalar R);
sasfit_common_DLLEXP scalar sasfit_gauss_fc(scalar Q, scalar R);
sasfit_common_DLLEXP scalar sasfit_g(scalar fp, scalar a);
sasfit_common_DLLEXP void sasfit_intdeo(double (*f)(double, void *), double a, double omega, double *aw, double *i, double *err, void *fparams);
sasfit_common_DLLEXP void sasfit_intdeoini(int lenaw, double tiny, double eps, double *aw);
sasfit_common_DLLEXP void sasfit_intdei(double (*f)(double, void *), double a, double *aw, double *i, double *err, void *fparams);
sasfit_common_DLLEXP void sasfit_intdeiini(int lenaw, double tiny, double eps, double *aw);
/**
* Returns the shared library prefix for the current platform.
@@ -137,6 +137,24 @@ sasfit_common_stubs_t sasfit_common_stubs = {
sasfit_3f2, /* 93 */
sasfit_2f1, /* 94 */
sasfit_pfq, /* 95 */
NULL, /* 96 */
NULL, /* 97 */
NULL, /* 98 */
NULL, /* 99 */
NULL, /* 100 */
NULL, /* 101 */
NULL, /* 102 */
NULL, /* 103 */
NULL, /* 104 */
NULL, /* 105 */
NULL, /* 106 */
NULL, /* 107 */
NULL, /* 108 */
NULL, /* 109 */
sasfit_intdeiini, /* 110 */
sasfit_intdei, /* 111 */
sasfit_intdeoini, /* 112 */
sasfit_intdeo, /* 113 */
};
/* !END!: Do not edit above this line, see sasfit_common.decls for modifications. */
@@ -1,3 +1,18 @@
/*
* The code is an adaption of a numerical package from Takuya OOURA
* (email: ooura@mmm.t.u-tokyo.ac.jp)
* under http://www.kurims.kyoto-u.ac.jp/~ooura/ he supplies several numerical packages in Fortran and C.
* The code below has been taken from his libarary package for numerical integration (Quadrature) -
* DE Formula (Almighty Quadrature).
*
* see also the References:
* [1] M.Mori, Developments in the double exponential formula for numerical integration, Proceedings of the International Congress of Mathematicians, Kyoto 1990, 1991, Springer-Verlag, 1585-1594.
* [2] H.Takahasi and M.Mori, Double exponential formulas for numerical integration, Pub. RIMS Kyoto Univ. 9, 1974, 721-741
* [3] T.Ooura and M.Mori, The Double exponential formula for oscillatory functions over the half infinite interval, Journal of Computational and Applied Mathematics 38, 1991, 353-360
* [4] M.Mori and T.Ooura, Double exponential formulas for Fourier type integrals with a divergent integrand, Contributions in Numerical Mathematics, ed. R.P.Agarwal, World Scientific Series in Applicable Analysis, 2, 1993, 301-308
* [5] H.Toda and H.Ono, Some remarks for efficient usage of the double exponential formulas(in Japanese), Kokyuroku, RIMS, Kyoto Univ. 339, 1978, 74-109.
*/
/*
DE-Quadrature
Numerical Automatic Integrator for Improper Integral
@@ -176,7 +191,31 @@ intdeo
#include <math.h>
#include "include/sasfit_common.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#ifdef MACOSX
#include <sys/malloc.h>
#else
#include <malloc.h>
#endif
#include <sasfit_sd.h>
#include <sasfit_sq.h>
#include <tcl.h>
#include <sasfit.h>
#include <SASFIT_nr.h>
#include <SASFIT_x_tcl.h>
#include <SASFIT_resolution.h>
#include <tcl_cmds.h>
#include <sasfit_common.h>
#include <sasfit_plugin.h>
#include <sasfit_core.h>
#include <sasfit_plugin_backend.h>
void intdeini(int lenaw, double tiny, double eps, double *aw)
{
@@ -572,12 +611,39 @@ void sasfit_intdeoini(int lenaw, double tiny, double eps, double *aw)
void sasfit_intdeo(double (*f)(double, void *), double a, double omega, double *aw,
double *i, double *err, void *fparams)
double *i, double *err, void *GIP)
{
int lenawm, nk0, noff0, nk, noff, lmax, m, k, j, jm, l;
double eps, per, perw, w02, ir, h, iback, irback, t, tk,
xa, fm, fp, errh, s0, s1, s2, errd;
Tcl_Interp *interp;
scalar z;
scalar *par;
scalar *Ifit;
scalar *Isub;
scalar *dydpar;
scalar *Intdydpar, *Intdydpar_r;
scalar IntIsub, IntIsub_r, IntIsub_err, IntIsub_back,IntIsub_errh, IntIsub_errd, IntIsub_rback, fmIsub, fpIsub, s2Isub, s0Isub, s1Isub;
int max_SD;
sasfit_analytpar *AP;
int error_type;
bool *error;
interp = (( sasfit_GzIntStruct *) GIP)->interp;
z = (( sasfit_GzIntStruct *) GIP)->z;
par = (( sasfit_GzIntStruct *) GIP)->par;
Ifit = (( sasfit_GzIntStruct *) GIP)->Ifit;
Isub = (( sasfit_GzIntStruct *) GIP)->Isub;
dydpar = (( sasfit_GzIntStruct *) GIP)->dydpar;
max_SD = (( sasfit_GzIntStruct *) GIP)->max_SD;
AP = (( sasfit_GzIntStruct *) GIP)->AP;
error_type = ((sasfit_GzIntStruct *) GIP)->error_type;
error = (( sasfit_GzIntStruct *) GIP)->error;
// Intdydpar = (scalar *)malloc((max_SD*(3*MAXPAR))*sizeof(scalar));
// Intdydpar_r = (scalar *)malloc((max_SD*(3*MAXPAR))*sizeof(scalar));
lenawm = (int) (aw[0] + 0.5);
nk0 = (int) (aw[1] + 0.5);
noff0 = 6;
@@ -588,16 +654,28 @@ void sasfit_intdeo(double (*f)(double, void *), double a, double omega, double *
per = 1 / fabs(omega);
w02 = 2 * aw[noff + 2];
perw = per * w02;
*i = (*f)(a + aw[noff] * per,fparams);
*i = (*f)(a + aw[noff] * per,GIP);
IntIsub = *Isub;
ir = *i * aw[noff + 1];
IntIsub_r = IntIsub* aw[noff + 1];
*i *= aw[noff + 2];
IntIsub *= aw[noff + 2];
*err = fabs(*i);
IntIsub_err = fabs(IntIsub);
h = 2;
m = 1;
k = noff;
do {
iback = *i;
irback = ir;
IntIsub_back = IntIsub;
IntIsub_rback = IntIsub_r;
t = h * 0.5;
do {
if (k == noff) {
@@ -607,78 +685,118 @@ void sasfit_intdeo(double (*f)(double, void *), double a, double omega, double *
do {
j += 3;
xa = per * aw[j];
fm = (*f)(a + xa,fparams);
fp = (*f)(a + xa + perw * tk,fparams);
fm = (*f)(a + xa,GIP);
fmIsub = *Isub;
fp = (*f)(a + xa + perw * tk,GIP);
fpIsub = *Isub;
ir += (fm + fp) * aw[j + 1];
IntIsub_r += (fmIsub + fpIsub) * aw[j + 1];
fm *= aw[j + 2];
fmIsub *= aw[j + 2];
fp *= w02 - aw[j + 2];
fpIsub *= w02 - aw[j + 2];
*i += fm + fp;
IntIsub += fm + fp;
*err += fabs(fm) + fabs(fp);
IntIsub_err += fabs(fm) + fabs(fp);
tk += 1;
} while (aw[j] > eps && j < k);
errh = *err * aw[5];
*err *= eps;
IntIsub_errh = IntIsub * aw[5];
IntIsub_err *=eps;
jm = j - noff;
} else {
tk = t;
for (j = k + 3; j <= k + jm; j += 3) {
xa = per * aw[j];
fm = (*f)(a + xa,fparams);
fp = (*f)(a + xa + perw * tk,fparams);
fm = (*f)(a + xa,GIP);
fmIsub = *Isub;
fp = (*f)(a + xa + perw * tk,GIP);
fpIsub = *Isub;
ir += (fm + fp) * aw[j + 1];
IntIsub_r += (fmIsub + fpIsub) * aw[j + 1];
fm *= aw[j + 2];
fmIsub *= aw[j + 2];
fp *= w02 - aw[j + 2];
fpIsub *= w02 - aw[j + 2];
*i += fm + fp;
IntIsub += fm + fp;
tk += 1;
}
j = k + jm;
k += nk;
}
while (fabs(fm) > *err && j < k) {
j += 3;
fm = (*f)(a + per * aw[j],fparams);
fm = (*f)(a + per * aw[j],GIP);
fmIsub = *Isub;
ir += fm * aw[j + 1];
IntIsub_r += fmIsub * aw[j + 1];
fm *= aw[j + 2];
fmIsub *= aw[j + 2];
*i += fm;
IntIsub += fm;
}
fm = (*f)(a + perw * tk,fparams);
fm = (*f)(a + perw * tk,GIP);
fmIsub = *Isub;
s2 = w02 * fm;
s2Isub = w02 * fmIsub;
*i += s2;
IntIsub += s2Isub;
if (fabs(fp) > *err || fabs(s2) > *err) {
l = 0;
for (;;) {
l++;
s0 = 0;
s1 = 0;
s2 = fm * aw[noff0 + 1];
s0Isub=0;
s1Isub=0;
s2Isub = fmIsub * aw[noff0 + 1];
for (j = noff0 + 2; j <= noff - 2; j += 2) {
tk += 1;
fm = (*f)(a + perw * tk,fparams);
fm = (*f)(a + perw * tk,GIP);
fmIsub = *Isub;
s0 += fm;
s1 += fm * aw[j];
s2 += fm * aw[j + 1];
s0Isub += fmIsub;
s1Isub += fmIsub * aw[j];
s2Isub += fmIsub * aw[j + 1];
}
if (s2 <= *err || l >= lmax) break;
*i += w02 * s0;
IntIsub += w02 * s0Isub;
}
*i += s1;
IntIsub += s1Isub;
if (s2 > *err) *err = s2;
if (s2Isub > IntIsub_err) IntIsub_err = s2Isub;
}
t += h;
} while (t < 1);
if (m == 1) {
errd = 1 + 2 * errh;
IntIsub_errd = 1 + 2 * IntIsub_errh;
} else {
errd = h * (fabs(*i - 2 * iback) + fabs(ir - 2 * irback));
IntIsub_errd = h * (fabs(*i - 2 * IntIsub_back) + fabs(IntIsub_r - 2 * IntIsub_rback));
}
h *= 0.5;
m *= 2;
} while (errd > errh && 2 * k - noff <= lenawm);
*i *= h * per;
IntIsub *= h * per;
if (errd > errh) {
*err = -errd * per;
} else {
*err *= per * m * 0.5;
};
if (IntIsub_errd > IntIsub_errh) {
IntIsub_err = -IntIsub_errd * per;
} else {
IntIsub_err *= per * m * 0.5;
}
}
Oops, something went wrong.

0 comments on commit b1d5c67

Please sign in to comment.