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

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since dbf1a60 was dbf1a60, checked in by butler, 16 months ago

Add comments to c code and clean documentatin

Added comments to c code in both parallelepiped and core shell
parallelepiped noting the change of integration varialbes in the
computation. Cleaned up and final corrections to the core shell
documentation and did some cleaning of parallelipiped. In particular
tried to bring a bit more consistency between the docs.

addresses #896

  • Property mode set to 100644
File size: 5.3 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
30Iq(double q,
31    double core_sld,
32    double arim_sld,
33    double brim_sld,
34    double crim_sld,
35    double solvent_sld,
36    double length_a,
37    double length_b,
38    double length_c,
39    double thick_rim_a,
40    double thick_rim_b,
41    double thick_rim_c)
42{
43    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c
44    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
45    // Code is rewritten, the code is compliant with Diva Singh's thesis now (Dirk Honecker)
46    // Code rewritten; cross checked against hollow rectangular prism and realspace (PAK)
47
48    const double half_q = 0.5*q;
49
50    const double tA = length_a + 2.0*thick_rim_a;
51    const double tB = length_b + 2.0*thick_rim_b;
52    const double tC = length_c + 2.0*thick_rim_c;
53
54    // Scale factors
55    const double dr0 = (core_sld-solvent_sld);
56    const double drA = (arim_sld-solvent_sld);
57    const double drB = (brim_sld-solvent_sld);
58    const double drC = (crim_sld-solvent_sld);
59
60    // outer integral (with gauss points), integration limits = 0, 1
61    // substitute d_cos_alpha for sin_alpha d_alpha
62    double outer_sum = 0; //initialize integral
63    for( int i=0; i<GAUSS_N; i++) {
64        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 );
65        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha);
66        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q);
67        const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q);
68
69        // inner integral (with gauss points), integration limits = 0, 1
70        // substitute beta = PI/2 u (so 2/PI * d_(PI/2 * beta) = d_beta)
71        double inner_sum = 0.0;
72        for(int j=0; j<GAUSS_N; j++) {
73            const double u = 0.5 * ( GAUSS_Z[j] + 1.0 );
74            double sin_beta, cos_beta;
75            SINCOS(M_PI_2*u, sin_beta, cos_beta);
76            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta);
77            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta);
78            const double siAt = tA * sas_sinx_x(tA * mu * sin_beta);
79            const double siBt = tB * sas_sinx_x(tB * mu * cos_beta);
80
81#if OVERLAPPING
82            const double f = dr0*siA*siB*siC
83                + drA*(siAt-siA)*siB*siC
84                + drB*siAt*(siBt-siB)*siC
85                + drC*siAt*siBt*(siCt-siC);
86#else
87            const double f = dr0*siA*siB*siC
88                + drA*(siAt-siA)*siB*siC
89                + drB*siA*(siBt-siB)*siC
90                + drC*siA*siB*(siCt-siC);
91#endif
92
93            inner_sum += GAUSS_W[j] * f * f;
94        }
95        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5
96        inner_sum *= 0.5;
97        // now sum up the outer integral
98        outer_sum += GAUSS_W[i] * inner_sum;
99    }
100    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5
101    outer_sum *= 0.5;
102
103    //convert from [1e-12 A-1] to [cm-1]
104    return 1.0e-4 * outer_sum;
105}
106
107static double
108Iqabc(double qa, double qb, double qc,
109    double core_sld,
110    double arim_sld,
111    double brim_sld,
112    double crim_sld,
113    double solvent_sld,
114    double length_a,
115    double length_b,
116    double length_c,
117    double thick_rim_a,
118    double thick_rim_b,
119    double thick_rim_c)
120{
121    // cspkernel in csparallelepiped recoded here
122    const double dr0 = core_sld-solvent_sld;
123    const double drA = arim_sld-solvent_sld;
124    const double drB = brim_sld-solvent_sld;
125    const double drC = crim_sld-solvent_sld;
126
127    const double tA = length_a + 2.0*thick_rim_a;
128    const double tB = length_b + 2.0*thick_rim_b;
129    const double tC = length_c + 2.0*thick_rim_c;
130    const double siA = length_a*sas_sinx_x(0.5*length_a*qa);
131    const double siB = length_b*sas_sinx_x(0.5*length_b*qb);
132    const double siC = length_c*sas_sinx_x(0.5*length_c*qc);
133    const double siAt = tA*sas_sinx_x(0.5*tA*qa);
134    const double siBt = tB*sas_sinx_x(0.5*tB*qb);
135    const double siCt = tC*sas_sinx_x(0.5*tC*qc);
136
137#if OVERLAPPING
138    const double f = dr0*siA*siB*siC
139        + drA*(siAt-siA)*siB*siC
140        + drB*siAt*(siBt-siB)*siC
141        + drC*siAt*siBt*(siCt-siC);
142#else
143    const double f = dr0*siA*siB*siC
144        + drA*(siAt-siA)*siB*siC
145        + drB*siA*(siBt-siB)*siC
146        + drC*siA*siB*(siCt-siC);
147#endif
148
149    return 1.0e-4 * f * f;
150}
Note: See TracBrowser for help on using the repository browser.