source: sasview/sansmodels/src/libigor/libSphere.c @ 31af069

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 31af069 was 67424cd, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

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