source: sasview/sansmodels/src/libigor/libSphere.c @ b9405cd

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 b9405cd was 6e93a02, checked in by Jae Cho <jhjcho@…>, 15 years ago

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