source: sasmodels/sasmodels/models/spherical_sld.c @ c3ccaec

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since c3ccaec was c3ccaec, checked in by GitHub <noreply@…>, 7 years ago

Merge branch 'master' into ticket-815

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