source: sasmodels/sasmodels/models/spherical_sld.c @ 925ad6e

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 925ad6e was 925ad6e, checked in by wojciech, 6 years ago

sph_j1c translated to sas_3j1x_x

  • Property mode set to 100644
File size: 3.1 KB
Line 
1static double form_volume(
2    int n_shells,
3    double thickness[],
4    double interface[])
5{
6    double r = 0.0;
7    for (int i=0; i < n_shells; i++) {
8        r += thickness[i] + interface[i];
9    }
10    return M_4PI_3*cube(r);
11}
12
13static double blend(int shape, double nu, double z)
14{
15    if (shape==0) {
16        const double num = sas_erf(nu * M_SQRT1_2 * (2.0*z - 1.0));
17        const double denom = 2.0 * sas_erf(nu * M_SQRT1_2);
18        return num/denom + 0.5;
19    } else if (shape==1) {
20        return pow(z, nu);
21    } else if (shape==2) {
22        return 1.0 - pow(1. - z, nu);
23    } else if (shape==3) {
24        return expm1(-nu*z)/expm1(-nu);
25    } else if (shape==4) {
26        return expm1(nu*z)/expm1(nu);
27    } else {
28        return NAN;
29    }
30}
31
32static double f_linear(double q, double r, double contrast, double slope)
33{
34    const double qr = q * r;
35    const double qrsq = qr * qr;
36    const double bes = sas_3j1x_x(qr);
37    double sinqr, cosqr;
38    SINCOS(qr, sinqr, cosqr);
39    const double fun = 3.0*r*(2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq);
40    const double vol = M_4PI_3 * cube(r);
41    return vol*(bes*contrast + fun*slope);
42}
43
44static double Iq(
45    double q,
46    int n_shells,
47    double sld_solvent,
48    double sld[],
49    double thickness[],
50    double interface[],
51    double shape[],
52    double nu[],
53    int n_steps)
54{
55    // iteration for # of shells + core + solvent
56    double f=0.0;
57    double r=0.0;
58    for (int shell=0; shell<n_shells; shell++){
59        const double sld_l = sld[shell];
60
61        // uniform shell; r=0 => r^3=0 => f=0, so works for core as well.
62        f -= M_4PI_3 * cube(r) * sld_l * sas_3j1x_x(q*r);
63        r += thickness[shell];
64        f += M_4PI_3 * cube(r) * sld_l * sas_3j1x_x(q*r);
65
66        // iterate over sub_shells in the interface
67        const double dr = interface[shell]/n_steps;
68        const double delta = (shell==n_shells-1 ? sld_solvent : sld[shell+1]) - sld_l;
69        const double nu_shell = fmax(fabs(nu[shell]), 1.e-14);
70        const int shape_shell = (int)(shape[shell]);
71
72        // if there is no interface the equations don't work
73        if (dr == 0.) continue;
74
75        double sld_in = sld_l;
76        for (int step=1; step <= n_steps; step++) {
77            // find sld_i at the outer boundary of sub-shell step
78            //const double z = (double)step/(double)n_steps;
79            const double z = (double)step/(double)n_steps;
80            const double fraction = blend(shape_shell, nu_shell, z);
81            const double sld_out = fraction*delta + sld_l;
82            // calculate slope
83            const double slope = (sld_out - sld_in)/dr;
84            const double contrast = sld_in - slope*r;
85
86            // iteration for the left and right boundary of the shells
87            f -= f_linear(q, r, contrast, slope);
88            r += dr;
89            f += f_linear(q, r, contrast, slope);
90            sld_in = sld_out;
91        }
92    }
93    // add in solvent effect
94    f -= M_4PI_3 * cube(r) * sld_solvent * sas_3j1x_x(q*r);
95
96    const double f2 = f * f * 1.0e-4;
97    return f2;
98}
99
Note: See TracBrowser for help on using the repository browser.