source: sans/XOP_Dev/SANSAnalysis/lib/libSphere.c @ 593

Last change on this file since 593 was 593, checked in by srkline, 13 years ago

Added new model functions (see ticket #239) to the XOP, also added to the library.

File size: 56.6 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
7#include "StandardHeaders.h"                    // Include ANSI headers, Mac headers
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];
20        radius = dp[1];
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){
28                f = 4.0/3.0*pi*radius*radius*radius*delrho*delrho*scale*1.0e8 + bkg;
29                return(f);
30        }
31       
32        bes = 3.0*(sin(q*radius)-q*radius*cos(q*radius))/(q*q*q)/(radius*radius*radius);
33        vol = 4.0*pi/3.0*radius*radius*radius;
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];
144        corrad = dp[1];
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.;
161        vpoly = 4*pi/3*zp3*zp2/zp1/zp1*pow((corrad+del),3);
162       
163       
164        // remember that h is the passed in value of q for the calculation
165        y = h *del;
166        x = h *corrad;
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       
227        rad = dp[0];            // radius (A)
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];
234        sigma = 2*rad;
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
455        double inten,h1,qw,qr,width,sig,averad3,Rg2,slds,sld;
456       
457        pi = 4.0*atan(1.0);
458        x= q;
459       
460        scale = dp[0];
461        rad = dp[1];            // radius (A)
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       
472        sig = pd*rad;
473        width = sqrt(3.0)*sig;
474       
475        //x is the q-value
476        qw = x*width;
477        qr = x*rad;
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) {
485                h1 = scale*cont*cont*1.e8*4.*pi/3*pow(rad,3);
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
488                Rg2 = 3.0/5.0*rad*rad*( 1.+28.*pow(pd,2)+126.*pow(pd,4)+108.*pow(pd,6)+27.*pow(pd,8) );
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
510        // h2 = 2*sin(x*rad)*sin(x*width)-x*rad*cos(x*rad)*sin(x*width)-x*width*sin(x*rad)*cos(x*width)
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       
516        averad3 =  rad*rad*rad*(1.0+3.0*pd*pd);
517        inten /= 4.0*pi/3.0*averad3;
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];
543        rad=dp[1];
544        pd=dp[2];
545        sig=pd*rad;
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        }
555        vb = 4.0*sig +rad;
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
573        inten /= (4.0*pi/3.0*rad*rad*rad)*(1.0+3.0*pd*pd);
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];
596        rad=dp[1];              //rad is the median radius
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;
665        double scale,corrad,thick,shlrad,pp,drho1,drho2,sig,zz,bkg;             //my local names
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];
673        corrad = dp[1];
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;   
682        shlrad = corrad + thick;
683        drho1 = sld1-sld2;              //core-shell
684        drho2 = sld2-sld3;              //shell-solvent
685        zp1 = zz + 1.;
686        zp2 = zz + 2.;
687        zp3 = zz + 3.;
688        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((corrad+thick),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;
696        pp=corrad/shlrad;
697        volume=pi43*shlrad*shlrad*shlrad;
698        c1=drho1*volume;
699        c2=drho2*volume;
700       
701        arg1 = x*shlrad*pp;
702        arg2 = x*shlrad;
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       
716        //add in the background
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;
1499        double answer,zp1,zp2,zp3,vpoly,range,temp_1sf[7],pi;
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
1545        answer = (vb-va)/2.0*summ;
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);
1552        answer /= vpoly;
1553//scale
1554        answer *= scale;
1555// add in the background
1556        answer += bkg;
1557               
1558    return(answer);
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;
1566        double answer,zp1,zp2,zp3,vpoly,range,temp_2sf[9],pi;
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
1618        answer = (vb-va)/2.0*summ;
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);
1625        answer /= vpoly;
1626//scale
1627        answer *= scale;
1628// add in the background
1629        answer += bkg;
1630               
1631    return(answer);
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;
1639        double answer,zp1,zp2,zp3,vpoly,range,temp_3sf[11],pi;
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
1695        answer = (vb-va)/2.0*summ;
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);
1702        answer /= vpoly;
1703//scale
1704        answer *= scale;
1705// add in the background
1706        answer += bkg;
1707               
1708    return(answer);
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;
1716        double answer,zp1,zp2,zp3,vpoly,range,temp_4sf[13],pi;
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
1776        answer = (vb-va)/2.0*summ;
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);
1783        answer /= vpoly;
1784//scale
1785        answer *= scale;
1786// add in the background
1787        answer += bkg;
1788               
1789    return(answer);
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
1821        Rad = w[3];                                     //Sphere radius
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
1829        latticeScale = 2.0*(4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn/sqrt(3.0/4.0),3);
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
1842                answer = (vbj-vaj)/2.0*summj;
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       
1849        answer = (vb-va)/2.0*summ;
1850        // Multiply by contrast^2
1851        answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale;
1852        // add in the background
1853        answer += background;
1854       
1855        return answer;
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) {
1904                f = 4.0/3.0*pi*radius*radius*radius*delrho*delrho*1.0e8;
1905                return(f);
1906        }
1907       
1908        bes = 3.0*(sin(x*radius)-x*radius*cos(x*radius))/(x*x*x)/(radius*radius*radius);
1909        vol = 4.0*pi/3.0*radius*radius*radius;
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
1945        Rad = w[3];                                     //Sphere radius
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
1952        latticeScale = 4.0*(4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn*sqrt(2.0),3);
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
1965                answer = (vbj-vaj)/2.0*summj;
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       
1972        answer = (vb-va)/2.0*summ;
1973        // Multiply by contrast^2
1974        answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale;
1975        // add in the background
1976        answer += background;
1977       
1978        return answer;
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 = Pi/2.0;            //orintational average, outer integral
2041        vaj = 0.0;
2042        vbj = Pi/2.0;           //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
2049        Rad = w[3];                                     //Sphere radius
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
2056        latticeScale = (4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn,3);
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
2069                answer = (vbj-vaj)/2.0*summj;
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       
2076        answer = (vb-va)/2.0*summ;
2077        // Multiply by contrast^2
2078        answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale;
2079        // add in the background
2080        answer += background;
2081       
2082        return answer;
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 *= 2.0/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
2126// scattering from a uniform sphere with a Gaussian size distribution
2127//
2128double
2129FuzzySpheres(double dp[], double q)
2130{
2131        double pi,x;
2132        double scale,rad,pd,sig,rho,rhos,bkg,delrho,sig_surf,f2,bes,vol,f;              //my local names
2133        double va,vb,zi,yy,summ,inten;
2134        int nord=20,ii;
2135       
2136        pi = 4.0*atan(1.0);
2137        x= q;
2138       
2139        scale=dp[0];
2140        rad=dp[1];
2141        pd=dp[2];
2142        sig=pd*rad;
2143        sig_surf = dp[3];
2144        rho=dp[4];
2145        rhos=dp[5];
2146        delrho=rho-rhos;
2147        bkg=dp[6];
2148       
2149                       
2150        va = -4.0*sig + rad;
2151        if (va<0) {
2152                va=0;           //to avoid numerical error when  va<0 (-ve q-value)
2153        }
2154        vb = 4.0*sig +rad;
2155       
2156        summ = 0.0;             // initialize integral
2157        for(ii=0;ii<nord;ii+=1) {
2158                // calculate Gauss points on integration interval (r-value for evaluation)
2159                zi = ( Gauss20Z[ii]*(vb-va) + vb + va )/2.0;
2160                // calculate sphere scattering
2161                //
2162                //handle q==0 separately
2163                if (x==0.0) {
2164                        f2 = 4.0/3.0*pi*zi*zi*zi*delrho*delrho*1.0e8;
2165                        f2 *= exp(-0.5*sig_surf*sig_surf*x*x);
2166                        f2 *= exp(-0.5*sig_surf*sig_surf*x*x);
2167                } else {
2168                        bes = 3.0*(sin(x*zi)-x*zi*cos(x*zi))/(x*x*x)/(zi*zi*zi);
2169                        vol = 4.0*pi/3.0*zi*zi*zi;
2170                        f = vol*bes*delrho;             // [=] A
2171                        f *= exp(-0.5*sig_surf*sig_surf*x*x);
2172                        // normalize to single particle volume, convert to 1/cm
2173                        f2 = f * f / vol * 1.0e8;               // [=] 1/cm
2174                }
2175       
2176                yy = Gauss20Wt[ii] *  Gauss_distr(sig,rad,zi) * f2;
2177                yy *= 4.0*pi/3.0*zi*zi*zi;              //un-normalize by current sphere volume
2178               
2179                summ += yy;             //add to the running total of the quadrature
2180               
2181               
2182        }
2183        // calculate value of integral to return
2184        inten = (vb-va)/2.0*summ;
2185       
2186        //re-normalize by polydisperse sphere volume
2187        inten /= (4.0*pi/3.0*rad*rad*rad)*(1.0+3.0*pd*pd);
2188       
2189        inten *= scale;
2190        inten += bkg;
2191       
2192    return(inten);      //scale, and add in the background
2193}
2194
2195
Note: See TracBrowser for help on using the repository browser.