diff --git a/explore/precision.py b/explore/precision.py index fdc009b7..05880f76 100755 --- a/explore/precision.py +++ b/explore/precision.py @@ -803,6 +803,21 @@ def np_2J1x_x(x): ocl_function=make_ocl("return sas_2J1x_x(q);", "sas_2J1x_x", ["lib/polevl.c", "lib/sas_J1.c"]), ) +add_function( + name="sin_arccos", + mp_function=lambda x: (mp.sin(mp.acos(x))), + np_function=lambda x: (np.sin(np.arccos(x))), + ocl_function=make_ocl("return sin(acos(q));", "sas_sin_arccos"), + limits=(0, 1), +) +add_function( + name="sin_from_cos", + mp_function=lambda x: (mp.sqrt(1 - x*x)), + np_function=lambda x: (np.sqrt(1-x*x)), + ocl_function=make_ocl("return sqrt(1.-q*q);", "sas_sin_from_cos"), + limits=(0, 1), +) + ALL_FUNCTIONS = set(FUNCTIONS.keys()) ALL_FUNCTIONS.discard("loggamma") # use cephes-based gammaln instead ALL_FUNCTIONS.discard("3j1/x:taylor") diff --git a/sasmodels/compare.py b/sasmodels/compare.py index 26b67ada..1c5c5ca5 100755 --- a/sasmodels/compare.py +++ b/sasmodels/compare.py @@ -75,7 +75,7 @@ def __call__(self, **par: float) -> np.ndarray: ... === model parameters === -preset*/-random[=seed] preset or random parameters - -sets=n generates n random datasets with the seed given by -random=seed + -sets=n generates n random datasets using -seed=seed -pars/-nopars* prints the parameter set or not -sphere[=150] set up spherical integration over theta/phi using n points -mono*/-poly suppress or allow polydispersity on generated parameters @@ -1049,7 +1049,7 @@ def plot_models(opts, result, limits=None, setnum=0): '2d', '1d', 'sesans', # Parameter set - 'preset', 'random', 'random=', 'sets=', + 'preset', 'random', 'random=', 'sets=', 'seed=', 'nopars', 'pars', 'sphere', 'sphere=', # integrate over a sphere in 2d with n points 'poly', 'mono', @@ -1253,6 +1253,7 @@ def parse_opts(argv): elif arg.startswith('-res='): opts['res'] = arg[5:] elif arg.startswith('-noise='): opts['noise'] = float(arg[7:]) elif arg.startswith('-sets='): opts['sets'] = int(arg[6:]) + elif arg.startswith('-seed='): opts['seed'] = int(arg[6:]) elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:] elif arg.startswith('-cutoff='): opts['cutoff'] = arg[8:] elif arg.startswith('-title='): opts['title'] = arg[7:] diff --git a/sasmodels/direct_model.py b/sasmodels/direct_model.py index 7c1004ff..e9c133c6 100644 --- a/sasmodels/direct_model.py +++ b/sasmodels/direct_model.py @@ -453,14 +453,18 @@ def test_reparameterize(): except Exception: pass -def _direct_calculate(model, data, pars): +def _direct_calculate(model, data, pars, ngauss=0): from .core import build_model, load_model_info + from .generate import set_integration_size + model_info = load_model_info(model) + if ngauss != 0: + set_integration_size(model_info, ngauss) kernel = build_model(model_info) calculator = DirectModel(data, kernel) return calculator(**pars) -def Iq(model, q, dq=None, ql=None, qw=None, **pars): +def Iq(model, q, dq=None, ql=None, qw=None, ngauss=0, **pars): """ Compute I(q) for *model*. Resolution is *dq* for pinhole or *ql* and *qw* for slit geometry. Use 0 or None for infinite slits. @@ -498,16 +502,16 @@ def broadcast(v): else np.full(len(q), v) if np.isscalar(v) else _as_numpy(v)) data.dxl, data.dxw = broadcast(ql), broadcast(qw) - return _direct_calculate(model, data, pars) + return _direct_calculate(model, data, pars, ngauss=ngauss) -def Iqxy(model, qx, qy, dqx=None, dqy=None, **pars): +def Iqxy(model, qx, qy, dqx=None, dqy=None, ngauss=0, **pars): """ Compute I(qx, qy) for *model*. Resolution is *dqx* and *dqy*. See :func:`Iq` for details on model and parameters. """ from .data import Data2D data = Data2D(x=qx, y=qy, dx=dqx, dy=dqy) - return _direct_calculate(model, data, pars) + return _direct_calculate(model, data, pars, ngauss=ngauss) def Gxi(model, xi, **pars): """ @@ -528,6 +532,8 @@ def main(): if len(sys.argv) < 3: print("usage: python -m sasmodels.direct_model modelname (q|qx,qy) par=val ...") sys.exit(1) + + ngauss = 0 model = sys.argv[1] call = sys.argv[2].upper() pars = dict((k, (float(v) if not k.endswith("_pd_type") else v)) @@ -542,13 +548,13 @@ def main(): dq = dqw = dql = None #dq = [q*0.05] # 5% pinhole resolution dqw, dql = [q*0.05], [1.0] # 5% horizontal slit resolution - print(Iq(model, [q], dq=dq, qw=dqw, ql=dql, **pars)[0]) + print(Iq(model, [q], dq=dq, qw=dqw, ql=dql, ngauss=ngauss, **pars)[0]) #print(Gxi(model, [q], **pars)[0]) elif len(values) == 2: qx, qy = values dq = None #dq = [0.005] # 5% pinhole resolution at q = 0.1 - print(Iqxy(model, [qx], [qy], dqx=dq, dqy=dq, **pars)[0]) + print(Iqxy(model, [qx], [qy], dqx=dq, dqy=dq, ngauss=ngauss, **pars)[0]) else: print("use q or qx,qy") sys.exit(1) diff --git a/sasmodels/generate.py b/sasmodels/generate.py index 0ae362d6..2bc58b18 100644 --- a/sasmodels/generate.py +++ b/sasmodels/generate.py @@ -291,13 +291,27 @@ def set_integration_size(info, n): Note: this really ought to be a method in modelinfo, but that leads to import loops. """ - if info.source and any(lib.startswith('lib/gauss') for lib in info.source): - from .gengauss import gengauss - path = joinpath(MODEL_PATH, "lib", "gauss%d.c"%n) - if not exists(path): - gengauss(n, path) - info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss') - else lib for lib in info.source] + from .gengauss import gengauss + + if not info.source: + return + + # Generate the integration points + path = joinpath(MODEL_PATH, "lib", f"gauss{n}.c") + if not exists(path): + # print(f"building Gaussian integration points of size {n} in {str(path)}") + gengauss(n, path) + + # Replace adaptive.c or lib/gauss.c + try: + index = info.source.index("lib/adaptive.c") + info.source[index:index+1] = [f"lib/gauss{n}.c", "lib/nonadaptive.c"] + except ValueError: + for index in range(len(info.source)-1, -1, -1): + if info.source[index].startswith("lib/gauss"): + info.source[index] = f"lib/gauss{n}.c" + break + # print("info.source is now", info.source) def format_units(units): # type: (str) -> str diff --git a/sasmodels/gengauss.py b/sasmodels/gengauss.py index a9ce027f..2ef06313 100755 --- a/sasmodels/gengauss.py +++ b/sasmodels/gengauss.py @@ -21,19 +21,18 @@ def gengauss(n, path): array_size = n with open(path, "w") as fid: - fid.write("""\ -// Generated by sasmodels.gengauss.gengauss(%d) + fid.write(f"""\ +// Generated by sasmodels.gengauss.gengauss({n}) #ifdef GAUSS_N # undef GAUSS_N # undef GAUSS_Z # undef GAUSS_W #endif -#define GAUSS_N %d -#define GAUSS_Z Gauss%dZ -#define GAUSS_W Gauss%dWt - -"""%(n, n, n, n)) +#define GAUSS_N {n} +#define GAUSS_Z Gauss{n}Z +#define GAUSS_W Gauss{n}Wt +""") if array_size != n: fid.write("// Note: using array size %d so that it is a multiple of 4\n\n"%array_size) diff --git a/sasmodels/model_test.py b/sasmodels/model_test.py index 390df745..80fec08e 100755 --- a/sasmodels/model_test.py +++ b/sasmodels/model_test.py @@ -579,6 +579,31 @@ def _build_test(test): for test in tests: yield _build_test(test) +def _generate_target_values(modelname, ngauss=0): + from .generate import set_integration_size + + model_info = load_model_info(modelname) + if ngauss != 0: + set_integration_size(model_info, ngauss) + model = build_model(model_info, platform="dll", dtype="d") + + for pars, q, Iq in model_info.tests: + qin = q + if isinstance(Iq, float): + q, Iq = [q], [Iq] + if isinstance(q[0], tuple): + qx, qy = zip(*q) + q_vectors = [np.array(qx), np.array(qy)] + else: + q_vectors = [np.array(q)] + kernel = model.make_kernel(q_vectors) + target = np.array(Iq) + actual = call_kernel(kernel, pars) + if True or (actual != target).any(): + print("Test:", modelname, pars) + print(" q = ", qin) + print(f" current => [{', '.join(f'{v:.15g}' for v in target)}]") + print(f" ngauss={ngauss} => [{', '.join(f'{v:.15g}' for v in actual)}]") def main(): # type: () -> int @@ -601,6 +626,11 @@ def main(): help="Engines on which to run the test. " "Valid values are opencl, cuda, dll, and all. " "Defaults to all if no value is given") + parser.add_argument("-t", "--targets", action="store_true", + help="Generate target values for test.") + parser.add_argument("--ngauss", type=int, default=10000, + help="Number of gauss points to use in integration for " + "target values. Warning: this is very slow the first time.") parser.add_argument("models", nargs="*", help="The names of the models to be tested. " "If the first model is 'all', then all but the listed " @@ -630,9 +660,13 @@ def main(): print("unknown engine " + opts.engine) return 1 - runner = TestRunner(verbosity=opts.verbose, **test_args) - result = runner.run(make_suite(loaders, opts.models)) - return 1 if result.failures or result.errors else 0 + if opts.targets: + for model in opts.models: + _generate_target_values(model, ngauss=opts.ngauss) + else: + runner = TestRunner(verbosity=opts.verbose, **test_args) + result = runner.run(make_suite(loaders, opts.models)) + return 1 if result.failures or result.errors else 0 if __name__ == "__main__": diff --git a/sasmodels/models/barbell.c b/sasmodels/models/barbell.c index 87f1553e..f0b4c5c5 100644 --- a/sasmodels/models/barbell.c +++ b/sasmodels/models/barbell.c @@ -19,13 +19,18 @@ _bell_kernel(double qab, double qc, double h, double radius_bell, const double m = radius_bell*qc; // cos argument slope const double b = (half_length+h)*qc; // cos argument intercept const double qab_r = radius_bell*qab; // Q*R*sin(theta) + + const double qr_max = fmax(qab_r, m+b); + constant double *w, *z; + int n = gauss_weights(qr_max, &w, &z); + double total = 0.0; - for (int i = 0; i < GAUSS_N; i++){ - const double t = GAUSS_Z[i]*zm + zb; + for (int i = 0; i < n; i++){ + const double t = z[i]*zm + zb; const double radical = 1.0 - t*t; const double bj = sas_2J1x_x(qab_r*sqrt(radical)); const double Fq = cos(m*t + b) * radical * bj; - total += GAUSS_W[i] * Fq; + total += w[i] * Fq; } // translate dx in [-1,1] to dx in [lower,upper] const double integral = total*zm; @@ -110,21 +115,30 @@ Fq(double q,double *F1, double *F2, double sld, double solvent_sld, const double h = sqrt(square(radius_bell) - square(radius)); const double half_length = 0.5*length; + // The term h comes from solving the right triangle with diagonal + // equal to the bell radius and horizontal equal to the bar radius. + // The result is the height of the equator above the end of the rod. + // To get the total length of bar+bell use bar length + 2*(bell radius + h). + // We want the radius, so divide that by two. + const double qr_max = q*(half_length + radius_bell + h); + constant double *w, *z; + int n = gauss_weights(qr_max, &w, &z); + // translate a point in [-1,1] to a point in [0, pi/2] const double zm = M_PI_4; const double zb = M_PI_4; double total_F1 = 0.0; double total_F2 = 0.0; - for (int i = 0; i < GAUSS_N; i++){ - const double theta = GAUSS_Z[i]*zm + zb; + for (int i = 0; i < n; i++){ + const double theta = z[i]*zm + zb; double sin_theta, cos_theta; // slots to hold sincos function output SINCOS(theta, sin_theta, cos_theta); const double qab = q*sin_theta; const double qc = q*cos_theta; const double Aq = _fq(qab, qc, h, radius_bell, radius, half_length); // scale by sin_theta for spherical coord integration - total_F1 += GAUSS_W[i] * Aq * sin_theta; - total_F2 += GAUSS_W[i] * Aq * Aq * sin_theta; + total_F1 += w[i] * Aq * sin_theta; + total_F2 += w[i] * Aq * Aq * sin_theta; } // translate dx in [-1,1] to dx in [lower,upper] const double form_avg = total_F1 * zm; diff --git a/sasmodels/models/barbell.py b/sasmodels/models/barbell.py index c29b7954..13b4be74 100644 --- a/sasmodels/models/barbell.py +++ b/sasmodels/models/barbell.py @@ -117,7 +117,7 @@ ] # pylint: enable=bad-whitespace, line-too-long -source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"] +source = ["lib/polevl.c", "lib/sas_J1.c", "lib/adaptive.c", "barbell.c"] valid = "radius_bell >= radius" have_Fq = True radius_effective_modes = [ diff --git a/sasmodels/models/capped_cylinder.c b/sasmodels/models/capped_cylinder.c index c792574c..3c369c86 100644 --- a/sasmodels/models/capped_cylinder.c +++ b/sasmodels/models/capped_cylinder.c @@ -7,7 +7,7 @@ // length is the cylinder length, or the separation between the lens halves // theta is the angle of the cylinder wrt q. static double -_cap_kernel(double qab, double qc, double h, double radius_cap, +_cap_kernel(double qab, double qc, double h, double radius_cap, double radius, double half_length) { // translate a point in [-1,1] to a point in [lower,upper] @@ -26,13 +26,25 @@ _cap_kernel(double qab, double qc, double h, double radius_cap, const double m = radius_cap*qc; // cos argument slope const double b = (half_length+h)*qc; // cos argument intercept const double qab_r = radius_cap*qab; // Q*R*sin(theta) + + // m+b = qc*(half_length + radius_cap + h). With h in [-radius_cap, 0] depending + // on cylinder radius, that means m+b is in qc*[length/2, length_2 + radius_cap]. + // The qab_r term will be very large for mostly flat caps. Since the bj term will + // oscillate at this frequency, it seems like we should increase the number of + // gauss points to accomodate. However, if we use the radius of the cylinder + // we seem to get good results, so use that to set the number of integration points. + //const double qr_max = fmax(qab_r, m+b); + const double qr_max = fmax(qab*radius, m+b); + constant double *w, *z; + int n = gauss_weights(qr_max, &w, &z); + double total = 0.0; - for (int i=0; i [0, 1] const double m = 0.5; const double b = 0.5; double total_F1 = 0.0; //initialize intergral double total_F2 = 0.0; //initialize intergral - for(int i=0;i 20, 120 => 76, 1200 => 500, 14500 => 5000 + // The error threshhold is abrupt: if n is too low the result is bad, but + // there is little benefit for having too large an n. + // These results are specific to the pringle inner integral and may not hold in general. + // Check that all models have reasonable cutoff. + // *w = Gauss5000Wt; *z = Gauss5000Z; return 5000; // max precision + // *w = Gauss76Wt; *z = Gauss76Z; return 76; // default + if (qr < 10) { + *w = Gauss20Wt; *z = Gauss20Z; return 20; + } else if (qr < 100) { + *w = Gauss76Wt; *z = Gauss76Z; return 76; + } else if (qr < 750) { + *w = Gauss500Wt; *z = Gauss500Z; return 500; + } else { + *w = Gauss5000Wt; *z = Gauss5000Z; return 5000; + } +} \ No newline at end of file diff --git a/sasmodels/models/lib/nonadaptive.c b/sasmodels/models/lib/nonadaptive.c new file mode 100644 index 00000000..993121ba --- /dev/null +++ b/sasmodels/models/lib/nonadaptive.c @@ -0,0 +1,11 @@ +// To force a fixed rather than adaptive integration scheme, replace [..., "lib/adaptive.c", ...] +// with [..., "lib/gauss.c", "lib/nonadaptive.c", ...] in your source lists. + +#define GAUSS_76 0.0 // Doesn't matter since qr is ignored +#define GAUSS_500 0.0 // Doesn't matter since qr is ignored + +static int +gauss_weights(double qr, constant double **w, constant double **z) +{ + *w = GAUSS_W; *z = GAUSS_Z; return GAUSS_N; +} \ No newline at end of file diff --git a/sasmodels/models/octahedron_truncated.c b/sasmodels/models/octahedron_truncated.c index 679abeff..093ef433 100644 --- a/sasmodels/models/octahedron_truncated.c +++ b/sasmodels/models/octahedron_truncated.c @@ -1,5 +1,5 @@ -#include -#include +//#include +//#include //truncated octahedron volume // NOTE: needs to be called form_volume() for a shape category @@ -11,96 +11,13 @@ form_volume(double length_a, double b2a_ratio, double c2a_ratio, double t) // length_c is the half height along the c axis of the octahedron without truncature // b2a_ratio is length_b divided by Length_a // c2a_ratio is Length_c divided by Length_a -// t varies from 0.5 (cuboctahedron) to 1 (octahedron) - return (4./3.) * cube(length_a) * b2a_ratio * c2a_ratio *(1.-3*cube(1.-t)); +// t varies from 0 (octahedron) to 0.5 (cuboctahedron) +// updated convention with t=0 for on truncation +// t and tinv are exchanged starting from the previous version of the code + const double tinv = 1.0 - t; + return (4./3.) * cube(length_a) * b2a_ratio * c2a_ratio *(1.-3*cube(t)); } -// remark: Iq() is generally not used because have_Fq is set to True in the Python file -static double -Iq(double q, - double sld, - double solvent_sld, - double length_a, - double b2a_ratio, - double c2a_ratio, - double t) -{ - const double length_b = length_a * b2a_ratio; - const double length_c = length_a * c2a_ratio; - - - //Integration limits to use in Gaussian quadrature - const double v1a = 0.0; - const double v1b = M_PI_2; //theta integration limits - const double v2a = 0.0; - const double v2b = M_PI_2; //phi integration limits - - double outer_sum = 0.0; - for(int i=0; i Lb > La so that the integration will go faster. + const double maybe_min = fmin(length_a, length_b); + const double maybe_max = fmax(length_a, length_b); + const double maybe_mid = fmax(maybe_min, length_c); + const double La = fmin(maybe_min, length_c); + const double Lb = fmin(maybe_max, maybe_mid); + const double Lc = fmax(maybe_max, maybe_mid); + + // Find the circumradius by truncating the line from Lc to Lb. This chops of + // the similar triangle with sides (t)Lc, (t)Lb, leaving coordinate ((1-t)Lc, tLb). + // The distance to the origin then follows. + const double qr_max_outer = q*sqrt(square(tinv*Lc) + square(t*Lb)); + constant double *z_outer, *w_outer; + int n_outer = gauss_weights(qr_max_outer, &w_outer, &z_outer); + + const double qr_max_inner = q*sqrt(square(tinv*Lb) + square(t*La)); + constant double *z_inner, *w_inner; + int n_inner = gauss_weights(qr_max_inner, &w_inner, &z_inner); + + //printf("La=%g Lb=%g Lc=%g npoints = %d x %d = %d\n", La, Lb, Lc, n_outer, n_inner, n_outer*n_inner); - - //Integration limits to use in Gaussian quadrature const double v1a = 0.0; const double v1b = M_PI_2; //theta integration limits const double v2a = 0.0; @@ -127,15 +64,19 @@ Fq(double q, double outer_sum_F1 = 0.0; double outer_sum_F2 = 0.0; - for(int i=0; i [0, 1] + const double sin_theta = sqrt(1.0 - square(cos_theta)); // = sin(acos(cos_theta)) + const double qc = q * cos_theta; + const double qz = qc * Lc; + const double qz2 = square(qz); double inner_sum_F1 = 0.0; double inner_sum_F2 = 0.0; - for(int j=0; j [0, pi/2] double sin_phi, cos_phi; SINCOS(phi, sin_phi, cos_phi); @@ -143,37 +84,47 @@ Fq(double q, // q is the modulus of the scattering vector in [A-1] // NOTE: capital QX QY QZ are the three components in [A-1] of the scattering vector // NOTE: qx qy qz are rescaled components (no unit) for computing AA, BB and CC terms - const double Qx = q * sin_theta * cos_phi; - const double Qy = q * sin_theta * sin_phi; - const double Qz = q * cos_theta; - const double qx = Qx * length_a; - const double qy = Qy * length_b; - const double qz = Qz * length_c; - const double AA = 1./(2*(qy*qy-qz*qz)*(qy*qy-qx*qx))*((qy-qx)*sin(qy*(1.-t)-qx*t)+(qy+qx)*sin(qy*(1.-t)+qx*t))+ - 1./(2*(qz*qz-qx*qx)*(qz*qz-qy*qy))*((qz-qx)*sin(qz*(1.-t)-qx*t)+(qz+qx)*sin(qz*(1.-t)+qx*t)); - - const double BB = 1./(2*(qz*qz-qx*qx)*(qz*qz-qy*qy))*((qz-qy)*sin(qz*(1.-t)-qy*t)+(qz+qy)*sin(qz*(1.-t)+qy*t))+ - 1./(2*(qx*qx-qy*qy)*(qx*qx-qz*qz))*((qx-qy)*sin(qx*(1.-t)-qy*t)+(qx+qy)*sin(qx*(1.-t)+qy*t)); - - const double CC = 1./(2*(qx*qx-qy*qy)*(qx*qx-qz*qz))*((qx-qz)*sin(qx*(1.-t)-qz*t)+(qx+qz)*sin(qx*(1.-t)+qz*t))+ - 1./(2*(qy*qy-qz*qz)*(qy*qy-qx*qx))*((qy-qz)*sin(qy*(1.-t)-qz*t)+(qy+qz)*sin(qy*(1.-t)+qz*t)); - - // normalisation to 1. of AP at q = 0. Division by a Factor 4/3. - const double AP = 6./(1.-3*(1.-t)*(1.-t)*(1.-t))*(AA+BB+CC); - - - inner_sum_F1 += GAUSS_W[j] * AP; - inner_sum_F2 += GAUSS_W[j] * AP * AP; + const double qa = q * sin_theta * cos_phi; + const double qb = q * sin_theta * sin_phi; + const double qx = qa * La; + const double qy = qb * Lb; + + // TODO: test for q=0 and return the limiting value. + // From the equations, lim q -> 0 seems to be O(1/q^2), so it diverges. And indeed, + // for q < 1e-8 the function begins to rise. + // TODO: calculation is unstable for small q. + // PAK: reordered the equations and moved factor of 1/2 to normalization. + // PAK: verified that t=1 matches Wei-Ren Chen et al. for a regular octahedron. + const double qx2 = square(qx); + const double qy2 = square(qy); + const double AA = + ((qy-qx)*sin(qy*t-qx*tinv) + (qy+qx)*sin(qy*t+qx*tinv)) / ((qy2-qz2)*(qy2-qx2)) + + ((qz-qx)*sin(qz*t-qx*tinv) + (qz+qx)*sin(qz*t+qx*tinv)) / ((qz2-qx2)*(qz2-qy2)); + + const double BB = + ((qz-qy)*sin(qz*t-qy*tinv) + (qz+qy)*sin(qz*t+qy*tinv)) / ((qz2-qx2)*(qz2-qy2)) + + ((qx-qy)*sin(qx*t-qy*tinv) + (qx+qy)*sin(qx*t+qy*tinv)) / ((qx2-qy2)*(qx2-qz2)); + + const double CC = + ((qx-qz)*sin(qx*t-qz*tinv) + (qx+qz)*sin(qx*t+qz*tinv)) / ((qx2-qy2)*(qx2-qz2)) + + ((qy-qz)*sin(qy*t-qz*tinv) + (qy+qz)*sin(qy*t+qz*tinv)) / ((qy2-qz2)*(qy2-qx2)); + + // normalisation to 1. of AP at q = 0. Division by a Factor 4/3. + const double AP = 3./(1. - 3.*cube(tinv)) * (AA+BB+CC); + + inner_sum_F1 += w_inner[j] * AP; + inner_sum_F2 += w_inner[j] * AP * AP; } - inner_sum_F1 = 0.5 * (v2b-v2a) * inner_sum_F1; - inner_sum_F2 = 0.5 * (v2b-v2a) * inner_sum_F2; - outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * sin_theta; - outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * sin_theta; + // Note: moving the pi/2 scaling outside the sum. No sin(theta) term because of u=cos(theta) + // TODO: With the u=cos(theta) substitution should the F1 integral be negative? + // TODO: check for this in all models. + outer_sum_F1 += w_outer[i] * inner_sum_F1; + outer_sum_F2 += w_outer[i] * inner_sum_F2; } - - outer_sum_F1 *= 0.5*(v1b-v1a); - outer_sum_F2 *= 0.5*(v1b-v1a); + // sum(w) = 2, so scaling + outer_sum_F1 *= 0.5*M_PI_4; + outer_sum_F2 *= 0.5*M_PI_4; // The factor 2 appears because the theta integral has been defined between // 0 and pi/2, instead of 0 to pi. @@ -181,15 +132,13 @@ Fq(double q, outer_sum_F2 /= M_PI_2; // Multiply by contrast and volume - // contrast - const double s = (sld-solvent_sld); - // volume - // s *= form_volume(length_a, b2a_ratio,c2a_ratio, t); - + const double s = (sld-solvent_sld) * form_volume(length_a, b2a_ratio,c2a_ratio, t); + // Convert from [1e-12 A-1] to [cm-1] - *F1 = 1e-2 * s * form_volume(length_a, b2a_ratio,c2a_ratio, t) * outer_sum_F1; - *F2 = 1e-4 * square(s * form_volume(length_a, b2a_ratio,c2a_ratio, t)) * outer_sum_F2; + *F1 = 1e-2 * s * outer_sum_F1; + *F2 = 1e-4 * square(s) * outer_sum_F2; + // TODO: inf and NaN shouldn't occur, but if they do they shouldn't be replaced by zero. if (isnan(*F1) || isinf(*F1)) { *F1 = 0.0; } @@ -210,25 +159,35 @@ Iqabc(double qa, double qb, double qc, { const double length_b = length_a * b2a_ratio; const double length_c = length_a * c2a_ratio; - + const double tinv = 1.0 - t; //HERE: Octahedron formula - // NOTE: qa qb qc are the three components in [A-1] of the scattering vector + // NOTE: qa qb qc are the three components in 1/Ang of the scattering vector // NOTE: qx qy qz are rescaled components (no unit) for computing AA, BB and CC terms const double qx = qa * length_a; const double qy = qb * length_b; const double qz = qc * length_c; - const double AA = 1./(2*(qy*qy-qz*qz)*(qy*qy-qx*qx))*((qy-qx)*sin(qy*(1.-t)-qx*t)+(qy+qx)*sin(qy*(1.-t)+qx*t))+ - 1./(2*(qz*qz-qx*qx)*(qz*qz-qy*qy))*((qz-qx)*sin(qz*(1.-t)-qx*t)+(qz+qx)*sin(qz*(1.-t)+qx*t)); - const double BB = 1./(2*(qz*qz-qx*qx)*(qz*qz-qy*qy))*((qz-qy)*sin(qz*(1.-t)-qy*t)+(qz+qy)*sin(qz*(1.-t)+qy*t))+ - 1./(2*(qx*qx-qy*qy)*(qx*qx-qz*qz))*((qx-qy)*sin(qx*(1.-t)-qy*t)+(qx+qy)*sin(qx*(1.-t)+qy*t)); + // TODO: calculation is unstable for small q. + // PAK: reordered the equations and moved factor of 1/2 to normalization. + const double qx2 = square(qx); + const double qy2 = square(qy); + const double qz2 = square(qz); + const double AA = + ((qy-qx)*sin(qy*t-qx*tinv) + (qy+qx)*sin(qy*t+qx*tinv)) / ((qy2-qz2)*(qy2-qx2)) + + ((qz-qx)*sin(qz*t-qx*tinv) + (qz+qx)*sin(qz*t+qx*tinv)) / ((qz2-qx2)*(qz2-qy2)); + + const double BB = + ((qz-qy)*sin(qz*t-qy*tinv) + (qz+qy)*sin(qz*t+qy*tinv)) / ((qz2-qx2)*(qz2-qy2)) + + ((qx-qy)*sin(qx*t-qy*tinv) + (qx+qy)*sin(qx*t+qy*tinv)) / ((qx2-qy2)*(qx2-qz2)); + + const double CC = + ((qx-qz)*sin(qx*t-qz*tinv) + (qx+qz)*sin(qx*t+qz*tinv)) / ((qx2-qy2)*(qx2-qz2)) + + ((qy-qz)*sin(qy*t-qz*tinv) + (qy+qz)*sin(qy*t+qz*tinv)) / ((qy2-qz2)*(qy2-qx2)); - const double CC = 1./(2*(qx*qx-qy*qy)*(qx*qx-qz*qz))*((qx-qz)*sin(qx*(1.-t)-qz*t)+(qx+qz)*sin(qx*(1.-t)+qz*t))+ - 1./(2*(qy*qy-qz*qz)*(qy*qy-qx*qx))*((qy-qz)*sin(qy*(1.-t)-qz*t)+(qy+qz)*sin(qy*(1.-t)+qz*t)); + // normalisation to 1. of AP at q = 0. Division by a Factor 4/3. + const double AP = 3./(1. - 3.*cube(t)) * (AA+BB+CC); - // normalisation to 1. of AP at q = 0. Division by a Factor 4/3. - const double AP = 6./(1.-3*(1.-t)*(1.-t)*(1.-t))*(AA+BB+CC); // Multiply by contrast and volume // contrast diff --git a/sasmodels/models/octahedron_truncated.py b/sasmodels/models/octahedron_truncated.py index d2f3c454..ad63b88d 100644 --- a/sasmodels/models/octahedron_truncated.py +++ b/sasmodels/models/octahedron_truncated.py @@ -253,7 +253,9 @@ "rotation about c axis"], ] -source = ["lib/gauss20.c", "octahedron_truncated.c"] +valid = "t >= 0.5 && t <= 1.0" +single = False +source = ["lib/adaptive.c", "octahedron_truncated.c"] # change to "lib/gauss76.c" or "lib/gauss150.c" to increase the number of integration points # for the orientational average. Note that it will increase calculation times. diff --git a/sasmodels/models/parallelepiped.c b/sasmodels/models/parallelepiped.c index 7b5262b6..876dce3b 100644 --- a/sasmodels/models/parallelepiped.c +++ b/sasmodels/models/parallelepiped.c @@ -75,34 +75,42 @@ Fq(double q, const double a_scaled = length_a / length_b; const double c_scaled = length_c / length_b; + const double qr_max_outer = q*fmax(sqrt(length_a*length_a + length_b*length_b), length_c)/2; + constant double *w_outer, *z_outer; + int n_outer = gauss_weights(qr_max_outer, &w_outer, &z_outer); + // outer integral (with gauss points), integration limits = 0, 1 double outer_total_F1 = 0.0; //initialize integral double outer_total_F2 = 0.0; //initialize integral - for( int i=0; i