# source:sans/XOP_Dev/SANSAnalysis/lib/libSphere.c@453

Last change on this file since 453 was 453, checked in by srkline, 14 years ago

Additions to the library of the 2008 model functions. Direct proting of the Igor code, duplicated by these XOPs

File size: 54.9 KB
Line
1/*      SimpleFit.c
2
3A simplified project designed to act as a template for your curve fitting function.
4The fitting function is a simple polynomial. It works but is of no practical use.
5*/
6
8#include "GaussWeights.h"
9#include "libSphere.h"
10
11// scattering from a sphere - hardly needs to be an XOP...
12double
13SphereForm(double dp[], double q)
14{
15        double scale,radius,delrho,bkg,sldSph,sldSolv;          //my local names
16        double bes,f,vol,f2,pi;
17
18        pi = 4.0*atan(1.0);
19        scale = dp[0];
21        sldSph = dp[2];
22        sldSolv = dp[3];
23        bkg = dp[4];
24
25        delrho = sldSph - sldSolv;
26        //handle q==0 separately
27        if(q==0){
29                return(f);
30        }
31
34        f = vol*bes*delrho;             // [=] 
35                                                        // normalize to single particle volume, convert to 1/cm
36        f2 = f * f / vol * 1.0e8;               // [=] 1/cm
37
38        return(scale*f2+bkg);   //scale, and add in the background
39}
40
41// scattering from a monodisperse core-shell sphere - hardly needs to be an XOP...
42double
43CoreShellForm(double dp[], double q)
44{
45        double x,pi;
46        double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg;           //my local names
47        double bes,f,vol,qr,contr,f2;
48
49        pi = 4.0*atan(1.0);
50        x=q;
51
52        scale = dp[0];
53        rcore = dp[1];
54        thick = dp[2];
55        rhocore = dp[3];
56        rhoshel = dp[4];
57        rhosolv = dp[5];
58        bkg = dp[6];
59        // core first, then add in shell
60        qr=x*rcore;
61        contr = rhocore-rhoshel;
62        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
63        vol = 4.0*pi/3.0*rcore*rcore*rcore;
64        f = vol*bes*contr;
65        //now the shell
66        qr=x*(rcore+thick);
67        contr = rhoshel-rhosolv;
68        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
69        vol = 4.0*pi/3.0*pow((rcore+thick),3);
70        f += vol*bes*contr;
71
72        // normalize to particle volume and rescale from [-1] to [cm-1]
73        f2 = f*f/vol*1.0e8;
74
75        //scale if desired
76        f2 *= scale;
77        // then add in the background
78        f2 += bkg;
79
80        return(f2);
81}
82
83// scattering from a unilamellar vesicle - hardly needs to be an XOP...
84// same functional form as the core-shell sphere, but more intuitive for a vesicle
85double
86VesicleForm(double dp[], double q)
87{
88        double x,pi;
89        double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg;           //my local names
90        double bes,f,vol,qr,contr,f2;
91        pi = 4.0*atan(1.0);
92        x= q;
93
94        scale = dp[0];
95        rcore = dp[1];
96        thick = dp[2];
97        rhocore = dp[3];
98        rhosolv = rhocore;
99        rhoshel = dp[4];
100        bkg = dp[5];
101        // core first, then add in shell
102        qr=x*rcore;
103        contr = rhocore-rhoshel;
104        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
105        vol = 4.0*pi/3.0*rcore*rcore*rcore;
106        f = vol*bes*contr;
107        //now the shell
108        qr=x*(rcore+thick);
109        contr = rhoshel-rhosolv;
110        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
111        vol = 4.0*pi/3.0*pow((rcore+thick),3);
112        f += vol*bes*contr;
113
114        // normalize to the particle volume and rescale from [-1] to [cm-1]
115        //note that for the vesicle model, the volume is ONLY the shell volume
116        vol = 4.0*pi/3.0*(pow((rcore+thick),3)-pow(rcore,3));
117        f2 = f*f/vol*1.0e8;
118
119        //scale if desired
120        f2 *= scale;
121        // then add in the background
122        f2 += bkg;
123
124        return(f2);
125}
126
127
128// scattering from a core shell sphere with a (Schulz) polydisperse core and constant shell thickness
129//
130double
131PolyCoreForm(double dp[], double q)
132{
133        double pi;
134        double scale,corrad,sig,zz,del,drho1,drho2,form,bkg;            //my local names
135        double d, g ,h;
136        double qq, x, y, c1, c2, c3, c4, c5, c6, c7, c8, c9, t1, t2, t3;
137        double t4, t5, tb, cy, sy, tb1, tb2, tb3, c2y, zp1, zp2;
138        double zp3,vpoly;
139        double s2y, arg1, arg2, arg3, drh1, drh2;
140
141        pi = 4.0*atan(1.0);
142        qq= q;
143        scale = dp[0];
145        sig = dp[2];
146        del = dp[3];
147        drho1 = dp[4]-dp[5];            //core-shell
148        drho2 = dp[5]-dp[6];            //shell-solvent
149        bkg = dp[7];
150
151        zz = (1.0/sig)*(1.0/sig) - 1.0;
152
153        h=qq;
154
155        drh1 = drho1;
156        drh2 = drho2;
157        g = drh2 * -1. / drh1;
158        zp1 = zz + 1.;
159        zp2 = zz + 2.;
160        zp3 = zz + 3.;
162
163
164        // remember that h is the passed in value of q for the calculation
165        y = h *del;
167        d = atan(x * 2. / zp1);
168        arg1 = zp1 * d;
169        arg2 = zp2 * d;
170        arg3 = zp3 * d;
171        sy = sin(y);
172        cy = cos(y);
173        s2y = sin(y * 2.);
174        c2y = cos(y * 2.);
175        c1 = .5 - g * (cy + y * sy) + g * g * .5 * (y * y + 1.);
176        c2 = g * y * (g - cy);
177        c3 = (g * g + 1.) * .5 - g * cy;
178        c4 = g * g * (y * cy - sy) * (y * cy - sy) - c1;
179        c5 = g * 2. * sy * (1. - g * (y * sy + cy)) + c2;
180        c6 = c3 - g * g * sy * sy;
181        c7 = g * sy - g * .5 * g * (y * y + 1.) * s2y - c5;
182        c8 = c4 - .5 + g * cy - g * .5 * g * (y * y + 1.) * c2y;
183        c9 = g * sy * (1. - g * cy);
184
185        tb = log(zp1 * zp1 / (zp1 * zp1 + x * 4. * x));
186        tb1 = exp(zp1 * .5 * tb);
187        tb2 = exp(zp2 * .5 * tb);
188        tb3 = exp(zp3 * .5 * tb);
189
190        t1 = c1 + c2 * x + c3 * x * x * zp2 / zp1;
191        t2 = tb1 * (c4 * cos(arg1) + c7 * sin(arg1));
192        t3 = x * tb2 * (c5 * cos(arg2) + c8 * sin(arg2));
193        t4 = zp2 / zp1 * x * x * tb3 * (c6 * cos(arg3) + c9 * sin(arg3));
194        t5 = t1 + t2 + t3 + t4;
195        form = t5 * 16. * pi * pi * drh1 * drh1 / pow(qq,6);
196        //      normalize by the average volume !!! corrected for polydispersity
197        // and convert to cm-1
198        form /= vpoly;
199        form *= 1.0e8;
200        //Scale
201        form *= scale;
202        // then add in the background
203        form += bkg;
204
205        return(form);
206}
207
208// scattering from a uniform sphere with a (Schulz) size distribution
209// structure factor effects are explicitly and correctly included.
210//
211double
212PolyHardSphereIntensity(double dp[], double q)
213{
214        double pi;
215        double rad,z2,phi,cont,bkg,sigma;               //my local names
216        double mu,mu1,d1,d2,d3,d4,d5,d6,capd,rho;
217        double ll,l1,bb,cc,chi,chi1,chi2,ee,t1,t2,t3,pp;
218        double ka,zz,v1,v2,p1,p2;
219        double h1,h2,h3,h4,e1,yy,y1,s1,s2,s3,hint1,hint2;
220        double capl,capl1,capmu,capmu1,r3,pq;
221        double ka2,r1,heff;
222        double hh,k,slds,sld;
223
224        pi = 4.0*atan(1.0);
225        k= q;
226
228        z2 = dp[1];             //polydispersity (0<z2<1)
229        phi = dp[2];            // volume fraction (0<phi<1)
230        slds = dp[3];
231        sld = dp[4];
232        cont = (slds - sld)*1.0e4;              // contrast (odd units)
233        bkg = dp[5];
235
236        zz=1.0/(z2*z2)-1.0;
237        bb = sigma/(zz+1.0);
238        cc = zz+1.0;
239
240        //*c   Compute the number density by <r-cubed>, not <r> cubed*/
241        r1 = sigma/2.0;
242        r3 = r1*r1*r1;
243        r3 *= (zz+2.0)*(zz+3.0)/((zz+1.0)*(zz+1.0));
244        rho=phi/(1.3333333333*pi*r3);
245        t1 = rho*bb*cc;
246        t2 = rho*bb*bb*cc*(cc+1.0);
247        t3 = rho*bb*bb*bb*cc*(cc+1.0)*(cc+2.0);
248        capd = 1.0-pi*t3/6.0;
249        //************
250        v1=1.0/(1.0+bb*bb*k*k);
251        v2=1.0/(4.0+bb*bb*k*k);
252        pp=pow(v1,(cc/2.0))*sin(cc*atan(bb*k));
253        p1=bb*cc*pow(v1,((cc+1.0)/2.0))*sin((cc+1.0)*atan(bb*k));
254        p2=cc*(cc+1.0)*bb*bb*pow(v1,((cc+2.0)/2.0))*sin((cc+2.0)*atan(bb*k));
255        mu=pow(2,cc)*pow(v2,(cc/2.0))*sin(cc*atan(bb*k/2.0));
256        mu1=pow(2,(cc+1.0))*bb*cc*pow(v2,((cc+1.0)/2.0))*sin((cc+1.0)*atan(k*bb/2.0));
257        s1=bb*cc;
258        s2=cc*(cc+1.0)*bb*bb;
259        s3=cc*(cc+1.0)*(cc+2.0)*bb*bb*bb;
260        chi=pow(v1,(cc/2.0))*cos(cc*atan(bb*k));
261        chi1=bb*cc*pow(v1,((cc+1.0)/2.0))*cos((cc+1.0)*atan(bb*k));
262        chi2=cc*(cc+1.0)*bb*bb*pow(v1,((cc+2.0)/2.0))*cos((cc+2.0)*atan(bb*k));
263        ll=pow(2,cc)*pow(v2,(cc/2.0))*cos(cc*atan(bb*k/2.0));
264        l1=pow(2,(cc+1.0))*bb*cc*pow(v2,((cc+1.0)/2.0))*cos((cc+1.0)*atan(k*bb/2.0));
265        d1=(pi/capd)*(2.0+(pi/capd)*(t3-(rho/k)*(k*s3-p2)));
266        d2=pow((pi/capd),2)*(rho/k)*(k*s2-p1);
267        d3=(-1.0)*pow((pi/capd),2)*(rho/k)*(k*s1-pp);
268        d4=(pi/capd)*(k-(pi/capd)*(rho/k)*(chi1-s1));
269        d5=pow((pi/capd),2)*((rho/k)*(chi-1.0)+0.5*k*t2);
270        d6=pow((pi/capd),2)*(rho/k)*(chi2-s2);
271
272        e1=pow((pi/capd),2)*pow((rho/k/k),2)*((chi-1.0)*(chi2-s2)-(chi1-s1)*(chi1-s1)-(k*s1-pp)*(k*s3-p2)+pow((k*s2-p1),2));
273        ee=1.0-(2.0*pi/capd)*(1.0+0.5*pi*t3/capd)*(rho/k/k/k)*(k*s1-pp)-(2.0*pi/capd)*rho/k/k*((chi1-s1)+(0.25*pi*t2/capd)*(chi2-s2))-e1;
274        y1=pow((pi/capd),2)*pow((rho/k/k),2)*((k*s1-pp)*(chi2-s2)-2.0*(k*s2-p1)*(chi1-s1)+(k*s3-p2)*(chi-1.0));
275        yy = (2.0*pi/capd)*(1.0+0.5*pi*t3/capd)*(rho/k/k/k)*(chi+0.5*k*k*s2-1.0)-(2.0*pi*rho/capd/k/k)*(k*s2-p1+(0.25*pi*t2/capd)*(k*s3-p2))-y1;
276
277        capl=2.0*pi*cont*rho/k/k/k*(pp-0.5*k*(s1+chi1));
278        capl1=2.0*pi*cont*rho/k/k/k*(p1-0.5*k*(s2+chi2));
279        capmu=2.0*pi*cont*rho/k/k/k*(1.0-chi-0.5*k*p1);
280        capmu1=2.0*pi*cont*rho/k/k/k*(s1-chi1-0.5*k*p2);
281
282        h1=capl*(capl*(yy*d1-ee*d6)+capl1*(yy*d2-ee*d4)+capmu*(ee*d1+yy*d6)+capmu1*(ee*d2+yy*d4));
283        h2=capl1*(capl*(yy*d2-ee*d4)+capl1*(yy*d3-ee*d5)+capmu*(ee*d2+yy*d4)+capmu1*(ee*d3+yy*d5));
284        h3=capmu*(capl*(ee*d1+yy*d6)+capl1*(ee*d2+yy*d4)+capmu*(ee*d6-yy*d1)+capmu1*(ee*d4-yy*d2));
285        h4=capmu1*(capl*(ee*d2+yy*d4)+capl1*(ee*d3+yy*d5)+capmu*(ee*d4-yy*d2)+capmu1*(ee*d5-yy*d3));
286
287        //*  This part computes the second integral in equation (1) of the paper.*/
288
289        hint1 = -2.0*(h1+h2+h3+h4)/(k*k*k*(ee*ee+yy*yy));
290
291        //*  This part computes the first integral in equation (1).  It also
292        // generates the KC approximated effective structure factor.*/
293
294        pq=4.0*pi*cont*(sin(k*sigma/2.0)-0.5*k*sigma*cos(k*sigma/2.0));
295        hint2=8.0*pi*pi*rho*cont*cont/(k*k*k*k*k*k)*(1.0-chi-k*p1+0.25*k*k*(s2+chi2));
296
297        ka=k*(sigma/2.0);
298        //
299        hh=hint1+hint2;         // this is the model intensity
300                                                //
301        heff=1.0+hint1/hint2;
302        ka2=ka*ka;
303        //*
304        //  heff is PY analytical solution for intensity divided by the
305        //   form factor.  happ is the KC approximated effective S(q)
306
307         //*******************
308         //  add in the background then return the intensity value
309
310         return(hh+bkg);        //scale, and add in the background
311}
312
313// scattering from a uniform sphere with a (Schulz) size distribution, bimodal population
314// NO CROSS TERM IS ACCOUNTED FOR == DILUTE SOLUTION!!
315//
316double
317BimodalSchulzSpheres(double dp[], double q)
318{
319        double x,pq;
320        double scale,ravg,pd,bkg,rho,rhos;              //my local names
321        double scale2,ravg2,pd2,rho2;           //my local names
322
323        x= q;
324
325        scale = dp[0];
326        ravg = dp[1];
327        pd = dp[2];
328        rho = dp[3];
329        scale2 = dp[4];
330        ravg2 = dp[5];
331        pd2 = dp[6];
332        rho2 = dp[7];
333        rhos = dp[8];
334        bkg = dp[9];
335
336        pq = SchulzSphere_Fn( scale,  ravg,  pd,  rho,  rhos, x);
337        pq += SchulzSphere_Fn( scale2,  ravg2,  pd2,  rho2,  rhos, x);
338        // add in the background
339        pq += bkg;
340
341        return (pq);
342}
343
344// scattering from a uniform sphere with a (Schulz) size distribution
345//
346double
347SchulzSpheres(double dp[], double q)
348{
349        double x,pq;
350        double scale,ravg,pd,bkg,rho,rhos;              //my local names
351
352        x= q;
353
354        scale = dp[0];
355        ravg = dp[1];
356        pd = dp[2];
357        rho = dp[3];
358        rhos = dp[4];
359        bkg = dp[5];
360        pq = SchulzSphere_Fn( scale,  ravg,  pd,  rho,  rhos, x);
361        // add in the background
362        pq += bkg;
363
364        return(pq);
365}
366
367// calculates everything but the background
368double
369SchulzSphere_Fn(double scale, double ravg, double pd, double rho, double rhos, double x)
370{
371        double zp1,zp2,zp3,zp4,zp5,zp6,zp7,vpoly;
372        double aa,at1,at2,rt1,rt2,rt3,t1,t2,t3;
373        double v1,v2,v3,g1,pq,pi,delrho,zz;
374        double i_zero,Rg2,zp8;
375
376        pi = 4.0*atan(1.0);
377        delrho = rho-rhos;
378        zz = (1.0/pd)*(1.0/pd) - 1.0;
379
380        zp1 = zz + 1.0;
381        zp2 = zz + 2.0;
382        zp3 = zz + 3.0;
383        zp4 = zz + 4.0;
384        zp5 = zz + 5.0;
385        zp6 = zz + 6.0;
386        zp7 = zz + 7.0;
387        //
388
389        //small QR limit - use Guinier approx
390        zp8 = zz+8;
391        if(x*ravg < 0.1) {
392                i_zero = scale*delrho*delrho*1.e8*4.*pi/3.*pow(ravg,3);
393                i_zero *= zp6*zp5*zp4/zp1/zp1/zp1;              //6th moment / 3rd moment
394                Rg2 = 3.*zp8*zp7/5./(zp1*zp1)*ravg*ravg;
395                pq = i_zero*exp(-x*x*Rg2/3.);
396                //pq += bkg;            //unlike the Igor code, the backgorund is added in the wrapper (above)
397                return(pq);
398        }
399        //
400
401        aa = (zz+1)/x/ravg;
402
403        at1 = atan(1.0/aa);
404        at2 = atan(2.0/aa);
405        //
406        //  calculations are performed to avoid  large # errors
407        // - trick is to propogate the a^(z+7) term through the g1
408        //
409        t1 = zp7*log10(aa) - zp1/2.0*log10(aa*aa+4.0);
410        t2 = zp7*log10(aa) - zp3/2.0*log10(aa*aa+4.0);
411        t3 = zp7*log10(aa) - zp2/2.0*log10(aa*aa+4.0);
412        //      print t1,t2,t3
413        rt1 = pow(10,t1);
414        rt2 = pow(10,t2);
415        rt3 = pow(10,t3);
416        v1 = pow(aa,6) - rt1*cos(zp1*at2);
417        v2 = zp1*zp2*( pow(aa,4) + rt2*cos(zp3*at2) );
418        v3 = -2.0*zp1*rt3*sin(zp2*at2);
419        g1 = (v1+v2+v3);
420
421        pq = log10(g1) - 6.0*log10(zp1) + 6.0*log10(ravg);
422        pq = pow(10,pq)*8*pi*pi*delrho*delrho;
423
424        //
425        // beta factor is not used here, but could be for the
426        // decoupling approximation
427        //
428        //      g11 = g1
429        //      gd = -zp7*log(aa)
430        //      g1 = log(g11) + gd
431        //
432        //      t1 = zp1*at1
433        //      t2 = zp2*at1
434        //      g2 = sin( t1 ) - zp1/sqrt(aa*aa+1)*cos( t2 )
435        //      g22 = g2*g2
436        //      beta = zp1*log(aa) - zp1*log(aa*aa+1) - g1 + log(g22)
437        //      beta = 2*alog(beta)
438
439        //re-normalize by the average volume
440        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*ravg*ravg*ravg;
441        pq /= vpoly;
442        //scale, convert to cm^-1
443        pq *= scale * 1.0e8;
444
445    return(pq);
446}
447
448// scattering from a uniform sphere with a rectangular size distribution
449//
450double
451PolyRectSpheres(double dp[], double q)
452{
453        double pi,x;
454        double scale,rad,pd,cont,bkg;           //my local names
456
457        pi = 4.0*atan(1.0);
458        x= q;
459
460        scale = dp[0];
462        pd = dp[2];             //polydispersity of rectangular distribution
463        slds = dp[3];
464        sld = dp[4];
465        cont = slds - sld;              // contrast (A^-2)
466        bkg = dp[5];
467
468        // as usual, poly = sig/ravg
469        // for the rectangular distribution, sig = width/sqrt(3)
470        // width is the HALF- WIDTH of the rectangular distrubution
471
473        width = sqrt(3.0)*sig;
474
475        //x is the q-value
476        qw = x*width;
478
479        // as for the numerical inatabilities at low QR, the function is calculating the sines and cosines
480        // just fine - the problem seems to be that the
481        // leading terms nearly cancel with the last term (the -6*qr... term), to within machine
482        // precision - the difference is on the order of 10^-20
483        // so just use the limiting Guiner value
484        if(qr<0.1) {
486                h1 *=   (1. + 15.*pow(pd,2) + 27.*pow(pd,4) +27./7.*pow(pd,6) );                                //6th moment
487                h1 /= (1.+3.*pd*pd);    //3rd moment
489                Rg2 /= (1.+15.*pow(pd,2)+27.*pow(pd,4)+27./7.*pow(pd,6));
490                h1 *= exp(-1./3.*Rg2*x*x);
491                h1 += bkg;
492                return(h1);
493        }
494
495        // normal calculation
496        h1 = -0.5*qw + qr*qr*qw + (qw*qw*qw)/3.0;
497        h1 -= 5.0/2.0*cos(2*qr)*sin(qw)*cos(qw);
498        h1 += 0.5*qr*qr*cos(2*qr)*sin(2*qw);
499        h1 += 0.5*qw*qw*cos(2*qr)*sin(2*qw);
500        h1 += qw*qr*sin(2*qr)*cos(2*qw);
501        h1 += 3.0*qw*(cos(qr)*cos(qw))*(cos(qr)*cos(qw));
502        h1+= 3.0*qw*(sin(qr)*sin(qw))*(sin(qr)*sin(qw));
503        h1 -= 6.0*qr*cos(qr)*sin(qr)*cos(qw)*sin(qw);
504
505        // calculate P(q) = <f^2>
506        inten = 8.0*pi*pi*cont*cont/width/pow(x,7)*h1;
507
508        // beta(q) would be calculated as 2/width/x/h1*h2*h2
509        // with
511
512        // normalize to the average volume
513        // <R^3> = ravg^3*(1+3*pd^2)
514        // or... "zf"  = (1 + 3*p^2), which will be greater than one
515
518        //resacle to 1/cm
519        inten *= 1.0e8;
520        //scale the result
521        inten *= scale;
522        // then add in the background
523        inten += bkg;
524
525        return(inten);
526}
527
528
529// scattering from a uniform sphere with a Gaussian size distribution
530//
531double
532GaussPolySphere(double dp[], double q)
533{
534        double pi,x;
535        double scale,rad,pd,sig,rho,rhos,bkg,delrho;            //my local names
536        double va,vb,zi,yy,summ,inten;
537        int nord=20,ii;
538
539        pi = 4.0*atan(1.0);
540        x= q;
541
542        scale=dp[0];
544        pd=dp[2];
546        rho=dp[3];
547        rhos=dp[4];
548        delrho=rho-rhos;
549        bkg=dp[5];
550
551        va = -4.0*sig + rad;
552        if (va<0) {
553                va=0;           //to avoid numerical error when  va<0 (-ve q-value)
554        }
556
557        summ = 0.0;             // initialize integral
558        for(ii=0;ii<nord;ii+=1) {
559                // calculate Gauss points on integration interval (r-value for evaluation)
560                zi = ( Gauss20Z[ii]*(vb-va) + vb + va )/2.0;
561                // calculate sphere scattering
562                //return(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr)); pass qr
563                yy = F_func(x*zi)*(4.0*pi/3.0*zi*zi*zi)*delrho;
564                yy *= yy;
565                yy *= Gauss20Wt[ii] *  Gauss_distr(sig,rad,zi);
566
567                summ += yy;             //add to the running total of the quadrature
568        }
569        // calculate value of integral to return
570        inten = (vb-va)/2.0*summ;
571
572        //re-normalize by polydisperse sphere volume
574
575        inten *= 1.0e8;
576        inten *= scale;
577        inten += bkg;
578
579    return(inten);      //scale, and add in the background
580}
581
582// scattering from a uniform sphere with a LogNormal size distribution
583//
584double
585LogNormalPolySphere(double dp[], double q)
586{
587        double pi,x;
588        double scale,rad,sig,rho,rhos,bkg,delrho,mu,r3;         //my local names
589        double va,vb,zi,yy,summ,inten;
590        int nord=76,ii;
591
592        pi = 4.0*atan(1.0);
593        x= q;
594
595        scale=dp[0];
597        mu = log(dp[1]);
598        sig=dp[2];
599        rho=dp[3];
600        rhos=dp[4];
601        delrho=rho-rhos;
602        bkg=dp[5];
603
604        va = -3.5*sig + mu;
605        va = exp(va);
606        if (va<0) {
607                va=0;           //to avoid numerical error when  va<0 (-ve q-value)
608        }
609        vb = 3.5*sig*(1.0+sig) +mu;
610        vb = exp(vb);
611
612        summ = 0.0;             // initialize integral
613        for(ii=0;ii<nord;ii+=1) {
614                // calculate Gauss points on integration interval (r-value for evaluation)
615                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
616                // calculate sphere scattering
617                //return(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr)); pass qr
618                yy = F_func(x*zi)*(4.0*pi/3.0*zi*zi*zi)*delrho;
619                yy *= yy;
620                yy *= Gauss76Wt[ii] *  LogNormal_distr(sig,mu,zi);
621
622                summ += yy;             //add to the running total of the quadrature
623        }
624        // calculate value of integral to return
625        inten = (vb-va)/2.0*summ;
626
627        //re-normalize by polydisperse sphere volume
628        r3 = exp(3.0*mu + 9.0/2.0*sig*sig);             // <R^3> directly
629        inten /= (4.0*pi/3.0*r3);               //polydisperse volume
630
631        inten *= 1.0e8;
632        inten *= scale;
633        inten += bkg;
634
635        return(inten);
636}
637
638static double
639LogNormal_distr(double sig, double mu, double pt)
640{
641        double retval,pi;
642
643        pi = 4.0*atan(1.0);
644        retval = (1/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig );
645        return(retval);
646}
647
648static double
649Gauss_distr(double sig, double avg, double pt)
650{
651        double retval,Pi;
652
653        Pi = 4.0*atan(1.0);
654        retval = (1.0/ (sig*sqrt(2.0*Pi)) )*exp(-(avg-pt)*(avg-pt)/sig/sig/2.0);
655        return(retval);
656}
657
658// scattering from a core shell sphere with a (Schulz) polydisperse core and constant ratio (shell thickness)/(core radius)
659// - the polydispersity is of the WHOLE sphere
660//
661double
662PolyCoreShellRatio(double dp[], double q)
663{
664        double pi,x;
666        double sld1,sld2,sld3,zp1,zp2,zp3,vpoly;
667        double pi43,c1,c2,form,volume,arg1,arg2;
668
669        pi = 4.0*atan(1.0);
670        x= q;
671
672        scale = dp[0];
674        thick = dp[2];
675        sig = dp[3];
676        sld1 = dp[4];
677        sld2 = dp[5];
678        sld3 = dp[6];
679        bkg = dp[7];
680
681        zz = (1.0/sig)*(1.0/sig) - 1.0;
683        drho1 = sld1-sld2;              //core-shell
684        drho2 = sld2-sld3;              //shell-solvent
685        zp1 = zz + 1.;
686        zp2 = zz + 2.;
687        zp3 = zz + 3.;
689
690        // the beta factor is not calculated
691        // the calculated form factor <f^2> has units [length^2]
692        // and must be multiplied by number density [l^-3] and the correct unit
693        // conversion to get to absolute scale
694
695        pi43=4.0/3.0*pi;
698        c1=drho1*volume;
699        c2=drho2*volume;
700
703
704        form=pow(pp,6)*c1*c1*fnt2(arg1,zz);
705        form += c2*c2*fnt2(arg2,zz);
706        form += 2.0*c1*c2*fnt3(arg2,pp,zz);
707
708        //convert the result to [cm^-1]
709
710        //scale the result
711        // - divide by the polydisperse volume, mult by 10^8
712        form  /= vpoly;
713        form *= 1.0e8;
714        form *= scale;
715
717        form += bkg;
718
719        return(form);
720}
721
722//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
723//c
724//c      function fnt2(y,z)
725//c
726double
727fnt2(double yy, double zz)
728{
729        double z1,z2,z3,u,ww,term1,term2,term3,ans;
730
731        z1=zz+1.0;
732        z2=zz+2.0;
733        z3=zz+3.0;
734        u=yy/z1;
735        ww=atan(2.0*u);
736        term1=cos(z1*ww)/pow((1.0+4.0*u*u),(z1/2.0));
737        term2=2.0*yy*sin(z2*ww)/pow((1.0+4.0*u*u),(z2/2.0));
738        term3=1.0+cos(z3*ww)/pow((1.0+4.0*u*u),(z3/2.0));
739        ans=(4.50/z1/pow(yy,6))*(z1*(1.0-term1-term2)+yy*yy*z2*term3);
740
741        return(ans);
742}
743
744//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
745//c
746//c      function fnt3(y,p,z)
747//c
748double
749fnt3(double yy, double pp, double zz)
750{
751        double z1,z2,z3,yp,yn,up,un,vp,vn,term1,term2,term3,term4,term5,term6,ans;
752
753        z1=zz+1.0;
754        z2=zz+2.0;
755        z3=zz+3.0;
756        yp=(1.0+pp)*yy;
757        yn=(1.0-pp)*yy;
758        up=yp/z1;
759        un=yn/z1;
760        vp=atan(up);
761        vn=atan(un);
762        term1=cos(z1*vn)/pow((1.0+un*un),(z1/2.0));
763        term2=cos(z1*vp)/pow((1.0+up*up),(z1/2.0));
764        term3=cos(z3*vn)/pow((1.0+un*un),(z3/2.0));
765        term4=cos(z3*vp)/pow((1.0+up*up),(z3/2.0));
766        term5=yn*sin(z2*vn)/pow((1.0+un*un),(z2/2.0));
767        term6=yp*sin(z2*vp)/pow((1.0+up*up),(z2/2.0));
768        ans=4.5/z1/pow(yy,6);
769        ans *=(z1*(term1-term2)+yy*yy*pp*z2*(term3+term4)+z1*(term5-term6));
770
771        return(ans);
772}
773
774// scattering from a a binary population of hard spheres, 3 partial structure factors
775// are properly accounted for...
776//       Input (fitting) variables are:
777//      larger sphere radius(angstroms) = guess[0]
778//      smaller sphere radius (A) = w[1]
779//      number fraction of larger spheres = guess[2]
780//      total volume fraction of spheres = guess[3]
781//      size ratio, alpha(0<a<1) = derived
782//      SLD(A-2) of larger particle = guess[4]
783//      SLD(A-2) of smaller particle = guess[5]
784//      SLD(A-2) of the solvent = guess[6]
785//      background = guess[7]
786double
787BinaryHS(double dp[], double q)
788{
789        double x,pi;
790        double r2,r1,nf2,phi,aa,rho2,rho1,rhos,inten,bgd;               //my local names
791        double psf11,psf12,psf22;
792        double phi1,phi2,phr,a3;
793        double v1,v2,n1,n2,qr1,qr2,b1,b2;
794        int err;
795
796        pi = 4.0*atan(1.0);
797        x= q;
798        r2 = dp[0];
799        r1 = dp[1];
800        phi2 = dp[2];
801        phi1 = dp[3];
802        rho2 = dp[4];
803        rho1 = dp[5];
804        rhos = dp[6];
805        bgd = dp[7];
806
807
808        phi = phi1 + phi2;
809        aa = r1/r2;
810        //calculate the number fraction of larger spheres (eqn 2 in reference)
811        a3=aa*aa*aa;
812        phr=phi2/phi;
813        nf2 = phr*a3/(1.0-phr+phr*a3);
814        // calculate the PSF's here
815        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
816
817        // /* do form factor calculations  */
818
819        v1 = 4.0*pi/3.0*r1*r1*r1;
820        v2 = 4.0*pi/3.0*r2*r2*r2;
821
822        n1 = phi1/v1;
823        n2 = phi2/v2;
824
825        qr1 = r1*x;
826        qr2 = r2*x;
827
828        b1 = r1*r1*r1*(rho1-rhos)*4.0*pi*(sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1;
829        b2 = r2*r2*r2*(rho2-rhos)*4.0*pi*(sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2;
830        inten = n1*b1*b1*psf11;
831        inten += sqrt(n1*n2)*2.0*b1*b2*psf12;
832        inten += n2*b2*b2*psf22;
833        ///* convert I(1/A) to (1/cm)  */
834        inten *= 1.0e8;
835
836        inten += bgd;
837
838        return(inten);
839}
840
841double
842BinaryHS_PSF11(double dp[], double q)
843{
844        double x,pi;
845        double r2,r1,nf2,phi,aa,rho2,rho1,rhos,bgd;             //my local names
846        double psf11,psf12,psf22;
847        double phi1,phi2,phr,a3;
848        int err;
849
850        pi = 4.0*atan(1.0);
851        x= q;
852        r2 = dp[0];
853        r1 = dp[1];
854        phi2 = dp[2];
855        phi1 = dp[3];
856        rho2 = dp[4];
857        rho1 = dp[5];
858        rhos = dp[6];
859        bgd = dp[7];
860        phi = phi1 + phi2;
861        aa = r1/r2;
862        //calculate the number fraction of larger spheres (eqn 2 in reference)
863        a3=aa*aa*aa;
864        phr=phi2/phi;
865        nf2 = phr*a3/(1.0-phr+phr*a3);
866        // calculate the PSF's here
867        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
868
869    return(psf11);      //scale, and add in the background
870}
871
872double
873BinaryHS_PSF12(double dp[], double q)
874{
875        double x,pi;
876        double r2,r1,nf2,phi,aa,rho2,rho1,rhos,bgd;             //my local names
877        double psf11,psf12,psf22;
878        double phi1,phi2,phr,a3;
879        int err;
880
881        pi = 4.0*atan(1.0);
882        x= q;
883        r2 = dp[0];
884        r1 = dp[1];
885        phi2 = dp[2];
886        phi1 = dp[3];
887        rho2 = dp[4];
888        rho1 = dp[5];
889        rhos = dp[6];
890        bgd = dp[7];
891        phi = phi1 + phi2;
892        aa = r1/r2;
893        //calculate the number fraction of larger spheres (eqn 2 in reference)
894        a3=aa*aa*aa;
895        phr=phi2/phi;
896        nf2 = phr*a3/(1.0-phr+phr*a3);
897        // calculate the PSF's here
898        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
899
900    return(psf12);      //scale, and add in the background
901}
902
903double
904BinaryHS_PSF22(double dp[], double q)
905{
906        double x,pi;
907        double r2,r1,nf2,phi,aa,rho2,rho1,rhos,bgd;             //my local names
908        double psf11,psf12,psf22;
909        double phi1,phi2,phr,a3;
910        int err;
911
912        pi = 4.0*atan(1.0);
913        x= q;
914
915        r2 = dp[0];
916        r1 = dp[1];
917        phi2 = dp[2];
918        phi1 = dp[3];
919        rho2 = dp[4];
920        rho1 = dp[5];
921        rhos = dp[6];
922        bgd = dp[7];
923        phi = phi1 + phi2;
924        aa = r1/r2;
925        //calculate the number fraction of larger spheres (eqn 2 in reference)
926        a3=aa*aa*aa;
927        phr=phi2/phi;
928        nf2 = phr*a3/(1.0-phr+phr*a3);
929        // calculate the PSF's here
930        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
931
932    return(psf22);      //scale, and add in the background
933}
934
935int
936ashcroft(double qval, double r2, double nf2, double aa, double phi, double *s11, double *s22, double *s12)
937{
938        //      variable qval,r2,nf2,aa,phi,&s11,&s22,&s12
939
940        //   calculate constant terms
941        double s1,s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4;
942        double a1,a2i,a2,b1,b2,b12,gm1,gm12;
943        double err=0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3;
944        double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13;
945        double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1;
946
947        s2 = 2.0*r2;
948        s1 = aa*s2;
949        v = phi;
950        a3 = aa*aa*aa;
951        v1=((1.-nf2)*a3/(nf2+(1.-nf2)*a3))*v;
952        v2=(nf2/(nf2+(1.-nf2)*a3))*v;
953        g11=((1.+.5*v)+1.5*v2*(aa-1.))/(1.-v)/(1.-v);
954        g22=((1.+.5*v)+1.5*v1*(1./aa-1.))/(1.-v)/(1.-v);
955        g12=((1.+.5*v)+1.5*(1.-aa)*(v1-v2)/(1.+aa))/(1.-v)/(1.-v);
956        wmv = 1/(1.-v);
957        wmv3 = wmv*wmv*wmv;
958        wmv4 = wmv*wmv3;
959        a1=3.*wmv4*((v1+a3*v2)*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2))) + ((v1+a3*v2)*(1.+2.*v)+(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)-3.*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*wmv3;
960        a2i=((v1+a3*v2)*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*3*wmv4 + ((v1+a3*v2)*(1.+2.*v)+a3*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*aa-3.*v1*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*wmv3;
961        a2=a2i/a3;
962        b1=-6.*(v1*g11*g11+.25*v2*(1.+aa)*(1.+aa)*aa*g12*g12);
963        b2=-6.*(v2*g22*g22+.25*v1/a3*(1.+aa)*(1.+aa)*g12*g12);
964        b12=-3.*aa*(1.+aa)*(v1*g11/aa/aa+v2*g22)*g12;
965        gm1=(v1*a1+a3*v2*a2)*.5;
966        gm12=2.*gm1*(1.-aa)/aa;
967        //c
968        //c   calculate the direct correlation functions and print results
969        //c
970        //      do 20 j=1,npts
971
972        yy=qval*s2;
973        //c   calculate direct correlation functions
974        //c   ----c11
975        ay=aa*yy;
976        ay2 = ay*ay;
977        ay3 = ay*ay*ay;
978        t1=a1*(sin(ay)-ay*cos(ay));
979        t2=b1*(2.*ay*sin(ay)-(ay2-2.)*cos(ay)-2.)/ay;
980        t3=gm1*((4.*ay*ay2-24.*ay)*sin(ay)-(ay2*ay2-12.*ay2+24.)*cos(ay)+24.)/ay3;
981        f11=24.*v1*(t1+t2+t3)/ay3;
982
983        //c ------c22
984        y2=yy*yy;
985        y3=yy*y2;
986        tt1=a2*(sin(yy)-yy*cos(yy));
987        tt2=b2*(2.*yy*sin(yy)-(y2-2.)*cos(yy)-2.)/yy;
988        tt3=gm1*((4.*y3-24.*yy)*sin(yy)-(y2*y2-12.*y2+24.)*cos(yy)+24.)/ay3;
989        f22=24.*v2*(tt1+tt2+tt3)/y3;
990
991        //c   -----c12
992        yl=.5*yy*(1.-aa);
993        yl3=yl*yl*yl;
994        wma3 = (1.-aa)*(1.-aa)*(1.-aa);
995        y1=aa*yy;
996        y13 = y1*y1*y1;
997        ttt1=3.*wma3*v*sqrt(nf2)*sqrt(1.-nf2)*a1*(sin(yl)-yl*cos(yl))/((nf2+(1.-nf2)*a3)*yl3);
998        t21=b12*(2.*y1*cos(y1)+(y1*y1-2.)*sin(y1));
999        t22=gm12*((3.*y1*y1-6.)*cos(y1)+(y1*y1*y1-6.*y1)*sin(y1)+6.)/y1;
1000        t23=gm1*((4.*y13-24.*y1)*cos(y1)+(y13*y1-12.*y1*y1+24.)*sin(y1))/(y1*y1);
1001        t31=b12*(2.*y1*sin(y1)-(y1*y1-2.)*cos(y1)-2.);
1002        t32=gm12*((3.*y1*y1-6.)*sin(y1)-(y1*y1*y1-6.*y1)*cos(y1))/y1;
1003        t33=gm1*((4.*y13-24.*y1)*sin(y1)-(y13*y1-12.*y1*y1+24.)*cos(y1)+24.)/(y1*y1);
1004        t41=cos(yl)*((sin(y1)-y1*cos(y1))/(y1*y1) + (1.-aa)/(2.*aa)*(1.-cos(y1))/y1);
1005        t42=sin(yl)*((cos(y1)+y1*sin(y1)-1.)/(y1*y1) + (1.-aa)/(2.*aa)*sin(y1)/y1);
1006        ttt2=sin(yl)*(t21+t22+t23)/(y13*y1);
1007        ttt3=cos(yl)*(t31+t32+t33)/(y13*y1);
1008        ttt4=a1*(t41+t42)/y1;
1009        f12=ttt1+24.*v*sqrt(nf2)*sqrt(1.-nf2)*a3*(ttt2+ttt3+ttt4)/(nf2+(1.-nf2)*a3);
1010
1011        c11=f11;
1012        c22=f22;
1013        c12=f12;
1014        *s11=1./(1.+c11-(c12)*c12/(1.+c22));
1015        *s22=1./(1.+c22-(c12)*c12/(1.+c11));
1016        *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));
1017
1018        return(err);
1019}
1020
1021
1022
1023/*
1024 // calculates the scattering from a spherical particle made up of a core (aqueous) surrounded
1025 // by N spherical layers, each of which is a PAIR of shells, solvent + surfactant since there
1026 //must always be a surfactant layer on the outside
1027 //
1028 // bragg peaks arise naturally from the periodicity of the sample
1029 // resolution smeared version gives he most appropriate view of the model
1030
1031        Warning:
1032 The call to WaveData() below returns a pointer to the middle
1033 of an unlocked Macintosh handle. In the unlikely event that your
1034 calculations could cause memory to move, you should copy the coefficient
1035 values to local variables or an array before such operations.
1036 */
1037double
1038MultiShell(double dp[], double q)
1039{
1040        double x;
1041        double scale,rcore,tw,ts,rhocore,rhoshel,num,bkg;               //my local names
1042        int ii;
1043        double fval,voli,ri,sldi;
1044        double pi;
1045
1046        pi = 4.0*atan(1.0);
1047
1048        x= q;
1049        scale = dp[0];
1050        rcore = dp[1];
1051        ts = dp[2];
1052        tw = dp[3];
1053        rhocore = dp[4];
1054        rhoshel = dp[5];
1055        num = dp[6];
1056        bkg = dp[7];
1057
1058        //calculate with a loop, two shells at a time
1059
1060        ii=0;
1061        fval=0;
1062
1063        do {
1064                ri = rcore + (double)ii*ts + (double)ii*tw;
1065                voli = 4*pi/3*ri*ri*ri;
1066                sldi = rhocore-rhoshel;
1067                fval += voli*sldi*F_func(ri*x);
1068                ri += ts;
1069                voli = 4*pi/3*ri*ri*ri;
1070                sldi = rhoshel-rhocore;
1071                fval += voli*sldi*F_func(ri*x);
1072                ii+=1;          //do 2 layers at a time
1073        } while(ii<=num-1);  //change to make 0 < num < 2 correspond to unilamellar vesicles (C. Glinka, 11/24/03)
1074
1075        fval *= fval;           //square it
1076        fval /= voli;           //normalize by the overall volume
1077        fval *= scale*1e8;
1078        fval += bkg;
1079
1080        return(fval);
1081}
1082
1083/*
1084 // calculates the scattering from a POLYDISPERSE spherical particle made up of a core (aqueous) surrounded
1085 // by N spherical layers, each of which is a PAIR of shells, solvent + surfactant since there
1086 //must always be a surfactant layer on the outside
1087 //
1088 // bragg peaks arise naturally from the periodicity of the sample
1089 // resolution smeared version gives he most appropriate view of the model
1090 //
1091 // Polydispersity is of the total (outer) radius. This is converted into a distribution of MLV's
1092 // with integer numbers of layers, with a minimum of one layer... a vesicle... depending
1093 // on the parameters, the "distribution" of MLV's that is used may be truncated
1094 //
1095        Warning:
1096 The call to WaveData() below returns a pointer to the middle
1097 of an unlocked Macintosh handle. In the unlikely event that your
1098 calculations could cause memory to move, you should copy the coefficient
1099 values to local variables or an array before such operations.
1100 */
1101double
1102PolyMultiShell(double dp[], double q)
1103{
1104        double x;
1105        double scale,rcore,tw,ts,rhocore,rhoshel,bkg;           //my local names
1106        int ii,minPairs,maxPairs,first;
1107        double fval,ri,pi;
1108        double avg,pd,zz,lo,hi,r1,r2,d1,d2,distr;
1109
1110        pi = 4.0*atan(1.0);
1111        x= q;
1112
1113        scale = dp[0];
1114        avg = dp[1];            // average (total) outer radius
1115        pd = dp[2];
1116        rcore = dp[3];
1117        ts = dp[4];
1118        tw = dp[5];
1119        rhocore = dp[6];
1120        rhoshel = dp[7];
1121        bkg = dp[8];
1122
1123        zz = (1.0/pd)*(1.0/pd)-1.0;
1124
1125        //max radius set to be 5 std deviations past mean
1126        hi = avg + pd*avg*5.0;
1127        lo = avg - pd*avg*5.0;
1128
1129        maxPairs = trunc( (hi-rcore+tw)/(ts+tw) );
1130        minPairs = trunc( (lo-rcore+tw)/(ts+tw) );
1131        minPairs = (minPairs < 1) ? 1 : minPairs;       // need a minimum of one
1132
1133        ii=minPairs;
1134        fval=0;
1135        d1 = 0;
1136        d2 = 0;
1137        r1 = 0;
1138        r2 = 0;
1139        distr = 0;
1140        first = 1;
1141        do {
1142                //make the current values old
1143                r1 = r2;
1144                d1 = d2;
1145
1146                ri = (double)ii*(ts+tw) - tw + rcore;
1147                fval += SchulzPoint(ri,avg,zz) * MultiShellGuts(x, rcore, ts, tw, rhocore, rhoshel, ii) * (4*pi/3*ri*ri*ri);
1148                // get a running integration of the fraction of the distribution used, but not the first time
1149                r2 = ri;
1150                d2 = SchulzPoint(ri,avg,zz);
1151                if( !first ) {
1152                        distr += 0.5*(d1+d2)*(r2-r1);           //cheap trapezoidal integration
1153                }
1154                ii+=1;
1155                first = 0;
1156        } while(ii<=maxPairs);
1157
1158        fval /= 4*pi/3*avg*avg*avg;             //normalize by the overall volume
1159        fval /= distr;
1160        fval *= scale;
1161        fval += bkg;
1162
1163        return(fval);
1164}
1165
1166double
1167MultiShellGuts(double x,double rcore,double ts,double tw,double rhocore,double rhoshel,int num) {
1168
1169    double ri,sldi,fval,voli,pi;
1170    int ii;
1171
1172        pi = 4.0*atan(1.0);
1173    ii=0;
1174    fval=0;
1175
1176    do {
1177        ri = rcore + (double)ii*ts + (double)ii*tw;
1178        voli = 4*pi/3*ri*ri*ri;
1179        sldi = rhocore-rhoshel;
1180        fval += voli*sldi*F_func(ri*x);
1181        ri += ts;
1182        voli = 4*pi/3*ri*ri*ri;
1183        sldi = rhoshel-rhocore;
1184        fval += voli*sldi*F_func(ri*x);
1185        ii+=1;          //do 2 layers at a time
1186    } while(ii<=num-1);  //change to make 0 < num < 2 correspond to unilamellar vesicles (C. Glinka, 11/24/03)
1187
1188    fval *= fval;
1189    fval /= voli;
1190    fval *= 1e8;
1191
1192    return(fval);       // this result still needs to be multiplied by scale and have background added
1193}
1194
1195static double
1196SchulzPoint(double x, double avg, double zz) {
1197
1198    double dr;
1199
1200    dr = zz*log(x) - gammln(zz+1)+(zz+1)*log((zz+1)/avg)-(x/avg*(zz+1));
1201    return (exp(dr));
1202}
1203
1204static double
1205gammln(double xx) {
1206
1207    double x,y,tmp,ser;
1208    static double cof[6]={76.18009172947146,-86.50532032941677,
1209                24.01409824083091,-1.231739572450155,
1210                0.1208650973866179e-2,-0.5395239384953e-5};
1211    int j;
1212
1213    y=x=xx;
1214    tmp=x+5.5;
1215    tmp -= (x+0.5)*log(tmp);
1216    ser=1.000000000190015;
1217    for (j=0;j<=5;j++) ser += cof[j]/++y;
1218    return -tmp+log(2.5066282746310005*ser/x);
1219}
1220
1221double
1222F_func(double qr) {
1223        return(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr));
1224}
1225
1226double
1227OneShell(double dp[], double q)
1228{
1229        // variables are:
1230        //[0] scale factor
1231        //[1] radius of core []
1232        //[2] SLD of the core   [-2]
1233        //[3] thickness of the shell    []
1234        //[4] SLD of the shell
1235        //[5] SLD of the solvent
1236        //[6] background        [cm-1]
1237
1238        double x,pi;
1239        double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg;           //my local names
1240        double bes,f,vol,qr,contr,f2;
1241
1242        pi = 4.0*atan(1.0);
1243        x=q;
1244
1245        scale = dp[0];
1246        rcore = dp[1];
1247        rhocore = dp[2];
1248        thick = dp[3];
1249        rhoshel = dp[4];
1250        rhosolv = dp[5];
1251        bkg = dp[6];
1252
1253        // core first, then add in shell
1254        qr=x*rcore;
1255        contr = rhocore-rhoshel;
1256        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1257        vol = 4.0*pi/3.0*rcore*rcore*rcore;
1258        f = vol*bes*contr;
1259        //now the shell
1260        qr=x*(rcore+thick);
1261        contr = rhoshel-rhosolv;
1262        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1263        vol = 4.0*pi/3.0*pow((rcore+thick),3);
1264        f += vol*bes*contr;
1265
1266        // normalize to particle volume and rescale from [-1] to [cm-1]
1267        f2 = f*f/vol*1.0e8;
1268
1269        //scale if desired
1270        f2 *= scale;
1271        // then add in the background
1272        f2 += bkg;
1273
1274        return(f2);
1275}
1276
1277double
1278TwoShell(double dp[], double q)
1279{
1280        // variables are:
1281        //[0] scale factor
1282        //[1] radius of core []
1283        //[2] SLD of the core   [-2]
1284        //[3] thickness of shell 1 []
1285        //[4] SLD of shell 1
1286        //[5] thickness of shell 2 []
1287        //[6] SLD of shell 2
1288        //[7] SLD of the solvent
1289        //[8] background        [cm-1]
1290
1291        double x,pi;
1292        double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;         //my local names
1293        double bes,f,vol,qr,contr,f2;
1294        double rhoshel2,thick2;
1295
1296        pi = 4.0*atan(1.0);
1297        x=q;
1298
1299        scale = dp[0];
1300        rcore = dp[1];
1301        rhocore = dp[2];
1302        thick1 = dp[3];
1303        rhoshel1 = dp[4];
1304        thick2 = dp[5];
1305        rhoshel2 = dp[6];
1306        rhosolv = dp[7];
1307        bkg = dp[8];
1308                // core first, then add in shells
1309        qr=x*rcore;
1310        contr = rhocore-rhoshel1;
1311        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1312        vol = 4.0*pi/3*rcore*rcore*rcore;
1313        f = vol*bes*contr;
1314        //now the shell (1)
1315        qr=x*(rcore+thick1);
1316        contr = rhoshel1-rhoshel2;
1317        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1318        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1);
1319        f += vol*bes*contr;
1320        //now the shell (2)
1321        qr=x*(rcore+thick1+thick2);
1322        contr = rhoshel2-rhosolv;
1323        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1324        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2);
1325        f += vol*bes*contr;
1326
1327
1328        // normalize to particle volume and rescale from [-1] to [cm-1]
1329        f2 = f*f/vol*1.0e8;
1330
1331        //scale if desired
1332        f2 *= scale;
1333        // then add in the background
1334        f2 += bkg;
1335
1336        return(f2);
1337}
1338
1339double
1340ThreeShell(double dp[], double q)
1341{
1342        // variables are:
1343        //[0] scale factor
1344        //[1] radius of core []
1345        //[2] SLD of the core   [-2]
1346        //[3] thickness of shell 1 []
1347        //[4] SLD of shell 1
1348        //[5] thickness of shell 2 []
1349        //[6] SLD of shell 2
1350        //[7] thickness of shell 3
1351        //[8] SLD of shell 3
1352        //[9] SLD of solvent
1353        //[10] background       [cm-1]
1354
1355        double x,pi;
1356        double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;         //my local names
1357        double bes,f,vol,qr,contr,f2;
1358        double rhoshel2,thick2,rhoshel3,thick3;
1359
1360        pi = 4.0*atan(1.0);
1361        x=q;
1362
1363        scale = dp[0];
1364        rcore = dp[1];
1365        rhocore = dp[2];
1366        thick1 = dp[3];
1367        rhoshel1 = dp[4];
1368        thick2 = dp[5];
1369        rhoshel2 = dp[6];
1370        thick3 = dp[7];
1371        rhoshel3 = dp[8];
1372        rhosolv = dp[9];
1373        bkg = dp[10];
1374
1375                // core first, then add in shells
1376        qr=x*rcore;
1377        contr = rhocore-rhoshel1;
1378        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1379        vol = 4.0*pi/3*rcore*rcore*rcore;
1380        f = vol*bes*contr;
1381        //now the shell (1)
1382        qr=x*(rcore+thick1);
1383        contr = rhoshel1-rhoshel2;
1384        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1385        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1);
1386        f += vol*bes*contr;
1387        //now the shell (2)
1388        qr=x*(rcore+thick1+thick2);
1389        contr = rhoshel2-rhoshel3;
1390        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1391        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2);
1392        f += vol*bes*contr;
1393        //now the shell (3)
1394        qr=x*(rcore+thick1+thick2+thick3);
1395        contr = rhoshel3-rhosolv;
1396        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1397        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3);
1398        f += vol*bes*contr;
1399
1400        // normalize to particle volume and rescale from [-1] to [cm-1]
1401        f2 = f*f/vol*1.0e8;
1402
1403        //scale if desired
1404        f2 *= scale;
1405        // then add in the background
1406        f2 += bkg;
1407
1408        return(f2);
1409}
1410
1411double
1412FourShell(double dp[], double q)
1413{
1414        // variables are:
1415        //[0] scale factor
1416        //[1] radius of core []
1417        //[2] SLD of the core   [-2]
1418        //[3] thickness of shell 1 []
1419        //[4] SLD of shell 1
1420        //[5] thickness of shell 2 []
1421        //[6] SLD of shell 2
1422        //[7] thickness of shell 3
1423        //[8] SLD of shell 3
1424        //[9] thickness of shell 3
1425        //[10] SLD of shell 3
1426        //[11] SLD of solvent
1427        //[12] background       [cm-1]
1428
1429        double x,pi;
1430        double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;         //my local names
1431        double bes,f,vol,qr,contr,f2;
1432        double rhoshel2,thick2,rhoshel3,thick3,rhoshel4,thick4;
1433
1434        pi = 4.0*atan(1.0);
1435        x=q;
1436
1437        scale = dp[0];
1438        rcore = dp[1];
1439        rhocore = dp[2];
1440        thick1 = dp[3];
1441        rhoshel1 = dp[4];
1442        thick2 = dp[5];
1443        rhoshel2 = dp[6];
1444        thick3 = dp[7];
1445        rhoshel3 = dp[8];
1446        thick4 = dp[9];
1447        rhoshel4 = dp[10];
1448        rhosolv = dp[11];
1449        bkg = dp[12];
1450
1451                // core first, then add in shells
1452        qr=x*rcore;
1453        contr = rhocore-rhoshel1;
1454        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1455        vol = 4.0*pi/3*rcore*rcore*rcore;
1456        f = vol*bes*contr;
1457        //now the shell (1)
1458        qr=x*(rcore+thick1);
1459        contr = rhoshel1-rhoshel2;
1460        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1461        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1);
1462        f += vol*bes*contr;
1463        //now the shell (2)
1464        qr=x*(rcore+thick1+thick2);
1465        contr = rhoshel2-rhoshel3;
1466        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1467        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2);
1468        f += vol*bes*contr;
1469        //now the shell (3)
1470        qr=x*(rcore+thick1+thick2+thick3);
1471        contr = rhoshel3-rhoshel4;
1472        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1473        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3);
1474        f += vol*bes*contr;
1475        //now the shell (4)
1476        qr=x*(rcore+thick1+thick2+thick3+thick4);
1477        contr = rhoshel4-rhosolv;
1478        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1479        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4);
1480        f += vol*bes*contr;
1481
1482
1483        // normalize to particle volume and rescale from [-1] to [cm-1]
1484        f2 = f*f/vol*1.0e8;
1485
1486        //scale if desired
1487        f2 *= scale;
1488        // then add in the background
1489        f2 += bkg;
1490
1491        return(f2);
1492}
1493
1494double
1495PolyOneShell(double dp[], double x)
1496{
1497        double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg,pd,zz;             //my local names
1498        double va,vb,summ,yyy,zi;
1500        int nord=76,ii;
1501
1502        pi = 4.0*atan(1.0);
1503
1504        scale = dp[0];
1505        rcore = dp[1];
1506        pd = dp[2];
1507        rhocore = dp[3];
1508        thick = dp[4];
1509        rhoshel = dp[5];
1510        rhosolv = dp[6];
1511        bkg = dp[7];
1512
1513        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only
1514
1515        range = 8.0;                    //std deviations for the integration
1516        va = rcore*(1.0-range*pd);
1517        if (va<0) {
1518                va=0;           //otherwise numerical error when pd >= 0.3, making a<0
1519        }
1520        if (pd>0.3) {
1521                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail
1522        }
1523        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius?
1524
1525        //temp set scale=1 and bkg=0 for quadrature calc
1526        temp_1sf[0] = 1.0;
1527        temp_1sf[1] = dp[1];    //the core radius will be changed in the loop
1528        temp_1sf[2] = dp[3];
1529        temp_1sf[3] = dp[4];
1530        temp_1sf[4] = dp[5];
1531        temp_1sf[5] = dp[6];
1532        temp_1sf[6] = 0.0;
1533
1534        summ = 0.0;             // initialize integral
1535        for(ii=0;ii<nord;ii+=1) {
1536                // calculate Gauss points on integration interval (r-value for evaluation)
1537                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
1538                temp_1sf[1] = zi;
1539                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * OneShell(temp_1sf,x);
1540                //un-normalize by volume
1541                yyy *= 4.0*pi/3.0*pow((zi+thick),3);
1542                summ += yyy;            //add to the running total of the quadrature
1543        }
1544        // calculate value of integral to return
1546
1547        //re-normalize by the average volume
1548        zp1 = zz + 1.0;
1549        zp2 = zz + 2.0;
1550        zp3 = zz + 3.0;
1551        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick),3);
1553//scale
1557
1559}
1560
1561double
1562PolyTwoShell(double dp[], double x)
1563{
1564        double scale,rcore,rhocore,rhosolv,bkg,pd,zz;           //my local names
1565        double va,vb,summ,yyy,zi;
1567        int nord=76,ii;
1568        double thick1,thick2;
1569        double rhoshel1,rhoshel2;
1570
1571        scale = dp[0];
1572        rcore = dp[1];
1573        pd = dp[2];
1574        rhocore = dp[3];
1575        thick1 = dp[4];
1576        rhoshel1 = dp[5];
1577        thick2 = dp[6];
1578        rhoshel2 = dp[7];
1579        rhosolv = dp[8];
1580        bkg = dp[9];
1581
1582        pi = 4.0*atan(1.0);
1583
1584        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only
1585
1586        range = 8.0;                    //std deviations for the integration
1587        va = rcore*(1.0-range*pd);
1588        if (va<0) {
1589                va=0;           //otherwise numerical error when pd >= 0.3, making a<0
1590        }
1591        if (pd>0.3) {
1592                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail
1593        }
1594        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius?
1595
1596        //temp set scale=1 and bkg=0 for quadrature calc
1597        temp_2sf[0] = 1.0;
1598        temp_2sf[1] = dp[1];            //the core radius will be changed in the loop
1599        temp_2sf[2] = dp[3];
1600        temp_2sf[3] = dp[4];
1601        temp_2sf[4] = dp[5];
1602        temp_2sf[5] = dp[6];
1603        temp_2sf[6] = dp[7];
1604        temp_2sf[7] = dp[8];
1605        temp_2sf[8] = 0.0;
1606
1607        summ = 0.0;             // initialize integral
1608        for(ii=0;ii<nord;ii+=1) {
1609                // calculate Gauss points on integration interval (r-value for evaluation)
1610                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
1611                temp_2sf[1] = zi;
1612                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * TwoShell(temp_2sf,x);
1613                //un-normalize by volume
1614                yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2),3);
1615                summ += yyy;            //add to the running total of the quadrature
1616        }
1617        // calculate value of integral to return
1619
1620        //re-normalize by the average volume
1621        zp1 = zz + 1.0;
1622        zp2 = zz + 2.0;
1623        zp3 = zz + 3.0;
1624        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2),3);
1626//scale
1630
1632}
1633
1634double
1635PolyThreeShell(double dp[], double x)
1636{
1637        double scale,rcore,rhocore,rhosolv,bkg,pd,zz;           //my local names
1638        double va,vb,summ,yyy,zi;
1640        int nord=76,ii;
1641        double thick1,thick2,thick3;
1642        double rhoshel1,rhoshel2,rhoshel3;
1643
1644        scale = dp[0];
1645        rcore = dp[1];
1646        pd = dp[2];
1647        rhocore = dp[3];
1648        thick1 = dp[4];
1649        rhoshel1 = dp[5];
1650        thick2 = dp[6];
1651        rhoshel2 = dp[7];
1652        thick3 = dp[8];
1653        rhoshel3 = dp[9];
1654        rhosolv = dp[10];
1655        bkg = dp[11];
1656
1657        pi = 4.0*atan(1.0);
1658
1659        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only
1660
1661        range = 8.0;                    //std deviations for the integration
1662        va = rcore*(1.0-range*pd);
1663        if (va<0) {
1664                va=0;           //otherwise numerical error when pd >= 0.3, making a<0
1665        }
1666        if (pd>0.3) {
1667                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail
1668        }
1669        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius?
1670
1671        //temp set scale=1 and bkg=0 for quadrature calc
1672        temp_3sf[0] = 1.0;
1673        temp_3sf[1] = dp[1];            //the core radius will be changed in the loop
1674        temp_3sf[2] = dp[3];
1675        temp_3sf[3] = dp[4];
1676        temp_3sf[4] = dp[5];
1677        temp_3sf[5] = dp[6];
1678        temp_3sf[6] = dp[7];
1679        temp_3sf[7] = dp[8];
1680        temp_3sf[8] = dp[9];
1681        temp_3sf[9] = dp[10];
1682        temp_3sf[10] = 0.0;
1683
1684        summ = 0.0;             // initialize integral
1685        for(ii=0;ii<nord;ii+=1) {
1686                // calculate Gauss points on integration interval (r-value for evaluation)
1687                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
1688                temp_3sf[1] = zi;
1689                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * ThreeShell(temp_3sf,x);
1690                //un-normalize by volume
1691                yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2+thick3),3);
1692                summ += yyy;            //add to the running total of the quadrature
1693        }
1694        // calculate value of integral to return
1696
1697        //re-normalize by the average volume
1698        zp1 = zz + 1.0;
1699        zp2 = zz + 2.0;
1700        zp3 = zz + 3.0;
1701        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2+thick3),3);
1703//scale
1707
1709}
1710
1711double
1712PolyFourShell(double dp[], double x)
1713{
1714        double scale,rcore,rhocore,rhosolv,bkg,pd,zz;           //my local names
1715        double va,vb,summ,yyy,zi;
1717        int nord=76,ii;
1718        double thick1,thick2,thick3,thick4;
1719        double rhoshel1,rhoshel2,rhoshel3,rhoshel4;
1720
1721        scale = dp[0];
1722        rcore = dp[1];
1723        pd = dp[2];
1724        rhocore = dp[3];
1725        thick1 = dp[4];
1726        rhoshel1 = dp[5];
1727        thick2 = dp[6];
1728        rhoshel2 = dp[7];
1729        thick3 = dp[8];
1730        rhoshel3 = dp[9];
1731        thick4 = dp[10];
1732        rhoshel4 = dp[11];
1733        rhosolv = dp[12];
1734        bkg = dp[13];
1735
1736        pi = 4.0*atan(1.0);
1737
1738        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only
1739
1740        range = 8.0;                    //std deviations for the integration
1741        va = rcore*(1.0-range*pd);
1742        if (va<0) {
1743                va=0;           //otherwise numerical error when pd >= 0.3, making a<0
1744        }
1745        if (pd>0.3) {
1746                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail
1747        }
1748        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius?
1749
1750        //temp set scale=1 and bkg=0 for quadrature calc
1751        temp_4sf[0] = 1.0;
1752        temp_4sf[1] = dp[1];            //the core radius will be changed in the loop
1753        temp_4sf[2] = dp[3];
1754        temp_4sf[3] = dp[4];
1755        temp_4sf[4] = dp[5];
1756        temp_4sf[5] = dp[6];
1757        temp_4sf[6] = dp[7];
1758        temp_4sf[7] = dp[8];
1759        temp_4sf[8] = dp[9];
1760        temp_4sf[9] = dp[10];
1761        temp_4sf[10] = dp[11];
1762        temp_4sf[11] = dp[12];
1763        temp_4sf[12] = 0.0;
1764
1765        summ = 0.0;             // initialize integral
1766        for(ii=0;ii<nord;ii+=1) {
1767                // calculate Gauss points on integration interval (r-value for evaluation)
1768                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
1769                temp_4sf[1] = zi;
1770                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * FourShell(temp_4sf,x);
1771                //un-normalize by volume
1772                yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2+thick3+thick4),3);
1773                summ += yyy;            //add to the running total of the quadrature
1774        }
1775        // calculate value of integral to return
1777
1778        //re-normalize by the average volume
1779        zp1 = zz + 1.0;
1780        zp2 = zz + 2.0;
1781        zp3 = zz + 3.0;
1782        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2+thick3+thick4),3);
1784//scale
1788
1790}
1791
1792
1793/*      BCC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x
1794
1795Uses 150 pt Gaussian quadrature for both integrals
1796
1797*/
1798double
1799BCC_ParaCrystal(double w[], double x)
1800{
1801        int i,j;
1802        double Pi;
1803        double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;           //local variables of coefficient wave
1804        int nordi=150;                  //order of integration
1805        int nordj=150;
1806        double va,vb;           //upper and lower integration limits
1807        double summ,zi,yyy,answer;                      //running tally of integration
1808        double summj,vaj,vbj,zij;                       //for the inner integration
1809
1810        Pi = 4.0*atan(1.0);
1811        va = 0.0;
1812        vb = 2.0*Pi;            //orintational average, outer integral
1813        vaj = 0.0;
1814        vbj = Pi;               //endpoints of inner integral
1815
1816        summ = 0.0;                     //initialize intergral
1817
1818        scale = w[0];
1819        Dnn = w[1];                                     //Nearest neighbor distance A
1820        gg = w[2];                                      //Paracrystal distortion factor
1822        sld = w[4];
1823        sldSolv = w[5];
1824        background = w[6];
1825
1826        contrast = sld - sldSolv;
1827
1828        //Volume fraction calculated from lattice symmetry and sphere radius
1830
1831        for(i=0;i<nordi;i++) {
1832                //setup inner integral over the ellipsoidal cross-section
1833                summj=0.0;
1834                zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;            //the outer dummy is phi
1835                for(j=0;j<nordj;j++) {
1836                        //20 gauss points for the inner integral
1837                        zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;               //the inner dummy is theta
1838                        yyy = Gauss150Wt[j] * BCC_Integrand(w,x,zi,zij);
1839                        summj += yyy;
1840                }
1841                //now calculate the value of the inner integral
1843
1844                //now calculate outer integral
1845                yyy = Gauss150Wt[i] * answer;
1846                summ += yyy;
1847        }               //final scaling is done at the end of the function, after the NT_FP64 case
1848
1850        // Multiply by contrast^2
1852        // add in the background
1854
1856}
1857
1858// xx is phi (outer)
1859// yy is theta (inner)
1860double
1861BCC_Integrand(double w[], double qq, double xx, double yy) {
1862
1863        double retVal,temp1,temp3,aa,Da,Dnn,gg,Pi;
1864
1865        Dnn = w[1]; //Nearest neighbor distance A
1866        gg = w[2]; //Paracrystal distortion factor
1867        aa = Dnn;
1868        Da = gg*aa;
1869
1870        Pi = 4.0*atan(1.0);
1871        temp1 = qq*qq*Da*Da;
1872        temp3 = qq*aa;
1873
1874        retVal = BCCeval(yy,xx,temp1,temp3);
1875        retVal /=4.0*Pi;
1876
1877        return(retVal);
1878}
1879
1880double
1881BCCeval(double Theta, double Phi, double temp1, double temp3) {
1882
1883        double temp6,temp7,temp8,temp9,temp10;
1884        double result;
1885
1886        temp6 = sin(Theta);
1887        temp7 = sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)+cos(Theta);
1888        temp8 = -1.0*sin(Theta)*cos(Phi)-sin(Theta)*sin(Phi)+cos(Theta);
1889        temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)-cos(Theta);
1890        temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
1891        result = pow(1.0-(temp10*temp10),3)*temp6/((1.0-2.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp9))+(temp10*temp10)));
1892
1893        return (result);
1894}
1895
1896double
1897SphereForm_Paracrystal(double radius, double delrho, double x) {
1898
1899        double bes,f,vol,f2,pi;
1900        pi = 4.0*atan(1.0);
1901        //
1902        //handle q==0 separately
1903        if(x==0) {
1905                return(f);
1906        }
1907
1910        f = vol*bes*delrho      ;       // [=] 
1911        // normalize to single particle volume, convert to 1/cm
1912        f2 = f * f / vol * 1.0e8;               // [=] 1/cm
1913
1914        return (f2);
1915}
1916
1917/*      FCC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x
1918
1919Uses 150 pt Gaussian quadrature for both integrals
1920
1921*/
1922double
1923FCC_ParaCrystal(double w[], double x)
1924{
1925        int i,j;
1926        double Pi;
1927        double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;           //local variables of coefficient wave
1928        int nordi=150;                  //order of integration
1929        int nordj=150;
1930        double va,vb;           //upper and lower integration limits
1931        double summ,zi,yyy,answer;                      //running tally of integration
1932        double summj,vaj,vbj,zij;                       //for the inner integration
1933
1934        Pi = 4.0*atan(1.0);
1935        va = 0.0;
1936        vb = 2.0*Pi;            //orintational average, outer integral
1937        vaj = 0.0;
1938        vbj = Pi;               //endpoints of inner integral
1939
1940        summ = 0.0;                     //initialize intergral
1941
1942        scale = w[0];
1943        Dnn = w[1];                                     //Nearest neighbor distance A
1944        gg = w[2];                                      //Paracrystal distortion factor
1946        sld = w[4];
1947        sldSolv = w[5];
1948        background = w[6];
1949
1950        contrast = sld - sldSolv;
1951        //Volume fraction calculated from lattice symmetry and sphere radius
1953
1954        for(i=0;i<nordi;i++) {
1955                //setup inner integral over the ellipsoidal cross-section
1956                summj=0.0;
1957                zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;            //the outer dummy is phi
1958                for(j=0;j<nordj;j++) {
1959                        //20 gauss points for the inner integral
1960                        zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;               //the inner dummy is theta
1961                        yyy = Gauss150Wt[j] * FCC_Integrand(w,x,zi,zij);
1962                        summj += yyy;
1963                }
1964                //now calculate the value of the inner integral
1966
1967                //now calculate outer integral
1968                yyy = Gauss150Wt[i] * answer;
1969                summ += yyy;
1970        }               //final scaling is done at the end of the function, after the NT_FP64 case
1971
1973        // Multiply by contrast^2
1975        // add in the background
1977
1979}
1980
1981
1982// xx is phi (outer)
1983// yy is theta (inner)
1984double
1985FCC_Integrand(double w[], double qq, double xx, double yy) {
1986
1987        double retVal,temp1,temp3,aa,Da,Dnn,gg,Pi;
1988
1989        Pi = 4.0*atan(1.0);
1990        Dnn = w[1]; //Nearest neighbor distance A
1991        gg = w[2]; //Paracrystal distortion factor
1992        aa = Dnn;
1993        Da = gg*aa;
1994
1995        temp1 = qq*qq*Da*Da;
1996        temp3 = qq*aa;
1997
1998        retVal = FCCeval(yy,xx,temp1,temp3);
1999        retVal /=4*Pi;
2000
2001        return(retVal);
2002}
2003
2004double
2005FCCeval(double Theta, double Phi, double temp1, double temp3) {
2006
2007        double temp6,temp7,temp8,temp9,temp10;
2008        double result;
2009
2010        temp6 = sin(Theta);
2011        temp7 = sin(Theta)*sin(Phi)+cos(Theta);
2012        temp8 = -1.0*sin(Theta)*cos(Phi)+cos(Theta);
2013        temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi);
2014        temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
2015        result = pow((1.0-(temp10*temp10)),3)*temp6/((1.0-2.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp9))+(temp10*temp10)));
2016
2017        return (result);
2018}
2019
2020
2021/*      SC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x
2022
2023Uses 150 pt Gaussian quadrature for both integrals
2024
2025*/
2026double
2027SC_ParaCrystal(double w[], double x)
2028{
2029        int i,j;
2030        double Pi;
2031        double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;           //local variables of coefficient wave
2032        int nordi=150;                  //order of integration
2033        int nordj=150;
2034        double va,vb;           //upper and lower integration limits
2035        double summ,zi,yyy,answer;                      //running tally of integration
2036        double summj,vaj,vbj,zij;                       //for the inner integration
2037
2038        Pi = 4.0*atan(1.0);
2039        va = 0.0;
2040        vb = 2.0*Pi;            //orintational average, outer integral
2041        vaj = 0.0;
2042        vbj = Pi;               //endpoints of inner integral
2043
2044        summ = 0.0;                     //initialize intergral
2045
2046        scale = w[0];
2047        Dnn = w[1];                                     //Nearest neighbor distance A
2048        gg = w[2];                                      //Paracrystal distortion factor
2050        sld = w[4];
2051        sldSolv = w[5];
2052        background = w[6];
2053
2054        contrast = sld - sldSolv;
2055        //Volume fraction calculated from lattice symmetry and sphere radius
2057
2058        for(i=0;i<nordi;i++) {
2059                //setup inner integral over the ellipsoidal cross-section
2060                summj=0.0;
2061                zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;            //the outer dummy is phi
2062                for(j=0;j<nordj;j++) {
2063                        //20 gauss points for the inner integral
2064                        zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;               //the inner dummy is theta
2065                        yyy = Gauss150Wt[j] * SC_Integrand(w,x,zi,zij);
2066                        summj += yyy;
2067                }
2068                //now calculate the value of the inner integral
2070
2071                //now calculate outer integral
2072                yyy = Gauss150Wt[i] * answer;
2073                summ += yyy;
2074        }               //final scaling is done at the end of the function, after the NT_FP64 case
2075
2077        // Multiply by contrast^2
2079        // add in the background
2081
2083}
2084
2085// xx is phi (outer)
2086// yy is theta (inner)
2087double
2088SC_Integrand(double w[], double qq, double xx, double yy) {
2089
2090        double retVal,temp1,temp2,temp3,temp4,temp5,aa,Da,Dnn,gg,Pi;
2091
2092        Pi = 4.0*atan(1.0);
2093        Dnn = w[1]; //Nearest neighbor distance A
2094        gg = w[2]; //Paracrystal distortion factor
2095        aa = Dnn;
2096        Da = gg*aa;
2097
2098        temp1 = qq*qq*Da*Da;
2099        temp2 = pow( 1.0-exp(-1.0*temp1) ,3);
2100        temp3 = qq*aa;
2101        temp4 = 2.0*exp(-0.5*temp1);
2102        temp5 = exp(-1.0*temp1);
2103
2104
2105        retVal = temp2*SCeval(yy,xx,temp3,temp4,temp5);
2106        retVal /= 4*Pi;
2107
2108        return(retVal);
2109}
2110
2111double
2112SCeval(double Theta, double Phi, double temp3, double temp4, double temp5) { //Function to calculate integrand values for simple cubic structure
2113
2114        double temp6,temp7,temp8,temp9; //Theta and phi dependent parts of the equation
2115        double result;
2116
2117        temp6 = sin(Theta);
2118        temp7 = -1.0*temp3*sin(Theta)*cos(Phi);
2119        temp8 = temp3*sin(Theta)*sin(Phi);
2120        temp9 = temp3*cos(Theta);
2121        result = temp6/((1.0-temp4*cos((temp7))+temp5)*(1.0-temp4*cos((temp8))+temp5)*(1.0-temp4*cos((temp9))+temp5));
2122
2123        return (result);
2124}
2125
Note: See TracBrowser for help on using the repository browser.