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

Last change on this file since 594 was 594, checked in by srkline, 13 years ago

Correcting my merge mistakes in CSParallelepiped

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