1 | double form_volume(double length_a, double length_b, double length_c, |
---|
2 | double thick_rim_a, double thick_rim_b, double thick_rim_c); |
---|
3 | double 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); |
---|
6 | double 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 | |
---|
11 | double 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 | |
---|
21 | double 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 = sas_sinx_x(mu_proj * sin_uu * a_scaled); |
---|
90 | const double si2 = sas_sinx_x(mu_proj * cos_uu); |
---|
91 | const double si3 = sas_sinx_x(mu_proj * sin_uu * ta); |
---|
92 | const double si4 = sas_sinx_x(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 = sas_sinx_x(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 | |
---|
120 | double 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, zhat, yhat, xhat; |
---|
137 | ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); |
---|
138 | const double qa = q*xhat; |
---|
139 | const double qb = q*yhat; |
---|
140 | const double qc = q*zhat; |
---|
141 | |
---|
142 | // cspkernel in csparallelepiped recoded here |
---|
143 | const double dr0 = core_sld-solvent_sld; |
---|
144 | const double drA = arim_sld-solvent_sld; |
---|
145 | const double drB = brim_sld-solvent_sld; |
---|
146 | const double drC = crim_sld-solvent_sld; |
---|
147 | |
---|
148 | double Vin = length_a * length_b * length_c; |
---|
149 | double V1 = 2.0 * thick_rim_a * length_b * length_c; // incorrect V1(aa*bb*cc+2*ta*bb*cc) |
---|
150 | double V2 = 2.0 * length_a * thick_rim_b * length_c; // incorrect V2(aa*bb*cc+2*aa*tb*cc) |
---|
151 | double V3 = 2.0 * length_a * length_b * thick_rim_c; |
---|
152 | // As for the 1D case, Vot is not used |
---|
153 | //double Vot = (length_a * length_b * length_c + |
---|
154 | // 2.0 * thick_rim_a * length_b * length_c + |
---|
155 | // 2.0 * length_a * thick_rim_b * length_c + |
---|
156 | // 2.0 * length_a * length_b * thick_rim_c); |
---|
157 | |
---|
158 | // The definitions of ta, tb, tc are not the same as in the 1D case because there is no |
---|
159 | // the scaling by B. The use of length_a for the 3 of them seems clearly a mistake to me, |
---|
160 | // but for the moment I let it like this until understanding better the code. |
---|
161 | double ta = length_a + 2.0*thick_rim_a; |
---|
162 | double tb = length_a + 2.0*thick_rim_b; |
---|
163 | double tc = length_a + 2.0*thick_rim_c; |
---|
164 | //handle arg=0 separately, as sin(t)/t -> 1 as t->0 |
---|
165 | double siA = sas_sinx_x(0.5*length_a*qa); |
---|
166 | double siB = sas_sinx_x(0.5*length_b*qb); |
---|
167 | double siC = sas_sinx_x(0.5*length_c*qc); |
---|
168 | double siAt = sas_sinx_x(0.5*ta*qa); |
---|
169 | double siBt = sas_sinx_x(0.5*tb*qb); |
---|
170 | double siCt = sas_sinx_x(0.5*tc*qc); |
---|
171 | |
---|
172 | |
---|
173 | // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed |
---|
174 | // in the 1D code, but should be checked! |
---|
175 | double f = ( dr0*siA*siB*siC*Vin |
---|
176 | + drA*(siAt-siA)*siB*siC*V1 |
---|
177 | + drB*siA*(siBt-siB)*siC*V2 |
---|
178 | + drC*siA*siB*(siCt*siCt-siC)*V3); |
---|
179 | |
---|
180 | return 1.0e-4 * f * f; |
---|
181 | } |
---|