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

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since a34b811 was a34b811, checked in by Paul Kienzle <pkienzle@…>, 5 months ago

use radius_effective/radius_effective_mode/radius_effective_modes consistently throughout the code

  • Property mode set to 100644
File size: 5.5 KB
Line 
1static double
2outer_radius(double fp_n_shells, double thickness[], double interface[])
3{
4    int n_shells= (int)(fp_n_shells + 0.5);
5    double r = 0.0;
6    for (int i=0; i < n_shells; i++) {
7        r += thickness[i] + interface[i];
8    }
9    return r;
10}
11
12static double form_volume(
13    double fp_n_shells,
14    double thickness[],
15    double interface[])
16{
17    return M_4PI_3*cube(outer_radius(fp_n_shells, thickness, interface));
18}
19
20static double
21radius_effective(int mode, double fp_n_shells, double thickness[], double interface[])
22{
23    // case 1: outer radius
24    return outer_radius(fp_n_shells, thickness, interface);
25}
26
27static double blend(int shape, double nu, double z)
28{
29    if (shape==0) {
30        const double num = sas_erf(nu * M_SQRT1_2 * (2.0*z - 1.0));
31        const double denom = 2.0 * sas_erf(nu * M_SQRT1_2);
32        return num/denom + 0.5;
33    } else if (shape==1) {
34        return pow(z, nu);
35    } else if (shape==2) {
36        return 1.0 - pow(1.0 - z, nu);
37    } else if (shape==3) {
38        return expm1(-nu*z)/expm1(-nu);
39    } else if (shape==4) {
40        return expm1(nu*z)/expm1(nu);
41    } else {
42        return NAN;
43    }
44}
45
46static double f_linear(double q, double r, double contrast, double slope)
47{
48    const double qr = q * r;
49    const double qrsq = qr * qr;
50    const double bes = sas_3j1x_x(qr);
51    double sinqr, cosqr;
52    SINCOS(qr, sinqr, cosqr);
53    const double fun = 3.0*r*(2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq);
54    const double vol = M_4PI_3 * cube(r);
55    return vol*(bes*contrast + fun*slope);
56}
57
58static double Iq(
59    double q,
60    double fp_n_shells,
61    double sld_solvent,
62    double sld[],
63    double thickness[],
64    double interface[],
65    double shape[],
66    double nu[],
67    double fp_n_steps)
68{
69    // iteration for # of shells + core + solvent
70    int n_shells = (int)(fp_n_shells + 0.5);
71    int n_steps = (int)(fp_n_steps + 0.5);
72    double f=0.0;
73    double r=0.0;
74    for (int shell=0; shell<n_shells; shell++){
75        const double sld_l = sld[shell];
76
77        // uniform shell; r=0 => r^3=0 => f=0, so works for core as well.
78        f -= M_4PI_3 * cube(r) * sld_l * sas_3j1x_x(q*r);
79        r += thickness[shell];
80        f += M_4PI_3 * cube(r) * sld_l * sas_3j1x_x(q*r);
81
82        // iterate over sub_shells in the interface
83        const double dr = interface[shell]/n_steps;
84        const double delta = (shell==n_shells-1 ? sld_solvent : sld[shell+1]) - sld_l;
85        const double nu_shell = fmax(fabs(nu[shell]), 1.e-14);
86        const int shape_shell = (int)(shape[shell]);
87
88        // if there is no interface the equations don't work
89        if (dr == 0.) continue;
90
91        double sld_in = sld_l;
92        for (int step=1; step <= n_steps; step++) {
93            // find sld_i at the outer boundary of sub-shell step
94            //const double z = (double)step/(double)n_steps;
95            const double z = (double)step/(double)n_steps;
96            const double fraction = blend(shape_shell, nu_shell, z);
97            const double sld_out = fraction*delta + sld_l;
98            // calculate slope
99            const double slope = (sld_out - sld_in)/dr;
100            const double contrast = sld_in - slope*r;
101
102            // iteration for the left and right boundary of the shells
103            f -= f_linear(q, r, contrast, slope);
104            r += dr;
105            f += f_linear(q, r, contrast, slope);
106            sld_in = sld_out;
107        }
108    }
109    // add in solvent effect
110    f -= M_4PI_3 * cube(r) * sld_solvent * sas_3j1x_x(q*r);
111
112    const double f2 = f * f * 1.0e-4;
113    return f2;
114}
115
116static void Fq(
117    double q,
118    double *F1,
119    double *F2,
120    double fp_n_shells,
121    double sld_solvent,
122    double sld[],
123    double thickness[],
124    double interface[],
125    double shape[],
126    double nu[],
127    double fp_n_steps)
128{
129    // iteration for # of shells + core + solvent
130    int n_shells = (int)(fp_n_shells + 0.5);
131    int n_steps = (int)(fp_n_steps + 0.5);
132    double f=0.0;
133    double r=0.0;
134    for (int shell=0; shell<n_shells; shell++){
135        const double sld_l = sld[shell];
136
137        // uniform shell; r=0 => r^3=0 => f=0, so works for core as well.
138        f -= M_4PI_3 * cube(r) * sld_l * sas_3j1x_x(q*r);
139        r += thickness[shell];
140        f += M_4PI_3 * cube(r) * sld_l * sas_3j1x_x(q*r);
141
142        // iterate over sub_shells in the interface
143        const double dr = interface[shell]/n_steps;
144        const double delta = (shell==n_shells-1 ? sld_solvent : sld[shell+1]) - sld_l;
145        const double nu_shell = fmax(fabs(nu[shell]), 1.e-14);
146        const int shape_shell = (int)(shape[shell]);
147
148        // if there is no interface the equations don't work
149        if (dr == 0.) continue;
150
151        double sld_in = sld_l;
152        for (int step=1; step <= n_steps; step++) {
153            // find sld_i at the outer boundary of sub-shell step
154            //const double z = (double)step/(double)n_steps;
155            const double z = (double)step/(double)n_steps;
156            const double fraction = blend(shape_shell, nu_shell, z);
157            const double sld_out = fraction*delta + sld_l;
158            // calculate slope
159            const double slope = (sld_out - sld_in)/dr;
160            const double contrast = sld_in - slope*r;
161
162            // iteration for the left and right boundary of the shells
163            f -= f_linear(q, r, contrast, slope);
164            r += dr;
165            f += f_linear(q, r, contrast, slope);
166            sld_in = sld_out;
167        }
168    }
169    // add in solvent effect
170    f -= M_4PI_3 * cube(r) * sld_solvent * sas_3j1x_x(q*r);
171
172    *F1 = 1e-2*f;
173    *F2 = 1e-4*f*f;
174}
Note: See TracBrowser for help on using the repository browser.