source: sasview/sansmodels/src/sans/models/libigor/libSphere.c @ 2605da22

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 2605da22 was 34c2649, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Re #4 Fixed warnings
  • Property mode set to 100644
File size: 57.6 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
655/*
656static double
657LogNormal_distr(double sig, double mu, double pt)
658{       
659        double retval,pi;
660       
661        pi = 4.0*atan(1.0);
662        retval = (1.0/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig );
663        return(retval);
664}
665
666static double
667Gauss_distr(double sig, double avg, double pt)
668{       
669        double retval,Pi;
670       
671        Pi = 4.0*atan(1.0);
672        retval = (1.0/ (sig*sqrt(2.0*Pi)) )*exp(-(avg-pt)*(avg-pt)/sig/sig/2.0);
673        return(retval);
674}
675*/
676
677// scattering from a core shell sphere with a (Schulz) polydisperse core and constant ratio (shell thickness)/(core radius)
678// - the polydispersity is of the WHOLE sphere
679//
680double
681PolyCoreShellRatio(double dp[], double q)
682{
683        double pi,x;
684        double scale,corrad,thick,shlrad,pp,drho1,drho2,sig,zz,bkg;             //my local names
685        double sld1,sld2,sld3,zp1,zp2,zp3,vpoly;
686        double pi43,c1,c2,form,volume,arg1,arg2;
687       
688        pi = 4.0*atan(1.0);
689        x= q;
690       
691        scale = dp[0];
692        corrad = dp[1];
693        thick = dp[2];
694        sig = dp[3];
695        sld1 = dp[4];
696        sld2 = dp[5];
697        sld3 = dp[6];
698        bkg = dp[7];
699       
700        zz = (1.0/sig)*(1.0/sig) - 1.0;   
701        shlrad = corrad + thick;
702        drho1 = sld1-sld2;              //core-shell
703        drho2 = sld2-sld3;              //shell-solvent
704        zp1 = zz + 1.;
705        zp2 = zz + 2.;
706        zp3 = zz + 3.;
707        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((corrad+thick),3);
708       
709        // the beta factor is not calculated
710        // the calculated form factor <f^2> has units [length^2]
711        // and must be multiplied by number density [l^-3] and the correct unit
712        // conversion to get to absolute scale
713       
714        pi43=4.0/3.0*pi;
715        pp=corrad/shlrad;
716        volume=pi43*shlrad*shlrad*shlrad;
717        c1=drho1*volume;
718        c2=drho2*volume;
719       
720        arg1 = x*shlrad*pp;
721        arg2 = x*shlrad;
722       
723        form=pow(pp,6)*c1*c1*fnt2(arg1,zz);
724        form += c2*c2*fnt2(arg2,zz);
725        form += 2.0*c1*c2*fnt3(arg2,pp,zz);
726       
727        //convert the result to [cm^-1]
728       
729        //scale the result
730        // - divide by the polydisperse volume, mult by 10^8
731        form  /= vpoly;
732        form *= 1.0e8;
733        form *= scale;
734       
735        //add in the background
736        form += bkg;
737       
738        return(form);
739}
740
741//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
742//c
743//c      function fnt2(y,z)
744//c
745double
746fnt2(double yy, double zz)
747{
748        double z1,z2,z3,u,ww,term1,term2,term3,ans;
749       
750        z1=zz+1.0;
751        z2=zz+2.0;
752        z3=zz+3.0;
753        u=yy/z1;
754        ww=atan(2.0*u);
755        term1=cos(z1*ww)/pow((1.0+4.0*u*u),(z1/2.0));
756        term2=2.0*yy*sin(z2*ww)/pow((1.0+4.0*u*u),(z2/2.0));
757        term3=1.0+cos(z3*ww)/pow((1.0+4.0*u*u),(z3/2.0));
758        ans=(4.50/z1/pow(yy,6))*(z1*(1.0-term1-term2)+yy*yy*z2*term3);
759       
760        return(ans);
761}
762
763//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
764//c
765//c      function fnt3(y,p,z)
766//c
767double
768fnt3(double yy, double pp, double zz)
769{       
770        double z1,z2,z3,yp,yn,up,un,vp,vn,term1,term2,term3,term4,term5,term6,ans;
771       
772        z1=zz+1.0;
773        z2=zz+2.0;
774        z3=zz+3.0;
775        yp=(1.0+pp)*yy;
776        yn=(1.0-pp)*yy;
777        up=yp/z1;
778        un=yn/z1;
779        vp=atan(up);
780        vn=atan(un);
781        term1=cos(z1*vn)/pow((1.0+un*un),(z1/2.0));
782        term2=cos(z1*vp)/pow((1.0+up*up),(z1/2.0));
783        term3=cos(z3*vn)/pow((1.0+un*un),(z3/2.0));
784        term4=cos(z3*vp)/pow((1.0+up*up),(z3/2.0));
785        term5=yn*sin(z2*vn)/pow((1.0+un*un),(z2/2.0));
786        term6=yp*sin(z2*vp)/pow((1.0+up*up),(z2/2.0));
787        ans=4.5/z1/pow(yy,6);
788        ans *=(z1*(term1-term2)+yy*yy*pp*z2*(term3+term4)+z1*(term5-term6));
789       
790        return(ans);
791}
792
793// scattering from a a binary population of hard spheres, 3 partial structure factors
794// are properly accounted for...
795//       Input (fitting) variables are:
796//      larger sphere radius(angstroms) = guess[0]
797//      smaller sphere radius (A) = w[1]
798//      number fraction of larger spheres = guess[2]
799//      total volume fraction of spheres = guess[3]
800//      size ratio, alpha(0<a<1) = derived
801//      SLD(A-2) of larger particle = guess[4]
802//      SLD(A-2) of smaller particle = guess[5]
803//      SLD(A-2) of the solvent = guess[6]
804//      background = guess[7]
805double
806BinaryHS(double dp[], double q)
807{
808        double x,pi;
809        double r2,r1,nf2,phi,aa,rho2,rho1,rhos,inten,bgd;               //my local names
810        double psf11,psf12,psf22;
811        double phi1,phi2,phr,a3;
812        double v1,v2,n1,n2,qr1,qr2,b1,b2,sc1,sc2;
813        int err;
814       
815        pi = 4.0*atan(1.0);
816        x= q;
817        r2 = dp[0];
818        r1 = dp[1];
819        phi2 = dp[2];
820        phi1 = dp[3];
821        rho2 = dp[4];
822        rho1 = dp[5];
823        rhos = dp[6];
824        bgd = dp[7];
825       
826       
827        phi = phi1 + phi2;
828        aa = r1/r2;
829        //calculate the number fraction of larger spheres (eqn 2 in reference)
830        a3=aa*aa*aa;
831        phr=phi2/phi;
832        nf2 = phr*a3/(1.0-phr+phr*a3);
833        // calculate the PSF's here
834        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
835       
836        // /* do form factor calculations  */
837       
838        v1 = 4.0*pi/3.0*r1*r1*r1;
839        v2 = 4.0*pi/3.0*r2*r2*r2;
840       
841        n1 = phi1/v1;
842        n2 = phi2/v2;
843       
844        qr1 = r1*x;
845        qr2 = r2*x;
846
847        if (qr1 == 0){
848                sc1 = 1.0/3.0;
849        }else{
850                sc1 = (sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1;
851        }
852        if (qr2 == 0){
853                sc2 = 1.0/3.0;
854        }else{
855                sc2 = (sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2;
856        }
857        b1 = r1*r1*r1*(rho1-rhos)*4.0*pi*sc1;
858        b2 = r2*r2*r2*(rho2-rhos)*4.0*pi*sc2;
859        inten = n1*b1*b1*psf11;
860        inten += sqrt(n1*n2)*2.0*b1*b2*psf12;
861        inten += n2*b2*b2*psf22;
862        ///* convert I(1/A) to (1/cm)  */
863        inten *= 1.0e8;
864       
865        inten += bgd;
866       
867        return(inten);
868}
869
870double
871BinaryHS_PSF11(double dp[], double q)
872{
873        double x,pi;
874        double r2,r1,nf2,phi,aa,rho2,rho1,rhos,bgd;             //my local names
875        double psf11,psf12,psf22;
876        double phi1,phi2,phr,a3;
877        int err;
878       
879        pi = 4.0*atan(1.0);
880        x= q;
881        r2 = dp[0];
882        r1 = dp[1];
883        phi2 = dp[2];
884        phi1 = dp[3];
885        rho2 = dp[4];
886        rho1 = dp[5];
887        rhos = dp[6];
888        bgd = dp[7];
889        phi = phi1 + phi2;
890        aa = r1/r2;
891        //calculate the number fraction of larger spheres (eqn 2 in reference)
892        a3=aa*aa*aa;
893        phr=phi2/phi;
894        nf2 = phr*a3/(1.0-phr+phr*a3);
895        // calculate the PSF's here
896        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
897       
898    return(psf11);      //scale, and add in the background
899}
900
901double
902BinaryHS_PSF12(double dp[], double q)
903{
904        double x,pi;
905        double r2,r1,nf2,phi,aa,rho2,rho1,rhos,bgd;             //my local names
906        double psf11,psf12,psf22;
907        double phi1,phi2,phr,a3;
908        int err;
909       
910        pi = 4.0*atan(1.0);
911        x= q;
912        r2 = dp[0];
913        r1 = dp[1];
914        phi2 = dp[2];
915        phi1 = dp[3];
916        rho2 = dp[4];
917        rho1 = dp[5];
918        rhos = dp[6];
919        bgd = dp[7];
920        phi = phi1 + phi2;
921        aa = r1/r2;
922        //calculate the number fraction of larger spheres (eqn 2 in reference)
923        a3=aa*aa*aa;
924        phr=phi2/phi;
925        nf2 = phr*a3/(1.0-phr+phr*a3);
926        // calculate the PSF's here
927        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
928       
929    return(psf12);      //scale, and add in the background
930}
931
932double
933BinaryHS_PSF22(double dp[], double q)
934{
935        double x,pi;
936        double r2,r1,nf2,phi,aa,rho2,rho1,rhos,bgd;             //my local names
937        double psf11,psf12,psf22;
938        double phi1,phi2,phr,a3;
939        int err;
940       
941        pi = 4.0*atan(1.0);
942        x= q;
943       
944        r2 = dp[0];
945        r1 = dp[1];
946        phi2 = dp[2];
947        phi1 = dp[3];
948        rho2 = dp[4];
949        rho1 = dp[5];
950        rhos = dp[6];
951        bgd = dp[7];
952        phi = phi1 + phi2;
953        aa = r1/r2;
954        //calculate the number fraction of larger spheres (eqn 2 in reference)
955        a3=aa*aa*aa;
956        phr=phi2/phi;
957        nf2 = phr*a3/(1.0-phr+phr*a3);
958        // calculate the PSF's here
959        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
960       
961    return(psf22);      //scale, and add in the background
962}
963
964int
965ashcroft(double qval, double r2, double nf2, double aa, double phi, double *s11, double *s22, double *s12)
966{
967        //      variable qval,r2,nf2,aa,phi,&s11,&s22,&s12
968       
969        //   calculate constant terms
970        double s1,s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4;
971        double a1,a2i,a2,b1,b2,b12,gm1,gm12;
972        double err=0.0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3;
973        double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13;
974        double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1;
975       
976        s2 = 2.0*r2;
977        s1 = aa*s2;
978        v = phi;
979        a3 = aa*aa*aa;
980        v1=((1.-nf2)*a3/(nf2+(1.-nf2)*a3))*v;
981        v2=(nf2/(nf2+(1.-nf2)*a3))*v;
982        g11=((1.+.5*v)+1.5*v2*(aa-1.))/(1.-v)/(1.-v);
983        g22=((1.+.5*v)+1.5*v1*(1./aa-1.))/(1.-v)/(1.-v);
984        g12=((1.+.5*v)+1.5*(1.-aa)*(v1-v2)/(1.+aa))/(1.-v)/(1.-v);
985        wmv = 1/(1.-v);
986        wmv3 = wmv*wmv*wmv;
987        wmv4 = wmv*wmv3;
988        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;
989        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;
990        a2=a2i/a3;
991        b1=-6.*(v1*g11*g11+.25*v2*(1.+aa)*(1.+aa)*aa*g12*g12);
992        b2=-6.*(v2*g22*g22+.25*v1/a3*(1.+aa)*(1.+aa)*g12*g12);
993        b12=-3.*aa*(1.+aa)*(v1*g11/aa/aa+v2*g22)*g12;
994        gm1=(v1*a1+a3*v2*a2)*.5;
995        gm12=2.*gm1*(1.-aa)/aa;
996        //c 
997        //c   calculate the direct correlation functions and print results
998        //c
999        //      do 20 j=1,npts
1000       
1001        yy=qval*s2;
1002        //c   calculate direct correlation functions
1003        //c   ----c11
1004        ay=aa*yy;
1005        ay2 = ay*ay;
1006        ay3 = ay*ay*ay;
1007        t1=a1*(sin(ay)-ay*cos(ay));
1008        t2=b1*(2.*ay*sin(ay)-(ay2-2.)*cos(ay)-2.)/ay;
1009        t3=gm1*((4.*ay*ay2-24.*ay)*sin(ay)-(ay2*ay2-12.*ay2+24.)*cos(ay)+24.)/ay3;
1010        f11=24.*v1*(t1+t2+t3)/ay3;
1011       
1012        //c ------c22
1013        y2=yy*yy;
1014        y3=yy*y2;
1015        tt1=a2*(sin(yy)-yy*cos(yy));
1016        tt2=b2*(2.*yy*sin(yy)-(y2-2.)*cos(yy)-2.)/yy;
1017        tt3=gm1*((4.*y3-24.*yy)*sin(yy)-(y2*y2-12.*y2+24.)*cos(yy)+24.)/ay3;
1018        f22=24.*v2*(tt1+tt2+tt3)/y3;
1019       
1020        //c   -----c12
1021        yl=.5*yy*(1.-aa);
1022        yl3=yl*yl*yl;
1023        wma3 = (1.-aa)*(1.-aa)*(1.-aa);
1024        y1=aa*yy;
1025        y13 = y1*y1*y1;
1026        ttt1=3.*wma3*v*sqrt(nf2)*sqrt(1.-nf2)*a1*(sin(yl)-yl*cos(yl))/((nf2+(1.-nf2)*a3)*yl3);
1027        t21=b12*(2.*y1*cos(y1)+(y1*y1-2.)*sin(y1));
1028        t22=gm12*((3.*y1*y1-6.)*cos(y1)+(y1*y1*y1-6.*y1)*sin(y1)+6.)/y1;
1029        t23=gm1*((4.*y13-24.*y1)*cos(y1)+(y13*y1-12.*y1*y1+24.)*sin(y1))/(y1*y1);
1030        t31=b12*(2.*y1*sin(y1)-(y1*y1-2.)*cos(y1)-2.);
1031        t32=gm12*((3.*y1*y1-6.)*sin(y1)-(y1*y1*y1-6.*y1)*cos(y1))/y1;
1032        t33=gm1*((4.*y13-24.*y1)*sin(y1)-(y13*y1-12.*y1*y1+24.)*cos(y1)+24.)/(y1*y1);
1033        t41=cos(yl)*((sin(y1)-y1*cos(y1))/(y1*y1) + (1.-aa)/(2.*aa)*(1.-cos(y1))/y1);
1034        t42=sin(yl)*((cos(y1)+y1*sin(y1)-1.)/(y1*y1) + (1.-aa)/(2.*aa)*sin(y1)/y1);
1035        ttt2=sin(yl)*(t21+t22+t23)/(y13*y1);
1036        ttt3=cos(yl)*(t31+t32+t33)/(y13*y1);
1037        ttt4=a1*(t41+t42)/y1;
1038        f12=ttt1+24.*v*sqrt(nf2)*sqrt(1.-nf2)*a3*(ttt2+ttt3+ttt4)/(nf2+(1.-nf2)*a3);
1039       
1040        c11=f11;
1041        c22=f22;
1042        c12=f12;
1043        *s11=1./(1.+c11-(c12)*c12/(1.+c22));
1044        *s22=1./(1.+c22-(c12)*c12/(1.+c11)); 
1045        *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));   
1046       
1047        return(err);
1048}
1049
1050
1051
1052/*
1053 // calculates the scattering from a spherical particle made up of a core (aqueous) surrounded
1054 // by N spherical layers, each of which is a PAIR of shells, solvent + surfactant since there
1055 //must always be a surfactant layer on the outside
1056 //
1057 // bragg peaks arise naturally from the periodicity of the sample
1058 // resolution smeared version gives he most appropriate view of the model
1059 
1060        Warning:
1061 The call to WaveData() below returns a pointer to the middle
1062 of an unlocked Macintosh handle. In the unlikely event that your
1063 calculations could cause memory to move, you should copy the coefficient
1064 values to local variables or an array before such operations.
1065 */
1066double
1067MultiShell(double dp[], double q)
1068{
1069        double x;
1070        double scale,rcore,tw,ts,rhocore,rhoshel,num,bkg;               //my local names
1071        int ii;
1072        double fval,voli,ri,sldi;
1073        double pi;
1074       
1075        pi = 4.0*atan(1.0);
1076       
1077        x= q;
1078        scale = dp[0];
1079        rcore = dp[1];
1080        ts = dp[2];
1081        tw = dp[3];
1082        rhocore = dp[4];
1083        rhoshel = dp[5];
1084        num = dp[6];
1085        bkg = dp[7];
1086       
1087        //calculate with a loop, two shells at a time
1088       
1089        ii=0;
1090        fval=0.0;
1091       
1092        do {
1093                ri = rcore + (double)ii*ts + (double)ii*tw;
1094                voli = 4.0*pi/3.0*ri*ri*ri;
1095                sldi = rhocore-rhoshel;
1096                fval += voli*sldi*F_func(ri*x);
1097                ri += ts;
1098                voli = 4.0*pi/3.0*ri*ri*ri;
1099                sldi = rhoshel-rhocore;
1100                fval += voli*sldi*F_func(ri*x);
1101                ii+=1;          //do 2 layers at a time
1102        } while(ii<=num-1);  //change to make 0 < num < 2 correspond to unilamellar vesicles (C. Glinka, 11/24/03)
1103       
1104        fval *= fval;           //square it
1105        fval /= voli;           //normalize by the overall volume
1106        fval *= scale*1.0e8;
1107        fval += bkg;
1108       
1109        return(fval);
1110}
1111
1112/*
1113 // calculates the scattering from a POLYDISPERSE spherical particle made up of a core (aqueous) surrounded
1114 // by N spherical layers, each of which is a PAIR of shells, solvent + surfactant since there
1115 //must always be a surfactant layer on the outside
1116 //
1117 // bragg peaks arise naturally from the periodicity of the sample
1118 // resolution smeared version gives he most appropriate view of the model
1119 //
1120 // Polydispersity is of the total (outer) radius. This is converted into a distribution of MLV's
1121 // with integer numbers of layers, with a minimum of one layer... a vesicle... depending
1122 // on the parameters, the "distribution" of MLV's that is used may be truncated
1123 //
1124        Warning:
1125 The call to WaveData() below returns a pointer to the middle
1126 of an unlocked Macintosh handle. In the unlikely event that your
1127 calculations could cause memory to move, you should copy the coefficient
1128 values to local variables or an array before such operations.
1129 */
1130double
1131PolyMultiShell(double dp[], double q)
1132{
1133        double x;
1134        double scale,rcore,tw,ts,rhocore,rhoshel,bkg;           //my local names
1135        int ii,minPairs,maxPairs,first;
1136        double fval,ri,pi;
1137        double avg,pd,zz,lo,hi,r1,r2,d1,d2,distr;
1138       
1139        pi = 4.0*atan(1.0);     
1140        x= q;
1141       
1142        scale = dp[0];
1143        avg = dp[1];            // average (total) outer radius
1144        pd = dp[2];
1145        rcore = dp[3];
1146        ts = dp[4];
1147        tw = dp[5];
1148        rhocore = dp[6];
1149        rhoshel = dp[7];
1150        bkg = dp[8];
1151       
1152        zz = (1.0/pd)*(1.0/pd)-1.0;
1153       
1154        //max radius set to be 5 std deviations past mean
1155        hi = avg + pd*avg*5.0;
1156        lo = avg - pd*avg*5.0;
1157       
1158        maxPairs = trunc( (hi-rcore+tw)/(ts+tw) );
1159        minPairs = trunc( (lo-rcore+tw)/(ts+tw) );
1160        minPairs = (minPairs < 1) ? 1 : minPairs;       // need a minimum of one
1161       
1162        ii=minPairs;
1163        fval=0.0;
1164        d1 = 0.0;
1165        d2 = 0.0;
1166        r1 = 0.0;
1167        r2 = 0.0;
1168        distr = 0.0;
1169        first = 1.0;
1170        do {
1171                //make the current values old
1172                r1 = r2;
1173                d1 = d2;
1174               
1175                ri = (double)ii*(ts+tw) - tw + rcore;
1176                fval += SchulzPoint(ri,avg,zz) * MultiShellGuts(x, rcore, ts, tw, rhocore, rhoshel, ii) * (4*pi/3*ri*ri*ri);
1177                // get a running integration of the fraction of the distribution used, but not the first time
1178                r2 = ri;
1179                d2 = SchulzPoint(ri,avg,zz);
1180                if( !first ) {
1181                        distr += 0.5*(d1+d2)*(r2-r1);           //cheap trapezoidal integration
1182                }
1183                ii+=1;
1184                first = 0;
1185        } while(ii<=maxPairs);
1186       
1187        fval /= 4.0*pi/3.0*avg*avg*avg;         //normalize by the overall volume
1188        fval /= distr;
1189        fval *= scale;
1190        fval += bkg;
1191       
1192        return(fval);
1193}
1194
1195double
1196MultiShellGuts(double x,double rcore,double ts,double tw,double rhocore,double rhoshel,int num) {
1197       
1198    double ri,sldi,fval,voli,pi;
1199    int ii;
1200   
1201        pi = 4.0*atan(1.0);
1202    ii=0;
1203    fval=0.0;
1204   
1205    do {
1206        ri = rcore + (double)ii*ts + (double)ii*tw;
1207        voli = 4.0*pi/3.0*ri*ri*ri;
1208        sldi = rhocore-rhoshel;
1209        fval += voli*sldi*F_func(ri*x);
1210        ri += ts;
1211        voli = 4.0*pi/3.0*ri*ri*ri;
1212        sldi = rhoshel-rhocore;
1213        fval += voli*sldi*F_func(ri*x);
1214        ii+=1;          //do 2 layers at a time
1215    } while(ii<=num-1);  //change to make 0 < num < 2 correspond to unilamellar vesicles (C. Glinka, 11/24/03)
1216   
1217    fval *= fval;
1218    fval /= voli;
1219    fval *= 1.0e8;
1220   
1221    return(fval);       // this result still needs to be multiplied by scale and have background added
1222}
1223
1224/*
1225static double
1226SchulzPoint(double x, double avg, double zz) {
1227       
1228    double dr;
1229   
1230    dr = zz*log(x) - gammln(zz+1.0)+(zz+1.0)*log((zz+1.0)/avg)-(x/avg*(zz+1.0));
1231    return (exp(dr));
1232}
1233
1234static double
1235gammln(double xx) {
1236       
1237    double x,y,tmp,ser;
1238    static double cof[6]={76.18009172947146,-86.50532032941677,
1239                24.01409824083091,-1.231739572450155,
1240                0.1208650973866179e-2,-0.5395239384953e-5};
1241    int j;
1242       
1243    y=x=xx;
1244    tmp=x+5.5;
1245    tmp -= (x+0.5)*log(tmp);
1246    ser=1.000000000190015;
1247    for (j=0;j<=5;j++) ser += cof[j]/++y;
1248    return -tmp+log(2.5066282746310005*ser/x);
1249}
1250*/
1251
1252double
1253F_func(double qr) {
1254        double sc;
1255        if (qr == 0.0){
1256                sc = 1.0;
1257        }else{
1258                sc=(3.0*(sin(qr) - qr*cos(qr))/(qr*qr*qr));
1259        }
1260        return sc;
1261}
1262
1263double
1264OneShell(double dp[], double q)
1265{
1266        // variables are:
1267        //[0] scale factor
1268        //[1] radius of core [ï¿œ]
1269        //[2] SLD of the core   [ï¿œ-2]
1270        //[3] thickness of the shell    [ï¿œ]
1271        //[4] SLD of the shell
1272        //[5] SLD of the solvent
1273        //[6] background        [cm-1]
1274       
1275        double x,pi;
1276        double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg;           //my local names
1277        double bes,f,vol,qr,contr,f2;
1278       
1279        pi = 4.0*atan(1.0);
1280        x=q;
1281       
1282        scale = dp[0];
1283        rcore = dp[1];
1284        rhocore = dp[2];
1285        thick = dp[3];
1286        rhoshel = dp[4];
1287        rhosolv = dp[5];
1288        bkg = dp[6];
1289       
1290        // core first, then add in shell
1291        qr=x*rcore;
1292        contr = rhocore-rhoshel;
1293        if(qr == 0){
1294                bes = 1.0;
1295        }else{
1296                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1297        }
1298        vol = 4.0*pi/3.0*rcore*rcore*rcore;
1299        f = vol*bes*contr;
1300        //now the shell
1301        qr=x*(rcore+thick);
1302        contr = rhoshel-rhosolv;
1303        if(qr == 0){
1304                bes = 1.0;
1305        }else{
1306                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1307        }
1308        vol = 4.0*pi/3.0*pow((rcore+thick),3);
1309        f += vol*bes*contr;
1310       
1311        // normalize to particle volume and rescale from [ï¿œ-1] to [cm-1]
1312        f2 = f*f/vol*1.0e8;
1313       
1314        //scale if desired
1315        f2 *= scale;
1316        // then add in the background
1317        f2 += bkg;
1318       
1319        return(f2);
1320}
1321
1322double
1323TwoShell(double dp[], double q)
1324{
1325        // variables are:
1326        //[0] scale factor
1327        //[1] radius of core [ï¿œ]
1328        //[2] SLD of the core   [ï¿œ-2]
1329        //[3] thickness of shell 1 [ï¿œ]
1330        //[4] SLD of shell 1
1331        //[5] thickness of shell 2 [ï¿œ]
1332        //[6] SLD of shell 2
1333        //[7] SLD of the solvent
1334        //[8] background        [cm-1]
1335       
1336        double x,pi;
1337        double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;         //my local names
1338        double bes,f,vol,qr,contr,f2;
1339        double rhoshel2,thick2;
1340       
1341        pi = 4.0*atan(1.0);
1342        x=q;
1343       
1344        scale = dp[0];
1345        rcore = dp[1];
1346        rhocore = dp[2];
1347        thick1 = dp[3];
1348        rhoshel1 = dp[4];
1349        thick2 = dp[5];
1350        rhoshel2 = dp[6];       
1351        rhosolv = dp[7];
1352        bkg = dp[8];
1353                // core first, then add in shells
1354        qr=x*rcore;
1355        contr = rhocore-rhoshel1;
1356        if(qr == 0){
1357                bes = 1.0;
1358        }else{
1359                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1360        }
1361        vol = 4.0*pi/3.0*rcore*rcore*rcore;
1362        f = vol*bes*contr;
1363        //now the shell (1)
1364        qr=x*(rcore+thick1);
1365        contr = rhoshel1-rhoshel2;
1366        if(qr == 0){
1367                bes = 1.0;
1368        }else{
1369                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1370        }
1371        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1);
1372        f += vol*bes*contr;
1373        //now the shell (2)
1374        qr=x*(rcore+thick1+thick2);
1375        contr = rhoshel2-rhosolv;
1376        if(qr == 0){
1377                bes = 1.0;
1378        }else{
1379                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1380        }
1381        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2);
1382        f += vol*bes*contr;
1383       
1384               
1385        // normalize to particle volume and rescale from [ï¿œ-1] to [cm-1]
1386        f2 = f*f/vol*1.0e8;
1387       
1388        //scale if desired
1389        f2 *= scale;
1390        // then add in the background
1391        f2 += bkg;
1392       
1393        return(f2);
1394}
1395
1396double
1397ThreeShell(double dp[], double q)
1398{
1399        // variables are:
1400        //[0] scale factor
1401        //[1] radius of core [ï¿œ]
1402        //[2] SLD of the core   [ï¿œ-2]
1403        //[3] thickness of shell 1 [ï¿œ]
1404        //[4] SLD of shell 1
1405        //[5] thickness of shell 2 [ï¿œ]
1406        //[6] SLD of shell 2
1407        //[7] thickness of shell 3
1408        //[8] SLD of shell 3
1409        //[9] SLD of solvent
1410        //[10] background       [cm-1]
1411       
1412        double x,pi;
1413        double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;         //my local names
1414        double bes,f,vol,qr,contr,f2;
1415        double rhoshel2,thick2,rhoshel3,thick3;
1416       
1417        pi = 4.0*atan(1.0);
1418        x=q;
1419       
1420        scale = dp[0];
1421        rcore = dp[1];
1422        rhocore = dp[2];
1423        thick1 = dp[3];
1424        rhoshel1 = dp[4];
1425        thick2 = dp[5];
1426        rhoshel2 = dp[6];       
1427        thick3 = dp[7];
1428        rhoshel3 = dp[8];       
1429        rhosolv = dp[9];
1430        bkg = dp[10];
1431       
1432                // core first, then add in shells
1433        qr=x*rcore;
1434        contr = rhocore-rhoshel1;
1435        if(qr == 0){
1436                bes = 1.0;
1437        }else{
1438                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1439        }
1440        vol = 4.0*pi/3.0*rcore*rcore*rcore;
1441        f = vol*bes*contr;
1442        //now the shell (1)
1443        qr=x*(rcore+thick1);
1444        contr = rhoshel1-rhoshel2;
1445        if(qr == 0){
1446                bes = 1.0;
1447        }else{
1448                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1449        }
1450        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1);
1451        f += vol*bes*contr;
1452        //now the shell (2)
1453        qr=x*(rcore+thick1+thick2);
1454        contr = rhoshel2-rhoshel3;
1455        if(qr == 0){
1456                bes = 1.0;
1457        }else{
1458                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1459        }
1460        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2);
1461        f += vol*bes*contr;
1462        //now the shell (3)
1463        qr=x*(rcore+thick1+thick2+thick3);
1464        contr = rhoshel3-rhosolv;
1465        if(qr == 0){
1466                bes = 1.0;
1467        }else{
1468                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1469        }
1470        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3);
1471        f += vol*bes*contr;
1472               
1473        // normalize to particle volume and rescale from [ï¿œ-1] to [cm-1]
1474        f2 = f*f/vol*1.0e8;
1475       
1476        //scale if desired
1477        f2 *= scale;
1478        // then add in the background
1479        f2 += bkg;
1480       
1481        return(f2);
1482}
1483
1484double
1485FourShell(double dp[], double q)
1486{
1487        // variables are:
1488        //[0] scale factor
1489        //[1] radius of core [ï¿œ]
1490        //[2] SLD of the core   [ï¿œ-2]
1491        //[3] thickness of shell 1 [ï¿œ]
1492        //[4] SLD of shell 1
1493        //[5] thickness of shell 2 [ï¿œ]
1494        //[6] SLD of shell 2
1495        //[7] thickness of shell 3
1496        //[8] SLD of shell 3
1497        //[9] thickness of shell 3
1498        //[10] SLD of shell 3
1499        //[11] SLD of solvent
1500        //[12] background       [cm-1]
1501       
1502        double x,pi;
1503        double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;         //my local names
1504        double bes,f,vol,qr,contr,f2;
1505        double rhoshel2,thick2,rhoshel3,thick3,rhoshel4,thick4;
1506       
1507        pi = 4.0*atan(1.0);
1508        x=q;
1509       
1510        scale = dp[0];
1511        rcore = dp[1];
1512        rhocore = dp[2];
1513        thick1 = dp[3];
1514        rhoshel1 = dp[4];
1515        thick2 = dp[5];
1516        rhoshel2 = dp[6];       
1517        thick3 = dp[7];
1518        rhoshel3 = dp[8];
1519        thick4 = dp[9];
1520        rhoshel4 = dp[10];     
1521        rhosolv = dp[11];
1522        bkg = dp[12];
1523       
1524                // core first, then add in shells
1525        qr=x*rcore;
1526        contr = rhocore-rhoshel1;
1527        if(qr == 0){
1528                bes = 1.0;
1529        }else{
1530                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1531        }
1532        vol = 4.0*pi/3.0*rcore*rcore*rcore;
1533        f = vol*bes*contr;
1534        //now the shell (1)
1535        qr=x*(rcore+thick1);
1536        contr = rhoshel1-rhoshel2;
1537        if(qr == 0){
1538                bes = 1.0;
1539        }else{
1540                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1541        }
1542        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1);
1543        f += vol*bes*contr;
1544        //now the shell (2)
1545        qr=x*(rcore+thick1+thick2);
1546        contr = rhoshel2-rhoshel3;
1547        if(qr == 0){
1548                bes = 1.0;
1549        }else{
1550                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1551        }
1552        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2);
1553        f += vol*bes*contr;
1554        //now the shell (3)
1555        qr=x*(rcore+thick1+thick2+thick3);
1556        contr = rhoshel3-rhoshel4;
1557        if(qr == 0){
1558                bes = 1.0;
1559        }else{
1560                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1561        }
1562        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3);
1563        f += vol*bes*contr;
1564        //now the shell (4)
1565        qr=x*(rcore+thick1+thick2+thick3+thick4);
1566        contr = rhoshel4-rhosolv;
1567        if(qr == 0){
1568                bes = 1.0;
1569        }else{
1570                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
1571        }
1572        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4);
1573        f += vol*bes*contr;
1574       
1575               
1576        // normalize to particle volume and rescale from [ï¿œ-1] to [cm-1]
1577        f2 = f*f/vol*1.0e8;
1578       
1579        //scale if desired
1580        f2 *= scale;
1581        // then add in the background
1582        f2 += bkg;
1583       
1584        return(f2);
1585}
1586
1587double
1588PolyOneShell(double dp[], double x)
1589{
1590        double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg,pd,zz;             //my local names
1591        double va,vb,summ,yyy,zi;
1592        double answer,zp1,zp2,zp3,vpoly,range,temp_1sf[7],pi;
1593        int nord=76,ii;
1594       
1595        pi = 4.0*atan(1.0);
1596       
1597        scale = dp[0];
1598        rcore = dp[1];
1599        pd = dp[2];
1600        rhocore = dp[3];
1601        thick = dp[4];
1602        rhoshel = dp[5];
1603        rhosolv = dp[6];
1604        bkg = dp[7];
1605               
1606        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only
1607       
1608        range = 8.0;                    //std deviations for the integration
1609        va = rcore*(1.0-range*pd);
1610        if (va<0.0) {
1611                va=0.0;         //otherwise numerical error when pd >= 0.3, making a<0
1612        }
1613        if (pd>0.3) {
1614                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail
1615        }
1616        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius?
1617       
1618        //temp set scale=1 and bkg=0 for quadrature calc
1619        temp_1sf[0] = 1.0;
1620        temp_1sf[1] = dp[1];    //the core radius will be changed in the loop
1621        temp_1sf[2] = dp[3];
1622        temp_1sf[3] = dp[4];
1623        temp_1sf[4] = dp[5];
1624        temp_1sf[5] = dp[6];
1625        temp_1sf[6] = 0.0;
1626       
1627        summ = 0.0;             // initialize integral
1628        for(ii=0;ii<nord;ii+=1) {
1629                // calculate Gauss points on integration interval (r-value for evaluation)
1630                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
1631                temp_1sf[1] = zi;
1632                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * OneShell(temp_1sf,x);
1633                //un-normalize by volume
1634                yyy *= 4.0*pi/3.0*pow((zi+thick),3);
1635                summ += yyy;            //add to the running total of the quadrature
1636        }
1637        // calculate value of integral to return
1638        answer = (vb-va)/2.0*summ;
1639       
1640        //re-normalize by the average volume
1641        zp1 = zz + 1.0;
1642        zp2 = zz + 2.0;
1643        zp3 = zz + 3.0;
1644        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick),3);
1645        answer /= vpoly;
1646//scale
1647        answer *= scale;
1648// add in the background
1649        answer += bkg;
1650               
1651    return(answer);
1652}
1653
1654double
1655PolyTwoShell(double dp[], double x)
1656{
1657        double scale,rcore,rhocore,rhosolv,bkg,pd,zz;           //my local names
1658        double va,vb,summ,yyy,zi;
1659        double answer,zp1,zp2,zp3,vpoly,range,temp_2sf[9],pi;
1660        int nord=76,ii;
1661        double thick1,thick2;
1662        double rhoshel1,rhoshel2;
1663       
1664        scale = dp[0];
1665        rcore = dp[1];
1666        pd = dp[2];
1667        rhocore = dp[3];
1668        thick1 = dp[4];
1669        rhoshel1 = dp[5];
1670        thick2 = dp[6];
1671        rhoshel2 = dp[7];
1672        rhosolv = dp[8];
1673        bkg = dp[9];
1674       
1675        pi = 4.0*atan(1.0);
1676               
1677        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only
1678       
1679        range = 8.0;                    //std deviations for the integration
1680        va = rcore*(1.0-range*pd);
1681        if (va<0.0) {
1682                va=0.0;         //otherwise numerical error when pd >= 0.3, making a<0
1683        }
1684        if (pd>0.3) {
1685                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail
1686        }
1687        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius?
1688       
1689        //temp set scale=1 and bkg=0 for quadrature calc
1690        temp_2sf[0] = 1.0;
1691        temp_2sf[1] = dp[1];            //the core radius will be changed in the loop
1692        temp_2sf[2] = dp[3];
1693        temp_2sf[3] = dp[4];
1694        temp_2sf[4] = dp[5];
1695        temp_2sf[5] = dp[6];
1696        temp_2sf[6] = dp[7];
1697        temp_2sf[7] = dp[8];
1698        temp_2sf[8] = 0.0;
1699       
1700        summ = 0.0;             // initialize integral
1701        for(ii=0;ii<nord;ii+=1) {
1702                // calculate Gauss points on integration interval (r-value for evaluation)
1703                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
1704                temp_2sf[1] = zi;
1705                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * TwoShell(temp_2sf,x);
1706                //un-normalize by volume
1707                yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2),3);
1708                summ += yyy;            //add to the running total of the quadrature
1709        }
1710        // calculate value of integral to return
1711        answer = (vb-va)/2.0*summ;
1712       
1713        //re-normalize by the average volume
1714        zp1 = zz + 1.0;
1715        zp2 = zz + 2.0;
1716        zp3 = zz + 3.0;
1717        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2),3);
1718        answer /= vpoly;
1719//scale
1720        answer *= scale;
1721// add in the background
1722        answer += bkg;
1723               
1724    return(answer);
1725}
1726
1727double
1728PolyThreeShell(double dp[], double x)
1729{
1730        double scale,rcore,rhocore,rhosolv,bkg,pd,zz;           //my local names
1731        double va,vb,summ,yyy,zi;
1732        double answer,zp1,zp2,zp3,vpoly,range,temp_3sf[11],pi;
1733        int nord=76,ii;
1734        double thick1,thick2,thick3;
1735        double rhoshel1,rhoshel2,rhoshel3;
1736       
1737        scale = dp[0];
1738        rcore = dp[1];
1739        pd = dp[2];
1740        rhocore = dp[3];
1741        thick1 = dp[4];
1742        rhoshel1 = dp[5];
1743        thick2 = dp[6];
1744        rhoshel2 = dp[7];
1745        thick3 = dp[8];
1746        rhoshel3 = dp[9];
1747        rhosolv = dp[10];
1748        bkg = dp[11];
1749       
1750        pi = 4.0*atan(1.0);
1751               
1752        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only
1753       
1754        range = 8.0;                    //std deviations for the integration
1755        va = rcore*(1.0-range*pd);
1756        if (va<0) {
1757                va=0;           //otherwise numerical error when pd >= 0.3, making a<0
1758        }
1759        if (pd>0.3) {
1760                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail
1761        }
1762        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius?
1763       
1764        //temp set scale=1 and bkg=0 for quadrature calc
1765        temp_3sf[0] = 1.0;
1766        temp_3sf[1] = dp[1];            //the core radius will be changed in the loop
1767        temp_3sf[2] = dp[3];
1768        temp_3sf[3] = dp[4];
1769        temp_3sf[4] = dp[5];
1770        temp_3sf[5] = dp[6];
1771        temp_3sf[6] = dp[7];
1772        temp_3sf[7] = dp[8];
1773        temp_3sf[8] = dp[9];
1774        temp_3sf[9] = dp[10];
1775        temp_3sf[10] = 0.0;
1776       
1777        summ = 0.0;             // initialize integral
1778        for(ii=0;ii<nord;ii+=1) {
1779                // calculate Gauss points on integration interval (r-value for evaluation)
1780                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
1781                temp_3sf[1] = zi;
1782                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * ThreeShell(temp_3sf,x);
1783                //un-normalize by volume
1784                yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2+thick3),3);
1785                summ += yyy;            //add to the running total of the quadrature
1786        }
1787        // calculate value of integral to return
1788        answer = (vb-va)/2.0*summ;
1789       
1790        //re-normalize by the average volume
1791        zp1 = zz + 1.0;
1792        zp2 = zz + 2.0;
1793        zp3 = zz + 3.0;
1794        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2+thick3),3);
1795        answer /= vpoly;
1796//scale
1797        answer *= scale;
1798// add in the background
1799        answer += bkg;
1800               
1801    return(answer);
1802}
1803
1804double
1805PolyFourShell(double dp[], double x)
1806{
1807        double scale,rcore,rhocore,rhosolv,bkg,pd,zz;           //my local names
1808        double va,vb,summ,yyy,zi;
1809        double answer,zp1,zp2,zp3,vpoly,range,temp_4sf[13],pi;
1810        int nord=76,ii;
1811        double thick1,thick2,thick3,thick4;
1812        double rhoshel1,rhoshel2,rhoshel3,rhoshel4;
1813       
1814        scale = dp[0];
1815        rcore = dp[1];
1816        pd = dp[2];
1817        rhocore = dp[3];
1818        thick1 = dp[4];
1819        rhoshel1 = dp[5];
1820        thick2 = dp[6];
1821        rhoshel2 = dp[7];
1822        thick3 = dp[8];
1823        rhoshel3 = dp[9];
1824        thick4 = dp[10];
1825        rhoshel4 = dp[11];
1826        rhosolv = dp[12];
1827        bkg = dp[13];
1828       
1829        pi = 4.0*atan(1.0);
1830               
1831        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only
1832       
1833        range = 8.0;                    //std deviations for the integration
1834        va = rcore*(1.0-range*pd);
1835        if (va<0) {
1836                va=0;           //otherwise numerical error when pd >= 0.3, making a<0
1837        }
1838        if (pd>0.3) {
1839                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail
1840        }
1841        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius?
1842       
1843        //temp set scale=1 and bkg=0 for quadrature calc
1844        temp_4sf[0] = 1.0;
1845        temp_4sf[1] = dp[1];            //the core radius will be changed in the loop
1846        temp_4sf[2] = dp[3];
1847        temp_4sf[3] = dp[4];
1848        temp_4sf[4] = dp[5];
1849        temp_4sf[5] = dp[6];
1850        temp_4sf[6] = dp[7];
1851        temp_4sf[7] = dp[8];
1852        temp_4sf[8] = dp[9];
1853        temp_4sf[9] = dp[10];
1854        temp_4sf[10] = dp[11];
1855        temp_4sf[11] = dp[12];
1856        temp_4sf[12] = 0.0;
1857               
1858        summ = 0.0;             // initialize integral
1859        for(ii=0;ii<nord;ii+=1) {
1860                // calculate Gauss points on integration interval (r-value for evaluation)
1861                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
1862                temp_4sf[1] = zi;
1863                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * FourShell(temp_4sf,x);
1864                //un-normalize by volume
1865                yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2+thick3+thick4),3);
1866                summ += yyy;            //add to the running total of the quadrature
1867        }
1868        // calculate value of integral to return
1869        answer = (vb-va)/2.0*summ;
1870       
1871        //re-normalize by the average volume
1872        zp1 = zz + 1.0;
1873        zp2 = zz + 2.0;
1874        zp3 = zz + 3.0;
1875        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2+thick3+thick4),3);
1876        answer /= vpoly;
1877//scale
1878        answer *= scale;
1879// add in the background
1880        answer += bkg;
1881               
1882    return(answer);
1883}
1884
1885
1886/*      BCC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x
1887
1888Uses 150 pt Gaussian quadrature for both integrals
1889
1890*/
1891double
1892BCC_ParaCrystal(double w[], double x)
1893{
1894        int i,j;
1895        double Pi;
1896        double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;           //local variables of coefficient wave
1897        int nordi=150;                  //order of integration
1898        int nordj=150;
1899        double va,vb;           //upper and lower integration limits
1900        double summ,zi,yyy,answer;                      //running tally of integration
1901        double summj,vaj,vbj,zij;                       //for the inner integration
1902       
1903        Pi = 4.0*atan(1.0);
1904        va = 0.0;
1905        vb = 2.0*Pi;            //orintational average, outer integral
1906        vaj = 0.0;
1907        vbj = Pi;               //endpoints of inner integral
1908       
1909        summ = 0.0;                     //initialize intergral
1910       
1911        scale = w[0];
1912        Dnn = w[1];                                     //Nearest neighbor distance A
1913        gg = w[2];                                      //Paracrystal distortion factor
1914        Rad = w[3];                                     //Sphere radius
1915        sld = w[4];
1916        sldSolv = w[5];
1917        background = w[6]; 
1918       
1919        contrast = sld - sldSolv;
1920       
1921        //Volume fraction calculated from lattice symmetry and sphere radius
1922        latticeScale = 2.0*(4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn/sqrt(3.0/4.0),3);
1923       
1924        for(i=0;i<nordi;i++) {
1925                //setup inner integral over the ellipsoidal cross-section
1926                summj=0.0;
1927                zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;            //the outer dummy is phi
1928                for(j=0;j<nordj;j++) {
1929                        //20 gauss points for the inner integral
1930                        zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;               //the inner dummy is theta
1931                        yyy = Gauss150Wt[j] * BCC_Integrand(w,x,zi,zij);
1932                        summj += yyy;
1933                }
1934                //now calculate the value of the inner integral
1935                answer = (vbj-vaj)/2.0*summj;
1936               
1937                //now calculate outer integral
1938                yyy = Gauss150Wt[i] * answer;
1939                summ += yyy;
1940        }               //final scaling is done at the end of the function, after the NT_FP64 case
1941       
1942        answer = (vb-va)/2.0*summ;
1943        // Multiply by contrast^2
1944        answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale;
1945        // add in the background
1946        answer += background;
1947       
1948        return answer;
1949}
1950
1951// xx is phi (outer)
1952// yy is theta (inner)
1953double
1954BCC_Integrand(double w[], double qq, double xx, double yy) {
1955       
1956        double retVal,temp1,temp3,aa,Da,Dnn,gg,Pi;
1957       
1958        Dnn = w[1]; //Nearest neighbor distance A
1959        gg = w[2]; //Paracrystal distortion factor
1960        aa = Dnn;
1961        Da = gg*aa;
1962       
1963        Pi = 4.0*atan(1.0);
1964        temp1 = qq*qq*Da*Da;
1965        temp3 = qq*aa; 
1966       
1967        retVal = BCCeval(yy,xx,temp1,temp3);
1968        retVal /=4.0*Pi;
1969       
1970        return(retVal);
1971}
1972
1973double
1974BCCeval(double Theta, double Phi, double temp1, double temp3) {
1975
1976        double temp6,temp7,temp8,temp9,temp10;
1977        double result;
1978       
1979        temp6 = sin(Theta);
1980        temp7 = sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)+cos(Theta);
1981        temp8 = -1.0*sin(Theta)*cos(Phi)-sin(Theta)*sin(Phi)+cos(Theta);
1982        temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)-cos(Theta);
1983        temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
1984        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)));
1985       
1986        return (result);
1987}
1988
1989double
1990SphereForm_Paracrystal(double radius, double delrho, double x) {                                       
1991       
1992        double bes,f,vol,f2,pi;
1993        pi = 4.0*atan(1.0);
1994        //
1995        //handle q==0 separately
1996        if(x==0) {
1997                f = 4.0/3.0*pi*radius*radius*radius*delrho*delrho*1.0e8;
1998                return(f);
1999        }
2000       
2001        bes = 3.0*(sin(x*radius)-x*radius*cos(x*radius))/(x*x*x)/(radius*radius*radius);
2002        vol = 4.0*pi/3.0*radius*radius*radius;
2003        f = vol*bes*delrho      ;       // [=] ï¿œ
2004        // normalize to single particle volume, convert to 1/cm
2005        f2 = f * f / vol * 1.0e8;               // [=] 1/cm
2006       
2007        return (f2);
2008}
2009
2010/*      FCC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x
2011
2012Uses 150 pt Gaussian quadrature for both integrals
2013
2014*/
2015double
2016FCC_ParaCrystal(double w[], double x)
2017{
2018        int i,j;
2019        double Pi;
2020        double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;           //local variables of coefficient wave
2021        int nordi=150;                  //order of integration
2022        int nordj=150;
2023        double va,vb;           //upper and lower integration limits
2024        double summ,zi,yyy,answer;                      //running tally of integration
2025        double summj,vaj,vbj,zij;                       //for the inner integration
2026       
2027        Pi = 4.0*atan(1.0);
2028        va = 0.0;
2029        vb = 2.0*Pi;            //orintational average, outer integral
2030        vaj = 0.0;
2031        vbj = Pi;               //endpoints of inner integral
2032       
2033        summ = 0.0;                     //initialize intergral
2034       
2035        scale = w[0];
2036        Dnn = w[1];                                     //Nearest neighbor distance A
2037        gg = w[2];                                      //Paracrystal distortion factor
2038        Rad = w[3];                                     //Sphere radius
2039        sld = w[4];
2040        sldSolv = w[5];
2041        background = w[6]; 
2042       
2043        contrast = sld - sldSolv;
2044        //Volume fraction calculated from lattice symmetry and sphere radius
2045        latticeScale = 4.0*(4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn*sqrt(2.0),3);
2046       
2047        for(i=0;i<nordi;i++) {
2048                //setup inner integral over the ellipsoidal cross-section
2049                summj=0.0;
2050                zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;            //the outer dummy is phi
2051                for(j=0;j<nordj;j++) {
2052                        //20 gauss points for the inner integral
2053                        zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;               //the inner dummy is theta
2054                        yyy = Gauss150Wt[j] * FCC_Integrand(w,x,zi,zij);
2055                        summj += yyy;
2056                }
2057                //now calculate the value of the inner integral
2058                answer = (vbj-vaj)/2.0*summj;
2059               
2060                //now calculate outer integral
2061                yyy = Gauss150Wt[i] * answer;
2062                summ += yyy;
2063        }               //final scaling is done at the end of the function, after the NT_FP64 case
2064       
2065        answer = (vb-va)/2.0*summ;
2066        // Multiply by contrast^2
2067        answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale;
2068        // add in the background
2069        answer += background;
2070       
2071        return answer;
2072}
2073
2074
2075// xx is phi (outer)
2076// yy is theta (inner)
2077double
2078FCC_Integrand(double w[], double qq, double xx, double yy) {
2079       
2080        double retVal,temp1,temp3,aa,Da,Dnn,gg,Pi;
2081       
2082        Pi = 4.0*atan(1.0);
2083        Dnn = w[1]; //Nearest neighbor distance A
2084        gg = w[2]; //Paracrystal distortion factor
2085        aa = Dnn;
2086        Da = gg*aa;
2087       
2088        temp1 = qq*qq*Da*Da;
2089        temp3 = qq*aa;
2090       
2091        retVal = FCCeval(yy,xx,temp1,temp3);
2092        retVal /=4*Pi;
2093       
2094        return(retVal);
2095}
2096
2097double
2098FCCeval(double Theta, double Phi, double temp1, double temp3) {
2099
2100        double temp6,temp7,temp8,temp9,temp10;
2101        double result;
2102       
2103        temp6 = sin(Theta);
2104        temp7 = sin(Theta)*sin(Phi)+cos(Theta);
2105        temp8 = -1.0*sin(Theta)*cos(Phi)+cos(Theta);
2106        temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi);
2107        temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
2108        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)));
2109       
2110        return (result);
2111}
2112
2113
2114/*      SC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x
2115
2116Uses 150 pt Gaussian quadrature for both integrals
2117
2118*/
2119double
2120SC_ParaCrystal(double w[], double x)
2121{
2122        int i,j;
2123        double Pi;
2124        double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;           //local variables of coefficient wave
2125        int nordi=150;                  //order of integration
2126        int nordj=150;
2127        double va,vb;           //upper and lower integration limits
2128        double summ,zi,yyy,answer;                      //running tally of integration
2129        double summj,vaj,vbj,zij;                       //for the inner integration
2130       
2131        Pi = 4.0*atan(1.0);
2132        va = 0.0;
2133        vb = Pi/2.0;            //orintational average, outer integral
2134        vaj = 0.0;
2135        vbj = Pi/2.0;           //endpoints of inner integral
2136       
2137        summ = 0.0;                     //initialize intergral
2138       
2139        scale = w[0];
2140        Dnn = w[1];                                     //Nearest neighbor distance A
2141        gg = w[2];                                      //Paracrystal distortion factor
2142        Rad = w[3];                                     //Sphere radius
2143        sld = w[4];
2144        sldSolv = w[5];
2145        background = w[6]; 
2146       
2147        contrast = sld - sldSolv;
2148        //Volume fraction calculated from lattice symmetry and sphere radius
2149        latticeScale = (4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn,3);
2150       
2151        for(i=0;i<nordi;i++) {
2152                //setup inner integral over the ellipsoidal cross-section
2153                summj=0.0;
2154                zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;            //the outer dummy is phi
2155                for(j=0;j<nordj;j++) {
2156                        //20 gauss points for the inner integral
2157                        zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;               //the inner dummy is theta
2158                        yyy = Gauss150Wt[j] * SC_Integrand(w,x,zi,zij);
2159                        summj += yyy;
2160                }
2161                //now calculate the value of the inner integral
2162                answer = (vbj-vaj)/2.0*summj;
2163               
2164                //now calculate outer integral
2165                yyy = Gauss150Wt[i] * answer;
2166                summ += yyy;
2167        }               //final scaling is done at the end of the function, after the NT_FP64 case
2168       
2169        answer = (vb-va)/2.0*summ;
2170        // Multiply by contrast^2
2171        answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale;
2172        // add in the background
2173        answer += background;
2174       
2175        return answer;
2176}
2177
2178// xx is phi (outer)
2179// yy is theta (inner)
2180double
2181SC_Integrand(double w[], double qq, double xx, double yy) {
2182       
2183        double retVal,temp1,temp2,temp3,temp4,temp5,aa,Da,Dnn,gg,Pi;
2184       
2185        Pi = 4.0*atan(1.0);
2186        Dnn = w[1]; //Nearest neighbor distance A
2187        gg = w[2]; //Paracrystal distortion factor
2188        aa = Dnn;
2189        Da = gg*aa;
2190       
2191        temp1 = qq*qq*Da*Da;
2192        temp2 = pow( 1.0-exp(-1.0*temp1) ,3);
2193        temp3 = qq*aa;
2194        temp4 = 2.0*exp(-0.5*temp1);
2195        temp5 = exp(-1.0*temp1);
2196       
2197       
2198        retVal = temp2*SCeval(yy,xx,temp3,temp4,temp5);
2199        retVal *= 2.0/Pi;
2200       
2201        return(retVal);
2202}
2203
2204double
2205SCeval(double Theta, double Phi, double temp3, double temp4, double temp5) { //Function to calculate integrand values for simple cubic structure
2206
2207        double temp6,temp7,temp8,temp9; //Theta and phi dependent parts of the equation
2208        double result;
2209       
2210        temp6 = sin(Theta);
2211        temp7 = -1.0*temp3*sin(Theta)*cos(Phi);
2212        temp8 = temp3*sin(Theta)*sin(Phi);
2213        temp9 = temp3*cos(Theta);
2214        result = temp6/((1.0-temp4*cos((temp7))+temp5)*(1.0-temp4*cos((temp8))+temp5)*(1.0-temp4*cos((temp9))+temp5));
2215       
2216        return (result);
2217}
2218
2219// scattering from a uniform sphere with a Gaussian size distribution
2220//
2221double
2222FuzzySpheres(double dp[], double q)
2223{
2224        double pi,x;
2225        double scale,rad,pd,sig,rho,rhos,bkg,delrho,sig_surf,f2,bes,vol,f;              //my local names
2226        double va,vb,zi,yy,summ,inten;
2227        int nord=20,ii;
2228       
2229        pi = 4.0*atan(1.0);
2230        x= q;
2231       
2232        scale=dp[0];
2233        rad=dp[1];
2234        pd=dp[2];
2235        sig=pd*rad;
2236        sig_surf = dp[3];
2237        rho=dp[4];
2238        rhos=dp[5];
2239        delrho=rho-rhos;
2240        bkg=dp[6];
2241       
2242                       
2243        va = -4.0*sig + rad;
2244        if (va<0) {
2245                va=0;           //to avoid numerical error when  va<0 (-ve q-value)
2246        }
2247        vb = 4.0*sig +rad;
2248       
2249        summ = 0.0;             // initialize integral
2250        for(ii=0;ii<nord;ii+=1) {
2251                // calculate Gauss points on integration interval (r-value for evaluation)
2252                zi = ( Gauss20Z[ii]*(vb-va) + vb + va )/2.0;
2253                // calculate sphere scattering
2254                //
2255                //handle q==0 separately
2256                if (x==0.0) {
2257                        f2 = 4.0/3.0*pi*zi*zi*zi*delrho*delrho*1.0e8;
2258                        f2 *= exp(-0.5*sig_surf*sig_surf*x*x);
2259                        f2 *= exp(-0.5*sig_surf*sig_surf*x*x);
2260                } else {
2261                        bes = 3.0*(sin(x*zi)-x*zi*cos(x*zi))/(x*x*x)/(zi*zi*zi);
2262                        vol = 4.0*pi/3.0*zi*zi*zi;
2263                        f = vol*bes*delrho;             // [=] A
2264                        f *= exp(-0.5*sig_surf*sig_surf*x*x);
2265                        // normalize to single particle volume, convert to 1/cm
2266                        f2 = f * f / vol * 1.0e8;               // [=] 1/cm
2267                }
2268       
2269                yy = Gauss20Wt[ii] *  Gauss_distr(sig,rad,zi) * f2;
2270                yy *= 4.0*pi/3.0*zi*zi*zi;              //un-normalize by current sphere volume
2271               
2272                summ += yy;             //add to the running total of the quadrature
2273               
2274               
2275        }
2276        // calculate value of integral to return
2277        inten = (vb-va)/2.0*summ;
2278       
2279        //re-normalize by polydisperse sphere volume
2280        inten /= (4.0*pi/3.0*rad*rad*rad)*(1.0+3.0*pd*pd);
2281       
2282        inten *= scale;
2283        inten += bkg;
2284       
2285    return(inten);      //scale, and add in the background
2286}
2287
2288
Note: See TracBrowser for help on using the repository browser.