1 | static double form_volume( |
---|
2 | int n_shells, |
---|
3 | double radius_core, |
---|
4 | double thick_inter0, |
---|
5 | double thick_flat[], |
---|
6 | double thick_inter[]) |
---|
7 | { |
---|
8 | int i; |
---|
9 | double r = radius_core; |
---|
10 | r += thick_inter0; |
---|
11 | for (i=0; i < n_shells; i++) { |
---|
12 | r += thick_inter[i]; |
---|
13 | r += thick_flat[i]; |
---|
14 | } |
---|
15 | return M_4PI_3*cube(r); |
---|
16 | } |
---|
17 | |
---|
18 | |
---|
19 | static double sphere_sld_kernel( |
---|
20 | double q, |
---|
21 | int n_shells, |
---|
22 | int npts_inter, |
---|
23 | double radius_core, |
---|
24 | double sld_core, |
---|
25 | double sld_solvent, |
---|
26 | double func_inter_core, |
---|
27 | double thick_inter_core, |
---|
28 | double nu_inter_core, |
---|
29 | double sld_flat[], |
---|
30 | double thick_flat[], |
---|
31 | double func_inter[], |
---|
32 | double thick_inter[], |
---|
33 | double nu_inter[] ) { |
---|
34 | |
---|
35 | int i,j,k; |
---|
36 | int n_s; |
---|
37 | |
---|
38 | double sld_i,sld_f,dz,bes,fun,f,vol,qr,r,contr,f2; |
---|
39 | double sign,slope=0.0; |
---|
40 | double pi; |
---|
41 | |
---|
42 | double total_thick=0.0; |
---|
43 | |
---|
44 | //TODO: This part can be further simplified |
---|
45 | int fun_type[12]; |
---|
46 | double sld[12]; |
---|
47 | double thick_internal[12]; |
---|
48 | double thick[12]; |
---|
49 | double fun_coef[12]; |
---|
50 | |
---|
51 | fun_type[0] = func_inter_core; |
---|
52 | fun_coef[0] = fabs(nu_inter_core); |
---|
53 | sld[0] = sld_core; |
---|
54 | thick[0] = radius_core; |
---|
55 | thick_internal[0] = thick_inter_core; |
---|
56 | |
---|
57 | for (i =1; i<=n_shells; i++){ |
---|
58 | sld[i] = sld_flat[i-1]; |
---|
59 | thick_internal[i]= thick_inter[i-1]; |
---|
60 | thick[i] = thick_flat[i-1]; |
---|
61 | fun_type[i] = func_inter[i-1]; |
---|
62 | fun_coef[i] = fabs(nu_inter[i-1]); |
---|
63 | total_thick += thick[i]; |
---|
64 | total_thick += thick_internal[i]; //doesn't account for core layer |
---|
65 | } |
---|
66 | |
---|
67 | sld[n_shells+1] = sld_solvent; |
---|
68 | thick[n_shells+1] = total_thick/5.0; |
---|
69 | thick_internal[n_shells+1] = 0.0; |
---|
70 | fun_coef[n_shells+1] = 0.0; |
---|
71 | fun_type[n_shells+1] = 0; |
---|
72 | |
---|
73 | pi = 4.0*atan(1.0); |
---|
74 | f = 0.0; |
---|
75 | r = 0.0; |
---|
76 | vol = 0.0; |
---|
77 | sld_f = sld_core; |
---|
78 | |
---|
79 | dz = 0.0; |
---|
80 | // iteration for # of shells + core + solvent |
---|
81 | for (i=0;i<=n_shells+1; i++){ |
---|
82 | // iteration for flat and interface |
---|
83 | for (j=0;j<2;j++){ |
---|
84 | // iteration for sub_shells in the interface |
---|
85 | // starts from #1 sub-layer |
---|
86 | for (n_s=1;n_s<=npts_inter; n_s++){ |
---|
87 | // for solvent, it doesn't have an interface |
---|
88 | if (i==n_shells+1 && j==1) |
---|
89 | break; |
---|
90 | // for flat layers |
---|
91 | if (j==0){ |
---|
92 | dz = thick[i]; |
---|
93 | sld_i = sld[i]; |
---|
94 | slope = 0.0; |
---|
95 | } |
---|
96 | // for interfacial sub_shells |
---|
97 | else{ |
---|
98 | dz = thick_internal[i]/npts_inter; |
---|
99 | // find sld_i at the outer boundary of sub-layer #n_s |
---|
100 | sld_i = intersldfunc(fun_type[i], npts_inter, n_s, |
---|
101 | fun_coef[i], sld[i], sld[i+1]); |
---|
102 | // calculate slope |
---|
103 | slope= (sld_i -sld_f)/dz; |
---|
104 | } |
---|
105 | contr = sld_f-slope*r; |
---|
106 | // iteration for the left and right boundary of the shells |
---|
107 | for (k=0; k<2; k++){ |
---|
108 | // At r=0, the contribution to I is zero, so skip it. |
---|
109 | if ( i == 0 && j == 0 && k == 0){ |
---|
110 | continue; |
---|
111 | } |
---|
112 | // On the top of slovent there is no interface; skip it. |
---|
113 | if (i == n_shells+1 && k == 1){ |
---|
114 | continue; |
---|
115 | } |
---|
116 | // At the right side (outer) boundary |
---|
117 | if ( k == 1){ |
---|
118 | sign = 1.0; |
---|
119 | r += dz; |
---|
120 | } |
---|
121 | // At the left side (inner) boundary |
---|
122 | else{ |
---|
123 | sign = -1.0; |
---|
124 | } |
---|
125 | |
---|
126 | qr = q * r; |
---|
127 | fun = 0.0; |
---|
128 | |
---|
129 | if(qr == 0.0){ |
---|
130 | bes = sign * 1.0; |
---|
131 | } |
---|
132 | else{ |
---|
133 | // for flat sub-layer |
---|
134 | //TODO: Single precision calculation fails here |
---|
135 | bes = sign * sph_j1c(qr); |
---|
136 | if (fabs(slope) > 0.0 ){ |
---|
137 | const double qrsq = qr*qr; |
---|
138 | double sinqr, cosqr; |
---|
139 | SINCOS(qr, sinqr, cosqr); |
---|
140 | fun = sign * 3.0 * r * |
---|
141 | (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq * qrsq); |
---|
142 | // In the onioon model Jae-He's formula is rederived |
---|
143 | // and gives following: |
---|
144 | //fun = 3.0 * sign * r * |
---|
145 | //(2.0*cosqr + qr*sinqr)/(qrsq*qrsq); |
---|
146 | //But this seems not to be working in this case... |
---|
147 | } |
---|
148 | } |
---|
149 | |
---|
150 | // update total volume |
---|
151 | vol = M_4PI_3 * cube(r); |
---|
152 | f += vol * (bes * contr + fun * slope); |
---|
153 | } |
---|
154 | sld_f = sld_i; |
---|
155 | // no sub-layer iteration (n_s loop) for the flat layer |
---|
156 | if (j==0) |
---|
157 | break; |
---|
158 | } |
---|
159 | } |
---|
160 | } |
---|
161 | |
---|
162 | f2 = f * f * 1.0e-4; |
---|
163 | return (f2); |
---|
164 | } |
---|
165 | |
---|
166 | |
---|
167 | /** |
---|
168 | * Function to evaluate 1D SphereSLD function |
---|
169 | * @param q: q-value |
---|
170 | * @return: function value |
---|
171 | */ |
---|
172 | static double Iq(double q, |
---|
173 | int n_shells, |
---|
174 | int npts_inter, |
---|
175 | double radius_core, |
---|
176 | double sld_core, |
---|
177 | double sld_solvent, |
---|
178 | double func_inter0, |
---|
179 | double thick_inter0, |
---|
180 | double nu_inter0, |
---|
181 | double sld_flat[], |
---|
182 | double thick_flat[], |
---|
183 | double func_inter[], |
---|
184 | double thick_inter[], |
---|
185 | double nu_inter[] ) { |
---|
186 | |
---|
187 | //printf("Number of points %d\n",npts_inter); |
---|
188 | double intensity; |
---|
189 | //TODO: Remove this container at later stage. |
---|
190 | /*double dp[60]; |
---|
191 | dp[0] = n_shells; |
---|
192 | //This is scale will also have to be removed at some stage |
---|
193 | dp[1] = 1.0; |
---|
194 | dp[2] = thick_inter0; |
---|
195 | dp[3] = func_inter0; |
---|
196 | dp[4] = sld_core; |
---|
197 | dp[5] = sld_solvent; |
---|
198 | dp[6] = 0.0; |
---|
199 | |
---|
200 | for (int i=0; i<n_shells; i++){ |
---|
201 | dp[i+7] = sld_flat[i]; |
---|
202 | dp[i+17] = thick_inter[i]; |
---|
203 | dp[i+27] = thick_flat[i]; |
---|
204 | dp[i+37] = func_inter[i]; |
---|
205 | dp[i+47] = nu_inter[i]; |
---|
206 | } |
---|
207 | |
---|
208 | dp[57] = npts_inter; |
---|
209 | dp[58] = nu_inter0; |
---|
210 | dp[59] = radius_core; |
---|
211 | */ |
---|
212 | intensity = sphere_sld_kernel(q, n_shells, npts_inter, radius_core, |
---|
213 | sld_core, sld_solvent, func_inter0, thick_inter0, nu_inter0, |
---|
214 | sld_flat, thick_flat, func_inter, thick_inter, nu_inter); |
---|
215 | //intensity *=1.0e-4; |
---|
216 | //printf("%10d\n",intensity); |
---|
217 | return intensity; |
---|
218 | } |
---|
219 | |
---|
220 | /** |
---|
221 | * Function to evaluate 2D SphereSLD function |
---|
222 | * @param q_x: value of Q along x |
---|
223 | * @param q_y: value of Q along y |
---|
224 | * @return: function value |
---|
225 | */ |
---|
226 | |
---|
227 | /*static double Iqxy(double qx, double qy, |
---|
228 | int n_shells, |
---|
229 | int npts_inter, |
---|
230 | double radius_core |
---|
231 | double sld_core, |
---|
232 | double sld_solvent, |
---|
233 | double sld_flat[], |
---|
234 | double thick_flat[], |
---|
235 | double func_inter[], |
---|
236 | double thick_inter[], |
---|
237 | double nu_inter[], |
---|
238 | ) { |
---|
239 | |
---|
240 | double q = sqrt(qx*qx + qy*qy); |
---|
241 | return Iq(q, n_shells, npts_inter, radius_core, sld_core, sld_solvent, |
---|
242 | sld_flat[], thick_flat[], func_inter[], thick_inter[], nu_inter[]) |
---|
243 | }*/ |
---|
244 | |
---|