source: sasmodels/sasmodels/models/core_shell_parallelepiped.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, 8 months ago

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

  • Property mode set to 100644
File size: 8.4 KB
Line 
1// Set OVERLAPPING to 1 in order to fill in the edges of the box, with
2// c endcaps and b overlapping a.  With the proper choice of parameters,
3// (setting rim slds to sld, core sld to solvent, rim thickness to thickness
4// and subtracting 2*thickness from length, this should match the hollow
5// rectangular prism.)  Set it to 0 for the documented behaviour.
6#define OVERLAPPING 0
7static double
8form_volume(double length_a, double length_b, double length_c,
9    double thick_rim_a, double thick_rim_b, double thick_rim_c)
10{
11    return
12#if OVERLAPPING
13        // Hollow rectangular prism only includes the volume of the shell
14        // so uncomment the next line when comparing.  Solid rectangular
15        // prism, or parallelepiped want filled cores, so comment when
16        // comparing.
17        //-length_a * length_b * length_c +
18        (length_a + 2.0*thick_rim_a) *
19        (length_b + 2.0*thick_rim_b) *
20        (length_c + 2.0*thick_rim_c);
21#else
22        length_a * length_b * length_c +
23        2.0 * thick_rim_a * length_b * length_c +
24        2.0 * length_a * thick_rim_b * length_c +
25        2.0 * length_a * length_b * thick_rim_c;
26#endif
27}
28
29static double
30radius_from_excluded_volume(double length_a, double length_b, double length_c,
31                   double thick_rim_a, double thick_rim_b, double thick_rim_c)
32{
33    double r_equiv, length;
34    double lengths[3] = {length_a+thick_rim_a, length_b+thick_rim_b, length_c+thick_rim_c};
35    double lengthmax = fmax(lengths[0],fmax(lengths[1],lengths[2]));
36    double length_1 = lengthmax;
37    double length_2 = lengthmax;
38    double length_3 = lengthmax;
39
40    for(int ilen=0; ilen<3; ilen++) {
41        if (lengths[ilen] < length_1) {
42            length_2 = length_1;
43            length_1 = lengths[ilen];
44            } else {
45                if (lengths[ilen] < length_2) {
46                        length_2 = lengths[ilen];
47                }
48            }
49    }
50    if(length_2-length_1 > length_3-length_2) {
51        r_equiv = sqrt(length_2*length_3/M_PI);
52        length  = length_1;
53    } else  {
54        r_equiv = sqrt(length_1*length_2/M_PI);
55        length  = length_3;
56    }
57
58    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length + (r_equiv + length)*(M_PI*r_equiv + length)));
59}
60
61static double
62radius_from_volume(double length_a, double length_b, double length_c,
63                   double thick_rim_a, double thick_rim_b, double thick_rim_c)
64{
65    const double volume = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
66    return cbrt(volume/M_4PI_3);
67}
68
69static double
70radius_from_crosssection(double length_a, double length_b, double thick_rim_a, double thick_rim_b)
71{
72    const double area_xsec_paral = length_a*length_b + 2.0*thick_rim_a*length_b + 2.0*thick_rim_b*length_a;
73    return sqrt(area_xsec_paral/M_PI);
74}
75
76static double
77effective_radius(int mode, double length_a, double length_b, double length_c,
78                 double thick_rim_a, double thick_rim_b, double thick_rim_c)
79{
80    switch (mode) {
81    default:
82    case 1: // equivalent cylinder excluded volume
83        return radius_from_excluded_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
84    case 2: // equivalent volume sphere
85        return radius_from_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
86    case 3: // half outer length a
87        return 0.5 * length_a + thick_rim_a;
88    case 4: // half outer length b
89        return 0.5 * length_b + thick_rim_b;
90    case 5: // half outer length c
91        return 0.5 * length_c + thick_rim_c;
92    case 6: // equivalent circular cross-section
93        return radius_from_crosssection(length_a, length_b, thick_rim_a, thick_rim_b);
94    case 7: // half outer ab diagonal
95        return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b));
96    case 8: // half outer diagonal
97        return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b) + square(length_c+ 2.0*thick_rim_c));
98    }
99}
100
101static void
102Fq(double q,
103    double *F1,
104    double *F2,
105    double core_sld,
106    double arim_sld,
107    double brim_sld,
108    double crim_sld,
109    double solvent_sld,
110    double length_a,
111    double length_b,
112    double length_c,
113    double thick_rim_a,
114    double thick_rim_b,
115    double thick_rim_c)
116{
117    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c
118    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
119    // Code is rewritten, the code is compliant with Diva Singh's thesis now (Dirk Honecker)
120    // Code rewritten; cross checked against hollow rectangular prism and realspace (PAK)
121
122    const double half_q = 0.5*q;
123
124    const double tA = length_a + 2.0*thick_rim_a;
125    const double tB = length_b + 2.0*thick_rim_b;
126    const double tC = length_c + 2.0*thick_rim_c;
127
128    // Scale factors
129    const double dr0 = (core_sld-solvent_sld);
130    const double drA = (arim_sld-solvent_sld);
131    const double drB = (brim_sld-solvent_sld);
132    const double drC = (crim_sld-solvent_sld);
133
134    // outer integral (with gauss points), integration limits = 0, 1
135    // substitute d_cos_alpha for sin_alpha d_alpha
136    double outer_sum_F1 = 0; //initialize integral
137    double outer_sum_F2 = 0; //initialize integral
138    for( int i=0; i<GAUSS_N; i++) {
139        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 );
140        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha);
141        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q);
142        const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q);
143
144        // inner integral (with gauss points), integration limits = 0, 1
145        // substitute beta = PI/2 u (so 2/PI * d_(PI/2 * beta) = d_beta)
146        double inner_sum_F1 = 0.0;
147        double inner_sum_F2 = 0.0;
148        for(int j=0; j<GAUSS_N; j++) {
149            const double u = 0.5 * ( GAUSS_Z[j] + 1.0 );
150            double sin_beta, cos_beta;
151            SINCOS(M_PI_2*u, sin_beta, cos_beta);
152            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta);
153            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta);
154            const double siAt = tA * sas_sinx_x(tA * mu * sin_beta);
155            const double siBt = tB * sas_sinx_x(tB * mu * cos_beta);
156
157#if OVERLAPPING
158            const double f = dr0*siA*siB*siC
159                + drA*(siAt-siA)*siB*siC
160                + drB*siAt*(siBt-siB)*siC
161                + drC*siAt*siBt*(siCt-siC);
162#else
163            const double f = dr0*siA*siB*siC
164                + drA*(siAt-siA)*siB*siC
165                + drB*siA*(siBt-siB)*siC
166                + drC*siA*siB*(siCt-siC);
167#endif
168
169            inner_sum_F1 += GAUSS_W[j] * f;
170            inner_sum_F2 += GAUSS_W[j] * f * f;
171        }
172        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5
173        // and sum up the outer integral
174        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * 0.5;
175        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * 0.5;
176    }
177    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5
178    outer_sum_F1 *= 0.5;
179    outer_sum_F2 *= 0.5;
180
181    //convert from [1e-12 A-1] to [cm-1]
182    *F1 = 1.0e-2 * outer_sum_F1;
183    *F2 = 1.0e-4 * outer_sum_F2;
184}
185
186static double
187Iqabc(double qa, double qb, double qc,
188    double core_sld,
189    double arim_sld,
190    double brim_sld,
191    double crim_sld,
192    double solvent_sld,
193    double length_a,
194    double length_b,
195    double length_c,
196    double thick_rim_a,
197    double thick_rim_b,
198    double thick_rim_c)
199{
200    // cspkernel in csparallelepiped recoded here
201    const double dr0 = core_sld-solvent_sld;
202    const double drA = arim_sld-solvent_sld;
203    const double drB = brim_sld-solvent_sld;
204    const double drC = crim_sld-solvent_sld;
205
206    const double tA = length_a + 2.0*thick_rim_a;
207    const double tB = length_b + 2.0*thick_rim_b;
208    const double tC = length_c + 2.0*thick_rim_c;
209    const double siA = length_a*sas_sinx_x(0.5*length_a*qa);
210    const double siB = length_b*sas_sinx_x(0.5*length_b*qb);
211    const double siC = length_c*sas_sinx_x(0.5*length_c*qc);
212    const double siAt = tA*sas_sinx_x(0.5*tA*qa);
213    const double siBt = tB*sas_sinx_x(0.5*tB*qb);
214    const double siCt = tC*sas_sinx_x(0.5*tC*qc);
215
216#if OVERLAPPING
217    const double f = dr0*siA*siB*siC
218        + drA*(siAt-siA)*siB*siC
219        + drB*siAt*(siBt-siB)*siC
220        + drC*siAt*siBt*(siCt-siC);
221#else
222    const double f = dr0*siA*siB*siC
223        + drA*(siAt-siA)*siB*siC
224        + drB*siA*(siBt-siB)*siC
225        + drC*siA*siB*(siCt-siC);
226#endif
227
228    return 1.0e-4 * f * f;
229}
Note: See TracBrowser for help on using the repository browser.