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

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since fc0b7aa was fc0b7aa, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

core_shell_parallelepiped: modify calculation to be consistent with hollow rectangular prism

  • Property mode set to 100644
File size: 6.7 KB
RevLine 
[fc0b7aa]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
[becded3]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)
[44bd2be]10{
[fc0b7aa]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
[44bd2be]27}
28
[becded3]29static double
30Iq(double q,
[44bd2be]31    double core_sld,
32    double arim_sld,
33    double brim_sld,
34    double crim_sld,
35    double solvent_sld,
[2222134]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)
[44bd2be]42{
[3a1fc7d]43    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c
[44bd2be]44    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
[c69d6d6]45    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker)
[2a0b2b1]46
[14838a3]47    const double mu = 0.5 * q * length_b;
[2a0b2b1]48
[44bd2be]49    // Scale sides by B
[3a1fc7d]50    const double a_over_b = length_a / length_b;
51    const double c_over_b = length_c / length_b;
[14838a3]52
[3a1fc7d]53    double tA_over_b = a_over_b + 2.0*thick_rim_a/length_b;
54    double tB_over_b = 1+ 2.0*thick_rim_b/length_b;
55    double tC_over_b = c_over_b + 2.0*thick_rim_c/length_b;
[14838a3]56
57    double Vin = length_a * length_b * length_c;
[fc0b7aa]58#if OVERLAPPING
59    const double capA_area = length_b*length_c;
60    const double capB_area = (length_a+2.*thick_rim_a)*length_c;
61    const double capC_area = (length_a+2.*thick_rim_a)*(length_b+2.*thick_rim_b);
62#else
63    const double capA_area = length_b*length_c;
64    const double capB_area = length_a*length_c;
65    const double capC_area = length_a*length_b;
66#endif
67    const double Va = length_a * capA_area;
68    const double Vb = length_b * capB_area;
69    const double Vc = length_c * capC_area;
70    const double Vat = Va + 2.0 * thick_rim_a * capA_area;
71    const double Vbt = Vb + 2.0 * thick_rim_b * capB_area;
72    const double Vct = Vc + 2.0 * thick_rim_c * capC_area;
[14838a3]73
74    // Scale factors (note that drC is not used later)
[3a1fc7d]75    const double dr0 = (core_sld-solvent_sld);
76    const double drA = (arim_sld-solvent_sld);
77    const double drB = (brim_sld-solvent_sld);
78    const double drC = (crim_sld-solvent_sld);
[14838a3]79
80    // outer integral (with gauss points), integration limits = 0, 1
[3a1fc7d]81    double outer_sum = 0; //initialize integral
[14838a3]82    for( int i=0; i<76; i++) {
83        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );
84        double mu_proj = mu * sqrt(1.0-sigma*sigma);
85
[3a1fc7d]86        // inner integral (with gauss points), integration limits = 0, pi/2
87        const double siC = sas_sinx_x(mu * sigma * c_over_b);
88        const double siCt = sas_sinx_x(mu * sigma * tC_over_b);
89        double inner_sum = 0.0;
[14838a3]90        for(int j=0; j<76; j++) {
91            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
92            double sin_uu, cos_uu;
93            SINCOS(M_PI_2*uu, sin_uu, cos_uu);
[3a1fc7d]94            const double siA = sas_sinx_x(mu_proj * sin_uu * a_over_b);
95            const double siB = sas_sinx_x(mu_proj * cos_uu );
96            const double siAt = sas_sinx_x(mu_proj * sin_uu * tA_over_b);
97            const double siBt = sas_sinx_x(mu_proj * cos_uu * tB_over_b);
98
[fc0b7aa]99#if OVERLAPPING
100            const double f = dr0*Vin*siA*siB*siC
101                + drA*(Vat*siAt-Va*siA)*siB*siC
102                + drB*siAt*(Vbt*siBt-Vb*siB)*siC
103                + drC*siAt*siBt*(Vct*siCt-Vc*siC);
104#else
105            const double f = dr0*Vin*siA*siB*siC
106                + drA*(Vat*siAt-Va*siA)*siB*siC
107                + drB*siA*(Vbt*siBt-Vb*siB)*siC
108                + drC*siA*siB*(Vct*siCt-Vc*siC);
109#endif
110
111            inner_sum += Gauss76Wt[j] * f * f;
[44bd2be]112        }
[3a1fc7d]113        inner_sum *= 0.5;
[14838a3]114        // now sum up the outer integral
[3a1fc7d]115        outer_sum += Gauss76Wt[i] * inner_sum;
[44bd2be]116    }
[3a1fc7d]117    outer_sum *= 0.5;
[44bd2be]118
[14838a3]119    //convert from [1e-12 A-1] to [cm-1]
[3a1fc7d]120    return 1.0e-4 * outer_sum;
[44bd2be]121}
122
[becded3]123static double
124Iqxy(double qa, double qb, double qc,
[44bd2be]125    double core_sld,
126    double arim_sld,
127    double brim_sld,
128    double crim_sld,
129    double solvent_sld,
[2222134]130    double length_a,
131    double length_b,
132    double length_c,
133    double thick_rim_a,
134    double thick_rim_b,
[becded3]135    double thick_rim_c)
[44bd2be]136{
[14838a3]137    // cspkernel in csparallelepiped recoded here
138    const double dr0 = core_sld-solvent_sld;
139    const double drA = arim_sld-solvent_sld;
140    const double drB = brim_sld-solvent_sld;
141    const double drC = crim_sld-solvent_sld;
142
143    double Vin = length_a * length_b * length_c;
[fc0b7aa]144#if OVERLAPPING
145    const double capA_area = length_b*length_c;
146    const double capB_area = (length_a+2.*thick_rim_a)*length_c;
147    const double capC_area = (length_a+2.*thick_rim_a)*(length_b+2.*thick_rim_b);
148#else
149    const double capA_area = length_b*length_c;
150    const double capB_area = length_a*length_c;
151    const double capC_area = length_a*length_b;
152#endif
153    const double Va = length_a * capA_area;
154    const double Vb = length_b * capB_area;
155    const double Vc = length_c * capC_area;
156    const double Vat = Va + 2.0 * thick_rim_a * capA_area;
157    const double Vbt = Vb + 2.0 * thick_rim_b * capB_area;
158    const double Vct = Vc + 2.0 * thick_rim_c * capC_area;
[14838a3]159
[44bd2be]160    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
[c69d6d6]161    // the scaling by B.
[fc0b7aa]162    const double tA = length_a + 2.0*thick_rim_a;
163    const double tB = length_b + 2.0*thick_rim_b;
164    const double tC = length_c + 2.0*thick_rim_c;
165    const double siA = sas_sinx_x(0.5*length_a*qa);
166    const double siB = sas_sinx_x(0.5*length_b*qb);
167    const double siC = sas_sinx_x(0.5*length_c*qc);
168    const double siAt = sas_sinx_x(0.5*tA*qa);
169    const double siBt = sas_sinx_x(0.5*tB*qb);
170    const double siCt = sas_sinx_x(0.5*tC*qc);
171
172#if OVERLAPPING
173    const double f = dr0*Vin*siA*siB*siC
174        + drA*(Vat*siAt-Va*siA)*siB*siC
175        + drB*siAt*(Vbt*siBt-Vb*siB)*siC
176        + drC*siAt*siBt*(Vct*siCt-Vc*siC);
177#else
178    const double f = dr0*Vin*siA*siB*siC
179        + drA*(Vat*siAt-Va*siA)*siB*siC
180        + drB*siA*(Vbt*siBt-Vb*siB)*siC
181        + drC*siA*siB*(Vct*siCt-Vc*siC);
182#endif
[2a0b2b1]183
[44bd2be]184    return 1.0e-4 * f * f;
185}
Note: See TracBrowser for help on using the repository browser.