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

Last change on this file since 834 was 834, checked in by srkline, 11 years ago

Changes to the XOP code to upgrade to ToolKit? v6. Changes are the ones outlined in the Appendix A of the TK6 manual. SOme of the XOP support routines and the #pragma for 2-byte structures have changed. Per Howard Rodstein, there is no need to change the c files to cpp. c should work and compile just fine.

These changes work correctly on my mac. Next is to make sure that they work correctly on the two build machines.

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