1 | /* CylinderFit.c |
---|
2 | |
---|
3 | A simplified project designed to act as a template for your curve fitting function. |
---|
4 | The fitting function is a Cylinder form factor. No resolution effects are included (yet) |
---|
5 | */ |
---|
6 | |
---|
7 | #include "StandardHeaders.h" // Include ANSI headers, Mac headers |
---|
8 | #include "GaussWeights.h" |
---|
9 | #include "libCylinder.h" |
---|
10 | /* CylinderForm : calculates the form factor of a cylinder at the give x-value p->x |
---|
11 | |
---|
12 | Warning: |
---|
13 | The call to WaveData() below returns a pointer to the middle |
---|
14 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
15 | calculations could cause memory to move, you should copy the coefficient |
---|
16 | values to local variables or an array before such operations. |
---|
17 | */ |
---|
18 | double |
---|
19 | CylinderForm(double dp[], double q) |
---|
20 | { |
---|
21 | int i; |
---|
22 | double Pi; |
---|
23 | double scale,radius,length,delrho,bkg,halfheight; //local variables of coefficient wave |
---|
24 | int nord=76; //order of integration |
---|
25 | double uplim,lolim; //upper and lower integration limits |
---|
26 | double summ,zi,yyy,answer,vcyl; //running tally of integration |
---|
27 | |
---|
28 | Pi = 4.0*atan(1.0); |
---|
29 | lolim = 0; |
---|
30 | uplim = Pi/2.0; |
---|
31 | |
---|
32 | summ = 0.0; //initialize intergral |
---|
33 | |
---|
34 | scale = dp[0]; //make local copies in case memory moves |
---|
35 | radius = dp[1]; |
---|
36 | length = dp[2]; |
---|
37 | delrho = dp[3]; |
---|
38 | bkg = dp[4]; |
---|
39 | halfheight = length/2.0; |
---|
40 | for(i=0;i<nord;i++) { |
---|
41 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
42 | yyy = Gauss76Wt[i] * CylKernel(q, radius, halfheight, zi); |
---|
43 | summ += yyy; |
---|
44 | } |
---|
45 | |
---|
46 | answer = (uplim-lolim)/2.0*summ; |
---|
47 | // Multiply by contrast^2 |
---|
48 | answer *= delrho*delrho; |
---|
49 | //normalize by cylinder volume |
---|
50 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
51 | vcyl=Pi*radius*radius*length; |
---|
52 | answer *= vcyl; |
---|
53 | //convert to [cm-1] |
---|
54 | answer *= 1.0e8; |
---|
55 | //Scale |
---|
56 | answer *= scale; |
---|
57 | // add in the background |
---|
58 | answer += bkg; |
---|
59 | |
---|
60 | return answer; |
---|
61 | } |
---|
62 | |
---|
63 | /* EllipCyl76X : calculates the form factor of a elliptical cylinder at the given x-value p->x |
---|
64 | |
---|
65 | Uses 76 pt Gaussian quadrature for both integrals |
---|
66 | |
---|
67 | Warning: |
---|
68 | The call to WaveData() below returns a pointer to the middle |
---|
69 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
70 | calculations could cause memory to move, you should copy the coefficient |
---|
71 | values to local variables or an array before such operations. |
---|
72 | */ |
---|
73 | double |
---|
74 | EllipCyl76(double dp[], double q) |
---|
75 | { |
---|
76 | int i,j; |
---|
77 | double Pi; |
---|
78 | double scale,ra,nu,length,delrho,bkg; //local variables of coefficient wave |
---|
79 | int nord=76; //order of integration |
---|
80 | double va,vb; //upper and lower integration limits |
---|
81 | double summ,zi,yyy,answer,vell; //running tally of integration |
---|
82 | double summj,vaj,vbj,zij,arg, si; //for the inner integration |
---|
83 | |
---|
84 | Pi = 4.0*atan(1.0); |
---|
85 | va = 0; |
---|
86 | vb = 1; //orintational average, outer integral |
---|
87 | vaj=0; |
---|
88 | vbj=Pi; //endpoints of inner integral |
---|
89 | |
---|
90 | summ = 0.0; //initialize intergral |
---|
91 | |
---|
92 | scale = dp[0]; //make local copies in case memory moves |
---|
93 | ra = dp[1]; |
---|
94 | nu = dp[2]; |
---|
95 | length = dp[3]; |
---|
96 | delrho = dp[4]; |
---|
97 | bkg = dp[5]; |
---|
98 | for(i=0;i<nord;i++) { |
---|
99 | //setup inner integral over the ellipsoidal cross-section |
---|
100 | summj=0; |
---|
101 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy |
---|
102 | arg = ra*sqrt(1-zi*zi); |
---|
103 | for(j=0;j<nord;j++) { |
---|
104 | //76 gauss points for the inner integral as well |
---|
105 | zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "y" dummy |
---|
106 | yyy = Gauss76Wt[j] * EllipCylKernel(q,arg,nu,zij); |
---|
107 | summj += yyy; |
---|
108 | } |
---|
109 | //now calculate the value of the inner integral |
---|
110 | answer = (vbj-vaj)/2.0*summj; |
---|
111 | //divide integral by Pi |
---|
112 | answer /=Pi; |
---|
113 | |
---|
114 | //now calculate outer integral |
---|
115 | arg = q*length*zi/2; |
---|
116 | if (arg == 0){ |
---|
117 | si = 1; |
---|
118 | }else{ |
---|
119 | si = sin(arg) * sin(arg) / arg / arg; |
---|
120 | } |
---|
121 | yyy = Gauss76Wt[i] * answer * si; |
---|
122 | summ += yyy; |
---|
123 | } |
---|
124 | answer = (vb-va)/2.0*summ; |
---|
125 | // Multiply by contrast^2 |
---|
126 | answer *= delrho*delrho; |
---|
127 | //normalize by cylinder volume |
---|
128 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
129 | vell = Pi*ra*(nu*ra)*length; |
---|
130 | answer *= vell; |
---|
131 | //convert to [cm-1] |
---|
132 | answer *= 1.0e8; |
---|
133 | //Scale |
---|
134 | answer *= scale; |
---|
135 | // add in the background |
---|
136 | answer += bkg; |
---|
137 | |
---|
138 | return answer; |
---|
139 | } |
---|
140 | |
---|
141 | /* EllipCyl20X : calculates the form factor of a elliptical cylinder at the given x-value p->x |
---|
142 | |
---|
143 | Uses 76 pt Gaussian quadrature for orientational integral |
---|
144 | Uses 20 pt quadrature for the inner integral over the elliptical cross-section |
---|
145 | |
---|
146 | Warning: |
---|
147 | The call to WaveData() below returns a pointer to the middle |
---|
148 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
149 | calculations could cause memory to move, you should copy the coefficient |
---|
150 | values to local variables or an array before such operations. |
---|
151 | */ |
---|
152 | double |
---|
153 | EllipCyl20(double dp[], double q) |
---|
154 | { |
---|
155 | int i,j; |
---|
156 | double Pi; |
---|
157 | double scale,ra,nu,length,delrho,bkg; //local variables of coefficient wave |
---|
158 | int nordi=76; //order of integration |
---|
159 | int nordj=20; |
---|
160 | double va,vb; //upper and lower integration limits |
---|
161 | double summ,zi,yyy,answer,vell; //running tally of integration |
---|
162 | double summj,vaj,vbj,zij,arg,si; //for the inner integration |
---|
163 | |
---|
164 | Pi = 4.0*atan(1.0); |
---|
165 | va = 0; |
---|
166 | vb = 1; //orintational average, outer integral |
---|
167 | vaj=0; |
---|
168 | vbj=Pi; //endpoints of inner integral |
---|
169 | |
---|
170 | summ = 0.0; //initialize intergral |
---|
171 | |
---|
172 | scale = dp[0]; //make local copies in case memory moves |
---|
173 | ra = dp[1]; |
---|
174 | nu = dp[2]; |
---|
175 | length = dp[3]; |
---|
176 | delrho = dp[4]; |
---|
177 | bkg = dp[5]; |
---|
178 | for(i=0;i<nordi;i++) { |
---|
179 | //setup inner integral over the ellipsoidal cross-section |
---|
180 | summj=0; |
---|
181 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy |
---|
182 | arg = ra*sqrt(1-zi*zi); |
---|
183 | for(j=0;j<nordj;j++) { |
---|
184 | //20 gauss points for the inner integral |
---|
185 | zij = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "y" dummy |
---|
186 | yyy = Gauss20Wt[j] * EllipCylKernel(q,arg,nu,zij); |
---|
187 | summj += yyy; |
---|
188 | } |
---|
189 | //now calculate the value of the inner integral |
---|
190 | answer = (vbj-vaj)/2.0*summj; |
---|
191 | //divide integral by Pi |
---|
192 | answer /=Pi; |
---|
193 | |
---|
194 | //now calculate outer integral |
---|
195 | arg = q*length*zi/2; |
---|
196 | if (arg == 0){ |
---|
197 | si = 1; |
---|
198 | }else{ |
---|
199 | si = sin(arg) * sin(arg) / arg / arg; |
---|
200 | } |
---|
201 | yyy = Gauss76Wt[i] * answer * si; |
---|
202 | summ += yyy; |
---|
203 | } |
---|
204 | |
---|
205 | answer = (vb-va)/2.0*summ; |
---|
206 | // Multiply by contrast^2 |
---|
207 | answer *= delrho*delrho; |
---|
208 | //normalize by cylinder volume |
---|
209 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
210 | vell = Pi*ra*(nu*ra)*length; |
---|
211 | answer *= vell; |
---|
212 | //convert to [cm-1] |
---|
213 | answer *= 1.0e8; |
---|
214 | //Scale |
---|
215 | answer *= scale; |
---|
216 | // add in the background |
---|
217 | answer += bkg; |
---|
218 | |
---|
219 | return answer; |
---|
220 | } |
---|
221 | |
---|
222 | /* TriaxialEllipsoidX : calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x |
---|
223 | |
---|
224 | Uses 76 pt Gaussian quadrature for both integrals |
---|
225 | |
---|
226 | Warning: |
---|
227 | The call to WaveData() below returns a pointer to the middle |
---|
228 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
229 | calculations could cause memory to move, you should copy the coefficient |
---|
230 | values to local variables or an array before such operations. |
---|
231 | */ |
---|
232 | double |
---|
233 | TriaxialEllipsoid(double dp[], double q) |
---|
234 | { |
---|
235 | int i,j; |
---|
236 | double Pi; |
---|
237 | double scale,aa,bb,cc,delrho,bkg; //local variables of coefficient wave |
---|
238 | int nordi=76; //order of integration |
---|
239 | int nordj=76; |
---|
240 | double va,vb; //upper and lower integration limits |
---|
241 | double summ,zi,yyy,answer; //running tally of integration |
---|
242 | double summj,vaj,vbj,zij; //for the inner integration |
---|
243 | |
---|
244 | Pi = 4.0*atan(1.0); |
---|
245 | va = 0; |
---|
246 | vb = 1; //orintational average, outer integral |
---|
247 | vaj = 0; |
---|
248 | vbj = 1; //endpoints of inner integral |
---|
249 | |
---|
250 | summ = 0.0; //initialize intergral |
---|
251 | |
---|
252 | scale = dp[0]; //make local copies in case memory moves |
---|
253 | aa = dp[1]; |
---|
254 | bb = dp[2]; |
---|
255 | cc = dp[3]; |
---|
256 | delrho = dp[4]; |
---|
257 | bkg = dp[5]; |
---|
258 | for(i=0;i<nordi;i++) { |
---|
259 | //setup inner integral over the ellipsoidal cross-section |
---|
260 | summj=0; |
---|
261 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the "x" dummy |
---|
262 | for(j=0;j<nordj;j++) { |
---|
263 | //20 gauss points for the inner integral |
---|
264 | zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "y" dummy |
---|
265 | yyy = Gauss76Wt[j] * TriaxialKernel(q,aa,bb,cc,zi,zij); |
---|
266 | summj += yyy; |
---|
267 | } |
---|
268 | //now calculate the value of the inner integral |
---|
269 | answer = (vbj-vaj)/2.0*summj; |
---|
270 | |
---|
271 | //now calculate outer integral |
---|
272 | yyy = Gauss76Wt[i] * answer; |
---|
273 | summ += yyy; |
---|
274 | } //final scaling is done at the end of the function, after the NT_FP64 case |
---|
275 | |
---|
276 | answer = (vb-va)/2.0*summ; |
---|
277 | // Multiply by contrast^2 |
---|
278 | answer *= delrho*delrho; |
---|
279 | //normalize by ellipsoid volume |
---|
280 | answer *= 4*Pi/3*aa*bb*cc; |
---|
281 | //convert to [cm-1] |
---|
282 | answer *= 1.0e8; |
---|
283 | //Scale |
---|
284 | answer *= scale; |
---|
285 | // add in the background |
---|
286 | answer += bkg; |
---|
287 | |
---|
288 | return answer; |
---|
289 | } |
---|
290 | |
---|
291 | /* ParallelepipedX : calculates the form factor of a Parallelepiped (a rectangular solid) |
---|
292 | at the given x-value p->x |
---|
293 | |
---|
294 | Uses 76 pt Gaussian quadrature for both integrals |
---|
295 | |
---|
296 | Warning: |
---|
297 | The call to WaveData() below returns a pointer to the middle |
---|
298 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
299 | calculations could cause memory to move, you should copy the coefficient |
---|
300 | values to local variables or an array before such operations. |
---|
301 | */ |
---|
302 | double |
---|
303 | Parallelepiped(double dp[], double q) |
---|
304 | { |
---|
305 | int i,j; |
---|
306 | double scale,aa,bb,cc,delrho,bkg; //local variables of coefficient wave |
---|
307 | int nordi=76; //order of integration |
---|
308 | int nordj=76; |
---|
309 | double va,vb; //upper and lower integration limits |
---|
310 | double summ,yyy,answer; //running tally of integration |
---|
311 | double summj,vaj,vbj; //for the inner integration |
---|
312 | double mu,mudum,arg,sigma,uu,vol; |
---|
313 | |
---|
314 | |
---|
315 | // Pi = 4.0*atan(1.0); |
---|
316 | va = 0; |
---|
317 | vb = 1; //orintational average, outer integral |
---|
318 | vaj = 0; |
---|
319 | vbj = 1; //endpoints of inner integral |
---|
320 | |
---|
321 | summ = 0.0; //initialize intergral |
---|
322 | |
---|
323 | scale = dp[0]; //make local copies in case memory moves |
---|
324 | aa = dp[1]; |
---|
325 | bb = dp[2]; |
---|
326 | cc = dp[3]; |
---|
327 | delrho = dp[4]; |
---|
328 | bkg = dp[5]; |
---|
329 | |
---|
330 | mu = q*bb; |
---|
331 | vol = aa*bb*cc; |
---|
332 | // normalize all WRT bb |
---|
333 | aa = aa/bb; |
---|
334 | cc = cc/bb; |
---|
335 | |
---|
336 | for(i=0;i<nordi;i++) { |
---|
337 | //setup inner integral over the ellipsoidal cross-section |
---|
338 | summj=0; |
---|
339 | sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; //the outer dummy |
---|
340 | |
---|
341 | for(j=0;j<nordj;j++) { |
---|
342 | //76 gauss points for the inner integral |
---|
343 | uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the inner dummy |
---|
344 | mudum = mu*sqrt(1-sigma*sigma); |
---|
345 | yyy = Gauss76Wt[j] * PPKernel(aa,mudum,uu); |
---|
346 | summj += yyy; |
---|
347 | } |
---|
348 | //now calculate the value of the inner integral |
---|
349 | answer = (vbj-vaj)/2.0*summj; |
---|
350 | |
---|
351 | arg = mu*cc*sigma/2; |
---|
352 | if ( arg == 0 ) { |
---|
353 | answer *= 1; |
---|
354 | } else { |
---|
355 | answer *= sin(arg)*sin(arg)/arg/arg; |
---|
356 | } |
---|
357 | |
---|
358 | //now sum up the outer integral |
---|
359 | yyy = Gauss76Wt[i] * answer; |
---|
360 | summ += yyy; |
---|
361 | } //final scaling is done at the end of the function, after the NT_FP64 case |
---|
362 | |
---|
363 | answer = (vb-va)/2.0*summ; |
---|
364 | // Multiply by contrast^2 |
---|
365 | answer *= delrho*delrho; |
---|
366 | //normalize by volume |
---|
367 | answer *= vol; |
---|
368 | //convert to [cm-1] |
---|
369 | answer *= 1.0e8; |
---|
370 | //Scale |
---|
371 | answer *= scale; |
---|
372 | // add in the background |
---|
373 | answer += bkg; |
---|
374 | |
---|
375 | return answer; |
---|
376 | } |
---|
377 | |
---|
378 | /* HollowCylinderX : calculates the form factor of a Hollow Cylinder |
---|
379 | at the given x-value p->x |
---|
380 | |
---|
381 | Uses 76 pt Gaussian quadrature for the single integral |
---|
382 | |
---|
383 | Warning: |
---|
384 | The call to WaveData() below returns a pointer to the middle |
---|
385 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
386 | calculations could cause memory to move, you should copy the coefficient |
---|
387 | values to local variables or an array before such operations. |
---|
388 | */ |
---|
389 | double |
---|
390 | HollowCylinder(double dp[], double q) |
---|
391 | { |
---|
392 | int i; |
---|
393 | double scale,rcore,rshell,length,delrho,bkg; //local variables of coefficient wave |
---|
394 | int nord=76; //order of integration |
---|
395 | double va,vb,zi; //upper and lower integration limits |
---|
396 | double summ,answer,pi; //running tally of integration |
---|
397 | |
---|
398 | pi = 4.0*atan(1.0); |
---|
399 | va = 0; |
---|
400 | vb = 1; //limits of numerical integral |
---|
401 | |
---|
402 | summ = 0.0; //initialize intergral |
---|
403 | |
---|
404 | scale = dp[0]; //make local copies in case memory moves |
---|
405 | rcore = dp[1]; |
---|
406 | rshell = dp[2]; |
---|
407 | length = dp[3]; |
---|
408 | delrho = dp[4]; |
---|
409 | bkg = dp[5]; |
---|
410 | |
---|
411 | for(i=0;i<nord;i++) { |
---|
412 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; |
---|
413 | summ += Gauss76Wt[i] * HolCylKernel(q, rcore, rshell, length, zi); |
---|
414 | } |
---|
415 | |
---|
416 | answer = (vb-va)/2.0*summ; |
---|
417 | // Multiply by contrast^2 |
---|
418 | answer *= delrho*delrho; |
---|
419 | //normalize by volume |
---|
420 | answer *= pi*(rshell*rshell-rcore*rcore)*length; |
---|
421 | //convert to [cm-1] |
---|
422 | answer *= 1.0e8; |
---|
423 | //Scale |
---|
424 | answer *= scale; |
---|
425 | // add in the background |
---|
426 | answer += bkg; |
---|
427 | |
---|
428 | return answer; |
---|
429 | } |
---|
430 | |
---|
431 | /* EllipsoidFormX : calculates the form factor of an ellipsoid of revolution with semiaxes a:a:nua |
---|
432 | at the given x-value p->x |
---|
433 | |
---|
434 | Uses 76 pt Gaussian quadrature for the single integral |
---|
435 | |
---|
436 | Warning: |
---|
437 | The call to WaveData() below returns a pointer to the middle |
---|
438 | of an unlocked Macintosh handle. In the unlikely event that your |
---|
439 | calculations could cause memory to move, you should copy the coefficient |
---|
440 | values to local variables or an array before such operations. |
---|
441 | */ |
---|
442 | double |
---|
443 | EllipsoidForm(double dp[], double q) |
---|
444 | { |
---|
445 | int i; |
---|
446 | double scale,a,nua,delrho,bkg; //local variables of coefficient wave |
---|
447 | int nord=76; //order of integration |
---|
448 | double va,vb,zi; //upper and lower integration limits |
---|
449 | double summ,answer,pi; //running tally of integration |
---|
450 | |
---|
451 | pi = 4.0*atan(1.0); |
---|
452 | va = 0; |
---|
453 | vb = 1; //limits of numerical integral |
---|
454 | |
---|
455 | summ = 0.0; //initialize intergral |
---|
456 | |
---|
457 | scale = dp[0]; //make local copies in case memory moves |
---|
458 | nua = dp[1]; |
---|
459 | a = dp[2]; |
---|
460 | delrho = dp[3]; |
---|
461 | bkg = dp[4]; |
---|
462 | |
---|
463 | for(i=0;i<nord;i++) { |
---|
464 | zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; |
---|
465 | summ += Gauss76Wt[i] * EllipsoidKernel(q, a, nua, zi); |
---|
466 | } |
---|
467 | |
---|
468 | answer = (vb-va)/2.0*summ; |
---|
469 | // Multiply by contrast^2 |
---|
470 | answer *= delrho*delrho; |
---|
471 | //normalize by volume |
---|
472 | answer *= 4*pi/3*a*a*nua; |
---|
473 | //convert to [cm-1] |
---|
474 | answer *= 1.0e8; |
---|
475 | //Scale |
---|
476 | answer *= scale; |
---|
477 | // add in the background |
---|
478 | answer += bkg; |
---|
479 | |
---|
480 | return answer; |
---|
481 | } |
---|
482 | |
---|
483 | |
---|
484 | /* Cyl_PolyRadiusX : calculates the form factor of a cylinder at the given x-value p->x |
---|
485 | the cylinder has a polydisperse cross section |
---|
486 | |
---|
487 | */ |
---|
488 | double |
---|
489 | Cyl_PolyRadius(double dp[], double q) |
---|
490 | { |
---|
491 | int i; |
---|
492 | double scale,radius,length,pd,delrho,bkg; //local variables of coefficient wave |
---|
493 | int nord=20; //order of integration |
---|
494 | double uplim,lolim; //upper and lower integration limits |
---|
495 | double summ,zi,yyy,answer,Vpoly; //running tally of integration |
---|
496 | double range,zz,Pi; |
---|
497 | |
---|
498 | Pi = 4.0*atan(1.0); |
---|
499 | range = 3.4; |
---|
500 | |
---|
501 | summ = 0.0; //initialize intergral |
---|
502 | |
---|
503 | scale = dp[0]; //make local copies in case memory moves |
---|
504 | radius = dp[1]; |
---|
505 | length = dp[2]; |
---|
506 | pd = dp[3]; |
---|
507 | delrho = dp[4]; |
---|
508 | bkg = dp[5]; |
---|
509 | |
---|
510 | zz = (1.0/pd)*(1.0/pd) - 1.0; |
---|
511 | |
---|
512 | lolim = radius*(1.0-range*pd); //set the upper/lower limits to cover the distribution |
---|
513 | if(lolim<0) { |
---|
514 | lolim = 0; |
---|
515 | } |
---|
516 | if(pd>0.3) { |
---|
517 | range = 3.4 + (pd-0.3)*18.0; |
---|
518 | } |
---|
519 | uplim = radius*(1.0+range*pd); |
---|
520 | |
---|
521 | for(i=0;i<nord;i++) { |
---|
522 | zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
523 | yyy = Gauss20Wt[i] * Cyl_PolyRadKernel(q, radius, length, zz, delrho, zi); |
---|
524 | summ += yyy; |
---|
525 | } |
---|
526 | |
---|
527 | answer = (uplim-lolim)/2.0*summ; |
---|
528 | //normalize by average cylinder volume |
---|
529 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
530 | Vpoly=Pi*radius*radius*length*(zz+2.0)/(zz+1.0); |
---|
531 | answer /= Vpoly; |
---|
532 | //convert to [cm-1] |
---|
533 | answer *= 1.0e8; |
---|
534 | //Scale |
---|
535 | answer *= scale; |
---|
536 | // add in the background |
---|
537 | answer += bkg; |
---|
538 | |
---|
539 | return answer; |
---|
540 | } |
---|
541 | |
---|
542 | /* Cyl_PolyLengthX : calculates the form factor of a cylinder at the given x-value p->x |
---|
543 | the cylinder has a polydisperse Length |
---|
544 | |
---|
545 | */ |
---|
546 | double |
---|
547 | Cyl_PolyLength(double dp[], double q) |
---|
548 | { |
---|
549 | int i; |
---|
550 | double scale,radius,length,pd,delrho,bkg; //local variables of coefficient wave |
---|
551 | int nord=20; //order of integration |
---|
552 | double uplim,lolim; //upper and lower integration limits |
---|
553 | double summ,zi,yyy,answer,Vpoly; //running tally of integration |
---|
554 | double range,zz,Pi; |
---|
555 | |
---|
556 | |
---|
557 | Pi = 4.0*atan(1.0); |
---|
558 | range = 3.4; |
---|
559 | |
---|
560 | summ = 0.0; //initialize intergral |
---|
561 | |
---|
562 | scale = dp[0]; //make local copies in case memory moves |
---|
563 | radius = dp[1]; |
---|
564 | length = dp[2]; |
---|
565 | pd = dp[3]; |
---|
566 | delrho = dp[4]; |
---|
567 | bkg = dp[5]; |
---|
568 | |
---|
569 | zz = (1.0/pd)*(1.0/pd) - 1.0; |
---|
570 | |
---|
571 | lolim = length*(1.0-range*pd); //set the upper/lower limits to cover the distribution |
---|
572 | if(lolim<0) { |
---|
573 | lolim = 0; |
---|
574 | } |
---|
575 | if(pd>0.3) { |
---|
576 | range = 3.4 + (pd-0.3)*18.0; |
---|
577 | } |
---|
578 | uplim = length*(1.0+range*pd); |
---|
579 | |
---|
580 | for(i=0;i<nord;i++) { |
---|
581 | zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
582 | yyy = Gauss20Wt[i] * Cyl_PolyLenKernel(q, radius, length, zz, delrho, zi); |
---|
583 | summ += yyy; |
---|
584 | } |
---|
585 | |
---|
586 | answer = (uplim-lolim)/2.0*summ; |
---|
587 | //normalize by average cylinder volume (first moment) |
---|
588 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
589 | Vpoly=Pi*radius*radius*length; |
---|
590 | answer /= Vpoly; |
---|
591 | //convert to [cm-1] |
---|
592 | answer *= 1.0e8; |
---|
593 | //Scale |
---|
594 | answer *= scale; |
---|
595 | // add in the background |
---|
596 | answer += bkg; |
---|
597 | |
---|
598 | return answer; |
---|
599 | } |
---|
600 | |
---|
601 | /* CoreShellCylinderX : calculates the form factor of a cylinder at the given x-value p->x |
---|
602 | the cylinder has a core-shell structure |
---|
603 | |
---|
604 | */ |
---|
605 | double |
---|
606 | CoreShellCylinder(double dp[], double q) |
---|
607 | { |
---|
608 | int i; |
---|
609 | double scale,rcore,length,bkg; //local variables of coefficient wave |
---|
610 | double thick,rhoc,rhos,rhosolv; |
---|
611 | int nord=76; //order of integration |
---|
612 | double uplim,lolim,halfheight; //upper and lower integration limits |
---|
613 | double summ,zi,yyy,answer,Vcyl; //running tally of integration |
---|
614 | double Pi; |
---|
615 | |
---|
616 | Pi = 4.0*atan(1.0); |
---|
617 | |
---|
618 | lolim = 0.0; |
---|
619 | uplim = Pi/2.0; |
---|
620 | |
---|
621 | summ = 0.0; //initialize intergral |
---|
622 | |
---|
623 | scale = dp[0]; //make local copies in case memory moves |
---|
624 | rcore = dp[1]; |
---|
625 | thick = dp[2]; |
---|
626 | length = dp[3]; |
---|
627 | rhoc = dp[4]; |
---|
628 | rhos = dp[5]; |
---|
629 | rhosolv = dp[6]; |
---|
630 | bkg = dp[7]; |
---|
631 | |
---|
632 | halfheight = length/2.0; |
---|
633 | |
---|
634 | for(i=0;i<nord;i++) { |
---|
635 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
636 | yyy = Gauss76Wt[i] * CoreShellCylKernel(q, rcore, thick, rhoc,rhos,rhosolv, halfheight, zi); |
---|
637 | summ += yyy; |
---|
638 | } |
---|
639 | |
---|
640 | answer = (uplim-lolim)/2.0*summ; |
---|
641 | // length is the total core length |
---|
642 | Vcyl=Pi*(rcore+thick)*(rcore+thick)*(length+2.0*thick); |
---|
643 | answer /= Vcyl; |
---|
644 | //convert to [cm-1] |
---|
645 | answer *= 1.0e8; |
---|
646 | //Scale |
---|
647 | answer *= scale; |
---|
648 | // add in the background |
---|
649 | answer += bkg; |
---|
650 | |
---|
651 | return answer; |
---|
652 | } |
---|
653 | |
---|
654 | |
---|
655 | /* PolyCoShCylinderX : calculates the form factor of a core-shell cylinder at the given x-value p->x |
---|
656 | the cylinder has a polydisperse CORE radius |
---|
657 | |
---|
658 | */ |
---|
659 | double |
---|
660 | PolyCoShCylinder(double dp[], double q) |
---|
661 | { |
---|
662 | int i; |
---|
663 | double scale,radius,length,sigma,bkg; //local variables of coefficient wave |
---|
664 | double rad,radthick,facthick,rhoc,rhos,rhosolv; |
---|
665 | int nord=20; //order of integration |
---|
666 | double uplim,lolim; //upper and lower integration limits |
---|
667 | double summ,yyy,answer,Vpoly; //running tally of integration |
---|
668 | double Pi,AR,Rsqrsumm,Rsqryyy,Rsqr; |
---|
669 | |
---|
670 | Pi = 4.0*atan(1.0); |
---|
671 | |
---|
672 | summ = 0.0; //initialize intergral |
---|
673 | Rsqrsumm = 0.0; |
---|
674 | |
---|
675 | scale = dp[0]; |
---|
676 | radius = dp[1]; |
---|
677 | sigma = dp[2]; //sigma is the standard mean deviation |
---|
678 | length = dp[3]; |
---|
679 | radthick = dp[4]; |
---|
680 | facthick= dp[5]; |
---|
681 | rhoc = dp[6]; |
---|
682 | rhos = dp[7]; |
---|
683 | rhosolv = dp[8]; |
---|
684 | bkg = dp[9]; |
---|
685 | |
---|
686 | lolim = exp(log(radius)-(4.*sigma)); |
---|
687 | if (lolim<0) { |
---|
688 | lolim=0; //to avoid numerical error when va<0 (-ve r value) |
---|
689 | } |
---|
690 | uplim = exp(log(radius)+(4.*sigma)); |
---|
691 | |
---|
692 | for(i=0;i<nord;i++) { |
---|
693 | rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
694 | AR=(1.0/(rad*sigma*sqrt(2.0*Pi)))*exp(-(0.5*((log(radius/rad))/sigma)*((log(radius/rad))/sigma))); |
---|
695 | yyy = AR* Gauss20Wt[i] * CSCylIntegration(q,rad,radthick,facthick,rhoc,rhos,rhosolv,length); |
---|
696 | Rsqryyy= Gauss20Wt[i] * AR * (rad+radthick)*(rad+radthick); //SRK normalize to total dimensions |
---|
697 | summ += yyy; |
---|
698 | Rsqrsumm += Rsqryyy; |
---|
699 | } |
---|
700 | |
---|
701 | answer = (uplim-lolim)/2.0*summ; |
---|
702 | Rsqr = (uplim-lolim)/2.0*Rsqrsumm; |
---|
703 | //normalize by average cylinder volume |
---|
704 | Vpoly = Pi*Rsqr*(length+2*facthick); |
---|
705 | answer /= Vpoly; |
---|
706 | //convert to [cm-1] |
---|
707 | answer *= 1.0e8; |
---|
708 | //Scale |
---|
709 | answer *= scale; |
---|
710 | // add in the background |
---|
711 | answer += bkg; |
---|
712 | |
---|
713 | return answer; |
---|
714 | } |
---|
715 | |
---|
716 | /* OblateFormX : calculates the form factor of a core-shell Oblate ellipsoid at the given x-value p->x |
---|
717 | the ellipsoid has a core-shell structure |
---|
718 | |
---|
719 | */ |
---|
720 | double |
---|
721 | OblateForm(double dp[], double q) |
---|
722 | { |
---|
723 | int i; |
---|
724 | double scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg; |
---|
725 | int nord=76; //order of integration |
---|
726 | double uplim,lolim; //upper and lower integration limits |
---|
727 | double summ,zi,yyy,answer,oblatevol; //running tally of integration |
---|
728 | double Pi; |
---|
729 | |
---|
730 | Pi = 4.0*atan(1.0); |
---|
731 | |
---|
732 | lolim = 0.0; |
---|
733 | uplim = 1.0; |
---|
734 | |
---|
735 | summ = 0.0; //initialize intergral |
---|
736 | |
---|
737 | |
---|
738 | scale = dp[0]; //make local copies in case memory moves |
---|
739 | crmaj = dp[1]; |
---|
740 | crmin = dp[2]; |
---|
741 | trmaj = dp[3]; |
---|
742 | trmin = dp[4]; |
---|
743 | delpc = dp[5]; |
---|
744 | delps = dp[6]; |
---|
745 | bkg = dp[7]; |
---|
746 | |
---|
747 | for(i=0;i<nord;i++) { |
---|
748 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
749 | yyy = Gauss76Wt[i] * gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,q); |
---|
750 | summ += yyy; |
---|
751 | } |
---|
752 | |
---|
753 | answer = (uplim-lolim)/2.0*summ; |
---|
754 | // normalize by particle volume |
---|
755 | oblatevol = 4*Pi/3*trmaj*trmaj*trmin; |
---|
756 | answer /= oblatevol; |
---|
757 | |
---|
758 | //convert to [cm-1] |
---|
759 | answer *= 1.0e8; |
---|
760 | //Scale |
---|
761 | answer *= scale; |
---|
762 | // add in the background |
---|
763 | answer += bkg; |
---|
764 | |
---|
765 | return answer; |
---|
766 | } |
---|
767 | |
---|
768 | /* ProlateFormX : calculates the form factor of a core-shell Prolate ellipsoid at the given x-value p->x |
---|
769 | the ellipsoid has a core-shell structure |
---|
770 | |
---|
771 | */ |
---|
772 | double |
---|
773 | ProlateForm(double dp[], double q) |
---|
774 | { |
---|
775 | int i; |
---|
776 | double scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg; |
---|
777 | int nord=76; //order of integration |
---|
778 | double uplim,lolim; //upper and lower integration limits |
---|
779 | double summ,zi,yyy,answer,prolatevol; //running tally of integration |
---|
780 | double Pi; |
---|
781 | |
---|
782 | Pi = 4.0*atan(1.0); |
---|
783 | |
---|
784 | lolim = 0.0; |
---|
785 | uplim = 1.0; |
---|
786 | |
---|
787 | summ = 0.0; //initialize intergral |
---|
788 | |
---|
789 | scale = dp[0]; //make local copies in case memory moves |
---|
790 | crmaj = dp[1]; |
---|
791 | crmin = dp[2]; |
---|
792 | trmaj = dp[3]; |
---|
793 | trmin = dp[4]; |
---|
794 | delpc = dp[5]; |
---|
795 | delps = dp[6]; |
---|
796 | bkg = dp[7]; |
---|
797 | |
---|
798 | for(i=0;i<nord;i++) { |
---|
799 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
800 | yyy = Gauss76Wt[i] * gfn2(zi,crmaj,crmin,trmaj,trmin,delpc,delps,q); |
---|
801 | summ += yyy; |
---|
802 | } |
---|
803 | |
---|
804 | answer = (uplim-lolim)/2.0*summ; |
---|
805 | // normalize by particle volume |
---|
806 | prolatevol = 4*Pi/3*trmaj*trmin*trmin; |
---|
807 | answer /= prolatevol; |
---|
808 | |
---|
809 | //convert to [cm-1] |
---|
810 | answer *= 1.0e8; |
---|
811 | //Scale |
---|
812 | answer *= scale; |
---|
813 | // add in the background |
---|
814 | answer += bkg; |
---|
815 | |
---|
816 | return answer; |
---|
817 | } |
---|
818 | |
---|
819 | |
---|
820 | /* StackedDiscsX : calculates the form factor of a stacked "tactoid" of core shell disks |
---|
821 | like clay platelets that are not exfoliated |
---|
822 | |
---|
823 | */ |
---|
824 | double |
---|
825 | StackedDiscs(double dp[], double q) |
---|
826 | { |
---|
827 | int i; |
---|
828 | double scale,length,bkg,rcore,thick,rhoc,rhol,rhosolv,N,gsd; //local variables of coefficient wave |
---|
829 | double va,vb,vcyl,summ,yyy,zi,halfheight,d,answer; |
---|
830 | int nord=76; //order of integration |
---|
831 | double Pi; |
---|
832 | |
---|
833 | |
---|
834 | Pi = 4.0*atan(1.0); |
---|
835 | |
---|
836 | va = 0.0; |
---|
837 | vb = Pi/2.0; |
---|
838 | |
---|
839 | summ = 0.0; //initialize intergral |
---|
840 | |
---|
841 | scale = dp[0]; |
---|
842 | rcore = dp[1]; |
---|
843 | length = dp[2]; |
---|
844 | thick = dp[3]; |
---|
845 | rhoc = dp[4]; |
---|
846 | rhol = dp[5]; |
---|
847 | rhosolv = dp[6]; |
---|
848 | N = dp[7]; |
---|
849 | gsd = dp[8]; |
---|
850 | bkg = dp[9]; |
---|
851 | |
---|
852 | d=2.0*thick+length; |
---|
853 | halfheight = length/2.0; |
---|
854 | |
---|
855 | for(i=0;i<nord;i++) { |
---|
856 | zi = ( Gauss76Z[i]*(vb-va) + vb + va )/2.0; |
---|
857 | yyy = Gauss76Wt[i] * Stackdisc_kern(q, rcore, rhoc,rhol,rhosolv, halfheight,thick,zi,gsd,d,N); |
---|
858 | summ += yyy; |
---|
859 | } |
---|
860 | |
---|
861 | answer = (vb-va)/2.0*summ; |
---|
862 | // length is the total core length |
---|
863 | vcyl=Pi*rcore*rcore*(2.0*thick+length)*N; |
---|
864 | answer /= vcyl; |
---|
865 | //Convert to [cm-1] |
---|
866 | answer *= 1.0e8; |
---|
867 | //Scale |
---|
868 | answer *= scale; |
---|
869 | // add in the background |
---|
870 | answer += bkg; |
---|
871 | |
---|
872 | return answer; |
---|
873 | } |
---|
874 | |
---|
875 | |
---|
876 | /* LamellarFFX : calculates the form factor of a lamellar structure - no S(q) effects included |
---|
877 | |
---|
878 | */ |
---|
879 | double |
---|
880 | LamellarFF(double dp[], double q) |
---|
881 | { |
---|
882 | double scale,del,sig,contr,bkg; //local variables of coefficient wave |
---|
883 | double inten, qval,Pq; |
---|
884 | double Pi; |
---|
885 | |
---|
886 | |
---|
887 | Pi = 4.0*atan(1.0); |
---|
888 | scale = dp[0]; |
---|
889 | del = dp[1]; |
---|
890 | sig = dp[2]*del; |
---|
891 | contr = dp[3]; |
---|
892 | bkg = dp[4]; |
---|
893 | qval = q; |
---|
894 | |
---|
895 | Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig)); |
---|
896 | |
---|
897 | inten = 2.0*Pi*scale*Pq/(qval*qval); //this is now dimensionless... |
---|
898 | |
---|
899 | inten /= del; //normalize by the thickness (in A) |
---|
900 | |
---|
901 | inten *= 1.0e8; // 1/A to 1/cm |
---|
902 | |
---|
903 | return(inten+bkg); |
---|
904 | } |
---|
905 | /* LamellarFFX : calculates the form factor of a lamellar structure - no S(q) effects included |
---|
906 | -NO polydispersion included |
---|
907 | */ |
---|
908 | double lamellar_kernel(double dp[], double q){ |
---|
909 | double scale,del,sld_bi,sld_sol,contr,bkg; //local variables of coefficient wave |
---|
910 | double inten, qval,Pq; |
---|
911 | double Pi; |
---|
912 | |
---|
913 | |
---|
914 | Pi = 4.0*atan(1.0); |
---|
915 | scale = dp[0]; |
---|
916 | del = dp[1]; |
---|
917 | sld_bi = dp[2]; |
---|
918 | sld_sol = dp[3]; |
---|
919 | bkg = dp[4]; |
---|
920 | qval = q; |
---|
921 | contr = sld_bi -sld_sol; |
---|
922 | |
---|
923 | Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)); |
---|
924 | |
---|
925 | inten = 2.0*Pi*scale*Pq/(qval*qval); //this is now dimensionless... |
---|
926 | |
---|
927 | inten /= del; //normalize by the thickness (in A) |
---|
928 | |
---|
929 | inten *= 1.0e8; // 1/A to 1/cm |
---|
930 | |
---|
931 | return(inten+bkg); |
---|
932 | } |
---|
933 | /* LamellarPSX : calculates the form factor of a lamellar structure - with S(q) effects included |
---|
934 | ------- |
---|
935 | ------- resolution effects ARE included, but only a CONSTANT default value, not the real q-dependent resolution!! |
---|
936 | |
---|
937 | */ |
---|
938 | double |
---|
939 | LamellarPS(double dp[], double q) |
---|
940 | { |
---|
941 | double scale,dd,del,sig,contr,NN,Cp,bkg; //local variables of coefficient wave |
---|
942 | double inten, qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ; |
---|
943 | double Pi,Euler,dQDefault,fii; |
---|
944 | int ii,NNint; |
---|
945 | |
---|
946 | Euler = 0.5772156649; // Euler's constant |
---|
947 | dQDefault = 0.0;//0.0025; //[=] 1/A, q-resolution, default value |
---|
948 | dQ = dQDefault; |
---|
949 | |
---|
950 | Pi = 4.0*atan(1.0); |
---|
951 | qval = q; |
---|
952 | |
---|
953 | scale = dp[0]; |
---|
954 | dd = dp[1]; |
---|
955 | del = dp[2]; |
---|
956 | sig = dp[3]*del; |
---|
957 | contr = dp[4]; |
---|
958 | NN = trunc(dp[5]); //be sure that NN is an integer |
---|
959 | Cp = dp[6]; |
---|
960 | bkg = dp[7]; |
---|
961 | |
---|
962 | Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig)); |
---|
963 | |
---|
964 | NNint = (int)NN; //cast to an integer for the loop |
---|
965 | ii=0; |
---|
966 | Sq = 0.0; |
---|
967 | for(ii=1;ii<(NNint-1);ii+=1) { |
---|
968 | |
---|
969 | fii = (double)ii; //do I really need to do this? |
---|
970 | |
---|
971 | temp = 0.0; |
---|
972 | alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler); |
---|
973 | t1 = 2.0*dQ*dQ*dd*dd*alpha; |
---|
974 | t2 = 2.0*qval*qval*dd*dd*alpha; |
---|
975 | t3 = dQ*dQ*dd*dd*ii*ii; |
---|
976 | |
---|
977 | temp = 1.0-ii/NN; |
---|
978 | temp *= cos(dd*qval*ii/(1.0+t1)); |
---|
979 | temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); |
---|
980 | temp /= sqrt(1.0+t1); |
---|
981 | |
---|
982 | Sq += temp; |
---|
983 | } |
---|
984 | |
---|
985 | Sq *= 2.0; |
---|
986 | Sq += 1.0; |
---|
987 | |
---|
988 | inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval); |
---|
989 | |
---|
990 | inten *= 1.0e8; // 1/A to 1/cm |
---|
991 | |
---|
992 | return(inten+bkg); |
---|
993 | } |
---|
994 | |
---|
995 | |
---|
996 | /* LamellarPS_HGX : calculates the form factor of a lamellar structure - with S(q) effects included |
---|
997 | ------- |
---|
998 | ------- resolution effects ARE included, but only a CONSTANT default value, not the real q-dependent resolution!! |
---|
999 | |
---|
1000 | */ |
---|
1001 | double |
---|
1002 | LamellarPS_HG(double dp[], double q) |
---|
1003 | { |
---|
1004 | double scale,dd,delT,delH,SLD_T,SLD_H,SLD_S,NN,Cp,bkg; //local variables of coefficient wave |
---|
1005 | double inten,qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ,drh,drt; |
---|
1006 | double Pi,Euler,dQDefault,fii; |
---|
1007 | int ii,NNint; |
---|
1008 | |
---|
1009 | |
---|
1010 | Euler = 0.5772156649; // Euler's constant |
---|
1011 | dQDefault = 0; //0.0025; //[=] 1/A, q-resolution, default value |
---|
1012 | dQ = dQDefault; |
---|
1013 | |
---|
1014 | Pi = 4.0*atan(1.0); |
---|
1015 | qval= q; |
---|
1016 | |
---|
1017 | scale = dp[0]; |
---|
1018 | dd = dp[1]; |
---|
1019 | delT = dp[2]; |
---|
1020 | delH = dp[3]; |
---|
1021 | SLD_T = dp[4]; |
---|
1022 | SLD_H = dp[5]; |
---|
1023 | SLD_S = dp[6]; |
---|
1024 | NN = trunc(dp[7]); //be sure that NN is an integer |
---|
1025 | Cp = dp[8]; |
---|
1026 | bkg = dp[9]; |
---|
1027 | |
---|
1028 | |
---|
1029 | drh = SLD_H - SLD_S; |
---|
1030 | drt = SLD_T - SLD_S; //correction 13FEB06 by L.Porcar |
---|
1031 | |
---|
1032 | Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT); |
---|
1033 | Pq *= Pq; |
---|
1034 | Pq *= 4.0/(qval*qval); |
---|
1035 | |
---|
1036 | NNint = (int)NN; //cast to an integer for the loop |
---|
1037 | ii=0; |
---|
1038 | Sq = 0.0; |
---|
1039 | for(ii=1;ii<(NNint-1);ii+=1) { |
---|
1040 | |
---|
1041 | fii = (double)ii; //do I really need to do this? |
---|
1042 | |
---|
1043 | temp = 0.0; |
---|
1044 | alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler); |
---|
1045 | t1 = 2.0*dQ*dQ*dd*dd*alpha; |
---|
1046 | t2 = 2.0*qval*qval*dd*dd*alpha; |
---|
1047 | t3 = dQ*dQ*dd*dd*ii*ii; |
---|
1048 | |
---|
1049 | temp = 1.0-ii/NN; |
---|
1050 | temp *= cos(dd*qval*ii/(1.0+t1)); |
---|
1051 | temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); |
---|
1052 | temp /= sqrt(1.0+t1); |
---|
1053 | |
---|
1054 | Sq += temp; |
---|
1055 | } |
---|
1056 | |
---|
1057 | Sq *= 2.0; |
---|
1058 | Sq += 1.0; |
---|
1059 | |
---|
1060 | inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval); |
---|
1061 | |
---|
1062 | inten *= 1.0e8; // 1/A to 1/cm |
---|
1063 | |
---|
1064 | return(inten+bkg); |
---|
1065 | } |
---|
1066 | |
---|
1067 | /* LamellarFF_HGX : calculates the form factor of a lamellar structure - no S(q) effects included |
---|
1068 | but extra SLD for head groups is included |
---|
1069 | |
---|
1070 | */ |
---|
1071 | double |
---|
1072 | LamellarFF_HG(double dp[], double q) |
---|
1073 | { |
---|
1074 | double scale,delT,delH,slds,sldh,sldt,bkg; //local variables of coefficient wave |
---|
1075 | double inten, qval,Pq,drh,drt; |
---|
1076 | double Pi; |
---|
1077 | |
---|
1078 | |
---|
1079 | Pi = 4.0*atan(1.0); |
---|
1080 | qval= q; |
---|
1081 | scale = dp[0]; |
---|
1082 | delT = dp[1]; |
---|
1083 | delH = dp[2]; |
---|
1084 | sldt = dp[3]; |
---|
1085 | sldh = dp[4]; |
---|
1086 | slds = dp[5]; |
---|
1087 | bkg = dp[6]; |
---|
1088 | |
---|
1089 | |
---|
1090 | drh = sldh - slds; |
---|
1091 | drt = sldt - slds; //correction 13FEB06 by L.Porcar |
---|
1092 | |
---|
1093 | Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT); |
---|
1094 | Pq *= Pq; |
---|
1095 | Pq *= 4.0/(qval*qval); |
---|
1096 | |
---|
1097 | inten = 2.0*Pi*scale*Pq/(qval*qval); //dimensionless... |
---|
1098 | |
---|
1099 | inten /= 2.0*(delT+delH); //normalize by the bilayer thickness |
---|
1100 | |
---|
1101 | inten *= 1.0e8; // 1/A to 1/cm |
---|
1102 | |
---|
1103 | return(inten+bkg); |
---|
1104 | } |
---|
1105 | |
---|
1106 | /* FlexExclVolCylX : calculates the form factor of a flexible cylinder with a circular cross section |
---|
1107 | -- incorporates Wei-Ren Chen's fixes - 2006 |
---|
1108 | |
---|
1109 | */ |
---|
1110 | double |
---|
1111 | FlexExclVolCyl(double dp[], double q) |
---|
1112 | { |
---|
1113 | double scale,L,B,bkg,rad,qr,cont; |
---|
1114 | double Pi,flex,crossSect,answer; |
---|
1115 | |
---|
1116 | |
---|
1117 | Pi = 4.0*atan(1.0); |
---|
1118 | |
---|
1119 | scale = dp[0]; //make local copies in case memory moves |
---|
1120 | L = dp[1]; |
---|
1121 | B = dp[2]; |
---|
1122 | rad = dp[3]; |
---|
1123 | cont = dp[4]; |
---|
1124 | bkg = dp[5]; |
---|
1125 | |
---|
1126 | |
---|
1127 | qr = q*rad; |
---|
1128 | |
---|
1129 | flex = Sk_WR(q,L,B); |
---|
1130 | |
---|
1131 | crossSect = (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); |
---|
1132 | flex *= crossSect; |
---|
1133 | flex *= Pi*rad*rad*L; |
---|
1134 | flex *= cont*cont; |
---|
1135 | flex *= 1.0e8; |
---|
1136 | answer = scale*flex + bkg; |
---|
1137 | |
---|
1138 | return answer; |
---|
1139 | } |
---|
1140 | |
---|
1141 | /* FlexCyl_EllipX : calculates the form factor of a flexible cylinder with an elliptical cross section |
---|
1142 | -- incorporates Wei-Ren Chen's fixes - 2006 |
---|
1143 | |
---|
1144 | */ |
---|
1145 | double |
---|
1146 | FlexCyl_Ellip(double dp[], double q) |
---|
1147 | { |
---|
1148 | double scale,L,B,bkg,rad,qr,cont,ellRatio; |
---|
1149 | double Pi,flex,crossSect,answer; |
---|
1150 | |
---|
1151 | |
---|
1152 | Pi = 4.0*atan(1.0); |
---|
1153 | scale = dp[0]; //make local copies in case memory moves |
---|
1154 | L = dp[1]; |
---|
1155 | B = dp[2]; |
---|
1156 | rad = dp[3]; |
---|
1157 | ellRatio = dp[4]; |
---|
1158 | cont = dp[5]; |
---|
1159 | bkg = dp[6]; |
---|
1160 | |
---|
1161 | qr = q*rad; |
---|
1162 | |
---|
1163 | flex = Sk_WR(q,L,B); |
---|
1164 | |
---|
1165 | crossSect = EllipticalCross_fn(q,rad,(rad*ellRatio)); |
---|
1166 | flex *= crossSect; |
---|
1167 | flex *= Pi*rad*rad*ellRatio*L; |
---|
1168 | flex *= cont*cont; |
---|
1169 | flex *= 1.0e8; |
---|
1170 | answer = scale*flex + bkg; |
---|
1171 | |
---|
1172 | return answer; |
---|
1173 | } |
---|
1174 | |
---|
1175 | double |
---|
1176 | EllipticalCross_fn(double qq, double a, double b) |
---|
1177 | { |
---|
1178 | double uplim,lolim,Pi,summ,arg,zi,yyy,answer; |
---|
1179 | int i,nord=76; |
---|
1180 | |
---|
1181 | Pi = 4.0*atan(1.0); |
---|
1182 | lolim=0.0; |
---|
1183 | uplim=Pi/2.0; |
---|
1184 | summ=0.0; |
---|
1185 | |
---|
1186 | for(i=0;i<nord;i++) { |
---|
1187 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
1188 | arg = qq*sqrt(a*a*sin(zi)*sin(zi)+b*b*cos(zi)*cos(zi)); |
---|
1189 | yyy = pow((2.0 * NR_BessJ1(arg) / arg),2); |
---|
1190 | yyy *= Gauss76Wt[i]; |
---|
1191 | summ += yyy; |
---|
1192 | } |
---|
1193 | answer = (uplim-lolim)/2.0*summ; |
---|
1194 | answer *= 2.0/Pi; |
---|
1195 | return(answer); |
---|
1196 | |
---|
1197 | } |
---|
1198 | /* FlexCyl_PolyLenX : calculates the form factor of a flecible cylinder at the given x-value p->x |
---|
1199 | the cylinder has a polydisperse Length |
---|
1200 | |
---|
1201 | */ |
---|
1202 | double |
---|
1203 | FlexCyl_PolyLen(double dp[], double q) |
---|
1204 | { |
---|
1205 | int i; |
---|
1206 | double scale,radius,length,pd,delrho,bkg,lb; //local variables of coefficient wave |
---|
1207 | int nord=20; //order of integration |
---|
1208 | double uplim,lolim; //upper and lower integration limits |
---|
1209 | double summ,zi,yyy,answer,Vpoly; //running tally of integration |
---|
1210 | double range,zz,Pi; |
---|
1211 | |
---|
1212 | Pi = 4.0*atan(1.0); |
---|
1213 | range = 3.4; |
---|
1214 | |
---|
1215 | summ = 0.0; //initialize intergral |
---|
1216 | scale = dp[0]; //make local copies in case memory moves |
---|
1217 | length = dp[1]; //radius |
---|
1218 | pd = dp[2]; // average length |
---|
1219 | lb = dp[3]; |
---|
1220 | radius = dp[4]; |
---|
1221 | delrho = dp[5]; |
---|
1222 | bkg = dp[6]; |
---|
1223 | |
---|
1224 | zz = (1.0/pd)*(1.0/pd) - 1.0; |
---|
1225 | |
---|
1226 | lolim = length*(1.0-range*pd); //set the upper/lower limits to cover the distribution |
---|
1227 | if(lolim<0) { |
---|
1228 | lolim = 0; |
---|
1229 | } |
---|
1230 | if(pd>0.3) { |
---|
1231 | range = 3.4 + (pd-0.3)*18.0; |
---|
1232 | } |
---|
1233 | uplim = length*(1.0+range*pd); |
---|
1234 | |
---|
1235 | for(i=0;i<nord;i++) { |
---|
1236 | zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
1237 | yyy = Gauss20Wt[i] * FlePolyLen_kernel(q,radius,length,lb,zz,delrho,zi); |
---|
1238 | summ += yyy; |
---|
1239 | } |
---|
1240 | |
---|
1241 | answer = (uplim-lolim)/2.0*summ; |
---|
1242 | //normalize by average cylinder volume (first moment), using the average length |
---|
1243 | Vpoly=Pi*radius*radius*length; |
---|
1244 | answer /= Vpoly; |
---|
1245 | |
---|
1246 | answer *=delrho*delrho; |
---|
1247 | |
---|
1248 | //convert to [cm-1] |
---|
1249 | answer *= 1.0e8; |
---|
1250 | //Scale |
---|
1251 | answer *= scale; |
---|
1252 | // add in the background |
---|
1253 | answer += bkg; |
---|
1254 | |
---|
1255 | return answer; |
---|
1256 | } |
---|
1257 | |
---|
1258 | /* FlexCyl_PolyLenX : calculates the form factor of a flexible cylinder at the given x-value p->x |
---|
1259 | the cylinder has a polydisperse cross sectional radius |
---|
1260 | |
---|
1261 | */ |
---|
1262 | double |
---|
1263 | FlexCyl_PolyRad(double dp[], double q) |
---|
1264 | { |
---|
1265 | int i; |
---|
1266 | double scale,radius,length,pd,delrho,bkg,lb; //local variables of coefficient wave |
---|
1267 | int nord=76; //order of integration |
---|
1268 | double uplim,lolim; //upper and lower integration limits |
---|
1269 | double summ,zi,yyy,answer,Vpoly; //running tally of integration |
---|
1270 | double range,zz,Pi; |
---|
1271 | |
---|
1272 | |
---|
1273 | Pi = 4.0*atan(1.0); |
---|
1274 | range = 3.4; |
---|
1275 | |
---|
1276 | summ = 0.0; //initialize intergral |
---|
1277 | |
---|
1278 | scale = dp[0]; //make local copies in case memory moves |
---|
1279 | length = dp[1]; //radius |
---|
1280 | lb = dp[2]; // average length |
---|
1281 | radius = dp[3]; |
---|
1282 | pd = dp[4]; |
---|
1283 | delrho = dp[5]; |
---|
1284 | bkg = dp[6]; |
---|
1285 | |
---|
1286 | zz = (1.0/pd)*(1.0/pd) - 1.0; |
---|
1287 | |
---|
1288 | lolim = radius*(1.0-range*pd); //set the upper/lower limits to cover the distribution |
---|
1289 | if(lolim<0) { |
---|
1290 | lolim = 0; |
---|
1291 | } |
---|
1292 | if(pd>0.3) { |
---|
1293 | range = 3.4 + (pd-0.3)*18.0; |
---|
1294 | } |
---|
1295 | uplim = radius*(1.0+range*pd); |
---|
1296 | |
---|
1297 | for(i=0;i<nord;i++) { |
---|
1298 | //zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
1299 | //yyy = Gauss20Wt[i] * FlePolyRad_kernel(q,radius,length,lb,zz,delrho,zi); |
---|
1300 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
1301 | yyy = Gauss76Wt[i] * FlePolyRad_kernel(q,radius,length,lb,zz,delrho,zi); |
---|
1302 | summ += yyy; |
---|
1303 | } |
---|
1304 | |
---|
1305 | answer = (uplim-lolim)/2.0*summ; |
---|
1306 | //normalize by average cylinder volume (second moment), using the average radius |
---|
1307 | Vpoly = Pi*radius*radius*length*(zz+2.0)/(zz+1.0); |
---|
1308 | answer /= Vpoly; |
---|
1309 | |
---|
1310 | answer *=delrho*delrho; |
---|
1311 | |
---|
1312 | //convert to [cm-1] |
---|
1313 | answer *= 1.0e8; |
---|
1314 | //Scale |
---|
1315 | answer *= scale; |
---|
1316 | // add in the background |
---|
1317 | answer += bkg; |
---|
1318 | |
---|
1319 | return answer; |
---|
1320 | } |
---|
1321 | |
---|
1322 | /////////functions for WRC implementation of flexible cylinders |
---|
1323 | static double |
---|
1324 | Sk_WR(double q, double L, double b) |
---|
1325 | { |
---|
1326 | // |
---|
1327 | double p1,p2,p1short,p2short,q0,qconnect; |
---|
1328 | double C,epsilon,ans,q0short,Sexvmodify,pi; |
---|
1329 | |
---|
1330 | pi = 4.0*atan(1.0); |
---|
1331 | |
---|
1332 | p1 = 4.12; |
---|
1333 | p2 = 4.42; |
---|
1334 | p1short = 5.36; |
---|
1335 | p2short = 5.62; |
---|
1336 | q0 = 3.1; |
---|
1337 | qconnect = q0/b; |
---|
1338 | // |
---|
1339 | q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0); |
---|
1340 | |
---|
1341 | // |
---|
1342 | if(L/b > 10.0) { |
---|
1343 | C = 3.06/pow((L/b),0.44); |
---|
1344 | epsilon = 0.176; |
---|
1345 | } else { |
---|
1346 | C = 1.0; |
---|
1347 | epsilon = 0.170; |
---|
1348 | } |
---|
1349 | // |
---|
1350 | |
---|
1351 | if( L > 4*b ) { // Longer Chains |
---|
1352 | if (q*b <= 3.1) { //Modified by Yun on Oct. 15, |
---|
1353 | Sexvmodify = Sexvnew(q, L, b); |
---|
1354 | ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) - (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L); |
---|
1355 | } else { //q(i)*b > 3.1 |
---|
1356 | ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) + a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + pi/(q*L); |
---|
1357 | } |
---|
1358 | } else { //L <= 4*b Shorter Chains |
---|
1359 | if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) { |
---|
1360 | if (q*b<=0.01) { |
---|
1361 | ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0; |
---|
1362 | } else { |
---|
1363 | ans = Sdebye1(q,L,b); |
---|
1364 | } |
---|
1365 | } else { //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3) |
---|
1366 | ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) + a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) + pi/(q*L); |
---|
1367 | } |
---|
1368 | } |
---|
1369 | |
---|
1370 | return(ans); |
---|
1371 | //return(a2long(q, L, b, p1, p2, q0)); |
---|
1372 | } |
---|
1373 | |
---|
1374 | //WR named this w (too generic) |
---|
1375 | static double |
---|
1376 | w_WR(double x) |
---|
1377 | { |
---|
1378 | double yy; |
---|
1379 | yy = 0.5*(1 + tanh((x - 1.523)/0.1477)); |
---|
1380 | |
---|
1381 | return (yy); |
---|
1382 | } |
---|
1383 | |
---|
1384 | // |
---|
1385 | static double |
---|
1386 | u1(double q, double L, double b) |
---|
1387 | { |
---|
1388 | double yy; |
---|
1389 | |
---|
1390 | yy = Rgsquareshort(q,L,b)*q*q; |
---|
1391 | |
---|
1392 | return (yy); |
---|
1393 | } |
---|
1394 | |
---|
1395 | // was named u |
---|
1396 | static double |
---|
1397 | u_WR(double q, double L, double b) |
---|
1398 | { |
---|
1399 | double yy; |
---|
1400 | yy = Rgsquare(q,L,b)*q*q; |
---|
1401 | return (yy); |
---|
1402 | } |
---|
1403 | |
---|
1404 | |
---|
1405 | |
---|
1406 | // |
---|
1407 | static double |
---|
1408 | Rgsquarezero(double q, double L, double b) |
---|
1409 | { |
---|
1410 | double yy; |
---|
1411 | yy = (L*b/6.0) * (1.0 - 1.5*(b/L) + 1.5*pow((b/L),2) - 0.75*pow((b/L),3)*(1.0 - exp(-2.0*(L/b)))); |
---|
1412 | |
---|
1413 | return (yy); |
---|
1414 | } |
---|
1415 | |
---|
1416 | // |
---|
1417 | static double |
---|
1418 | Rgsquareshort(double q, double L, double b) |
---|
1419 | { |
---|
1420 | double yy; |
---|
1421 | yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b); |
---|
1422 | |
---|
1423 | return (yy); |
---|
1424 | } |
---|
1425 | |
---|
1426 | // |
---|
1427 | static double |
---|
1428 | Rgsquare(double q, double L, double b) |
---|
1429 | { |
---|
1430 | double yy; |
---|
1431 | yy = AlphaSquare(L/b)*L*b/6.0; |
---|
1432 | |
---|
1433 | return (yy); |
---|
1434 | } |
---|
1435 | |
---|
1436 | // |
---|
1437 | static double |
---|
1438 | AlphaSquare(double x) |
---|
1439 | { |
---|
1440 | double yy; |
---|
1441 | yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); |
---|
1442 | |
---|
1443 | return (yy); |
---|
1444 | } |
---|
1445 | |
---|
1446 | // ?? funciton is not used - but should the log actually be log10??? |
---|
1447 | static double |
---|
1448 | miu(double x) |
---|
1449 | { |
---|
1450 | double yy; |
---|
1451 | yy = (1.0/8.0)*(9.0*x - 2.0 + 2.0*log(1.0 + x)/x)*exp(1.0/2.565*(1.0/x + (1.0 - 1.0/(x*x))*log(1.0 + x))); |
---|
1452 | |
---|
1453 | return (yy); |
---|
1454 | } |
---|
1455 | |
---|
1456 | // |
---|
1457 | static double |
---|
1458 | Sdebye(double q, double L, double b) |
---|
1459 | { |
---|
1460 | double yy; |
---|
1461 | yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2)); |
---|
1462 | |
---|
1463 | return (yy); |
---|
1464 | } |
---|
1465 | |
---|
1466 | // |
---|
1467 | static double |
---|
1468 | Sdebye1(double q, double L, double b) |
---|
1469 | { |
---|
1470 | double yy; |
---|
1471 | yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) ); |
---|
1472 | |
---|
1473 | return (yy); |
---|
1474 | } |
---|
1475 | |
---|
1476 | // |
---|
1477 | static double |
---|
1478 | Sexv(double q, double L, double b) |
---|
1479 | { |
---|
1480 | double yy,C1,C2,C3,miu,Rg2; |
---|
1481 | C1=1.22; |
---|
1482 | C2=0.4288; |
---|
1483 | C3=-1.651; |
---|
1484 | miu = 0.585; |
---|
1485 | |
---|
1486 | Rg2 = Rgsquare(q,L,b); |
---|
1487 | |
---|
1488 | yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) + C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu))); |
---|
1489 | |
---|
1490 | return (yy); |
---|
1491 | } |
---|
1492 | |
---|
1493 | // this must be WR modified version |
---|
1494 | static double |
---|
1495 | Sexvnew(double q, double L, double b) |
---|
1496 | { |
---|
1497 | double yy,C1,C2,C3,miu; |
---|
1498 | double del=1.05,C_star2,Rg2; |
---|
1499 | |
---|
1500 | C1=1.22; |
---|
1501 | C2=0.4288; |
---|
1502 | C3=-1.651; |
---|
1503 | miu = 0.585; |
---|
1504 | |
---|
1505 | //calculating the derivative to decide on the corection (cutoff) term? |
---|
1506 | // I have modified this from WRs original code |
---|
1507 | |
---|
1508 | if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) { |
---|
1509 | C_star2 = 0.0; |
---|
1510 | } else { |
---|
1511 | C_star2 = 1.0; |
---|
1512 | } |
---|
1513 | |
---|
1514 | Rg2 = Rgsquare(q,L,b); |
---|
1515 | |
---|
1516 | yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + C_star2*w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) + C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu))); |
---|
1517 | |
---|
1518 | return (yy); |
---|
1519 | } |
---|
1520 | |
---|
1521 | // these are the messy ones |
---|
1522 | static double |
---|
1523 | a2short(double q, double L, double b, double p1short, double p2short, double q0) |
---|
1524 | { |
---|
1525 | double yy,Rg2_sh; |
---|
1526 | double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p; |
---|
1527 | double pi; |
---|
1528 | |
---|
1529 | E = 2.718281828459045091; |
---|
1530 | pi = 4.0*atan(1.0); |
---|
1531 | Rg2_sh = Rgsquareshort(q,L,b); |
---|
1532 | Rg2_sh2 = Rg2_sh*Rg2_sh; |
---|
1533 | t1 = ((q0*q0*Rg2_sh)/(b*b)); |
---|
1534 | Et1 = pow(E,t1); |
---|
1535 | Emt1 =pow(E,-t1); |
---|
1536 | q02 = q0*q0; |
---|
1537 | q0p = pow(q0,(-4.0 + p2short) ); |
---|
1538 | |
---|
1539 | //E is the number e |
---|
1540 | yy = ((-(1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b*b*b*L - 8.0*b*b*b*Et1*L - 2.0*b*b*b*L*p1short + 2.0*b*b*b*Et1*L*p1short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p1short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p1short*pi*q02*q0*Rg2_sh2))))))); |
---|
1541 | |
---|
1542 | return (yy); |
---|
1543 | } |
---|
1544 | |
---|
1545 | // |
---|
1546 | static double |
---|
1547 | a1short(double q, double L, double b, double p1short, double p2short, double q0) |
---|
1548 | { |
---|
1549 | double yy,Rg2_sh; |
---|
1550 | double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3; |
---|
1551 | double pi; |
---|
1552 | |
---|
1553 | E = 2.718281828459045091; |
---|
1554 | pi = 4.0*atan(1.0); |
---|
1555 | Rg2_sh = Rgsquareshort(q,L,b); |
---|
1556 | Rg2_sh2 = Rg2_sh*Rg2_sh; |
---|
1557 | b3 = b*b*b; |
---|
1558 | t1 = ((q0*q0*Rg2_sh)/(b*b)); |
---|
1559 | Et1 = pow(E,t1); |
---|
1560 | Emt1 =pow(E,-t1); |
---|
1561 | q02 = q0*q0; |
---|
1562 | q0p = pow(q0,(-4.0 + p1short) ); |
---|
1563 | |
---|
1564 | yy = ((1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short + 2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p2short*pi*q02*q0*Rg2_sh2)))))); |
---|
1565 | |
---|
1566 | return(yy); |
---|
1567 | } |
---|
1568 | |
---|
1569 | // this one will be lots of trouble |
---|
1570 | static double |
---|
1571 | a2long(double q, double L, double b, double p1, double p2, double q0) |
---|
1572 | { |
---|
1573 | double yy,C1,C2,C3,C4,C5,miu,C,Rg2; |
---|
1574 | double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi; |
---|
1575 | double E,b2,b3,b4,q02,q03,q04,q05,Rg22; |
---|
1576 | |
---|
1577 | pi = 4.0*atan(1.0); |
---|
1578 | E = 2.718281828459045091; |
---|
1579 | if( L/b > 10.0) { |
---|
1580 | C = 3.06/pow((L/b),0.44); |
---|
1581 | } else { |
---|
1582 | C = 1.0; |
---|
1583 | } |
---|
1584 | |
---|
1585 | C1 = 1.22; |
---|
1586 | C2 = 0.4288; |
---|
1587 | C3 = -1.651; |
---|
1588 | C4 = 1.523; |
---|
1589 | C5 = 0.1477; |
---|
1590 | miu = 0.585; |
---|
1591 | |
---|
1592 | Rg2 = Rgsquare(q,L,b); |
---|
1593 | Rg22 = Rg2*Rg2; |
---|
1594 | b2 = b*b; |
---|
1595 | b3 = b*b*b; |
---|
1596 | b4 = b3*b; |
---|
1597 | q02 = q0*q0; |
---|
1598 | q03 = q0*q0*q0; |
---|
1599 | q04 = q03*q0; |
---|
1600 | q05 = q04*q0; |
---|
1601 | |
---|
1602 | t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) )); |
---|
1603 | |
---|
1604 | t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + (14*b3*pow(E,(-((q02*Rg2)/b2))))/(15*q03*Rg2) + (2*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7*b2)/(15*q02*Rg2)))*Rg2)/b)))/L; |
---|
1605 | |
---|
1606 | t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(2*C5); |
---|
1607 | |
---|
1608 | t4 = (b4*sqrt(Rg2)*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22); |
---|
1609 | |
---|
1610 | t5 = (2*b4*(((2*q0*Rg2)/b - (2*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); |
---|
1611 | |
---|
1612 | t6 = (8*b4*b*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22); |
---|
1613 | |
---|
1614 | t7 = (((-((3*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 3/miu)))/miu)) - (2*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 2/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 1/miu)))/miu)); |
---|
1615 | |
---|
1616 | t8 = ((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); |
---|
1617 | |
---|
1618 | t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7*b2)/(15*q02*Rg2))) + (7*b2)/(15*q02*Rg2))))/L; |
---|
1619 | |
---|
1620 | t10 = (2*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); |
---|
1621 | |
---|
1622 | |
---|
1623 | yy = ((-1*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1)/miu))))*((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))); |
---|
1624 | |
---|
1625 | return (yy); |
---|
1626 | } |
---|
1627 | |
---|
1628 | //need to define this on my own |
---|
1629 | static double |
---|
1630 | sech_WR(double x) |
---|
1631 | { |
---|
1632 | return(1/cosh(x)); |
---|
1633 | } |
---|
1634 | |
---|
1635 | // |
---|
1636 | static double |
---|
1637 | a1long(double q, double L, double b, double p1, double p2, double q0) |
---|
1638 | { |
---|
1639 | double yy,C,C1,C2,C3,C4,C5,miu,Rg2; |
---|
1640 | double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15; |
---|
1641 | double E,pi; |
---|
1642 | double b2,b3,b4,q02,q03,q04,q05,Rg22; |
---|
1643 | |
---|
1644 | pi = 4.0*atan(1.0); |
---|
1645 | E = 2.718281828459045091; |
---|
1646 | |
---|
1647 | if( L/b > 10.0) { |
---|
1648 | C = 3.06/pow((L/b),0.44); |
---|
1649 | } else { |
---|
1650 | C = 1.0; |
---|
1651 | } |
---|
1652 | |
---|
1653 | C1 = 1.22; |
---|
1654 | C2 = 0.4288; |
---|
1655 | C3 = -1.651; |
---|
1656 | C4 = 1.523; |
---|
1657 | C5 = 0.1477; |
---|
1658 | miu = 0.585; |
---|
1659 | |
---|
1660 | Rg2 = Rgsquare(q,L,b); |
---|
1661 | Rg22 = Rg2*Rg2; |
---|
1662 | b2 = b*b; |
---|
1663 | b3 = b*b*b; |
---|
1664 | b4 = b3*b; |
---|
1665 | q02 = q0*q0; |
---|
1666 | q03 = q0*q0*q0; |
---|
1667 | q04 = q03*q0; |
---|
1668 | q05 = q04*q0; |
---|
1669 | |
---|
1670 | t1 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2)))); |
---|
1671 | |
---|
1672 | t2 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); |
---|
1673 | |
---|
1674 | t3 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu)))); |
---|
1675 | |
---|
1676 | t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); |
---|
1677 | |
---|
1678 | t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)))); |
---|
1679 | |
---|
1680 | t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2)))*Rg2)/b))); |
---|
1681 | |
---|
1682 | t7 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); |
---|
1683 | |
---|
1684 | t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); |
---|
1685 | |
---|
1686 | t9 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); |
---|
1687 | |
---|
1688 | t10 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); |
---|
1689 | |
---|
1690 | t11 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu)); |
---|
1691 | |
---|
1692 | t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); |
---|
1693 | |
---|
1694 | t13 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02* Rg2))) + (7.0*b2)/(15.0*q02*Rg2)))); |
---|
1695 | |
---|
1696 | t14 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); |
---|
1697 | |
---|
1698 | t15 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu)))); |
---|
1699 | |
---|
1700 | |
---|
1701 | yy = (pow(q0,p1)*(((-((b*pi)/(L*q0))) +t1/L +t2/(q04*Rg22) + 1.0/2.0*t3*t4)) + (t5*((pow(q0,(p1 - p2))*(((-pow(q0,(-p1)))*(((b2*pi)/(L*q02) +t6/L +t7/(2.0*C5) -t8/(C5*q04*Rg22) +t9/(q04*Rg22) -t10/(q05*Rg22) + 1.0/2.0*t11*t12)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) +t13/L +t14/(q04*Rg22) + 1.0/2.0*t15*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))))); |
---|
1702 | |
---|
1703 | return (yy); |
---|
1704 | } |
---|
1705 | |
---|
1706 | |
---|
1707 | |
---|
1708 | /////////////// |
---|
1709 | |
---|
1710 | // |
---|
1711 | // FUNCTION gfn2: CONTAINS F(Q,A,B,mu)**2 AS GIVEN |
---|
1712 | // BY (53) AND (56,57) IN CHEN AND |
---|
1713 | // KOTLARCHYK REFERENCE |
---|
1714 | // |
---|
1715 | // <PROLATE ELLIPSOIDS> |
---|
1716 | // |
---|
1717 | double |
---|
1718 | gfn2(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq) |
---|
1719 | { |
---|
1720 | // local variables |
---|
1721 | double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,gfn2,pi43,gfn,Pi; |
---|
1722 | |
---|
1723 | Pi = 4.0*atan(1.0); |
---|
1724 | |
---|
1725 | pi43=4.0/3.0*Pi; |
---|
1726 | aa = crmaj; |
---|
1727 | bb = crmin; |
---|
1728 | u2 = (aa*aa*xx*xx + bb*bb*(1.0-xx*xx)); |
---|
1729 | ut2 = (trmaj*trmaj*xx*xx + trmin*trmin*(1.0-xx*xx)); |
---|
1730 | uq = sqrt(u2)*qq; |
---|
1731 | ut= sqrt(ut2)*qq; |
---|
1732 | vc = pi43*aa*bb*bb; |
---|
1733 | vt = pi43*trmaj*trmin*trmin; |
---|
1734 | if (uq == 0){ |
---|
1735 | siq = 1/3; |
---|
1736 | }else{ |
---|
1737 | siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq; |
---|
1738 | } |
---|
1739 | if (ut == 0){ |
---|
1740 | sit = 1/3; |
---|
1741 | }else{ |
---|
1742 | sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut; |
---|
1743 | } |
---|
1744 | gfnc = 3.0*siq*vc*delpc; |
---|
1745 | gfnt = 3.0*sit*vt*delps; |
---|
1746 | gfn = gfnc+gfnt; |
---|
1747 | gfn2 = gfn*gfn; |
---|
1748 | |
---|
1749 | return (gfn2); |
---|
1750 | } |
---|
1751 | |
---|
1752 | // |
---|
1753 | // FUNCTION gfn4: CONTAINS F(Q,A,B,MU)**2 AS GIVEN |
---|
1754 | // BY (53) & (58-59) IN CHEN AND |
---|
1755 | // KOTLARCHYK REFERENCE |
---|
1756 | // |
---|
1757 | // <OBLATE ELLIPSOID> |
---|
1758 | // function gfn4 for oblate ellipsoids |
---|
1759 | double |
---|
1760 | gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq) |
---|
1761 | { |
---|
1762 | // local variables |
---|
1763 | double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,tgfn,gfn4,pi43,Pi; |
---|
1764 | |
---|
1765 | Pi = 4.0*atan(1.0); |
---|
1766 | pi43=4.0/3.0*Pi; |
---|
1767 | aa = crmaj; |
---|
1768 | bb = crmin; |
---|
1769 | u2 = (bb*bb*xx*xx + aa*aa*(1.0-xx*xx)); |
---|
1770 | ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx)); |
---|
1771 | uq = sqrt(u2)*qq; |
---|
1772 | ut= sqrt(ut2)*qq; |
---|
1773 | vc = pi43*aa*aa*bb; |
---|
1774 | vt = pi43*trmaj*trmaj*trmin; |
---|
1775 | if (uq == 0){ |
---|
1776 | siq = 1/3; |
---|
1777 | }else{ |
---|
1778 | siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq; |
---|
1779 | } |
---|
1780 | if (ut == 0){ |
---|
1781 | sit = 1/3; |
---|
1782 | }else{ |
---|
1783 | sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut; |
---|
1784 | } |
---|
1785 | gfnc = 3.0*siq*vc*delpc; |
---|
1786 | gfnt = 3.0*sit*vt*delps; |
---|
1787 | tgfn = gfnc+gfnt; |
---|
1788 | gfn4 = tgfn*tgfn; |
---|
1789 | |
---|
1790 | return (gfn4); |
---|
1791 | } |
---|
1792 | |
---|
1793 | double |
---|
1794 | FlePolyLen_kernel(double q, double radius, double length, double lb, double zz, double delrho, double zi) |
---|
1795 | { |
---|
1796 | double Pq,vcyl,dl; |
---|
1797 | double Pi,qr; |
---|
1798 | |
---|
1799 | Pi = 4.0*atan(1.0); |
---|
1800 | qr = q*radius; |
---|
1801 | |
---|
1802 | Pq = Sk_WR(q,zi,lb); //does not have cross section term |
---|
1803 | if (qr !=0){ |
---|
1804 | Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); |
---|
1805 | } //else Pk *=1; |
---|
1806 | vcyl=Pi*radius*radius*zi; |
---|
1807 | Pq *= vcyl*vcyl; |
---|
1808 | |
---|
1809 | dl = SchulzPoint_cpr(zi,length,zz); |
---|
1810 | return (Pq*dl); |
---|
1811 | |
---|
1812 | } |
---|
1813 | |
---|
1814 | double |
---|
1815 | FlePolyRad_kernel(double q, double ravg, double Lc, double Lb, double zz, double delrho, double zi) |
---|
1816 | { |
---|
1817 | double Pq,vcyl,dr; |
---|
1818 | double Pi,qr; |
---|
1819 | |
---|
1820 | Pi = 4.0*atan(1.0); |
---|
1821 | qr = q*zi; |
---|
1822 | |
---|
1823 | Pq = Sk_WR(q,Lc,Lb); //does not have cross section term |
---|
1824 | if (qr !=0){ |
---|
1825 | Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); |
---|
1826 | } |
---|
1827 | |
---|
1828 | vcyl=Pi*zi*zi*Lc; |
---|
1829 | Pq *= vcyl*vcyl; |
---|
1830 | |
---|
1831 | dr = SchulzPoint_cpr(zi,ravg,zz); |
---|
1832 | return (Pq*dr); |
---|
1833 | |
---|
1834 | } |
---|
1835 | |
---|
1836 | double |
---|
1837 | CSCylIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length) |
---|
1838 | { |
---|
1839 | double answer,halfheight,Pi; |
---|
1840 | double lolim,uplim,summ,yyy,zi; |
---|
1841 | int nord,i; |
---|
1842 | |
---|
1843 | // set up the integration end points |
---|
1844 | Pi = 4.0*atan(1.0); |
---|
1845 | nord = 76; |
---|
1846 | lolim = 0; |
---|
1847 | uplim = Pi/2; |
---|
1848 | halfheight = length/2.0; |
---|
1849 | |
---|
1850 | summ = 0.0; // initialize integral |
---|
1851 | i=0; |
---|
1852 | for(i=0;i<nord;i++) { |
---|
1853 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
1854 | yyy = Gauss76Wt[i] * CScyl(qq, rad, radthick, facthick, rhoc,rhos,rhosolv, halfheight, zi); |
---|
1855 | summ += yyy; |
---|
1856 | } |
---|
1857 | |
---|
1858 | // calculate value of integral to return |
---|
1859 | answer = (uplim-lolim)/2.0*summ; |
---|
1860 | return (answer); |
---|
1861 | } |
---|
1862 | |
---|
1863 | double |
---|
1864 | CScyl(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length, double dum) |
---|
1865 | { |
---|
1866 | // qq is the q-value for the calculation (1/A) |
---|
1867 | // radius is the core radius of the cylinder (A) |
---|
1868 | // radthick and facthick are the radial and face layer thicknesses |
---|
1869 | // rho(n) are the respective SLD's |
---|
1870 | // length is the *Half* CORE-LENGTH of the cylinder |
---|
1871 | // dum is the dummy variable for the integration (theta) |
---|
1872 | |
---|
1873 | double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval; |
---|
1874 | double Pi; |
---|
1875 | |
---|
1876 | Pi = 4.0*atan(1.0); |
---|
1877 | |
---|
1878 | dr1 = rhoc-rhos; |
---|
1879 | dr2 = rhos-rhosolv; |
---|
1880 | vol1 = Pi*rad*rad*(2.0*length); |
---|
1881 | vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick); |
---|
1882 | |
---|
1883 | besarg1 = qq*rad*sin(dum); |
---|
1884 | besarg2 = qq*(rad+radthick)*sin(dum); |
---|
1885 | sinarg1 = qq*length*cos(dum); |
---|
1886 | sinarg2 = qq*(length+facthick)*cos(dum); |
---|
1887 | if (besarg1 == 0){ |
---|
1888 | be1 = 0.5; |
---|
1889 | }else{ |
---|
1890 | be1 = NR_BessJ1(besarg1)/besarg1; |
---|
1891 | } |
---|
1892 | if (besarg2 == 0){ |
---|
1893 | be2 = 0.5; |
---|
1894 | }else{ |
---|
1895 | be2 = NR_BessJ1(besarg2)/besarg2; |
---|
1896 | } |
---|
1897 | if (sinarg1 == 0){ |
---|
1898 | si1 = 1; |
---|
1899 | }else{ |
---|
1900 | si1 = sin(sinarg1)/sinarg1; |
---|
1901 | } |
---|
1902 | if (besarg2 == 0){ |
---|
1903 | si2 = 1; |
---|
1904 | }else{ |
---|
1905 | si2 = sin(sinarg2)/sinarg2; |
---|
1906 | } |
---|
1907 | |
---|
1908 | t1 = 2.0*vol1*dr1*si1*be1; |
---|
1909 | t2 = 2.0*vol2*dr2*si2*be2; |
---|
1910 | |
---|
1911 | retval = ((t1+t2)*(t1+t2))*sin(dum); |
---|
1912 | return (retval); |
---|
1913 | |
---|
1914 | } |
---|
1915 | |
---|
1916 | |
---|
1917 | double |
---|
1918 | CoreShellCylKernel(double qq, double rcore, double thick, double rhoc, double rhos, double rhosolv, double length, double dum) |
---|
1919 | { |
---|
1920 | |
---|
1921 | double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval; |
---|
1922 | double Pi; |
---|
1923 | |
---|
1924 | Pi = 4.0*atan(1.0); |
---|
1925 | |
---|
1926 | dr1 = rhoc-rhos; |
---|
1927 | dr2 = rhos-rhosolv; |
---|
1928 | vol1 = Pi*rcore*rcore*(2.0*length); |
---|
1929 | vol2 = Pi*(rcore+thick)*(rcore+thick)*(2.0*length+2.0*thick); |
---|
1930 | |
---|
1931 | besarg1 = qq*rcore*sin(dum); |
---|
1932 | besarg2 = qq*(rcore+thick)*sin(dum); |
---|
1933 | sinarg1 = qq*length*cos(dum); |
---|
1934 | sinarg2 = qq*(length+thick)*cos(dum); |
---|
1935 | |
---|
1936 | if (besarg1 == 0){ |
---|
1937 | be1 = 0.5; |
---|
1938 | }else{ |
---|
1939 | be1 = NR_BessJ1(besarg1)/besarg1; |
---|
1940 | } |
---|
1941 | if (besarg2 == 0){ |
---|
1942 | be2 = 0.5; |
---|
1943 | }else{ |
---|
1944 | be2 = NR_BessJ1(besarg2)/besarg2; |
---|
1945 | } |
---|
1946 | if (sinarg1 == 0){ |
---|
1947 | si1 = 1; |
---|
1948 | }else{ |
---|
1949 | si1 = sin(sinarg1)/sinarg1; |
---|
1950 | } |
---|
1951 | if (besarg2 == 0){ |
---|
1952 | si2 = 1; |
---|
1953 | }else{ |
---|
1954 | si2 = sin(sinarg2)/sinarg2; |
---|
1955 | } |
---|
1956 | |
---|
1957 | t1 = 2.0*vol1*dr1*si1*be1; |
---|
1958 | t2 = 2.0*vol2*dr2*si2*be2; |
---|
1959 | |
---|
1960 | retval = ((t1+t2)*(t1+t2))*sin(dum); |
---|
1961 | |
---|
1962 | return (retval); |
---|
1963 | } |
---|
1964 | |
---|
1965 | double |
---|
1966 | Cyl_PolyLenKernel(double q, double radius, double len_avg, double zz, double delrho, double dumLen) |
---|
1967 | { |
---|
1968 | |
---|
1969 | double halfheight,uplim,lolim,zi,summ,yyy,Pi; |
---|
1970 | double answer,dr,Vcyl; |
---|
1971 | int i,nord; |
---|
1972 | |
---|
1973 | Pi = 4.0*atan(1.0); |
---|
1974 | lolim = 0; |
---|
1975 | uplim = Pi/2.0; |
---|
1976 | halfheight = dumLen/2.0; |
---|
1977 | nord=20; |
---|
1978 | summ = 0.0; |
---|
1979 | |
---|
1980 | //do the cylinder orientational average |
---|
1981 | for(i=0;i<nord;i++) { |
---|
1982 | zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
1983 | yyy = Gauss20Wt[i] * CylKernel(q, radius, halfheight, zi); |
---|
1984 | summ += yyy; |
---|
1985 | } |
---|
1986 | answer = (uplim-lolim)/2.0*summ; |
---|
1987 | // Multiply by contrast^2 |
---|
1988 | answer *= delrho*delrho; |
---|
1989 | // don't do the normal scaling to volume here |
---|
1990 | // instead, multiply by VCyl^2 to get rid of the normalization for this radius of cylinder |
---|
1991 | Vcyl = Pi*radius*radius*dumLen; |
---|
1992 | answer *= Vcyl*Vcyl; |
---|
1993 | |
---|
1994 | dr = SchulzPoint_cpr(dumLen,len_avg,zz); |
---|
1995 | return(dr*answer); |
---|
1996 | } |
---|
1997 | |
---|
1998 | |
---|
1999 | double |
---|
2000 | Stackdisc_kern(double qq, double rcore, double rhoc, double rhol, double rhosolv, double length, double thick, double dum, double gsd, double d, double N) |
---|
2001 | { |
---|
2002 | // qq is the q-value for the calculation (1/A) |
---|
2003 | // rcore is the core radius of the cylinder (A) |
---|
2004 | // rho(n) are the respective SLD's |
---|
2005 | // length is the *Half* CORE-LENGTH of the cylinder = L (A) |
---|
2006 | // dum is the dummy variable for the integration (x in Feigin's notation) |
---|
2007 | |
---|
2008 | //Local variables |
---|
2009 | double totald,dr1,dr2,besarg1,besarg2,be1,be2,si1,si2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt; |
---|
2010 | double Pi; |
---|
2011 | int kk; |
---|
2012 | |
---|
2013 | Pi = 4.0*atan(1.0); |
---|
2014 | |
---|
2015 | dr1 = rhoc-rhosolv; |
---|
2016 | dr2 = rhol-rhosolv; |
---|
2017 | area = Pi*rcore*rcore; |
---|
2018 | totald=2.0*(thick+length); |
---|
2019 | |
---|
2020 | besarg1 = qq*rcore*sin(dum); |
---|
2021 | besarg2 = qq*rcore*sin(dum); |
---|
2022 | |
---|
2023 | sinarg1 = qq*length*cos(dum); |
---|
2024 | sinarg2 = qq*(length+thick)*cos(dum); |
---|
2025 | |
---|
2026 | if (besarg1 == 0){ |
---|
2027 | be1 = 0.5; |
---|
2028 | }else{ |
---|
2029 | be1 = NR_BessJ1(besarg1)/besarg1; |
---|
2030 | } |
---|
2031 | if (besarg2 == 0){ |
---|
2032 | be2 = 0.5; |
---|
2033 | }else{ |
---|
2034 | be2 = NR_BessJ1(besarg2)/besarg2; |
---|
2035 | } |
---|
2036 | if (sinarg1 == 0){ |
---|
2037 | si1 = 1; |
---|
2038 | }else{ |
---|
2039 | si1 = sin(sinarg1)/sinarg1; |
---|
2040 | } |
---|
2041 | if (besarg2 == 0){ |
---|
2042 | si2 = 1; |
---|
2043 | }else{ |
---|
2044 | si2 = sin(sinarg2)/sinarg2; |
---|
2045 | } |
---|
2046 | |
---|
2047 | t1 = 2*area*(2*length)*dr1*(si1)*(be1); |
---|
2048 | t2 = 2*area*dr2*(totald*si2-2*length*si1)*(be2); |
---|
2049 | |
---|
2050 | retval =((t1+t2)*(t1+t2))*sin(dum); |
---|
2051 | |
---|
2052 | // loop for the structure facture S(q) |
---|
2053 | sqq=0.0; |
---|
2054 | for(kk=1;kk<N;kk+=1) { |
---|
2055 | dexpt=qq*cos(dum)*qq*cos(dum)*d*d*gsd*gsd*kk/2.0; |
---|
2056 | sqq=sqq+(N-kk)*cos(qq*cos(dum)*d*kk)*exp(-1.*dexpt); |
---|
2057 | } |
---|
2058 | |
---|
2059 | // end of loop for S(q) |
---|
2060 | sqq=1.0+2.0*sqq/N; |
---|
2061 | |
---|
2062 | retval *= sqq; |
---|
2063 | |
---|
2064 | return(retval); |
---|
2065 | } |
---|
2066 | |
---|
2067 | |
---|
2068 | double |
---|
2069 | Cyl_PolyRadKernel(double q, double radius, double length, double zz, double delrho, double dumRad) |
---|
2070 | { |
---|
2071 | |
---|
2072 | double halfheight,uplim,lolim,zi,summ,yyy,Pi; |
---|
2073 | double answer,dr,Vcyl; |
---|
2074 | int i,nord; |
---|
2075 | |
---|
2076 | Pi = 4.0*atan(1.0); |
---|
2077 | lolim = 0; |
---|
2078 | uplim = Pi/2.0; |
---|
2079 | halfheight = length/2.0; |
---|
2080 | // nord=20; |
---|
2081 | nord=76; |
---|
2082 | summ = 0.0; |
---|
2083 | |
---|
2084 | //do the cylinder orientational average |
---|
2085 | // for(i=0;i<nord;i++) { |
---|
2086 | // zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
2087 | // yyy = Gauss20Wt[i] * CylKernel(q, dumRad, halfheight, zi); |
---|
2088 | // summ += yyy; |
---|
2089 | // } |
---|
2090 | for(i=0;i<nord;i++) { |
---|
2091 | zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; |
---|
2092 | yyy = Gauss76Wt[i] * CylKernel(q, dumRad, halfheight, zi); |
---|
2093 | summ += yyy; |
---|
2094 | } |
---|
2095 | answer = (uplim-lolim)/2.0*summ; |
---|
2096 | // Multiply by contrast^2 |
---|
2097 | answer *= delrho*delrho; |
---|
2098 | // don't do the normal scaling to volume here |
---|
2099 | // instead, multiply by VCyl^2 to get rid of the normalization for this radius of cylinder |
---|
2100 | Vcyl = Pi*dumRad*dumRad*length; |
---|
2101 | answer *= Vcyl*Vcyl; |
---|
2102 | |
---|
2103 | dr = SchulzPoint_cpr(dumRad,radius,zz); |
---|
2104 | return(dr*answer); |
---|
2105 | } |
---|
2106 | |
---|
2107 | double |
---|
2108 | SchulzPoint_cpr(double dumRad, double radius, double zz) |
---|
2109 | { |
---|
2110 | double dr; |
---|
2111 | |
---|
2112 | dr = zz*log(dumRad) - gammaln(zz+1.0) + (zz+1.0)*log((zz+1.0)/radius)-(dumRad/radius*(zz+1.0)); |
---|
2113 | return(exp(dr)); |
---|
2114 | } |
---|
2115 | |
---|
2116 | static double |
---|
2117 | gammaln(double xx) |
---|
2118 | { |
---|
2119 | double x,y,tmp,ser; |
---|
2120 | static double cof[6]={76.18009172947146,-86.50532032941677, |
---|
2121 | 24.01409824083091,-1.231739572450155, |
---|
2122 | 0.1208650973866179e-2,-0.5395239384953e-5}; |
---|
2123 | int j; |
---|
2124 | |
---|
2125 | y=x=xx; |
---|
2126 | tmp=x+5.5; |
---|
2127 | tmp -= (x+0.5)*log(tmp); |
---|
2128 | ser=1.000000000190015; |
---|
2129 | for (j=0;j<=5;j++) ser += cof[j]/++y; |
---|
2130 | return -tmp+log(2.5066282746310005*ser/x); |
---|
2131 | } |
---|
2132 | |
---|
2133 | |
---|
2134 | double |
---|
2135 | EllipsoidKernel(double qq, double a, double nua, double dum) |
---|
2136 | { |
---|
2137 | double arg,nu,retval; //local variables |
---|
2138 | |
---|
2139 | nu = nua/a; |
---|
2140 | arg = qq*a*sqrt(1+dum*dum*(nu*nu-1)); |
---|
2141 | if (arg == 0){ |
---|
2142 | retval =1/3; |
---|
2143 | }else{ |
---|
2144 | retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg); |
---|
2145 | } |
---|
2146 | retval *= retval; |
---|
2147 | retval *= 9; |
---|
2148 | |
---|
2149 | return(retval); |
---|
2150 | }//Function EllipsoidKernel() |
---|
2151 | |
---|
2152 | double |
---|
2153 | HolCylKernel(double qq, double rcore, double rshell, double length, double dum) |
---|
2154 | { |
---|
2155 | double gamma,arg1,arg2,lam1,lam2,psi,sinarg,t2,retval; //local variables |
---|
2156 | |
---|
2157 | gamma = rcore/rshell; |
---|
2158 | arg1 = qq*rshell*sqrt(1-dum*dum); //1=shell (outer radius) |
---|
2159 | arg2 = qq*rcore*sqrt(1-dum*dum); //2=core (inner radius) |
---|
2160 | if (arg1 == 0){ |
---|
2161 | lam1 = 1; |
---|
2162 | }else{ |
---|
2163 | lam1 = 2*NR_BessJ1(arg1)/arg1; |
---|
2164 | } |
---|
2165 | if (arg2 == 0){ |
---|
2166 | lam2 = 1; |
---|
2167 | }else{ |
---|
2168 | lam2 = 2*NR_BessJ1(arg2)/arg2; |
---|
2169 | } |
---|
2170 | //Todo: Need to check psi behavior as gamma goes to 1. |
---|
2171 | psi = 1/(1-gamma*gamma)*(lam1 - gamma*gamma*lam2); //SRK 10/19/00 |
---|
2172 | sinarg = qq*length*dum/2; |
---|
2173 | if (sinarg == 0){ |
---|
2174 | t2 = 1; |
---|
2175 | }else{ |
---|
2176 | t2 = sin(sinarg)/sinarg; |
---|
2177 | } |
---|
2178 | |
---|
2179 | retval = psi*psi*t2*t2; |
---|
2180 | |
---|
2181 | return(retval); |
---|
2182 | }//Function HolCylKernel() |
---|
2183 | |
---|
2184 | double |
---|
2185 | PPKernel(double aa, double mu, double uu) |
---|
2186 | { |
---|
2187 | // mu passed in is really mu*sqrt(1-sig^2) |
---|
2188 | double arg1,arg2,Pi,tmp1,tmp2; //local variables |
---|
2189 | |
---|
2190 | Pi = 4.0*atan(1.0); |
---|
2191 | |
---|
2192 | //handle arg=0 separately, as sin(t)/t -> 1 as t->0 |
---|
2193 | arg1 = (mu/2)*cos(Pi*uu/2); |
---|
2194 | arg2 = (mu*aa/2)*sin(Pi*uu/2); |
---|
2195 | if(arg1==0) { |
---|
2196 | tmp1 = 1; |
---|
2197 | } else { |
---|
2198 | tmp1 = sin(arg1)*sin(arg1)/arg1/arg1; |
---|
2199 | } |
---|
2200 | |
---|
2201 | if (arg2==0) { |
---|
2202 | tmp2 = 1; |
---|
2203 | } else { |
---|
2204 | tmp2 = sin(arg2)*sin(arg2)/arg2/arg2; |
---|
2205 | } |
---|
2206 | |
---|
2207 | return (tmp1*tmp2); |
---|
2208 | |
---|
2209 | }//Function PPKernel() |
---|
2210 | |
---|
2211 | |
---|
2212 | double |
---|
2213 | TriaxialKernel(double q, double aa, double bb, double cc, double dx, double dy) |
---|
2214 | { |
---|
2215 | |
---|
2216 | double arg,val,pi; //local variables |
---|
2217 | |
---|
2218 | pi = 4.0*atan(1.0); |
---|
2219 | |
---|
2220 | arg = aa*aa*cos(pi*dx/2)*cos(pi*dx/2); |
---|
2221 | arg += bb*bb*sin(pi*dx/2)*sin(pi*dx/2)*(1-dy*dy); |
---|
2222 | arg += cc*cc*dy*dy; |
---|
2223 | arg = q*sqrt(arg); |
---|
2224 | if (arg == 0){ |
---|
2225 | val = 1; // as arg --> 0, val should go to 1.0 |
---|
2226 | }else{ |
---|
2227 | val = 9 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ); |
---|
2228 | } |
---|
2229 | return (val); |
---|
2230 | |
---|
2231 | }//Function TriaxialKernel() |
---|
2232 | |
---|
2233 | |
---|
2234 | double |
---|
2235 | CylKernel(double qq, double rr,double h, double theta) |
---|
2236 | { |
---|
2237 | |
---|
2238 | // qq is the q-value for the calculation (1/A) |
---|
2239 | // rr is the radius of the cylinder (A) |
---|
2240 | // h is the HALF-LENGTH of the cylinder = L/2 (A) |
---|
2241 | |
---|
2242 | double besarg,bj,retval,d1,t1,b1,t2,b2,siarg,be,si; //Local variables |
---|
2243 | |
---|
2244 | |
---|
2245 | besarg = qq*rr*sin(theta); |
---|
2246 | siarg = qq * h * cos(theta); |
---|
2247 | bj =NR_BessJ1(besarg); |
---|
2248 | |
---|
2249 | //* Computing 2nd power */ |
---|
2250 | d1 = sin(siarg); |
---|
2251 | t1 = d1 * d1; |
---|
2252 | //* Computing 2nd power */ |
---|
2253 | d1 = bj; |
---|
2254 | t2 = d1 * d1 * 4.0 * sin(theta); |
---|
2255 | //* Computing 2nd power */ |
---|
2256 | d1 = siarg; |
---|
2257 | b1 = d1 * d1; |
---|
2258 | //* Computing 2nd power */ |
---|
2259 | d1 = qq * rr * sin(theta); |
---|
2260 | b2 = d1 * d1; |
---|
2261 | if (besarg == 0){ |
---|
2262 | be = sin(theta); |
---|
2263 | }else{ |
---|
2264 | be = t2 / b2; |
---|
2265 | } |
---|
2266 | if (siarg == 0){ |
---|
2267 | si = 1; |
---|
2268 | }else{ |
---|
2269 | si = t1 / b1; |
---|
2270 | } |
---|
2271 | retval = be * si; |
---|
2272 | |
---|
2273 | return (retval); |
---|
2274 | |
---|
2275 | }//Function CylKernel() |
---|
2276 | |
---|
2277 | double |
---|
2278 | EllipCylKernel(double qq, double ra,double nu, double theta) |
---|
2279 | { |
---|
2280 | //this is the function LAMBDA1^2 in Feigin's notation |
---|
2281 | // qq is the q-value for the calculation (1/A) |
---|
2282 | // ra is the transformed radius"a" in Feigin's notation |
---|
2283 | // nu is the ratio (major radius)/(minor radius) of the Ellipsoid [=] --- |
---|
2284 | // theta is the dummy variable of the integration |
---|
2285 | |
---|
2286 | double retval,arg; //Local variables |
---|
2287 | |
---|
2288 | arg = qq*ra*sqrt((1+nu*nu)/2+(1-nu*nu)*cos(theta)/2); |
---|
2289 | if (arg == 0){ |
---|
2290 | retval = 1; |
---|
2291 | }else{ |
---|
2292 | retval = 2*NR_BessJ1(arg)/arg; |
---|
2293 | } |
---|
2294 | |
---|
2295 | //square it |
---|
2296 | retval *= retval; |
---|
2297 | |
---|
2298 | return(retval); |
---|
2299 | |
---|
2300 | }//Function EllipCylKernel() |
---|
2301 | |
---|
2302 | double NR_BessJ1(double x) |
---|
2303 | { |
---|
2304 | double ax,z; |
---|
2305 | double xx,y,ans,ans1,ans2; |
---|
2306 | |
---|
2307 | if ((ax=fabs(x)) < 8.0) { |
---|
2308 | y=x*x; |
---|
2309 | ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1 |
---|
2310 | +y*(-2972611.439+y*(15704.48260+y*(-30.16036606)))))); |
---|
2311 | ans2=144725228442.0+y*(2300535178.0+y*(18583304.74 |
---|
2312 | +y*(99447.43394+y*(376.9991397+y*1.0)))); |
---|
2313 | ans=ans1/ans2; |
---|
2314 | } else { |
---|
2315 | z=8.0/ax; |
---|
2316 | y=z*z; |
---|
2317 | xx=ax-2.356194491; |
---|
2318 | ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4 |
---|
2319 | +y*(0.2457520174e-5+y*(-0.240337019e-6)))); |
---|
2320 | ans2=0.04687499995+y*(-0.2002690873e-3 |
---|
2321 | +y*(0.8449199096e-5+y*(-0.88228987e-6 |
---|
2322 | +y*0.105787412e-6))); |
---|
2323 | ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2); |
---|
2324 | if (x < 0.0) ans = -ans; |
---|
2325 | } |
---|
2326 | |
---|
2327 | return(ans); |
---|
2328 | } |
---|