source: sans/XOP_Dev/SANSAnalysis/lib/libCylinder.c @ 632

Last change on this file since 632 was 632, checked in by srkline, 12 years ago

incorporating Jae-Hie's changes to the library functions to explicitly handle the cases at qr==0. Also a lot of changes to explicitly declare constant values as floating point, rather than chance the interpretation as integer. (probably not necessary)

  • Property svn:executable set to *
File size: 79.9 KB
RevLine 
[97]1/*      CylinderFit.c
2
3A simplified project designed to act as a template for your curve fitting function.
4The fitting function is a Cylinder form factor. No resolution effects are included (yet)
5*/
6
7#include "StandardHeaders.h"                    // Include ANSI headers, Mac headers
8#include "GaussWeights.h"
9#include "libCylinder.h"
10
11/*      CylinderForm  :  calculates the form factor of a cylinder at the give x-value p->x
12
13Warning:
14The call to WaveData() below returns a pointer to the middle
15of an unlocked Macintosh handle. In the unlikely event that your
16calculations could cause memory to move, you should copy the coefficient
17values to local variables or an array before such operations.
18*/
19double
20CylinderForm(double dp[], double q)
21{
22        int i;
23        double Pi;
[235]24        double scale,radius,length,delrho,bkg,halfheight,sldCyl,sldSolv;                //local variables of coefficient wave
[97]25        int nord=76;                    //order of integration
26        double uplim,lolim;             //upper and lower integration limits
27        double summ,zi,yyy,answer,vcyl;                 //running tally of integration
28       
29        Pi = 4.0*atan(1.0);
[632]30        lolim = 0.0;
[97]31        uplim = Pi/2.0;
32       
33        summ = 0.0;                     //initialize intergral
34       
35        scale = dp[0];                  //make local copies in case memory moves
36        radius = dp[1];
37        length = dp[2];
[235]38        sldCyl = dp[3];
39        sldSolv = dp[4];
40        bkg = dp[5];
41       
42        delrho = sldCyl-sldSolv;
[97]43        halfheight = length/2.0;
44        for(i=0;i<nord;i++) {
45                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
46                yyy = Gauss76Wt[i] * CylKernel(q, radius, halfheight, zi);
47                summ += yyy;
48        }
49       
50        answer = (uplim-lolim)/2.0*summ;
51        // Multiply by contrast^2
52        answer *= delrho*delrho;
53        //normalize by cylinder volume
54        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
55        vcyl=Pi*radius*radius*length;
56        answer *= vcyl;
57        //convert to [cm-1]
58        answer *= 1.0e8;
59        //Scale
60        answer *= scale;
61        // add in the background
62        answer += bkg;
63       
64        return answer;
65}
66
67/*      EllipCyl76X  :  calculates the form factor of a elliptical cylinder at the given x-value p->x
68
69Uses 76 pt Gaussian quadrature for both integrals
70
71Warning:
72The call to WaveData() below returns a pointer to the middle
73of an unlocked Macintosh handle. In the unlikely event that your
74calculations could cause memory to move, you should copy the coefficient
75values to local variables or an array before such operations.
76*/
77double
78EllipCyl76(double dp[], double q)
79{
80        int i,j;
[235]81        double Pi,slde,sld;
[97]82        double scale,ra,nu,length,delrho,bkg;           //local variables of coefficient wave
83        int nord=76;                    //order of integration
84        double va,vb;           //upper and lower integration limits
85        double summ,zi,yyy,answer,vell;                 //running tally of integration
[632]86        double summj,vaj,vbj,zij,arg, si;                       //for the inner integration
87
[97]88        Pi = 4.0*atan(1.0);
[632]89        va = 0.0;
90        vb = 1.0;               //orintational average, outer integral
91        vaj=0.0;
[97]92        vbj=Pi;         //endpoints of inner integral
93       
94        summ = 0.0;                     //initialize intergral
95       
96        scale = dp[0];                  //make local copies in case memory moves
97        ra = dp[1];
98        nu = dp[2];
99        length = dp[3];
[235]100        slde = dp[4];
101        sld = dp[5];
102        delrho = slde - sld;
103        bkg = dp[6];
104       
[97]105        for(i=0;i<nord;i++) {
106                //setup inner integral over the ellipsoidal cross-section
107                summj=0;
108                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy
[632]109                arg = ra*sqrt(1.0-zi*zi);
[97]110                for(j=0;j<nord;j++) {
111                        //76 gauss points for the inner integral as well
112                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "y" dummy
113                        yyy = Gauss76Wt[j] * EllipCylKernel(q,arg,nu,zij);
114                        summj += yyy;
115                }
116                //now calculate the value of the inner integral
117                answer = (vbj-vaj)/2.0*summj;
118                //divide integral by Pi
119                answer /=Pi;
120               
121                //now calculate outer integral
[632]122                arg = q*length*zi/2.0;
123                if (arg == 0.0){
124                        si = 1.0;
125                }else{
126                        si = sin(arg) * sin(arg) / arg / arg;
127                }
128                yyy = Gauss76Wt[i] * answer * si;
[97]129                summ += yyy;
130        }
131        answer = (vb-va)/2.0*summ;
132        // Multiply by contrast^2
133        answer *= delrho*delrho;
134        //normalize by cylinder volume
135        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
136        vell = Pi*ra*(nu*ra)*length;
137        answer *= vell;
138        //convert to [cm-1]
139        answer *= 1.0e8;
140        //Scale
141        answer *= scale;
142        // add in the background
143        answer += bkg;
144       
145        return answer;
146}
147
148/*      EllipCyl20X  :  calculates the form factor of a elliptical cylinder at the given x-value p->x
149
150Uses 76 pt Gaussian quadrature for orientational integral
151Uses 20 pt quadrature for the inner integral over the elliptical cross-section
152
153Warning:
154The call to WaveData() below returns a pointer to the middle
155of an unlocked Macintosh handle. In the unlikely event that your
156calculations could cause memory to move, you should copy the coefficient
157values to local variables or an array before such operations.
158*/
159double
160EllipCyl20(double dp[], double q)
161{
162        int i,j;
[235]163        double Pi,slde,sld;
[97]164        double scale,ra,nu,length,delrho,bkg;           //local variables of coefficient wave
165        int nordi=76;                   //order of integration
166        int nordj=20;
167        double va,vb;           //upper and lower integration limits
168        double summ,zi,yyy,answer,vell;                 //running tally of integration
[632]169        double summj,vaj,vbj,zij,arg,si;                        //for the inner integration
170
[97]171        Pi = 4.0*atan(1.0);
[632]172        va = 0.0;
173        vb = 1.0;               //orintational average, outer integral
174        vaj=0.0;
[97]175        vbj=Pi;         //endpoints of inner integral
176       
177        summ = 0.0;                     //initialize intergral
178       
179        scale = dp[0];                  //make local copies in case memory moves
180        ra = dp[1];
181        nu = dp[2];
182        length = dp[3];
[235]183        slde = dp[4];
184        sld = dp[5];
185        delrho = slde - sld;
186        bkg = dp[6];
187       
[97]188        for(i=0;i<nordi;i++) {
189                //setup inner integral over the ellipsoidal cross-section
190                summj=0;
191                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy
[632]192                arg = ra*sqrt(1.0-zi*zi);
[97]193                for(j=0;j<nordj;j++) {
194                        //20 gauss points for the inner integral
195                        zij = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "y" dummy
196                        yyy = Gauss20Wt[j] * EllipCylKernel(q,arg,nu,zij);
197                        summj += yyy;
198                }
199                //now calculate the value of the inner integral
200                answer = (vbj-vaj)/2.0*summj;
201                //divide integral by Pi
202                answer /=Pi;
203               
204                //now calculate outer integral
[632]205                arg = q*length*zi/2.0;
206                if (arg == 0.0){
207                        si = 1.0;
208                }else{
209                        si = sin(arg) * sin(arg) / arg / arg;
210                }
211                yyy = Gauss76Wt[i] * answer * si;
[97]212                summ += yyy;
213        }
214       
215        answer = (vb-va)/2.0*summ;
216        // Multiply by contrast^2
217        answer *= delrho*delrho;
218        //normalize by cylinder volume
219        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
220        vell = Pi*ra*(nu*ra)*length;
221        answer *= vell;
222        //convert to [cm-1]
223        answer *= 1.0e8;
224        //Scale
225        answer *= scale;
226        // add in the background
[632]227        answer += bkg;
228
[97]229        return answer;
230}
231
232/*      TriaxialEllipsoidX  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x
233
234Uses 76 pt Gaussian quadrature for both integrals
235
236Warning:
237The call to WaveData() below returns a pointer to the middle
238of an unlocked Macintosh handle. In the unlikely event that your
239calculations could cause memory to move, you should copy the coefficient
240values to local variables or an array before such operations.
241*/
242double
243TriaxialEllipsoid(double dp[], double q)
244{
245        int i,j;
246        double Pi;
247        double scale,aa,bb,cc,delrho,bkg;               //local variables of coefficient wave
248        int nordi=76;                   //order of integration
249        int nordj=76;
250        double va,vb;           //upper and lower integration limits
251        double summ,zi,yyy,answer;                      //running tally of integration
[235]252        double summj,vaj,vbj,zij,slde,sld;                      //for the inner integration
[97]253       
254        Pi = 4.0*atan(1.0);
[453]255        va = 0.0;
256        vb = 1.0;               //orintational average, outer integral
257        vaj = 0.0;
258        vbj = 1.0;              //endpoints of inner integral
[97]259       
260        summ = 0.0;                     //initialize intergral
261       
262        scale = dp[0];                  //make local copies in case memory moves
263        aa = dp[1];
264        bb = dp[2];
265        cc = dp[3];
[235]266        slde = dp[4];
267        sld = dp[5];
268        delrho = slde - sld;
269        bkg = dp[6];
[97]270        for(i=0;i<nordi;i++) {
271                //setup inner integral over the ellipsoidal cross-section
[453]272                summj=0.0;
[97]273                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy
274                for(j=0;j<nordj;j++) {
275                        //20 gauss points for the inner integral
276                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "y" dummy
277                        yyy = Gauss76Wt[j] * TriaxialKernel(q,aa,bb,cc,zi,zij);
278                        summj += yyy;
279                }
280                //now calculate the value of the inner integral
281                answer = (vbj-vaj)/2.0*summj;
282               
283                //now calculate outer integral
284                yyy = Gauss76Wt[i] * answer;
285                summ += yyy;
286        }               //final scaling is done at the end of the function, after the NT_FP64 case
287       
288        answer = (vb-va)/2.0*summ;
289        // Multiply by contrast^2
290        answer *= delrho*delrho;
291        //normalize by ellipsoid volume
[453]292        answer *= 4.0*Pi/3.0*aa*bb*cc;
[97]293        //convert to [cm-1]
294        answer *= 1.0e8;
295        //Scale
296        answer *= scale;
297        // add in the background
298        answer += bkg;
299       
300        return answer;
301}
302
303/*      ParallelepipedX  :  calculates the form factor of a Parallelepiped (a rectangular solid)
304at the given x-value p->x
305
306Uses 76 pt Gaussian quadrature for both integrals
307
308Warning:
309The call to WaveData() below returns a pointer to the middle
310of an unlocked Macintosh handle. In the unlikely event that your
311calculations could cause memory to move, you should copy the coefficient
312values to local variables or an array before such operations.
313*/
314double
315Parallelepiped(double dp[], double q)
316{
317        int i,j;
318        double scale,aa,bb,cc,delrho,bkg;               //local variables of coefficient wave
319        int nordi=76;                   //order of integration
320        int nordj=76;
321        double va,vb;           //upper and lower integration limits
322        double summ,yyy,answer;                 //running tally of integration
323        double summj,vaj,vbj;                   //for the inner integration
[235]324        double mu,mudum,arg,sigma,uu,vol,sldp,sld;
[97]325       
326       
327        //      Pi = 4.0*atan(1.0);
[632]328        va = 0.0;
329        vb = 1.0;               //orintational average, outer integral
330        vaj = 0.0;
331        vbj = 1.0;              //endpoints of inner integral
332
[97]333        summ = 0.0;                     //initialize intergral
334       
335        scale = dp[0];                  //make local copies in case memory moves
336        aa = dp[1];
337        bb = dp[2];
338        cc = dp[3];
[235]339        sldp = dp[4];
340        sld = dp[5];
341        delrho = sldp - sld;
342        bkg = dp[6];
[97]343       
344        mu = q*bb;
345        vol = aa*bb*cc;
346        // normalize all WRT bb
347        aa = aa/bb;
348        cc = cc/bb;
349       
350        for(i=0;i<nordi;i++) {
351                //setup inner integral over the ellipsoidal cross-section
[632]352                summj=0.0;
[97]353                sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;          //the outer dummy
354               
355                for(j=0;j<nordj;j++) {
356                        //76 gauss points for the inner integral
357                        uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;         //the inner dummy
[632]358                        mudum = mu*sqrt(1.0-sigma*sigma);
[97]359                        yyy = Gauss76Wt[j] * PPKernel(aa,mudum,uu);
360                        summj += yyy;
361                }
362                //now calculate the value of the inner integral
363                answer = (vbj-vaj)/2.0*summj;
[632]364
365                arg = mu*cc*sigma/2.0;
366                if ( arg == 0.0 ) {
367                        answer *= 1.0;
[97]368                } else {
369                        answer *= sin(arg)*sin(arg)/arg/arg;
370                }
371               
372                //now sum up the outer integral
373                yyy = Gauss76Wt[i] * answer;
374                summ += yyy;
375        }               //final scaling is done at the end of the function, after the NT_FP64 case
376       
377        answer = (vb-va)/2.0*summ;
378        // Multiply by contrast^2
379        answer *= delrho*delrho;
380        //normalize by volume
381        answer *= vol;
382        //convert to [cm-1]
383        answer *= 1.0e8;
384        //Scale
385        answer *= scale;
386        // add in the background
387        answer += bkg;
388       
389        return answer;
390}
391
392/*      HollowCylinderX  :  calculates the form factor of a Hollow Cylinder
393at the given x-value p->x
394
395Uses 76 pt Gaussian quadrature for the single integral
396
397Warning:
398The call to WaveData() below returns a pointer to the middle
399of an unlocked Macintosh handle. In the unlikely event that your
400calculations could cause memory to move, you should copy the coefficient
401values to local variables or an array before such operations.
402*/
403double
404HollowCylinder(double dp[], double q)
405{
406        int i;
407        double scale,rcore,rshell,length,delrho,bkg;            //local variables of coefficient wave
408        int nord=76;                    //order of integration
409        double va,vb,zi;                //upper and lower integration limits
[235]410        double summ,answer,pi,sldc,sld;                 //running tally of integration
[97]411       
412        pi = 4.0*atan(1.0);
[632]413        va = 0.0;
414        vb = 1.0;               //limits of numerical integral
415
[97]416        summ = 0.0;                     //initialize intergral
417       
418        scale = dp[0];          //make local copies in case memory moves
419        rcore = dp[1];
420        rshell = dp[2];
421        length = dp[3];
[235]422        sldc = dp[4];
423        sld = dp[5];
424        delrho = sldc - sld;
425        bkg = dp[6];
[97]426       
427        for(i=0;i<nord;i++) {
428                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
429                summ += Gauss76Wt[i] * HolCylKernel(q, rcore, rshell, length, zi);
430        }
431       
432        answer = (vb-va)/2.0*summ;
433        // Multiply by contrast^2
434        answer *= delrho*delrho;
435        //normalize by volume
436        answer *= pi*(rshell*rshell-rcore*rcore)*length;
437        //convert to [cm-1]
438        answer *= 1.0e8;
439        //Scale
440        answer *= scale;
441        // add in the background
442        answer += bkg;
443       
444        return answer;
445}
446
447/*      EllipsoidFormX  :  calculates the form factor of an ellipsoid of revolution with semiaxes a:a:nua
448at the given x-value p->x
449
450Uses 76 pt Gaussian quadrature for the single integral
451
452Warning:
453The call to WaveData() below returns a pointer to the middle
454of an unlocked Macintosh handle. In the unlikely event that your
455calculations could cause memory to move, you should copy the coefficient
456values to local variables or an array before such operations.
457*/
458double
459EllipsoidForm(double dp[], double q)
460{
461        int i;
462        double scale,a,nua,delrho,bkg;          //local variables of coefficient wave
463        int nord=76;                    //order of integration
464        double va,vb,zi;                //upper and lower integration limits
[235]465        double summ,answer,pi,slde,sld;                 //running tally of integration
[97]466       
467        pi = 4.0*atan(1.0);
[632]468        va = 0.0;
469        vb = 1.0;               //limits of numerical integral
470
[97]471        summ = 0.0;                     //initialize intergral
472       
473        scale = dp[0];                  //make local copies in case memory moves
474        nua = dp[1];
475        a = dp[2];
[235]476        slde = dp[3];
477        sld = dp[4];
478        delrho = slde - sld;
479        bkg = dp[5];
[97]480       
481        for(i=0;i<nord;i++) {
482                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
483                summ += Gauss76Wt[i] * EllipsoidKernel(q, a, nua, zi);
484        }
485       
486        answer = (vb-va)/2.0*summ;
487        // Multiply by contrast^2
488        answer *= delrho*delrho;
489        //normalize by volume
[632]490        answer *= 4.0*pi/3.0*a*a*nua;
[97]491        //convert to [cm-1]
492        answer *= 1.0e8;
493        //Scale
494        answer *= scale;
495        // add in the background
496        answer += bkg;
497       
498        return answer;
499}
500
501
502/*      Cyl_PolyRadiusX  :  calculates the form factor of a cylinder at the given x-value p->x
503the cylinder has a polydisperse cross section
504
505*/
506double
507Cyl_PolyRadius(double dp[], double q)
508{
509        int i;
510        double scale,radius,length,pd,delrho,bkg;               //local variables of coefficient wave
511        int nord=20;                    //order of integration
512        double uplim,lolim;             //upper and lower integration limits
513        double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
[235]514        double range,zz,Pi,sldc,sld;
[97]515       
516        Pi = 4.0*atan(1.0);
517        range = 3.4;
518       
519        summ = 0.0;                     //initialize intergral
520       
521        scale = dp[0];                  //make local copies in case memory moves
522        radius = dp[1];
523        length = dp[2];
524        pd = dp[3];
[235]525        sldc = dp[4];
526        sld = dp[5];
527        delrho = sldc - sld;
528        bkg = dp[6];
[97]529       
530        zz = (1.0/pd)*(1.0/pd) - 1.0;
531       
532        lolim = radius*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
[632]533        if(lolim<0.0) {
534                lolim = 0.0;
[97]535        }
536        if(pd>0.3) {
537                range = 3.4 + (pd-0.3)*18.0;
538        }
539        uplim = radius*(1.0+range*pd);
540       
541        for(i=0;i<nord;i++) {
542                zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
543                yyy = Gauss20Wt[i] * Cyl_PolyRadKernel(q, radius, length, zz, delrho, zi);
544                summ += yyy;
545        }
546       
547        answer = (uplim-lolim)/2.0*summ;
548        //normalize by average cylinder volume
549        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
550        Vpoly=Pi*radius*radius*length*(zz+2.0)/(zz+1.0);
551        answer /= Vpoly;
552        //convert to [cm-1]
553        answer *= 1.0e8;
554        //Scale
555        answer *= scale;
556        // add in the background
557        answer += bkg;
558       
559        return answer;
560}
561
562/*      Cyl_PolyLengthX  :  calculates the form factor of a cylinder at the given x-value p->x
563the cylinder has a polydisperse Length
564
565*/
566double
567Cyl_PolyLength(double dp[], double q)
568{
569        int i;
570        double scale,radius,length,pd,delrho,bkg;               //local variables of coefficient wave
571        int nord=20;                    //order of integration
572        double uplim,lolim;             //upper and lower integration limits
573        double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
[235]574        double range,zz,Pi,sldc,sld;
[97]575       
576       
577        Pi = 4.0*atan(1.0);
578        range = 3.4;
579       
580        summ = 0.0;                     //initialize intergral
581       
582        scale = dp[0];                  //make local copies in case memory moves
583        radius = dp[1];
584        length = dp[2];
585        pd = dp[3];
[235]586        sldc = dp[4];
587        sld = dp[5];
588        delrho = sldc - sld;
589        bkg = dp[6];
[97]590       
591        zz = (1.0/pd)*(1.0/pd) - 1.0;
592       
593        lolim = length*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
[632]594        if(lolim<0.0) {
595                lolim = 0.0;
[97]596        }
597        if(pd>0.3) {
598                range = 3.4 + (pd-0.3)*18.0;
599        }
600        uplim = length*(1.0+range*pd);
601       
602        for(i=0;i<nord;i++) {
603                zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
604                yyy = Gauss20Wt[i] * Cyl_PolyLenKernel(q, radius, length, zz, delrho, zi);
605                summ += yyy;
606        }
607       
608        answer = (uplim-lolim)/2.0*summ;
609        //normalize by average cylinder volume (first moment)
610        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
611        Vpoly=Pi*radius*radius*length;
612        answer /= Vpoly;
613        //convert to [cm-1]
614        answer *= 1.0e8;
615        //Scale
616        answer *= scale;
617        // add in the background
618        answer += bkg;
619       
620        return answer;
621}
622
623/*      CoreShellCylinderX  :  calculates the form factor of a cylinder at the given x-value p->x
624the cylinder has a core-shell structure
625
626*/
627double
628CoreShellCylinder(double dp[], double q)
629{
630        int i;
631        double scale,rcore,length,bkg;          //local variables of coefficient wave
632        double thick,rhoc,rhos,rhosolv;
633        int nord=76;                    //order of integration
634        double uplim,lolim,halfheight;          //upper and lower integration limits
635        double summ,zi,yyy,answer,Vcyl;                 //running tally of integration
636        double Pi;
637       
638        Pi = 4.0*atan(1.0);
639       
640        lolim = 0.0;
641        uplim = Pi/2.0;
642       
643        summ = 0.0;                     //initialize intergral
644       
645        scale = dp[0];          //make local copies in case memory moves
646        rcore = dp[1];
647        thick = dp[2];
648        length = dp[3];
649        rhoc = dp[4];
650        rhos = dp[5];
651        rhosolv = dp[6];
652        bkg = dp[7];
653       
654        halfheight = length/2.0;
655       
656        for(i=0;i<nord;i++) {
657                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
658                yyy = Gauss76Wt[i] * CoreShellCylKernel(q, rcore, thick, rhoc,rhos,rhosolv, halfheight, zi);
659                summ += yyy;
660        }
661       
662        answer = (uplim-lolim)/2.0*summ;
663        // length is the total core length
664        Vcyl=Pi*(rcore+thick)*(rcore+thick)*(length+2.0*thick);
665        answer /= Vcyl;
666        //convert to [cm-1]
667        answer *= 1.0e8;
668        //Scale
669        answer *= scale;
670        // add in the background
671        answer += bkg;
672       
673        return answer;
674}
675
676
677/*      PolyCoShCylinderX  :  calculates the form factor of a core-shell cylinder at the given x-value p->x
678the cylinder has a polydisperse CORE radius
679
680*/
681double
682PolyCoShCylinder(double dp[], double q)
683{
684        int i;
685        double scale,radius,length,sigma,bkg;           //local variables of coefficient wave
686        double rad,radthick,facthick,rhoc,rhos,rhosolv;
687        int nord=20;                    //order of integration
688        double uplim,lolim;             //upper and lower integration limits
689        double summ,yyy,answer,Vpoly;                   //running tally of integration
690        double Pi,AR,Rsqrsumm,Rsqryyy,Rsqr;
[501]691               
[97]692        Pi = 4.0*atan(1.0);
693       
694        summ = 0.0;                     //initialize intergral
695        Rsqrsumm = 0.0;
696       
697        scale = dp[0];
698        radius = dp[1];
699        sigma = dp[2];                          //sigma is the standard mean deviation
700        length = dp[3];
701        radthick = dp[4];
702        facthick= dp[5];
703        rhoc = dp[6];
704        rhos = dp[7];
705        rhosolv = dp[8];
706        bkg = dp[9];
707       
708        lolim = exp(log(radius)-(4.*sigma));
[632]709        if (lolim<0.0) {
710                lolim=0.0;              //to avoid numerical error when  va<0 (-ve r value)
[97]711        }
712        uplim = exp(log(radius)+(4.*sigma));
713       
714        for(i=0;i<nord;i++) {
715                rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
716                AR=(1.0/(rad*sigma*sqrt(2.0*Pi)))*exp(-(0.5*((log(radius/rad))/sigma)*((log(radius/rad))/sigma)));
717                yyy = AR* Gauss20Wt[i] * CSCylIntegration(q,rad,radthick,facthick,rhoc,rhos,rhosolv,length);
718                Rsqryyy= Gauss20Wt[i] * AR * (rad+radthick)*(rad+radthick);             //SRK normalize to total dimensions
719                summ += yyy;
720                Rsqrsumm += Rsqryyy;
721        }
722       
723        answer = (uplim-lolim)/2.0*summ;
724        Rsqr = (uplim-lolim)/2.0*Rsqrsumm;
725        //normalize by average cylinder volume
726        Vpoly = Pi*Rsqr*(length+2*facthick);
727        answer /= Vpoly;
728        //convert to [cm-1]
729        answer *= 1.0e8;
730        //Scale
731        answer *= scale;
732        // add in the background
733        answer += bkg;
734       
735        return answer;
736}
737
738/*      OblateFormX  :  calculates the form factor of a core-shell Oblate ellipsoid at the given x-value p->x
739the ellipsoid has a core-shell structure
740
741*/
742double
743OblateForm(double dp[], double q)
744{
745        int i;
746        double scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg;
747        int nord=76;                    //order of integration
748        double uplim,lolim;             //upper and lower integration limits
749        double summ,zi,yyy,answer,oblatevol;                    //running tally of integration
[235]750        double Pi,sldc,slds,sld;
[97]751       
752        Pi = 4.0*atan(1.0);
753       
754        lolim = 0.0;
755        uplim = 1.0;
756       
757        summ = 0.0;                     //initialize intergral
758       
759       
760        scale = dp[0];          //make local copies in case memory moves
761        crmaj = dp[1];
762        crmin = dp[2];
763        trmaj = dp[3];
764        trmin = dp[4];
[235]765        sldc = dp[5];
766        slds = dp[6];
767        sld = dp[7];
768        delpc = sldc - slds;    //core - shell
769        delps = slds - sld;             //shell - solvent
[632]770        bkg = dp[8];
771
[97]772        for(i=0;i<nord;i++) {
773                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
774                yyy = Gauss76Wt[i] * gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,q);
775                summ += yyy;
776        }
777       
778        answer = (uplim-lolim)/2.0*summ;
779        // normalize by particle volume
[632]780        oblatevol = 4.0*Pi/3.0*trmaj*trmaj*trmin;
[97]781        answer /= oblatevol;
782       
783        //convert to [cm-1]
784        answer *= 1.0e8;
785        //Scale
786        answer *= scale;
787        // add in the background
788        answer += bkg;
789       
790        return answer;
791}
792
793/*      ProlateFormX  :  calculates the form factor of a core-shell Prolate ellipsoid at the given x-value p->x
794the ellipsoid has a core-shell structure
795
796*/
797double
798ProlateForm(double dp[], double q)
799{
800        int i;
801        double scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg;
802        int nord=76;                    //order of integration
803        double uplim,lolim;             //upper and lower integration limits
804        double summ,zi,yyy,answer,prolatevol;                   //running tally of integration
[235]805        double Pi,sldc,slds,sld;
[97]806       
807        Pi = 4.0*atan(1.0);
808       
809        lolim = 0.0;
810        uplim = 1.0;
811       
812        summ = 0.0;                     //initialize intergral
813       
814        scale = dp[0];          //make local copies in case memory moves
815        crmaj = dp[1];
816        crmin = dp[2];
817        trmaj = dp[3];
818        trmin = dp[4];
[235]819        sldc = dp[5];
820        slds = dp[6];
821        sld = dp[7];
822        delpc = sldc - slds;            //core - shell
823        delps = slds - sld;             //shell  - sovent
[632]824        bkg = dp[8];
825
[97]826        for(i=0;i<nord;i++) {
827                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
828                yyy = Gauss76Wt[i] * gfn2(zi,crmaj,crmin,trmaj,trmin,delpc,delps,q);
829                summ += yyy;
830        }
831       
832        answer = (uplim-lolim)/2.0*summ;
833        // normalize by particle volume
[632]834        prolatevol = 4.0*Pi/3.0*trmaj*trmin*trmin;
[97]835        answer /= prolatevol;
836       
837        //convert to [cm-1]
838        answer *= 1.0e8;
839        //Scale
840        answer *= scale;
841        // add in the background
842        answer += bkg;
843       
844        return answer;
845}
846
847
848/*      StackedDiscsX  :  calculates the form factor of a stacked "tactoid" of core shell disks
849like clay platelets that are not exfoliated
850
851*/
852double
853StackedDiscs(double dp[], double q)
854{
855        int i;
856        double scale,length,bkg,rcore,thick,rhoc,rhol,rhosolv,N,gsd;            //local variables of coefficient wave
857        double va,vb,vcyl,summ,yyy,zi,halfheight,d,answer;
858        int nord=76;                    //order of integration
859        double Pi;
860       
861       
862        Pi = 4.0*atan(1.0);
863       
864        va = 0.0;
865        vb = Pi/2.0;
866       
867        summ = 0.0;                     //initialize intergral
868       
869        scale = dp[0];
870        rcore = dp[1];
871        length = dp[2];
872        thick = dp[3];
873        rhoc = dp[4];
874        rhol = dp[5];
875        rhosolv = dp[6];
876        N = dp[7];
877        gsd = dp[8];
878        bkg = dp[9];
879       
880        d=2.0*thick+length;
881        halfheight = length/2.0;
882       
883        for(i=0;i<nord;i++) {
884                zi = ( Gauss76Z[i]*(vb-va) + vb + va )/2.0;
885                yyy = Gauss76Wt[i] * Stackdisc_kern(q, rcore, rhoc,rhol,rhosolv, halfheight,thick,zi,gsd,d,N);
886                summ += yyy;
887        }
888       
889        answer = (vb-va)/2.0*summ;
890        // length is the total core length
891        vcyl=Pi*rcore*rcore*(2.0*thick+length)*N;
892        answer /= vcyl;
893        //Convert to [cm-1]
894        answer *= 1.0e8;
895        //Scale
896        answer *= scale;
897        // add in the background
898        answer += bkg;
899       
900        return answer;
901}
902
903
904/*      LamellarFFX  :  calculates the form factor of a lamellar structure - no S(q) effects included
905
906*/
907double
908LamellarFF(double dp[], double q)
909{
910        double scale,del,sig,contr,bkg;         //local variables of coefficient wave
911        double inten, qval,Pq;
[235]912        double Pi,sldb,sld;
[97]913       
914       
915        Pi = 4.0*atan(1.0);
916        scale = dp[0];
917        del = dp[1];
918        sig = dp[2]*del;
[235]919        sldb = dp[3];
920        sld = dp[4];
921        contr = sldb - sld;
922        bkg = dp[5];
[156]923        qval=q;
[97]924       
925        Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig));
926       
927        inten = 2.0*Pi*scale*Pq/(qval*qval);            //this is now dimensionless...
928       
929        inten /= del;                   //normalize by the thickness (in A)
930       
931        inten *= 1.0e8;         // 1/A to 1/cm
932       
933        return(inten+bkg);
934}
935
936/*      LamellarPSX  :  calculates the form factor of a lamellar structure - with S(q) effects included
[356]937--- now the proper resolution effects are used - the "default" resolution is turned off (= 0) and the
938model is smeared just like any other function
939*/
[97]940double
941LamellarPS(double dp[], double q)
942{
943        double scale,dd,del,sig,contr,NN,Cp,bkg;                //local variables of coefficient wave
[235]944        double inten,qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ;
945        double Pi,Euler,dQDefault,fii,sldb,sld;
[97]946        int ii,NNint;
[235]947//      char buf[256];
948
[97]949       
950        Euler = 0.5772156649;           // Euler's constant
[356]951//      dQDefault = 0.0025;             //[=] 1/A, q-resolution, default value
952        dQDefault = 0.0;
[97]953        dQ = dQDefault;
954       
955        Pi = 4.0*atan(1.0);
956        qval = q;
957       
958        scale = dp[0];
959        dd = dp[1];
960        del = dp[2];
961        sig = dp[3]*del;
[235]962        sldb = dp[4];
963        sld = dp[5];
964        contr = sldb - sld;
965        NN = trunc(dp[6]);              //be sure that NN is an integer
966        Cp = dp[7];
967        bkg = dp[8];
[97]968       
969        Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig));
970       
971        NNint = (int)NN;                //cast to an integer for the loop
[235]972       
973//      sprintf(buf, "qval = %g\r", qval);
974//      XOPNotice(buf);
975       
[97]976        ii=0;
977        Sq = 0.0;
978        for(ii=1;ii<(NNint-1);ii+=1) {
979       
980                fii = (double)ii;               //do I really need to do this?
981               
982                temp = 0.0;
[235]983                alpha = Cp/4.0/Pi/Pi*(log(Pi*fii) + Euler);
[97]984                t1 = 2.0*dQ*dQ*dd*dd*alpha;
985                t2 = 2.0*qval*qval*dd*dd*alpha;
[235]986                t3 = dQ*dQ*dd*dd*fii*fii;
[97]987               
[235]988                temp = 1.0-fii/NN;
989                temp *= cos(dd*qval*fii/(1.0+t1));
[97]990                temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
991                temp /= sqrt(1.0+t1);
992               
993                Sq += temp;
994        }
995       
996        Sq *= 2.0;
997        Sq += 1.0;
998       
999        inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval);
1000       
1001        inten *= 1.0e8;         // 1/A to 1/cm
1002       
1003    return(inten+bkg);
1004}
1005
1006
1007/*      LamellarPS_HGX  :  calculates the form factor of a lamellar structure - with S(q) effects included
[356]1008--- now the proper resolution effects are used - the "default" resolution is turned off (= 0) and the
1009model is smeared just like any other function
1010*/
[97]1011double
1012LamellarPS_HG(double dp[], double q)
1013{
1014        double scale,dd,delT,delH,SLD_T,SLD_H,SLD_S,NN,Cp,bkg;          //local variables of coefficient wave
1015        double inten,qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ,drh,drt;
1016        double Pi,Euler,dQDefault,fii;
1017        int ii,NNint;
1018       
1019       
1020        Euler = 0.5772156649;           // Euler's constant
[356]1021//      dQDefault = 0.0025;             //[=] 1/A, q-resolution, default value
1022        dQDefault = 0.0;
[97]1023        dQ = dQDefault;
1024       
1025        Pi = 4.0*atan(1.0);
1026        qval= q;
1027       
1028        scale = dp[0];
1029        dd = dp[1];
1030        delT = dp[2];
1031        delH = dp[3];
1032        SLD_T = dp[4];
1033        SLD_H = dp[5];
1034        SLD_S = dp[6];
1035        NN = trunc(dp[7]);              //be sure that NN is an integer
1036        Cp = dp[8];
1037        bkg = dp[9];
1038       
1039       
1040        drh = SLD_H - SLD_S;
1041        drt = SLD_T - SLD_S;    //correction 13FEB06 by L.Porcar
1042       
1043        Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT);
1044        Pq *= Pq;
1045        Pq *= 4.0/(qval*qval);
1046       
1047        NNint = (int)NN;                //cast to an integer for the loop
1048        ii=0;
1049        Sq = 0.0;
1050        for(ii=1;ii<(NNint-1);ii+=1) {
1051       
1052                fii = (double)ii;               //do I really need to do this?
1053               
1054                temp = 0.0;
1055                alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler);
1056                t1 = 2.0*dQ*dQ*dd*dd*alpha;
1057                t2 = 2.0*qval*qval*dd*dd*alpha;
1058                t3 = dQ*dQ*dd*dd*ii*ii;
1059               
1060                temp = 1.0-ii/NN;
1061                temp *= cos(dd*qval*ii/(1.0+t1));
1062                temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
1063                temp /= sqrt(1.0+t1);
1064               
1065                Sq += temp;
1066        }
1067       
1068        Sq *= 2.0;
1069        Sq += 1.0;
1070       
1071        inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval);
1072       
1073        inten *= 1.0e8;         // 1/A to 1/cm
1074       
1075        return(inten+bkg);
1076}
1077
1078/*      LamellarFF_HGX  :  calculates the form factor of a lamellar structure - no S(q) effects included
1079but extra SLD for head groups is included
1080
1081*/
1082double
1083LamellarFF_HG(double dp[], double q)
1084{
1085        double scale,delT,delH,slds,sldh,sldt,bkg;              //local variables of coefficient wave
1086        double inten, qval,Pq,drh,drt;
1087        double Pi;
1088       
1089       
1090        Pi = 4.0*atan(1.0);
1091        qval= q;
1092        scale = dp[0];
1093        delT = dp[1];
1094        delH = dp[2];
1095        sldt = dp[3];
1096        sldh = dp[4];
1097        slds = dp[5];
1098        bkg = dp[6];
1099       
1100       
1101        drh = sldh - slds;
1102        drt = sldt - slds;              //correction 13FEB06 by L.Porcar
1103       
1104        Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT);
1105        Pq *= Pq;
1106        Pq *= 4.0/(qval*qval);
1107       
1108        inten = 2.0*Pi*scale*Pq/(qval*qval);            //dimensionless...
1109       
1110        inten /= 2.0*(delT+delH);                       //normalize by the bilayer thickness
1111       
1112        inten *= 1.0e8;         // 1/A to 1/cm
1113       
1114        return(inten+bkg);
1115}
1116
1117/*      FlexExclVolCylX  :  calculates the form factor of a flexible cylinder with a circular cross section
1118-- incorporates Wei-Ren Chen's fixes - 2006
1119
1120        */
1121double
1122FlexExclVolCyl(double dp[], double q)
1123{
[235]1124        double scale,L,B,bkg,rad,qr,cont,sldc,slds;
[97]1125        double Pi,flex,crossSect,answer;
1126       
1127       
1128        Pi = 4.0*atan(1.0);
1129       
1130        scale = dp[0];          //make local copies in case memory moves
1131        L = dp[1];
1132        B = dp[2];
1133        rad = dp[3];
[235]1134        sldc = dp[4];
1135        slds = dp[5];
1136        cont = sldc-slds;
1137        bkg = dp[6];
[97]1138       
1139   
1140        qr = q*rad;
1141       
1142        flex = Sk_WR(q,L,B);
1143       
1144        crossSect = (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
1145        flex *= crossSect;
1146        flex *= Pi*rad*rad*L;
1147        flex *= cont*cont;
1148        flex *= 1.0e8;
1149        answer = scale*flex + bkg;
1150       
1151        return answer;
1152}
1153
1154/*      FlexCyl_EllipX  :  calculates the form factor of a flexible cylinder with an elliptical cross section
1155-- incorporates Wei-Ren Chen's fixes - 2006
1156
1157        */
1158double
1159FlexCyl_Ellip(double dp[], double q)
1160{
[235]1161        double scale,L,B,bkg,rad,qr,cont,ellRatio,slds,sldc;
[97]1162        double Pi,flex,crossSect,answer;
1163       
1164       
1165        Pi = 4.0*atan(1.0);
1166        scale = dp[0];          //make local copies in case memory moves
1167        L = dp[1];
1168        B = dp[2];
1169        rad = dp[3];
1170        ellRatio = dp[4];
[235]1171        sldc = dp[5];
1172        slds = dp[6];
1173        bkg = dp[7];
[97]1174       
[235]1175        cont = sldc - slds;
[97]1176        qr = q*rad;
1177       
1178        flex = Sk_WR(q,L,B);
1179       
1180        crossSect = EllipticalCross_fn(q,rad,(rad*ellRatio));
1181        flex *= crossSect;
1182        flex *= Pi*rad*rad*ellRatio*L;
1183        flex *= cont*cont;
1184        flex *= 1.0e8;
1185        answer = scale*flex + bkg;
1186       
1187        return answer;
1188}
1189
1190double
1191EllipticalCross_fn(double qq, double a, double b)
1192{
1193    double uplim,lolim,Pi,summ,arg,zi,yyy,answer;
1194    int i,nord=76;
1195   
1196    Pi = 4.0*atan(1.0);
1197    lolim=0.0;
1198    uplim=Pi/2.0;
1199    summ=0.0;
1200   
1201    for(i=0;i<nord;i++) {
1202                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1203                arg = qq*sqrt(a*a*sin(zi)*sin(zi)+b*b*cos(zi)*cos(zi));
1204                yyy = pow((2.0 * NR_BessJ1(arg) / arg),2);
1205                yyy *= Gauss76Wt[i];
1206                summ += yyy;
1207    }
1208    answer = (uplim-lolim)/2.0*summ;
1209    answer *= 2.0/Pi;
1210    return(answer);
1211       
1212}
1213/*      FlexCyl_PolyLenX  :  calculates the form factor of a flecible cylinder at the given x-value p->x
1214the cylinder has a polydisperse Length
1215
1216*/
1217double
1218FlexCyl_PolyLen(double dp[], double q)
1219{
1220        int i;
[235]1221        double scale,radius,length,pd,bkg,lb,delrho,sldc,slds;          //local variables of coefficient wave
[97]1222        int nord=20;                    //order of integration
1223        double uplim,lolim;             //upper and lower integration limits
1224        double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
1225        double range,zz,Pi;
1226       
1227        Pi = 4.0*atan(1.0);
1228        range = 3.4;
1229       
1230        summ = 0.0;                     //initialize intergral
1231        scale = dp[0];                  //make local copies in case memory moves
1232        length = dp[1];                 //radius
1233        pd = dp[2];                     // average length
1234        lb = dp[3];
1235        radius = dp[4];
[235]1236        sldc = dp[5];
1237        slds = dp[6];
1238        bkg = dp[7];
[97]1239       
[235]1240        delrho = sldc - slds;
[97]1241        zz = (1.0/pd)*(1.0/pd) - 1.0;
1242       
1243        lolim = length*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
[632]1244        if(lolim<0.0) {
1245                lolim = 0.0;
[97]1246        }
1247        if(pd>0.3) {
1248                range = 3.4 + (pd-0.3)*18.0;
1249        }
1250        uplim = length*(1.0+range*pd);
1251       
1252        for(i=0;i<nord;i++) {
1253                zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1254                yyy = Gauss20Wt[i] * FlePolyLen_kernel(q,radius,length,lb,zz,delrho,zi);
1255                summ += yyy;
1256        }
1257       
1258        answer = (uplim-lolim)/2.0*summ;
1259        //normalize by average cylinder volume (first moment), using the average length
1260        Vpoly=Pi*radius*radius*length;
1261        answer /= Vpoly;
1262       
1263        answer *=delrho*delrho;
1264       
1265        //convert to [cm-1]
1266        answer *= 1.0e8;
1267        //Scale
1268        answer *= scale;
1269        // add in the background
1270        answer += bkg;
1271       
1272        return answer;
1273}
1274
1275/*      FlexCyl_PolyLenX  :  calculates the form factor of a flexible cylinder at the given x-value p->x
1276the cylinder has a polydisperse cross sectional radius
1277
1278*/
1279double
1280FlexCyl_PolyRad(double dp[], double q)
1281{
1282        int i;
[235]1283        double scale,radius,length,pd,delrho,bkg,lb,sldc,slds;          //local variables of coefficient wave
[97]1284        int nord=76;                    //order of integration
1285        double uplim,lolim;             //upper and lower integration limits
1286        double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
1287        double range,zz,Pi;
1288       
1289       
1290        Pi = 4.0*atan(1.0);
1291        range = 3.4;
1292       
1293        summ = 0.0;                     //initialize intergral
1294       
1295        scale = dp[0];                  //make local copies in case memory moves
1296        length = dp[1];                 //radius
1297        lb = dp[2];                     // average length
1298        radius = dp[3];
1299        pd = dp[4];
[235]1300        sldc = dp[5];
1301        slds = dp[6];
1302        bkg = dp[7];
[97]1303       
[235]1304        delrho = sldc-slds;
[97]1305        zz = (1.0/pd)*(1.0/pd) - 1.0;
1306       
1307        lolim = radius*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
[632]1308        if(lolim<0.0) {
1309                lolim = 0.0;
[97]1310        }
1311        if(pd>0.3) {
1312                range = 3.4 + (pd-0.3)*18.0;
1313        }
1314        uplim = radius*(1.0+range*pd);
1315       
1316        for(i=0;i<nord;i++) {
1317                //zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1318                //yyy = Gauss20Wt[i] * FlePolyRad_kernel(q,radius,length,lb,zz,delrho,zi);
1319                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1320                yyy = Gauss76Wt[i] * FlePolyRad_kernel(q,radius,length,lb,zz,delrho,zi);
1321                summ += yyy;
1322        }
1323       
1324        answer = (uplim-lolim)/2.0*summ;
1325        //normalize by average cylinder volume (second moment), using the average radius
1326        Vpoly = Pi*radius*radius*length*(zz+2.0)/(zz+1.0);
1327        answer /= Vpoly;
1328       
1329        answer *=delrho*delrho;
1330       
1331        //convert to [cm-1]
1332        answer *= 1.0e8;
1333        //Scale
1334        answer *= scale;
1335        // add in the background
1336        answer += bkg;
1337       
1338        return answer;
1339}
1340
1341/////////functions for WRC implementation of flexible cylinders
1342static double
1343Sk_WR(double q, double L, double b)
1344{
1345        //
1346        double p1,p2,p1short,p2short,q0,qconnect;
1347        double C,epsilon,ans,q0short,Sexvmodify,pi;
1348       
1349        pi = 4.0*atan(1.0);
1350       
1351        p1 = 4.12;
1352        p2 = 4.42;
1353        p1short = 5.36;
1354        p2short = 5.62;
1355        q0 = 3.1;
1356        qconnect = q0/b;
1357        //     
1358        q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0);
1359       
1360        //
1361        if(L/b > 10.0) {
1362                C = 3.06/pow((L/b),0.44);
1363                epsilon = 0.176;
1364        } else {
1365                C = 1.0;
1366                epsilon = 0.170;
1367        }
1368        //
1369       
1370        if( L > 4*b ) { // Longer Chains
1371                if (q*b <= 3.1) {               //Modified by Yun on Oct. 15,
1372                        Sexvmodify = Sexvnew(q, L, b);
1373                        ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) - (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L);
1374                } else { //q(i)*b > 3.1
1375                        ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) + a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + pi/(q*L);
1376                } 
1377        } else { //L <= 4*b Shorter Chains
1378                if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) {
1379                        if (q*b<=0.01) {
1380                                ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0;
1381                        } else {
1382                                ans = Sdebye1(q,L,b);
1383                        }
1384                } else {        //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3)
1385                        ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) + a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) + pi/(q*L);
1386                }
1387        }
1388       
1389        return(ans);
1390        //return(a2long(q, L, b, p1, p2, q0));
1391}
1392
1393//WR named this w (too generic)
1394static double
1395w_WR(double x)
1396{
1397    double yy;
1398    yy = 0.5*(1 + tanh((x - 1.523)/0.1477));
1399       
1400    return (yy);
1401}
1402
1403//
1404static double
1405u1(double q, double L, double b)
1406{
1407    double yy;
1408       
1409    yy = Rgsquareshort(q,L,b)*q*q;
1410   
1411    return (yy);
1412}
1413
1414// was named u
1415static double
1416u_WR(double q, double L, double b)
1417{
1418    double yy;
1419    yy = Rgsquare(q,L,b)*q*q;
1420    return (yy);
1421}
1422
1423
1424
1425//
1426static double
1427Rgsquarezero(double q, double L, double b)
1428{   
1429    double yy;
1430    yy = (L*b/6.0) * (1.0 - 1.5*(b/L) + 1.5*pow((b/L),2) - 0.75*pow((b/L),3)*(1.0 - exp(-2.0*(L/b))));
1431   
1432    return (yy);
1433}
1434
1435//
1436static double
1437Rgsquareshort(double q, double L, double b)
1438{   
1439    double yy;
1440    yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b);
1441   
1442    return (yy);
1443}
1444
1445//
1446static double
1447Rgsquare(double q, double L, double b)
1448{
1449    double yy;
1450    yy = AlphaSquare(L/b)*L*b/6.0;
1451   
1452    return (yy);
1453}
1454
1455//
1456static double
1457AlphaSquare(double x)
1458{   
1459    double yy;
1460    yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) );
1461       
1462    return (yy);
1463}
1464
1465// ?? funciton is not used - but should the log actually be log10???
1466static double
1467miu(double x)
1468{   
1469    double yy;
1470    yy = (1.0/8.0)*(9.0*x - 2.0 + 2.0*log(1.0 + x)/x)*exp(1.0/2.565*(1.0/x + (1.0 - 1.0/(x*x))*log(1.0 + x)));
1471   
1472    return (yy);
1473}
1474
1475//
1476static double
1477Sdebye(double q, double L, double b)
1478{   
1479    double yy;
1480    yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2));
1481       
1482    return (yy);
1483}
1484
1485//
1486static double
1487Sdebye1(double q, double L, double b)
1488{   
1489    double yy;
1490    yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) );
1491   
1492    return (yy);
1493}
1494
1495//
1496static double
1497Sexv(double q, double L, double b)
1498{   
1499    double yy,C1,C2,C3,miu,Rg2;
1500    C1=1.22;
1501    C2=0.4288;
1502    C3=-1.651;
1503    miu = 0.585;
1504       
1505    Rg2 = Rgsquare(q,L,b);
1506   
1507    yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) +  C2*pow((q*sqrt(Rg2)),(-2.0/miu)) +    C3*pow((q*sqrt(Rg2)),(-3.0/miu)));
1508   
1509    return (yy);
1510}
1511
1512// this must be WR modified version
1513static double
1514Sexvnew(double q, double L, double b)
1515{   
1516    double yy,C1,C2,C3,miu;
1517    double del=1.05,C_star2,Rg2;
1518   
1519    C1=1.22;
1520    C2=0.4288;
1521    C3=-1.651;
1522    miu = 0.585;
1523       
1524    //calculating the derivative to decide on the corection (cutoff) term?
1525    // I have modified this from WRs original code
1526   
1527    if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) {
1528        C_star2 = 0.0;
1529    } else {
1530        C_star2 = 1.0;
1531    }
1532       
1533    Rg2 = Rgsquare(q,L,b);
1534   
1535        yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + C_star2*w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) + C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu)));
1536       
1537    return (yy);
1538}
1539
1540// these are the messy ones
1541static double
1542a2short(double q, double L, double b, double p1short, double p2short, double q0)
1543{   
1544    double yy,Rg2_sh;
1545    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p;
1546    double pi;
1547       
1548    E = 2.718281828459045091;
1549        pi = 4.0*atan(1.0);
1550    Rg2_sh = Rgsquareshort(q,L,b);
1551    Rg2_sh2 = Rg2_sh*Rg2_sh;
1552    t1 = ((q0*q0*Rg2_sh)/(b*b));
1553    Et1 = pow(E,t1);
1554    Emt1 =pow(E,-t1); 
1555    q02 = q0*q0;
1556    q0p = pow(q0,(-4.0 + p2short) );
1557   
1558    //E is the number e
1559    yy = ((-(1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b*b*b*L - 8.0*b*b*b*Et1*L - 2.0*b*b*b*L*p1short + 2.0*b*b*b*Et1*L*p1short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p1short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p1short*pi*q02*q0*Rg2_sh2)))))));
1560       
1561    return (yy);
1562}
1563
1564//
1565static double
1566a1short(double q, double L, double b, double p1short, double p2short, double q0)
1567{   
1568    double yy,Rg2_sh;
1569    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3;
1570        double pi;
1571   
1572    E = 2.718281828459045091;
1573        pi = 4.0*atan(1.0);
1574    Rg2_sh = Rgsquareshort(q,L,b);
1575    Rg2_sh2 = Rg2_sh*Rg2_sh;
1576    b3 = b*b*b;
1577    t1 = ((q0*q0*Rg2_sh)/(b*b));
1578    Et1 = pow(E,t1);
1579    Emt1 =pow(E,-t1); 
1580    q02 = q0*q0;
1581    q0p = pow(q0,(-4.0 + p1short) );
[632]1582
1583    yy = ((1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short + 2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p2short*pi*q02*q0*Rg2_sh2))))));
1584
[97]1585    return(yy);
1586}
1587
1588// this one will be lots of trouble
1589static double
1590a2long(double q, double L, double b, double p1, double p2, double q0)
1591{
1592    double yy,C1,C2,C3,C4,C5,miu,C,Rg2;
1593    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi;
1594    double E,b2,b3,b4,q02,q03,q04,q05,Rg22;
1595   
1596    pi = 4.0*atan(1.0);
1597    E = 2.718281828459045091;
1598    if( L/b > 10.0) {
1599        C = 3.06/pow((L/b),0.44);
1600    } else {
1601        C = 1.0;
1602    }
1603       
1604    C1 = 1.22;
1605    C2 = 0.4288;
1606    C3 = -1.651;
1607    C4 = 1.523;
1608    C5 = 0.1477;
1609    miu = 0.585;
1610       
1611    Rg2 = Rgsquare(q,L,b);
1612    Rg22 = Rg2*Rg2;
1613    b2 = b*b;
1614    b3 = b*b*b;
1615    b4 = b3*b;
1616    q02 = q0*q0;
1617    q03 = q0*q0*q0;
1618    q04 = q03*q0;
1619    q05 = q04*q0;
1620   
1621    t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) ));
[632]1622
1623    t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7*b2)/(15.0*q02*Rg2)))*Rg2)/b)))/L;
1624
1625    t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2.0))/(2.0*C5);
1626
1627    t4 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22);
1628
1629    t5 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);
1630
1631    t6 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22);
1632
1633    t7 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu));
1634
1635    t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
1636
1637    t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))))/L;
1638
1639    t10 = (2.0*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);
1640
1641
1642    yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))));
1643
[97]1644    return (yy);
1645}
1646
1647//need to define this on my own
1648static double
1649sech_WR(double x)
1650{       
1651        return(1/cosh(x));
1652}
1653
1654//
1655static double
1656a1long(double q, double L, double b, double p1, double p2, double q0)
1657{   
1658    double yy,C,C1,C2,C3,C4,C5,miu,Rg2;
1659    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15;
1660    double E,pi;
1661    double b2,b3,b4,q02,q03,q04,q05,Rg22;
1662       
1663    pi = 4.0*atan(1.0);
1664    E = 2.718281828459045091;
1665   
1666    if( L/b > 10.0) {
1667        C = 3.06/pow((L/b),0.44);
1668    } else {
1669        C = 1.0;
1670    }
1671       
1672    C1 = 1.22;
1673    C2 = 0.4288;
1674    C3 = -1.651;
1675    C4 = 1.523;
1676    C5 = 0.1477;
1677    miu = 0.585;
1678       
1679    Rg2 = Rgsquare(q,L,b);
1680    Rg22 = Rg2*Rg2;
1681    b2 = b*b;
1682    b3 = b*b*b;
1683    b4 = b3*b;
1684    q02 = q0*q0;
1685    q03 = q0*q0*q0;
1686    q04 = q03*q0;
1687    q05 = q04*q0;
1688   
1689    t1 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))));
1690   
1691    t2 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
1692   
1693    t3 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))));
1694   
1695    t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
1696       
1697    t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2))));
1698   
1699    t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2)))*Rg2)/b)));
1700   
1701    t7 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));
1702   
1703    t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));
1704   
1705    t9 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
1706   
1707    t10 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
1708   
1709    t11 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu));
1710   
1711    t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
1712   
1713    t13 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02* Rg2))) + (7.0*b2)/(15.0*q02*Rg2))));
1714   
1715    t14 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
1716   
1717    t15 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))));
1718       
1719   
1720    yy = (pow(q0,p1)*(((-((b*pi)/(L*q0))) +t1/L +t2/(q04*Rg22) + 1.0/2.0*t3*t4)) + (t5*((pow(q0,(p1 - p2))*(((-pow(q0,(-p1)))*(((b2*pi)/(L*q02) +t6/L +t7/(2.0*C5) -t8/(C5*q04*Rg22) +t9/(q04*Rg22) -t10/(q05*Rg22) + 1.0/2.0*t11*t12)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) +t13/L +t14/(q04*Rg22) + 1.0/2.0*t15*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))));
1721       
1722    return (yy);
1723}
1724
1725
1726
1727///////////////
1728
1729//
1730//     FUNCTION gfn2:    CONTAINS F(Q,A,B,mu)**2  AS GIVEN
1731//                       BY (53) AND (56,57) IN CHEN AND
1732//                       KOTLARCHYK REFERENCE
1733//
1734//     <PROLATE ELLIPSOIDS>
1735//
1736double
1737gfn2(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq)
1738{
1739        // local variables
[632]1740        double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,gfn2,pi43,gfn,Pi;
1741
[97]1742        Pi = 4.0*atan(1.0);
1743       
1744        pi43=4.0/3.0*Pi;
1745        aa = crmaj;
1746        bb = crmin;
1747        u2 = (aa*aa*xx*xx + bb*bb*(1.0-xx*xx));
1748        ut2 = (trmaj*trmaj*xx*xx + trmin*trmin*(1.0-xx*xx));
1749        uq = sqrt(u2)*qq;
1750        ut= sqrt(ut2)*qq;
1751        vc = pi43*aa*bb*bb;
1752        vt = pi43*trmaj*trmin*trmin;
[632]1753        if (uq == 0.0){
1754                siq = 1.0/3.0;
1755        }else{
1756                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq;
1757        }
1758        if (ut == 0.0){
1759                sit = 1.0/3.0;
1760        }else{
1761                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut;
1762        }
1763        gfnc = 3.0*siq*vc*delpc;
1764        gfnt = 3.0*sit*vt*delps;
[97]1765        gfn = gfnc+gfnt;
1766        gfn2 = gfn*gfn;
1767       
1768        return (gfn2);
1769}
1770
1771//
1772//     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN
1773//                       BY (53) & (58-59) IN CHEN AND
1774//                       KOTLARCHYK REFERENCE
1775//
1776//       <OBLATE ELLIPSOID>
1777// function gfn4 for oblate ellipsoids
1778double
1779gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq)
1780{
1781        // local variables
[632]1782        double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,tgfn,gfn4,pi43,Pi;
1783
[97]1784        Pi = 4.0*atan(1.0);
1785        pi43=4.0/3.0*Pi;
1786        aa = crmaj;
1787        bb = crmin;
1788        u2 = (bb*bb*xx*xx + aa*aa*(1.0-xx*xx));
1789        ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx));
1790        uq = sqrt(u2)*qq;
1791        ut= sqrt(ut2)*qq;
1792        vc = pi43*aa*aa*bb;
1793        vt = pi43*trmaj*trmaj*trmin;
[632]1794        if (uq == 0.0){
1795                siq = 1.0/3.0;
1796        }else{
1797                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq;
1798        }
1799        if (ut == 0.0){
1800                sit = 1.0/3.0;
1801        }else{
1802                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut;
1803        }
1804        gfnc = 3.0*siq*vc*delpc;
1805        gfnt = 3.0*sit*vt*delps;
[97]1806        tgfn = gfnc+gfnt;
1807        gfn4 = tgfn*tgfn;
1808       
1809        return (gfn4);
1810}
1811
1812double
1813FlePolyLen_kernel(double q, double radius, double length, double lb, double zz, double delrho, double zi)
1814{
1815    double Pq,vcyl,dl;
1816    double Pi,qr;
1817   
1818    Pi = 4.0*atan(1.0);
1819    qr = q*radius;
1820   
1821    Pq = Sk_WR(q,zi,lb);                //does not have cross section term
[632]1822    if (qr !=0){
1823        Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
1824    } 
[97]1825    vcyl=Pi*radius*radius*zi;
1826    Pq *= vcyl*vcyl;
1827   
1828    dl = SchulzPoint_cpr(zi,length,zz);
1829    return (Pq*dl);     
1830       
1831}
1832
1833double
1834FlePolyRad_kernel(double q, double ravg, double Lc, double Lb, double zz, double delrho, double zi)
1835{
1836    double Pq,vcyl,dr;
1837    double Pi,qr;
1838   
1839    Pi = 4.0*atan(1.0);
1840    qr = q*zi;
1841   
1842    Pq = Sk_WR(q,Lc,Lb);                //does not have cross section term
[632]1843    if (qr !=0){
1844        Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
1845    }
1846
[97]1847    vcyl=Pi*zi*zi*Lc;
1848    Pq *= vcyl*vcyl;
1849   
1850    dr = SchulzPoint_cpr(zi,ravg,zz);
1851    return (Pq*dr);     
1852       
1853}
1854
1855double
1856CSCylIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length)
1857{
1858        double answer,halfheight,Pi;
1859        double lolim,uplim,summ,yyy,zi;
1860        int nord,i;
1861       
1862        // set up the integration end points
1863        Pi = 4.0*atan(1.0);
1864        nord = 76;
[632]1865        lolim = 0.0;
1866        uplim = Pi/2.0;
[97]1867        halfheight = length/2.0;
1868       
1869        summ = 0.0;                             // initialize integral
1870        i=0;
1871        for(i=0;i<nord;i++) {
1872                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1873                yyy = Gauss76Wt[i] * CScyl(qq, rad, radthick, facthick, rhoc,rhos,rhosolv, halfheight, zi);
1874                summ += yyy;
1875        }
1876       
1877        // calculate value of integral to return
1878        answer = (uplim-lolim)/2.0*summ;
1879        return (answer);
1880}
1881
1882double
1883CScyl(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length, double dum)
1884{       
1885        // qq is the q-value for the calculation (1/A)
1886        // radius is the core radius of the cylinder (A)
1887        //  radthick and facthick are the radial and face layer thicknesses
1888        // rho(n) are the respective SLD's
1889        // length is the *Half* CORE-LENGTH of the cylinder
1890        // dum is the dummy variable for the integration (theta)
[632]1891
1892        double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval;
[97]1893        double Pi;
1894       
1895        Pi = 4.0*atan(1.0); 
1896       
1897        dr1 = rhoc-rhos;
1898        dr2 = rhos-rhosolv;
1899        vol1 = Pi*rad*rad*(2.0*length);
1900        vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick);
1901       
1902        besarg1 = qq*rad*sin(dum);
1903        besarg2 = qq*(rad+radthick)*sin(dum);
1904        sinarg1 = qq*length*cos(dum);
1905        sinarg2 = qq*(length+facthick)*cos(dum);
[632]1906        if (besarg1 == 0.0){
1907                be1 = 0.5;
1908        }else{
1909                be1 = NR_BessJ1(besarg1)/besarg1;
1910        }
1911        if (besarg2 == 0.0){
1912                be2 = 0.5;
1913        }else{
1914                be2 = NR_BessJ1(besarg2)/besarg2;
1915        }
1916        if (sinarg1 == 0.0){
1917                si1 = 1.0;
1918        }else{
1919                si1 = sin(sinarg1)/sinarg1;
1920        }
1921        if (besarg2 == 0.0){
1922                si2 = 1.0;
1923        }else{
1924                si2 = sin(sinarg2)/sinarg2;
1925        }
1926
1927        t1 = 2.0*vol1*dr1*si1*be1;
1928        t2 = 2.0*vol2*dr2*si2*be2;
1929
[97]1930        retval = ((t1+t2)*(t1+t2))*sin(dum);
1931        return (retval);
1932   
1933}
1934
1935
1936double
1937CoreShellCylKernel(double qq, double rcore, double thick, double rhoc, double rhos, double rhosolv, double length, double dum)
1938{
[632]1939
1940    double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval;
[97]1941    double Pi;
1942   
1943    Pi = 4.0*atan(1.0);
1944   
1945    dr1 = rhoc-rhos;
1946    dr2 = rhos-rhosolv;
1947    vol1 = Pi*rcore*rcore*(2.0*length);
1948    vol2 = Pi*(rcore+thick)*(rcore+thick)*(2.0*length+2.0*thick);
1949   
1950    besarg1 = qq*rcore*sin(dum);
1951    besarg2 = qq*(rcore+thick)*sin(dum);
1952    sinarg1 = qq*length*cos(dum);
1953    sinarg2 = qq*(length+thick)*cos(dum);
[632]1954
1955        if (besarg1 == 0.0){
1956                be1 = 0.5;
1957        }else{
1958                be1 = NR_BessJ1(besarg1)/besarg1;
1959        }
1960        if (besarg2 == 0.0){
1961                be2 = 0.5;
1962        }else{
1963                be2 = NR_BessJ1(besarg2)/besarg2;
1964        }
1965        if (sinarg1 == 0.0){
1966                si1 = 1.0;
1967        }else{
1968                si1 = sin(sinarg1)/sinarg1;
1969        }
1970        if (besarg2 == 0.0){
1971                si2 = 1.0;
1972        }else{
1973                si2 = sin(sinarg2)/sinarg2;
1974        }
1975
1976    t1 = 2.0*vol1*dr1*si1*be1;
1977    t2 = 2.0*vol2*dr2*si2*be2;
1978
[97]1979    retval = ((t1+t2)*(t1+t2))*sin(dum);
1980       
1981    return (retval);
1982}
1983
1984double
1985Cyl_PolyLenKernel(double q, double radius, double len_avg, double zz, double delrho, double dumLen)
1986{
1987       
1988    double halfheight,uplim,lolim,zi,summ,yyy,Pi;
1989    double answer,dr,Vcyl;
1990    int i,nord;
1991   
1992    Pi = 4.0*atan(1.0);
[632]1993    lolim = 0.0;
[97]1994    uplim = Pi/2.0;
1995    halfheight = dumLen/2.0;
1996    nord=20;
1997    summ = 0.0;
1998   
1999    //do the cylinder orientational average
2000    for(i=0;i<nord;i++) {
2001                zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
2002                yyy = Gauss20Wt[i] * CylKernel(q, radius, halfheight, zi);
2003                summ += yyy;
2004    }
2005    answer = (uplim-lolim)/2.0*summ;
2006    // Multiply by contrast^2
2007    answer *= delrho*delrho;
2008    // don't do the normal scaling to volume here
2009    // instead, multiply by VCyl^2 to get rid of the normalization for this radius of cylinder
2010    Vcyl = Pi*radius*radius*dumLen;
2011    answer *= Vcyl*Vcyl;
2012   
2013    dr = SchulzPoint_cpr(dumLen,len_avg,zz);
2014    return(dr*answer);
2015}
2016
2017
2018double
2019Stackdisc_kern(double qq, double rcore, double rhoc, double rhol, double rhosolv, double length, double thick, double dum, double gsd, double d, double N)
2020{               
2021        // qq is the q-value for the calculation (1/A)
2022        // rcore is the core radius of the cylinder (A)
2023        // rho(n) are the respective SLD's
2024        // length is the *Half* CORE-LENGTH of the cylinder = L (A)
2025        // dum is the dummy variable for the integration (x in Feigin's notation)
[632]2026
2027        //Local variables
2028        double totald,dr1,dr2,besarg1,besarg2,be1,be2,si1,si2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt;
[97]2029        double Pi;
2030        int kk;
2031       
2032        Pi = 4.0*atan(1.0);
2033       
2034        dr1 = rhoc-rhosolv;
2035        dr2 = rhol-rhosolv;
2036        area = Pi*rcore*rcore;
2037        totald=2.0*(thick+length);
2038       
2039        besarg1 = qq*rcore*sin(dum);
2040        besarg2 = qq*rcore*sin(dum);
2041       
2042        sinarg1 = qq*length*cos(dum);
2043        sinarg2 = qq*(length+thick)*cos(dum);
[632]2044
2045        if (besarg1 == 0.0){
2046                be1 = 0.5;
2047        }else{
2048                be1 = NR_BessJ1(besarg1)/besarg1;
2049        }
2050        if (besarg2 == 0.0){
2051                be2 = 0.5;
2052        }else{
2053                be2 = NR_BessJ1(besarg2)/besarg2;
2054        }
2055        if (sinarg1 == 0.0){
2056                si1 = 1.0;
2057        }else{
2058                si1 = sin(sinarg1)/sinarg1;
2059        }
2060        if (besarg2 == 0.0){
2061                si2 = 1.0;
2062        }else{
2063                si2 = sin(sinarg2)/sinarg2;
2064        }
2065
2066        t1 = 2.0*area*(2.0*length)*dr1*(si1)*(be1);
2067        t2 = 2.0*area*dr2*(totald*si2-2.0*length*si1)*(be2);
2068
[97]2069        retval =((t1+t2)*(t1+t2))*sin(dum);
2070       
2071        // loop for the structure facture S(q)
2072        sqq=0.0;
2073        for(kk=1;kk<N;kk+=1) {
2074                dexpt=qq*cos(dum)*qq*cos(dum)*d*d*gsd*gsd*kk/2.0;
2075                sqq=sqq+(N-kk)*cos(qq*cos(dum)*d*kk)*exp(-1.*dexpt);
2076        }                       
2077       
2078        // end of loop for S(q)
2079        sqq=1.0+2.0*sqq/N;
2080        retval *= sqq;
2081   
2082        return(retval);
2083}
2084
2085
2086double
2087Cyl_PolyRadKernel(double q, double radius, double length, double zz, double delrho, double dumRad)
2088{
2089       
2090    double halfheight,uplim,lolim,zi,summ,yyy,Pi;
2091    double answer,dr,Vcyl;
2092    int i,nord;
2093   
2094    Pi = 4.0*atan(1.0);
[632]2095    lolim = 0.0;
[97]2096    uplim = Pi/2.0;
2097    halfheight = length/2.0;
2098        //    nord=20;
2099    nord=76;
2100    summ = 0.0;
2101   
2102    //do the cylinder orientational average
2103        //    for(i=0;i<nord;i++) {
2104        //            zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
2105        //            yyy = Gauss20Wt[i] * CylKernel(q, dumRad, halfheight, zi);
2106        //            summ += yyy;
2107        //    }
2108    for(i=0;i<nord;i++) {
2109                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
2110                yyy = Gauss76Wt[i] * CylKernel(q, dumRad, halfheight, zi);
2111                summ += yyy;
2112    }
2113    answer = (uplim-lolim)/2.0*summ;
2114    // Multiply by contrast^2
2115    answer *= delrho*delrho;
2116    // don't do the normal scaling to volume here
2117    // instead, multiply by VCyl^2 to get rid of the normalization for this radius of cylinder
2118    Vcyl = Pi*dumRad*dumRad*length;
2119    answer *= Vcyl*Vcyl;
2120   
2121    dr = SchulzPoint_cpr(dumRad,radius,zz);
2122    return(dr*answer);
2123}
2124
2125double
2126SchulzPoint_cpr(double dumRad, double radius, double zz)
2127{
2128    double dr;
2129   
2130    dr = zz*log(dumRad) - gammaln(zz+1.0) + (zz+1.0)*log((zz+1.0)/radius)-(dumRad/radius*(zz+1.0));
2131    return(exp(dr));
2132}
2133
2134static double
2135gammaln(double xx)
2136{
2137        double x,y,tmp,ser;
2138        static double cof[6]={76.18009172947146,-86.50532032941677,
2139                24.01409824083091,-1.231739572450155,
2140                0.1208650973866179e-2,-0.5395239384953e-5};
2141        int j;
2142       
2143        y=x=xx;
2144        tmp=x+5.5;
2145        tmp -= (x+0.5)*log(tmp);
2146        ser=1.000000000190015;
2147        for (j=0;j<=5;j++) ser += cof[j]/++y;
2148        return -tmp+log(2.5066282746310005*ser/x);
2149}
2150
2151
2152double
2153EllipsoidKernel(double qq, double a, double nua, double dum)
2154{
2155    double arg,nu,retval;               //local variables
2156       
2157    nu = nua/a;
[632]2158    arg = qq*a*sqrt(1.0+dum*dum*(nu*nu-1.0));
2159    if (arg == 0.0){
2160        retval =1.0/3.0;
2161    }else{
2162        retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg);
2163    }
[97]2164    retval *= retval;
[632]2165    retval *= 9.0;
2166
[97]2167    return(retval);
2168}//Function EllipsoidKernel()
2169
2170double
2171HolCylKernel(double qq, double rcore, double rshell, double length, double dum)
2172{
2173    double gamma,arg1,arg2,lam1,lam2,psi,sinarg,t2,retval;              //local variables
2174   
2175    gamma = rcore/rshell;
[632]2176    arg1 = qq*rshell*sqrt(1.0-dum*dum);         //1=shell (outer radius)
2177    arg2 = qq*rcore*sqrt(1.0-dum*dum);                  //2=core (inner radius)
2178    if (arg1 == 0.0){
2179        lam1 = 1.0;
2180    }else{
2181        lam1 = 2.0*NR_BessJ1(arg1)/arg1;
2182    }
2183    if (arg2 == 0.0){
2184        lam2 = 1.0;
2185    }else{
2186        lam2 = 2.0*NR_BessJ1(arg2)/arg2;
2187    }
2188    //Todo: Need to check psi behavior as gamma goes to 1.
2189    psi = 1.0/(1.0-gamma*gamma)*(lam1 -  gamma*gamma*lam2);             //SRK 10/19/00
2190    sinarg = qq*length*dum/2.0;
2191    if (sinarg == 0.0){
2192        t2 = 1.0;
2193    }else{
2194        t2 = sin(sinarg)/sinarg;
2195    }
2196
[97]2197    retval = psi*psi*t2*t2;
2198   
2199    return(retval);
2200}//Function HolCylKernel()
2201
2202double
2203PPKernel(double aa, double mu, double uu)
2204{
2205    // mu passed in is really mu*sqrt(1-sig^2)
2206    double arg1,arg2,Pi,tmp1,tmp2;                      //local variables
2207       
2208    Pi = 4.0*atan(1.0);
2209   
2210    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
[632]2211    arg1 = (mu/2.0)*cos(Pi*uu/2.0);
2212    arg2 = (mu*aa/2.0)*sin(Pi*uu/2.0);
2213    if(arg1==0.0) {
2214                tmp1 = 1.0;
[97]2215        } else {
2216                tmp1 = sin(arg1)*sin(arg1)/arg1/arg1;
2217    }
[632]2218
2219    if (arg2==0.0) {
2220                tmp2 = 1.0;
[97]2221        } else {
2222                tmp2 = sin(arg2)*sin(arg2)/arg2/arg2;
2223    }
2224       
2225    return (tmp1*tmp2);
2226       
2227}//Function PPKernel()
2228
2229
2230double
2231TriaxialKernel(double q, double aa, double bb, double cc, double dx, double dy)
2232{
2233       
2234    double arg,val,pi;                  //local variables
2235       
2236    pi = 4.0*atan(1.0);
2237   
2238    arg = aa*aa*cos(pi*dx/2)*cos(pi*dx/2);
2239    arg += bb*bb*sin(pi*dx/2)*sin(pi*dx/2)*(1-dy*dy);
2240    arg += cc*cc*dy*dy;
2241    arg = q*sqrt(arg);
[632]2242    if (arg == 0.0){
2243        val = 1.0; // as arg --> 0, val should go to 1.0
2244    }else{
2245        val = 9.0 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) );
2246    }
[97]2247    return (val);
2248       
2249}//Function TriaxialKernel()
2250
2251
2252double
2253CylKernel(double qq, double rr,double h, double theta)
2254{
2255       
2256        // qq is the q-value for the calculation (1/A)
2257        // rr is the radius of the cylinder (A)
2258        // h is the HALF-LENGTH of the cylinder = L/2 (A)
[632]2259
2260    double besarg,bj,retval,d1,t1,b1,t2,b2,siarg,be,si;          //Local variables
2261
2262
[97]2263    besarg = qq*rr*sin(theta);
[632]2264    siarg = qq * h * cos(theta);
[97]2265    bj =NR_BessJ1(besarg);
2266       
2267        //* Computing 2nd power */
[632]2268    d1 = sin(siarg);
[97]2269    t1 = d1 * d1;
2270        //* Computing 2nd power */
2271    d1 = bj;
2272    t2 = d1 * d1 * 4.0 * sin(theta);
2273        //* Computing 2nd power */
[632]2274    d1 = siarg;
[97]2275    b1 = d1 * d1;
2276        //* Computing 2nd power */
2277    d1 = qq * rr * sin(theta);
2278    b2 = d1 * d1;
[632]2279    if (besarg == 0.0){
2280        be = sin(theta);
2281    }else{
2282        be = t2 / b2;
2283    }
2284    if (siarg == 0.0){
2285        si = 1.0;
2286    }else{
2287        si = t1 / b1;
2288    }
2289    retval = be * si;
2290
[97]2291    return (retval);
2292       
2293}//Function CylKernel()
2294
2295double
2296EllipCylKernel(double qq, double ra,double nu, double theta)
2297{
2298        //this is the function LAMBDA1^2 in Feigin's notation   
2299        // qq is the q-value for the calculation (1/A)
2300        // ra is the transformed radius"a" in Feigin's notation
2301        // nu is the ratio (major radius)/(minor radius) of the Ellipsoid [=] ---
2302        // theta is the dummy variable of the integration
2303       
2304        double retval,arg;               //Local variables
2305       
[632]2306        arg = qq*ra*sqrt((1.0+nu*nu)/2+(1.0-nu*nu)*cos(theta)/2);
2307        if (arg == 0.0){
2308                retval = 1.0;
2309        }else{
2310                retval = 2.0*NR_BessJ1(arg)/arg;
2311        }
2312
[97]2313        //square it
2314        retval *= retval;
2315       
2316    return(retval);
2317       
2318}//Function EllipCylKernel()
2319
2320double NR_BessJ1(double x)
2321{
2322        double ax,z;
2323        double xx,y,ans,ans1,ans2;
2324       
2325        if ((ax=fabs(x)) < 8.0) {
2326                y=x*x;
2327                ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
2328                                                                                                  +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
2329                ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
2330                                                                                           +y*(99447.43394+y*(376.9991397+y*1.0))));
2331                ans=ans1/ans2;
2332        } else {
2333                z=8.0/ax;
2334                y=z*z;
2335                xx=ax-2.356194491;
2336                ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
2337                                                                   +y*(0.2457520174e-5+y*(-0.240337019e-6))));
2338                ans2=0.04687499995+y*(-0.2002690873e-3
2339                                                          +y*(0.8449199096e-5+y*(-0.88228987e-6
2340                                                                                                         +y*0.105787412e-6)));
2341                ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
2342                if (x < 0.0) ans = -ans;
2343        }
2344       
2345        return(ans);
2346}
[453]2347
2348/*      Lamellar_ParaCrystal - Pedersen's model
2349
2350*/
2351double
2352Lamellar_ParaCrystal(double w[], double q)
2353{
2354//       Input (fitting) variables are:
2355        //[0] scale factor
2356        //[1]   thickness
2357        //[2]   number of layers
2358        //[3]   spacing between layers
2359        //[4]   polydispersity of spacing
2360        //[5] SLD lamellar
2361        //[6] SLD solvent
2362        //[7] incoherent background
2363//      give them nice names
2364        double inten,qval,scale,th,nl,davg,pd,contr,bkg,xn;
2365        double xi,ww,Pbil,Znq,Snq,an,sldLayer,sldSolvent,pi;
2366        long n1,n2;
2367       
2368        pi = 4.0*atan(1.0);
2369        scale = w[0];
2370        th = w[1];
2371        nl = w[2];
2372        davg = w[3];
2373        pd = w[4];
2374        sldLayer = w[5];
2375        sldSolvent = w[6];
2376        bkg = w[7];
2377       
2378        contr = w[5] - w[6];
2379        qval = q;
2380               
2381        //get the fractional part of nl, to determine the "mixing" of N's
2382       
2383        n1 = trunc(nl);         //rounds towards zero
2384        n2 = n1 + 1;
2385        xn = (double)n2 - nl;                   //fractional contribution of n1
2386       
2387        ww = exp(-qval*qval*pd*pd*davg*davg/2.0);
2388       
2389        //calculate the n1 contribution
2390        an = paraCryst_an(ww,qval,davg,n1);
2391        Snq = paraCryst_sn(ww,qval,davg,n1,an);
2392       
2393        Znq = xn*Snq;
2394       
2395        //calculate the n2 contribution
2396        an = paraCryst_an(ww,qval,davg,n2);
2397        Snq = paraCryst_sn(ww,qval,davg,n2,an);
2398       
2399        Znq += (1.0-xn)*Snq;
2400       
2401        //and the independent contribution
2402        Znq += (1.0-ww*ww)/(1.0+ww*ww-2.0*ww*cos(qval*davg));
2403       
2404        //the limit when NL approaches infinity
2405//      Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg))
2406       
2407        xi = th/2.0;            //use 1/2 the bilayer thickness
2408        Pbil = (sin(qval*xi)/(qval*xi))*(sin(qval*xi)/(qval*xi));
2409       
2410        inten = 2.0*pi*contr*contr*Pbil*Znq/(qval*qval);
2411        inten *= 1.0e8;
2412       
2413        return(scale*inten+bkg);
2414}
2415
2416// functions for the lamellar paracrystal model
2417double
2418paraCryst_sn(double ww, double qval, double davg, long nl, double an) {
2419       
2420        double Snq;
2421       
2422        Snq = an/( (double)nl*pow((1.0+ww*ww-2.0*ww*cos(qval*davg)),2) );
2423       
2424        return(Snq);
2425}
2426
2427
2428double
2429paraCryst_an(double ww, double qval, double davg, long nl) {
2430       
2431        double an;
2432       
2433        an = 4.0*ww*ww - 2.0*(ww*ww*ww+ww)*cos(qval*davg);
2434        an -= 4.0*pow(ww,(nl+2))*cos((double)nl*qval*davg);
2435        an += 2.0*pow(ww,(nl+3))*cos((double)(nl-1)*qval*davg);
2436        an += 2.0*pow(ww,(nl+1))*cos((double)(nl+1)*qval*davg);
2437       
2438        return(an);
2439}
2440
2441
2442/*      Spherocylinder  :
2443
2444Uses 76 pt Gaussian quadrature for both integrals
2445*/
2446double
2447Spherocylinder(double w[], double x)
2448{
2449        int i,j;
2450        double Pi;
2451        double scale,contr,bkg,sldc,slds;
2452        double len,rad,hDist,endRad;
2453        int nordi=76;                   //order of integration
2454        int nordj=76;
2455        double va,vb;           //upper and lower integration limits
2456        double summ,zi,yyy,answer;                      //running tally of integration
2457        double summj,vaj,vbj,zij;                       //for the inner integration
2458        double SphCyl_tmp[7],arg1,arg2,inner;
2459       
2460       
2461        scale = w[0];
2462        rad = w[1];
2463        len = w[2];
2464        sldc = w[3];
2465        slds = w[4];
2466        bkg = w[5];
2467       
2468        SphCyl_tmp[0] = w[0];
2469        SphCyl_tmp[1] = w[1];
2470        SphCyl_tmp[2] = w[2];
2471        SphCyl_tmp[3] = w[1];           //end radius is same as cylinder radius
2472        SphCyl_tmp[4] = w[3];
2473        SphCyl_tmp[5] = w[4];
2474        SphCyl_tmp[6] = w[5];
2475       
2476        hDist = 0;              //by definition for this model
2477        endRad = rad;
2478       
2479        contr = sldc-slds;
2480       
2481        Pi = 4.0*atan(1.0);
2482        va = 0.0;
2483        vb = Pi/2.0;            //orintational average, outer integral
2484        vaj = -1.0*hDist/endRad;
2485        vbj = 1.0;              //endpoints of inner integral
2486
2487        summ = 0.0;                     //initialize intergral
2488
2489        for(i=0;i<nordi;i++) {
2490                //setup inner integral over the ellipsoidal cross-section
2491                summj=0.0;
2492                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
2493               
2494                for(j=0;j<nordj;j++) {
2495                        //20 gauss points for the inner integral
2496                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
2497                        yyy = Gauss76Wt[j] * SphCyl_kernel(SphCyl_tmp,x,zij,zi);
2498                        summj += yyy;
2499                }
2500                //now calculate the value of the inner integral
2501                inner = (vbj-vaj)/2.0*summj;
2502                inner *= 4.0*Pi*endRad*endRad*endRad;
2503               
2504                //now calculate outer integrand
2505                arg1 = x*len/2.0*cos(zi);
2506                arg2 = x*rad*sin(zi);
2507                yyy = inner;
2508               
[632]2509                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
[453]2510                        yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2;
2511                } else {
2512                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2;
2513                }
2514                yyy *= yyy;
2515                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
2516                yyy *= Gauss76Wt[i];
2517                summ += yyy;
2518        }               //final scaling is done at the end of the function, after the NT_FP64 case
2519       
2520        answer = (vb-va)/2.0*summ;
2521
2522        answer /= Pi*rad*rad*len + Pi*4.0*endRad*endRad*endRad/3.0;             //divide by volume
2523        answer *= 1.0e8;                //convert to cm^-1
2524        answer *= contr*contr;
2525        answer *= scale;
2526        answer += bkg;
2527               
2528        return answer;
2529}
2530
2531
2532// inner integral of the sphereocylinder model, special case of hDist = 0
2533//
2534double
2535SphCyl_kernel(double w[], double x, double tt, double theta) {
2536
2537        double val,arg1,arg2;
2538        double scale,bkg,sldc,slds;
2539        double len,rad,hDist,endRad;
2540        scale = w[0];
2541        rad = w[1];
2542        len = w[2];
2543        endRad = w[3];
2544        sldc = w[4];
2545        slds = w[5];
2546        bkg = w[6];
2547       
2548        hDist = 0.0;
2549               
2550        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0);
2551        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt);
2552       
2553        val = cos(arg1)*(1.0-tt*tt)*NR_BessJ1(arg2)/arg2;
2554       
2555        return(val);
2556}
2557
2558
2559/*      Convex Lens  : special case where L ~ 0 and hDist < 0
2560
2561Uses 76 pt Gaussian quadrature for both integrals
2562*/
2563double
2564ConvexLens(double w[], double x)
2565{
2566        int i,j;
2567        double Pi;
2568        double scale,contr,bkg,sldc,slds;
2569        double len,rad,hDist,endRad;
2570        int nordi=76;                   //order of integration
2571        int nordj=76;
2572        double va,vb;           //upper and lower integration limits
2573        double summ,zi,yyy,answer;                      //running tally of integration
2574        double summj,vaj,vbj,zij;                       //for the inner integration
2575        double CLens_tmp[7],arg1,arg2,inner,hh;
2576       
2577       
2578        scale = w[0];
2579        rad = w[1];
2580//      len = w[2]
2581        endRad = w[2];
2582        sldc = w[3];
2583        slds = w[4];
2584        bkg = w[5];
2585       
2586        len = 0.01;
2587       
2588        CLens_tmp[0] = w[0];
2589        CLens_tmp[1] = w[1];
2590        CLens_tmp[2] = len;                     //length is some small number, essentially zero
2591        CLens_tmp[3] = w[2];
2592        CLens_tmp[4] = w[3];
2593        CLens_tmp[5] = w[4];
2594        CLens_tmp[6] = w[5];
2595               
2596        hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad));         //by definition for this model
2597       
2598        contr = sldc-slds;
2599       
2600        Pi = 4.0*atan(1.0);
2601        va = 0.0;
2602        vb = Pi/2.0;            //orintational average, outer integral
2603        vaj = -1.0*hDist/endRad;
2604        vbj = 1.0;              //endpoints of inner integral
2605
2606        summ = 0.0;                     //initialize intergral
2607
2608        for(i=0;i<nordi;i++) {
2609                //setup inner integral over the ellipsoidal cross-section
2610                summj=0.0;
2611                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
2612               
2613                for(j=0;j<nordj;j++) {
2614                        //20 gauss points for the inner integral
2615                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
2616                        yyy = Gauss76Wt[j] * ConvLens_kernel(CLens_tmp,x,zij,zi);
2617                        summj += yyy;
2618                }
2619                //now calculate the value of the inner integral
2620                inner = (vbj-vaj)/2.0*summj;
2621                inner *= 4.0*Pi*endRad*endRad*endRad;
2622               
2623                //now calculate outer integrand
2624                arg1 = x*len/2.0*cos(zi);
2625                arg2 = x*rad*sin(zi);
2626                yyy = inner;
2627               
[632]2628                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
[453]2629                        yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2;
2630                } else {
2631                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2;
2632                }
2633                yyy *= yyy;
2634                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
2635                yyy *= Gauss76Wt[i];
2636                summ += yyy;
2637        }               //final scaling is done at the end of the function, after the NT_FP64 case
2638       
2639        answer = (vb-va)/2.0*summ;
2640
2641        hh = fabs(hDist);               //need positive value for spherical cap volume
2642        answer /= 2.0*(1.0/3.0*Pi*(endRad-hh)*(endRad-hh)*(2.0*endRad+hh));             //divide by volume
2643        answer *= 1.0e8;                //convert to cm^-1
2644        answer *= contr*contr;
2645        answer *= scale;
2646        answer += bkg;
2647               
2648        return answer;
2649}
2650
2651/*      Capped Cylinder  : special case where L is nonzero and hDist < 0
2652
2653 -- uses the same Kernel as the Convex Lens
2654 
2655Uses 76 pt Gaussian quadrature for both integrals
2656*/
2657double
2658CappedCylinder(double w[], double x)
2659{
2660        int i,j;
2661        double Pi;
2662        double scale,contr,bkg,sldc,slds;
2663        double len,rad,hDist,endRad;
2664        int nordi=76;                   //order of integration
2665        int nordj=76;
2666        double va,vb;           //upper and lower integration limits
2667        double summ,zi,yyy,answer;                      //running tally of integration
2668        double summj,vaj,vbj,zij;                       //for the inner integration
2669        double arg1,arg2,inner,hh;
2670       
2671       
2672        scale = w[0];
2673        rad = w[1];
2674        len = w[2];
2675        endRad = w[3];
2676        sldc = w[4];
2677        slds = w[5];
2678        bkg = w[6];
2679               
2680        hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad));         //by definition for this model
2681       
2682        contr = sldc-slds;
2683       
2684        Pi = 4.0*atan(1.0);
2685        va = 0.0;
2686        vb = Pi/2.0;            //orintational average, outer integral
2687        vaj = -1.0*hDist/endRad;
2688        vbj = 1.0;              //endpoints of inner integral
2689
2690        summ = 0.0;                     //initialize intergral
2691
2692        for(i=0;i<nordi;i++) {
2693                //setup inner integral over the ellipsoidal cross-section
2694                summj=0.0;
2695                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
2696               
2697                for(j=0;j<nordj;j++) {
2698                        //20 gauss points for the inner integral
2699                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
2700                        yyy = Gauss76Wt[j] * ConvLens_kernel(w,x,zij,zi);               //uses the same kernel as ConvexLens, except here L != 0
2701                        summj += yyy;
2702                }
2703                //now calculate the value of the inner integral
2704                inner = (vbj-vaj)/2.0*summj;
2705                inner *= 4.0*Pi*endRad*endRad*endRad;
2706               
2707                //now calculate outer integrand
2708                arg1 = x*len/2.0*cos(zi);
2709                arg2 = x*rad*sin(zi);
2710                yyy = inner;
2711               
[632]2712                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
[453]2713                        yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2;
2714                } else {
2715                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2;
2716                }
2717                yyy *= yyy;
2718                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
2719                yyy *= Gauss76Wt[i];
2720                summ += yyy;
2721        }               //final scaling is done at the end of the function, after the NT_FP64 case
2722       
2723        answer = (vb-va)/2.0*summ;
2724
2725        hh = fabs(hDist);               //need positive value for spherical cap volume
2726        answer /= Pi*rad*rad*len + 2.0*(1.0/3.0*Pi*(endRad-hh)*(endRad-hh)*(2.0*endRad+hh));            //divide by volume
2727        answer *= 1.0e8;                //convert to cm^-1
2728        answer *= contr*contr;
2729        answer *= scale;
2730        answer += bkg;
2731               
2732        return answer;
2733}
2734
2735
2736
2737// inner integral of the ConvexLens model, special case where L ~ 0 and hDist < 0
2738//
2739double
2740ConvLens_kernel(double w[], double x, double tt, double theta) {
2741
2742        double val,arg1,arg2;
2743        double scale,bkg,sldc,slds;
2744        double len,rad,hDist,endRad;
2745        scale = w[0];
2746        rad = w[1];
2747        len = w[2];
2748        endRad = w[3];
2749        sldc = w[4];
2750        slds = w[5];
2751        bkg = w[6];
2752       
2753        hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad));
2754               
2755        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0);
2756        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt);
2757       
2758        val = cos(arg1)*(1.0-tt*tt)*NR_BessJ1(arg2)/arg2;
2759       
2760        return(val);
2761}
2762
2763
2764/*      Dumbbell  : special case where L ~ 0 and hDist > 0
2765
2766Uses 76 pt Gaussian quadrature for both integrals
2767*/
2768double
2769Dumbbell(double w[], double x)
2770{
2771        int i,j;
2772        double Pi;
2773        double scale,contr,bkg,sldc,slds;
2774        double len,rad,hDist,endRad;
2775        int nordi=76;                   //order of integration
2776        int nordj=76;
2777        double va,vb;           //upper and lower integration limits
2778        double summ,zi,yyy,answer;                      //running tally of integration
2779        double summj,vaj,vbj,zij;                       //for the inner integration
2780        double Dumb_tmp[7],arg1,arg2,inner;
2781       
2782       
2783        scale = w[0];
2784        rad = w[1];
2785//      len = w[2]
2786        endRad = w[2];
2787        sldc = w[3];
2788        slds = w[4];
2789        bkg = w[5];
2790       
2791        len = 0.01;
2792       
2793        Dumb_tmp[0] = w[0];
2794        Dumb_tmp[1] = w[1];
2795        Dumb_tmp[2] = len;              //length is some small number, essentially zero
2796        Dumb_tmp[3] = w[2];
2797        Dumb_tmp[4] = w[3];
2798        Dumb_tmp[5] = w[4];
2799        Dumb_tmp[6] = w[5];
2800                       
2801        hDist = sqrt(fabs(endRad*endRad-rad*rad));              //by definition for this model
2802       
2803        contr = sldc-slds;
2804       
2805        Pi = 4.0*atan(1.0);
2806        va = 0.0;
2807        vb = Pi/2.0;            //orintational average, outer integral
2808        vaj = -1.0*hDist/endRad;
2809        vbj = 1.0;              //endpoints of inner integral
2810
2811        summ = 0.0;                     //initialize intergral
2812
2813        for(i=0;i<nordi;i++) {
2814                //setup inner integral over the ellipsoidal cross-section
2815                summj=0.0;
2816                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
2817               
2818                for(j=0;j<nordj;j++) {
2819                        //20 gauss points for the inner integral
2820                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
2821                        yyy = Gauss76Wt[j] * Dumb_kernel(Dumb_tmp,x,zij,zi);
2822                        summj += yyy;
2823                }
2824                //now calculate the value of the inner integral
2825                inner = (vbj-vaj)/2.0*summj;
2826                inner *= 4.0*Pi*endRad*endRad*endRad;
2827               
2828                //now calculate outer integrand
2829                arg1 = x*len/2.0*cos(zi);
2830                arg2 = x*rad*sin(zi);
2831                yyy = inner;
2832               
[632]2833                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
[453]2834                        yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2;
2835                } else {
2836                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2;
2837                }
2838                yyy *= yyy;
2839                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
2840                yyy *= Gauss76Wt[i];
2841                summ += yyy;
2842        }               //final scaling is done at the end of the function, after the NT_FP64 case
2843       
2844        answer = (vb-va)/2.0*summ;
2845
2846        answer /= 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*hDist*hDist/3.0);              //divide by volume
2847        answer *= 1.0e8;                //convert to cm^-1
2848        answer *= contr*contr;
2849        answer *= scale;
2850        answer += bkg;
2851               
2852        return answer;
2853}
2854
2855
2856/*      Barbell  : "normal" case where L is nonzero 0 and hDist > 0
2857
2858-- uses the same kernel as the Dumbbell case
2859
2860Uses 76 pt Gaussian quadrature for both integrals
2861*/
2862double
2863Barbell(double w[], double x)
2864{
2865        int i,j;
2866        double Pi;
2867        double scale,contr,bkg,sldc,slds;
2868        double len,rad,hDist,endRad;
2869        int nordi=76;                   //order of integration
2870        int nordj=76;
2871        double va,vb;           //upper and lower integration limits
2872        double summ,zi,yyy,answer;                      //running tally of integration
2873        double summj,vaj,vbj,zij;                       //for the inner integration
2874        double arg1,arg2,inner;
2875       
2876       
2877        scale = w[0];
2878        rad = w[1];
2879        len = w[2];
2880        endRad = w[3];
2881        sldc = w[4];
2882        slds = w[5];
2883        bkg = w[6];
2884                       
2885        hDist = sqrt(fabs(endRad*endRad-rad*rad));              //by definition for this model
2886       
2887        contr = sldc-slds;
2888       
2889        Pi = 4.0*atan(1.0);
2890        va = 0.0;
2891        vb = Pi/2.0;            //orintational average, outer integral
2892        vaj = -1.0*hDist/endRad;
2893        vbj = 1.0;              //endpoints of inner integral
2894
2895        summ = 0.0;                     //initialize intergral
2896
2897        for(i=0;i<nordi;i++) {
2898                //setup inner integral over the ellipsoidal cross-section
2899                summj=0.0;
2900                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
2901               
2902                for(j=0;j<nordj;j++) {
2903                        //20 gauss points for the inner integral
2904                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
2905                        yyy = Gauss76Wt[j] * Dumb_kernel(w,x,zij,zi);           //uses the same Kernel as the Dumbbell, here L>0
2906                        summj += yyy;
2907                }
2908                //now calculate the value of the inner integral
2909                inner = (vbj-vaj)/2.0*summj;
2910                inner *= 4.0*Pi*endRad*endRad*endRad;
2911               
2912                //now calculate outer integrand
2913                arg1 = x*len/2.0*cos(zi);
2914                arg2 = x*rad*sin(zi);
2915                yyy = inner;
2916               
[632]2917                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
[453]2918                        yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2;
2919                } else {
2920                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2;
2921                }
2922                yyy *= yyy;
2923                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
2924                yyy *= Gauss76Wt[i];
2925                summ += yyy;
2926        }               //final scaling is done at the end of the function, after the NT_FP64 case
2927       
2928        answer = (vb-va)/2.0*summ;
2929
2930        answer /= Pi*rad*rad*len + 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*hDist*hDist/3.0);             //divide by volume
2931        answer *= 1.0e8;                //convert to cm^-1
2932        answer *= contr*contr;
2933        answer *= scale;
2934        answer += bkg;
2935               
2936        return answer;
2937}
2938
2939
2940
2941// inner integral of the Dumbbell model, special case where L ~ 0 and hDist > 0
2942//
2943// inner integral of the Barbell model if L is nonzero
2944//
2945double
2946Dumb_kernel(double w[], double x, double tt, double theta) {
2947
2948        double val,arg1,arg2;
2949        double scale,bkg,sldc,slds;
2950        double len,rad,hDist,endRad;
2951        scale = w[0];
2952        rad = w[1];
2953        len = w[2];
2954        endRad = w[3];
2955        sldc = w[4];
2956        slds = w[5];
2957        bkg = w[6];
2958       
2959        hDist = sqrt(fabs(endRad*endRad-rad*rad));
2960               
2961        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0);
2962        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt);
2963       
2964        val = cos(arg1)*(1.0-tt*tt)*NR_BessJ1(arg2)/arg2;
2965       
2966        return(val);
2967}
[500]2968
[501]2969double PolyCoreBicelle(double dp[], double q)
2970{
[500]2971        int i;
2972        int nord = 20;
2973        double scale, length, sigma, bkg, radius, radthick, facthick;
2974        double rhoc, rhoh, rhor, rhosolv;
2975        double answer, Vpoly;
[501]2976        double Pi,lolim,uplim,summ,yyy,rad,AR,Rsqr,Rsqrsumm,Rsqryyy;
[500]2977       
[501]2978        scale = dp[0];
2979        radius = dp[1];
2980        sigma = dp[2];                          //sigma is the standard mean deviation
2981        length = dp[3];
2982        radthick = dp[4];
2983        facthick= dp[5];
2984        rhoc = dp[6];
2985        rhoh = dp[7];
2986        rhor=dp[8];
2987        rhosolv = dp[9];
2988        bkg = dp[10];
[500]2989       
[501]2990        Pi = 4.0*atan(1.0);
2991       
[500]2992        lolim = exp(log(radius)-(4.*sigma));
[632]2993        if (lolim<0.0) {
2994                lolim=0.0;              //to avoid numerical error when  va<0 (-ve r value)
[500]2995        }
2996        uplim = exp(log(radius)+(4.*sigma));
2997       
[501]2998        summ = 0.0;
2999        Rsqrsumm = 0.0;
3000       
[500]3001        for(i=0;i<nord;i++) {
3002                rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
3003                AR=(1.0/(rad*sigma*sqrt(2.0*Pi)))*exp(-(0.5*((log(radius/rad))/sigma)*((log(radius/rad))/sigma)));
3004                yyy = AR* Gauss20Wt[i] * BicelleIntegration(q,rad,radthick,facthick,rhoc,rhoh,rhor,rhosolv,length);
3005                Rsqryyy= Gauss20Wt[i] * AR * (rad+radthick)*(rad+radthick);             //SRK normalize to total dimensions
3006                summ += yyy;
3007                Rsqrsumm += Rsqryyy;
3008        }
3009       
3010        answer = (uplim-lolim)/2.0*summ;
3011        Rsqr = (uplim-lolim)/2.0*Rsqrsumm;
3012        //normalize by average cylinder volume
3013        Vpoly = Pi*Rsqr*(length+2*facthick);
3014        answer /= Vpoly;
3015        //convert to [cm-1]
3016        answer *= 1.0e8;
3017        //Scale
3018        answer *= scale;
3019        // add in the background
3020        answer += bkg;
[501]3021               
3022        return answer;
[500]3023       
3024}
3025
3026double
3027BicelleIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length){
3028
3029        double answer,halfheight,Pi;
3030        double lolim,uplim,summ,yyy,zi;
3031        int nord,i;
3032       
3033        // set up the integration end points
3034        Pi = 4.0*atan(1.0);
3035        nord = 76;
[632]3036        lolim = 0.0;
[500]3037        uplim = Pi/2;
3038        halfheight = length/2.0;
3039       
3040        summ = 0.0;                             // initialize integral
3041        i=0;
3042        for(i=0;i<nord;i++) {
3043                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
3044                yyy = Gauss76Wt[i] * BicelleKernel(qq, rad, radthick, facthick, rhoc, rhoh, rhor,rhosolv, halfheight, zi);
3045                summ += yyy;
3046        }
3047       
3048        // calculate value of integral to return
3049        answer = (uplim-lolim)/2.0*summ;
3050        return(answer); 
3051}
3052
3053double
[594]3054BicelleKernel(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length, double dum)
3055{
[500]3056        double dr1,dr2,dr3;
3057        double besarg1,besarg2;
3058        double vol1,vol2,vol3;
3059        double sinarg1,sinarg2;
3060        double t1,t2,t3;
3061        double retval;
3062       
3063        double Pi = 4.0*atan(1.0);
3064       
3065        dr1 = rhoc-rhoh;
3066        dr2 = rhor-rhosolv;
3067        dr3=  rhoh-rhor;
[632]3068        vol1 = Pi*rad*rad*(2.0*length);
3069        vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick);
3070        vol3= Pi*(rad)*(rad)*(2.0*length+2.0*facthick);
[500]3071        besarg1 = qq*rad*sin(dum);
3072        besarg2 = qq*(rad+radthick)*sin(dum);
3073        sinarg1 = qq*length*cos(dum);
3074        sinarg2 = qq*(length+facthick)*cos(dum);
3075       
[632]3076        t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1;
3077        t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2;
3078        t3 = 2.0*vol3*dr3*sin(sinarg2)/sinarg2*NR_BessJ1(besarg1)/besarg1;
[500]3079       
3080        retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum);
[501]3081        return(retval);
[500]3082       
3083}
[501]3084
3085
3086double
[594]3087CSPPKernel(double dp[], double mu, double uu)
3088{       
[501]3089        double aa,bb,cc, ta,tb,tc; 
3090        double Vin,Vot,V1,V2;
3091        double rhoA,rhoB,rhoC, rhoP, rhosolv;
[594]3092        double dr0, drA,drB, drC;
[501]3093        double arg1,arg2,arg3,arg4,t1,t2, t3, t4;
[594]3094        double Pi,retVal;
[501]3095
3096        aa = dp[1];
3097        bb = dp[2];
3098        cc = dp[3];
3099        ta = dp[4];
3100        tb = dp[5];
3101        tc = dp[6];
3102        rhoA=dp[7];
3103        rhoB=dp[8];
3104        rhoC=dp[9];
3105        rhoP=dp[10];
3106        rhosolv=dp[11];
3107        dr0=rhoP-rhosolv;
3108        drA=rhoA-rhosolv;
3109        drB=rhoB-rhosolv;
3110        drC=rhoC-rhosolv; 
3111        Vin=(aa*bb*cc);
[594]3112        Vot=(aa*bb*cc+2.0*ta*bb*cc+2.0*aa*tb*cc+2.0*aa*bb*tc);
3113        V1=(2.0*ta*bb*cc);   //  incorrect V1 (aa*bb*cc+2*ta*bb*cc)
3114        V2=(2.0*aa*tb*cc);  // incorrect V2(aa*bb*cc+2*aa*tb*cc)
[501]3115        aa = aa/bb;
[594]3116        ta=(aa+2.0*ta)/bb;
3117        tb=(aa+2.0*tb)/bb;
[501]3118       
3119        Pi = 4.0*atan(1.0);
3120       
[594]3121        arg1 = (mu*aa/2.0)*sin(Pi*uu/2.0);
3122        arg2 = (mu/2.0)*cos(Pi*uu/2.0);
3123        arg3=  (mu*ta/2.0)*sin(Pi*uu/2.0);
3124        arg4=  (mu*tb/2.0)*cos(Pi*uu/2.0);
[501]3125                         
[632]3126        if(arg1==0.0){
[594]3127                t1 = 1.0;
[501]3128        } else {
3129                t1 = (sin(arg1)/arg1);                //defn for CSPP model sin(arg1)/arg1    test:  (sin(arg1)/arg1)*(sin(arg1)/arg1)   
3130        }
[632]3131        if(arg2==0.0){
[594]3132                t2 = 1.0;
[501]3133        } else {
3134                t2 = (sin(arg2)/arg2);           //defn for CSPP model sin(arg2)/arg2   test: (sin(arg2)/arg2)*(sin(arg2)/arg2)   
3135        }       
[632]3136        if(arg3==0.0){
[594]3137                t3 = 1.0;
[501]3138        } else {
3139                t3 = sin(arg3)/arg3;
3140        }
[632]3141        if(arg4==0.0){
[594]3142                t4 = 1.0;
[501]3143        } else {
3144                t4 = sin(arg4)/arg4;
3145        }
[594]3146        retVal =( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 )*( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 );   //  correct FF : square of sum of phase factors
[501]3147        return(retVal); 
[594]3148
3149}
3150
3151/*      CSParallelepiped  :  calculates the form factor of a Parallelepiped with a core-shell structure
3152 -- different SLDs can be used for the face and rim
3153
3154Uses 76 pt Gaussian quadrature for both integrals
3155*/
3156double
3157CSParallelepiped(double dp[], double q)
3158{
3159        int i,j;
3160        double scale,aa,bb,cc,ta,tb,tc,rhoA,rhoB,rhoC,rhoP,rhosolv,bkg;         //local variables of coefficient wave
3161        int nordi=76;                   //order of integration
3162        int nordj=76;
3163        double va,vb;           //upper and lower integration limits
3164        double summ,yyy,answer;                 //running tally of integration
3165        double summj,vaj,vbj;                   //for the inner integration
3166        double mu,mudum,arg,sigma,uu,vol;
[501]3167       
[594]3168       
3169        //      Pi = 4.0*atan(1.0);
3170        va = 0.0;
3171        vb = 1.0;               //orintational average, outer integral
3172        vaj = 0.0;
3173        vbj = 1.0;              //endpoints of inner integral
3174       
3175        summ = 0.0;                     //initialize intergral
3176       
3177        scale = dp[0];
3178        aa = dp[1];
3179        bb = dp[2];
3180        cc = dp[3];
3181        ta  = dp[4];
3182        tb  = dp[5];
3183        tc  = dp[6];   // is 0 at the moment 
3184        rhoA=dp[7];   //rim A SLD
3185        rhoB=dp[8];   //rim B SLD
3186        rhoC=dp[9];    //rim C SLD
3187        rhoP = dp[10];   //Parallelpiped core SLD
3188        rhosolv=dp[11];  // Solvent SLD
3189        bkg = dp[12];
3190       
3191        mu = q*bb;
3192        vol = aa*bb*cc+2.0*ta*bb*cc+2.0*aa*tb*cc+2.0*aa*bb*tc;          //calculate volume before rescaling
3193       
3194        // do the rescaling here, not in the kernel
3195        // normalize all WRT bb
3196        aa = aa/bb;
3197        cc = cc/bb;
3198       
3199        for(i=0;i<nordi;i++) {
3200                //setup inner integral over the ellipsoidal cross-section
[632]3201                summj=0.0;
[594]3202                sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;          //the outer dummy
3203               
3204                for(j=0;j<nordj;j++) {
3205                        //76 gauss points for the inner integral
3206                        uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;         //the inner dummy
3207                        mudum = mu*sqrt(1.0-sigma*sigma);
3208                        yyy = Gauss76Wt[j] * CSPPKernel(dp,mudum,uu);
3209                        summj += yyy;
3210                }
3211                //now calculate the value of the inner integral
3212                answer = (vbj-vaj)/2.0*summj;
3213               
3214                //finish the outer integral cc already scaled
3215                arg = mu*cc*sigma/2.0;
[632]3216                if ( arg == 0.0 ) {
[594]3217                        answer *= 1.0;
3218                } else {
3219                        answer *= sin(arg)*sin(arg)/arg/arg;
3220                }
3221               
3222                //now sum up the outer integral
3223                yyy = Gauss76Wt[i] * answer;
3224                summ += yyy;
3225        }               //final scaling is done at the end of the function, after the NT_FP64 case
3226       
3227        answer = (vb-va)/2.0*summ;
[501]3228
[594]3229        //normalize by volume
3230        answer /= vol;
3231        //convert to [cm-1]
3232        answer *= 1.0e8;
3233        //Scale
3234        answer *= scale;
3235        // add in the background
3236        answer += bkg;
3237       
3238        return answer;
3239}
[501]3240
Note: See TracBrowser for help on using the repository browser.