source: sasmodels/sasmodels/models/core_shell_parallelepiped.c @ d42dd4a

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since d42dd4a was d42dd4a, checked in by pkienzle, 8 months ago

fix compiler warnings for CUDA

  • Property mode set to 100644
File size: 7.2 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_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    const double volume = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
34    return cbrt(volume/M_4PI_3);
35}
36
37static double
38radius_from_crosssection(double length_a, double length_b, double thick_rim_a, double thick_rim_b)
39{
40    const double area_xsec_paral = length_a*length_b + 2.0*thick_rim_a*length_b + 2.0*thick_rim_b*length_a;
41    return sqrt(area_xsec_paral/M_PI);
42}
43
44static double
45effective_radius(int mode, double length_a, double length_b, double length_c,
46                 double thick_rim_a, double thick_rim_b, double thick_rim_c)
47{
48    switch (mode) {
49    default:
50    case 1: // equivalent sphere
51        return radius_from_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
52    case 2: // half outer length a
53        return 0.5 * length_a + thick_rim_a;
54    case 3: // half outer length b
55        return 0.5 * length_b + thick_rim_b;
56    case 4: // half outer length c
57        return 0.5 * length_c + thick_rim_c;
58    case 5: // equivalent circular cross-section
59        return radius_from_crosssection(length_a, length_b, thick_rim_a, thick_rim_b);
60    case 6: // half outer ab diagonal
61        return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b));
62    case 7: // half outer diagonal
63        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));
64    }
65}
66
67static void
68Fq(double q,
69    double *F1,
70    double *F2,
71    double core_sld,
72    double arim_sld,
73    double brim_sld,
74    double crim_sld,
75    double solvent_sld,
76    double length_a,
77    double length_b,
78    double length_c,
79    double thick_rim_a,
80    double thick_rim_b,
81    double thick_rim_c)
82{
83    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c
84    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
85    // Code is rewritten, the code is compliant with Diva Singh's thesis now (Dirk Honecker)
86    // Code rewritten; cross checked against hollow rectangular prism and realspace (PAK)
87
88    const double half_q = 0.5*q;
89
90    const double tA = length_a + 2.0*thick_rim_a;
91    const double tB = length_b + 2.0*thick_rim_b;
92    const double tC = length_c + 2.0*thick_rim_c;
93
94    // Scale factors
95    const double dr0 = (core_sld-solvent_sld);
96    const double drA = (arim_sld-solvent_sld);
97    const double drB = (brim_sld-solvent_sld);
98    const double drC = (crim_sld-solvent_sld);
99
100    // outer integral (with gauss points), integration limits = 0, 1
101    // substitute d_cos_alpha for sin_alpha d_alpha
102    double outer_sum_F1 = 0; //initialize integral
103    double outer_sum_F2 = 0; //initialize integral
104    for( int i=0; i<GAUSS_N; i++) {
105        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 );
106        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha);
107        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q);
108        const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q);
109
110        // inner integral (with gauss points), integration limits = 0, 1
111        // substitute beta = PI/2 u (so 2/PI * d_(PI/2 * beta) = d_beta)
112        double inner_sum_F1 = 0.0;
113        double inner_sum_F2 = 0.0;
114        for(int j=0; j<GAUSS_N; j++) {
115            const double u = 0.5 * ( GAUSS_Z[j] + 1.0 );
116            double sin_beta, cos_beta;
117            SINCOS(M_PI_2*u, sin_beta, cos_beta);
118            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta);
119            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta);
120            const double siAt = tA * sas_sinx_x(tA * mu * sin_beta);
121            const double siBt = tB * sas_sinx_x(tB * mu * cos_beta);
122
123#if OVERLAPPING
124            const double f = dr0*siA*siB*siC
125                + drA*(siAt-siA)*siB*siC
126                + drB*siAt*(siBt-siB)*siC
127                + drC*siAt*siBt*(siCt-siC);
128#else
129            const double f = dr0*siA*siB*siC
130                + drA*(siAt-siA)*siB*siC
131                + drB*siA*(siBt-siB)*siC
132                + drC*siA*siB*(siCt-siC);
133#endif
134
135            inner_sum_F1 += GAUSS_W[j] * f;
136            inner_sum_F2 += GAUSS_W[j] * f * f;
137        }
138        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5
139        // and sum up the outer integral
140        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * 0.5;
141        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * 0.5;
142    }
143    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5
144    outer_sum_F1 *= 0.5;
145    outer_sum_F2 *= 0.5;
146
147    //convert from [1e-12 A-1] to [cm-1]
148    *F1 = 1.0e-2 * outer_sum_F1;
149    *F2 = 1.0e-4 * outer_sum_F2;
150}
151
152static double
153Iqabc(double qa, double qb, double qc,
154    double core_sld,
155    double arim_sld,
156    double brim_sld,
157    double crim_sld,
158    double solvent_sld,
159    double length_a,
160    double length_b,
161    double length_c,
162    double thick_rim_a,
163    double thick_rim_b,
164    double thick_rim_c)
165{
166    // cspkernel in csparallelepiped recoded here
167    const double dr0 = core_sld-solvent_sld;
168    const double drA = arim_sld-solvent_sld;
169    const double drB = brim_sld-solvent_sld;
170    const double drC = crim_sld-solvent_sld;
171
172    const double tA = length_a + 2.0*thick_rim_a;
173    const double tB = length_b + 2.0*thick_rim_b;
174    const double tC = length_c + 2.0*thick_rim_c;
175    const double siA = length_a*sas_sinx_x(0.5*length_a*qa);
176    const double siB = length_b*sas_sinx_x(0.5*length_b*qb);
177    const double siC = length_c*sas_sinx_x(0.5*length_c*qc);
178    const double siAt = tA*sas_sinx_x(0.5*tA*qa);
179    const double siBt = tB*sas_sinx_x(0.5*tB*qb);
180    const double siCt = tC*sas_sinx_x(0.5*tC*qc);
181
182#if OVERLAPPING
183    const double f = dr0*siA*siB*siC
184        + drA*(siAt-siA)*siB*siC
185        + drB*siAt*(siBt-siB)*siC
186        + drC*siAt*siBt*(siCt-siC);
187#else
188    const double f = dr0*siA*siB*siC
189        + drA*(siAt-siA)*siB*siC
190        + drB*siA*(siBt-siB)*siC
191        + drC*siA*siB*(siCt-siC);
192#endif
193
194    return 1.0e-4 * f * f;
195}
Note: See TracBrowser for help on using the repository browser.