source: sans/Analysis/branches/ajj_23APR07/XOPs/SANSAnalysis/Sphere.c @ 94

Last change on this file since 94 was 93, checked in by ajj, 16 years ago

Add combined XOP code. Currently only contains XCode project file to build Universal binary suitable for Igor 6.
File size: 54.1 KB
Line 
1/*      SimpleFit.c
2
3        A simplified project designed to act as a template for your curve fitting function.
4        The fitting function is a simple polynomial. It works but is of no practical use.
5*/
6
7#pragma XOP_SET_STRUCT_PACKING                  // All structures are 2-byte-aligned.
8
9#include "XOPStandardHeaders.h"                 // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h
10#include "SANSAnalysis.h"
11#include "Sphere.h"
12
13// scattering from a sphere - hardly needs to be an XOP...
14int
15SphereFormX(FitParamsPtr p)
16{
17        DOUBLE *dp;                             // Pointer to double precision wave data.
18        float *fp;                              // Pointer to single precision wave data.
19        DOUBLE q;
20        DOUBLE scale,radius,delrho,bkg;         //my local names
21        DOUBLE bes,f,vol,f2,pi;
22       
23        if (p->waveHandle == NIL) {
24                SetNaN64(&p->result);
25                return NON_EXISTENT_WAVE;
26        }
27       
28        pi = 4.0*atan(1.0);
29        q= p->x;
30       
31        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
32                case NT_FP32:
33                        fp= WaveData(p->waveHandle);
34                        scale = fp[0];
35                        radius = fp[1];
36                        delrho = fp[2];
37                        bkg = fp[3];
38                   
39                        break;
40                case NT_FP64:
41                        dp= WaveData(p->waveHandle);
42                        scale = dp[0];
43                        radius = dp[1];
44                        delrho = dp[2];
45                        bkg = dp[3];
46                                               
47                        break;
48                default:                                                                // We can't handle this wave data type.
49                        SetNaN64(&p->result);
50                        return REQUIRES_SP_OR_DP_WAVE;
51        }
52        //handle q==0 separately
53        if(q==0){
54            f = 4.0/3.0*pi*radius*radius*radius*delrho*delrho*scale*1.0e8 + bkg;
55            p->result= f;
56            return(0);
57        }
58       
59        bes = 3.0*(sin(q*radius)-q*radius*cos(q*radius))/(q*q*q)/(radius*radius*radius);
60        vol = 4.0*pi/3.0*radius*radius*radius;
61        f = vol*bes*delrho;             // [=] 
62        // normalize to single particle volume, convert to 1/cm
63        f2 = f * f / vol * 1.0e8;               // [=] 1/cm
64       
65        p->result= (scale*f2+bkg);      //scale, and add in the background
66        return 0;
67}
68
69// scattering from a monodisperse core-shell sphere - hardly needs to be an XOP...
70int
71CoreShellFormX(FitParamsPtr p)
72{
73        DOUBLE *dp;                             // Pointer to double precision wave data.
74        float *fp;                              // Pointer to single precision wave data.
75        DOUBLE x,pi;
76        DOUBLE scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg;           //my local names
77        DOUBLE bes,f,vol,qr,contr,f2;
78       
79        if (p->waveHandle == NIL) {
80                SetNaN64(&p->result);
81                return NON_EXISTENT_WAVE;
82        }
83       
84        pi = 4.0*atan(1.0);
85        x= p->x;
86       
87        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
88                case NT_FP32:
89                        fp= WaveData(p->waveHandle);
90                        scale = fp[0];
91                        rcore = fp[1];
92                        thick = fp[2];
93                        rhocore = fp[3];
94                        rhoshel = fp[4];
95                        rhosolv = fp[5];
96                        bkg = fp[6];
97                   
98                        break;
99                case NT_FP64:
100                        dp= WaveData(p->waveHandle);
101                        scale = dp[0];
102                        rcore = dp[1];
103                        thick = dp[2];
104                        rhocore = dp[3];
105                        rhoshel = dp[4];
106                        rhosolv = dp[5];
107                        bkg = dp[6];
108                                               
109                        break;
110                default:                                                                // We can't handle this wave data type.
111                        SetNaN64(&p->result);
112                        return REQUIRES_SP_OR_DP_WAVE;
113        }
114       // core first, then add in shell
115        qr=x*rcore;
116        contr = rhocore-rhoshel;
117        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
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        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
124        vol = 4.0*pi/3.0*pow((rcore+thick),3);
125        f += vol*bes*contr;
126       
127        // normalize to particle volume and rescale from [-1] to [cm-1]
128        f2 = f*f/vol*1.0e8;
129       
130        //scale if desired
131        f2 *= scale;
132        // then add in the background
133        f2 += bkg;
134               
135        p->result= f2;
136        return 0;
137}
138
139// scattering from a unilamellar vesicle - hardly needs to be an XOP...
140// same functional form as the core-shell sphere, but more intuitive for a vesicle
141int
142VesicleFormX(FitParamsPtr p)
143{
144        DOUBLE *dp;                             // Pointer to double precision wave data.
145        float *fp;                              // Pointer to single precision wave data.
146        DOUBLE x,pi;
147        DOUBLE scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg;           //my local names
148        DOUBLE bes,f,vol,qr,contr,f2;
149       
150        if (p->waveHandle == NIL) {
151                SetNaN64(&p->result);
152                return NON_EXISTENT_WAVE;
153        }
154       
155        pi = 4.0*atan(1.0);
156        x= p->x;
157       
158        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
159                case NT_FP32:
160                        fp= WaveData(p->waveHandle);
161                        scale = fp[0];
162                        rcore = fp[1];
163                        thick = fp[2];
164                        rhocore = fp[3];
165                        rhosolv = rhocore;
166                        rhoshel = fp[4];
167                        bkg = fp[5];
168                   
169                        break;
170                case NT_FP64:
171                        dp= WaveData(p->waveHandle);
172                        scale = dp[0];
173                        rcore = dp[1];
174                        thick = dp[2];
175                        rhocore = dp[3];
176                        rhosolv = rhocore;
177                        rhoshel = dp[4];
178                        bkg = dp[5];
179                                               
180                        break;
181                default:                                                                // We can't handle this wave data type.
182                        SetNaN64(&p->result);
183                        return REQUIRES_SP_OR_DP_WAVE;
184        }
185        // core first, then add in shell
186        qr=x*rcore;
187        contr = rhocore-rhoshel;
188        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
189        vol = 4.0*pi/3.0*rcore*rcore*rcore;
190        f = vol*bes*contr;
191        //now the shell
192        qr=x*(rcore+thick);
193        contr = rhoshel-rhosolv;
194        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
195        vol = 4.0*pi/3.0*pow((rcore+thick),3);
196        f += vol*bes*contr;
197       
198        // normalize to the particle volume and rescale from [-1] to [cm-1]
199        //note that for the vesicle model, the volume is ONLY the shell volume
200        vol = 4.0*pi/3.0*(pow((rcore+thick),3)-pow(rcore,3));
201        f2 = f*f/vol*1.0e8;
202       
203        //scale if desired
204        f2 *= scale;
205        // then add in the background
206        f2 += bkg;
207       
208        p->result= f2;  //scale, and add in the background
209        return 0;
210}
211
212
213// scattering from a core shell sphere with a (Schulz) polydisperse core and constant shell thickness
214//
215int
216PolyCoreFormX(FitParamsPtr p)
217{
218        DOUBLE *dp;                             // Pointer to double precision wave data.
219        float *fp;                              // Pointer to single precision wave data.
220        DOUBLE pi;
221        DOUBLE scale,corrad,sig,zz,del,drho1,drho2,form,bkg;            //my local names
222        DOUBLE d, g ,h;
223        DOUBLE qq, x, y, c1, c2, c3, c4, c5, c6, c7, c8, c9, t1, t2, t3;
224        DOUBLE t4, t5, tb, cy, sy, tb1, tb2, tb3, c2y, zp1, zp2;
225        DOUBLE zp3,vpoly;
226        DOUBLE s2y, arg1, arg2, arg3, drh1, drh2;
227
228        if (p->waveHandle == NIL) {
229                SetNaN64(&p->result);
230                return NON_EXISTENT_WAVE;
231        }
232       
233        pi = 4.0*atan(1.0);
234        qq= p->x;
235       
236        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
237                case NT_FP32:
238                        fp= WaveData(p->waveHandle);
239                        scale = fp[0];
240                        corrad = fp[1];
241                        sig = fp[2];
242                        del = fp[3];
243                        drho1 = fp[4]-fp[5];            //core-shell
244                        drho2 = fp[5]-fp[6];            //shell-solvent
245                        bkg = fp[7];
246                                           
247                        break;
248                case NT_FP64:
249                        dp= WaveData(p->waveHandle);
250                        scale = dp[0];
251                        corrad = dp[1];
252                        sig = dp[2];
253                        del = dp[3];
254                        drho1 = dp[4]-dp[5];            //core-shell
255                        drho2 = dp[5]-dp[6];            //shell-solvent
256                        bkg = dp[7];
257                   
258                        break;
259                default:                                                                // We can't handle this wave data type.
260                        SetNaN64(&p->result);
261                        return REQUIRES_SP_OR_DP_WAVE;
262        }
263       
264        zz = (1.0/sig)*(1.0/sig) - 1.0;   
265       
266        h=qq;
267       
268        drh1 = drho1;
269        drh2 = drho2;
270        g = drh2 * -1. / drh1;
271        zp1 = zz + 1.;
272        zp2 = zz + 2.;
273        zp3 = zz + 3.;
274        vpoly = 4*pi/3*zp3*zp2/zp1/zp1*pow((corrad+del),3);
275
276
277        // remember that h is the passed in value of q for the calculation
278        y = h *del;
279        x = h *corrad;
280        d = atan(x * 2. / zp1);
281        arg1 = zp1 * d;
282        arg2 = zp2 * d;
283        arg3 = zp3 * d;
284        sy = sin(y);
285        cy = cos(y);
286        s2y = sin(y * 2.);
287        c2y = cos(y * 2.);
288        c1 = .5 - g * (cy + y * sy) + g * g * .5 * (y * y + 1.);
289        c2 = g * y * (g - cy);
290        c3 = (g * g + 1.) * .5 - g * cy;
291        c4 = g * g * (y * cy - sy) * (y * cy - sy) - c1;
292        c5 = g * 2. * sy * (1. - g * (y * sy + cy)) + c2;
293        c6 = c3 - g * g * sy * sy;
294        c7 = g * sy - g * .5 * g * (y * y + 1.) * s2y - c5;
295        c8 = c4 - .5 + g * cy - g * .5 * g * (y * y + 1.) * c2y;
296        c9 = g * sy * (1. - g * cy);
297
298        tb = log(zp1 * zp1 / (zp1 * zp1 + x * 4. * x));
299        tb1 = exp(zp1 * .5 * tb);
300        tb2 = exp(zp2 * .5 * tb);
301        tb3 = exp(zp3 * .5 * tb);
302
303        t1 = c1 + c2 * x + c3 * x * x * zp2 / zp1;
304        t2 = tb1 * (c4 * cos(arg1) + c7 * sin(arg1));
305        t3 = x * tb2 * (c5 * cos(arg2) + c8 * sin(arg2));
306        t4 = zp2 / zp1 * x * x * tb3 * (c6 * cos(arg3) + c9 * sin(arg3));
307        t5 = t1 + t2 + t3 + t4;
308        form = t5 * 16. * pi * pi * drh1 * drh1 / pow(qq,6);
309//      normalize by the average volume !!! corrected for polydispersity
310// and convert to cm-1
311        form /= vpoly;
312        form *= 1.0e8;
313        //Scale
314        form *= scale;
315        // then add in the background
316        form += bkg;
317       
318        p->result= (form);      //scale, and add in the background
319        return 0;
320}
321
322// scattering from a uniform sphere with a (Schulz) size distribution
323// structure factor effects are explicitly and correctly included.
324//
325int
326PolyHardSphereIntensityX(FitParamsPtr p)
327{
328        DOUBLE *dp;                             // Pointer to double precision wave data.
329        float *fp;                              // Pointer to single precision wave data.
330        DOUBLE pi;
331        DOUBLE rad,z2,phi,cont,bkg,sigma;               //my local names
332        DOUBLE mu,mu1,d1,d2,d3,d4,d5,d6,capd,rho;
333        DOUBLE ll,l1,bb,cc,chi,chi1,chi2,ee,t1,t2,t3,pp;
334        DOUBLE ka,zz,v1,v2,p1,p2;
335        DOUBLE h1,h2,h3,h4,e1,yy,y1,s1,s2,s3,hint1,hint2;
336        DOUBLE capl,capl1,capmu,capmu1,r3,pq;
337        DOUBLE ka2,r1,heff;
338        DOUBLE hh,k;
339
340        if (p->waveHandle == NIL) {
341                SetNaN64(&p->result);
342                return NON_EXISTENT_WAVE;
343        }
344       
345        pi = 4.0*atan(1.0);
346        k= p->x;
347       
348        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
349                case NT_FP32:
350                        fp= WaveData(p->waveHandle);
351                        rad = fp[0];            // radius (A)
352                        z2 = fp[1];             //polydispersity (0<z2<1)
353                        phi = fp[2];            // volume fraction (0<phi<1)
354                        cont = fp[3]*1.0e4;             // contrast (odd units)
355                        bkg = fp[4];
356                                                   
357                        break;
358                case NT_FP64:
359                        dp= WaveData(p->waveHandle);
360                        rad = dp[0];            // radius (A)
361                        z2 = dp[1];             //polydispersity (0<z2<1)
362                        phi = dp[2];            // volume fraction (0<phi<1)
363                        cont = dp[3]*1.0e4;             // contrast (odd units)
364                        bkg = dp[4];
365
366                        break;
367                default:                                                                // We can't handle this wave data type.
368                        SetNaN64(&p->result);
369                        return REQUIRES_SP_OR_DP_WAVE;
370        }
371       
372        sigma = 2*rad;
373       
374        zz=1.0/(z2*z2)-1.0;
375        bb = sigma/(zz+1.0);
376        cc = zz+1.0;
377
378//*c   Compute the number density by <r-cubed>, not <r> cubed*/
379        r1 = sigma/2.0;
380        r3 = r1*r1*r1;
381        r3 *= (zz+2.0)*(zz+3.0)/((zz+1.0)*(zz+1.0));
382        rho=phi/(1.3333333333*pi*r3);
383        t1 = rho*bb*cc;
384        t2 = rho*bb*bb*cc*(cc+1.0);
385        t3 = rho*bb*bb*bb*cc*(cc+1.0)*(cc+2.0);
386        capd = 1.0-pi*t3/6.0;
387//************
388        v1=1.0/(1.0+bb*bb*k*k);
389        v2=1.0/(4.0+bb*bb*k*k);
390        pp=pow(v1,(cc/2.0))*sin(cc*atan(bb*k));
391        p1=bb*cc*pow(v1,((cc+1.0)/2.0))*sin((cc+1.0)*atan(bb*k));
392        p2=cc*(cc+1.0)*bb*bb*pow(v1,((cc+2.0)/2.0))*sin((cc+2.0)*atan(bb*k));
393        mu=pow(2,cc)*pow(v2,(cc/2.0))*sin(cc*atan(bb*k/2.0));
394        mu1=pow(2,(cc+1.0))*bb*cc*pow(v2,((cc+1.0)/2.0))*sin((cc+1.0)*atan(k*bb/2.0));
395        s1=bb*cc;
396        s2=cc*(cc+1.0)*bb*bb;
397        s3=cc*(cc+1.0)*(cc+2.0)*bb*bb*bb;
398        chi=pow(v1,(cc/2.0))*cos(cc*atan(bb*k));
399        chi1=bb*cc*pow(v1,((cc+1.0)/2.0))*cos((cc+1.0)*atan(bb*k));
400        chi2=cc*(cc+1.0)*bb*bb*pow(v1,((cc+2.0)/2.0))*cos((cc+2.0)*atan(bb*k));
401        ll=pow(2,cc)*pow(v2,(cc/2.0))*cos(cc*atan(bb*k/2.0));
402        l1=pow(2,(cc+1.0))*bb*cc*pow(v2,((cc+1.0)/2.0))*cos((cc+1.0)*atan(k*bb/2.0));
403        d1=(pi/capd)*(2.0+(pi/capd)*(t3-(rho/k)*(k*s3-p2)));
404        d2=pow((pi/capd),2)*(rho/k)*(k*s2-p1);
405        d3=(-1.0)*pow((pi/capd),2)*(rho/k)*(k*s1-pp);
406        d4=(pi/capd)*(k-(pi/capd)*(rho/k)*(chi1-s1));
407        d5=pow((pi/capd),2)*((rho/k)*(chi-1.0)+0.5*k*t2);
408        d6=pow((pi/capd),2)*(rho/k)*(chi2-s2);
409       
410        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));
411        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;
412        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));
413        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;   
414
415        capl=2.0*pi*cont*rho/k/k/k*(pp-0.5*k*(s1+chi1));
416        capl1=2.0*pi*cont*rho/k/k/k*(p1-0.5*k*(s2+chi2));
417        capmu=2.0*pi*cont*rho/k/k/k*(1.0-chi-0.5*k*p1);
418        capmu1=2.0*pi*cont*rho/k/k/k*(s1-chi1-0.5*k*p2);
419 
420        h1=capl*(capl*(yy*d1-ee*d6)+capl1*(yy*d2-ee*d4)+capmu*(ee*d1+yy*d6)+capmu1*(ee*d2+yy*d4));
421        h2=capl1*(capl*(yy*d2-ee*d4)+capl1*(yy*d3-ee*d5)+capmu*(ee*d2+yy*d4)+capmu1*(ee*d3+yy*d5));
422        h3=capmu*(capl*(ee*d1+yy*d6)+capl1*(ee*d2+yy*d4)+capmu*(ee*d6-yy*d1)+capmu1*(ee*d4-yy*d2));
423        h4=capmu1*(capl*(ee*d2+yy*d4)+capl1*(ee*d3+yy*d5)+capmu*(ee*d4-yy*d2)+capmu1*(ee*d5-yy*d3));
424
425//*  This part computes the second integral in equation (1) of the paper.*/
426
427        hint1 = -2.0*(h1+h2+h3+h4)/(k*k*k*(ee*ee+yy*yy));
428 
429//*  This part computes the first integral in equation (1).  It also
430// generates the KC approximated effective structure factor.*/
431
432        pq=4.0*pi*cont*(sin(k*sigma/2.0)-0.5*k*sigma*cos(k*sigma/2.0));
433        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));
434 
435        ka=k*(sigma/2.0);
436//
437        hh=hint1+hint2;         // this is the model intensity
438//
439        heff=1.0+hint1/hint2;
440        ka2=ka*ka;
441//*
442//  heff is PY analytical solution for intensity divided by the
443//   form factor.  happ is the KC approximated effective S(q)
444 
445//*******************
446//  add in the background then return the intensity value
447               
448        p->result= (hh+bkg);    //scale, and add in the background
449        return 0;
450}
451
452// scattering from a uniform sphere with a (Schulz) size distribution, bimodal population
453// NO CROSS TERM IS ACCOUNTED FOR == DILUTE SOLUTION!!
454//
455int
456BimodalSchulzSpheresX(FitParamsPtr p)
457{
458        DOUBLE *dp;                             // Pointer to double precision wave data.
459        float *fp;                              // Pointer to single precision wave data.
460        DOUBLE x,pq;
461        DOUBLE scale,ravg,pd,bkg,rho,rhos;              //my local names
462        DOUBLE scale2,ravg2,pd2,rho2;           //my local names
463       
464        if (p->waveHandle == NIL) {
465                SetNaN64(&p->result);
466                return NON_EXISTENT_WAVE;
467        }
468       
469        x= p->x;
470       
471        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
472                case NT_FP32:
473                        fp= WaveData(p->waveHandle);
474                        scale = fp[0];
475                        ravg = fp[1];
476                        pd = fp[2];
477                        rho = fp[3];
478                        scale2 = fp[4];
479                        ravg2 = fp[5];
480                        pd2 = fp[6];
481                        rho2 = fp[7];
482                        rhos = fp[8];
483                        bkg = fp[9];
484                                           
485                        break;
486                case NT_FP64:
487                        dp= WaveData(p->waveHandle);
488                        scale = dp[0];
489                        ravg = dp[1];
490                        pd = dp[2];
491                        rho = dp[3];
492                        scale2 = dp[4];
493                        ravg2 = dp[5];
494                        pd2 = dp[6];
495                        rho2 = dp[7];
496                        rhos = dp[8];
497                        bkg = dp[9];
498                   
499                        break;
500                default:                                                                // We can't handle this wave data type.
501                        SetNaN64(&p->result);
502                        return REQUIRES_SP_OR_DP_WAVE;
503        }
504
505        pq = SchulzSphere_Fn( scale,  ravg,  pd,  rho,  rhos, x);
506        pq += SchulzSphere_Fn( scale2,  ravg2,  pd2,  rho2,  rhos, x);
507// add in the background
508        pq += bkg;
509               
510        p->result= pq;  //scale, and add in the background
511        return 0;
512}
513
514// scattering from a uniform sphere with a (Schulz) size distribution
515//
516int
517SchulzSpheresX(FitParamsPtr p)
518{
519        DOUBLE *dp;                             // Pointer to double precision wave data.
520        float *fp;                              // Pointer to single precision wave data.
521        DOUBLE x,pq;
522        DOUBLE scale,ravg,pd,bkg,rho,rhos;              //my local names
523       
524        if (p->waveHandle == NIL) {
525                SetNaN64(&p->result);
526                return NON_EXISTENT_WAVE;
527        }
528       
529        x= p->x;
530       
531        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
532                case NT_FP32:
533                        fp= WaveData(p->waveHandle);
534                        scale = fp[0];
535                        ravg = fp[1];
536                        pd = fp[2];
537                        rho = fp[3];
538                        rhos = fp[4];
539                        bkg = fp[5];
540                                           
541                        break;
542                case NT_FP64:
543                        dp= WaveData(p->waveHandle);
544                        scale = dp[0];
545                        ravg = dp[1];
546                        pd = dp[2];
547                        rho = dp[3];
548                        rhos = dp[4];
549                        bkg = dp[5];
550                   
551                        break;
552                default:                                                                // We can't handle this wave data type.
553                        SetNaN64(&p->result);
554                        return REQUIRES_SP_OR_DP_WAVE;
555        }
556        pq = SchulzSphere_Fn( scale,  ravg,  pd,  rho,  rhos, x);
557// add in the background
558        pq += bkg;
559               
560        p->result= pq;  //scale, and add in the background
561        return 0;
562}
563
564// calculates everything but the background
565DOUBLE
566SchulzSphere_Fn(DOUBLE scale, DOUBLE ravg, DOUBLE pd, DOUBLE rho, DOUBLE rhos, DOUBLE x)
567{
568        DOUBLE zp1,zp2,zp3,zp4,zp5,zp6,zp7,vpoly;
569        DOUBLE aa,at1,at2,rt1,rt2,rt3,t1,t2,t3;
570        DOUBLE v1,v2,v3,g1,pq,pi,delrho,zz;
571       
572        pi = 4.0*atan(1.0);
573        delrho = rho-rhos;
574        zz = (1.0/pd)*(1.0/pd) - 1.0;   
575       
576        zp1 = zz + 1.0;
577        zp2 = zz + 2.0;
578        zp3 = zz + 3.0;
579        zp4 = zz + 4.0;
580        zp5 = zz + 5.0;
581        zp6 = zz + 6.0;
582        zp7 = zz + 7.0;
583//
584        aa = (zz+1)/x/ravg;
585
586        at1 = atan(1.0/aa);
587        at2 = atan(2.0/aa);
588//
589//  calculations are performed to avoid  large # errors
590// - trick is to propogate the a^(z+7) term through the g1
591//
592        t1 = zp7*log10(aa) - zp1/2.0*log10(aa*aa+4.0);
593        t2 = zp7*log10(aa) - zp3/2.0*log10(aa*aa+4.0);
594        t3 = zp7*log10(aa) - zp2/2.0*log10(aa*aa+4.0);
595//      print t1,t2,t3
596        rt1 = pow(10,t1);
597        rt2 = pow(10,t2);
598        rt3 = pow(10,t3);
599        v1 = pow(aa,6) - rt1*cos(zp1*at2);
600        v2 = zp1*zp2*( pow(aa,4) + rt2*cos(zp3*at2) );
601        v3 = -2.0*zp1*rt3*sin(zp2*at2);
602        g1 = (v1+v2+v3);
603       
604        pq = log10(g1) - 6.0*log10(zp1) + 6.0*log10(ravg);
605        pq = pow(10,pq)*8*pi*pi*delrho*delrho;
606       
607//
608// beta factor is not used here, but could be for the
609// decoupling approximation
610//
611//      g11 = g1
612//      gd = -zp7*log(aa)
613//      g1 = log(g11) + gd
614//                       
615//      t1 = zp1*at1
616//      t2 = zp2*at1
617//      g2 = sin( t1 ) - zp1/sqrt(aa*aa+1)*cos( t2 )
618//      g22 = g2*g2
619//      beta = zp1*log(aa) - zp1*log(aa*aa+1) - g1 + log(g22)
620//      beta = 2*alog(beta)
621       
622//re-normalize by the average volume
623        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*ravg*ravg*ravg;
624        pq /= vpoly;
625//scale, convert to cm^-1
626        pq *= scale * 1.0e8;
627   
628    return(pq);
629}
630
631// scattering from a uniform sphere with a rectangular size distribution
632//
633int
634PolyRectSpheresX(FitParamsPtr p)
635{
636        DOUBLE *dp;                             // Pointer to double precision wave data.
637        float *fp;                              // Pointer to single precision wave data.
638        DOUBLE pi,x;
639        DOUBLE scale,rad,pd,cont,bkg;           //my local names
640        DOUBLE inten,h1,qw,qr,width,sig,averad3;
641       
642        if (p->waveHandle == NIL) {
643                SetNaN64(&p->result);
644                return NON_EXISTENT_WAVE;
645        }
646       
647        pi = 4.0*atan(1.0);
648        x= p->x;
649       
650        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
651                case NT_FP32:
652                        fp= WaveData(p->waveHandle);
653                        scale = fp[0];
654                        rad = fp[1];            // radius (A)
655                        pd = fp[2];             //polydispersity of rectangular distribution
656                        cont = fp[3];           // contrast (A^-2)
657                        bkg = fp[4];
658                                           
659                        break;
660                case NT_FP64:
661                        dp= WaveData(p->waveHandle);
662                        scale = dp[0];
663                        rad = dp[1];            // radius (A)
664                        pd = dp[2];             //polydispersity of rectangular distribution
665                        cont = dp[3];           // contrast (A^-2)
666                        bkg = dp[4];
667
668                   
669                        break;
670                default:                                                                // We can't handle this wave data type.
671                        SetNaN64(&p->result);
672                        return REQUIRES_SP_OR_DP_WAVE;
673        }
674       
675       // as usual, poly = sig/ravg
676        // for the rectangular distribution, sig = width/sqrt(3)
677        // width is the HALF- WIDTH of the rectangular distrubution
678       
679        sig = pd*rad;
680        width = sqrt(3.0)*sig;
681       
682        //x is the q-value
683        qw = x*width;
684        qr = x*rad;
685        h1 = -0.5*qw + qr*qr*qw + (qw*qw*qw)/3.0;
686        h1 -= 5.0/2.0*cos(2*qr)*sin(qw)*cos(qw);
687        h1 += 0.5*qr*qr*cos(2*qr)*sin(2*qw);
688        h1 += 0.5*qw*qw*cos(2*qr)*sin(2*qw);
689        h1 += qw*qr*sin(2*qr)*cos(2*qw);
690        h1 += 3.0*qw*(cos(qr)*cos(qw))*(cos(qr)*cos(qw));
691        h1+= 3.0*qw*(sin(qr)*sin(qw))*(sin(qr)*sin(qw));
692        h1 -= 6.0*qr*cos(qr)*sin(qr)*cos(qw)*sin(qw);
693       
694        // calculate P(q) = <f^2>
695        inten = 8.0*pi*pi*cont*cont/width/pow(x,7)*h1;
696       
697        // beta(q) would be calculated as 2/width/x/h1*h2*h2
698        // with
699        // h2 = 2*sin(x*rad)*sin(x*width)-x*rad*cos(x*rad)*sin(x*width)-x*width*sin(x*rad)*cos(x*width)
700
701        // normalize to the average volume
702        // <R^3> = ravg^3*(1+3*pd^2)
703        // or... "zf"  = (1 + 3*p^2), which will be greater than one
704       
705        averad3 =  rad*rad*rad*(1.0+3.0*pd*pd);
706        inten /= 4.0*pi/3.0*averad3;
707        //resacle to 1/cm
708        inten *= 1.0e8;
709        //scale the result
710        inten *= scale;
711        // then add in the background
712        inten += bkg;
713
714        p->result= inten;       //scale, and add in the background
715        return 0;
716}
717
718
719// scattering from a uniform sphere with a Gaussian size distribution
720//
721int
722GaussPolySphereX(FitParamsPtr p)
723{
724        DOUBLE *dp;                             // Pointer to double precision wave data.
725        float *fp;                              // Pointer to single precision wave data.
726        DOUBLE pi,x;
727        DOUBLE scale,rad,pd,sig,rho,rhos,bkg,delrho;            //my local names
728        DOUBLE va,vb,zi,yy,summ,inten;
729        int nord=20,ii;
730       
731        if (p->waveHandle == NIL) {
732                SetNaN64(&p->result);
733                return NON_EXISTENT_WAVE;
734        }
735       
736        pi = 4.0*atan(1.0);
737        x= p->x;
738       
739        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
740                case NT_FP32:
741                        fp= WaveData(p->waveHandle);
742                        scale=fp[0];
743                        rad=fp[1];
744                        pd=fp[2];
745                        sig=pd*rad;
746                        rho=fp[3];
747                        rhos=fp[4];
748                        delrho=rho-rhos;
749                        bkg=fp[5];
750                                                   
751                        break;
752                case NT_FP64:
753                        dp= WaveData(p->waveHandle);
754                        scale=dp[0];
755                        rad=dp[1];
756                        pd=dp[2];
757                        sig=pd*rad;
758                        rho=dp[3];
759                        rhos=dp[4];
760                        delrho=rho-rhos;
761                        bkg=dp[5];
762                   
763                        break;
764                default:                                                                // We can't handle this wave data type.
765                        SetNaN64(&p->result);
766                        return REQUIRES_SP_OR_DP_WAVE;
767        }
768       
769        va = -4.0*sig + rad;
770        if (va<0) {
771                va=0;           //to avoid numerical error when  va<0 (-ve q-value)
772        }
773        vb = 4.0*sig +rad;
774       
775        summ = 0.0;             // initialize integral
776        for(ii=0;ii<nord;ii+=1) {
777                // calculate Gauss points on integration interval (r-value for evaluation)
778                zi = ( Gauss20Z[ii]*(vb-va) + vb + va )/2.0;
779                // calculate sphere scattering
780                //return(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr)); pass qr
781                yy = F_func(x*zi)*(4.0*pi/3.0*zi*zi*zi)*delrho;
782                yy *= yy;
783                yy *= Gauss20Wt[ii] *  Gauss_distr(sig,rad,zi);
784               
785                summ += yy;             //add to the running total of the quadrature
786        }
787// calculate value of integral to return
788        inten = (vb-va)/2.0*summ;
789       
790        //re-normalize by polydisperse sphere volume
791        inten /= (4.0*pi/3.0*rad*rad*rad)*(1.0+3.0*pd*pd);
792       
793        inten *= 1.0e8;
794        inten *= scale;
795        inten += bkg;
796               
797        p->result= inten;       //scale, and add in the background
798        return 0;
799}
800
801// scattering from a uniform sphere with a LogNormal size distribution
802//
803int
804LogNormalPolySphereX(FitParamsPtr p)
805{
806        DOUBLE *dp;                             // Pointer to double precision wave data.
807        float *fp;                              // Pointer to single precision wave data.
808        DOUBLE pi,x;
809        DOUBLE scale,rad,sig,rho,rhos,bkg,delrho,mu,r3;         //my local names
810        DOUBLE va,vb,zi,yy,summ,inten;
811        int nord=76,ii;
812       
813        if (p->waveHandle == NIL) {
814                SetNaN64(&p->result);
815                return NON_EXISTENT_WAVE;
816        }
817       
818        pi = 4.0*atan(1.0);
819        x= p->x;
820       
821        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
822                case NT_FP32:
823                        fp= WaveData(p->waveHandle);
824                        scale=fp[0];
825                        rad=fp[1];              //rad is the median radius
826                        mu = log(fp[1]);
827                        sig=fp[2];
828                        rho=fp[3];
829                        rhos=fp[4];
830                        delrho=rho-rhos;
831                        bkg=fp[5];
832                                                   
833                        break;
834                case NT_FP64:
835                        dp= WaveData(p->waveHandle);
836                        scale=dp[0];
837                        rad=dp[1];              //rad is the median radius
838                        mu = log(dp[1]);
839                        sig=dp[2];
840                        rho=dp[3];
841                        rhos=dp[4];
842                        delrho=rho-rhos;
843                        bkg=dp[5];
844                   
845                        break;
846                default:                                                                // We can't handle this wave data type.
847                        SetNaN64(&p->result);
848                        return REQUIRES_SP_OR_DP_WAVE;
849        }
850       
851        va = -3.5*sig + mu;
852        va = exp(va);
853        if (va<0) {
854                va=0;           //to avoid numerical error when  va<0 (-ve q-value)
855        }
856        vb = 3.5*sig*(1.0+sig) +mu;
857        vb = exp(vb);
858       
859        summ = 0.0;             // initialize integral
860        for(ii=0;ii<nord;ii+=1) {
861                // calculate Gauss points on integration interval (r-value for evaluation)
862                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0;
863                // calculate sphere scattering
864                //return(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr)); pass qr
865                yy = F_func(x*zi)*(4.0*pi/3.0*zi*zi*zi)*delrho;
866                yy *= yy;
867                yy *= Gauss76Wt[ii] *  LogNormal_distr(sig,mu,zi);
868               
869                summ += yy;             //add to the running total of the quadrature
870        }
871// calculate value of integral to return
872        inten = (vb-va)/2.0*summ;
873       
874        //re-normalize by polydisperse sphere volume
875        r3 = exp(3.0*mu + 9.0/2.0*sig*sig);             // <R^3> directly
876        inten /= (4.0*pi/3.0*r3);               //polydisperse volume
877               
878        inten *= 1.0e8;
879        inten *= scale;
880        inten += bkg;
881               
882        p->result= inten;       //scale, and add in the background
883        return 0;
884}
885
886static DOUBLE
887LogNormal_distr(DOUBLE sig, DOUBLE mu, DOUBLE pt)
888{       
889        DOUBLE retval,pi;
890       
891        pi = 4.0*atan(1.0);
892        retval = (1/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig );
893        return(retval);
894}
895
896static DOUBLE
897Gauss_distr(DOUBLE sig, DOUBLE avg, DOUBLE pt)
898{       
899        DOUBLE retval,Pi;
900       
901        Pi = 4.0*atan(1.0);
902        retval = (1.0/ (sig*sqrt(2.0*Pi)) )*exp(-(avg-pt)*(avg-pt)/sig/sig/2.0);
903        return(retval);
904}
905
906// scattering from a core shell sphere with a (Schulz) polydisperse core and constant ratio (shell thickness)/(core radius)
907// - the polydispersity is of the WHOLE sphere
908//
909int
910PolyCoreShellRatioX(FitParamsPtr p)
911{
912        DOUBLE *dp;                             // Pointer to double precision wave data.
913        float *fp;                              // Pointer to single precision wave data.
914        DOUBLE pi,x;
915        DOUBLE scale,corrad,thick,shlrad,pp,drho1,drho2,sig,zz,bkg;             //my local names
916        DOUBLE sld1,sld2,sld3,zp1,zp2,zp3,vpoly;
917        DOUBLE pi43,c1,c2,form,volume,arg1,arg2;
918
919        if (p->waveHandle == NIL) {
920                SetNaN64(&p->result);
921                return NON_EXISTENT_WAVE;
922        }
923       
924        pi = 4.0*atan(1.0);
925        x= p->x;
926       
927        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
928                case NT_FP32:
929                        fp= WaveData(p->waveHandle);
930                        scale = fp[0];
931                        corrad = fp[1];
932                        thick = fp[2];
933                        sig = fp[3];
934                        sld1 = fp[4];
935                        sld2 = fp[5];
936                        sld3 = fp[6];
937                        bkg = fp[7];
938                                           
939                        break;
940                case NT_FP64:
941                        dp= WaveData(p->waveHandle);
942                        scale = dp[0];
943                        corrad = dp[1];
944                        thick = dp[2];
945                        sig = dp[3];
946                        sld1 = dp[4];
947                        sld2 = dp[5];
948                        sld3 = dp[6];
949                        bkg = dp[7];
950                   
951                        break;
952                default:                                                                // We can't handle this wave data type.
953                        SetNaN64(&p->result);
954                        return REQUIRES_SP_OR_DP_WAVE;
955        }
956       
957        zz = (1.0/sig)*(1.0/sig) - 1.0;   
958        shlrad = corrad + thick;
959        drho1 = sld1-sld2;              //core-shell
960        drho2 = sld2-sld3;              //shell-solvent
961        zp1 = zz + 1.;
962        zp2 = zz + 2.;
963        zp3 = zz + 3.;
964        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((corrad+thick),3);
965       
966        // the beta factor is not calculated
967        // the calculated form factor <f^2> has units [length^2]
968        // and must be multiplied by number density [l^-3] and the correct unit
969        // conversion to get to absolute scale
970       
971        pi43=4.0/3.0*pi;
972        pp=corrad/shlrad;
973        volume=pi43*shlrad*shlrad*shlrad;
974        c1=drho1*volume;
975        c2=drho2*volume;
976       
977        arg1 = x*shlrad*pp;
978        arg2 = x*shlrad;
979               
980        form=pow(pp,6)*c1*c1*fnt2(arg1,zz);
981        form += c2*c2*fnt2(arg2,zz);
982        form += 2.0*c1*c2*fnt3(arg2,pp,zz);
983       
984        //convert the result to [cm^-1]
985       
986        //scale the result
987        // - divide by the polydisperse volume, mult by 10^8
988        form  /= vpoly;
989        form *= 1.0e8;
990        form *= scale;
991
992        //add in the background
993        form += bkg;
994       
995        p->result= form;        //scale, and add in the background
996        return 0;
997}
998
999//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1000//c
1001//c      function fnt2(y,z)
1002//c
1003DOUBLE
1004fnt2(DOUBLE yy, DOUBLE zz)
1005{
1006        DOUBLE z1,z2,z3,u,ww,term1,term2,term3,ans;
1007       
1008        z1=zz+1.0;
1009        z2=zz+2.0;
1010        z3=zz+3.0;
1011        u=yy/z1;
1012        ww=atan(2.0*u);
1013        term1=cos(z1*ww)/pow((1.0+4.0*u*u),(z1/2.0));
1014        term2=2.0*yy*sin(z2*ww)/pow((1.0+4.0*u*u),(z2/2.0));
1015        term3=1.0+cos(z3*ww)/pow((1.0+4.0*u*u),(z3/2.0));
1016        ans=(4.50/z1/pow(yy,6))*(z1*(1.0-term1-term2)+yy*yy*z2*term3);
1017
1018        return(ans);
1019}
1020
1021//cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1022//c
1023//c      function fnt3(y,p,z)
1024//c
1025DOUBLE
1026fnt3(DOUBLE yy, DOUBLE pp, DOUBLE zz)
1027{       
1028        DOUBLE z1,z2,z3,yp,yn,up,un,vp,vn,term1,term2,term3,term4,term5,term6,ans;
1029       
1030        z1=zz+1.0;
1031        z2=zz+2.0;
1032        z3=zz+3.0;
1033        yp=(1.0+pp)*yy;
1034        yn=(1.0-pp)*yy;
1035        up=yp/z1;
1036        un=yn/z1;
1037        vp=atan(up);
1038        vn=atan(un);
1039        term1=cos(z1*vn)/pow((1.0+un*un),(z1/2.0));
1040        term2=cos(z1*vp)/pow((1.0+up*up),(z1/2.0));
1041        term3=cos(z3*vn)/pow((1.0+un*un),(z3/2.0));
1042        term4=cos(z3*vp)/pow((1.0+up*up),(z3/2.0));
1043        term5=yn*sin(z2*vn)/pow((1.0+un*un),(z2/2.0));
1044        term6=yp*sin(z2*vp)/pow((1.0+up*up),(z2/2.0));
1045        ans=4.5/z1/pow(yy,6);
1046        ans *=(z1*(term1-term2)+yy*yy*pp*z2*(term3+term4)+z1*(term5-term6));
1047     
1048        return(ans);
1049}
1050
1051// scattering from a a binary population of hard spheres, 3 partial structure factors
1052// are properly accounted for...
1053//       Input (fitting) variables are:
1054//      larger sphere radius(angstroms) = guess[0]
1055//      smaller sphere radius (A) = w[1]
1056//      number fraction of larger spheres = guess[2]
1057//      total volume fraction of spheres = guess[3]
1058//      size ratio, alpha(0<a<1) = derived
1059//      SLD(A-2) of larger particle = guess[4]
1060//      SLD(A-2) of smaller particle = guess[5]
1061//      SLD(A-2) of the solvent = guess[6]
1062//      background = guess[7]
1063int
1064BinaryHSX(FitParamsPtr p)
1065{
1066        DOUBLE *dp;                             // Pointer to double precision wave data.
1067        float *fp;                              // Pointer to single precision wave data.
1068        DOUBLE x,pi;
1069        DOUBLE r2,r1,nf2,phi,aa,rho2,rho1,rhos,inten,bgd;               //my local names
1070        DOUBLE psf11,psf12,psf22;
1071        DOUBLE phi1,phi2,phr,a3;
1072        DOUBLE v1,v2,n1,n2,qr1,qr2,b1,b2;
1073        int err;
1074       
1075        if (p->waveHandle == NIL) {
1076                SetNaN64(&p->result);
1077                return NON_EXISTENT_WAVE;
1078        }
1079       
1080        pi = 4.0*atan(1.0);
1081        x= p->x;
1082       
1083        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
1084                case NT_FP32:
1085                        fp= WaveData(p->waveHandle);
1086                        r2 = fp[0];
1087                        r1 = fp[1];
1088                        phi2 = fp[2];
1089                        phi1 = fp[3];
1090                        rho2 = fp[4];
1091                        rho1 = fp[5];
1092                        rhos = fp[6];
1093                        bgd = fp[7];
1094                   
1095                        break;
1096                case NT_FP64:
1097                        dp= WaveData(p->waveHandle);
1098                        r2 = dp[0];
1099                        r1 = dp[1];
1100                        phi2 = dp[2];
1101                        phi1 = dp[3];
1102                        rho2 = dp[4];
1103                        rho1 = dp[5];
1104                        rhos = dp[6];
1105                        bgd = dp[7];
1106                                               
1107                        break;
1108                default:                                                                // We can't handle this wave data type.
1109                        SetNaN64(&p->result);
1110                        return REQUIRES_SP_OR_DP_WAVE;
1111        }
1112        phi = phi1 + phi2;
1113        aa = r1/r2;
1114        //calculate the number fraction of larger spheres (eqn 2 in reference)
1115        a3=aa*aa*aa;
1116        phr=phi2/phi;
1117        nf2 = phr*a3/(1.0-phr+phr*a3);
1118        // calculate the PSF's here
1119        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
1120       
1121        // /* do form factor calculations  */
1122       
1123        v1 = 4.0*pi/3.0*r1*r1*r1;
1124        v2 = 4.0*pi/3.0*r2*r2*r2;
1125
1126        n1 = phi1/v1;
1127        n2 = phi2/v2;
1128
1129        qr1 = r1*x;
1130        qr2 = r2*x;
1131       
1132        b1 = r1*r1*r1*(rho1-rhos)*4.0*pi*(sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1;
1133        b2 = r2*r2*r2*(rho2-rhos)*4.0*pi*(sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2;
1134        inten = n1*b1*b1*psf11;
1135        inten += sqrt(n1*n2)*2.0*b1*b2*psf12;
1136        inten += n2*b2*b2*psf22;
1137///* convert I(1/A) to (1/cm)  */
1138        inten *= 1.0e8;
1139       
1140        inten += bgd;
1141       
1142        p->result= inten;       //scale, and add in the background
1143        return 0;
1144}
1145
1146int
1147BinaryHS_PSF11X(FitParamsPtr p)
1148{
1149        DOUBLE *dp;                             // Pointer to double precision wave data.
1150        float *fp;                              // Pointer to single precision wave data.
1151        DOUBLE x,pi;
1152        DOUBLE r2,r1,nf2,phi,aa,rho2,rho1,rhos,bgd;             //my local names
1153        DOUBLE psf11,psf12,psf22;
1154        DOUBLE phi1,phi2,phr,a3;
1155        int err;
1156       
1157        if (p->waveHandle == NIL) {
1158                SetNaN64(&p->result);
1159                return NON_EXISTENT_WAVE;
1160        }
1161       
1162        pi = 4.0*atan(1.0);
1163        x= p->x;
1164       
1165        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
1166                case NT_FP32:
1167                        fp= WaveData(p->waveHandle);
1168                        r2 = fp[0];
1169                        r1 = fp[1];
1170                        phi2 = fp[2];
1171                        phi1 = fp[3];
1172                        rho2 = fp[4];
1173                        rho1 = fp[5];
1174                        rhos = fp[6];
1175                        bgd = fp[7];
1176                   
1177                        break;
1178                case NT_FP64:
1179                        dp= WaveData(p->waveHandle);
1180                        r2 = dp[0];
1181                        r1 = dp[1];
1182                        phi2 = dp[2];
1183                        phi1 = dp[3];
1184                        rho2 = dp[4];
1185                        rho1 = dp[5];
1186                        rhos = dp[6];
1187                        bgd = dp[7];
1188                                               
1189                        break;
1190                default:                                                                // We can't handle this wave data type.
1191                        SetNaN64(&p->result);
1192                        return REQUIRES_SP_OR_DP_WAVE;
1193        }
1194        phi = phi1 + phi2;
1195        aa = r1/r2;
1196        //calculate the number fraction of larger spheres (eqn 2 in reference)
1197        a3=aa*aa*aa;
1198        phr=phi2/phi;
1199        nf2 = phr*a3/(1.0-phr+phr*a3);
1200        // calculate the PSF's here
1201        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
1202       
1203        p->result= psf11;       //scale, and add in the background
1204        return 0;
1205}
1206
1207int
1208BinaryHS_PSF12X(FitParamsPtr p)
1209{
1210        DOUBLE *dp;                             // Pointer to double precision wave data.
1211        float *fp;                              // Pointer to single precision wave data.
1212        DOUBLE x,pi;
1213        DOUBLE r2,r1,nf2,phi,aa,rho2,rho1,rhos,bgd;             //my local names
1214        DOUBLE psf11,psf12,psf22;
1215        DOUBLE phi1,phi2,phr,a3;
1216        int err;
1217       
1218        if (p->waveHandle == NIL) {
1219                SetNaN64(&p->result);
1220                return NON_EXISTENT_WAVE;
1221        }
1222       
1223        pi = 4.0*atan(1.0);
1224        x= p->x;
1225       
1226        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
1227                case NT_FP32:
1228                        fp= WaveData(p->waveHandle);
1229                        r2 = fp[0];
1230                        r1 = fp[1];
1231                        phi2 = fp[2];
1232                        phi1 = fp[3];
1233                        rho2 = fp[4];
1234                        rho1 = fp[5];
1235                        rhos = fp[6];
1236                        bgd = fp[7];
1237                   
1238                        break;
1239                case NT_FP64:
1240                        dp= WaveData(p->waveHandle);
1241                        r2 = dp[0];
1242                        r1 = dp[1];
1243                        phi2 = dp[2];
1244                        phi1 = dp[3];
1245                        rho2 = dp[4];
1246                        rho1 = dp[5];
1247                        rhos = dp[6];
1248                        bgd = dp[7];
1249                                               
1250                        break;
1251                default:                                                                // We can't handle this wave data type.
1252                        SetNaN64(&p->result);
1253                        return REQUIRES_SP_OR_DP_WAVE;
1254        }
1255        phi = phi1 + phi2;
1256        aa = r1/r2;
1257        //calculate the number fraction of larger spheres (eqn 2 in reference)
1258        a3=aa*aa*aa;
1259        phr=phi2/phi;
1260        nf2 = phr*a3/(1.0-phr+phr*a3);
1261        // calculate the PSF's here
1262        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
1263       
1264        p->result= psf12;       //scale, and add in the background
1265        return 0;
1266}
1267
1268int
1269BinaryHS_PSF22X(FitParamsPtr p)
1270{
1271        DOUBLE *dp;                             // Pointer to double precision wave data.
1272        float *fp;                              // Pointer to single precision wave data.
1273        DOUBLE x,pi;
1274        DOUBLE r2,r1,nf2,phi,aa,rho2,rho1,rhos,bgd;             //my local names
1275        DOUBLE psf11,psf12,psf22;
1276        DOUBLE phi1,phi2,phr,a3;
1277        int err;
1278       
1279        if (p->waveHandle == NIL) {
1280                SetNaN64(&p->result);
1281                return NON_EXISTENT_WAVE;
1282        }
1283       
1284        pi = 4.0*atan(1.0);
1285        x= p->x;
1286       
1287        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
1288                case NT_FP32:
1289                        fp= WaveData(p->waveHandle);
1290                        r2 = fp[0];
1291                        r1 = fp[1];
1292                        phi2 = fp[2];
1293                        phi1 = fp[3];
1294                        rho2 = fp[4];
1295                        rho1 = fp[5];
1296                        rhos = fp[6];
1297                        bgd = fp[7];
1298                   
1299                        break;
1300                case NT_FP64:
1301                        dp= WaveData(p->waveHandle);
1302                        r2 = dp[0];
1303                        r1 = dp[1];
1304                        phi2 = dp[2];
1305                        phi1 = dp[3];
1306                        rho2 = dp[4];
1307                        rho1 = dp[5];
1308                        rhos = dp[6];
1309                        bgd = dp[7];
1310                                               
1311                        break;
1312                default:                                                                // We can't handle this wave data type.
1313                        SetNaN64(&p->result);
1314                        return REQUIRES_SP_OR_DP_WAVE;
1315        }
1316        phi = phi1 + phi2;
1317        aa = r1/r2;
1318        //calculate the number fraction of larger spheres (eqn 2 in reference)
1319        a3=aa*aa*aa;
1320        phr=phi2/phi;
1321        nf2 = phr*a3/(1.0-phr+phr*a3);
1322        // calculate the PSF's here
1323        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
1324       
1325        p->result= psf22;       //scale, and add in the background
1326        return 0;
1327}
1328
1329int
1330ashcroft(DOUBLE qval, DOUBLE r2, DOUBLE nf2, DOUBLE aa, DOUBLE phi, DOUBLE *s11, DOUBLE *s22, DOUBLE *s12)
1331{
1332//      variable qval,r2,nf2,aa,phi,&s11,&s22,&s12
1333
1334//   calculate constant terms
1335        DOUBLE s1,s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4;
1336        DOUBLE a1,a2i,a2,b1,b2,b12,gm1,gm12;
1337        DOUBLE err=0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3;
1338        DOUBLE c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13;
1339        DOUBLE t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1;
1340
1341        s2 = 2.0*r2;
1342        s1 = aa*s2;
1343        v = phi;
1344        a3 = aa*aa*aa;
1345        v1=((1.-nf2)*a3/(nf2+(1.-nf2)*a3))*v;
1346        v2=(nf2/(nf2+(1.-nf2)*a3))*v;
1347        g11=((1.+.5*v)+1.5*v2*(aa-1.))/(1.-v)/(1.-v);
1348        g22=((1.+.5*v)+1.5*v1*(1./aa-1.))/(1.-v)/(1.-v);
1349        g12=((1.+.5*v)+1.5*(1.-aa)*(v1-v2)/(1.+aa))/(1.-v)/(1.-v);
1350        wmv = 1/(1.-v);
1351        wmv3 = wmv*wmv*wmv;
1352        wmv4 = wmv*wmv3;
1353        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;
1354        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;
1355        a2=a2i/a3;
1356        b1=-6.*(v1*g11*g11+.25*v2*(1.+aa)*(1.+aa)*aa*g12*g12);
1357        b2=-6.*(v2*g22*g22+.25*v1/a3*(1.+aa)*(1.+aa)*g12*g12);
1358        b12=-3.*aa*(1.+aa)*(v1*g11/aa/aa+v2*g22)*g12;
1359        gm1=(v1*a1+a3*v2*a2)*.5;
1360        gm12=2.*gm1*(1.-aa)/aa;
1361//c 
1362//c   calculate the direct correlation functions and print results
1363//c
1364//      do 20 j=1,npts
1365
1366        yy=qval*s2;
1367//c   calculate direct correlation functions
1368//c   ----c11
1369        ay=aa*yy;
1370        ay2 = ay*ay;
1371        ay3 = ay*ay*ay;
1372        t1=a1*(sin(ay)-ay*cos(ay));
1373        t2=b1*(2.*ay*sin(ay)-(ay2-2.)*cos(ay)-2.)/ay;
1374        t3=gm1*((4.*ay*ay2-24.*ay)*sin(ay)-(ay2*ay2-12.*ay2+24.)*cos(ay)+24.)/ay3;
1375        f11=24.*v1*(t1+t2+t3)/ay3;
1376
1377//c ------c22
1378        y2=yy*yy;
1379        y3=yy*y2;
1380        tt1=a2*(sin(yy)-yy*cos(yy));
1381        tt2=b2*(2.*yy*sin(yy)-(y2-2.)*cos(yy)-2.)/yy;
1382        tt3=gm1*((4.*y3-24.*yy)*sin(yy)-(y2*y2-12.*y2+24.)*cos(yy)+24.)/ay3;
1383        f22=24.*v2*(tt1+tt2+tt3)/y3;
1384
1385//c   -----c12
1386        yl=.5*yy*(1.-aa);
1387        yl3=yl*yl*yl;
1388        wma3 = (1.-aa)*(1.-aa)*(1.-aa);
1389        y1=aa*yy;
1390        y13 = y1*y1*y1;
1391        ttt1=3.*wma3*v*sqrt(nf2)*sqrt(1.-nf2)*a1*(sin(yl)-yl*cos(yl))/((nf2+(1.-nf2)*a3)*yl3);
1392        t21=b12*(2.*y1*cos(y1)+(y1*y1-2.)*sin(y1));
1393        t22=gm12*((3.*y1*y1-6.)*cos(y1)+(y1*y1*y1-6.*y1)*sin(y1)+6.)/y1;
1394        t23=gm1*((4.*y13-24.*y1)*cos(y1)+(y13*y1-12.*y1*y1+24.)*sin(y1))/(y1*y1);
1395        t31=b12*(2.*y1*sin(y1)-(y1*y1-2.)*cos(y1)-2.);
1396        t32=gm12*((3.*y1*y1-6.)*sin(y1)-(y1*y1*y1-6.*y1)*cos(y1))/y1;
1397        t33=gm1*((4.*y13-24.*y1)*sin(y1)-(y13*y1-12.*y1*y1+24.)*cos(y1)+24.)/(y1*y1);
1398        t41=cos(yl)*((sin(y1)-y1*cos(y1))/(y1*y1) + (1.-aa)/(2.*aa)*(1.-cos(y1))/y1);
1399        t42=sin(yl)*((cos(y1)+y1*sin(y1)-1.)/(y1*y1) + (1.-aa)/(2.*aa)*sin(y1)/y1);
1400        ttt2=sin(yl)*(t21+t22+t23)/(y13*y1);
1401        ttt3=cos(yl)*(t31+t32+t33)/(y13*y1);
1402        ttt4=a1*(t41+t42)/y1;
1403        f12=ttt1+24.*v*sqrt(nf2)*sqrt(1.-nf2)*a3*(ttt2+ttt3+ttt4)/(nf2+(1.-nf2)*a3);
1404
1405        c11=f11;
1406        c22=f22;
1407        c12=f12;
1408        *s11=1./(1.+c11-(c12)*c12/(1.+c22));
1409        *s22=1./(1.+c22-(c12)*c12/(1.+c11)); 
1410        *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));   
1411
1412        return(err);
1413}
1414
1415
1416
1417/*
1418// calculates the scattering from a spherical particle made up of a core (aqueous) surrounded
1419// by N spherical layers, each of which is a PAIR of shells, solvent + surfactant since there
1420//must always be a surfactant layer on the outside
1421//
1422// bragg peaks arise naturally from the periodicity of the sample
1423// resolution smeared version gives he most appropriate view of the model
1424
1425        Warning:
1426                The call to WaveData() below returns a pointer to the middle
1427                of an unlocked Macintosh handle. In the unlikely event that your
1428                calculations could cause memory to move, you should copy the coefficient
1429                values to local variables or an array before such operations.
1430*/
1431int
1432MultiShellX(FitParamsPtr p)
1433{
1434        DOUBLE *dp;                             // Pointer to double precision wave data.
1435        float *fp;                              // Pointer to single precision wave data.
1436        DOUBLE x;
1437        DOUBLE scale,rcore,tw,ts,rhocore,rhoshel,num,bkg;               //my local names
1438        int ii;
1439        DOUBLE fval,voli,ri,sldi;
1440       
1441        if (p->waveHandle == NIL) {
1442                SetNaN64(&p->result);
1443                return NON_EXISTENT_WAVE;
1444        }
1445               
1446        x= p->x;
1447        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
1448                case NT_FP32:
1449                        fp= WaveData(p->waveHandle);
1450                        scale = fp[0];
1451                        rcore = fp[1];
1452                        ts = fp[2];
1453                        tw = fp[3];
1454                        rhocore = fp[4];
1455                        rhoshel = fp[5];
1456                        num = fp[6];
1457                        bkg = fp[7];
1458                       
1459                        //calculate with a loop, two shells at a time
1460                       
1461                        ii=0;
1462                        fval=0;
1463                       
1464                        do {
1465                            ri = rcore + (DOUBLE)ii*ts + (DOUBLE)ii*tw;
1466                            voli = 4*pi/3*ri*ri*ri;
1467                            sldi = rhocore-rhoshel;
1468                            fval += voli*sldi*F_func(ri*x);
1469                            ri += ts;
1470                            voli = 4*pi/3*ri*ri*ri;
1471                            sldi = rhoshel-rhocore;
1472                            fval += voli*sldi*F_func(ri*x);
1473                            ii+=1;              //do 2 layers at a time
1474                         } while(ii<=num-1);  //change to make 0 < num < 2 correspond to unilamellar vesicles (C. Glinka, 11/24/03)
1475                       
1476 
1477                        break;
1478                case NT_FP64:
1479                        dp= WaveData(p->waveHandle);
1480                        scale = dp[0];
1481                        rcore = dp[1];
1482                        ts = dp[2];
1483                        tw = dp[3];
1484                        rhocore = dp[4];
1485                        rhoshel = dp[5];
1486                        num = dp[6];
1487                        bkg = dp[7];
1488                       
1489                        //calculate with a loop, two shells at a time
1490                       
1491                        ii=0;
1492                        fval=0;
1493                       
1494                        do {
1495                            ri = rcore + (DOUBLE)ii*ts + (DOUBLE)ii*tw;
1496                            voli = 4*pi/3*ri*ri*ri;
1497                            sldi = rhocore-rhoshel;
1498                            fval += voli*sldi*F_func(ri*x);
1499                            ri += ts;
1500                            voli = 4*pi/3*ri*ri*ri;
1501                            sldi = rhoshel-rhocore;
1502                            fval += voli*sldi*F_func(ri*x);
1503                            ii+=1;              //do 2 layers at a time
1504                         } while(ii<=num-1);  //change to make 0 < num < 2 correspond to unilamellar vesicles (C. Glinka, 11/24/03)
1505                       
1506                        break;
1507                default:                                                                // We can't handle this wave data type.
1508                        SetNaN64(&p->result);
1509                        return REQUIRES_SP_OR_DP_WAVE;
1510        }
1511       
1512        fval *= fval;           //square it
1513        fval /= voli;           //normalize by the overall volume
1514        fval *= scale*1e8;
1515        fval += bkg;
1516       
1517        p->result= fval;
1518       
1519        return 0;
1520}
1521
1522/*
1523// calculates the scattering from a POLYDISPERSE spherical particle made up of a core (aqueous) surrounded
1524// by N spherical layers, each of which is a PAIR of shells, solvent + surfactant since there
1525//must always be a surfactant layer on the outside
1526//
1527// bragg peaks arise naturally from the periodicity of the sample
1528// resolution smeared version gives he most appropriate view of the model
1529//
1530// Polydispersity is of the total (outer) radius. This is converted into a distribution of MLV's
1531// with integer numbers of layers, with a minimum of one layer... a vesicle... depending
1532// on the parameters, the "distribution" of MLV's that is used may be truncated
1533//
1534        Warning:
1535                The call to WaveData() below returns a pointer to the middle
1536                of an unlocked Macintosh handle. In the unlikely event that your
1537                calculations could cause memory to move, you should copy the coefficient
1538                values to local variables or an array before such operations.
1539*/
1540int
1541PolyMultiShellX(FitParamsPtr p)
1542{
1543        DOUBLE *dp;                             // Pointer to double precision wave data.
1544        float *fp;                              // Pointer to single precision wave data.
1545        DOUBLE x;
1546        DOUBLE scale,rcore,tw,ts,rhocore,rhoshel,bkg;           //my local names
1547        int ii,minPairs,maxPairs,first;
1548        DOUBLE fval,ri,pi;
1549        DOUBLE avg,pd,zz,lo,hi,r1,r2,d1,d2,distr;
1550       
1551        if (p->waveHandle == NIL) {
1552                SetNaN64(&p->result);
1553                return NON_EXISTENT_WAVE;
1554        }
1555       
1556        pi = 4.0*atan(1.0);     
1557        x= p->x;
1558       
1559        switch(WaveType(p->waveHandle)){                        // We can handle single and double precision coefficient waves.
1560                case NT_FP32:
1561                        fp= WaveData(p->waveHandle);
1562                        scale = fp[0];
1563                        avg = fp[1];            // average (total) outer radius
1564                        pd = fp[2];
1565                        rcore = fp[3];
1566                        ts = fp[4];
1567                        tw = fp[5];
1568                        rhocore = fp[6];
1569                        rhoshel = fp[7];
1570                        bkg = fp[8];
1571                       
1572                        zz = (1.0/pd)*(1.0/pd)-1.0;
1573                       
1574                        //max radius set to be 5 std deviations past mean
1575                        hi = avg + pd*avg*5.0;
1576                        lo = avg - pd*avg*5.0;
1577                       
1578                        maxPairs = trunc( (hi-rcore+tw)/(ts+tw) );
1579                        minPairs = trunc( (lo-rcore+tw)/(ts+tw) );
1580                        minPairs = (minPairs < 1) ? 1 : minPairs;       // need a minimum of one
1581                       
1582                        ii=minPairs;
1583                        fval=0;
1584                        d1 = 0;
1585                        d2 = 0;
1586                        r1 = 0;
1587                        r2 = 0;
1588                        distr = 0;
1589                        first = 1;
1590                        do {
1591                            //make the current values old
1592                            r1 = r2;
1593                            d1 = d2;
1594                           
1595                            ri = (DOUBLE)ii*(ts+tw) - tw + rcore;
1596                            fval += SchulzPoint(ri,avg,zz) * MultiShellGuts(x, rcore, ts, tw, rhocore, rhoshel, ii) * (4*pi/3*ri*ri*ri);
1597                            // get a running integration of the fraction of the distribution used, but not the first time
1598                            r2 = ri;
1599                            d2 = SchulzPoint(ri,avg,zz);
1600                            if( !first ) {
1601                                distr += 0.5*(d1+d2)*(r2-r1);           //cheap trapezoidal integration
1602                            }
1603                            ii+=1;
1604                            first = 0;
1605                         } while(ii<=maxPairs);
1606                       
1607 
1608                        break;
1609                case NT_FP64:
1610                        dp= WaveData(p->waveHandle);
1611                        scale = dp[0];
1612                        avg = dp[1];            // average (total) outer radius
1613                        pd = dp[2];
1614                        rcore = dp[3];
1615                        ts = dp[4];
1616                        tw = dp[5];
1617                        rhocore = dp[6];
1618                        rhoshel = dp[7];
1619                        bkg = dp[8];
1620                       
1621                        zz = (1.0/pd)*(1.0/pd)-1.0;
1622                       
1623                        //max radius set to be 5 std deviations past mean
1624                        hi = avg + pd*avg*5.0;
1625                        lo = avg - pd*avg*5.0;
1626                       
1627                        maxPairs = trunc( (hi-rcore+tw)/(ts+tw) );
1628                        minPairs = trunc( (lo-rcore+tw)/(ts+tw) );
1629                        minPairs = (minPairs < 1) ? 1 : minPairs;       // need a minimum of one
1630                       
1631                        ii=minPairs;
1632                        fval=0;
1633                        d1 = 0;
1634                        d2 = 0;
1635                        r1 = 0;
1636                        r2 = 0;
1637                        distr = 0;
1638                        first = 1;
1639                        do {
1640                            //make the current values old
1641                            r1 = r2;
1642                            d1 = d2;
1643                           
1644                            ri = (DOUBLE)ii*(ts+tw) - tw + rcore;
1645                            fval += SchulzPoint(ri,avg,zz) * MultiShellGuts(x, rcore, ts, tw, rhocore, rhoshel, ii) * (4*pi/3*ri*ri*ri);
1646                            // get a running integration of the fraction of the distribution used, but not the first time
1647                            r2 = ri;
1648                            d2 = SchulzPoint(ri,avg,zz);
1649                            if( !first ) {
1650                                distr += 0.5*(d1+d2)*(r2-r1);           //cheap trapezoidal integration
1651                            }
1652                            ii+=1;
1653                            first = 0;
1654                         } while(ii<=maxPairs);
1655                       
1656                        break;
1657                default:                                                                // We can't handle this wave data type.
1658                        SetNaN64(&p->result);
1659                        return REQUIRES_SP_OR_DP_WAVE;
1660        }
1661       
1662        fval /= 4*pi/3*avg*avg*avg;             //normalize by the overall volume
1663        fval /= distr;
1664        fval *= scale;
1665        fval += bkg;
1666       
1667        p->result= fval;
1668       
1669        return 0;
1670}
1671
1672DOUBLE
1673MultiShellGuts(DOUBLE x,DOUBLE rcore,DOUBLE ts,DOUBLE tw,DOUBLE rhocore,DOUBLE rhoshel,int num) {
1674
1675    DOUBLE ri,sldi,fval,voli;
1676    int ii;
1677   
1678    ii=0;
1679    fval=0;
1680   
1681    do {
1682        ri = rcore + (DOUBLE)ii*ts + (DOUBLE)ii*tw;
1683        voli = 4*pi/3*ri*ri*ri;
1684        sldi = rhocore-rhoshel;
1685        fval += voli*sldi*F_func(ri*x);
1686        ri += ts;
1687        voli = 4*pi/3*ri*ri*ri;
1688        sldi = rhoshel-rhocore;
1689        fval += voli*sldi*F_func(ri*x);
1690        ii+=1;          //do 2 layers at a time
1691    } while(ii<=num-1);  //change to make 0 < num < 2 correspond to unilamellar vesicles (C. Glinka, 11/24/03)
1692   
1693    fval *= fval;
1694    fval /= voli;
1695    fval *= 1e8;
1696   
1697    return(fval);       // this result still needs to be multiplied by scale and have background added
1698}
1699
1700static DOUBLE
1701SchulzPoint(DOUBLE x, DOUBLE avg, DOUBLE zz) {
1702
1703    DOUBLE dr;
1704   
1705    dr = zz*log(x) - gammln(zz+1)+(zz+1)*log((zz+1)/avg)-(x/avg*(zz+1));
1706    return (exp(dr));
1707}
1708
1709static DOUBLE
1710gammln(double xx) {
1711
1712    double x,y,tmp,ser;
1713    static double cof[6]={76.18009172947146,-86.50532032941677,
1714            24.01409824083091,-1.231739572450155,
1715            0.1208650973866179e-2,-0.5395239384953e-5};
1716    int j;
1717
1718    y=x=xx;
1719    tmp=x+5.5;
1720    tmp -= (x+0.5)*log(tmp);
1721    ser=1.000000000190015;
1722    for (j=0;j<=5;j++) ser += cof[j]/++y;
1723    return -tmp+log(2.5066282746310005*ser/x);
1724}
1725
1726DOUBLE
1727F_func(double qr) {
1728        return(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr));
1729}
1730
1731
1732
1733#pragma XOP_RESET_STRUCT_PACKING                        // All structures are 2-byte-aligned.
1734
1735///////////end of XOP
1736
1737
Note: See TracBrowser for help on using the repository browser.