source: sasview/src/sas/models/c_extension/libigor/libCylinder.c @ 8d302cd

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 8d302cd was 8d302cd, checked in by richardh, 9 years ago

this time with capital C in libCylinder.c for new lammellarPS and lamellarPS_HG

  • Property mode set to 100644
File size: 89.0 KB
Line 
1/*      CylinderFit.c
2
3A simplified project designed to act as a template for your curve fitting function.
4The 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
11/////////functions for WRC implementation of flexible cylinders
12static double
13gammaln(double xx)
14{
15  double x,y,tmp,ser;
16  static double cof[6]={76.18009172947146,-86.50532032941677,
17    24.01409824083091,-1.231739572450155,
18    0.1208650973866179e-2,-0.5395239384953e-5};
19  int j;
20
21  y=x=xx;
22  tmp=x+5.5;
23  tmp -= (x+0.5)*log(tmp);
24  ser=1.000000000190015;
25  for (j=0;j<=5;j++) ser += cof[j]/++y;
26  return -tmp+log(2.5066282746310005*ser/x);
27}
28/////////functions for WRC implementation of flexible cylinders
29
30//
31static double
32AlphaSquare(double x)
33{
34    double yy;
35    yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) );
36
37    return (yy);
38}
39
40//
41static double
42Rgsquarezero(double q, double L, double b)
43{
44    double yy;
45    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))));
46
47    return (yy);
48}
49
50//
51static double
52Rgsquareshort(double q, double L, double b)
53{
54    double yy;
55    yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b);
56
57    return (yy);
58}
59
60//
61static double
62Rgsquare(double q, double L, double b)
63{
64    double yy;
65    yy = AlphaSquare(L/b)*L*b/6.0;
66
67    return (yy);
68}
69
70// ?? funciton is not used - but should the log actually be log10???
71/*
72static double
73miu(double x)
74{
75    double yy;
76    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)));
77
78    return (yy);
79}
80*/
81
82//WR named this w (too generic)
83static double
84w_WR(double x)
85{
86    double yy;
87    yy = 0.5*(1 + tanh((x - 1.523)/0.1477));
88
89    return (yy);
90}
91
92//
93static double
94u1(double q, double L, double b)
95{
96    double yy;
97
98    yy = Rgsquareshort(q,L,b)*q*q;
99
100    return (yy);
101}
102
103// was named u
104static double
105u_WR(double q, double L, double b)
106{
107    double yy;
108    yy = Rgsquare(q,L,b)*q*q;
109    return (yy);
110}
111
112
113//
114static double
115Sdebye1(double q, double L, double b)
116{
117    double yy;
118    yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) );
119
120    return (yy);
121}
122
123
124//
125static double
126Sdebye(double q, double L, double b)
127{
128    double yy;
129    yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2));
130
131    return (yy);
132}
133
134//
135static double
136Sexv(double q, double L, double b)
137{
138    double yy,C1,C2,C3,miu,Rg2;
139    C1=1.22;
140    C2=0.4288;
141    C3=-1.651;
142    miu = 0.585;
143
144    Rg2 = Rgsquare(q,L,b);
145
146    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)));
147
148    return (yy);
149}
150
151// this must be WR modified version
152static double
153Sexvnew(double q, double L, double b)
154{
155    double yy,C1,C2,C3,miu;
156    double del=1.05,C_star2,Rg2;
157
158    C1=1.22;
159    C2=0.4288;
160    C3=-1.651;
161    miu = 0.585;
162
163    //calculating the derivative to decide on the corection (cutoff) term?
164    // I have modified this from WRs original code
165
166    if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) {
167        C_star2 = 0.0;
168    } else {
169        C_star2 = 1.0;
170    }
171
172    Rg2 = Rgsquare(q,L,b);
173
174  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)));
175
176    return (yy);
177}
178
179// these are the messy ones
180static double
181a2short(double q, double L, double b, double p1short, double p2short, double q0)
182{
183    double yy,Rg2_sh;
184    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p;
185    double pi;
186
187    E = 2.718281828459045091;
188  pi = 4.0*atan(1.0);
189    Rg2_sh = Rgsquareshort(q,L,b);
190    Rg2_sh2 = Rg2_sh*Rg2_sh;
191    t1 = ((q0*q0*Rg2_sh)/(b*b));
192    Et1 = pow(E,t1);
193    Emt1 =pow(E,-t1);
194    q02 = q0*q0;
195    q0p = pow(q0,(-4.0 + p2short) );
196
197    //E is the number e
198    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)))))));
199
200    return (yy);
201}
202
203//
204static double
205a1short(double q, double L, double b, double p1short, double p2short, double q0)
206{
207    double yy,Rg2_sh;
208    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3;
209  double pi;
210
211    E = 2.718281828459045091;
212  pi = 4.0*atan(1.0);
213    Rg2_sh = Rgsquareshort(q,L,b);
214    Rg2_sh2 = Rg2_sh*Rg2_sh;
215    b3 = b*b*b;
216    t1 = ((q0*q0*Rg2_sh)/(b*b));
217    Et1 = pow(E,t1);
218    Emt1 =pow(E,-t1);
219    q02 = q0*q0;
220    q0p = pow(q0,(-4.0 + p1short) );
221
222    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))))));
223
224    return(yy);
225}
226
227
228//need to define this on my own
229static double
230sech_WR(double x)
231{
232  return(1/cosh(x));
233}
234
235// this one will be lots of trouble
236static double
237a2long(double q, double L, double b, double p1, double p2, double q0)
238{
239    double yy,C1,C2,C3,C4,C5,miu,C,Rg2;
240    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi;
241    double E,b2,b3,b4,q02,q03,q04,q05,Rg22;
242
243    pi = 4.0*atan(1.0);
244    E = 2.718281828459045091;
245    if( L/b > 10.0) {
246        C = 3.06/pow((L/b),0.44);
247    } else {
248        C = 1.0;
249    }
250
251    C1 = 1.22;
252    C2 = 0.4288;
253    C3 = -1.651;
254    C4 = 1.523;
255    C5 = 0.1477;
256    miu = 0.585;
257
258    Rg2 = Rgsquare(q,L,b);
259    Rg22 = Rg2*Rg2;
260    b2 = b*b;
261    b3 = b*b*b;
262    b4 = b3*b;
263    q02 = q0*q0;
264    q03 = q0*q0*q0;
265    q04 = q03*q0;
266    q05 = q04*q0;
267
268    t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) ));
269
270    t2 = (b*C*(((-1.0*((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*b2)/(15.0*q02*Rg2)))*Rg2)/b)))/L;
271
272    t3 = (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.0))/(2.0*C5);
273
274    t4 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22);
275
276    t5 = (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))))))/(q04*Rg22);
277
278    t6 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22);
279
280    t7 = (((-((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));
281
282    t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
283
284    t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))))/L;
285
286    t10 = (2.0*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);
287
288
289    yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((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))))*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))));
290
291    return (yy);
292}
293//
294static double
295a1long(double q, double L, double b, double p1, double p2, double q0)
296{
297    double yy,C,C1,C2,C3,C4,C5,miu,Rg2;
298    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15;
299    double E,pi;
300    double b2,b3,b4,q02,q03,q04,q05,Rg22;
301
302    pi = 4.0*atan(1.0);
303    E = 2.718281828459045091;
304
305    if( L/b > 10.0) {
306        C = 3.06/pow((L/b),0.44);
307    } else {
308        C = 1.0;
309    }
310
311    C1 = 1.22;
312    C2 = 0.4288;
313    C3 = -1.651;
314    C4 = 1.523;
315    C5 = 0.1477;
316    miu = 0.585;
317
318    Rg2 = Rgsquare(q,L,b);
319    Rg22 = Rg2*Rg2;
320    b2 = b*b;
321    b3 = b*b*b;
322    b4 = b3*b;
323    q02 = q0*q0;
324    q03 = q0*q0*q0;
325    q04 = q03*q0;
326    q05 = q04*q0;
327
328    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))));
329
330    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))))));
331
332    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))));
333
334    t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
335
336    t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2))));
337
338    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)));
339
340    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));
341
342    t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));
343
344    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))))));
345
346    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))))));
347
348    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));
349
350    t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
351
352    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))));
353
354    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))))));
355
356    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))));
357
358
359    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)))))))))));
360
361    return (yy);
362}
363
364
365static double
366Sk_WR(double q, double L, double b)
367{
368  //
369  double p1,p2,p1short,p2short,q0;
370  double C,ans,q0short,Sexvmodify,pi;
371
372  pi = 4.0*atan(1.0);
373
374  p1 = 4.12;
375  p2 = 4.42;
376  p1short = 5.36;
377  p2short = 5.62;
378  q0 = 3.1;
379  //
380  q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0);
381
382  //
383  if(L/b > 10.0) {
384    C = 3.06/pow((L/b),0.44);
385  } else {
386    C = 1.0;
387  }
388  //
389
390  if( L > 4*b ) { // Longer Chains
391    if (q*b <= 3.1) {   //Modified by Yun on Oct. 15,
392      Sexvmodify = Sexvnew(q, L, b);
393      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);
394    } else { //q(i)*b > 3.1
395      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);
396    }
397  } else { //L <= 4*b Shorter Chains
398    if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) {
399      if (q*b<=0.01) {
400        ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0;
401      } else {
402        ans = Sdebye1(q,L,b);
403      }
404    } else {  //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3)
405      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);
406    }
407  }
408
409  return(ans);
410  //return(a2long(q, L, b, p1, p2, q0));
411}
412
413/*      CylinderForm  :  calculates the form factor of a cylinder at the give x-value p->x
414
415Warning:
416The call to WaveData() below returns a pointer to the middle
417of an unlocked Macintosh handle. In the unlikely event that your
418calculations could cause memory to move, you should copy the coefficient
419values to local variables or an array before such operations.
420*/
421double
422CylinderForm(double dp[], double q)
423{
424
425        int i;
426        double Pi;
427        double scale,radius,length,delrho,bkg,halfheight,sldCyl,sldSolv;                //local variables of coefficient wave
428        int nord=76;                    //order of integration
429        double uplim,lolim;             //upper and lower integration limits
430        double summ,zi,yyy,answer,vcyl;                 //running tally of integration
431       
432        Pi = 4.0*atan(1.0);
433        lolim = 0.0;
434        uplim = Pi/2.0;
435       
436        summ = 0.0;                     //initialize intergral
437       
438        scale = dp[0];                  //make local copies in case memory moves
439        radius = dp[1];
440        length = dp[2];
441        sldCyl = dp[3];
442        sldSolv = dp[4];
443        bkg = dp[5];
444       
445        delrho = sldCyl-sldSolv;
446        halfheight = length/2.0;
447        for(i=0;i<nord;i++) {
448                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
449                yyy = Gauss76Wt[i] * CylKernel(q, radius, halfheight, zi);
450                summ += yyy;
451        }
452       
453        answer = (uplim-lolim)/2.0*summ;
454        // Multiply by contrast^2
455        answer *= delrho*delrho;
456        //normalize by cylinder volume
457        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
458        vcyl=Pi*radius*radius*length;
459        answer *= vcyl;
460        //convert to [cm-1]
461        answer *= 1.0e8;
462        //Scale
463        answer *= scale;
464        // add in the background
465        answer += bkg;
466       
467        return answer;
468}
469
470/*      EllipCyl76X  :  calculates the form factor of a elliptical cylinder at the given x-value p->x
471
472Uses 76 pt Gaussian quadrature for both integrals
473
474Warning:
475The call to WaveData() below returns a pointer to the middle
476of an unlocked Macintosh handle. In the unlikely event that your
477calculations could cause memory to move, you should copy the coefficient
478values to local variables or an array before such operations.
479*/
480double
481EllipCyl76(double dp[], double q)
482{
483        int i,j;
484        double Pi,slde,sld;
485        double scale,ra,nu,length,delrho,bkg;           //local variables of coefficient wave
486        int nord=76;                    //order of integration
487        double va,vb;           //upper and lower integration limits
488        double summ,zi,yyy,answer,vell;                 //running tally of integration
489        double summj,vaj,vbj,zij,arg, si;                       //for the inner integration
490
491        Pi = 4.0*atan(1.0);
492        va = 0.0;
493        vb = 1.0;               //orintational average, outer integral
494        vaj=0.0;
495        vbj=Pi;         //endpoints of inner integral
496       
497        summ = 0.0;                     //initialize intergral
498       
499        scale = dp[0];                  //make local copies in case memory moves
500        ra = dp[1];
501        nu = dp[2];
502        length = dp[3];
503        slde = dp[4];
504        sld = dp[5];
505        delrho = slde - sld;
506        bkg = dp[6];
507       
508        for(i=0;i<nord;i++) {
509                //setup inner integral over the ellipsoidal cross-section
510                summj=0;
511                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy
512                arg = ra*sqrt(1.0-zi*zi);
513                for(j=0;j<nord;j++) {
514                        //76 gauss points for the inner integral as well
515                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "y" dummy
516                        yyy = Gauss76Wt[j] * EllipCylKernel(q,arg,nu,zij);
517                        summj += yyy;
518                }
519                //now calculate the value of the inner integral
520                answer = (vbj-vaj)/2.0*summj;
521                //divide integral by Pi
522                answer /=Pi;
523               
524                //now calculate outer integral
525                arg = q*length*zi/2.0;
526                if (arg == 0.0){
527                        si = 1.0;
528                }else{
529                        si = sin(arg) * sin(arg) / arg / arg;
530                }
531                yyy = Gauss76Wt[i] * answer * si;
532                summ += yyy;
533        }
534        answer = (vb-va)/2.0*summ;
535        // Multiply by contrast^2
536        answer *= delrho*delrho;
537        //normalize by cylinder volume
538        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
539        vell = Pi*ra*(nu*ra)*length;
540        answer *= vell;
541        //convert to [cm-1]
542        answer *= 1.0e8;
543        //Scale
544        answer *= scale;
545        // add in the background
546        answer += bkg;
547       
548        return answer;
549}
550
551/*      EllipCyl20X  :  calculates the form factor of a elliptical cylinder at the given x-value p->x
552
553Uses 76 pt Gaussian quadrature for orientational integral
554Uses 20 pt quadrature for the inner integral over the elliptical cross-section
555
556Warning:
557The call to WaveData() below returns a pointer to the middle
558of an unlocked Macintosh handle. In the unlikely event that your
559calculations could cause memory to move, you should copy the coefficient
560values to local variables or an array before such operations.
561*/
562double
563EllipCyl20(double dp[], double q)
564{
565        int i,j;
566        double Pi,slde,sld;
567        double scale,ra,nu,length,delrho,bkg;           //local variables of coefficient wave
568        int nordi=76;                   //order of integration
569        int nordj=20;
570        double va,vb;           //upper and lower integration limits
571        double summ,zi,yyy,answer,vell;                 //running tally of integration
572        double summj,vaj,vbj,zij,arg,si;                        //for the inner integration
573
574        Pi = 4.0*atan(1.0);
575        va = 0.0;
576        vb = 1.0;               //orintational average, outer integral
577        vaj=0.0;
578        vbj=Pi;         //endpoints of inner integral
579       
580        summ = 0.0;                     //initialize intergral
581       
582        scale = dp[0];                  //make local copies in case memory moves
583        ra = dp[1];
584        nu = dp[2];
585        length = dp[3];
586        slde = dp[4];
587        sld = dp[5];
588        delrho = slde - sld;
589        bkg = dp[6];
590       
591        for(i=0;i<nordi;i++) {
592                //setup inner integral over the ellipsoidal cross-section
593                summj=0;
594                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy
595                arg = ra*sqrt(1.0-zi*zi);
596                for(j=0;j<nordj;j++) {
597                        //20 gauss points for the inner integral
598                        zij = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "y" dummy
599                        yyy = Gauss20Wt[j] * EllipCylKernel(q,arg,nu,zij);
600                        summj += yyy;
601                }
602                //now calculate the value of the inner integral
603                answer = (vbj-vaj)/2.0*summj;
604                //divide integral by Pi
605                answer /=Pi;
606               
607                //now calculate outer integral
608                arg = q*length*zi/2.0;
609                if (arg == 0.0){
610                        si = 1.0;
611                }else{
612                        si = sin(arg) * sin(arg) / arg / arg;
613                }
614                yyy = Gauss76Wt[i] * answer * si;
615                summ += yyy;
616        }
617       
618        answer = (vb-va)/2.0*summ;
619        // Multiply by contrast^2
620        answer *= delrho*delrho;
621        //normalize by cylinder volume
622        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
623        vell = Pi*ra*(nu*ra)*length;
624        answer *= vell;
625        //convert to [cm-1]
626        answer *= 1.0e8;
627        //Scale
628        answer *= scale;
629        // add in the background
630        answer += bkg;
631
632        return answer;
633}
634
635/*      TriaxialEllipsoidX  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x
636
637Uses 76 pt Gaussian quadrature for both integrals
638
639Warning:
640The call to WaveData() below returns a pointer to the middle
641of an unlocked Macintosh handle. In the unlikely event that your
642calculations could cause memory to move, you should copy the coefficient
643values to local variables or an array before such operations.
644*/
645double
646TriaxialEllipsoid(double dp[], double q)
647{
648        int i,j;
649        double Pi;
650        double scale,aa,bb,cc,delrho,bkg;               //local variables of coefficient wave
651        int nordi=76;                   //order of integration
652        int nordj=76;
653        double va,vb;           //upper and lower integration limits
654        double summ,zi,yyy,answer;                      //running tally of integration
655        double summj,vaj,vbj,zij,slde,sld;                      //for the inner integration
656       
657        Pi = 4.0*atan(1.0);
658        va = 0.0;
659        vb = 1.0;               //orintational average, outer integral
660        vaj = 0.0;
661        vbj = 1.0;              //endpoints of inner integral
662       
663        summ = 0.0;                     //initialize intergral
664       
665        scale = dp[0];                  //make local copies in case memory moves
666        aa = dp[1];
667        bb = dp[2];
668        cc = dp[3];
669        slde = dp[4];
670        sld = dp[5];
671        delrho = slde - sld;
672        bkg = dp[6];
673        for(i=0;i<nordi;i++) {
674                //setup inner integral over the ellipsoidal cross-section
675                summj=0.0;
676                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy
677                for(j=0;j<nordj;j++) {
678                        //20 gauss points for the inner integral
679                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "y" dummy
680                        yyy = Gauss76Wt[j] * TriaxialKernel(q,aa,bb,cc,zi,zij);
681                        summj += yyy;
682                }
683                //now calculate the value of the inner integral
684                answer = (vbj-vaj)/2.0*summj;
685               
686                //now calculate outer integral
687                yyy = Gauss76Wt[i] * answer;
688                summ += yyy;
689        }               //final scaling is done at the end of the function, after the NT_FP64 case
690       
691        answer = (vb-va)/2.0*summ;
692        // Multiply by contrast^2
693        answer *= delrho*delrho;
694        //normalize by ellipsoid volume
695        answer *= 4.0*Pi/3.0*aa*bb*cc;
696        //convert to [cm-1]
697        answer *= 1.0e8;
698        //Scale
699        answer *= scale;
700        // add in the background
701        answer += bkg;
702       
703        return answer;
704}
705
706/*      ParallelepipedX  :  calculates the form factor of a Parallelepiped (a rectangular solid)
707at the given x-value p->x
708
709Uses 76 pt Gaussian quadrature for both integrals
710
711Warning:
712The call to WaveData() below returns a pointer to the middle
713of an unlocked Macintosh handle. In the unlikely event that your
714calculations could cause memory to move, you should copy the coefficient
715values to local variables or an array before such operations.
716*/
717double
718Parallelepiped(double dp[], double q)
719{
720        int i,j;
721        double scale,aa,bb,cc,delrho,bkg;               //local variables of coefficient wave
722        int nordi=76;                   //order of integration
723        int nordj=76;
724        double va,vb;           //upper and lower integration limits
725        double summ,yyy,answer;                 //running tally of integration
726        double summj,vaj,vbj;                   //for the inner integration
727        double mu,mudum,arg,sigma,uu,vol,sldp,sld;
728       
729       
730        //      Pi = 4.0*atan(1.0);
731        va = 0.0;
732        vb = 1.0;               //orintational average, outer integral
733        vaj = 0.0;
734        vbj = 1.0;              //endpoints of inner integral
735
736        summ = 0.0;                     //initialize intergral
737       
738        scale = dp[0];                  //make local copies in case memory moves
739        aa = dp[1];
740        bb = dp[2];
741        cc = dp[3];
742        sldp = dp[4];
743        sld = dp[5];
744        delrho = sldp - sld;
745        bkg = dp[6];
746       
747        mu = q*bb;
748        vol = aa*bb*cc;
749        // normalize all WRT bb
750        aa = aa/bb;
751        cc = cc/bb;
752       
753        for(i=0;i<nordi;i++) {
754                //setup inner integral over the ellipsoidal cross-section
755                summj=0.0;
756                sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;          //the outer dummy
757               
758                for(j=0;j<nordj;j++) {
759                        //76 gauss points for the inner integral
760                        uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;         //the inner dummy
761                        mudum = mu*sqrt(1.0-sigma*sigma);
762                        yyy = Gauss76Wt[j] * PPKernel(aa,mudum,uu);
763                        summj += yyy;
764                }
765                //now calculate the value of the inner integral
766                answer = (vbj-vaj)/2.0*summj;
767
768                arg = mu*cc*sigma/2.0;
769                if ( arg == 0.0 ) {
770                        answer *= 1.0;
771                } else {
772                        answer *= sin(arg)*sin(arg)/arg/arg;
773                }
774               
775                //now sum up the outer integral
776                yyy = Gauss76Wt[i] * answer;
777                summ += yyy;
778        }               //final scaling is done at the end of the function, after the NT_FP64 case
779       
780        answer = (vb-va)/2.0*summ;
781        // Multiply by contrast^2
782        answer *= delrho*delrho;
783        //normalize by volume
784        answer *= vol;
785        //convert to [cm-1]
786        answer *= 1.0e8;
787        //Scale
788        answer *= scale;
789        // add in the background
790        answer += bkg;
791       
792        return answer;
793}
794
795/*      HollowCylinderX  :  calculates the form factor of a Hollow Cylinder
796at the given x-value p->x
797
798Uses 76 pt Gaussian quadrature for the single integral
799
800Warning:
801The call to WaveData() below returns a pointer to the middle
802of an unlocked Macintosh handle. In the unlikely event that your
803calculations could cause memory to move, you should copy the coefficient
804values to local variables or an array before such operations.
805*/
806double
807HollowCylinder(double dp[], double q)
808{
809        int i;
810        double scale,rcore,rshell,length,delrho,bkg;            //local variables of coefficient wave
811        int nord=76;                    //order of integration
812        double va,vb,zi;                //upper and lower integration limits
813        double summ,answer,pi,sldc,sld;                 //running tally of integration
814       
815        pi = 4.0*atan(1.0);
816        va = 0.0;
817        vb = 1.0;               //limits of numerical integral
818
819        summ = 0.0;                     //initialize intergral
820       
821        scale = dp[0];          //make local copies in case memory moves
822        rcore = dp[1];
823        rshell = dp[2];
824        length = dp[3];
825        sldc = dp[4];
826        sld = dp[5];
827        delrho = sldc - sld;
828        bkg = dp[6];
829       
830        for(i=0;i<nord;i++) {
831                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
832                summ += Gauss76Wt[i] * HolCylKernel(q, rcore, rshell, length, zi);
833        }
834       
835        answer = (vb-va)/2.0*summ;
836        // Multiply by contrast^2
837        answer *= delrho*delrho;
838        //normalize by volume
839        answer *= pi*(rshell*rshell-rcore*rcore)*length;
840        //convert to [cm-1]
841        answer *= 1.0e8;
842        //Scale
843        answer *= scale;
844        // add in the background
845        answer += bkg;
846       
847        return answer;
848}
849
850/*      EllipsoidFormX  :  calculates the form factor of an ellipsoid of revolution with semiaxes a:a:nua
851at the given x-value p->x
852
853Uses 76 pt Gaussian quadrature for the single integral
854
855Warning:
856The call to WaveData() below returns a pointer to the middle
857of an unlocked Macintosh handle. In the unlikely event that your
858calculations could cause memory to move, you should copy the coefficient
859values to local variables or an array before such operations.
860*/
861double
862EllipsoidForm(double dp[], double q)
863{
864        int i;
865        double scale,a,nua,delrho,bkg;          //local variables of coefficient wave
866        int nord=76;                    //order of integration
867        double va,vb,zi;                //upper and lower integration limits
868        double summ,answer,pi,slde,sld;                 //running tally of integration
869       
870        pi = 4.0*atan(1.0);
871        va = 0.0;
872        vb = 1.0;               //limits of numerical integral
873
874        summ = 0.0;                     //initialize intergral
875       
876        scale = dp[0];                  //make local copies in case memory moves
877        nua = dp[1];
878        a = dp[2];
879        slde = dp[3];
880        sld = dp[4];
881        delrho = slde - sld;
882        bkg = dp[5];
883       
884        for(i=0;i<nord;i++) {
885                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
886                summ += Gauss76Wt[i] * EllipsoidKernel(q, a, nua, zi);
887        }
888       
889        answer = (vb-va)/2.0*summ;
890        // Multiply by contrast^2
891        answer *= delrho*delrho;
892        //normalize by volume
893        answer *= 4.0*pi/3.0*a*a*nua;
894        //convert to [cm-1]
895        answer *= 1.0e8;
896        //Scale
897        answer *= scale;
898        // add in the background
899        answer += bkg;
900       
901        return answer;
902}
903
904
905/*      Cyl_PolyRadiusX  :  calculates the form factor of a cylinder at the given x-value p->x
906the cylinder has a polydisperse cross section
907
908*/
909double
910Cyl_PolyRadius(double dp[], double q)
911{
912        int i;
913        double scale,radius,length,pd,delrho,bkg;               //local variables of coefficient wave
914        int nord=20;                    //order of integration
915        double uplim,lolim;             //upper and lower integration limits
916        double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
917        double range,zz,Pi,sldc,sld;
918       
919        Pi = 4.0*atan(1.0);
920        range = 3.4;
921       
922        summ = 0.0;                     //initialize intergral
923       
924        scale = dp[0];                  //make local copies in case memory moves
925        radius = dp[1];
926        length = dp[2];
927        pd = dp[3];
928        sldc = dp[4];
929        sld = dp[5];
930        delrho = sldc - sld;
931        bkg = dp[6];
932       
933        zz = (1.0/pd)*(1.0/pd) - 1.0;
934       
935        lolim = radius*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
936        if(lolim<0.0) {
937                lolim = 0.0;
938        }
939        if(pd>0.3) {
940                range = 3.4 + (pd-0.3)*18.0;
941        }
942        uplim = radius*(1.0+range*pd);
943       
944        for(i=0;i<nord;i++) {
945                zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
946                yyy = Gauss20Wt[i] * Cyl_PolyRadKernel(q, radius, length, zz, delrho, zi);
947                summ += yyy;
948        }
949       
950        answer = (uplim-lolim)/2.0*summ;
951        //normalize by average cylinder volume
952        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
953        Vpoly=Pi*radius*radius*length*(zz+2.0)/(zz+1.0);
954        answer /= Vpoly;
955        //convert to [cm-1]
956        answer *= 1.0e8;
957        //Scale
958        answer *= scale;
959        // add in the background
960        answer += bkg;
961       
962        return answer;
963}
964
965/*      Cyl_PolyLengthX  :  calculates the form factor of a cylinder at the given x-value p->x
966the cylinder has a polydisperse Length
967
968*/
969double
970Cyl_PolyLength(double dp[], double q)
971{
972        int i;
973        double scale,radius,length,pd,delrho,bkg;               //local variables of coefficient wave
974        int nord=20;                    //order of integration
975        double uplim,lolim;             //upper and lower integration limits
976        double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
977        double range,zz,Pi,sldc,sld;
978       
979       
980        Pi = 4.0*atan(1.0);
981        range = 3.4;
982       
983        summ = 0.0;                     //initialize intergral
984       
985        scale = dp[0];                  //make local copies in case memory moves
986        radius = dp[1];
987        length = dp[2];
988        pd = dp[3];
989        sldc = dp[4];
990        sld = dp[5];
991        delrho = sldc - sld;
992        bkg = dp[6];
993       
994        zz = (1.0/pd)*(1.0/pd) - 1.0;
995       
996        lolim = length*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
997        if(lolim<0.0) {
998                lolim = 0.0;
999        }
1000        if(pd>0.3) {
1001                range = 3.4 + (pd-0.3)*18.0;
1002        }
1003        uplim = length*(1.0+range*pd);
1004       
1005        for(i=0;i<nord;i++) {
1006                zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1007                yyy = Gauss20Wt[i] * Cyl_PolyLenKernel(q, radius, length, zz, delrho, zi);
1008                summ += yyy;
1009        }
1010       
1011        answer = (uplim-lolim)/2.0*summ;
1012        //normalize by average cylinder volume (first moment)
1013        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
1014        Vpoly=Pi*radius*radius*length;
1015        answer /= Vpoly;
1016        //convert to [cm-1]
1017        answer *= 1.0e8;
1018        //Scale
1019        answer *= scale;
1020        // add in the background
1021        answer += bkg;
1022       
1023        return answer;
1024}
1025
1026/*      CoreShellCylinderX  :  calculates the form factor of a cylinder at the given x-value p->x
1027the cylinder has a core-shell structure
1028
1029*/
1030double
1031CoreShellCylinder(double dp[], double q)
1032{
1033        int i;
1034        double scale,rcore,length,bkg;          //local variables of coefficient wave
1035        double thick,rhoc,rhos,rhosolv;
1036        int nord=76;                    //order of integration
1037        double uplim,lolim,halfheight;          //upper and lower integration limits
1038        double summ,zi,yyy,answer,Vcyl;                 //running tally of integration
1039        double Pi;
1040       
1041        Pi = 4.0*atan(1.0);
1042       
1043        lolim = 0.0;
1044        uplim = Pi/2.0;
1045       
1046        summ = 0.0;                     //initialize intergral
1047       
1048        scale = dp[0];          //make local copies in case memory moves
1049        rcore = dp[1];
1050        thick = dp[2];
1051        length = dp[3];
1052        rhoc = dp[4];
1053        rhos = dp[5];
1054        rhosolv = dp[6];
1055        bkg = dp[7];
1056       
1057        halfheight = length/2.0;
1058       
1059        for(i=0;i<nord;i++) {
1060                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1061                yyy = Gauss76Wt[i] * CoreShellCylKernel(q, rcore, thick, rhoc,rhos,rhosolv, halfheight, zi);
1062                summ += yyy;
1063        }
1064       
1065        answer = (uplim-lolim)/2.0*summ;
1066        // length is the total core length
1067        Vcyl=Pi*(rcore+thick)*(rcore+thick)*(length+2.0*thick);
1068        answer /= Vcyl;
1069        //convert to [cm-1]
1070        answer *= 1.0e8;
1071        //Scale
1072        answer *= scale;
1073        // add in the background
1074        answer += bkg;
1075       
1076        return answer;
1077}
1078
1079
1080/*      PolyCoShCylinderX  :  calculates the form factor of a core-shell cylinder at the given x-value p->x
1081the cylinder has a polydisperse CORE radius
1082
1083*/
1084double
1085PolyCoShCylinder(double dp[], double q)
1086{
1087        int i;
1088        double scale,radius,length,sigma,bkg;           //local variables of coefficient wave
1089        double rad,radthick,facthick,rhoc,rhos,rhosolv;
1090        int nord=20;                    //order of integration
1091        double uplim,lolim;             //upper and lower integration limits
1092        double summ,yyy,answer,Vpoly;                   //running tally of integration
1093        double Pi,AR,Rsqrsumm,Rsqryyy,Rsqr;
1094               
1095        Pi = 4.0*atan(1.0);
1096       
1097        summ = 0.0;                     //initialize intergral
1098        Rsqrsumm = 0.0;
1099       
1100        scale = dp[0];
1101        radius = dp[1];
1102        sigma = dp[2];                          //sigma is the standard mean deviation
1103        length = dp[3];
1104        radthick = dp[4];
1105        facthick= dp[5];
1106        rhoc = dp[6];
1107        rhos = dp[7];
1108        rhosolv = dp[8];
1109        bkg = dp[9];
1110       
1111        lolim = exp(log(radius)-(4.*sigma));
1112        if (lolim<0.0) {
1113                lolim=0.0;              //to avoid numerical error when  va<0 (-ve r value)
1114        }
1115        uplim = exp(log(radius)+(4.*sigma));
1116       
1117        for(i=0;i<nord;i++) {
1118                rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1119                AR=(1.0/(rad*sigma*sqrt(2.0*Pi)))*exp(-(0.5*((log(radius/rad))/sigma)*((log(radius/rad))/sigma)));
1120                yyy = AR* Gauss20Wt[i] * CSCylIntegration(q,rad,radthick,facthick,rhoc,rhos,rhosolv,length);
1121                Rsqryyy= Gauss20Wt[i] * AR * (rad+radthick)*(rad+radthick);             //SRK normalize to total dimensions
1122                summ += yyy;
1123                Rsqrsumm += Rsqryyy;
1124        }
1125       
1126        answer = (uplim-lolim)/2.0*summ;
1127        Rsqr = (uplim-lolim)/2.0*Rsqrsumm;
1128        //normalize by average cylinder volume
1129        Vpoly = Pi*Rsqr*(length+2*facthick);
1130        answer /= Vpoly;
1131        //convert to [cm-1]
1132        answer *= 1.0e8;
1133        //Scale
1134        answer *= scale;
1135        // add in the background
1136        answer += bkg;
1137       
1138        return answer;
1139}
1140
1141/*      OblateFormX  :  calculates the form factor of a core-shell Oblate ellipsoid at the given x-value p->x
1142the ellipsoid has a core-shell structure
1143
1144*/
1145double
1146OblateForm(double dp[], double q)
1147{
1148        int i;
1149        double scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg;
1150        int nord=76;                    //order of integration
1151        double uplim,lolim;             //upper and lower integration limits
1152        double summ,zi,yyy,answer,oblatevol;                    //running tally of integration
1153        double Pi,sldc,slds,sld;
1154       
1155        Pi = 4.0*atan(1.0);
1156       
1157        lolim = 0.0;
1158        uplim = 1.0;
1159       
1160        summ = 0.0;                     //initialize intergral
1161       
1162       
1163        scale = dp[0];          //make local copies in case memory moves
1164        crmaj = dp[1];
1165        crmin = dp[2];
1166        trmaj = dp[3];
1167        trmin = dp[4];
1168        sldc = dp[5];
1169        slds = dp[6];
1170        sld = dp[7];
1171        delpc = sldc - slds;    //core - shell
1172        delps = slds - sld;             //shell - solvent
1173        bkg = dp[8];
1174
1175        for(i=0;i<nord;i++) {
1176                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1177                yyy = Gauss76Wt[i] * gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,q);
1178                summ += yyy;
1179        }
1180       
1181        answer = (uplim-lolim)/2.0*summ;
1182        // normalize by particle volume
1183        oblatevol = 4.0*Pi/3.0*trmaj*trmaj*trmin;
1184        answer /= oblatevol;
1185       
1186        //convert to [cm-1]
1187        answer *= 1.0e8;
1188        //Scale
1189        answer *= scale;
1190        // add in the background
1191        answer += bkg;
1192       
1193        return answer;
1194}
1195
1196/*      ProlateFormX  :  calculates the form factor of a core-shell Prolate ellipsoid at the given x-value p->x
1197the ellipsoid has a core-shell structure
1198
1199*/
1200double
1201ProlateForm(double dp[], double q)
1202{
1203        int i;
1204        double scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg;
1205        int nord=76;                    //order of integration
1206        double uplim,lolim;             //upper and lower integration limits
1207        double summ,zi,yyy,answer,prolatevol;                   //running tally of integration
1208        double Pi,sldc,slds,sld;
1209       
1210        Pi = 4.0*atan(1.0);
1211       
1212        lolim = 0.0;
1213        uplim = 1.0;
1214       
1215        summ = 0.0;                     //initialize intergral
1216       
1217        scale = dp[0];          //make local copies in case memory moves
1218        crmaj = dp[1];
1219        crmin = dp[2];
1220        trmaj = dp[3];
1221        trmin = dp[4];
1222        sldc = dp[5];
1223        slds = dp[6];
1224        sld = dp[7];
1225        delpc = sldc - slds;            //core - shell
1226        delps = slds - sld;             //shell  - sovent
1227        bkg = dp[8];
1228
1229        for(i=0;i<nord;i++) {
1230                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1231                yyy = Gauss76Wt[i] * gfn2(zi,crmaj,crmin,trmaj,trmin,delpc,delps,q);
1232                summ += yyy;
1233        }
1234       
1235        answer = (uplim-lolim)/2.0*summ;
1236        // normalize by particle volume
1237        prolatevol = 4.0*Pi/3.0*trmaj*trmin*trmin;
1238        answer /= prolatevol;
1239       
1240        //convert to [cm-1]
1241        answer *= 1.0e8;
1242        //Scale
1243        answer *= scale;
1244        // add in the background
1245        answer += bkg;
1246       
1247        return answer;
1248}
1249
1250
1251/*      StackedDiscsX  :  calculates the form factor of a stacked "tactoid" of core shell disks
1252like clay platelets that are not exfoliated
1253
1254*/
1255double
1256StackedDiscs(double dp[], double q)
1257{
1258        int i;
1259        double scale,length,bkg,rcore,thick,rhoc,rhol,rhosolv,N,gsd;            //local variables of coefficient wave
1260        double va,vb,vcyl,summ,yyy,zi,halfheight,d,answer;
1261        int nord=76;                    //order of integration
1262        double Pi;
1263       
1264       
1265        Pi = 4.0*atan(1.0);
1266       
1267        va = 0.0;
1268        vb = Pi/2.0;
1269       
1270        summ = 0.0;                     //initialize intergral
1271       
1272        scale = dp[0];
1273        rcore = dp[1];
1274        length = dp[2];
1275        thick = dp[3];
1276        rhoc = dp[4];
1277        rhol = dp[5];
1278        rhosolv = dp[6];
1279        N = dp[7];
1280        gsd = dp[8];
1281        bkg = dp[9];
1282       
1283        d=2.0*thick+length;
1284        halfheight = length/2.0;
1285       
1286        for(i=0;i<nord;i++) {
1287                zi = ( Gauss76Z[i]*(vb-va) + vb + va )/2.0;
1288                yyy = Gauss76Wt[i] * Stackdisc_kern(q, rcore, rhoc,rhol,rhosolv, halfheight,thick,zi,gsd,d,N);
1289                summ += yyy;
1290        }
1291       
1292        answer = (vb-va)/2.0*summ;
1293        // length is the total core length
1294        vcyl=Pi*rcore*rcore*(2.0*thick+length)*N;
1295        answer /= vcyl;
1296        //Convert to [cm-1]
1297        answer *= 1.0e8;
1298        //Scale
1299        answer *= scale;
1300        // add in the background
1301        answer += bkg;
1302       
1303        return answer;
1304}
1305
1306
1307/*      LamellarFFX  :  calculates the form factor of a lamellar structure - no S(q) effects included
1308
1309*/
1310double
1311LamellarFF(double dp[], double q)
1312{
1313        double scale,del,sig,contr,bkg;         //local variables of coefficient wave
1314        double inten, qval,Pq;
1315        double Pi,sldb,sld;
1316       
1317       
1318        Pi = 4.0*atan(1.0);
1319        scale = dp[0];
1320        del = dp[1];
1321        sig = dp[2]*del;
1322        sldb = dp[3];
1323        sld = dp[4];
1324        contr = sldb - sld;
1325        bkg = dp[5];
1326        qval=q;
1327       
1328        Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig));
1329       
1330        inten = 2.0*Pi*scale*Pq/(qval*qval);            //this is now dimensionless...
1331       
1332        inten /= del;                   //normalize by the thickness (in A)
1333       
1334        inten *= 1.0e8;         // 1/A to 1/cm
1335       
1336        return(inten+bkg);
1337}
1338
1339/*      LamellarPSX  :  calculates the form factor of a lamellar structure - with S(q) effects included
1340--- now the proper resolution effects are used - the "default" resolution is turned off (= 0) and the
1341model is smeared just like any other function
1342*/
1343double
1344LamellarPS(double dp[], double q)
1345{
1346        double scale,dd,del,sig,contr,NN,Cp,bkg;                //local variables of coefficient wave
1347        double inten,qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ;
1348        double Pi,Euler,dQDefault,fii,sldb,sld;
1349        int ii,NNint;
1350//      char buf[256];
1351
1352       
1353        Euler = 0.5772156649;           // Euler's constant
1354//      dQDefault = 0.0025;             //[=] 1/A, q-resolution, default value
1355        dQDefault = 0.0;
1356        dQ = dQDefault;
1357       
1358        Pi = 4.0*atan(1.0);
1359        qval = q;
1360       
1361        scale = dp[0];
1362        dd = dp[1];
1363        del = dp[2];
1364        sig = dp[3]*del;
1365        sldb = dp[4];
1366        sld = dp[5];
1367        contr = sldb - sld;
1368        NN = trunc(dp[6]);              //be sure that NN is an integer
1369        Cp = dp[7];
1370        bkg = dp[8];
1371       
1372        Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig));
1373       
1374        NNint = (int)NN;                //cast to an integer for the loop
1375       
1376//      sprintf(buf, "qval = %g\r", qval);
1377//      XOPNotice(buf);
1378       
1379        ii=0;
1380        Sq = 0.0;
1381        for(ii=1;ii<=(NNint-1);ii+=1) {
1382       
1383                fii = (double)ii;               //do I really need to do this?
1384               
1385                temp = 0.0;
1386                alpha = Cp/4.0/Pi/Pi*(log(Pi*fii) + Euler);
1387                t1 = 2.0*dQ*dQ*dd*dd*alpha;
1388                t2 = 2.0*qval*qval*dd*dd*alpha;
1389                t3 = dQ*dQ*dd*dd*fii*fii;
1390               
1391                temp = 1.0-fii/NN;
1392                temp *= cos(dd*qval*fii/(1.0+t1));
1393                temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
1394                temp /= sqrt(1.0+t1);
1395               
1396                Sq += temp;
1397        }
1398       
1399        Sq *= 2.0;
1400        Sq += 1.0;
1401       
1402        inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval);
1403       
1404        inten *= 1.0e8;         // 1/A to 1/cm
1405       
1406    return(inten+bkg);
1407}
1408
1409
1410/*      LamellarPS_HGX  :  calculates the form factor of a lamellar structure - with S(q) effects included
1411--- now the proper resolution effects are used - the "default" resolution is turned off (= 0) and the
1412model is smeared just like any other function
1413*/
1414double
1415LamellarPS_HG(double dp[], double q)
1416{
1417        double scale,dd,delT,delH,SLD_T,SLD_H,SLD_S,NN,Cp,bkg;          //local variables of coefficient wave
1418        double inten,qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ,drh,drt;
1419        double Pi,Euler,dQDefault,fii;
1420        int ii,NNint;
1421       
1422       
1423        Euler = 0.5772156649;           // Euler's constant
1424//      dQDefault = 0.0025;             //[=] 1/A, q-resolution, default value
1425        dQDefault = 0.0;
1426        dQ = dQDefault;
1427       
1428        Pi = 4.0*atan(1.0);
1429        qval= q;
1430       
1431        scale = dp[0];
1432        dd = dp[1];
1433        delT = dp[2];
1434        delH = dp[3];
1435        SLD_T = dp[4];
1436        SLD_H = dp[5];
1437        SLD_S = dp[6];
1438        NN = trunc(dp[7]);              //be sure that NN is an integer
1439        Cp = dp[8];
1440        bkg = dp[9];
1441       
1442       
1443        drh = SLD_H - SLD_S;
1444        drt = SLD_T - SLD_S;    //correction 13FEB06 by L.Porcar
1445       
1446        Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT);
1447        Pq *= Pq;
1448        Pq *= 4.0/(qval*qval);
1449       
1450        NNint = (int)NN;                //cast to an integer for the loop
1451        ii=0;
1452        Sq = 0.0;
1453        for(ii=1;ii<=(NNint-1);ii+=1) {
1454       
1455                fii = (double)ii;               //do I really need to do this?
1456               
1457                temp = 0.0;
1458                alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler);
1459                t1 = 2.0*dQ*dQ*dd*dd*alpha;
1460                t2 = 2.0*qval*qval*dd*dd*alpha;
1461                t3 = dQ*dQ*dd*dd*ii*ii;
1462               
1463                temp = 1.0-ii/NN;
1464                temp *= cos(dd*qval*ii/(1.0+t1));
1465                temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
1466                temp /= sqrt(1.0+t1);
1467               
1468                Sq += temp;
1469        }
1470       
1471        Sq *= 2.0;
1472        Sq += 1.0;
1473       
1474        inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval);
1475       
1476        inten *= 1.0e8;         // 1/A to 1/cm
1477       
1478        return(inten+bkg);
1479}
1480
1481/*      LamellarFF_HGX  :  calculates the form factor of a lamellar structure - no S(q) effects included
1482but extra SLD for head groups is included
1483
1484*/
1485double
1486LamellarFF_HG(double dp[], double q)
1487{
1488        double scale,delT,delH,slds,sldh,sldt,bkg;              //local variables of coefficient wave
1489        double inten, qval,Pq,drh,drt;
1490        double Pi;
1491       
1492       
1493        Pi = 4.0*atan(1.0);
1494        qval= q;
1495        scale = dp[0];
1496        delT = dp[1];
1497        delH = dp[2];
1498        sldt = dp[3];
1499        sldh = dp[4];
1500        slds = dp[5];
1501        bkg = dp[6];
1502       
1503       
1504        drh = sldh - slds;
1505        drt = sldt - slds;              //correction 13FEB06 by L.Porcar
1506       
1507        Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT);
1508        Pq *= Pq;
1509        Pq *= 4.0/(qval*qval);
1510       
1511        inten = 2.0*Pi*scale*Pq/(qval*qval);            //dimensionless...
1512       
1513        inten /= 2.0*(delT+delH);                       //normalize by the bilayer thickness
1514       
1515        inten *= 1.0e8;         // 1/A to 1/cm
1516       
1517        return(inten+bkg);
1518}
1519
1520/*      FlexExclVolCylX  :  calculates the form factor of a flexible cylinder with a circular cross section
1521-- incorporates Wei-Ren Chen's fixes - 2006
1522
1523        */
1524double
1525FlexExclVolCyl(double dp[], double q)
1526{
1527        double scale,L,B,bkg,rad,qr,cont,sldc,slds;
1528        double Pi,flex,crossSect,answer;
1529       
1530       
1531        Pi = 4.0*atan(1.0);
1532       
1533        scale = dp[0];          //make local copies in case memory moves
1534        L = dp[1];
1535        B = dp[2];
1536        rad = dp[3];
1537        sldc = dp[4];
1538        slds = dp[5];
1539        cont = sldc-slds;
1540        bkg = dp[6];
1541       
1542   
1543        qr = q*rad;
1544       
1545        flex = Sk_WR(q,L,B);
1546       
1547        crossSect = (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
1548        flex *= crossSect;
1549        flex *= Pi*rad*rad*L;
1550        flex *= cont*cont;
1551        flex *= 1.0e8;
1552        answer = scale*flex + bkg;
1553       
1554        return answer;
1555}
1556
1557/*      FlexCyl_EllipX  :  calculates the form factor of a flexible cylinder with an elliptical cross section
1558-- incorporates Wei-Ren Chen's fixes - 2006
1559
1560        */
1561double
1562FlexCyl_Ellip(double dp[], double q)
1563{
1564        double scale,L,B,bkg,rad,qr,cont,ellRatio,slds,sldc;
1565        double Pi,flex,crossSect,answer;
1566       
1567       
1568        Pi = 4.0*atan(1.0);
1569        scale = dp[0];          //make local copies in case memory moves
1570        L = dp[1];
1571        B = dp[2];
1572        rad = dp[3];
1573        ellRatio = dp[4];
1574        sldc = dp[5];
1575        slds = dp[6];
1576        bkg = dp[7];
1577       
1578        cont = sldc - slds;
1579        qr = q*rad;
1580       
1581        flex = Sk_WR(q,L,B);
1582       
1583        crossSect = EllipticalCross_fn(q,rad,(rad*ellRatio));
1584        flex *= crossSect;
1585        flex *= Pi*rad*rad*ellRatio*L;
1586        flex *= cont*cont;
1587        flex *= 1.0e8;
1588        answer = scale*flex + bkg;
1589       
1590        return answer;
1591}
1592
1593double
1594EllipticalCross_fn(double qq, double a, double b)
1595{
1596    double uplim,lolim,Pi,summ,arg,zi,yyy,answer;
1597    int i,nord=76;
1598   
1599    Pi = 4.0*atan(1.0);
1600    lolim=0.0;
1601    uplim=Pi/2.0;
1602    summ=0.0;
1603   
1604    for(i=0;i<nord;i++) {
1605                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1606                arg = qq*sqrt(a*a*sin(zi)*sin(zi)+b*b*cos(zi)*cos(zi));
1607                yyy = pow((2.0 * NR_BessJ1(arg) / arg),2);
1608                yyy *= Gauss76Wt[i];
1609                summ += yyy;
1610    }
1611    answer = (uplim-lolim)/2.0*summ;
1612    answer *= 2.0/Pi;
1613    return(answer);
1614       
1615}
1616/*      FlexCyl_PolyLenX  :  calculates the form factor of a flecible cylinder at the given x-value p->x
1617the cylinder has a polydisperse Length
1618
1619*/
1620double
1621FlexCyl_PolyLen(double dp[], double q)
1622{
1623        int i;
1624        double scale,radius,length,pd,bkg,lb,delrho,sldc,slds;          //local variables of coefficient wave
1625        int nord=20;                    //order of integration
1626        double uplim,lolim;             //upper and lower integration limits
1627        double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
1628        double range,zz,Pi;
1629       
1630        Pi = 4.0*atan(1.0);
1631        range = 3.4;
1632       
1633        summ = 0.0;                     //initialize intergral
1634        scale = dp[0];                  //make local copies in case memory moves
1635        length = dp[1];                 //radius
1636        pd = dp[2];                     // average length
1637        lb = dp[3];
1638        radius = dp[4];
1639        sldc = dp[5];
1640        slds = dp[6];
1641        bkg = dp[7];
1642       
1643        delrho = sldc - slds;
1644        zz = (1.0/pd)*(1.0/pd) - 1.0;
1645       
1646        lolim = length*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
1647        if(lolim<0.0) {
1648                lolim = 0.0;
1649        }
1650        if(pd>0.3) {
1651                range = 3.4 + (pd-0.3)*18.0;
1652        }
1653        uplim = length*(1.0+range*pd);
1654       
1655        for(i=0;i<nord;i++) {
1656                zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1657                yyy = Gauss20Wt[i] * FlePolyLen_kernel(q,radius,length,lb,zz,delrho,zi);
1658                summ += yyy;
1659        }
1660       
1661        answer = (uplim-lolim)/2.0*summ;
1662        //normalize by average cylinder volume (first moment), using the average length
1663        Vpoly=Pi*radius*radius*length;
1664        answer /= Vpoly;
1665       
1666        answer *=delrho*delrho;
1667       
1668        //convert to [cm-1]
1669        answer *= 1.0e8;
1670        //Scale
1671        answer *= scale;
1672        // add in the background
1673        answer += bkg;
1674       
1675        return answer;
1676}
1677
1678/*      FlexCyl_PolyLenX  :  calculates the form factor of a flexible cylinder at the given x-value p->x
1679the cylinder has a polydisperse cross sectional radius
1680
1681*/
1682double
1683FlexCyl_PolyRad(double dp[], double q)
1684{
1685        int i;
1686        double scale,radius,length,pd,delrho,bkg,lb,sldc,slds;          //local variables of coefficient wave
1687        int nord=76;                    //order of integration
1688        double uplim,lolim;             //upper and lower integration limits
1689        double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
1690        double range,zz,Pi;
1691       
1692       
1693        Pi = 4.0*atan(1.0);
1694        range = 3.4;
1695       
1696        summ = 0.0;                     //initialize intergral
1697       
1698        scale = dp[0];                  //make local copies in case memory moves
1699        length = dp[1];                 //radius
1700        lb = dp[2];                     // average length
1701        radius = dp[3];
1702        pd = dp[4];
1703        sldc = dp[5];
1704        slds = dp[6];
1705        bkg = dp[7];
1706       
1707        delrho = sldc-slds;
1708        zz = (1.0/pd)*(1.0/pd) - 1.0;
1709       
1710        lolim = radius*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
1711        if(lolim<0.0) {
1712                lolim = 0.0;
1713        }
1714        if(pd>0.3) {
1715                range = 3.4 + (pd-0.3)*18.0;
1716        }
1717        uplim = radius*(1.0+range*pd);
1718       
1719        for(i=0;i<nord;i++) {
1720                //zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1721                //yyy = Gauss20Wt[i] * FlePolyRad_kernel(q,radius,length,lb,zz,delrho,zi);
1722                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1723                yyy = Gauss76Wt[i] * FlePolyRad_kernel(q,radius,length,lb,zz,delrho,zi);
1724                summ += yyy;
1725        }
1726       
1727        answer = (uplim-lolim)/2.0*summ;
1728        //normalize by average cylinder volume (second moment), using the average radius
1729        Vpoly = Pi*radius*radius*length*(zz+2.0)/(zz+1.0);
1730        answer /= Vpoly;
1731       
1732        answer *=delrho*delrho;
1733       
1734        //convert to [cm-1]
1735        answer *= 1.0e8;
1736        //Scale
1737        answer *= scale;
1738        // add in the background
1739        answer += bkg;
1740       
1741        return answer;
1742}
1743
1744
1745
1746///////////////
1747
1748//
1749//     FUNCTION gfn2:    CONTAINS F(Q,A,B,mu)**2  AS GIVEN
1750//                       BY (53) AND (56,57) IN CHEN AND
1751//                       KOTLARCHYK REFERENCE
1752//
1753//     <PROLATE ELLIPSOIDS>
1754//
1755double
1756gfn2(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq)
1757{
1758        // local variables
1759        double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,gfn2,pi43,gfn,Pi;
1760
1761        Pi = 4.0*atan(1.0);
1762       
1763        pi43=4.0/3.0*Pi;
1764        aa = crmaj;
1765        bb = crmin;
1766        u2 = (aa*aa*xx*xx + bb*bb*(1.0-xx*xx));
1767        ut2 = (trmaj*trmaj*xx*xx + trmin*trmin*(1.0-xx*xx));
1768        uq = sqrt(u2)*qq;
1769        ut= sqrt(ut2)*qq;
1770        vc = pi43*aa*bb*bb;
1771        vt = pi43*trmaj*trmin*trmin;
1772        if (uq == 0.0){
1773                siq = 1.0/3.0;
1774        }else{
1775                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq;
1776        }
1777        if (ut == 0.0){
1778                sit = 1.0/3.0;
1779        }else{
1780                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut;
1781        }
1782        gfnc = 3.0*siq*vc*delpc;
1783        gfnt = 3.0*sit*vt*delps;
1784        gfn = gfnc+gfnt;
1785        gfn2 = gfn*gfn;
1786       
1787        return (gfn2);
1788}
1789
1790//
1791//     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN
1792//                       BY (53) & (58-59) IN CHEN AND
1793//                       KOTLARCHYK REFERENCE
1794//
1795//       <OBLATE ELLIPSOID>
1796// function gfn4 for oblate ellipsoids
1797double
1798gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq)
1799{
1800        // local variables
1801        double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,tgfn,gfn4,pi43,Pi;
1802
1803        Pi = 4.0*atan(1.0);
1804        pi43=4.0/3.0*Pi;
1805        aa = crmaj;
1806        bb = crmin;
1807        u2 = (bb*bb*xx*xx + aa*aa*(1.0-xx*xx));
1808        ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx));
1809        uq = sqrt(u2)*qq;
1810        ut= sqrt(ut2)*qq;
1811        vc = pi43*aa*aa*bb;
1812        vt = pi43*trmaj*trmaj*trmin;
1813        if (uq == 0.0){
1814                siq = 1.0/3.0;
1815        }else{
1816                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq;
1817        }
1818        if (ut == 0.0){
1819                sit = 1.0/3.0;
1820        }else{
1821                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut;
1822        }
1823        gfnc = 3.0*siq*vc*delpc;
1824        gfnt = 3.0*sit*vt*delps;
1825        tgfn = gfnc+gfnt;
1826        gfn4 = tgfn*tgfn;
1827       
1828        return (gfn4);
1829}
1830
1831double
1832FlePolyLen_kernel(double q, double radius, double length, double lb, double zz, double delrho, double zi)
1833{
1834    double Pq,vcyl,dl;
1835    double Pi,qr;
1836   
1837    Pi = 4.0*atan(1.0);
1838    qr = q*radius;
1839   
1840    Pq = Sk_WR(q,zi,lb);                //does not have cross section term
1841    if (qr !=0){
1842        Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
1843    } 
1844    vcyl=Pi*radius*radius*zi;
1845    Pq *= vcyl*vcyl;
1846   
1847    dl = SchulzPoint_cpr(zi,length,zz);
1848    return (Pq*dl);     
1849       
1850}
1851
1852double
1853FlePolyRad_kernel(double q, double ravg, double Lc, double Lb, double zz, double delrho, double zi)
1854{
1855    double Pq,vcyl,dr;
1856    double Pi,qr;
1857   
1858    Pi = 4.0*atan(1.0);
1859    qr = q*zi;
1860   
1861    Pq = Sk_WR(q,Lc,Lb);                //does not have cross section term
1862    if (qr !=0){
1863        Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
1864    }
1865
1866    vcyl=Pi*zi*zi*Lc;
1867    Pq *= vcyl*vcyl;
1868   
1869    dr = SchulzPoint_cpr(zi,ravg,zz);
1870    return (Pq*dr);     
1871       
1872}
1873
1874double
1875CSCylIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length)
1876{
1877        double answer,halfheight,Pi;
1878        double lolim,uplim,summ,yyy,zi;
1879        int nord,i;
1880       
1881        // set up the integration end points
1882        Pi = 4.0*atan(1.0);
1883        nord = 76;
1884        lolim = 0.0;
1885        uplim = Pi/2.0;
1886        halfheight = length/2.0;
1887       
1888        summ = 0.0;                             // initialize integral
1889        i=0;
1890        for(i=0;i<nord;i++) {
1891                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1892                yyy = Gauss76Wt[i] * CScyl(qq, rad, radthick, facthick, rhoc,rhos,rhosolv, halfheight, zi);
1893                summ += yyy;
1894        }
1895       
1896        // calculate value of integral to return
1897        answer = (uplim-lolim)/2.0*summ;
1898        return (answer);
1899}
1900
1901double
1902CScyl(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length, double dum)
1903{       
1904        // qq is the q-value for the calculation (1/A)
1905        // radius is the core radius of the cylinder (A)
1906        //  radthick and facthick are the radial and face layer thicknesses
1907        // rho(n) are the respective SLD's
1908        // length is the *Half* CORE-LENGTH of the cylinder
1909        // dum is the dummy variable for the integration (theta)
1910
1911        double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval;
1912        double Pi;
1913       
1914        Pi = 4.0*atan(1.0); 
1915       
1916        dr1 = rhoc-rhos;
1917        dr2 = rhos-rhosolv;
1918        vol1 = Pi*rad*rad*(2.0*length);
1919        vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick);
1920       
1921        besarg1 = qq*rad*sin(dum);
1922        besarg2 = qq*(rad+radthick)*sin(dum);
1923        sinarg1 = qq*length*cos(dum);
1924        sinarg2 = qq*(length+facthick)*cos(dum);
1925        if (besarg1 == 0.0){
1926                be1 = 0.5;
1927        }else{
1928                be1 = NR_BessJ1(besarg1)/besarg1;
1929        }
1930        if (besarg2 == 0.0){
1931                be2 = 0.5;
1932        }else{
1933                be2 = NR_BessJ1(besarg2)/besarg2;
1934        }
1935        if (sinarg1 == 0.0){
1936                si1 = 1.0;
1937        }else{
1938                si1 = sin(sinarg1)/sinarg1;
1939        }
1940        if (sinarg2 == 0.0){
1941                si2 = 1.0;
1942        }else{
1943                si2 = sin(sinarg2)/sinarg2;
1944        }
1945
1946        t1 = 2.0*vol1*dr1*si1*be1;
1947        t2 = 2.0*vol2*dr2*si2*be2;
1948
1949        retval = ((t1+t2)*(t1+t2))*sin(dum);
1950        return (retval);
1951   
1952}
1953
1954
1955double
1956CoreShellCylKernel(double qq, double rcore, double thick, double rhoc, double rhos, double rhosolv, double length, double dum)
1957{
1958
1959    double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval;
1960    double Pi;
1961   
1962    Pi = 4.0*atan(1.0);
1963   
1964    dr1 = rhoc-rhos;
1965    dr2 = rhos-rhosolv;
1966    vol1 = Pi*rcore*rcore*(2.0*length);
1967    vol2 = Pi*(rcore+thick)*(rcore+thick)*(2.0*length+2.0*thick);
1968   
1969    besarg1 = qq*rcore*sin(dum);
1970    besarg2 = qq*(rcore+thick)*sin(dum);
1971    sinarg1 = qq*length*cos(dum);
1972    sinarg2 = qq*(length+thick)*cos(dum);
1973
1974        if (besarg1 == 0.0){
1975                be1 = 0.5;
1976        }else{
1977                be1 = NR_BessJ1(besarg1)/besarg1;
1978        }
1979        if (besarg2 == 0.0){
1980                be2 = 0.5;
1981        }else{
1982                be2 = NR_BessJ1(besarg2)/besarg2;
1983        }
1984        if (sinarg1 == 0.0){
1985                si1 = 1.0;
1986        }else{
1987                si1 = sin(sinarg1)/sinarg1;
1988        }
1989        if (sinarg2 == 0.0){
1990                si2 = 1.0;
1991        }else{
1992                si2 = sin(sinarg2)/sinarg2;
1993        }
1994
1995    t1 = 2.0*vol1*dr1*si1*be1;
1996    t2 = 2.0*vol2*dr2*si2*be2;
1997
1998    retval = ((t1+t2)*(t1+t2))*sin(dum);
1999       
2000    return (retval);
2001}
2002
2003double
2004Cyl_PolyLenKernel(double q, double radius, double len_avg, double zz, double delrho, double dumLen)
2005{
2006       
2007    double halfheight,uplim,lolim,zi,summ,yyy,Pi;
2008    double answer,dr,Vcyl;
2009    int i,nord;
2010   
2011    Pi = 4.0*atan(1.0);
2012    lolim = 0.0;
2013    uplim = Pi/2.0;
2014    halfheight = dumLen/2.0;
2015    nord=20;
2016    summ = 0.0;
2017   
2018    //do the cylinder orientational average
2019    for(i=0;i<nord;i++) {
2020                zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
2021                yyy = Gauss20Wt[i] * CylKernel(q, radius, halfheight, zi);
2022                summ += yyy;
2023    }
2024    answer = (uplim-lolim)/2.0*summ;
2025    // Multiply by contrast^2
2026    answer *= delrho*delrho;
2027    // don't do the normal scaling to volume here
2028    // instead, multiply by VCyl^2 to get rid of the normalization for this radius of cylinder
2029    Vcyl = Pi*radius*radius*dumLen;
2030    answer *= Vcyl*Vcyl;
2031   
2032    dr = SchulzPoint_cpr(dumLen,len_avg,zz);
2033    return(dr*answer);
2034}
2035
2036
2037double
2038Stackdisc_kern(double qq, double rcore, double rhoc, double rhol, double rhosolv, double length, double thick, double dum, double gsd, double d, double N)
2039{               
2040        // qq is the q-value for the calculation (1/A)
2041        // rcore is the core radius of the cylinder (A)
2042        // rho(n) are the respective SLD's
2043        // length is the *Half* CORE-LENGTH of the cylinder = L (A)
2044        // dum is the dummy variable for the integration (x in Feigin's notation)
2045
2046        //Local variables
2047        double totald,dr1,dr2,besarg1,besarg2,be1,be2,si1,si2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt;
2048        double Pi;
2049        int kk;
2050       
2051        Pi = 4.0*atan(1.0);
2052       
2053        dr1 = rhoc-rhosolv;
2054        dr2 = rhol-rhosolv;
2055        area = Pi*rcore*rcore;
2056        totald=2.0*(thick+length);
2057       
2058        besarg1 = qq*rcore*sin(dum);
2059        besarg2 = qq*rcore*sin(dum);
2060       
2061        sinarg1 = qq*length*cos(dum);
2062        sinarg2 = qq*(length+thick)*cos(dum);
2063
2064        if (besarg1 == 0.0){
2065                be1 = 0.5;
2066        }else{
2067                be1 = NR_BessJ1(besarg1)/besarg1;
2068        }
2069        if (besarg2 == 0.0){
2070                be2 = 0.5;
2071        }else{
2072                be2 = NR_BessJ1(besarg2)/besarg2;
2073        }
2074        if (sinarg1 == 0.0){
2075                si1 = 1.0;
2076        }else{
2077                si1 = sin(sinarg1)/sinarg1;
2078        }
2079        if (sinarg2 == 0.0){
2080                si2 = 1.0;
2081        }else{
2082                si2 = sin(sinarg2)/sinarg2;
2083        }
2084
2085        t1 = 2.0*area*(2.0*length)*dr1*(si1)*(be1);
2086        t2 = 2.0*area*dr2*(totald*si2-2.0*length*si1)*(be2);
2087
2088        retval =((t1+t2)*(t1+t2))*sin(dum);
2089       
2090        // loop for the structure facture S(q)
2091        sqq=0.0;
2092        for(kk=1;kk<N;kk+=1) {
2093                dexpt=qq*cos(dum)*qq*cos(dum)*d*d*gsd*gsd*kk/2.0;
2094                sqq=sqq+(N-kk)*cos(qq*cos(dum)*d*kk)*exp(-1.*dexpt);
2095        }                       
2096       
2097        // end of loop for S(q)
2098        sqq=1.0+2.0*sqq/N;
2099        retval *= sqq;
2100   
2101        return(retval);
2102}
2103
2104
2105double
2106Cyl_PolyRadKernel(double q, double radius, double length, double zz, double delrho, double dumRad)
2107{
2108       
2109    double halfheight,uplim,lolim,zi,summ,yyy,Pi;
2110    double answer,dr,Vcyl;
2111    int i,nord;
2112   
2113    Pi = 4.0*atan(1.0);
2114    lolim = 0.0;
2115    uplim = Pi/2.0;
2116    halfheight = length/2.0;
2117        //    nord=20;
2118    nord=76;
2119    summ = 0.0;
2120   
2121    //do the cylinder orientational average
2122        //    for(i=0;i<nord;i++) {
2123        //            zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
2124        //            yyy = Gauss20Wt[i] * CylKernel(q, dumRad, halfheight, zi);
2125        //            summ += yyy;
2126        //    }
2127    for(i=0;i<nord;i++) {
2128                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
2129                yyy = Gauss76Wt[i] * CylKernel(q, dumRad, halfheight, zi);
2130                summ += yyy;
2131    }
2132    answer = (uplim-lolim)/2.0*summ;
2133    // Multiply by contrast^2
2134    answer *= delrho*delrho;
2135    // don't do the normal scaling to volume here
2136    // instead, multiply by VCyl^2 to get rid of the normalization for this radius of cylinder
2137    Vcyl = Pi*dumRad*dumRad*length;
2138    answer *= Vcyl*Vcyl;
2139   
2140    dr = SchulzPoint_cpr(dumRad,radius,zz);
2141    return(dr*answer);
2142}
2143
2144double
2145SchulzPoint_cpr(double dumRad, double radius, double zz)
2146{
2147    double dr;
2148   
2149    dr = zz*log(dumRad) - gammaln(zz+1.0) + (zz+1.0)*log((zz+1.0)/radius)-(dumRad/radius*(zz+1.0));
2150    return(exp(dr));
2151}
2152
2153
2154double
2155EllipsoidKernel(double qq, double a, double nua, double dum)
2156{
2157    double arg,nu,retval;               //local variables
2158       
2159    nu = nua/a;
2160    arg = qq*a*sqrt(1.0+dum*dum*(nu*nu-1.0));
2161    if (arg == 0.0){
2162        retval =1.0/3.0;
2163    }else{
2164        retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg);
2165    }
2166    retval *= retval;
2167    retval *= 9.0;
2168
2169    return(retval);
2170}//Function EllipsoidKernel()
2171
2172double
2173HolCylKernel(double qq, double rcore, double rshell, double length, double dum)
2174{
2175    double gamma,arg1,arg2,lam1,lam2,psi,sinarg,t2,retval;              //local variables
2176   
2177    gamma = rcore/rshell;
2178    arg1 = qq*rshell*sqrt(1.0-dum*dum);         //1=shell (outer radius)
2179    arg2 = qq*rcore*sqrt(1.0-dum*dum);                  //2=core (inner radius)
2180    if (arg1 == 0.0){
2181        lam1 = 1.0;
2182    }else{
2183        lam1 = 2.0*NR_BessJ1(arg1)/arg1;
2184    }
2185    if (arg2 == 0.0){
2186        lam2 = 1.0;
2187    }else{
2188        lam2 = 2.0*NR_BessJ1(arg2)/arg2;
2189    }
2190    //Todo: Need to check psi behavior as gamma goes to 1.
2191    psi = 1.0/(1.0-gamma*gamma)*(lam1 -  gamma*gamma*lam2);             //SRK 10/19/00
2192    sinarg = qq*length*dum/2.0;
2193    if (sinarg == 0.0){
2194        t2 = 1.0;
2195    }else{
2196        t2 = sin(sinarg)/sinarg;
2197    }
2198
2199    retval = psi*psi*t2*t2;
2200   
2201    return(retval);
2202}//Function HolCylKernel()
2203
2204double
2205PPKernel(double aa, double mu, double uu)
2206{
2207    // mu passed in is really mu*sqrt(1-sig^2)
2208    double arg1,arg2,Pi,tmp1,tmp2;                      //local variables
2209       
2210    Pi = 4.0*atan(1.0);
2211   
2212    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
2213    arg1 = (mu/2.0)*cos(Pi*uu/2.0);
2214    arg2 = (mu*aa/2.0)*sin(Pi*uu/2.0);
2215    if(arg1==0.0) {
2216                tmp1 = 1.0;
2217        } else {
2218                tmp1 = sin(arg1)*sin(arg1)/arg1/arg1;
2219    }
2220
2221    if (arg2==0.0) {
2222                tmp2 = 1.0;
2223        } else {
2224                tmp2 = sin(arg2)*sin(arg2)/arg2/arg2;
2225    }
2226       
2227    return (tmp1*tmp2);
2228       
2229}//Function PPKernel()
2230
2231
2232double
2233TriaxialKernel(double q, double aa, double bb, double cc, double dx, double dy)
2234{
2235       
2236    double arg,val,pi;                  //local variables
2237       
2238    pi = 4.0*atan(1.0);
2239   
2240    arg = aa*aa*cos(pi*dx/2)*cos(pi*dx/2);
2241    arg += bb*bb*sin(pi*dx/2)*sin(pi*dx/2)*(1-dy*dy);
2242    arg += cc*cc*dy*dy;
2243    arg = q*sqrt(arg);
2244    if (arg == 0.0){
2245        val = 1.0; // as arg --> 0, val should go to 1.0
2246    }else{
2247        val = 9.0 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) );
2248    }
2249    return (val);
2250       
2251}//Function TriaxialKernel()
2252
2253
2254double
2255CylKernel(double qq, double rr,double h, double theta)
2256{
2257       
2258        // qq is the q-value for the calculation (1/A)
2259        // rr is the radius of the cylinder (A)
2260        // h is the HALF-LENGTH of the cylinder = L/2 (A)
2261
2262    double besarg,bj,retval,d1,t1,b1,t2,b2,siarg,be,si;          //Local variables
2263
2264
2265    besarg = qq*rr*sin(theta);
2266    siarg = qq * h * cos(theta);
2267    bj =NR_BessJ1(besarg);
2268       
2269        //* Computing 2nd power */
2270    d1 = sin(siarg);
2271    t1 = d1 * d1;
2272        //* Computing 2nd power */
2273    d1 = bj;
2274    t2 = d1 * d1 * 4.0 * sin(theta);
2275        //* Computing 2nd power */
2276    d1 = siarg;
2277    b1 = d1 * d1;
2278        //* Computing 2nd power */
2279    d1 = qq * rr * sin(theta);
2280    b2 = d1 * d1;
2281    if (besarg == 0.0){
2282        be = sin(theta);
2283    }else{
2284        be = t2 / b2;
2285    }
2286    if (siarg == 0.0){
2287        si = 1.0;
2288    }else{
2289        si = t1 / b1;
2290    }
2291    retval = be * si;
2292
2293    return (retval);
2294       
2295}//Function CylKernel()
2296
2297double
2298EllipCylKernel(double qq, double ra,double nu, double theta)
2299{
2300        //this is the function LAMBDA1^2 in Feigin's notation   
2301        // qq is the q-value for the calculation (1/A)
2302        // ra is the transformed radius"a" in Feigin's notation
2303        // nu is the ratio (major radius)/(minor radius) of the Ellipsoid [=] ---
2304        // theta is the dummy variable of the integration
2305       
2306        double retval,arg;               //Local variables
2307       
2308        arg = qq*ra*sqrt((1.0+nu*nu)/2+(1.0-nu*nu)*cos(theta)/2);
2309        if (arg == 0.0){
2310                retval = 1.0;
2311        }else{
2312                retval = 2.0*NR_BessJ1(arg)/arg;
2313        }
2314
2315        //square it
2316        retval *= retval;
2317       
2318    return(retval);
2319       
2320}//Function EllipCylKernel()
2321
2322double NR_BessJ1(double x)
2323{
2324        double ax,z;
2325        double xx,y,ans,ans1,ans2;
2326       
2327        if ((ax=fabs(x)) < 8.0) {
2328                y=x*x;
2329                ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
2330                                                                                                  +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
2331                ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
2332                                                                                           +y*(99447.43394+y*(376.9991397+y*1.0))));
2333                ans=ans1/ans2;
2334        } else {
2335                z=8.0/ax;
2336                y=z*z;
2337                xx=ax-2.356194491;
2338                ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
2339                                                                   +y*(0.2457520174e-5+y*(-0.240337019e-6))));
2340                ans2=0.04687499995+y*(-0.2002690873e-3
2341                                                          +y*(0.8449199096e-5+y*(-0.88228987e-6
2342                                                                                                         +y*0.105787412e-6)));
2343                ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
2344                if (x < 0.0) ans = -ans;
2345        }
2346       
2347        return(ans);
2348}
2349
2350/*      Lamellar_ParaCrystal - Pedersen's model
2351
2352*/
2353double
2354Lamellar_ParaCrystal(double w[], double q)
2355{
2356//       Input (fitting) variables are:
2357        //[0] scale factor
2358        //[1]   thickness
2359        //[2]   number of layers
2360        //[3]   spacing between layers
2361        //[4]   polydispersity of spacing
2362        //[5] SLD lamellar
2363        //[6] SLD solvent
2364        //[7] incoherent background
2365//      give them nice names
2366        double inten,qval,scale,th,nl,davg,pd,contr,bkg,xn;
2367        double xi,ww,Pbil,Znq,Snq,an,sldLayer,sldSolvent,pi;
2368        long n1,n2;
2369       
2370        pi = 4.0*atan(1.0);
2371        scale = w[0];
2372        th = w[1];
2373        nl = w[2];
2374        davg = w[3];
2375        pd = w[4];
2376        sldLayer = w[5];
2377        sldSolvent = w[6];
2378        bkg = w[7];
2379       
2380        contr = w[5] - w[6];
2381        qval = q;
2382               
2383        //get the fractional part of nl, to determine the "mixing" of N's
2384       
2385        n1 = trunc(nl);         //rounds towards zero
2386        n2 = n1 + 1;
2387        xn = (double)n2 - nl;                   //fractional contribution of n1
2388       
2389        ww = exp(-qval*qval*pd*pd*davg*davg/2.0);
2390       
2391        //calculate the n1 contribution
2392        an = paraCryst_an(ww,qval,davg,n1);
2393        Snq = paraCryst_sn(ww,qval,davg,n1,an);
2394       
2395        Znq = xn*Snq;
2396       
2397        //calculate the n2 contribution
2398        an = paraCryst_an(ww,qval,davg,n2);
2399        Snq = paraCryst_sn(ww,qval,davg,n2,an);
2400       
2401        Znq += (1.0-xn)*Snq;
2402       
2403        //and the independent contribution
2404        Znq += (1.0-ww*ww)/(1.0+ww*ww-2.0*ww*cos(qval*davg));
2405       
2406        //the limit when NL approaches infinity
2407//      Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg))
2408       
2409        xi = th/2.0;            //use 1/2 the bilayer thickness
2410        Pbil = (sin(qval*xi)/(qval*xi))*(sin(qval*xi)/(qval*xi));
2411       
2412        inten = 2.0*pi*contr*contr*Pbil*Znq/(qval*qval);
2413        inten *= 1.0e8;
2414       
2415        return(scale*inten+bkg);
2416}
2417
2418// functions for the lamellar paracrystal model
2419double
2420paraCryst_sn(double ww, double qval, double davg, long nl, double an) {
2421       
2422        double Snq;
2423       
2424        Snq = an/( (double)nl*pow((1.0+ww*ww-2.0*ww*cos(qval*davg)),2) );
2425       
2426        return(Snq);
2427}
2428
2429
2430double
2431paraCryst_an(double ww, double qval, double davg, long nl) {
2432       
2433        double an;
2434       
2435        an = 4.0*ww*ww - 2.0*(ww*ww*ww+ww)*cos(qval*davg);
2436        an -= 4.0*pow(ww,(nl+2))*cos((double)nl*qval*davg);
2437        an += 2.0*pow(ww,(nl+3))*cos((double)(nl-1)*qval*davg);
2438        an += 2.0*pow(ww,(nl+1))*cos((double)(nl+1)*qval*davg);
2439       
2440        return(an);
2441}
2442
2443
2444/*      Spherocylinder  :
2445
2446Uses 76 pt Gaussian quadrature for both integrals
2447*/
2448double
2449Spherocylinder(double w[], double x)
2450{
2451        int i,j;
2452        double Pi;
2453        double scale,contr,bkg,sldc,slds;
2454        double len,rad,hDist,endRad;
2455        int nordi=76;                   //order of integration
2456        int nordj=76;
2457        double va,vb;           //upper and lower integration limits
2458        double summ,zi,yyy,answer;                      //running tally of integration
2459        double summj,vaj,vbj,zij;                       //for the inner integration
2460        double SphCyl_tmp[7],arg1,arg2,inner,be;
2461       
2462       
2463        scale = w[0];
2464        rad = w[1];
2465        len = w[2];
2466        sldc = w[3];
2467        slds = w[4];
2468        bkg = w[5];
2469       
2470        SphCyl_tmp[0] = w[0];
2471        SphCyl_tmp[1] = w[1];
2472        SphCyl_tmp[2] = w[2];
2473        SphCyl_tmp[3] = w[1];           //end radius is same as cylinder radius
2474        SphCyl_tmp[4] = w[3];
2475        SphCyl_tmp[5] = w[4];
2476        SphCyl_tmp[6] = w[5];
2477       
2478        hDist = 0;              //by definition for this model
2479        endRad = rad;
2480       
2481        contr = sldc-slds;
2482       
2483        Pi = 4.0*atan(1.0);
2484        va = 0.0;
2485        vb = Pi/2.0;            //orintational average, outer integral
2486        vaj = -1.0*hDist/endRad;
2487        vbj = 1.0;              //endpoints of inner integral
2488
2489        summ = 0.0;                     //initialize intergral
2490
2491        for(i=0;i<nordi;i++) {
2492                //setup inner integral over the ellipsoidal cross-section
2493                summj=0.0;
2494                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
2495               
2496                for(j=0;j<nordj;j++) {
2497                        //20 gauss points for the inner integral
2498                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
2499                        yyy = Gauss76Wt[j] * SphCyl_kernel(SphCyl_tmp,x,zij,zi);
2500                        summj += yyy;
2501                }
2502                //now calculate the value of the inner integral
2503                inner = (vbj-vaj)/2.0*summj;
2504                inner *= 4.0*Pi*endRad*endRad*endRad;
2505               
2506                //now calculate outer integrand
2507                arg1 = x*len/2.0*cos(zi);
2508                arg2 = x*rad*sin(zi);
2509                yyy = inner;
2510
2511                if(arg2 == 0) {
2512                        be = 0.5;
2513                } else {
2514                        be = NR_BessJ1(arg2)/arg2;
2515                }
2516               
2517                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
2518                        yyy += Pi*rad*rad*len*2.0*be;
2519                } else {
2520                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
2521                }
2522                yyy *= yyy;
2523                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
2524                yyy *= Gauss76Wt[i];
2525                summ += yyy;
2526        }               //final scaling is done at the end of the function, after the NT_FP64 case
2527       
2528        answer = (vb-va)/2.0*summ;
2529
2530        answer /= Pi*rad*rad*len + Pi*4.0*endRad*endRad*endRad/3.0;             //divide by volume
2531        answer *= 1.0e8;                //convert to cm^-1
2532        answer *= contr*contr;
2533        answer *= scale;
2534        answer += bkg;
2535               
2536        return answer;
2537}
2538
2539
2540// inner integral of the sphereocylinder model, special case of hDist = 0
2541//
2542double
2543SphCyl_kernel(double w[], double x, double tt, double theta) {
2544
2545        double val,arg1,arg2;
2546        double scale,bkg,sldc,slds;
2547        double len,rad,hDist,endRad,be;
2548        scale = w[0];
2549        rad = w[1];
2550        len = w[2];
2551        endRad = w[3];
2552        sldc = w[4];
2553        slds = w[5];
2554        bkg = w[6];
2555       
2556        hDist = 0.0;
2557               
2558        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0);
2559        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt);
2560       
2561        if(arg2 == 0) {
2562                be = 0.5;
2563        } else {
2564                be = NR_BessJ1(arg2)/arg2;
2565        }
2566        val = cos(arg1)*(1.0-tt*tt)*be;
2567       
2568        return(val);
2569}
2570
2571
2572/*      Convex Lens  : special case where L ~ 0 and hDist < 0
2573
2574Uses 76 pt Gaussian quadrature for both integrals
2575*/
2576double
2577ConvexLens(double w[], double x)
2578{
2579        int i,j;
2580        double Pi;
2581        double scale,contr,bkg,sldc,slds;
2582        double len,rad,hDist,endRad;
2583        int nordi=76;                   //order of integration
2584        int nordj=76;
2585        double va,vb;           //upper and lower integration limits
2586        double summ,zi,yyy,answer;                      //running tally of integration
2587        double summj,vaj,vbj,zij;                       //for the inner integration
2588        double CLens_tmp[7],arg1,arg2,inner,hh,be;
2589       
2590       
2591        scale = w[0];
2592        rad = w[1];
2593//      len = w[2]
2594        endRad = w[2];
2595        sldc = w[3];
2596        slds = w[4];
2597        bkg = w[5];
2598       
2599        len = 0.01;
2600       
2601        CLens_tmp[0] = w[0];
2602        CLens_tmp[1] = w[1];
2603        CLens_tmp[2] = len;                     //length is some small number, essentially zero
2604        CLens_tmp[3] = w[2];
2605        CLens_tmp[4] = w[3];
2606        CLens_tmp[5] = w[4];
2607        CLens_tmp[6] = w[5];
2608               
2609        hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad));         //by definition for this model
2610       
2611        contr = sldc-slds;
2612       
2613        Pi = 4.0*atan(1.0);
2614        va = 0.0;
2615        vb = Pi/2.0;            //orintational average, outer integral
2616        vaj = -1.0*hDist/endRad;
2617        vbj = 1.0;              //endpoints of inner integral
2618
2619        summ = 0.0;                     //initialize intergral
2620
2621        for(i=0;i<nordi;i++) {
2622                //setup inner integral over the ellipsoidal cross-section
2623                summj=0.0;
2624                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
2625               
2626                for(j=0;j<nordj;j++) {
2627                        //20 gauss points for the inner integral
2628                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
2629                        yyy = Gauss76Wt[j] * ConvLens_kernel(CLens_tmp,x,zij,zi);
2630                        summj += yyy;
2631                }
2632                //now calculate the value of the inner integral
2633                inner = (vbj-vaj)/2.0*summj;
2634                inner *= 4.0*Pi*endRad*endRad*endRad;
2635               
2636                //now calculate outer integrand
2637                arg1 = x*len/2.0*cos(zi);
2638                arg2 = x*rad*sin(zi);
2639                yyy = inner;
2640               
2641                if(arg2 == 0) {
2642                        be = 0.5;
2643                } else {
2644                        be = NR_BessJ1(arg2)/arg2;
2645                }
2646               
2647                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
2648                        yyy += Pi*rad*rad*len*2.0*be;
2649                } else {
2650                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
2651                }
2652                yyy *= yyy;
2653                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
2654                yyy *= Gauss76Wt[i];
2655                summ += yyy;
2656        }               //final scaling is done at the end of the function, after the NT_FP64 case
2657       
2658        answer = (vb-va)/2.0*summ;
2659
2660        hh = fabs(hDist);               //need positive value for spherical cap volume
2661        answer /= 2.0*(1.0/3.0*Pi*(endRad-hh)*(endRad-hh)*(2.0*endRad+hh));             //divide by volume
2662        answer *= 1.0e8;                //convert to cm^-1
2663        answer *= contr*contr;
2664        answer *= scale;
2665        answer += bkg;
2666               
2667        return answer;
2668}
2669
2670/*      Capped Cylinder  : special case where L is nonzero and hDist < 0
2671
2672 -- uses the same Kernel as the Convex Lens
2673 
2674Uses 76 pt Gaussian quadrature for both integrals
2675*/
2676double
2677CappedCylinder(double w[], double x)
2678{
2679        int i,j;
2680        double Pi;
2681        double scale,contr,bkg,sldc,slds;
2682        double len,rad,hDist,endRad;
2683        int nordi=76;                   //order of integration
2684        int nordj=76;
2685        double va,vb;           //upper and lower integration limits
2686        double summ,zi,yyy,answer;                      //running tally of integration
2687        double summj,vaj,vbj,zij;                       //for the inner integration
2688        double arg1,arg2,inner,hh,be;
2689       
2690       
2691        scale = w[0];
2692        rad = w[1];
2693        len = w[2];
2694        endRad = w[3];
2695        sldc = w[4];
2696        slds = w[5];
2697        bkg = w[6];
2698               
2699        hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad));         //by definition for this model
2700       
2701        contr = sldc-slds;
2702       
2703        Pi = 4.0*atan(1.0);
2704        va = 0.0;
2705        vb = Pi/2.0;            //orintational average, outer integral
2706        vaj = -1.0*hDist/endRad;
2707        vbj = 1.0;              //endpoints of inner integral
2708
2709        summ = 0.0;                     //initialize intergral
2710
2711        for(i=0;i<nordi;i++) {
2712                //setup inner integral over the ellipsoidal cross-section
2713                summj=0.0;
2714                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
2715               
2716                for(j=0;j<nordj;j++) {
2717                        //20 gauss points for the inner integral
2718                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
2719                        yyy = Gauss76Wt[j] * ConvLens_kernel(w,x,zij,zi);               //uses the same kernel as ConvexLens, except here L != 0
2720                        summj += yyy;
2721                }
2722                //now calculate the value of the inner integral
2723                inner = (vbj-vaj)/2.0*summj;
2724                inner *= 4.0*Pi*endRad*endRad*endRad;
2725               
2726                //now calculate outer integrand
2727                arg1 = x*len/2.0*cos(zi);
2728                arg2 = x*rad*sin(zi);
2729                yyy = inner;
2730               
2731                if(arg2 == 0) {
2732                        be = 0.5;
2733                } else {
2734                        be = NR_BessJ1(arg2)/arg2;
2735                }
2736               
2737                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
2738                        yyy += Pi*rad*rad*len*2.0*be;
2739                } else {
2740                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
2741                }
2742               
2743               
2744               
2745                yyy *= yyy;
2746                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
2747                yyy *= Gauss76Wt[i];
2748                summ += yyy;
2749        }               //final scaling is done at the end of the function, after the NT_FP64 case
2750       
2751        answer = (vb-va)/2.0*summ;
2752
2753        hh = fabs(hDist);               //need positive value for spherical cap volume
2754        answer /= Pi*rad*rad*len + 2.0*(1.0/3.0*Pi*(endRad-hh)*(endRad-hh)*(2.0*endRad+hh));            //divide by volume
2755        answer *= 1.0e8;                //convert to cm^-1
2756        answer *= contr*contr;
2757        answer *= scale;
2758        answer += bkg;
2759               
2760        return answer;
2761}
2762
2763
2764
2765// inner integral of the ConvexLens model, special case where L ~ 0 and hDist < 0
2766//
2767double
2768ConvLens_kernel(double w[], double x, double tt, double theta) {
2769
2770        double val,arg1,arg2;
2771        double scale,bkg,sldc,slds;
2772        double len,rad,hDist,endRad,be;
2773        scale = w[0];
2774        rad = w[1];
2775        len = w[2];
2776        endRad = w[3];
2777        sldc = w[4];
2778        slds = w[5];
2779        bkg = w[6];
2780       
2781        hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad));
2782               
2783        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0);
2784        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt);
2785       
2786        if(arg2 == 0) {
2787                be = 0.5;
2788        } else {
2789                be = NR_BessJ1(arg2)/arg2;
2790        }
2791       
2792        val = cos(arg1)*(1.0-tt*tt)*be;
2793       
2794        return(val);
2795}
2796
2797
2798/*      Dumbbell  : special case where L ~ 0 and hDist > 0
2799
2800Uses 76 pt Gaussian quadrature for both integrals
2801*/
2802double
2803Dumbbell(double w[], double x)
2804{
2805        int i,j;
2806        double Pi;
2807        double scale,contr,bkg,sldc,slds;
2808        double len,rad,hDist,endRad;
2809        int nordi=76;                   //order of integration
2810        int nordj=76;
2811        double va,vb;           //upper and lower integration limits
2812        double summ,zi,yyy,answer;                      //running tally of integration
2813        double summj,vaj,vbj,zij;                       //for the inner integration
2814        double Dumb_tmp[7],arg1,arg2,inner,be;
2815       
2816       
2817        scale = w[0];
2818        rad = w[1];
2819//      len = w[2]
2820        endRad = w[2];
2821        sldc = w[3];
2822        slds = w[4];
2823        bkg = w[5];
2824       
2825        len = 0.01;
2826       
2827        Dumb_tmp[0] = w[0];
2828        Dumb_tmp[1] = w[1];
2829        Dumb_tmp[2] = len;              //length is some small number, essentially zero
2830        Dumb_tmp[3] = w[2];
2831        Dumb_tmp[4] = w[3];
2832        Dumb_tmp[5] = w[4];
2833        Dumb_tmp[6] = w[5];
2834                       
2835        hDist = sqrt(fabs(endRad*endRad-rad*rad));              //by definition for this model
2836       
2837        contr = sldc-slds;
2838       
2839        Pi = 4.0*atan(1.0);
2840        va = 0.0;
2841        vb = Pi/2.0;            //orintational average, outer integral
2842        vaj = -1.0*hDist/endRad;
2843        vbj = 1.0;              //endpoints of inner integral
2844
2845        summ = 0.0;                     //initialize intergral
2846
2847        for(i=0;i<nordi;i++) {
2848                //setup inner integral over the ellipsoidal cross-section
2849                summj=0.0;
2850                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
2851               
2852                for(j=0;j<nordj;j++) {
2853                        //20 gauss points for the inner integral
2854                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
2855                        yyy = Gauss76Wt[j] * Dumb_kernel(Dumb_tmp,x,zij,zi);
2856                        summj += yyy;
2857                }
2858                //now calculate the value of the inner integral
2859                inner = (vbj-vaj)/2.0*summj;
2860                inner *= 4.0*Pi*endRad*endRad*endRad;
2861               
2862                //now calculate outer integrand
2863                arg1 = x*len/2.0*cos(zi);
2864                arg2 = x*rad*sin(zi);
2865                yyy = inner;
2866               
2867                if(arg2 == 0) {
2868                        be = 0.5;
2869                } else {
2870                        be = NR_BessJ1(arg2)/arg2;
2871                }
2872               
2873                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
2874                        yyy += Pi*rad*rad*len*2.0*be;
2875                } else {
2876                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
2877                }
2878                yyy *= yyy;
2879                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
2880                yyy *= Gauss76Wt[i];
2881                summ += yyy;
2882        }               //final scaling is done at the end of the function, after the NT_FP64 case
2883       
2884        answer = (vb-va)/2.0*summ;
2885
2886        answer /= 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*hDist*hDist/3.0);              //divide by volume
2887        answer *= 1.0e8;                //convert to cm^-1
2888        answer *= contr*contr;
2889        answer *= scale;
2890        answer += bkg;
2891               
2892        return answer;
2893}
2894
2895
2896/*      Barbell  : "normal" case where L is nonzero 0 and hDist > 0
2897
2898-- uses the same kernel as the Dumbbell case
2899
2900Uses 76 pt Gaussian quadrature for both integrals
2901*/
2902double
2903Barbell(double w[], double x)
2904{
2905        int i,j;
2906        double Pi;
2907        double scale,contr,bkg,sldc,slds;
2908        double len,rad,hDist,endRad;
2909        int nordi=76;                   //order of integration
2910        int nordj=76;
2911        double va,vb;           //upper and lower integration limits
2912        double summ,zi,yyy,answer;                      //running tally of integration
2913        double summj,vaj,vbj,zij;                       //for the inner integration
2914        double arg1,arg2,inner,be;
2915       
2916       
2917        scale = w[0];
2918        rad = w[1];
2919        len = w[2];
2920        endRad = w[3];
2921        sldc = w[4];
2922        slds = w[5];
2923        bkg = w[6];
2924                       
2925        hDist = sqrt(fabs(endRad*endRad-rad*rad));              //by definition for this model
2926       
2927        contr = sldc-slds;
2928       
2929        Pi = 4.0*atan(1.0);
2930        va = 0.0;
2931        vb = Pi/2.0;            //orintational average, outer integral
2932        vaj = -1.0*hDist/endRad;
2933        vbj = 1.0;              //endpoints of inner integral
2934
2935        summ = 0.0;                     //initialize intergral
2936
2937        for(i=0;i<nordi;i++) {
2938                //setup inner integral over the ellipsoidal cross-section
2939                summj=0.0;
2940                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
2941               
2942                for(j=0;j<nordj;j++) {
2943                        //20 gauss points for the inner integral
2944                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
2945                        yyy = Gauss76Wt[j] * Dumb_kernel(w,x,zij,zi);           //uses the same Kernel as the Dumbbell, here L>0
2946                        summj += yyy;
2947                }
2948                //now calculate the value of the inner integral
2949                inner = (vbj-vaj)/2.0*summj;
2950                inner *= 4.0*Pi*endRad*endRad*endRad;
2951               
2952                //now calculate outer integrand
2953                arg1 = x*len/2.0*cos(zi);
2954                arg2 = x*rad*sin(zi);
2955                yyy = inner;
2956               
2957                if(arg2 == 0) {
2958                        be = 0.5;
2959                } else {
2960                        be = NR_BessJ1(arg2)/arg2;
2961                }
2962               
2963                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
2964                        yyy += Pi*rad*rad*len*2.0*be;
2965                } else {
2966                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
2967                }
2968                yyy *= yyy;
2969                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
2970                yyy *= Gauss76Wt[i];
2971                summ += yyy;
2972        }               //final scaling is done at the end of the function, after the NT_FP64 case
2973       
2974        answer = (vb-va)/2.0*summ;
2975
2976        answer /= Pi*rad*rad*len + 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*hDist*hDist/3.0);             //divide by volume
2977        answer *= 1.0e8;                //convert to cm^-1
2978        answer *= contr*contr;
2979        answer *= scale;
2980        answer += bkg;
2981               
2982        return answer;
2983}
2984
2985
2986
2987// inner integral of the Dumbbell model, special case where L ~ 0 and hDist > 0
2988//
2989// inner integral of the Barbell model if L is nonzero
2990//
2991double
2992Dumb_kernel(double w[], double x, double tt, double theta) {
2993
2994        double val,arg1,arg2;
2995        double scale,bkg,sldc,slds;
2996        double len,rad,hDist,endRad,be;
2997        scale = w[0];
2998        rad = w[1];
2999        len = w[2];
3000        endRad = w[3];
3001        sldc = w[4];
3002        slds = w[5];
3003        bkg = w[6];
3004       
3005        hDist = sqrt(fabs(endRad*endRad-rad*rad));
3006               
3007        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0);
3008        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt);
3009       
3010        if(arg2 == 0) {
3011                be = 0.5;
3012        } else {
3013                be = NR_BessJ1(arg2)/arg2;
3014        }
3015        val = cos(arg1)*(1.0-tt*tt)*be;
3016       
3017        return(val);
3018}
3019
3020double PolyCoreBicelle(double dp[], double q)
3021{
3022        int i;
3023        int nord = 20;
3024        double scale, length, sigma, bkg, radius, radthick, facthick;
3025        double rhoc, rhoh, rhor, rhosolv;
3026        double answer, Vpoly;
3027        double Pi,lolim,uplim,summ,yyy,rad,AR,Rsqr,Rsqrsumm,Rsqryyy;
3028       
3029        scale = dp[0];
3030        radius = dp[1];
3031        sigma = dp[2];                          //sigma is the standard mean deviation
3032        length = dp[3];
3033        radthick = dp[4];
3034        facthick= dp[5];
3035        rhoc = dp[6];
3036        rhoh = dp[7];
3037        rhor=dp[8];
3038        rhosolv = dp[9];
3039        bkg = dp[10];
3040       
3041        Pi = 4.0*atan(1.0);
3042       
3043        lolim = exp(log(radius)-(4.*sigma));
3044        if (lolim<0.0) {
3045                lolim=0.0;              //to avoid numerical error when  va<0 (-ve r value)
3046        }
3047        uplim = exp(log(radius)+(4.*sigma));
3048       
3049        summ = 0.0;
3050        Rsqrsumm = 0.0;
3051       
3052        for(i=0;i<nord;i++) {
3053                rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
3054                AR=(1.0/(rad*sigma*sqrt(2.0*Pi)))*exp(-(0.5*((log(radius/rad))/sigma)*((log(radius/rad))/sigma)));
3055                yyy = AR* Gauss20Wt[i] * BicelleIntegration(q,rad,radthick,facthick,rhoc,rhoh,rhor,rhosolv,length);
3056                Rsqryyy= Gauss20Wt[i] * AR * (rad+radthick)*(rad+radthick);             //SRK normalize to total dimensions
3057                summ += yyy;
3058                Rsqrsumm += Rsqryyy;
3059        }
3060       
3061        answer = (uplim-lolim)/2.0*summ;
3062        Rsqr = (uplim-lolim)/2.0*Rsqrsumm;
3063        //normalize by average cylinder volume
3064        Vpoly = Pi*Rsqr*(length+2*facthick);
3065        answer /= Vpoly;
3066        //convert to [cm-1]
3067        answer *= 1.0e8;
3068        //Scale
3069        answer *= scale;
3070        // add in the background
3071        answer += bkg;
3072               
3073        return answer;
3074       
3075}
3076
3077double
3078BicelleIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length){
3079
3080        double answer,halfheight,Pi;
3081        double lolim,uplim,summ,yyy,zi;
3082        int nord,i;
3083       
3084        // set up the integration end points
3085        Pi = 4.0*atan(1.0);
3086        nord = 76;
3087        lolim = 0.0;
3088        uplim = Pi/2;
3089        halfheight = length/2.0;
3090       
3091        summ = 0.0;                             // initialize integral
3092        i=0;
3093        for(i=0;i<nord;i++) {
3094                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
3095                yyy = Gauss76Wt[i] * BicelleKernel(qq, rad, radthick, facthick, rhoc, rhoh, rhor,rhosolv, halfheight, zi);
3096                summ += yyy;
3097        }
3098       
3099        // calculate value of integral to return
3100        answer = (uplim-lolim)/2.0*summ;
3101        return(answer); 
3102}
3103
3104double
3105BicelleKernel(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length, double dum)
3106{
3107        double dr1,dr2,dr3;
3108        double besarg1,besarg2;
3109        double vol1,vol2,vol3;
3110        double sinarg1,sinarg2;
3111        double t1,t2,t3;
3112        double retval,si1,si2,be1,be2;
3113       
3114        double Pi = 4.0*atan(1.0);
3115       
3116        dr1 = rhoc-rhoh;
3117        dr2 = rhor-rhosolv;
3118        dr3=  rhoh-rhor;
3119        vol1 = Pi*rad*rad*(2.0*length);
3120        vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick);
3121        vol3= Pi*(rad)*(rad)*(2.0*length+2.0*facthick);
3122        besarg1 = qq*rad*sin(dum);
3123        besarg2 = qq*(rad+radthick)*sin(dum);
3124        sinarg1 = qq*length*cos(dum);
3125        sinarg2 = qq*(length+facthick)*cos(dum);
3126       
3127        if(besarg1 == 0) {
3128                be1 = 0.5;
3129        } else {
3130                be1 = NR_BessJ1(besarg1)/besarg1;
3131        }
3132        if(besarg2 == 0) {
3133                be2 = 0.5;
3134        } else {
3135                be2 = NR_BessJ1(besarg2)/besarg2;
3136        }       
3137        if(sinarg1 == 0) {
3138                si1 = 1.0;
3139        } else {
3140                si1 = sin(sinarg1)/sinarg1;
3141        }
3142        if(sinarg2 == 0) {
3143                si2 = 1.0;
3144        } else {
3145                si2 = sin(sinarg2)/sinarg2;
3146        }
3147        t1 = 2.0*vol1*dr1*si1*be1;
3148        t2 = 2.0*vol2*dr2*si2*be2;
3149        t3 = 2.0*vol3*dr3*si2*be1;
3150       
3151        retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum);
3152        return(retval);
3153       
3154}
3155
3156
3157double
3158CSPPKernel(double dp[], double mu, double uu)
3159{       
3160        double aa,bb,cc, ta,tb,tc; 
3161        double Vin,Vot,V1,V2;
3162        double rhoA,rhoB,rhoC, rhoP, rhosolv;
3163        double dr0, drA,drB, drC;
3164        double arg1,arg2,arg3,arg4,t1,t2, t3, t4;
3165        double Pi,retVal;
3166
3167        aa = dp[1];
3168        bb = dp[2];
3169        cc = dp[3];
3170        ta = dp[4];
3171        tb = dp[5];
3172        tc = dp[6];
3173        rhoA=dp[7];
3174        rhoB=dp[8];
3175        rhoC=dp[9];
3176        rhoP=dp[10];
3177        rhosolv=dp[11];
3178        dr0=rhoP-rhosolv;
3179        drA=rhoA-rhosolv;
3180        drB=rhoB-rhosolv;
3181        drC=rhoC-rhosolv; 
3182        Vin=(aa*bb*cc);
3183        Vot=(aa*bb*cc+2.0*ta*bb*cc+2.0*aa*tb*cc+2.0*aa*bb*tc);
3184        V1=(2.0*ta*bb*cc);   //  incorrect V1 (aa*bb*cc+2*ta*bb*cc)
3185        V2=(2.0*aa*tb*cc);  // incorrect V2(aa*bb*cc+2*aa*tb*cc)
3186        aa = aa/bb;
3187        ta=(aa+2.0*ta)/bb;
3188        tb=(aa+2.0*tb)/bb;
3189       
3190        Pi = 4.0*atan(1.0);
3191       
3192        arg1 = (mu*aa/2.0)*sin(Pi*uu/2.0);
3193        arg2 = (mu/2.0)*cos(Pi*uu/2.0);
3194        arg3=  (mu*ta/2.0)*sin(Pi*uu/2.0);
3195        arg4=  (mu*tb/2.0)*cos(Pi*uu/2.0);
3196                         
3197        if(arg1==0.0){
3198                t1 = 1.0;
3199        } else {
3200                t1 = (sin(arg1)/arg1);                //defn for CSPP model sin(arg1)/arg1    test:  (sin(arg1)/arg1)*(sin(arg1)/arg1)   
3201        }
3202        if(arg2==0.0){
3203                t2 = 1.0;
3204        } else {
3205                t2 = (sin(arg2)/arg2);           //defn for CSPP model sin(arg2)/arg2   test: (sin(arg2)/arg2)*(sin(arg2)/arg2)   
3206        }       
3207        if(arg3==0.0){
3208                t3 = 1.0;
3209        } else {
3210                t3 = sin(arg3)/arg3;
3211        }
3212        if(arg4==0.0){
3213                t4 = 1.0;
3214        } else {
3215                t4 = sin(arg4)/arg4;
3216        }
3217        retVal =( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 )*( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 );   //  correct FF : square of sum of phase factors
3218        return(retVal); 
3219
3220}
3221
3222/*      CSParallelepiped  :  calculates the form factor of a Parallelepiped with a core-shell structure
3223 -- different SLDs can be used for the face and rim
3224
3225Uses 76 pt Gaussian quadrature for both integrals
3226*/
3227double
3228CSParallelepiped(double dp[], double q)
3229{
3230        int i,j;
3231        double scale,aa,bb,cc,ta,tb,tc,rhoA,rhoB,rhoC,rhoP,rhosolv,bkg;         //local variables of coefficient wave
3232        int nordi=76;                   //order of integration
3233        int nordj=76;
3234        double va,vb;           //upper and lower integration limits
3235        double summ,yyy,answer;                 //running tally of integration
3236        double summj,vaj,vbj;                   //for the inner integration
3237        double mu,mudum,arg,sigma,uu,vol;
3238       
3239       
3240        //      Pi = 4.0*atan(1.0);
3241        va = 0.0;
3242        vb = 1.0;               //orintational average, outer integral
3243        vaj = 0.0;
3244        vbj = 1.0;              //endpoints of inner integral
3245       
3246        summ = 0.0;                     //initialize intergral
3247       
3248        scale = dp[0];
3249        aa = dp[1];
3250        bb = dp[2];
3251        cc = dp[3];
3252        ta  = dp[4];
3253        tb  = dp[5];
3254        tc  = dp[6];   // is 0 at the moment 
3255        rhoA=dp[7];   //rim A SLD
3256        rhoB=dp[8];   //rim B SLD
3257        rhoC=dp[9];    //rim C SLD
3258        rhoP = dp[10];   //Parallelpiped core SLD
3259        rhosolv=dp[11];  // Solvent SLD
3260        bkg = dp[12];
3261       
3262        mu = q*bb;
3263        vol = aa*bb*cc+2.0*ta*bb*cc+2.0*aa*tb*cc+2.0*aa*bb*tc;          //calculate volume before rescaling
3264       
3265        // do the rescaling here, not in the kernel
3266        // normalize all WRT bb
3267        aa = aa/bb;
3268        cc = cc/bb;
3269       
3270        for(i=0;i<nordi;i++) {
3271                //setup inner integral over the ellipsoidal cross-section
3272                summj=0.0;
3273                sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;          //the outer dummy
3274               
3275                for(j=0;j<nordj;j++) {
3276                        //76 gauss points for the inner integral
3277                        uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;         //the inner dummy
3278                        mudum = mu*sqrt(1.0-sigma*sigma);
3279                        yyy = Gauss76Wt[j] * CSPPKernel(dp,mudum,uu);
3280                        summj += yyy;
3281                }
3282                //now calculate the value of the inner integral
3283                answer = (vbj-vaj)/2.0*summj;
3284               
3285                //finish the outer integral cc already scaled
3286                arg = mu*cc*sigma/2.0;
3287                if ( arg == 0.0 ) {
3288                        answer *= 1.0;
3289                } else {
3290                        answer *= sin(arg)*sin(arg)/arg/arg;
3291                }
3292               
3293                //now sum up the outer integral
3294                yyy = Gauss76Wt[i] * answer;
3295                summ += yyy;
3296        }               //final scaling is done at the end of the function, after the NT_FP64 case
3297       
3298        answer = (vb-va)/2.0*summ;
3299
3300        //normalize by volume
3301        answer /= vol;
3302        //convert to [cm-1]
3303        answer *= 1.0e8;
3304        //Scale
3305        answer *= scale;
3306        // add in the background
3307        answer += bkg;
3308       
3309        return answer;
3310}
3311
3312
3313/*
3314 * Massive rectangular prism (a <= b <=c)
3315 * using eqns. (12)+(16) from Nayuk & Huber, Z. Phys. Chem. 226, 837 (2012).
3316 * It is totally equivalent to the Parallelepiped model in this library,
3317 * but instead of using as input the a, b, and c sides of the prism,
3318 * it uses a and the ratios (b/a) and (c/a).
3319 * This allows to keep the shape when using the polydispersity.
3320 */
3321
3322double
3323RectangularPrism(double dp[], double q)
3324{
3325    int i, j;
3326    double scale, aa, bb, cc, delrho, bkg;             
3327    int nordi=76;                               //order of integration
3328    int nordj=76;
3329    double Pi;
3330    double sumi, sumj;
3331    double v1a, v1b, v2a, v2b;          // upper and lower integration limits
3332    double answer;
3333    double theta, phi, arg, termA, termB, termC, AP, ahalf, bhalf, chalf;
3334    double vol, sldp, sld;
3335
3336    Pi = 4.0*atan(1.0);
3337
3338    scale = dp[0];     
3339    aa = dp[1];          // parameter short_side
3340    bb = aa * dp[2];     // parameter b/a ratio
3341    cc = aa * dp[3];     // parameter c/a ratio
3342    sldp = dp[4];        // scattering length density of the object
3343    sld = dp[5];         // scattering length density of the solvent
3344    delrho = sldp - sld;
3345    bkg = dp[6];
3346
3347    vol = aa*bb*cc;
3348    ahalf = 0.5 * aa;
3349    bhalf = 0.5 * bb;
3350    chalf = 0.5 * cc;
3351
3352   //Integration limits to use in Gaussian quadrature
3353    v1a = 0.;
3354    v1b = Pi/2.;  //theta integration limits
3355    v2a = 0.;
3356    v2b = Pi/2.;  //phi integration limits
3357
3358    sumi = 0.0;
3359    for(i=0;i<nordi;i++) {
3360
3361            theta = ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b )/2.0;  //z values for gaussian quadrature
3362
3363            arg = q*chalf*cos(theta);
3364            if (fabs(arg) > 1.e-16) {termC = sin(arg) / arg;} else {termC = 1.0;} 
3365
3366            sumj = 0.0;
3367            for(j=0;j<nordj;j++) {
3368
3369                phi = ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b )/2.0; //z values for gaussian quadrature
3370
3371                // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0
3372
3373                arg = q*ahalf*sin(theta)*sin(phi); 
3374                if (fabs(arg) > 1.e-16) {termA = sin(arg) / arg;} else {termA = 1.0;}
3375               
3376                arg = q*bhalf*sin(theta)*cos(phi); 
3377                if (fabs(arg) > 1.e-16) {termB = sin(arg) / arg;} else {termB = 1.0;}     
3378               
3379                AP = termA * termB * termC; 
3380
3381                sumj += Gauss76Wt[j] * (AP*AP);
3382
3383            }
3384
3385            sumj = (v2b-v2a) * sumj / 2.0;
3386            sumi += Gauss76Wt[i] * sumj * sin(theta);
3387
3388    }
3389
3390    answer = (v1b-v1a)/2.0*sumi;
3391
3392    // Normalize by Pi (Eqn. 16).
3393    // The term (ABC)^2 does not appear because it was introduced before from the defs of termA, termB, termC
3394    // The factor 2 appears because the theta integral has been defined between 0 and pi/2, instead of 0 to pi
3395    answer = 2.0 * answer / Pi; //Form factor P(q)
3396
3397    // Multiply by contrast^2
3398    answer *= delrho*delrho;
3399
3400    //normalize by volume
3401    answer *= vol;
3402       
3403    //convert to [cm-1]
3404    answer *= 1.0e8;
3405
3406    //Scale
3407    answer *= scale;
3408
3409    // add in the background
3410    answer += bkg;
3411
3412    return answer;
3413}
3414
3415
3416/*
3417 * Hollow rectangular prism (a <= b <=c) with infinitely thin walls
3418 * using eqn. (7-11) from Nayuk & Huber, Z. Phys. Chem. 226, 837 (2012)
3419 */
3420
3421double
3422RectangularHollowPrismInfThinWalls(double dp[], double q)
3423{
3424    int i,j;
3425    double scale,aa,bb,cc,delrho,bkg;
3426    int nordi=76;
3427    int nordj=76;
3428    double Pi;
3429    double sumi, sumj;
3430    double v1a, v1b, v2a, v2b;          //upper and lower integration limits
3431    double answer;
3432    double theta, phi, termAL_theta, termAT_theta, AL, AT, ahalf, bhalf, chalf, vol;
3433    double sldp, sld ;
3434
3435    Pi = 4.0*atan(1.0);
3436
3437    scale = dp[0];                //make local copies in case memory moves
3438    aa = dp[1];
3439    bb = aa * dp[2];
3440    cc = aa * dp[3];
3441    sldp = dp[4];        // scattering length density of the object
3442    sld = dp[5];         // scattering length density of the solvent
3443    delrho = sldp - sld;
3444    bkg = dp[6];
3445
3446    ahalf = 0.5 * aa;
3447    bhalf = 0.5 * bb;
3448    chalf = 0.5 * cc;
3449
3450    //Integration limits to use in Gaussian quadrature
3451    v1a = 0.;
3452    v1b = Pi/2;   //theta integration limits
3453    v2a = 0.;
3454    v2b = Pi/2.;  //phi integration limits
3455
3456    sumi = 0.0;
3457    for(i=0;i<nordi;i++) {
3458
3459        theta = ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b )/2.0;      //z values for gaussian quadrature
3460
3461        // To check potential problems if denominator goes to zero here !!!
3462        termAL_theta = 8.0*cos(q*chalf*cos(theta)) / (q*q*sin(theta)*sin(theta));
3463        termAT_theta = 8.0*sin(q*chalf*cos(theta)) / (q*q*sin(theta)*cos(theta));
3464
3465        sumj = 0.0;
3466        for(j=0;j<nordj;j++) {
3467
3468            phi = ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b )/2.0;    //z values for gaussian quadrature
3469
3470            // Amplitude AL from eqn. (7c)
3471            AL = termAL_theta * sin(q*ahalf*sin(theta)*sin(phi)) * sin(q*bhalf*sin(theta)*cos(phi)) / (sin(phi)*cos(phi));
3472
3473            // Amplitude AT from eqn. (9)
3474            AT = termAT_theta * (  (cos(q*ahalf*sin(theta)*sin(phi))*sin(q*bhalf*sin(theta)*cos(phi))/cos(phi)) + (cos(q*bhalf*sin(theta)*cos(phi))*sin(q*ahalf*sin(theta)*sin(phi))/sin(phi)) );
3475
3476            sumj += Gauss76Wt[j] * (AL+AT)*(AL+AT);
3477
3478        }
3479
3480        sumj = (v2b-v2a) * sumj / 2.0;
3481        sumi += Gauss76Wt[i] * sumj * sin(theta);
3482
3483    }
3484
3485    answer = (v1b-v1a)/2.0*sumi;
3486
3487    // Normalize as in Eqn. (11). The factor 2 is due to the different theta integration limit (pi/2 instead of pi)
3488    vol = 2.0*aa*bb+2.0*aa*cc+2.0*bb*cc;
3489    answer = 2.0 * answer / Pi / vol / vol;
3490
3491    // Multiply by contrast^2
3492    answer *= delrho*delrho;
3493
3494    //normalize by volume
3495    answer *= vol;
3496
3497    //convert to [cm-1]
3498    answer *= 1.0e8;
3499
3500    //Scale
3501    answer *= scale;
3502
3503    // add in the background
3504    answer += bkg;
3505
3506    return answer;
3507}
3508
3509
3510/*
3511 * Hollow rectangular prism (a <= b <=c)
3512 * using eqn. (13-15) from Nayuk & Huber, Z. Phys. Chem. 226, 837 (2012)
3513 */
3514 
3515double
3516RectangularHollowPrism(double dp[], double q)
3517{
3518    int i, j;
3519    double scale, aa, bb, cc, thickness, delrho, bkg;
3520    int nordi=76;
3521    int nordj=76;
3522    double Pi;
3523    double sumi, sumj;
3524    double v1a, v1b, v2a, v2b;          //upper and lower integration limits
3525    double answer;
3526    double theta, phi, arg, termA1, termB1, termC1, termA2, termB2, termC2, AP1, AP2, AP, ahalf, bhalf, chalf, vol;
3527    double sldp, sld ;
3528
3529    Pi = 4.0*atan(1.0);
3530
3531    scale = dp[0];
3532    aa = dp[1];
3533    bb = aa * dp[2];
3534    cc = aa * dp[3];
3535    thickness = dp[4];
3536    sldp = dp[5];        // scattering length density of the object
3537    sld = dp[6];         // scattering length density of the solvent
3538    delrho = sldp - sld;
3539    bkg = dp[7];
3540
3541    ahalf = 0.5 * aa;
3542    bhalf = 0.5 * bb;
3543    chalf = 0.5 * cc;
3544
3545    //Integration limits to use in Gaussian quadrature
3546    v1a = 0.;
3547    v1b = Pi/2;    //theta integration limits
3548    v2a = 0.;
3549    v2b = Pi/2.;   //phi integration limits
3550
3551    sumi = 0.0;
3552    for(i=0;i<nordi;i++) {
3553
3554        theta = ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b )/2.0;      //z values for gaussian quadrature
3555
3556            arg = q*chalf*cos(theta);
3557            if (fabs(arg) > 1.e-16) {termC1 = sin(arg) / arg;} else {termC1 = 1.0;}
3558            arg = q*(chalf-thickness)*cos(theta);
3559            if (fabs(arg) > 1.e-16) {termC2 = sin(arg) / arg;} else {termC2 = 1.0;}
3560
3561            sumj = 0.0;
3562            for(j=0;j<nordj;j++) {
3563
3564            phi = ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b )/2.0;    //z values for gaussian quadrature
3565
3566            // Amplitude AP from eqn. (13), rewritten to avoid round-off effects when arg=0
3567
3568                arg = q*ahalf*sin(theta)*sin(phi);
3569                if (fabs(arg) > 1.e-16) {termA1 = sin(arg) / arg;} else {termA1 = 1.0;}
3570                arg = q*(ahalf-thickness)*sin(theta)*sin(phi);
3571                if (fabs(arg) > 1.e-16) {termA2 = sin(arg) / arg;} else {termA2 = 1.0;}
3572
3573                arg = q*bhalf*sin(theta)*cos(phi);
3574                if (fabs(arg) > 1.e-16) {termB1 = sin(arg) / arg;} else {termB1 = 1.0;}
3575                arg = q*(bhalf-thickness)*sin(theta)*cos(phi);
3576                if (fabs(arg) > 1.e-16) {termB2 = sin(arg) / arg;} else {termB2 = 1.0;}
3577
3578            AP1 = (aa*bb*cc) * termA1 * termB1 * termC1;
3579            AP2 = 8 * (ahalf-thickness) * (bhalf-thickness) * (chalf-thickness) * termA2 * termB2 * termC2;
3580            AP = AP1 - AP2;
3581
3582            sumj += Gauss76Wt[j] * (AP*AP);
3583
3584            }
3585
3586            sumj = (v2b-v2a) * sumj / 2.0;
3587            sumi += Gauss76Wt[i] * sumj * sin(theta);
3588
3589    }
3590
3591    answer = (v1b-v1a)/2.0*sumi;
3592
3593    // Normalize as in Eqn. (15).
3594    // The factor 2 is due to the different theta integration limit (pi/2 instead of pi)
3595    vol = (aa*bb*cc-((aa-2*thickness)*(bb-2*thickness)*(cc-2*thickness)));
3596    answer = 2.0 * answer / Pi / vol / vol;
3597
3598    // Multiply by contrast^2
3599    answer *= delrho*delrho;
3600
3601    //normalize by volume
3602    answer *= vol;
3603
3604    //convert to [cm-1]
3605    answer *= 1.0e8;
3606
3607    //Scale
3608    answer *= scale;
3609
3610    // add in the background
3611    answer += bkg;
3612
3613    return answer;
3614}
3615
3616
Note: See TracBrowser for help on using the repository browser.