source: sasmodels/sasmodels/models/core_shell_parallelepiped.c @ 904cd9c

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

Merge branch 'master' into ticket-776-orientation

  • Property mode set to 100644
File size: 6.5 KB
Line 
1static double
2form_volume(double length_a, double length_b, double length_c,
3    double thick_rim_a, double thick_rim_b, double thick_rim_c)
4{
5    //return length_a * length_b * length_c;
6    return length_a * length_b * length_c +
7           2.0 * thick_rim_a * length_b * length_c +
8           2.0 * thick_rim_b * length_a * length_c +
9           2.0 * thick_rim_c * length_a * length_b;
10}
11
12static double
13Iq(double q,
14    double core_sld,
15    double arim_sld,
16    double brim_sld,
17    double crim_sld,
18    double solvent_sld,
19    double length_a,
20    double length_b,
21    double length_c,
22    double thick_rim_a,
23    double thick_rim_b,
24    double thick_rim_c)
25{
26    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled
27    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
28    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker)
29
30    const double mu = 0.5 * q * length_b;
31
32    //calculate volume before rescaling (in original code, but not used)
33    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
34    //double vol = length_a * length_b * length_c +
35    //       2.0 * thick_rim_a * length_b * length_c +
36    //       2.0 * thick_rim_b * length_a * length_c +
37    //       2.0 * thick_rim_c * length_a * length_b;
38
39    // Scale sides by B
40    const double a_scaled = length_a / length_b;
41    const double c_scaled = length_c / length_b;
42
43    double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b;
44    double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b;
45    double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present
46
47    double Vin = length_a * length_b * length_c;
48    //double Vot = (length_a * length_b * length_c +
49    //            2.0 * thick_rim_a * length_b * length_c +
50    //            2.0 * length_a * thick_rim_b * length_c +
51    //            2.0 * length_a * length_b * thick_rim_c);
52    double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc)
53    double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
54    double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present
55
56    // Scale factors (note that drC is not used later)
57    const double drho0 = (core_sld-solvent_sld);
58    const double drhoA = (arim_sld-solvent_sld);
59    const double drhoB = (brim_sld-solvent_sld);
60    const double drhoC = (crim_sld-solvent_sld);  // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot;
61
62
63    // Precompute scale factors for combining cross terms from the shape
64    const double scale23 = drhoA*V1;
65    const double scale14 = drhoB*V2;
66    const double scale24 = drhoC*V3;
67    const double scale11 = drho0*Vin;
68    const double scale12 = drho0*Vin - scale23 - scale14 - scale24;
69
70    // outer integral (with gauss points), integration limits = 0, 1
71    double outer_total = 0; //initialize integral
72
73    for( int i=0; i<76; i++) {
74        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );
75        double mu_proj = mu * sqrt(1.0-sigma*sigma);
76
77        // inner integral (with gauss points), integration limits = 0, 1
78        double inner_total = 0.0;
79        double inner_total_crim = 0.0;
80        for(int j=0; j<76; j++) {
81            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
82            double sin_uu, cos_uu;
83            SINCOS(M_PI_2*uu, sin_uu, cos_uu);
84            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
85            const double si2 = sas_sinx_x(mu_proj * cos_uu );
86            const double si3 = sas_sinx_x(mu_proj * sin_uu * ta);
87            const double si4 = sas_sinx_x(mu_proj * cos_uu * tb);
88
89            // Expression in libCylinder.c (neither drC nor Vot are used)
90            const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4;
91            const double form_crim = scale11*si1*si2;
92
93            //  correct FF : sum of square of phase factors
94            inner_total += Gauss76Wt[j] * form * form;
95            inner_total_crim += Gauss76Wt[j] * form_crim * form_crim;
96        }
97        inner_total *= 0.5;
98        inner_total_crim *= 0.5;
99        // now sum up the outer integral
100        const double si = sas_sinx_x(mu * c_scaled * sigma);
101        const double si_crim = sas_sinx_x(mu * tc * sigma);
102        outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim);
103    }
104    outer_total *= 0.5;
105
106    //convert from [1e-12 A-1] to [cm-1]
107    return 1.0e-4 * outer_total;
108}
109
110static double
111Iqxy(double qa, double qb, double qc,
112    double core_sld,
113    double arim_sld,
114    double brim_sld,
115    double crim_sld,
116    double solvent_sld,
117    double length_a,
118    double length_b,
119    double length_c,
120    double thick_rim_a,
121    double thick_rim_b,
122    double thick_rim_c)
123{
124    // cspkernel in csparallelepiped recoded here
125    const double dr0 = core_sld-solvent_sld;
126    const double drA = arim_sld-solvent_sld;
127    const double drB = brim_sld-solvent_sld;
128    const double drC = crim_sld-solvent_sld;
129
130    double Vin = length_a * length_b * length_c;
131    double V1 = 2.0 * thick_rim_a * length_b * length_c;    // incorrect V1(aa*bb*cc+2*ta*bb*cc)
132    double V2 = 2.0 * length_a * thick_rim_b * length_c;    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
133    double V3 = 2.0 * length_a * length_b * thick_rim_c;
134    // As for the 1D case, Vot is not used
135    //double Vot = (length_a * length_b * length_c +
136    //              2.0 * thick_rim_a * length_b * length_c +
137    //              2.0 * length_a * thick_rim_b * length_c +
138    //              2.0 * length_a * length_b * thick_rim_c);
139
140    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
141    // the scaling by B.
142    double ta = length_a + 2.0*thick_rim_a;
143    double tb = length_b + 2.0*thick_rim_b;
144    double tc = length_c + 2.0*thick_rim_c;
145    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
146    double siA = sas_sinx_x(0.5*length_a*qa);
147    double siB = sas_sinx_x(0.5*length_b*qb);
148    double siC = sas_sinx_x(0.5*length_c*qc);
149    double siAt = sas_sinx_x(0.5*ta*qa);
150    double siBt = sas_sinx_x(0.5*tb*qb);
151    double siCt = sas_sinx_x(0.5*tc*qc);
152
153
154    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed
155    // in the 1D code, but should be checked!
156    double f = ( dr0*siA*siB*siC*Vin
157               + drA*(siAt-siA)*siB*siC*V1
158               + drB*siA*(siBt-siB)*siC*V2
159               + drC*siA*siB*(siCt-siC)*V3);
160
161    return 1.0e-4 * f * f;
162}
Note: See TracBrowser for help on using the repository browser.