source: sasmodels/sasmodels/models/core_shell_parallelepiped.c @ 14838a3

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

update parallelepiped and core_shell_parallelepiped to use ORIENT macros

  • Property mode set to 100644
File size: 7.9 KB
Line 
1double form_volume(double length_a, double length_b, double length_c, 
2                   double thick_rim_a, double thick_rim_b, double thick_rim_c);
3double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld,
4          double solvent_sld, double length_a, double length_b, double length_c,
5          double thick_rim_a, double thick_rim_b, double thick_rim_c);
6double Iqxy(double qx, double qy, double core_sld, double arim_sld, double brim_sld,
7            double crim_sld, double solvent_sld, double length_a, double length_b,
8            double length_c, double thick_rim_a, double thick_rim_b,
9            double thick_rim_c, double theta, double phi, double psi);
10
11double form_volume(double length_a, double length_b, double length_c, 
12                   double thick_rim_a, double thick_rim_b, double thick_rim_c)
13{
14    //return length_a * length_b * length_c;
15    return length_a * length_b * length_c + 
16           2.0 * thick_rim_a * length_b * length_c + 
17           2.0 * thick_rim_b * length_a * length_c +
18           2.0 * thick_rim_c * length_a * length_b;
19}
20
21double Iq(double q,
22    double core_sld,
23    double arim_sld,
24    double brim_sld,
25    double crim_sld,
26    double solvent_sld,
27    double length_a,
28    double length_b,
29    double length_c,
30    double thick_rim_a,
31    double thick_rim_b,
32    double thick_rim_c)
33{
34    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled
35    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
36   
37    const double mu = 0.5 * q * length_b;
38   
39    //calculate volume before rescaling (in original code, but not used)
40    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);       
41    //double vol = length_a * length_b * length_c +
42    //       2.0 * thick_rim_a * length_b * length_c +
43    //       2.0 * thick_rim_b * length_a * length_c +
44    //       2.0 * thick_rim_c * length_a * length_b;
45   
46    // Scale sides by B
47    const double a_scaled = length_a / length_b;
48    const double c_scaled = length_c / length_b;
49
50    // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG)
51    // the a_scaled in the definition of tb was present in CSPPKernel in libCylinder.c,
52    // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions
53    // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!)
54    double ta = (a_scaled + 2.0*thick_rim_a)/length_b;
55    double tb = (a_scaled + 2.0*thick_rim_b)/length_b;
56
57    double Vin = length_a * length_b * length_c;
58    //double Vot = (length_a * length_b * length_c +
59    //            2.0 * thick_rim_a * length_b * length_c +
60    //            2.0 * length_a * thick_rim_b * length_c +
61    //            2.0 * length_a * length_b * thick_rim_c);
62    double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc)
63    double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
64
65    // Scale factors (note that drC is not used later)
66    const double drho0 = (core_sld-solvent_sld);
67    const double drhoA = (arim_sld-solvent_sld);
68    const double drhoB = (brim_sld-solvent_sld);
69    //const double drC_Vot = (crim_sld-solvent_sld)*Vot;
70
71    // Precompute scale factors for combining cross terms from the shape
72    const double scale23 = drhoA*V1;
73    const double scale14 = drhoB*V2;
74    const double scale12 = drho0*Vin - scale23 - scale14;
75
76    // outer integral (with gauss points), integration limits = 0, 1
77    double outer_total = 0; //initialize integral
78
79    for( int i=0; i<76; i++) {
80        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );
81        double mu_proj = mu * sqrt(1.0-sigma*sigma);
82
83        // inner integral (with gauss points), integration limits = 0, 1
84        double inner_total = 0.0;
85        for(int j=0; j<76; j++) {
86            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
87            double sin_uu, cos_uu;
88            SINCOS(M_PI_2*uu, sin_uu, cos_uu);
89            const double si1 = sinc(mu_proj * sin_uu * a_scaled);
90            const double si2 = sinc(mu_proj * cos_uu);
91            const double si3 = sinc(mu_proj * sin_uu * ta);
92            const double si4 = sinc(mu_proj * cos_uu * tb);
93
94            // Expression in libCylinder.c (neither drC nor Vot are used)
95            const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4;
96
97            // To note also that in csparallelepiped.cpp, there is a function called
98            // cspkernel, which seems to make more sense and has the following comment:
99            //   This expression is different from NIST/IGOR package (I strongly believe the IGOR is wrong!!!). 10/15/2010.
100            //   tmp =( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3)*
101            //   ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors
102            // This is the function called by csparallelepiped_analytical_2D_scaled,
103            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c       
104           
105            //  correct FF : sum of square of phase factors
106            inner_total += Gauss76Wt[j] * form * form;
107        }
108        inner_total *= 0.5;
109
110        // now sum up the outer integral
111        const double si = sinc(mu * c_scaled * sigma);
112        outer_total += Gauss76Wt[i] * inner_total * si * si;
113    }
114    outer_total *= 0.5;
115
116    //convert from [1e-12 A-1] to [cm-1]
117    return 1.0e-4 * outer_total;
118}
119
120double Iqxy(double qx, double qy,
121    double core_sld,
122    double arim_sld,
123    double brim_sld,
124    double crim_sld,
125    double solvent_sld,
126    double length_a,
127    double length_b,
128    double length_c,
129    double thick_rim_a,
130    double thick_rim_b,
131    double thick_rim_c,
132    double theta,
133    double phi,
134    double psi)
135{
136    double q, cos_val_a, cos_val_b, cos_val_c;
137    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a);
138
139    // cspkernel in csparallelepiped recoded here
140    const double dr0 = core_sld-solvent_sld;
141    const double drA = arim_sld-solvent_sld;
142    const double drB = brim_sld-solvent_sld;
143    const double drC = crim_sld-solvent_sld;
144
145    double Vin = length_a * length_b * length_c;
146    double V1 = 2.0 * thick_rim_a * length_b * length_c;    // incorrect V1(aa*bb*cc+2*ta*bb*cc)
147    double V2 = 2.0 * length_a * thick_rim_b * length_c;    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
148    double V3 = 2.0 * length_a * length_b * thick_rim_c;
149    // As for the 1D case, Vot is not used
150    //double Vot = (length_a * length_b * length_c +
151    //              2.0 * thick_rim_a * length_b * length_c +
152    //              2.0 * length_a * thick_rim_b * length_c +
153    //              2.0 * length_a * length_b * thick_rim_c);
154
155    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
156    // the scaling by B. The use of length_a for the 3 of them seems clearly a mistake to me,
157    // but for the moment I let it like this until understanding better the code.
158    double ta = length_a + 2.0*thick_rim_a;
159    double tb = length_a + 2.0*thick_rim_b;
160    double tc = length_a + 2.0*thick_rim_c;
161    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
162    double siA = sinc(0.5*q*length_a*cos_val_a);
163    double siB = sinc(0.5*q*length_b*cos_val_b);
164    double siC = sinc(0.5*q*length_c*cos_val_c);
165    double siAt = sinc(0.5*q*ta*cos_val_a);
166    double siBt = sinc(0.5*q*tb*cos_val_b);
167    double siCt = sinc(0.5*q*tc*cos_val_c);
168   
169
170    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed
171    // in the 1D code, but should be checked!
172    double f = ( dr0*siA*siB*siC*Vin
173               + drA*(siAt-siA)*siB*siC*V1
174               + drB*siA*(siBt-siB)*siC*V2
175               + drC*siA*siB*(siCt*siCt-siC)*V3);
176   
177    return 1.0e-4 * f * f;
178}
Note: See TracBrowser for help on using the repository browser.