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

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