1 | static double |
---|
2 | form_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 | 2.0 * thick_rim_a * length_b * length_c + |
---|
7 | 2.0 * thick_rim_b * length_a * length_c + |
---|
8 | 2.0 * thick_rim_c * length_a * length_b; |
---|
9 | } |
---|
10 | |
---|
11 | static double |
---|
12 | Iq(double q, |
---|
13 | double core_sld, |
---|
14 | double arim_sld, |
---|
15 | double brim_sld, |
---|
16 | double crim_sld, |
---|
17 | double solvent_sld, |
---|
18 | double length_a, |
---|
19 | double length_b, |
---|
20 | double length_c, |
---|
21 | double thick_rim_a, |
---|
22 | double thick_rim_b, |
---|
23 | double thick_rim_c) |
---|
24 | { |
---|
25 | // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c |
---|
26 | // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) |
---|
27 | //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) |
---|
28 | |
---|
29 | const double mu = 0.5 * q * length_b; |
---|
30 | |
---|
31 | // Scale sides by B |
---|
32 | const double a_over_b = length_a / length_b; |
---|
33 | const double c_over_b = length_c / length_b; |
---|
34 | |
---|
35 | double tA_over_b = a_over_b + 2.0*thick_rim_a/length_b; |
---|
36 | double tB_over_b = 1+ 2.0*thick_rim_b/length_b; |
---|
37 | double tC_over_b = c_over_b + 2.0*thick_rim_c/length_b; |
---|
38 | |
---|
39 | double Vin = length_a * length_b * length_c; |
---|
40 | double VtA = (2.0 * thick_rim_a * length_b * length_c); |
---|
41 | double VtB = (2.0 * length_a * thick_rim_b * length_c); |
---|
42 | double VtC = (2.0 * length_a * length_b * thick_rim_c); |
---|
43 | |
---|
44 | // Scale factors (note that drC is not used later) |
---|
45 | const double dr0 = (core_sld-solvent_sld); |
---|
46 | const double drA = (arim_sld-solvent_sld); |
---|
47 | const double drB = (brim_sld-solvent_sld); |
---|
48 | const double drC = (crim_sld-solvent_sld); |
---|
49 | |
---|
50 | // Precompute scale factors for combining cross terms from the shape |
---|
51 | const double dr0_Vin = dr0*Vin; |
---|
52 | const double drA_VtA = drA*VtA; |
---|
53 | const double drB_VtB = drB*VtB; |
---|
54 | const double drC_VtC = drC*VtC; |
---|
55 | const double drV_delta = dr0_Vin - drA_VtA - drB_VtB - drC_VtC; |
---|
56 | |
---|
57 | /* *************** algorithm description ****************** |
---|
58 | |
---|
59 | // Rewrite f as x*siC + y*siCt to move the siC/siCt calculation out |
---|
60 | // of the inner loop. That is: |
---|
61 | |
---|
62 | f = (di-da-db-dc) sa sb sc + da sa' sb sc + db sa sb' sc + dc sa sb sc' |
---|
63 | = [ (di-da-db-dc) sa sb + da sa' sb + db sa sb' ] sc + [dc sa sb] sc' |
---|
64 | = x sc + y sc' |
---|
65 | |
---|
66 | // where: |
---|
67 | di = delta rho_core V_core |
---|
68 | da = delta rho_rimA V_rimA |
---|
69 | db = delta rho_rimB V_rimB |
---|
70 | dc = delta rho_rimC V_rimC |
---|
71 | sa = j_0 (q_a a/2) // siA the code |
---|
72 | sb = j_0 (q_b b/2) |
---|
73 | sc = j_0 (q_c c/2) |
---|
74 | sa' = j_0(q_a a_rim/2) // siAt the code |
---|
75 | sb' = j_0(q_b b_rim/2) |
---|
76 | sc' = j_0(q_c c_rim/2) |
---|
77 | |
---|
78 | // qa, qb, and qc are generated using polar coordinates, with the |
---|
79 | // outer loop integrating over [0,1] after the u-substitution |
---|
80 | // sigma = cos(theta), sqrt(1-sigma^2) = sin(theta) |
---|
81 | // and inner loop integrating over [0, pi/2] as |
---|
82 | // uu = phi |
---|
83 | |
---|
84 | ************************************************************ */ |
---|
85 | |
---|
86 | // outer integral (with gauss points), integration limits = 0, 1 |
---|
87 | double outer_sum = 0; //initialize integral |
---|
88 | for( int i=0; i<76; i++) { |
---|
89 | double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); |
---|
90 | double mu_proj = mu * sqrt(1.0-sigma*sigma); |
---|
91 | |
---|
92 | // inner integral (with gauss points), integration limits = 0, pi/2 |
---|
93 | const double siC = sas_sinx_x(mu * sigma * c_over_b); |
---|
94 | const double siCt = sas_sinx_x(mu * sigma * tC_over_b); |
---|
95 | double inner_sum = 0.0; |
---|
96 | for(int j=0; j<76; j++) { |
---|
97 | const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); |
---|
98 | double sin_uu, cos_uu; |
---|
99 | SINCOS(M_PI_2*uu, sin_uu, cos_uu); |
---|
100 | const double siA = sas_sinx_x(mu_proj * sin_uu * a_over_b); |
---|
101 | const double siB = sas_sinx_x(mu_proj * cos_uu ); |
---|
102 | const double siAt = sas_sinx_x(mu_proj * sin_uu * tA_over_b); |
---|
103 | const double siBt = sas_sinx_x(mu_proj * cos_uu * tB_over_b); |
---|
104 | |
---|
105 | const double x = drV_delta*siA*siB + drA_VtA*siB*siAt + drB_VtB*siA*siBt; |
---|
106 | const double form = x*siC + drC_VtC*siA*siB*siCt; |
---|
107 | |
---|
108 | inner_sum += Gauss76Wt[j] * form * form; |
---|
109 | } |
---|
110 | inner_sum *= 0.5; |
---|
111 | // now sum up the outer integral |
---|
112 | outer_sum += Gauss76Wt[i] * inner_sum; |
---|
113 | } |
---|
114 | outer_sum *= 0.5; |
---|
115 | |
---|
116 | //convert from [1e-12 A-1] to [cm-1] |
---|
117 | return 1.0e-4 * outer_sum; |
---|
118 | } |
---|
119 | |
---|
120 | static double |
---|
121 | Iqxy(double qa, double qb, double qc, |
---|
122 | double core_sld, |
---|
123 | double arim_sld, |
---|
124 | double brim_sld, |
---|
125 | double crim_sld, |
---|
126 | double solvent_sld, |
---|
127 | double length_a, |
---|
128 | double length_b, |
---|
129 | double length_c, |
---|
130 | double thick_rim_a, |
---|
131 | double thick_rim_b, |
---|
132 | double thick_rim_c) |
---|
133 | { |
---|
134 | // cspkernel in csparallelepiped recoded here |
---|
135 | const double dr0 = core_sld-solvent_sld; |
---|
136 | const double drA = arim_sld-solvent_sld; |
---|
137 | const double drB = brim_sld-solvent_sld; |
---|
138 | const double drC = crim_sld-solvent_sld; |
---|
139 | |
---|
140 | double Vin = length_a * length_b * length_c; |
---|
141 | double VtA = 2.0 * thick_rim_a * length_b * length_c; |
---|
142 | double VtB = 2.0 * length_a * thick_rim_b * length_c; |
---|
143 | double VtC = 2.0 * length_a * length_b * thick_rim_c; |
---|
144 | |
---|
145 | // The definitions of ta, tb, tc are not the same as in the 1D case because there is no |
---|
146 | // the scaling by B. |
---|
147 | double tA = length_a + 2.0*thick_rim_a; |
---|
148 | double tB = length_b + 2.0*thick_rim_b; |
---|
149 | double tC = length_c + 2.0*thick_rim_c; |
---|
150 | //handle arg=0 separately, as sin(t)/t -> 1 as t->0 |
---|
151 | double siA = sas_sinx_x(0.5*length_a*qa); |
---|
152 | double siB = sas_sinx_x(0.5*length_b*qb); |
---|
153 | double siC = sas_sinx_x(0.5*length_c*qc); |
---|
154 | double siAt = sas_sinx_x(0.5*tA*qa); |
---|
155 | double siBt = sas_sinx_x(0.5*tB*qb); |
---|
156 | double siCt = sas_sinx_x(0.5*tC*qc); |
---|
157 | |
---|
158 | |
---|
159 | // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed |
---|
160 | // in the 1D code, but should be checked! |
---|
161 | double f = ( dr0*Vin*siA*siB*siC |
---|
162 | + drA*VtA*(siAt-siA)*siB*siC |
---|
163 | + drB*VtB*siA*(siBt-siB)*siC |
---|
164 | + drC*VtC*siA*siB*(siCt-siC)); |
---|
165 | |
---|
166 | return 1.0e-4 * f * f; |
---|
167 | } |
---|