source: sasmodels/sasmodels/models/rectangular_prism.c @ 99658f6

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 99658f6 was 99658f6, checked in by grethevj, 13 months ago

updated ER functions including cylinder excluded volume, to match 4.x

  • Property mode set to 100644
File size: 6.7 KB
Line 
1static double
2form_volume(double length_a, double b2a_ratio, double c2a_ratio)
3{
4    return length_a * (length_a*b2a_ratio) * (length_a*c2a_ratio);
5}
6
7static double
8radius_from_excluded_volume(double length_a, double b2a_ratio, double c2a_ratio)
9{
10    double const r_equiv   = sqrt(length_a*length_a*b2a_ratio/M_PI);
11    double const length_c  = c2a_ratio*length_a;
12    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length_c + (r_equiv + length_c)*(M_PI*r_equiv + length_c)));
13}
14
15static double
16effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio)
17{
18    switch (mode) {
19    default:
20    case 1: // equivalent cylinder excluded volume
21        return radius_from_excluded_volume(length_a,b2a_ratio,c2a_ratio);
22    case 2: // equivalent volume sphere
23        return cbrt(cube(length_a)*b2a_ratio*c2a_ratio/M_4PI_3);
24    case 3: // half length_a
25        return 0.5 * length_a;
26    case 4: // half length_b
27        return 0.5 * length_a*b2a_ratio;
28    case 5: // half length_c
29        return 0.5 * length_a*c2a_ratio;
30    case 6: // equivalent circular cross-section
31        return length_a*sqrt(b2a_ratio/M_PI);
32    case 7: // half ab diagonal
33        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio)));
34    case 8: // half diagonal
35        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio)));
36    }
37}
38
39static double
40Iq(double q,
41    double sld,
42    double solvent_sld,
43    double length_a,
44    double b2a_ratio,
45    double c2a_ratio)
46{
47    const double length_b = length_a * b2a_ratio;
48    const double length_c = length_a * c2a_ratio;
49    const double a_half = 0.5 * length_a;
50    const double b_half = 0.5 * length_b;
51    const double c_half = 0.5 * length_c;
52
53   //Integration limits to use in Gaussian quadrature
54    const double v1a = 0.0;
55    const double v1b = M_PI_2;  //theta integration limits
56    const double v2a = 0.0;
57    const double v2b = M_PI_2;  //phi integration limits
58
59    double outer_sum = 0.0;
60    for(int i=0; i<GAUSS_N; i++) {
61        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b );
62        double sin_theta, cos_theta;
63        SINCOS(theta, sin_theta, cos_theta);
64
65        const double termC = sas_sinx_x(q * c_half * cos_theta);
66
67        double inner_sum = 0.0;
68        for(int j=0; j<GAUSS_N; j++) {
69            double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b );
70            double sin_phi, cos_phi;
71            SINCOS(phi, sin_phi, cos_phi);
72
73            // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0
74            const double termA = sas_sinx_x(q * a_half * sin_theta * sin_phi);
75            const double termB = sas_sinx_x(q * b_half * sin_theta * cos_phi);
76            const double AP = termA * termB * termC;
77            inner_sum += GAUSS_W[j] * AP * AP;
78        }
79        inner_sum = 0.5 * (v2b-v2a) * inner_sum;
80        outer_sum += GAUSS_W[i] * inner_sum * sin_theta;
81    }
82
83    double answer = 0.5*(v1b-v1a)*outer_sum;
84
85    // Normalize by Pi (Eqn. 16).
86    // The term (ABC)^2 does not appear because it was introduced before on
87    // the definitions of termA, termB, termC.
88    // The factor 2 appears because the theta integral has been defined between
89    // 0 and pi/2, instead of 0 to pi.
90    answer /= M_PI_2; //Form factor P(q)
91
92    // Multiply by contrast^2 and volume^2
93    const double volume = length_a * length_b * length_c;
94    answer *= square((sld-solvent_sld)*volume);
95
96    // Convert from [1e-12 A-1] to [cm-1]
97    answer *= 1.0e-4;
98
99    return answer;
100}
101
102static void
103Fq(double q,
104    double *F1,
105    double *F2,
106    double sld,
107    double solvent_sld,
108    double length_a,
109    double b2a_ratio,
110    double c2a_ratio)
111{
112    const double length_b = length_a * b2a_ratio;
113    const double length_c = length_a * c2a_ratio;
114    const double a_half = 0.5 * length_a;
115    const double b_half = 0.5 * length_b;
116    const double c_half = 0.5 * length_c;
117
118   //Integration limits to use in Gaussian quadrature
119    const double v1a = 0.0;
120    const double v1b = M_PI_2;  //theta integration limits
121    const double v2a = 0.0;
122    const double v2b = M_PI_2;  //phi integration limits
123
124    double outer_sum_F1 = 0.0;
125    double outer_sum_F2 = 0.0;
126    for(int i=0; i<GAUSS_N; i++) {
127        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b );
128        double sin_theta, cos_theta;
129        SINCOS(theta, sin_theta, cos_theta);
130
131        const double termC = sas_sinx_x(q * c_half * cos_theta);
132
133        double inner_sum_F1 = 0.0;
134        double inner_sum_F2 = 0.0;
135        for(int j=0; j<GAUSS_N; j++) {
136            double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b );
137            double sin_phi, cos_phi;
138            SINCOS(phi, sin_phi, cos_phi);
139
140            // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0
141            const double termA = sas_sinx_x(q * a_half * sin_theta * sin_phi);
142            const double termB = sas_sinx_x(q * b_half * sin_theta * cos_phi);
143            const double AP = termA * termB * termC;
144            inner_sum_F1 += GAUSS_W[j] * AP;
145            inner_sum_F2 += GAUSS_W[j] * AP * AP;
146        }
147        inner_sum_F1 = 0.5 * (v2b-v2a) * inner_sum_F1;
148        inner_sum_F2 = 0.5 * (v2b-v2a) * inner_sum_F2;
149        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * sin_theta;
150        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * sin_theta;
151    }
152
153    outer_sum_F1 *= 0.5*(v1b-v1a);
154    outer_sum_F2 *= 0.5*(v1b-v1a);
155
156    // Normalize by Pi (Eqn. 16).
157    // The term (ABC)^2 does not appear because it was introduced before on
158    // the definitions of termA, termB, termC.
159    // The factor 2 appears because the theta integral has been defined between
160    // 0 and pi/2, instead of 0 to pi.
161    outer_sum_F1 /= M_PI_2;
162    outer_sum_F2 /= M_PI_2;
163
164    // Multiply by contrast and volume
165    const double s = (sld-solvent_sld) * (length_a * length_b * length_c);
166
167    // Convert from [1e-12 A-1] to [cm-1]
168    *F1 = 1e-2 * s * outer_sum_F1;
169    *F2 = 1e-4 * s * s * outer_sum_F2;
170}
171
172
173static double
174Iqabc(double qa, double qb, double qc,
175    double sld,
176    double solvent_sld,
177    double length_a,
178    double b2a_ratio,
179    double c2a_ratio)
180{
181    const double length_b = length_a * b2a_ratio;
182    const double length_c = length_a * c2a_ratio;
183    const double a_half = 0.5 * length_a;
184    const double b_half = 0.5 * length_b;
185    const double c_half = 0.5 * length_c;
186
187    // Amplitude AP from eqn. (13)
188    const double termA = sas_sinx_x(qa * a_half);
189    const double termB = sas_sinx_x(qb * b_half);
190    const double termC = sas_sinx_x(qc * c_half);
191    const double AP = termA * termB * termC;
192
193    // Multiply by contrast and volume
194    const double s = (sld-solvent_sld) * (length_a * length_b * length_c);
195
196    // Convert from [1e-12 A-1] to [cm-1]
197    return 1.0e-4 * square(s * AP);
198}
Note: See TracBrowser for help on using the repository browser.