source: sans/Analysis/branches/ajj_23APR07/XOPs/SANSAnalysis/lib/libSphere.c @ 235

Last change on this file since 235 was 235, checked in by srkline, 15 years ago

Changes to XOP library functions to convert models to take individual SLD's rather than contrast.

Individual SLD's are easier to work with - since they are experimental values, contrast is not.

This, unfortunately, means that what Mathieu has done may need to be tweaked, since our library has been tweaked.

File size: 31.8 KB
Line 
1/*      SimpleFit.c
2
3A simplified project designed to act as a template for your curve fitting function.
4The fitting function is a simple polynomial. It works but is of no practical use.
5*/
6
7#include "StandardHeaders.h"                    // Include ANSI headers, Mac headers
8#include "GaussWeights.h"
9#include "libSphere.h"
10
11// scattering from a sphere - hardly needs to be an XOP...
12double
13SphereForm(double dp[], double q)
14{
15        double scale,radius,delrho,bkg,sldSph,sldSolv;          //my local names
16        double bes,f,vol,f2,pi;
17       
18        pi = 4.0*atan(1.0);
19        scale = dp[0];
20        radius = dp[1];
21        sldSph = dp[2];
22        sldSolv = dp[3];
23        bkg = dp[4];
24       
25        delrho = sldSph - sldSolv;
26        //handle q==0 separately
27        if(q==0){
28                f = 4.0/3.0*pi*radius*radius*radius*delrho*delrho*scale*1.0e8 + bkg;
29                return(f);
30        }
31       
32        bes = 3.0*(sin(q*radius)-q*radius*cos(q*radius))/(q*q*q)/(radius*radius*radius);
33        vol = 4.0*pi/3.0*radius*radius*radius;
34        f = vol*bes*delrho;             // [=] 
35                                                        // normalize to single particle volume, convert to 1/cm
36        f2 = f * f / vol * 1.0e8;               // [=] 1/cm
37       
38        return(scale*f2+bkg);   //scale, and add in the background
39}
40
41// scattering from a monodisperse core-shell sphere - hardly needs to be an XOP...
42double
43CoreShellForm(double dp[], double q)
44{
45        double x,pi;
46        double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg;           //my local names
47        double bes,f,vol,qr,contr,f2;
48       
49        pi = 4.0*atan(1.0);
50        x=q;
51       
52        scale = dp[0];
53        rcore = dp[1];
54        thick = dp[2];
55        rhocore = dp[3];
56        rhoshel = dp[4];
57        rhosolv = dp[5];
58        bkg = dp[6];
59        // core first, then add in shell
60        qr=x*rcore;
61        contr = rhocore-rhoshel;
62        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
63        vol = 4.0*pi/3.0*rcore*rcore*rcore;
64        f = vol*bes*contr;
65        //now the shell
66        qr=x*(rcore+thick);
67        contr = rhoshel-rhosolv;
68        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
69        vol = 4.0*pi/3.0*pow((rcore+thick),3);
70        f += vol*bes*contr;
71       
72        // normalize to particle volume and rescale from [-1] to [cm-1]
73        f2 = f*f/vol*1.0e8;
74       
75        //scale if desired
76        f2 *= scale;
77        // then add in the background
78        f2 += bkg;
79       
80        return(f2);
81}
82
83// scattering from a unilamellar vesicle - hardly needs to be an XOP...
84// same functional form as the core-shell sphere, but more intuitive for a vesicle
85double
86VesicleForm(double dp[], double q)
87{
88        double x,pi;
89        double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg;           //my local names
90        double bes,f,vol,qr,contr,f2;
91        pi = 4.0*atan(1.0);
92        x= q;
93       
94        scale = dp[0];
95        rcore = dp[1];
96        thick = dp[2];
97        rhocore = dp[3];
98        rhosolv = rhocore;
99        rhoshel = dp[4];
100        bkg = dp[5];
101        // core first, then add in shell
102        qr=x*rcore;
103        contr = rhocore-rhoshel;
104        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
105        vol = 4.0*pi/3.0*rcore*rcore*rcore;
106        f = vol*bes*contr;
107        //now the shell
108        qr=x*(rcore+thick);
109        contr = rhoshel-rhosolv;
110        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
111        vol = 4.0*pi/3.0*pow((rcore+thick),3);
112        f += vol*bes*contr;
113       
114        // normalize to the particle volume and rescale from [-1] to [cm-1]
115        //note that for the vesicle model, the volume is ONLY the shell volume
116        vol = 4.0*pi/3.0*(pow((rcore+thick),3)-pow(rcore,3));
117        f2 = f*f/vol*1.0e8;
118       
119        //scale if desired
120        f2 *= scale;
121        // then add in the background
122        f2 += bkg;
123       
124        return(f2);
125}
126
127
128// scattering from a core shell sphere with a (Schulz) polydisperse core and constant shell thickness
129//
130double
131PolyCoreForm(double dp[], double q)
132{
133        double pi;
134        double scale,corrad,sig,zz,del,drho1,drho2,form,bkg;            //my local names
135        double d, g ,h;
136        double qq, x, y, c1, c2, c3, c4, c5, c6, c7, c8, c9, t1, t2, t3;
137        double t4, t5, tb, cy, sy, tb1, tb2, tb3, c2y, zp1, zp2;
138        double zp3,vpoly;
139        double s2y, arg1, arg2, arg3, drh1, drh2;
140       
141        pi = 4.0*atan(1.0);
142        qq= q;
143        scale = dp[0];
144        corrad = dp[1];
145        sig = dp[2];
146        del = dp[3];
147        drho1 = dp[4]-dp[5];            //core-shell
148        drho2 = dp[5]-dp[6];            //shell-solvent
149        bkg = dp[7];
150       
151        zz = (1.0/sig)*(1.0/sig) - 1.0;   
152       
153        h=qq;
154       
155        drh1 = drho1;
156        drh2 = drho2;
157        g = drh2 * -1. / drh1;
158        zp1 = zz + 1.;
159        zp2 = zz + 2.;
160        zp3 = zz + 3.;
161        vpoly = 4*pi/3*zp3*zp2/zp1/zp1*pow((corrad+del),3);
162       
163       
164        // remember that h is the passed in value of q for the calculation
165        y = h *del;
166        x = h *corrad;
167        d = atan(x * 2. / zp1);
168        arg1 = zp1 * d;
169        arg2 = zp2 * d;
170        arg3 = zp3 * d;
171        sy = sin(y);
172        cy = cos(y);
173        s2y = sin(y * 2.);
174        c2y = cos(y * 2.);
175        c1 = .5 - g * (cy + y * sy) + g * g * .5 * (y * y + 1.);
176        c2 = g * y * (g - cy);
177        c3 = (g * g + 1.) * .5 - g * cy;
178        c4 = g * g * (y * cy - sy) * (y * cy - sy) - c1;
179        c5 = g * 2. * sy * (1. - g * (y * sy + cy)) + c2;
180        c6 = c3 - g * g * sy * sy;
181        c7 = g * sy - g * .5 * g * (y * y + 1.) * s2y - c5;
182        c8 = c4 - .5 + g * cy - g * .5 * g * (y * y + 1.) * c2y;
183        c9 = g * sy * (1. - g * cy);
184       
185        tb = log(zp1 * zp1 / (zp1 * zp1 + x * 4. * x));
186        tb1 = exp(zp1 * .5 * tb);
187        tb2 = exp(zp2 * .5 * tb);
188        tb3 = exp(zp3 * .5 * tb);
189       
190        t1 = c1 + c2 * x + c3 * x * x * zp2 / zp1;
191        t2 = tb1 * (c4 * cos(arg1) + c7 * sin(arg1));
192        t3 = x * tb2 * (c5 * cos(arg2) + c8 * sin(arg2));
193        t4 = zp2 / zp1 * x * x * tb3 * (c6 * cos(arg3) + c9 * sin(arg3));
194        t5 = t1 + t2 + t3 + t4;
195        form = t5 * 16. * pi * pi * drh1 * drh1 / pow(qq,6);
196        //      normalize by the average volume !!! corrected for polydispersity
197        // and convert to cm-1
198        form /= vpoly;
199        form *= 1.0e8;
200        //Scale
201        form *= scale;
202        // then add in the background
203        form += bkg;
204       
205        return(form);
206}
207
208// scattering from a uniform sphere with a (Schulz) size distribution
209// structure factor effects are explicitly and correctly included.
210//
211double
212PolyHardSphereIntensity(double dp[], double q)
213{
214        double pi;
215        double rad,z2,phi,cont,bkg,sigma;               //my local names
216        double mu,mu1,d1,d2,d3,d4,d5,d6,capd,rho;
217        double ll,l1,bb,cc,chi,chi1,chi2,ee,t1,t2,t3,pp;
218        double ka,zz,v1,v2,p1,p2;
219        double h1,h2,h3,h4,e1,yy,y1,s1,s2,s3,hint1,hint2;
220        double capl,capl1,capmu,capmu1,r3,pq;
221        double ka2,r1,heff;
222        double hh,k,slds,sld;
223       
224        pi = 4.0*atan(1.0);
225        k= q;
226       
227        rad = dp[0];            // radius (A)
228        z2 = dp[1];             //polydispersity (0<z2<1)
229        phi = dp[2];            // volume fraction (0<phi<1)
230        slds = dp[3];
231        sld = dp[4];
232        cont = (slds - sld)*1.0e4;              // contrast (odd units)
233        bkg = dp[5];
234        sigma = 2*rad;
235       
236        zz=1.0/(z2*z2)-1.0;
237        bb = sigma/(zz+1.0);
238        cc = zz+1.0;
239       
240        //*c   Compute the number density by <r-cubed>, not <r> cubed*/
241        r1 = sigma/2.0;
242        r3 = r1*r1*r1;
243        r3 *= (zz+2.0)*(zz+3.0)/((zz+1.0)*(zz+1.0));
244        rho=phi/(1.3333333333*pi*r3);
245        t1 = rho*bb*cc;
246        t2 = rho*bb*bb*cc*(cc+1.0);
247        t3 = rho*bb*bb*bb*cc*(cc+1.0)*(cc+2.0);
248        capd = 1.0-pi*t3/6.0;
249        //************
250        v1=1.0/(1.0+bb*bb*k*k);
251        v2=1.0/(4.0+bb*bb*k*k);
252        pp=pow(v1,(cc/2.0))*sin(cc*atan(bb*k));
253        p1=bb*cc*pow(v1,((cc+1.0)/2.0))*sin((cc+1.0)*atan(bb*k));
254        p2=cc*(cc+1.0)*bb*bb*pow(v1,((cc+2.0)/2.0))*sin((cc+2.0)*atan(bb*k));
255        mu=pow(2,cc)*pow(v2,(cc/2.0))*sin(cc*atan(bb*k/2.0));
256        mu1=pow(2,(cc+1.0))*bb*cc*pow(v2,((cc+1.0)/2.0))*sin((cc+1.0)*atan(k*bb/2.0));
257        s1=bb*cc;
258        s2=cc*(cc+1.0)*bb*bb;
259        s3=cc*(cc+1.0)*(cc+2.0)*bb*bb*bb;
260        chi=pow(v1,(cc/2.0))*cos(cc*atan(bb*k));
261        chi1=bb*cc*pow(v1,((cc+1.0)/2.0))*cos((cc+1.0)*atan(bb*k));
262        chi2=cc*(cc+1.0)*bb*bb*pow(v1,((cc+2.0)/2.0))*cos((cc+2.0)*atan(bb*k));
263        ll=pow(2,cc)*pow(v2,(cc/2.0))*cos(cc*atan(bb*k/2.0));
264        l1=pow(2,(cc+1.0))*bb*cc*pow(v2,((cc+1.0)/2.0))*cos((cc+1.0)*atan(k*bb/2.0));
265        d1=(pi/capd)*(2.0+(pi/capd)*(t3-(rho/k)*(k*s3-p2)));
266        d2=pow((pi/capd),2)*(rho/k)*(k*s2-p1);
267        d3=(-1.0)*pow((pi/capd),2)*(rho/k)*(k*s1-pp);
268        d4=(pi/capd)*(k-(pi/capd)*(rho/k)*(chi1-s1));
269        d5=pow((pi/capd),2)*((rho/k)*(chi-1.0)+0.5*k*t2);
270        d6=pow((pi/capd),2)*(rho/k)*(chi2-s2);
271       
272        e1=pow((pi/capd),2)*pow((rho/k/k),2)*((chi-1.0)*(chi2-s2)-(chi1-s1)*(chi1-s1)-(k*s1-pp)*(k*s3-p2)+pow((k*s2-p1),2));
273        ee=1.0-(2.0*pi/capd)*(1.0+0.5*pi*t3/capd)*(rho/k/k/k)*(k*s1-pp)-(2.0*pi/capd)*rho/k/k*((chi1-s1)+(0.25*pi*t2/capd)*(chi2-s2))-e1;
274        y1=pow((pi/capd),2)*pow((rho/k/k),2)*((k*s1-pp)*(chi2-s2)-2.0*(k*s2-p1)*(chi1-s1)+(k*s3-p2)*(chi-1.0));
275        yy = (2.0*pi/capd)*(1.0+0.5*pi*t3/capd)*(rho/k/k/k)*(chi+0.5*k*k*s2-1.0)-(2.0*pi*rho/capd/k/k)*(k*s2-p1+(0.25*pi*t2/capd)*(k*s3-p2))-y1;   
276       
277        capl=2.0*pi*cont*rho/k/k/k*(pp-0.5*k*(s1+chi1));
278        capl1=2.0*pi*cont*rho/k/k/k*(p1-0.5*k*(s2+chi2));
279        capmu=2.0*pi*cont*rho/k/k/k*(1.0-chi-0.5*k*p1);
280        capmu1=2.0*pi*cont*rho/k/k/k*(s1-chi1-0.5*k*p2);
281       
282        h1=capl*(capl*(yy*d1-ee*d6)+capl1*(yy*d2-ee*d4)+capmu*(ee*d1+yy*d6)+capmu1*(ee*d2+yy*d4));
283        h2=capl1*(capl*(yy*d2-ee*d4)+capl1*(yy*d3-ee*d5)+capmu*(ee*d2+yy*d4)+capmu1*(ee*d3+yy*d5));
284        h3=capmu*(capl*(ee*d1+yy*d6)+capl1*(ee*d2+yy*d4)+capmu*(ee*d6-yy*d1)+capmu1*(ee*d4-yy*d2));
285        h4=capmu1*(capl*(ee*d2+yy*d4)+capl1*(ee*d3+yy*d5)+capmu*(ee*d4-yy*d2)+capmu1*(ee*d5-yy*d3));
286       
287        //*  This part computes the second integral in equation (1) of the paper.*/
288       
289        hint1 = -2.0*(h1+h2+h3+h4)/(k*k*k*(ee*ee+yy*yy));
290       
291        //*  This part computes the first integral in equation (1).  It also
292        // generates the KC approximated effective structure factor.*/
293       
294        pq=4.0*pi*cont*(sin(k*sigma/2.0)-0.5*k*sigma*cos(k*sigma/2.0));
295        hint2=8.0*pi*pi*rho*cont*cont/(k*k*k*k*k*k)*(1.0-chi-k*p1+0.25*k*k*(s2+chi2));
296       
297        ka=k*(sigma/2.0);
298        //
299        hh=hint1+hint2;         // this is the model intensity
300                                                //
301        heff=1.0+hint1/hint2;
302        ka2=ka*ka;
303        //*
304        //  heff is PY analytical solution for intensity divided by the
305        //   form factor.  happ is the KC approximated effective S(q)
306         
307         //*******************
308         //  add in the background then return the intensity value
309         
310         return(hh+bkg);        //scale, and add in the background
311}
312
313// scattering from a uniform sphere with a (Schulz) size distribution, bimodal population
314// NO CROSS TERM IS ACCOUNTED FOR == DILUTE SOLUTION!!
315//
316double
317BimodalSchulzSpheres(double dp[], double q)
318{
319        double x,pq;
320        double scale,ravg,pd,bkg,rho,rhos;              //my local names
321        double scale2,ravg2,pd2,rho2;           //my local names
322       
323        x= q;
324       
325        scale = dp[0];
326        ravg = dp[1];
327        pd = dp[2];
328        rho = dp[3];
329        scale2 = dp[4];
330        ravg2 = dp[5];
331        pd2 = dp[6];
332        rho2 = dp[7];
333        rhos = dp[8];
334        bkg = dp[9];
335       
336        pq = SchulzSphere_Fn( scale,  ravg,  pd,  rho,  rhos, x);
337        pq += SchulzSphere_Fn( scale2,  ravg2,  pd2,  rho2,  rhos, x);
338        // add in the background
339        pq += bkg;
340       
341        return (pq);
342}
343
344// scattering from a uniform sphere with a (Schulz) size distribution
345//
346double
347SchulzSpheres(double dp[], double q)
348{
349        double x,pq;
350        double scale,ravg,pd,bkg,rho,rhos;              //my local names
351       
352        x= q;
353       
354        scale = dp[0];
355        ravg = dp[1];
356        pd = dp[2];
357        rho = dp[3];
358        rhos = dp[4];
359        bkg = dp[5];
360        pq = SchulzSphere_Fn( scale,  ravg,  pd,  rho,  rhos, x);
361        // add in the background
362        pq += bkg;
363       
364        return(pq);
365}
366
367// calculates everything but the background
368double
369SchulzSphere_Fn(double scale, double ravg, double pd, double rho, double rhos, double x)
370{
371        double zp1,zp2,zp3,zp4,zp5,zp6,zp7,vpoly;
372        double aa,at1,at2,rt1,rt2,rt3,t1,t2,t3;
373        double v1,v2,v3,g1,pq,pi,delrho,zz;
374        double i_zero,Rg2,zp8;
375       
376        pi = 4.0*atan(1.0);
377        delrho = rho-rhos;
378        zz = (1.0/pd)*(1.0/pd) - 1.0;   
379       
380        zp1 = zz + 1.0;
381        zp2 = zz + 2.0;
382        zp3 = zz + 3.0;
383        zp4 = zz + 4.0;
384        zp5 = zz + 5.0;
385        zp6 = zz + 6.0;
386        zp7 = zz + 7.0;
387        //
388       
389        //small QR limit - use Guinier approx
390        zp8 = zz+8;
391        if(x*ravg < 0.1) {
392                i_zero = scale*delrho*delrho*1.e8*4.*pi/3.*pow(ravg,3);
393                i_zero *= zp6*zp5*zp4/zp1/zp1/zp1;              //6th moment / 3rd moment
394                Rg2 = 3.*zp8*zp7/5./(zp1*zp1)*ravg*ravg;
395                pq = i_zero*exp(-x*x*Rg2/3.);
396                //pq += bkg;            //unlike the Igor code, the backgorund is added in the wrapper (above)
397                return(pq);
398        }
399        //
400
401        aa = (zz+1)/x/ravg;
402       
403        at1 = atan(1.0/aa);
404        at2 = atan(2.0/aa);
405        //
406        //  calculations are performed to avoid  large # errors
407        // - trick is to propogate the a^(z+7) term through the g1
408        //
409        t1 = zp7*log10(aa) - zp1/2.0*log10(aa*aa+4.0);
410        t2 = zp7*log10(aa) - zp3/2.0*log10(aa*aa+4.0);
411        t3 = zp7*log10(aa) - zp2/2.0*log10(aa*aa+4.0);
412        //      print t1,t2,t3
413        rt1 = pow(10,t1);
414        rt2 = pow(10,t2);
415        rt3 = pow(10,t3);
416        v1 = pow(aa,6) - rt1*cos(zp1*at2);
417        v2 = zp1*zp2*( pow(aa,4) + rt2*cos(zp3*at2) );
418        v3 = -2.0*zp1*rt3*sin(zp2*at2);
419        g1 = (v1+v2+v3);
420       
421        pq = log10(g1) - 6.0*log10(zp1) + 6.0*log10(ravg);
422        pq = pow(10,pq)*8*pi*pi*delrho*delrho;
423       
424        //
425        // beta factor is not used here, but could be for the
426        // decoupling approximation
427        //
428        //      g11 = g1
429        //      gd = -zp7*log(aa)
430        //      g1 = log(g11) + gd
431        //                       
432        //      t1 = zp1*at1
433        //      t2 = zp2*at1
434        //      g2 = sin( t1 ) - zp1/sqrt(aa*aa+1)*cos( t2 )
435        //      g22 = g2*g2
436        //      beta = zp1*log(aa) - zp1*log(aa*aa+1) - g1 + log(g22)
437        //      beta = 2*alog(beta)
438       
439        //re-normalize by the average volume
440        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*ravg*ravg*ravg;
441        pq /= vpoly;
442        //scale, convert to cm^-1
443        pq *= scale * 1.0e8;
444   
445    return(pq);
446}
447
448// scattering from a uniform sphere with a rectangular size distribution
449//
450double
451PolyRectSpheres(double dp[], double q)
452{
453        double pi,x;
454        double scale,rad,pd,cont,bkg;           //my local names
455        double inten,h1,qw,qr,width,sig,averad3,Rg2,slds,sld;
456       
457        pi = 4.0*atan(1.0);
458        x= q;
459       
460        scale = dp[0];
461        rad = dp[1];            // radius (A)
462        pd = dp[2];             //polydispersity of rectangular distribution
463        slds = dp[3];
464        sld = dp[4];
465        cont = slds - sld;              // contrast (A^-2)
466        bkg = dp[5];
467       
468        // as usual, poly = sig/ravg
469        // for the rectangular distribution, sig = width/sqrt(3)
470        // width is the HALF- WIDTH of the rectangular distrubution
471       
472        sig = pd*rad;
473        width = sqrt(3.0)*sig;
474       
475        //x is the q-value
476        qw = x*width;
477        qr = x*rad;
478       
479        // as for the numerical inatabilities at low QR, the function is calculating the sines and cosines
480        // just fine - the problem seems to be that the
481        // leading terms nearly cancel with the last term (the -6*qr... term), to within machine
482        // precision - the difference is on the order of 10^-20
483        // so just use the limiting Guiner value
484        if(qr<0.1) {
485                h1 = scale*cont*cont*1.e8*4.*pi/3*pow(rad,3);
486                h1 *=   (1. + 15.*pow(pd,2) + 27.*pow(pd,4) +27./7.*pow(pd,6) );                                //6th moment
487                h1 /= (1.+3.*pd*pd);    //3rd moment
488                Rg2 = 3.0/5.0*rad*rad*( 1.+28.*pow(pd,2)+126.*pow(pd,4)+108.*pow(pd,6)+27.*pow(pd,8) );
489                Rg2 /= (1.+15.*pow(pd,2)+27.*pow(pd,4)+27./7.*pow(pd,6));
490                h1 *= exp(-1./3.*Rg2*x*x);
491                h1 += bkg;
492                return(h1);
493        }
494       
495        // normal calculation
496        h1 = -0.5*qw + qr*qr*qw + (qw*qw*qw)/3.0;
497        h1 -= 5.0/2.0*cos(2*qr)*sin(qw)*cos(qw);
498        h1 += 0.5*qr*qr*cos(2*qr)*sin(2*qw);
499        h1 += 0.5*qw*qw*cos(2*qr)*sin(2*qw);
500        h1 += qw*qr*sin(2*qr)*cos(2*qw);
501        h1 += 3.0*qw*(cos(qr)*cos(qw))*(cos(qr)*cos(qw));
502        h1+= 3.0*qw*(sin(qr)*sin(qw))*(sin(qr)*sin(qw));
503        h1 -= 6.0*qr*cos(qr)*sin(qr)*cos(qw)*sin(qw);
504       
505        // calculate P(q) = <f^2>
506        inten = 8.0*pi*pi*cont*cont/width/pow(x,7)*h1;
507       
508        // beta(q) would be calculated as 2/width/x/h1*h2*h2
509        // with
510        // h2 = 2*sin(x*rad)*sin(x*width)-x*rad*cos(x*rad)*sin(x*width)-x*width*sin(x*rad)*cos(x*width)
511       
512        // normalize to the average volume
513        // <R^3> = ravg^3*(1+3*pd^2)
514        // or... "zf"  = (1 + 3*p^2), which will be greater than one
515       
516        averad3 =  rad*rad*rad*(1.0+3.0*pd*pd);
517        inten /= 4.0*pi/3.0*averad3;
518        //resacle to 1/cm
519        inten *= 1.0e8;
520        //scale the result
521        inten *= scale;
522        // then add in the background
523        inten += bkg;
524       
525        return(inten);
526}
527
528
529// scattering from a uniform sphere with a Gaussian size distribution
530//
531double
532GaussPolySphere(double dp[], double q)
533{
534        double pi,x;
535        double scale,rad,pd,sig,rho,rhos,bkg,delrho;            //my local names
536        double va,vb,zi,yy,summ,inten;
537        int nord=20,ii;
538       
539        pi = 4.0*atan(1.0);
540        x= q;
541       
542        scale=dp[0];
543        rad=dp[1];
544        pd=dp[2];
545        sig=pd*rad;
546        rho=dp[3];
547        rhos=dp[4];
548        delrho=rho-rhos;
549        bkg=dp[5];
550       
551        va = -4.0*sig + rad;
552        if (va<0) {
553                va=0;           //to avoid numerical error when  va<0 (-ve q-value)
554        }
555        vb = 4.0*sig +rad;
556       
557        summ = 0.0;             // initialize integral
558        for(ii=0;ii<nord;ii+=1) {
559                // calculate Gauss points on integration interval (r-value for evaluation)
560                zi = ( Gauss20Z[ii]*(vb-va) + vb + va )/2.0;
561                // calculate sphere scattering
562                //return(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr)); pass qr
563                yy = F_func(x*zi)*(4.0*pi/3.0*zi*zi*zi)*delrho;
564                yy *= yy;
565                yy *= Gauss20Wt[ii] *  Gauss_distr(sig,rad,zi);
566               
567                summ += yy;             //add to the running total of the quadrature
568        }
569        // calculate value of integral to return
570        inten = (vb-va)/2.0*summ;
571       
572        //re-normalize by polydisperse sphere volume
573        inten /= (4.0*pi/3.0*rad*rad*rad)*(1.0+3.0*pd*pd);
574       
575        inten *= 1.0e8;
576        inten *= scale;
577        inten += bkg;
578       
579    return(inten);      //scale, and add in the background
580}
581
582// scattering from a uniform sphere with a LogNormal size distribution
583//
584double
585LogNormalPolySphere(double dp[], double q)
586{
587        double pi,x;
588        double scale,rad,sig,rho,rhos,bkg,delrho,mu,r3;         //my local names
589        double va,vb,zi,yy,summ,inten;
590        int nord=76,ii;
591       
592        pi = 4.0*atan(1.0);
593        x= q;
594       
595        scale=dp[0];
596        rad=dp[1];              //rad is the median radius
597        mu = log(dp[1]);
598        sig=dp[2];
599        rho=dp[3];
600        rhos=dp[4];
601        delrho=rho-rhos;
602        bkg=dp[5];
603       
604        va = -3.5*sig + mu;
605        va = exp(va);
606        if (va<0) {
607                va=0;           //to avoid numerical error when  va<0 (-ve q-value)
608        }
609        vb = 3.5*sig*(1.0+sig) +mu;
610        vb = exp(vb);
611       
612        summ = 0.0;             // initialize integral
613        for(ii=0;ii<nord;ii+=1) {
614                // calculate Gauss points on integration interval (r-value for evaluation)
615                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
616                // calculate sphere scattering
617                //return(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr)); pass qr
618                yy = F_func(x*zi)*(4.0*pi/3.0*zi*zi*zi)*delrho;
619                yy *= yy;
620                yy *= Gauss76Wt[ii] *  LogNormal_distr(sig,mu,zi);
621               
622                summ += yy;             //add to the running total of the quadrature
623        }
624        // calculate value of integral to return
625        inten = (vb-va)/2.0*summ;
626       
627        //re-normalize by polydisperse sphere volume
628        r3 = exp(3.0*mu + 9.0/2.0*sig*sig);             // <R^3> directly
629        inten /= (4.0*pi/3.0*r3);               //polydisperse volume
630       
631        inten *= 1.0e8;
632        inten *= scale;
633        inten += bkg;
634       
635        return(inten);
636}
637
638static double
639LogNormal_distr(double sig, double mu, double pt)
640{       
641        double retval,pi;
642       
643        pi = 4.0*atan(1.0);
644        retval = (1/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig );
645        return(retval);
646}
647
648static double
649Gauss_distr(double sig, double avg, double pt)
650{       
651        double retval,Pi;
652       
653        Pi = 4.0*atan(1.0);
654        retval = (1.0/ (sig*sqrt(2.0*Pi)) )*exp(-(avg-pt)*(avg-pt)/sig/sig/2.0);
655        return(retval);
656}
657
658// scattering from a core shell sphere with a (Schulz) polydisperse core and constant ratio (shell thickness)/(core radius)
659// - the polydispersity is of the WHOLE sphere
660//
661double
662PolyCoreShellRatio(double dp[], double q)
663{
664        double pi,x;
665        double scale,corrad,thick,shlrad,pp,drho1,drho2,sig,zz,bkg;             //my local names
666        double sld1,sld2,sld3,zp1,zp2,zp3,vpoly;
667        double pi43,c1,c2,form,volume,arg1,arg2;
668       
669        pi = 4.0*atan(1.0);
670        x= q;
671       
672        scale = dp[0];
673        corrad = dp[1];
674        thick = dp[2];
675        sig = dp[3];
676        sld1 = dp[4];
677        sld2 = dp[5];
678        sld3 = dp[6];
679        bkg = dp[7];
680       
681        zz = (1.0/sig)*(1.0/sig) - 1.0;   
682        shlrad = corrad + thick;
683        drho1 = sld1-sld2;              //core-shell
684        drho2 = sld2-sld3;              //shell-solvent
685        zp1 = zz + 1.;
686        zp2 = zz + 2.;
687        zp3 = zz + 3.;
688        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((corrad+thick),3);
689       
690        // the beta factor is not calculated
691        // the calculated form factor <f^2> has units [length^2]
692        // and must be multiplied by number density [l^-3] and the correct unit
693        // conversion to get to absolute scale
694       
695        pi43=4.0/3.0*pi;
696        pp=corrad/shlrad;
697        volume=pi43*shlrad*shlrad*shlrad;
698        c1=drho1*volume;
699        c2=drho2*volume;
700       
701        arg1 = x*shlrad*pp;
702        arg2 = x*shlrad;
703       
704        form=pow(pp,6)*c1*c1*fnt2(arg1,zz);
705        form += c2*c2*fnt2(arg2,zz);
706        form += 2.0*c1*c2*fnt3(arg2,pp,zz);
707       
708        //convert the result to [cm^-1]
709       
710        //scale the result
711        // - divide by the polydisperse volume, mult by 10^8
712        form  /= vpoly;
713        form *= 1.0e8;
714        form *= scale;
715       
716        //add in the background
717        form += bkg;
718       
719        return(form);
720}
721
722//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
723//c
724//c      function fnt2(y,z)
725//c
726double
727fnt2(double yy, double zz)
728{
729        double z1,z2,z3,u,ww,term1,term2,term3,ans;
730       
731        z1=zz+1.0;
732        z2=zz+2.0;
733        z3=zz+3.0;
734        u=yy/z1;
735        ww=atan(2.0*u);
736        term1=cos(z1*ww)/pow((1.0+4.0*u*u),(z1/2.0));
737        term2=2.0*yy*sin(z2*ww)/pow((1.0+4.0*u*u),(z2/2.0));
738        term3=1.0+cos(z3*ww)/pow((1.0+4.0*u*u),(z3/2.0));
739        ans=(4.50/z1/pow(yy,6))*(z1*(1.0-term1-term2)+yy*yy*z2*term3);
740       
741        return(ans);
742}
743
744//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
745//c
746//c      function fnt3(y,p,z)
747//c
748double
749fnt3(double yy, double pp, double zz)
750{       
751        double z1,z2,z3,yp,yn,up,un,vp,vn,term1,term2,term3,term4,term5,term6,ans;
752       
753        z1=zz+1.0;
754        z2=zz+2.0;
755        z3=zz+3.0;
756        yp=(1.0+pp)*yy;
757        yn=(1.0-pp)*yy;
758        up=yp/z1;
759        un=yn/z1;
760        vp=atan(up);
761        vn=atan(un);
762        term1=cos(z1*vn)/pow((1.0+un*un),(z1/2.0));
763        term2=cos(z1*vp)/pow((1.0+up*up),(z1/2.0));
764        term3=cos(z3*vn)/pow((1.0+un*un),(z3/2.0));
765        term4=cos(z3*vp)/pow((1.0+up*up),(z3/2.0));
766        term5=yn*sin(z2*vn)/pow((1.0+un*un),(z2/2.0));
767        term6=yp*sin(z2*vp)/pow((1.0+up*up),(z2/2.0));
768        ans=4.5/z1/pow(yy,6);
769        ans *=(z1*(term1-term2)+yy*yy*pp*z2*(term3+term4)+z1*(term5-term6));
770       
771        return(ans);
772}
773
774// scattering from a a binary population of hard spheres, 3 partial structure factors
775// are properly accounted for...
776//       Input (fitting) variables are:
777//      larger sphere radius(angstroms) = guess[0]
778//      smaller sphere radius (A) = w[1]
779//      number fraction of larger spheres = guess[2]
780//      total volume fraction of spheres = guess[3]
781//      size ratio, alpha(0<a<1) = derived
782//      SLD(A-2) of larger particle = guess[4]
783//      SLD(A-2) of smaller particle = guess[5]
784//      SLD(A-2) of the solvent = guess[6]
785//      background = guess[7]
786double
787BinaryHS(double dp[], double q)
788{
789        double x,pi;
790        double r2,r1,nf2,phi,aa,rho2,rho1,rhos,inten,bgd;               //my local names
791        double psf11,psf12,psf22;
792        double phi1,phi2,phr,a3;
793        double v1,v2,n1,n2,qr1,qr2,b1,b2;
794        int err;
795       
796        pi = 4.0*atan(1.0);
797        x= q;
798        r2 = dp[0];
799        r1 = dp[1];
800        phi2 = dp[2];
801        phi1 = dp[3];
802        rho2 = dp[4];
803        rho1 = dp[5];
804        rhos = dp[6];
805        bgd = dp[7];
806       
807       
808        phi = phi1 + phi2;
809        aa = r1/r2;
810        //calculate the number fraction of larger spheres (eqn 2 in reference)
811        a3=aa*aa*aa;
812        phr=phi2/phi;
813        nf2 = phr*a3/(1.0-phr+phr*a3);
814        // calculate the PSF's here
815        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
816       
817        // /* do form factor calculations  */
818       
819        v1 = 4.0*pi/3.0*r1*r1*r1;
820        v2 = 4.0*pi/3.0*r2*r2*r2;
821       
822        n1 = phi1/v1;
823        n2 = phi2/v2;
824       
825        qr1 = r1*x;
826        qr2 = r2*x;
827       
828        b1 = r1*r1*r1*(rho1-rhos)*4.0*pi*(sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1;
829        b2 = r2*r2*r2*(rho2-rhos)*4.0*pi*(sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2;
830        inten = n1*b1*b1*psf11;
831        inten += sqrt(n1*n2)*2.0*b1*b2*psf12;
832        inten += n2*b2*b2*psf22;
833        ///* convert I(1/A) to (1/cm)  */
834        inten *= 1.0e8;
835       
836        inten += bgd;
837       
838        return(inten);
839}
840
841double
842BinaryHS_PSF11(double dp[], double q)
843{
844        double x,pi;
845        double r2,r1,nf2,phi,aa,rho2,rho1,rhos,bgd;             //my local names
846        double psf11,psf12,psf22;
847        double phi1,phi2,phr,a3;
848        int err;
849       
850        pi = 4.0*atan(1.0);
851        x= q;
852        r2 = dp[0];
853        r1 = dp[1];
854        phi2 = dp[2];
855        phi1 = dp[3];
856        rho2 = dp[4];
857        rho1 = dp[5];
858        rhos = dp[6];
859        bgd = dp[7];
860        phi = phi1 + phi2;
861        aa = r1/r2;
862        //calculate the number fraction of larger spheres (eqn 2 in reference)
863        a3=aa*aa*aa;
864        phr=phi2/phi;
865        nf2 = phr*a3/(1.0-phr+phr*a3);
866        // calculate the PSF's here
867        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
868       
869    return(psf11);      //scale, and add in the background
870}
871
872double
873BinaryHS_PSF12(double dp[], double q)
874{
875        double x,pi;
876        double r2,r1,nf2,phi,aa,rho2,rho1,rhos,bgd;             //my local names
877        double psf11,psf12,psf22;
878        double phi1,phi2,phr,a3;
879        int err;
880       
881        pi = 4.0*atan(1.0);
882        x= q;
883        r2 = dp[0];
884        r1 = dp[1];
885        phi2 = dp[2];
886        phi1 = dp[3];
887        rho2 = dp[4];
888        rho1 = dp[5];
889        rhos = dp[6];
890        bgd = dp[7];
891        phi = phi1 + phi2;
892        aa = r1/r2;
893        //calculate the number fraction of larger spheres (eqn 2 in reference)
894        a3=aa*aa*aa;
895        phr=phi2/phi;
896        nf2 = phr*a3/(1.0-phr+phr*a3);
897        // calculate the PSF's here
898        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
899       
900    return(psf12);      //scale, and add in the background
901}
902
903double
904BinaryHS_PSF22(double dp[], double q)
905{
906        double x,pi;
907        double r2,r1,nf2,phi,aa,rho2,rho1,rhos,bgd;             //my local names
908        double psf11,psf12,psf22;
909        double phi1,phi2,phr,a3;
910        int err;
911       
912        pi = 4.0*atan(1.0);
913        x= q;
914       
915        r2 = dp[0];
916        r1 = dp[1];
917        phi2 = dp[2];
918        phi1 = dp[3];
919        rho2 = dp[4];
920        rho1 = dp[5];
921        rhos = dp[6];
922        bgd = dp[7];
923        phi = phi1 + phi2;
924        aa = r1/r2;
925        //calculate the number fraction of larger spheres (eqn 2 in reference)
926        a3=aa*aa*aa;
927        phr=phi2/phi;
928        nf2 = phr*a3/(1.0-phr+phr*a3);
929        // calculate the PSF's here
930        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
931       
932    return(psf22);      //scale, and add in the background
933}
934
935int
936ashcroft(double qval, double r2, double nf2, double aa, double phi, double *s11, double *s22, double *s12)
937{
938        //      variable qval,r2,nf2,aa,phi,&s11,&s22,&s12
939       
940        //   calculate constant terms
941        double s1,s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4;
942        double a1,a2i,a2,b1,b2,b12,gm1,gm12;
943        double err=0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3;
944        double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13;
945        double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1;
946       
947        s2 = 2.0*r2;
948        s1 = aa*s2;
949        v = phi;
950        a3 = aa*aa*aa;
951        v1=((1.-nf2)*a3/(nf2+(1.-nf2)*a3))*v;
952        v2=(nf2/(nf2+(1.-nf2)*a3))*v;
953        g11=((1.+.5*v)+1.5*v2*(aa-1.))/(1.-v)/(1.-v);
954        g22=((1.+.5*v)+1.5*v1*(1./aa-1.))/(1.-v)/(1.-v);
955        g12=((1.+.5*v)+1.5*(1.-aa)*(v1-v2)/(1.+aa))/(1.-v)/(1.-v);
956        wmv = 1/(1.-v);
957        wmv3 = wmv*wmv*wmv;
958        wmv4 = wmv*wmv3;
959        a1=3.*wmv4*((v1+a3*v2)*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2))) + ((v1+a3*v2)*(1.+2.*v)+(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)-3.*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*wmv3;
960        a2i=((v1+a3*v2)*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*3*wmv4 + ((v1+a3*v2)*(1.+2.*v)+a3*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*aa-3.*v1*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*wmv3;
961        a2=a2i/a3;
962        b1=-6.*(v1*g11*g11+.25*v2*(1.+aa)*(1.+aa)*aa*g12*g12);
963        b2=-6.*(v2*g22*g22+.25*v1/a3*(1.+aa)*(1.+aa)*g12*g12);
964        b12=-3.*aa*(1.+aa)*(v1*g11/aa/aa+v2*g22)*g12;
965        gm1=(v1*a1+a3*v2*a2)*.5;
966        gm12=2.*gm1*(1.-aa)/aa;
967        //c 
968        //c   calculate the direct correlation functions and print results
969        //c
970        //      do 20 j=1,npts
971       
972        yy=qval*s2;
973        //c   calculate direct correlation functions
974        //c   ----c11
975        ay=aa*yy;
976        ay2 = ay*ay;
977        ay3 = ay*ay*ay;
978        t1=a1*(sin(ay)-ay*cos(ay));
979        t2=b1*(2.*ay*sin(ay)-(ay2-2.)*cos(ay)-2.)/ay;
980        t3=gm1*((4.*ay*ay2-24.*ay)*sin(ay)-(ay2*ay2-12.*ay2+24.)*cos(ay)+24.)/ay3;
981        f11=24.*v1*(t1+t2+t3)/ay3;
982       
983        //c ------c22
984        y2=yy*yy;
985        y3=yy*y2;
986        tt1=a2*(sin(yy)-yy*cos(yy));
987        tt2=b2*(2.*yy*sin(yy)-(y2-2.)*cos(yy)-2.)/yy;
988        tt3=gm1*((4.*y3-24.*yy)*sin(yy)-(y2*y2-12.*y2+24.)*cos(yy)+24.)/ay3;
989        f22=24.*v2*(tt1+tt2+tt3)/y3;
990       
991        //c   -----c12
992        yl=.5*yy*(1.-aa);
993        yl3=yl*yl*yl;
994        wma3 = (1.-aa)*(1.-aa)*(1.-aa);
995        y1=aa*yy;
996        y13 = y1*y1*y1;
997        ttt1=3.*wma3*v*sqrt(nf2)*sqrt(1.-nf2)*a1*(sin(yl)-yl*cos(yl))/((nf2+(1.-nf2)*a3)*yl3);
998        t21=b12*(2.*y1*cos(y1)+(y1*y1-2.)*sin(y1));
999        t22=gm12*((3.*y1*y1-6.)*cos(y1)+(y1*y1*y1-6.*y1)*sin(y1)+6.)/y1;
1000        t23=gm1*((4.*y13-24.*y1)*cos(y1)+(y13*y1-12.*y1*y1+24.)*sin(y1))/(y1*y1);
1001        t31=b12*(2.*y1*sin(y1)-(y1*y1-2.)*cos(y1)-2.);
1002        t32=gm12*((3.*y1*y1-6.)*sin(y1)-(y1*y1*y1-6.*y1)*cos(y1))/y1;
1003        t33=gm1*((4.*y13-24.*y1)*sin(y1)-(y13*y1-12.*y1*y1+24.)*cos(y1)+24.)/(y1*y1);
1004        t41=cos(yl)*((sin(y1)-y1*cos(y1))/(y1*y1) + (1.-aa)/(2.*aa)*(1.-cos(y1))/y1);
1005        t42=sin(yl)*((cos(y1)+y1*sin(y1)-1.)/(y1*y1) + (1.-aa)/(2.*aa)*sin(y1)/y1);
1006        ttt2=sin(yl)*(t21+t22+t23)/(y13*y1);
1007        ttt3=cos(yl)*(t31+t32+t33)/(y13*y1);
1008        ttt4=a1*(t41+t42)/y1;
1009        f12=ttt1+24.*v*sqrt(nf2)*sqrt(1.-nf2)*a3*(ttt2+ttt3+ttt4)/(nf2+(1.-nf2)*a3);
1010       
1011        c11=f11;
1012        c22=f22;
1013        c12=f12;
1014        *s11=1./(1.+c11-(c12)*c12/(1.+c22));
1015        *s22=1./(1.+c22-(c12)*c12/(1.+c11)); 
1016        *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));   
1017       
1018        return(err);
1019}
1020
1021
1022
1023/*
1024 // calculates the scattering from a spherical particle made up of a core (aqueous) surrounded
1025 // by N spherical layers, each of which is a PAIR of shells, solvent + surfactant since there
1026 //must always be a surfactant layer on the outside
1027 //
1028 // bragg peaks arise naturally from the periodicity of the sample
1029 // resolution smeared version gives he most appropriate view of the model
1030 
1031        Warning:
1032 The call to WaveData() below returns a pointer to the middle
1033 of an unlocked Macintosh handle. In the unlikely event that your
1034 calculations could cause memory to move, you should copy the coefficient
1035 values to local variables or an array before such operations.
1036 */
1037double
1038MultiShell(double dp[], double q)
1039{
1040        double x;
1041        double scale,rcore,tw,ts,rhocore,rhoshel,num,bkg;               //my local names
1042        int ii;
1043        double fval,voli,ri,sldi;
1044        double pi;
1045       
1046        pi = 4.0*atan(1.0);
1047       
1048        x= q;
1049        scale = dp[0];
1050        rcore = dp[1];
1051        ts = dp[2];
1052        tw = dp[3];
1053        rhocore = dp[4];
1054        rhoshel = dp[5];
1055        num = dp[6];
1056        bkg = dp[7];
1057       
1058        //calculate with a loop, two shells at a time
1059       
1060        ii=0;
1061        fval=0;
1062       
1063        do {
1064                ri = rcore + (double)ii*ts + (double)ii*tw;
1065                voli = 4*pi/3*ri*ri*ri;
1066                sldi = rhocore-rhoshel;
1067                fval += voli*sldi*F_func(ri*x);
1068                ri += ts;
1069                voli = 4*pi/3*ri*ri*ri;
1070                sldi = rhoshel-rhocore;
1071                fval += voli*sldi*F_func(ri*x);
1072                ii+=1;          //do 2 layers at a time
1073        } while(ii<=num-1);  //change to make 0 < num < 2 correspond to unilamellar vesicles (C. Glinka, 11/24/03)
1074       
1075        fval *= fval;           //square it
1076        fval /= voli;           //normalize by the overall volume
1077        fval *= scale*1e8;
1078        fval += bkg;
1079       
1080        return(fval);
1081}
1082
1083/*
1084 // calculates the scattering from a POLYDISPERSE spherical particle made up of a core (aqueous) surrounded
1085 // by N spherical layers, each of which is a PAIR of shells, solvent + surfactant since there
1086 //must always be a surfactant layer on the outside
1087 //
1088 // bragg peaks arise naturally from the periodicity of the sample
1089 // resolution smeared version gives he most appropriate view of the model
1090 //
1091 // Polydispersity is of the total (outer) radius. This is converted into a distribution of MLV's
1092 // with integer numbers of layers, with a minimum of one layer... a vesicle... depending
1093 // on the parameters, the "distribution" of MLV's that is used may be truncated
1094 //
1095        Warning:
1096 The call to WaveData() below returns a pointer to the middle
1097 of an unlocked Macintosh handle. In the unlikely event that your
1098 calculations could cause memory to move, you should copy the coefficient
1099 values to local variables or an array before such operations.
1100 */
1101double
1102PolyMultiShell(double dp[], double q)
1103{
1104        double x;
1105        double scale,rcore,tw,ts,rhocore,rhoshel,bkg;           //my local names
1106        int ii,minPairs,maxPairs,first;
1107        double fval,ri,pi;
1108        double avg,pd,zz,lo,hi,r1,r2,d1,d2,distr;
1109       
1110        pi = 4.0*atan(1.0);     
1111        x= q;
1112       
1113        scale = dp[0];
1114        avg = dp[1];            // average (total) outer radius
1115        pd = dp[2];
1116        rcore = dp[3];
1117        ts = dp[4];
1118        tw = dp[5];
1119        rhocore = dp[6];
1120        rhoshel = dp[7];
1121        bkg = dp[8];
1122       
1123        zz = (1.0/pd)*(1.0/pd)-1.0;
1124       
1125        //max radius set to be 5 std deviations past mean
1126        hi = avg + pd*avg*5.0;
1127        lo = avg - pd*avg*5.0;
1128       
1129        maxPairs = trunc( (hi-rcore+tw)/(ts+tw) );
1130        minPairs = trunc( (lo-rcore+tw)/(ts+tw) );
1131        minPairs = (minPairs < 1) ? 1 : minPairs;       // need a minimum of one
1132       
1133        ii=minPairs;
1134        fval=0;
1135        d1 = 0;
1136        d2 = 0;
1137        r1 = 0;
1138        r2 = 0;
1139        distr = 0;
1140        first = 1;
1141        do {
1142                //make the current values old
1143                r1 = r2;
1144                d1 = d2;
1145               
1146                ri = (double)ii*(ts+tw) - tw + rcore;
1147                fval += SchulzPoint(ri,avg,zz) * MultiShellGuts(x, rcore, ts, tw, rhocore, rhoshel, ii) * (4*pi/3*ri*ri*ri);
1148                // get a running integration of the fraction of the distribution used, but not the first time
1149                r2 = ri;
1150                d2 = SchulzPoint(ri,avg,zz);
1151                if( !first ) {
1152                        distr += 0.5*(d1+d2)*(r2-r1);           //cheap trapezoidal integration
1153                }
1154                ii+=1;
1155                first = 0;
1156        } while(ii<=maxPairs);
1157       
1158        fval /= 4*pi/3*avg*avg*avg;             //normalize by the overall volume
1159        fval /= distr;
1160        fval *= scale;
1161        fval += bkg;
1162       
1163        return(fval);
1164}
1165
1166double
1167MultiShellGuts(double x,double rcore,double ts,double tw,double rhocore,double rhoshel,int num) {
1168       
1169    double ri,sldi,fval,voli,pi;
1170    int ii;
1171   
1172        pi = 4.0*atan(1.0);
1173    ii=0;
1174    fval=0;
1175   
1176    do {
1177        ri = rcore + (double)ii*ts + (double)ii*tw;
1178        voli = 4*pi/3*ri*ri*ri;
1179        sldi = rhocore-rhoshel;
1180        fval += voli*sldi*F_func(ri*x);
1181        ri += ts;
1182        voli = 4*pi/3*ri*ri*ri;
1183        sldi = rhoshel-rhocore;
1184        fval += voli*sldi*F_func(ri*x);
1185        ii+=1;          //do 2 layers at a time
1186    } while(ii<=num-1);  //change to make 0 < num < 2 correspond to unilamellar vesicles (C. Glinka, 11/24/03)
1187   
1188    fval *= fval;
1189    fval /= voli;
1190    fval *= 1e8;
1191   
1192    return(fval);       // this result still needs to be multiplied by scale and have background added
1193}
1194
1195static double
1196SchulzPoint(double x, double avg, double zz) {
1197       
1198    double dr;
1199   
1200    dr = zz*log(x) - gammln(zz+1)+(zz+1)*log((zz+1)/avg)-(x/avg*(zz+1));
1201    return (exp(dr));
1202}
1203
1204static double
1205gammln(double xx) {
1206       
1207    double x,y,tmp,ser;
1208    static double cof[6]={76.18009172947146,-86.50532032941677,
1209                24.01409824083091,-1.231739572450155,
1210                0.1208650973866179e-2,-0.5395239384953e-5};
1211    int j;
1212       
1213    y=x=xx;
1214    tmp=x+5.5;
1215    tmp -= (x+0.5)*log(tmp);
1216    ser=1.000000000190015;
1217    for (j=0;j<=5;j++) ser += cof[j]/++y;
1218    return -tmp+log(2.5066282746310005*ser/x);
1219}
1220
1221double
1222F_func(double qr) {
1223        return(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr));
1224}
1225
Note: See TracBrowser for help on using the repository browser.