source: sasview/sansmodels/src/libigor/libSphere.c @ 72a05d0

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 72a05d0 was 975ec8e, checked in by Jae Cho <jhjcho@…>, 15 years ago

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