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

Last change on this file since 632 was 632, checked in by srkline, 12 years ago

incorporating Jae-Hie's changes to the library functions to explicitly handle the cases at qr==0. Also a lot of changes to explicitly declare constant values as floating point, rather than chance the interpretation as integer. (probably not necessary)

File size: 57.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        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1290        vol = 4.0*pi/3.0*rcore*rcore*rcore;
1291        f = vol*bes*contr;
1292        //now the shell
1293        qr=x*(rcore+thick);
1294        contr = rhoshel-rhosolv;
1295        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1296        vol = 4.0*pi/3.0*pow((rcore+thick),3);
1297        f += vol*bes*contr;
1298       
1299        // normalize to particle volume and rescale from [-1] to [cm-1]
1300        f2 = f*f/vol*1.0e8;
1301       
1302        //scale if desired
1303        f2 *= scale;
1304        // then add in the background
1305        f2 += bkg;
1306       
1307        return(f2);
1308}
1309
1310double
1311TwoShell(double dp[], double q)
1312{
1313        // variables are:
1314        //[0] scale factor
1315        //[1] radius of core []
1316        //[2] SLD of the core   [-2]
1317        //[3] thickness of shell 1 []
1318        //[4] SLD of shell 1
1319        //[5] thickness of shell 2 []
1320        //[6] SLD of shell 2
1321        //[7] SLD of the solvent
1322        //[8] background        [cm-1]
1323       
1324        double x,pi;
1325        double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;         //my local names
1326        double bes,f,vol,qr,contr,f2;
1327        double rhoshel2,thick2;
1328       
1329        pi = 4.0*atan(1.0);
1330        x=q;
1331       
1332        scale = dp[0];
1333        rcore = dp[1];
1334        rhocore = dp[2];
1335        thick1 = dp[3];
1336        rhoshel1 = dp[4];
1337        thick2 = dp[5];
1338        rhoshel2 = dp[6];       
1339        rhosolv = dp[7];
1340        bkg = dp[8];
1341                // core first, then add in shells
1342        qr=x*rcore;
1343        contr = rhocore-rhoshel1;
1344        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1345        vol = 4.0*pi/3.0*rcore*rcore*rcore;
1346        f = vol*bes*contr;
1347        //now the shell (1)
1348        qr=x*(rcore+thick1);
1349        contr = rhoshel1-rhoshel2;
1350        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1351        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1);
1352        f += vol*bes*contr;
1353        //now the shell (2)
1354        qr=x*(rcore+thick1+thick2);
1355        contr = rhoshel2-rhosolv;
1356        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1357        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2);
1358        f += vol*bes*contr;
1359       
1360               
1361        // normalize to particle volume and rescale from [-1] to [cm-1]
1362        f2 = f*f/vol*1.0e8;
1363       
1364        //scale if desired
1365        f2 *= scale;
1366        // then add in the background
1367        f2 += bkg;
1368       
1369        return(f2);
1370}
1371
1372double
1373ThreeShell(double dp[], double q)
1374{
1375        // variables are:
1376        //[0] scale factor
1377        //[1] radius of core []
1378        //[2] SLD of the core   [-2]
1379        //[3] thickness of shell 1 []
1380        //[4] SLD of shell 1
1381        //[5] thickness of shell 2 []
1382        //[6] SLD of shell 2
1383        //[7] thickness of shell 3
1384        //[8] SLD of shell 3
1385        //[9] SLD of solvent
1386        //[10] background       [cm-1]
1387       
1388        double x,pi;
1389        double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;         //my local names
1390        double bes,f,vol,qr,contr,f2;
1391        double rhoshel2,thick2,rhoshel3,thick3;
1392       
1393        pi = 4.0*atan(1.0);
1394        x=q;
1395       
1396        scale = dp[0];
1397        rcore = dp[1];
1398        rhocore = dp[2];
1399        thick1 = dp[3];
1400        rhoshel1 = dp[4];
1401        thick2 = dp[5];
1402        rhoshel2 = dp[6];       
1403        thick3 = dp[7];
1404        rhoshel3 = dp[8];       
1405        rhosolv = dp[9];
1406        bkg = dp[10];
1407       
1408                // core first, then add in shells
1409        qr=x*rcore;
1410        contr = rhocore-rhoshel1;
1411        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1412        vol = 4.0*pi/3.0*rcore*rcore*rcore;
1413        f = vol*bes*contr;
1414        //now the shell (1)
1415        qr=x*(rcore+thick1);
1416        contr = rhoshel1-rhoshel2;
1417        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1418        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1);
1419        f += vol*bes*contr;
1420        //now the shell (2)
1421        qr=x*(rcore+thick1+thick2);
1422        contr = rhoshel2-rhoshel3;
1423        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1424        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2);
1425        f += vol*bes*contr;
1426        //now the shell (3)
1427        qr=x*(rcore+thick1+thick2+thick3);
1428        contr = rhoshel3-rhosolv;
1429        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1430        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3);
1431        f += vol*bes*contr;
1432               
1433        // normalize to particle volume and rescale from [-1] to [cm-1]
1434        f2 = f*f/vol*1.0e8;
1435       
1436        //scale if desired
1437        f2 *= scale;
1438        // then add in the background
1439        f2 += bkg;
1440       
1441        return(f2);
1442}
1443
1444double
1445FourShell(double dp[], double q)
1446{
1447        // variables are:
1448        //[0] scale factor
1449        //[1] radius of core []
1450        //[2] SLD of the core   [-2]
1451        //[3] thickness of shell 1 []
1452        //[4] SLD of shell 1
1453        //[5] thickness of shell 2 []
1454        //[6] SLD of shell 2
1455        //[7] thickness of shell 3
1456        //[8] SLD of shell 3
1457        //[9] thickness of shell 3
1458        //[10] SLD of shell 3
1459        //[11] SLD of solvent
1460        //[12] background       [cm-1]
1461       
1462        double x,pi;
1463        double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;         //my local names
1464        double bes,f,vol,qr,contr,f2;
1465        double rhoshel2,thick2,rhoshel3,thick3,rhoshel4,thick4;
1466       
1467        pi = 4.0*atan(1.0);
1468        x=q;
1469       
1470        scale = dp[0];
1471        rcore = dp[1];
1472        rhocore = dp[2];
1473        thick1 = dp[3];
1474        rhoshel1 = dp[4];
1475        thick2 = dp[5];
1476        rhoshel2 = dp[6];       
1477        thick3 = dp[7];
1478        rhoshel3 = dp[8];
1479        thick4 = dp[9];
1480        rhoshel4 = dp[10];     
1481        rhosolv = dp[11];
1482        bkg = dp[12];
1483       
1484                // core first, then add in shells
1485        qr=x*rcore;
1486        contr = rhocore-rhoshel1;
1487        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1488        vol = 4.0*pi/3.0*rcore*rcore*rcore;
1489        f = vol*bes*contr;
1490        //now the shell (1)
1491        qr=x*(rcore+thick1);
1492        contr = rhoshel1-rhoshel2;
1493        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1494        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1);
1495        f += vol*bes*contr;
1496        //now the shell (2)
1497        qr=x*(rcore+thick1+thick2);
1498        contr = rhoshel2-rhoshel3;
1499        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1500        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2);
1501        f += vol*bes*contr;
1502        //now the shell (3)
1503        qr=x*(rcore+thick1+thick2+thick3);
1504        contr = rhoshel3-rhoshel4;
1505        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1506        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3);
1507        f += vol*bes*contr;
1508        //now the shell (4)
1509        qr=x*(rcore+thick1+thick2+thick3+thick4);
1510        contr = rhoshel4-rhosolv;
1511        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1512        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4);
1513        f += vol*bes*contr;
1514       
1515               
1516        // normalize to particle volume and rescale from [-1] to [cm-1]
1517        f2 = f*f/vol*1.0e8;
1518       
1519        //scale if desired
1520        f2 *= scale;
1521        // then add in the background
1522        f2 += bkg;
1523       
1524        return(f2);
1525}
1526
1527double
1528PolyOneShell(double dp[], double x)
1529{
1530        double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg,pd,zz;             //my local names
1531        double va,vb,summ,yyy,zi;
1532        double answer,zp1,zp2,zp3,vpoly,range,temp_1sf[7],pi;
1533        int nord=76,ii;
1534       
1535        pi = 4.0*atan(1.0);
1536       
1537        scale = dp[0];
1538        rcore = dp[1];
1539        pd = dp[2];
1540        rhocore = dp[3];
1541        thick = dp[4];
1542        rhoshel = dp[5];
1543        rhosolv = dp[6];
1544        bkg = dp[7];
1545               
1546        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only
1547       
1548        range = 8.0;                    //std deviations for the integration
1549        va = rcore*(1.0-range*pd);
1550        if (va<0.0) {
1551                va=0.0;         //otherwise numerical error when pd >= 0.3, making a<0
1552        }
1553        if (pd>0.3) {
1554                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail
1555        }
1556        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius?
1557       
1558        //temp set scale=1 and bkg=0 for quadrature calc
1559        temp_1sf[0] = 1.0;
1560        temp_1sf[1] = dp[1];    //the core radius will be changed in the loop
1561        temp_1sf[2] = dp[3];
1562        temp_1sf[3] = dp[4];
1563        temp_1sf[4] = dp[5];
1564        temp_1sf[5] = dp[6];
1565        temp_1sf[6] = 0.0;
1566       
1567        summ = 0.0;             // initialize integral
1568        for(ii=0;ii<nord;ii+=1) {
1569                // calculate Gauss points on integration interval (r-value for evaluation)
1570                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
1571                temp_1sf[1] = zi;
1572                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * OneShell(temp_1sf,x);
1573                //un-normalize by volume
1574                yyy *= 4.0*pi/3.0*pow((zi+thick),3);
1575                summ += yyy;            //add to the running total of the quadrature
1576        }
1577        // calculate value of integral to return
1578        answer = (vb-va)/2.0*summ;
1579       
1580        //re-normalize by the average volume
1581        zp1 = zz + 1.0;
1582        zp2 = zz + 2.0;
1583        zp3 = zz + 3.0;
1584        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick),3);
1585        answer /= vpoly;
1586//scale
1587        answer *= scale;
1588// add in the background
1589        answer += bkg;
1590               
1591    return(answer);
1592}
1593
1594double
1595PolyTwoShell(double dp[], double x)
1596{
1597        double scale,rcore,rhocore,rhosolv,bkg,pd,zz;           //my local names
1598        double va,vb,summ,yyy,zi;
1599        double answer,zp1,zp2,zp3,vpoly,range,temp_2sf[9],pi;
1600        int nord=76,ii;
1601        double thick1,thick2;
1602        double rhoshel1,rhoshel2;
1603       
1604        scale = dp[0];
1605        rcore = dp[1];
1606        pd = dp[2];
1607        rhocore = dp[3];
1608        thick1 = dp[4];
1609        rhoshel1 = dp[5];
1610        thick2 = dp[6];
1611        rhoshel2 = dp[7];
1612        rhosolv = dp[8];
1613        bkg = dp[9];
1614       
1615        pi = 4.0*atan(1.0);
1616               
1617        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only
1618       
1619        range = 8.0;                    //std deviations for the integration
1620        va = rcore*(1.0-range*pd);
1621        if (va<0.0) {
1622                va=0.0;         //otherwise numerical error when pd >= 0.3, making a<0
1623        }
1624        if (pd>0.3) {
1625                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail
1626        }
1627        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius?
1628       
1629        //temp set scale=1 and bkg=0 for quadrature calc
1630        temp_2sf[0] = 1.0;
1631        temp_2sf[1] = dp[1];            //the core radius will be changed in the loop
1632        temp_2sf[2] = dp[3];
1633        temp_2sf[3] = dp[4];
1634        temp_2sf[4] = dp[5];
1635        temp_2sf[5] = dp[6];
1636        temp_2sf[6] = dp[7];
1637        temp_2sf[7] = dp[8];
1638        temp_2sf[8] = 0.0;
1639       
1640        summ = 0.0;             // initialize integral
1641        for(ii=0;ii<nord;ii+=1) {
1642                // calculate Gauss points on integration interval (r-value for evaluation)
1643                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
1644                temp_2sf[1] = zi;
1645                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * TwoShell(temp_2sf,x);
1646                //un-normalize by volume
1647                yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2),3);
1648                summ += yyy;            //add to the running total of the quadrature
1649        }
1650        // calculate value of integral to return
1651        answer = (vb-va)/2.0*summ;
1652       
1653        //re-normalize by the average volume
1654        zp1 = zz + 1.0;
1655        zp2 = zz + 2.0;
1656        zp3 = zz + 3.0;
1657        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2),3);
1658        answer /= vpoly;
1659//scale
1660        answer *= scale;
1661// add in the background
1662        answer += bkg;
1663               
1664    return(answer);
1665}
1666
1667double
1668PolyThreeShell(double dp[], double x)
1669{
1670        double scale,rcore,rhocore,rhosolv,bkg,pd,zz;           //my local names
1671        double va,vb,summ,yyy,zi;
1672        double answer,zp1,zp2,zp3,vpoly,range,temp_3sf[11],pi;
1673        int nord=76,ii;
1674        double thick1,thick2,thick3;
1675        double rhoshel1,rhoshel2,rhoshel3;
1676       
1677        scale = dp[0];
1678        rcore = dp[1];
1679        pd = dp[2];
1680        rhocore = dp[3];
1681        thick1 = dp[4];
1682        rhoshel1 = dp[5];
1683        thick2 = dp[6];
1684        rhoshel2 = dp[7];
1685        thick3 = dp[8];
1686        rhoshel3 = dp[9];
1687        rhosolv = dp[10];
1688        bkg = dp[11];
1689       
1690        pi = 4.0*atan(1.0);
1691               
1692        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only
1693       
1694        range = 8.0;                    //std deviations for the integration
1695        va = rcore*(1.0-range*pd);
1696        if (va<0) {
1697                va=0;           //otherwise numerical error when pd >= 0.3, making a<0
1698        }
1699        if (pd>0.3) {
1700                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail
1701        }
1702        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius?
1703       
1704        //temp set scale=1 and bkg=0 for quadrature calc
1705        temp_3sf[0] = 1.0;
1706        temp_3sf[1] = dp[1];            //the core radius will be changed in the loop
1707        temp_3sf[2] = dp[3];
1708        temp_3sf[3] = dp[4];
1709        temp_3sf[4] = dp[5];
1710        temp_3sf[5] = dp[6];
1711        temp_3sf[6] = dp[7];
1712        temp_3sf[7] = dp[8];
1713        temp_3sf[8] = dp[9];
1714        temp_3sf[9] = dp[10];
1715        temp_3sf[10] = 0.0;
1716       
1717        summ = 0.0;             // initialize integral
1718        for(ii=0;ii<nord;ii+=1) {
1719                // calculate Gauss points on integration interval (r-value for evaluation)
1720                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
1721                temp_3sf[1] = zi;
1722                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * ThreeShell(temp_3sf,x);
1723                //un-normalize by volume
1724                yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2+thick3),3);
1725                summ += yyy;            //add to the running total of the quadrature
1726        }
1727        // calculate value of integral to return
1728        answer = (vb-va)/2.0*summ;
1729       
1730        //re-normalize by the average volume
1731        zp1 = zz + 1.0;
1732        zp2 = zz + 2.0;
1733        zp3 = zz + 3.0;
1734        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2+thick3),3);
1735        answer /= vpoly;
1736//scale
1737        answer *= scale;
1738// add in the background
1739        answer += bkg;
1740               
1741    return(answer);
1742}
1743
1744double
1745PolyFourShell(double dp[], double x)
1746{
1747        double scale,rcore,rhocore,rhosolv,bkg,pd,zz;           //my local names
1748        double va,vb,summ,yyy,zi;
1749        double answer,zp1,zp2,zp3,vpoly,range,temp_4sf[13],pi;
1750        int nord=76,ii;
1751        double thick1,thick2,thick3,thick4;
1752        double rhoshel1,rhoshel2,rhoshel3,rhoshel4;
1753       
1754        scale = dp[0];
1755        rcore = dp[1];
1756        pd = dp[2];
1757        rhocore = dp[3];
1758        thick1 = dp[4];
1759        rhoshel1 = dp[5];
1760        thick2 = dp[6];
1761        rhoshel2 = dp[7];
1762        thick3 = dp[8];
1763        rhoshel3 = dp[9];
1764        thick4 = dp[10];
1765        rhoshel4 = dp[11];
1766        rhosolv = dp[12];
1767        bkg = dp[13];
1768       
1769        pi = 4.0*atan(1.0);
1770               
1771        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only
1772       
1773        range = 8.0;                    //std deviations for the integration
1774        va = rcore*(1.0-range*pd);
1775        if (va<0) {
1776                va=0;           //otherwise numerical error when pd >= 0.3, making a<0
1777        }
1778        if (pd>0.3) {
1779                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail
1780        }
1781        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius?
1782       
1783        //temp set scale=1 and bkg=0 for quadrature calc
1784        temp_4sf[0] = 1.0;
1785        temp_4sf[1] = dp[1];            //the core radius will be changed in the loop
1786        temp_4sf[2] = dp[3];
1787        temp_4sf[3] = dp[4];
1788        temp_4sf[4] = dp[5];
1789        temp_4sf[5] = dp[6];
1790        temp_4sf[6] = dp[7];
1791        temp_4sf[7] = dp[8];
1792        temp_4sf[8] = dp[9];
1793        temp_4sf[9] = dp[10];
1794        temp_4sf[10] = dp[11];
1795        temp_4sf[11] = dp[12];
1796        temp_4sf[12] = 0.0;
1797               
1798        summ = 0.0;             // initialize integral
1799        for(ii=0;ii<nord;ii+=1) {
1800                // calculate Gauss points on integration interval (r-value for evaluation)
1801                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
1802                temp_4sf[1] = zi;
1803                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * FourShell(temp_4sf,x);
1804                //un-normalize by volume
1805                yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2+thick3+thick4),3);
1806                summ += yyy;            //add to the running total of the quadrature
1807        }
1808        // calculate value of integral to return
1809        answer = (vb-va)/2.0*summ;
1810       
1811        //re-normalize by the average volume
1812        zp1 = zz + 1.0;
1813        zp2 = zz + 2.0;
1814        zp3 = zz + 3.0;
1815        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2+thick3+thick4),3);
1816        answer /= vpoly;
1817//scale
1818        answer *= scale;
1819// add in the background
1820        answer += bkg;
1821               
1822    return(answer);
1823}
1824
1825
1826/*      BCC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x
1827
1828Uses 150 pt Gaussian quadrature for both integrals
1829
1830*/
1831double
1832BCC_ParaCrystal(double w[], double x)
1833{
1834        int i,j;
1835        double Pi;
1836        double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;           //local variables of coefficient wave
1837        int nordi=150;                  //order of integration
1838        int nordj=150;
1839        double va,vb;           //upper and lower integration limits
1840        double summ,zi,yyy,answer;                      //running tally of integration
1841        double summj,vaj,vbj,zij;                       //for the inner integration
1842       
1843        Pi = 4.0*atan(1.0);
1844        va = 0.0;
1845        vb = 2.0*Pi;            //orintational average, outer integral
1846        vaj = 0.0;
1847        vbj = Pi;               //endpoints of inner integral
1848       
1849        summ = 0.0;                     //initialize intergral
1850       
1851        scale = w[0];
1852        Dnn = w[1];                                     //Nearest neighbor distance A
1853        gg = w[2];                                      //Paracrystal distortion factor
1854        Rad = w[3];                                     //Sphere radius
1855        sld = w[4];
1856        sldSolv = w[5];
1857        background = w[6]; 
1858       
1859        contrast = sld - sldSolv;
1860       
1861        //Volume fraction calculated from lattice symmetry and sphere radius
1862        latticeScale = 2.0*(4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn/sqrt(3.0/4.0),3);
1863       
1864        for(i=0;i<nordi;i++) {
1865                //setup inner integral over the ellipsoidal cross-section
1866                summj=0.0;
1867                zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;            //the outer dummy is phi
1868                for(j=0;j<nordj;j++) {
1869                        //20 gauss points for the inner integral
1870                        zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;               //the inner dummy is theta
1871                        yyy = Gauss150Wt[j] * BCC_Integrand(w,x,zi,zij);
1872                        summj += yyy;
1873                }
1874                //now calculate the value of the inner integral
1875                answer = (vbj-vaj)/2.0*summj;
1876               
1877                //now calculate outer integral
1878                yyy = Gauss150Wt[i] * answer;
1879                summ += yyy;
1880        }               //final scaling is done at the end of the function, after the NT_FP64 case
1881       
1882        answer = (vb-va)/2.0*summ;
1883        // Multiply by contrast^2
1884        answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale;
1885        // add in the background
1886        answer += background;
1887       
1888        return answer;
1889}
1890
1891// xx is phi (outer)
1892// yy is theta (inner)
1893double
1894BCC_Integrand(double w[], double qq, double xx, double yy) {
1895       
1896        double retVal,temp1,temp3,aa,Da,Dnn,gg,Pi;
1897       
1898        Dnn = w[1]; //Nearest neighbor distance A
1899        gg = w[2]; //Paracrystal distortion factor
1900        aa = Dnn;
1901        Da = gg*aa;
1902       
1903        Pi = 4.0*atan(1.0);
1904        temp1 = qq*qq*Da*Da;
1905        temp3 = qq*aa; 
1906       
1907        retVal = BCCeval(yy,xx,temp1,temp3);
1908        retVal /=4.0*Pi;
1909       
1910        return(retVal);
1911}
1912
1913double
1914BCCeval(double Theta, double Phi, double temp1, double temp3) {
1915
1916        double temp6,temp7,temp8,temp9,temp10;
1917        double result;
1918       
1919        temp6 = sin(Theta);
1920        temp7 = sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)+cos(Theta);
1921        temp8 = -1.0*sin(Theta)*cos(Phi)-sin(Theta)*sin(Phi)+cos(Theta);
1922        temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)-cos(Theta);
1923        temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
1924        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)));
1925       
1926        return (result);
1927}
1928
1929double
1930SphereForm_Paracrystal(double radius, double delrho, double x) {                                       
1931       
1932        double bes,f,vol,f2,pi;
1933        pi = 4.0*atan(1.0);
1934        //
1935        //handle q==0 separately
1936        if(x==0) {
1937                f = 4.0/3.0*pi*radius*radius*radius*delrho*delrho*1.0e8;
1938                return(f);
1939        }
1940       
1941        bes = 3.0*(sin(x*radius)-x*radius*cos(x*radius))/(x*x*x)/(radius*radius*radius);
1942        vol = 4.0*pi/3.0*radius*radius*radius;
1943        f = vol*bes*delrho      ;       // [=] 
1944        // normalize to single particle volume, convert to 1/cm
1945        f2 = f * f / vol * 1.0e8;               // [=] 1/cm
1946       
1947        return (f2);
1948}
1949
1950/*      FCC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x
1951
1952Uses 150 pt Gaussian quadrature for both integrals
1953
1954*/
1955double
1956FCC_ParaCrystal(double w[], double x)
1957{
1958        int i,j;
1959        double Pi;
1960        double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;           //local variables of coefficient wave
1961        int nordi=150;                  //order of integration
1962        int nordj=150;
1963        double va,vb;           //upper and lower integration limits
1964        double summ,zi,yyy,answer;                      //running tally of integration
1965        double summj,vaj,vbj,zij;                       //for the inner integration
1966       
1967        Pi = 4.0*atan(1.0);
1968        va = 0.0;
1969        vb = 2.0*Pi;            //orintational average, outer integral
1970        vaj = 0.0;
1971        vbj = Pi;               //endpoints of inner integral
1972       
1973        summ = 0.0;                     //initialize intergral
1974       
1975        scale = w[0];
1976        Dnn = w[1];                                     //Nearest neighbor distance A
1977        gg = w[2];                                      //Paracrystal distortion factor
1978        Rad = w[3];                                     //Sphere radius
1979        sld = w[4];
1980        sldSolv = w[5];
1981        background = w[6]; 
1982       
1983        contrast = sld - sldSolv;
1984        //Volume fraction calculated from lattice symmetry and sphere radius
1985        latticeScale = 4.0*(4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn*sqrt(2.0),3);
1986       
1987        for(i=0;i<nordi;i++) {
1988                //setup inner integral over the ellipsoidal cross-section
1989                summj=0.0;
1990                zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;            //the outer dummy is phi
1991                for(j=0;j<nordj;j++) {
1992                        //20 gauss points for the inner integral
1993                        zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;               //the inner dummy is theta
1994                        yyy = Gauss150Wt[j] * FCC_Integrand(w,x,zi,zij);
1995                        summj += yyy;
1996                }
1997                //now calculate the value of the inner integral
1998                answer = (vbj-vaj)/2.0*summj;
1999               
2000                //now calculate outer integral
2001                yyy = Gauss150Wt[i] * answer;
2002                summ += yyy;
2003        }               //final scaling is done at the end of the function, after the NT_FP64 case
2004       
2005        answer = (vb-va)/2.0*summ;
2006        // Multiply by contrast^2
2007        answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale;
2008        // add in the background
2009        answer += background;
2010       
2011        return answer;
2012}
2013
2014
2015// xx is phi (outer)
2016// yy is theta (inner)
2017double
2018FCC_Integrand(double w[], double qq, double xx, double yy) {
2019       
2020        double retVal,temp1,temp3,aa,Da,Dnn,gg,Pi;
2021       
2022        Pi = 4.0*atan(1.0);
2023        Dnn = w[1]; //Nearest neighbor distance A
2024        gg = w[2]; //Paracrystal distortion factor
2025        aa = Dnn;
2026        Da = gg*aa;
2027       
2028        temp1 = qq*qq*Da*Da;
2029        temp3 = qq*aa;
2030       
2031        retVal = FCCeval(yy,xx,temp1,temp3);
2032        retVal /=4*Pi;
2033       
2034        return(retVal);
2035}
2036
2037double
2038FCCeval(double Theta, double Phi, double temp1, double temp3) {
2039
2040        double temp6,temp7,temp8,temp9,temp10;
2041        double result;
2042       
2043        temp6 = sin(Theta);
2044        temp7 = sin(Theta)*sin(Phi)+cos(Theta);
2045        temp8 = -1.0*sin(Theta)*cos(Phi)+cos(Theta);
2046        temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi);
2047        temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
2048        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)));
2049       
2050        return (result);
2051}
2052
2053
2054/*      SC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x
2055
2056Uses 150 pt Gaussian quadrature for both integrals
2057
2058*/
2059double
2060SC_ParaCrystal(double w[], double x)
2061{
2062        int i,j;
2063        double Pi;
2064        double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;           //local variables of coefficient wave
2065        int nordi=150;                  //order of integration
2066        int nordj=150;
2067        double va,vb;           //upper and lower integration limits
2068        double summ,zi,yyy,answer;                      //running tally of integration
2069        double summj,vaj,vbj,zij;                       //for the inner integration
2070       
2071        Pi = 4.0*atan(1.0);
2072        va = 0.0;
2073        vb = Pi/2.0;            //orintational average, outer integral
2074        vaj = 0.0;
2075        vbj = Pi/2.0;           //endpoints of inner integral
2076       
2077        summ = 0.0;                     //initialize intergral
2078       
2079        scale = w[0];
2080        Dnn = w[1];                                     //Nearest neighbor distance A
2081        gg = w[2];                                      //Paracrystal distortion factor
2082        Rad = w[3];                                     //Sphere radius
2083        sld = w[4];
2084        sldSolv = w[5];
2085        background = w[6]; 
2086       
2087        contrast = sld - sldSolv;
2088        //Volume fraction calculated from lattice symmetry and sphere radius
2089        latticeScale = (4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn,3);
2090       
2091        for(i=0;i<nordi;i++) {
2092                //setup inner integral over the ellipsoidal cross-section
2093                summj=0.0;
2094                zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;            //the outer dummy is phi
2095                for(j=0;j<nordj;j++) {
2096                        //20 gauss points for the inner integral
2097                        zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;               //the inner dummy is theta
2098                        yyy = Gauss150Wt[j] * SC_Integrand(w,x,zi,zij);
2099                        summj += yyy;
2100                }
2101                //now calculate the value of the inner integral
2102                answer = (vbj-vaj)/2.0*summj;
2103               
2104                //now calculate outer integral
2105                yyy = Gauss150Wt[i] * answer;
2106                summ += yyy;
2107        }               //final scaling is done at the end of the function, after the NT_FP64 case
2108       
2109        answer = (vb-va)/2.0*summ;
2110        // Multiply by contrast^2
2111        answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale;
2112        // add in the background
2113        answer += background;
2114       
2115        return answer;
2116}
2117
2118// xx is phi (outer)
2119// yy is theta (inner)
2120double
2121SC_Integrand(double w[], double qq, double xx, double yy) {
2122       
2123        double retVal,temp1,temp2,temp3,temp4,temp5,aa,Da,Dnn,gg,Pi;
2124       
2125        Pi = 4.0*atan(1.0);
2126        Dnn = w[1]; //Nearest neighbor distance A
2127        gg = w[2]; //Paracrystal distortion factor
2128        aa = Dnn;
2129        Da = gg*aa;
2130       
2131        temp1 = qq*qq*Da*Da;
2132        temp2 = pow( 1.0-exp(-1.0*temp1) ,3);
2133        temp3 = qq*aa;
2134        temp4 = 2.0*exp(-0.5*temp1);
2135        temp5 = exp(-1.0*temp1);
2136       
2137       
2138        retVal = temp2*SCeval(yy,xx,temp3,temp4,temp5);
2139        retVal *= 2.0/Pi;
2140       
2141        return(retVal);
2142}
2143
2144double
2145SCeval(double Theta, double Phi, double temp3, double temp4, double temp5) { //Function to calculate integrand values for simple cubic structure
2146
2147        double temp6,temp7,temp8,temp9; //Theta and phi dependent parts of the equation
2148        double result;
2149       
2150        temp6 = sin(Theta);
2151        temp7 = -1.0*temp3*sin(Theta)*cos(Phi);
2152        temp8 = temp3*sin(Theta)*sin(Phi);
2153        temp9 = temp3*cos(Theta);
2154        result = temp6/((1.0-temp4*cos((temp7))+temp5)*(1.0-temp4*cos((temp8))+temp5)*(1.0-temp4*cos((temp9))+temp5));
2155       
2156        return (result);
2157}
2158
2159// scattering from a uniform sphere with a Gaussian size distribution
2160//
2161double
2162FuzzySpheres(double dp[], double q)
2163{
2164        double pi,x;
2165        double scale,rad,pd,sig,rho,rhos,bkg,delrho,sig_surf,f2,bes,vol,f;              //my local names
2166        double va,vb,zi,yy,summ,inten;
2167        int nord=20,ii;
2168       
2169        pi = 4.0*atan(1.0);
2170        x= q;
2171       
2172        scale=dp[0];
2173        rad=dp[1];
2174        pd=dp[2];
2175        sig=pd*rad;
2176        sig_surf = dp[3];
2177        rho=dp[4];
2178        rhos=dp[5];
2179        delrho=rho-rhos;
2180        bkg=dp[6];
2181       
2182                       
2183        va = -4.0*sig + rad;
2184        if (va<0) {
2185                va=0;           //to avoid numerical error when  va<0 (-ve q-value)
2186        }
2187        vb = 4.0*sig +rad;
2188       
2189        summ = 0.0;             // initialize integral
2190        for(ii=0;ii<nord;ii+=1) {
2191                // calculate Gauss points on integration interval (r-value for evaluation)
2192                zi = ( Gauss20Z[ii]*(vb-va) + vb + va )/2.0;
2193                // calculate sphere scattering
2194                //
2195                //handle q==0 separately
2196                if (x==0.0) {
2197                        f2 = 4.0/3.0*pi*zi*zi*zi*delrho*delrho*1.0e8;
2198                        f2 *= exp(-0.5*sig_surf*sig_surf*x*x);
2199                        f2 *= exp(-0.5*sig_surf*sig_surf*x*x);
2200                } else {
2201                        bes = 3.0*(sin(x*zi)-x*zi*cos(x*zi))/(x*x*x)/(zi*zi*zi);
2202                        vol = 4.0*pi/3.0*zi*zi*zi;
2203                        f = vol*bes*delrho;             // [=] A
2204                        f *= exp(-0.5*sig_surf*sig_surf*x*x);
2205                        // normalize to single particle volume, convert to 1/cm
2206                        f2 = f * f / vol * 1.0e8;               // [=] 1/cm
2207                }
2208       
2209                yy = Gauss20Wt[ii] *  Gauss_distr(sig,rad,zi) * f2;
2210                yy *= 4.0*pi/3.0*zi*zi*zi;              //un-normalize by current sphere volume
2211               
2212                summ += yy;             //add to the running total of the quadrature
2213               
2214               
2215        }
2216        // calculate value of integral to return
2217        inten = (vb-va)/2.0*summ;
2218       
2219        //re-normalize by polydisperse sphere volume
2220        inten /= (4.0*pi/3.0*rad*rad*rad)*(1.0+3.0*pd*pd);
2221       
2222        inten *= scale;
2223        inten += bkg;
2224       
2225    return(inten);      //scale, and add in the background
2226}
2227
2228
Note: See TracBrowser for help on using the repository browser.