source: sasview/sansmodels/src/libigor/libSphere.c @ 5ee2306

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 5ee2306 was ae3ce4e, checked in by Mathieu Doucet <doucetm@…>, 17 years ago

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