source: sans/Analysis/branches/ajj_23APR07/XOPs/SANSAnalysis/lib/libCylinder.c @ 97

Last change on this file since 97 was 97, checked in by ajj, 15 years ago

Now committing the code correctly - having relied on XCode's SVN interface which doesn't behave (quelle surprise).

I hope this isn't too screwed up.

  • Property svn:executable set to *
File size: 55.2 KB
Line 
1/*      CylinderFit.c
2
3A simplified project designed to act as a template for your curve fitting function.
4The fitting function is a Cylinder form factor. No resolution effects are included (yet)
5*/
6
7#include "StandardHeaders.h"                    // Include ANSI headers, Mac headers
8#include "GaussWeights.h"
9#include "libCylinder.h"
10
11/*      CylinderForm  :  calculates the form factor of a cylinder at the give x-value p->x
12
13Warning:
14The call to WaveData() below returns a pointer to the middle
15of an unlocked Macintosh handle. In the unlikely event that your
16calculations could cause memory to move, you should copy the coefficient
17values to local variables or an array before such operations.
18*/
19double
20CylinderForm(double dp[], double q)
21{
22        int i;
23        double Pi;
24        double scale,radius,length,delrho,bkg,halfheight;               //local variables of coefficient wave
25        int nord=76;                    //order of integration
26        double uplim,lolim;             //upper and lower integration limits
27        double summ,zi,yyy,answer,vcyl;                 //running tally of integration
28       
29        Pi = 4.0*atan(1.0);
30        lolim = 0;
31        uplim = Pi/2.0;
32       
33        summ = 0.0;                     //initialize intergral
34       
35        scale = dp[0];                  //make local copies in case memory moves
36        radius = dp[1];
37        length = dp[2];
38        delrho = dp[3];
39        bkg = dp[4];
40        halfheight = length/2.0;
41        for(i=0;i<nord;i++) {
42                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
43                yyy = Gauss76Wt[i] * CylKernel(q, radius, halfheight, zi);
44                summ += yyy;
45        }
46       
47        answer = (uplim-lolim)/2.0*summ;
48        // Multiply by contrast^2
49        answer *= delrho*delrho;
50        //normalize by cylinder volume
51        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
52        vcyl=Pi*radius*radius*length;
53        answer *= vcyl;
54        //convert to [cm-1]
55        answer *= 1.0e8;
56        //Scale
57        answer *= scale;
58        // add in the background
59        answer += bkg;
60       
61        return answer;
62}
63
64/*      EllipCyl76X  :  calculates the form factor of a elliptical cylinder at the given x-value p->x
65
66Uses 76 pt Gaussian quadrature for both integrals
67
68Warning:
69The call to WaveData() below returns a pointer to the middle
70of an unlocked Macintosh handle. In the unlikely event that your
71calculations could cause memory to move, you should copy the coefficient
72values to local variables or an array before such operations.
73*/
74double
75EllipCyl76(double dp[], double q)
76{
77        int i,j;
78        double Pi;
79        double scale,ra,nu,length,delrho,bkg;           //local variables of coefficient wave
80        int nord=76;                    //order of integration
81        double va,vb;           //upper and lower integration limits
82        double summ,zi,yyy,answer,vell;                 //running tally of integration
83        double summj,vaj,vbj,zij,arg;                   //for the inner integration
84       
85        Pi = 4.0*atan(1.0);
86        va = 0;
87        vb = 1;         //orintational average, outer integral
88        vaj=0;
89        vbj=Pi;         //endpoints of inner integral
90       
91        summ = 0.0;                     //initialize intergral
92       
93        scale = dp[0];                  //make local copies in case memory moves
94        ra = dp[1];
95        nu = dp[2];
96        length = dp[3];
97        delrho = dp[4];
98        bkg = dp[5];
99        for(i=0;i<nord;i++) {
100                //setup inner integral over the ellipsoidal cross-section
101                summj=0;
102                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy
103                arg = ra*sqrt(1-zi*zi);
104                for(j=0;j<nord;j++) {
105                        //76 gauss points for the inner integral as well
106                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "y" dummy
107                        yyy = Gauss76Wt[j] * EllipCylKernel(q,arg,nu,zij);
108                        summj += yyy;
109                }
110                //now calculate the value of the inner integral
111                answer = (vbj-vaj)/2.0*summj;
112                //divide integral by Pi
113                answer /=Pi;
114               
115                //now calculate outer integral
116                arg = q*length*zi/2;
117                yyy = Gauss76Wt[i] * answer * sin(arg) * sin(arg) / arg / arg;
118                summ += yyy;
119        }
120        answer = (vb-va)/2.0*summ;
121        // Multiply by contrast^2
122        answer *= delrho*delrho;
123        //normalize by cylinder volume
124        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
125        vell = Pi*ra*(nu*ra)*length;
126        answer *= vell;
127        //convert to [cm-1]
128        answer *= 1.0e8;
129        //Scale
130        answer *= scale;
131        // add in the background
132        answer += bkg;
133       
134        return answer;
135}
136
137/*      EllipCyl20X  :  calculates the form factor of a elliptical cylinder at the given x-value p->x
138
139Uses 76 pt Gaussian quadrature for orientational integral
140Uses 20 pt quadrature for the inner integral over the elliptical cross-section
141
142Warning:
143The call to WaveData() below returns a pointer to the middle
144of an unlocked Macintosh handle. In the unlikely event that your
145calculations could cause memory to move, you should copy the coefficient
146values to local variables or an array before such operations.
147*/
148double
149EllipCyl20(double dp[], double q)
150{
151        int i,j;
152        double Pi;
153        double scale,ra,nu,length,delrho,bkg;           //local variables of coefficient wave
154        int nordi=76;                   //order of integration
155        int nordj=20;
156        double va,vb;           //upper and lower integration limits
157        double summ,zi,yyy,answer,vell;                 //running tally of integration
158        double summj,vaj,vbj,zij,arg;                   //for the inner integration
159       
160        Pi = 4.0*atan(1.0);
161        va = 0;
162        vb = 1;         //orintational average, outer integral
163        vaj=0;
164        vbj=Pi;         //endpoints of inner integral
165       
166        summ = 0.0;                     //initialize intergral
167       
168        scale = dp[0];                  //make local copies in case memory moves
169        ra = dp[1];
170        nu = dp[2];
171        length = dp[3];
172        delrho = dp[4];
173        bkg = dp[5];
174        for(i=0;i<nordi;i++) {
175                //setup inner integral over the ellipsoidal cross-section
176                summj=0;
177                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy
178                arg = ra*sqrt(1-zi*zi);
179                for(j=0;j<nordj;j++) {
180                        //20 gauss points for the inner integral
181                        zij = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "y" dummy
182                        yyy = Gauss20Wt[j] * EllipCylKernel(q,arg,nu,zij);
183                        summj += yyy;
184                }
185                //now calculate the value of the inner integral
186                answer = (vbj-vaj)/2.0*summj;
187                //divide integral by Pi
188                answer /=Pi;
189               
190                //now calculate outer integral
191                arg = q*length*zi/2;
192                yyy = Gauss76Wt[i] * answer * sin(arg) * sin(arg) / arg / arg;
193                summ += yyy;
194        }
195       
196        answer = (vb-va)/2.0*summ;
197        // Multiply by contrast^2
198        answer *= delrho*delrho;
199        //normalize by cylinder volume
200        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
201        vell = Pi*ra*(nu*ra)*length;
202        answer *= vell;
203        //convert to [cm-1]
204        answer *= 1.0e8;
205        //Scale
206        answer *= scale;
207        // add in the background
208        answer += bkg; 
209       
210        return answer;
211}
212
213/*      TriaxialEllipsoidX  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x
214
215Uses 76 pt Gaussian quadrature for both integrals
216
217Warning:
218The call to WaveData() below returns a pointer to the middle
219of an unlocked Macintosh handle. In the unlikely event that your
220calculations could cause memory to move, you should copy the coefficient
221values to local variables or an array before such operations.
222*/
223double
224TriaxialEllipsoid(double dp[], double q)
225{
226        int i,j;
227        double Pi;
228        double scale,aa,bb,cc,delrho,bkg;               //local variables of coefficient wave
229        int nordi=76;                   //order of integration
230        int nordj=76;
231        double va,vb;           //upper and lower integration limits
232        double summ,zi,yyy,answer;                      //running tally of integration
233        double summj,vaj,vbj,zij;                       //for the inner integration
234       
235        Pi = 4.0*atan(1.0);
236        va = 0;
237        vb = 1;         //orintational average, outer integral
238        vaj = 0;
239        vbj = 1;                //endpoints of inner integral
240       
241        summ = 0.0;                     //initialize intergral
242       
243        scale = dp[0];                  //make local copies in case memory moves
244        aa = dp[1];
245        bb = dp[2];
246        cc = dp[3];
247        delrho = dp[4];
248        bkg = dp[5];
249        for(i=0;i<nordi;i++) {
250                //setup inner integral over the ellipsoidal cross-section
251                summj=0;
252                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy
253                for(j=0;j<nordj;j++) {
254                        //20 gauss points for the inner integral
255                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "y" dummy
256                        yyy = Gauss76Wt[j] * TriaxialKernel(q,aa,bb,cc,zi,zij);
257                        summj += yyy;
258                }
259                //now calculate the value of the inner integral
260                answer = (vbj-vaj)/2.0*summj;
261               
262                //now calculate outer integral
263                yyy = Gauss76Wt[i] * answer;
264                summ += yyy;
265        }               //final scaling is done at the end of the function, after the NT_FP64 case
266       
267        answer = (vb-va)/2.0*summ;
268        // Multiply by contrast^2
269        answer *= delrho*delrho;
270        //normalize by ellipsoid volume
271        answer *= 4*Pi/3*aa*bb*cc;
272        //convert to [cm-1]
273        answer *= 1.0e8;
274        //Scale
275        answer *= scale;
276        // add in the background
277        answer += bkg;
278       
279        return answer;
280}
281
282/*      ParallelepipedX  :  calculates the form factor of a Parallelepiped (a rectangular solid)
283at the given x-value p->x
284
285Uses 76 pt Gaussian quadrature for both integrals
286
287Warning:
288The call to WaveData() below returns a pointer to the middle
289of an unlocked Macintosh handle. In the unlikely event that your
290calculations could cause memory to move, you should copy the coefficient
291values to local variables or an array before such operations.
292*/
293double
294Parallelepiped(double dp[], double q)
295{
296        int i,j;
297        double scale,aa,bb,cc,delrho,bkg;               //local variables of coefficient wave
298        int nordi=76;                   //order of integration
299        int nordj=76;
300        double va,vb;           //upper and lower integration limits
301        double summ,yyy,answer;                 //running tally of integration
302        double summj,vaj,vbj;                   //for the inner integration
303        double mu,mudum,arg,sigma,uu,vol;
304       
305       
306        //      Pi = 4.0*atan(1.0);
307        va = 0;
308        vb = 1;         //orintational average, outer integral
309        vaj = 0;
310        vbj = 1;                //endpoints of inner integral
311       
312        summ = 0.0;                     //initialize intergral
313       
314        scale = dp[0];                  //make local copies in case memory moves
315        aa = dp[1];
316        bb = dp[2];
317        cc = dp[3];
318        delrho = dp[4];
319        bkg = dp[5];
320       
321        mu = q*bb;
322        vol = aa*bb*cc;
323        // normalize all WRT bb
324        aa = aa/bb;
325        cc = cc/bb;
326       
327        for(i=0;i<nordi;i++) {
328                //setup inner integral over the ellipsoidal cross-section
329                summj=0;
330                sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;          //the outer dummy
331               
332                for(j=0;j<nordj;j++) {
333                        //76 gauss points for the inner integral
334                        uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;         //the inner dummy
335                        mudum = mu*sqrt(1-sigma*sigma);
336                        yyy = Gauss76Wt[j] * PPKernel(aa,mudum,uu);
337                        summj += yyy;
338                }
339                //now calculate the value of the inner integral
340                answer = (vbj-vaj)/2.0*summj;
341               
342                arg = mu*cc*sigma/2;
343                if ( arg == 0 ) {
344                        answer *= 1;
345                } else {
346                        answer *= sin(arg)*sin(arg)/arg/arg;
347                }
348               
349                //now sum up the outer integral
350                yyy = Gauss76Wt[i] * answer;
351                summ += yyy;
352        }               //final scaling is done at the end of the function, after the NT_FP64 case
353       
354        answer = (vb-va)/2.0*summ;
355        // Multiply by contrast^2
356        answer *= delrho*delrho;
357        //normalize by volume
358        answer *= vol;
359        //convert to [cm-1]
360        answer *= 1.0e8;
361        //Scale
362        answer *= scale;
363        // add in the background
364        answer += bkg;
365       
366        return answer;
367}
368
369/*      HollowCylinderX  :  calculates the form factor of a Hollow Cylinder
370at the given x-value p->x
371
372Uses 76 pt Gaussian quadrature for the single integral
373
374Warning:
375The call to WaveData() below returns a pointer to the middle
376of an unlocked Macintosh handle. In the unlikely event that your
377calculations could cause memory to move, you should copy the coefficient
378values to local variables or an array before such operations.
379*/
380double
381HollowCylinder(double dp[], double q)
382{
383        int i;
384        double scale,rcore,rshell,length,delrho,bkg;            //local variables of coefficient wave
385        int nord=76;                    //order of integration
386        double va,vb,zi;                //upper and lower integration limits
387        double summ,answer,pi;                  //running tally of integration
388       
389        pi = 4.0*atan(1.0);
390        va = 0;
391        vb = 1;         //limits of numerical integral
392       
393        summ = 0.0;                     //initialize intergral
394       
395        scale = dp[0];          //make local copies in case memory moves
396        rcore = dp[1];
397        rshell = dp[2];
398        length = dp[3];
399        delrho = dp[4];
400        bkg = dp[5];
401       
402        for(i=0;i<nord;i++) {
403                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
404                summ += Gauss76Wt[i] * HolCylKernel(q, rcore, rshell, length, zi);
405        }
406       
407        answer = (vb-va)/2.0*summ;
408        // Multiply by contrast^2
409        answer *= delrho*delrho;
410        //normalize by volume
411        answer *= pi*(rshell*rshell-rcore*rcore)*length;
412        //convert to [cm-1]
413        answer *= 1.0e8;
414        //Scale
415        answer *= scale;
416        // add in the background
417        answer += bkg;
418       
419        return answer;
420}
421
422/*      EllipsoidFormX  :  calculates the form factor of an ellipsoid of revolution with semiaxes a:a:nua
423at the given x-value p->x
424
425Uses 76 pt Gaussian quadrature for the single integral
426
427Warning:
428The call to WaveData() below returns a pointer to the middle
429of an unlocked Macintosh handle. In the unlikely event that your
430calculations could cause memory to move, you should copy the coefficient
431values to local variables or an array before such operations.
432*/
433double
434EllipsoidForm(double dp[], double q)
435{
436        int i;
437        double scale,a,nua,delrho,bkg;          //local variables of coefficient wave
438        int nord=76;                    //order of integration
439        double va,vb,zi;                //upper and lower integration limits
440        double summ,answer,pi;                  //running tally of integration
441       
442        pi = 4.0*atan(1.0);
443        va = 0;
444        vb = 1;         //limits of numerical integral
445       
446        summ = 0.0;                     //initialize intergral
447       
448        scale = dp[0];                  //make local copies in case memory moves
449        nua = dp[1];
450        a = dp[2];
451        delrho = dp[3];
452        bkg = dp[4];
453       
454        for(i=0;i<nord;i++) {
455                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
456                summ += Gauss76Wt[i] * EllipsoidKernel(q, a, nua, zi);
457        }
458       
459        answer = (vb-va)/2.0*summ;
460        // Multiply by contrast^2
461        answer *= delrho*delrho;
462        //normalize by volume
463        answer *= 4*pi/3*a*a*nua;
464        //convert to [cm-1]
465        answer *= 1.0e8;
466        //Scale
467        answer *= scale;
468        // add in the background
469        answer += bkg;
470       
471        return answer;
472}
473
474
475/*      Cyl_PolyRadiusX  :  calculates the form factor of a cylinder at the given x-value p->x
476the cylinder has a polydisperse cross section
477
478*/
479double
480Cyl_PolyRadius(double dp[], double q)
481{
482        int i;
483        double scale,radius,length,pd,delrho,bkg;               //local variables of coefficient wave
484        int nord=20;                    //order of integration
485        double uplim,lolim;             //upper and lower integration limits
486        double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
487        double range,zz,Pi;
488       
489        Pi = 4.0*atan(1.0);
490        range = 3.4;
491       
492        summ = 0.0;                     //initialize intergral
493       
494        scale = dp[0];                  //make local copies in case memory moves
495        radius = dp[1];
496        length = dp[2];
497        pd = dp[3];
498        delrho = dp[4];
499        bkg = dp[5];
500       
501        zz = (1.0/pd)*(1.0/pd) - 1.0;
502       
503        lolim = radius*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
504        if(lolim<0) {
505                lolim = 0;
506        }
507        if(pd>0.3) {
508                range = 3.4 + (pd-0.3)*18.0;
509        }
510        uplim = radius*(1.0+range*pd);
511       
512        for(i=0;i<nord;i++) {
513                zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
514                yyy = Gauss20Wt[i] * Cyl_PolyRadKernel(q, radius, length, zz, delrho, zi);
515                summ += yyy;
516        }
517       
518        answer = (uplim-lolim)/2.0*summ;
519        //normalize by average cylinder volume
520        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
521        Vpoly=Pi*radius*radius*length*(zz+2.0)/(zz+1.0);
522        answer /= Vpoly;
523        //convert to [cm-1]
524        answer *= 1.0e8;
525        //Scale
526        answer *= scale;
527        // add in the background
528        answer += bkg;
529       
530        return answer;
531}
532
533/*      Cyl_PolyLengthX  :  calculates the form factor of a cylinder at the given x-value p->x
534the cylinder has a polydisperse Length
535
536*/
537double
538Cyl_PolyLength(double dp[], double q)
539{
540        int i;
541        double scale,radius,length,pd,delrho,bkg;               //local variables of coefficient wave
542        int nord=20;                    //order of integration
543        double uplim,lolim;             //upper and lower integration limits
544        double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
545        double range,zz,Pi;
546       
547       
548        Pi = 4.0*atan(1.0);
549        range = 3.4;
550       
551        summ = 0.0;                     //initialize intergral
552       
553        scale = dp[0];                  //make local copies in case memory moves
554        radius = dp[1];
555        length = dp[2];
556        pd = dp[3];
557        delrho = dp[4];
558        bkg = dp[5];
559       
560        zz = (1.0/pd)*(1.0/pd) - 1.0;
561       
562        lolim = length*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
563        if(lolim<0) {
564                lolim = 0;
565        }
566        if(pd>0.3) {
567                range = 3.4 + (pd-0.3)*18.0;
568        }
569        uplim = length*(1.0+range*pd);
570       
571        for(i=0;i<nord;i++) {
572                zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
573                yyy = Gauss20Wt[i] * Cyl_PolyLenKernel(q, radius, length, zz, delrho, zi);
574                summ += yyy;
575        }
576       
577        answer = (uplim-lolim)/2.0*summ;
578        //normalize by average cylinder volume (first moment)
579        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
580        Vpoly=Pi*radius*radius*length;
581        answer /= Vpoly;
582        //convert to [cm-1]
583        answer *= 1.0e8;
584        //Scale
585        answer *= scale;
586        // add in the background
587        answer += bkg;
588       
589        return answer;
590}
591
592/*      CoreShellCylinderX  :  calculates the form factor of a cylinder at the given x-value p->x
593the cylinder has a core-shell structure
594
595*/
596double
597CoreShellCylinder(double dp[], double q)
598{
599        int i;
600        double scale,rcore,length,bkg;          //local variables of coefficient wave
601        double thick,rhoc,rhos,rhosolv;
602        int nord=76;                    //order of integration
603        double uplim,lolim,halfheight;          //upper and lower integration limits
604        double summ,zi,yyy,answer,Vcyl;                 //running tally of integration
605        double Pi;
606       
607        Pi = 4.0*atan(1.0);
608       
609        lolim = 0.0;
610        uplim = Pi/2.0;
611       
612        summ = 0.0;                     //initialize intergral
613       
614        scale = dp[0];          //make local copies in case memory moves
615        rcore = dp[1];
616        thick = dp[2];
617        length = dp[3];
618        rhoc = dp[4];
619        rhos = dp[5];
620        rhosolv = dp[6];
621        bkg = dp[7];
622       
623        halfheight = length/2.0;
624       
625        for(i=0;i<nord;i++) {
626                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
627                yyy = Gauss76Wt[i] * CoreShellCylKernel(q, rcore, thick, rhoc,rhos,rhosolv, halfheight, zi);
628                summ += yyy;
629        }
630       
631        answer = (uplim-lolim)/2.0*summ;
632        // length is the total core length
633        Vcyl=Pi*(rcore+thick)*(rcore+thick)*(length+2.0*thick);
634        answer /= Vcyl;
635        //convert to [cm-1]
636        answer *= 1.0e8;
637        //Scale
638        answer *= scale;
639        // add in the background
640        answer += bkg;
641       
642        return answer;
643}
644
645
646/*      PolyCoShCylinderX  :  calculates the form factor of a core-shell cylinder at the given x-value p->x
647the cylinder has a polydisperse CORE radius
648
649*/
650double
651PolyCoShCylinder(double dp[], double q)
652{
653        int i;
654        double scale,radius,length,sigma,bkg;           //local variables of coefficient wave
655        double rad,radthick,facthick,rhoc,rhos,rhosolv;
656        int nord=20;                    //order of integration
657        double uplim,lolim;             //upper and lower integration limits
658        double summ,yyy,answer,Vpoly;                   //running tally of integration
659        double Pi,AR,Rsqrsumm,Rsqryyy,Rsqr;
660       
661        Pi = 4.0*atan(1.0);
662       
663        summ = 0.0;                     //initialize intergral
664        Rsqrsumm = 0.0;
665       
666        scale = dp[0];
667        radius = dp[1];
668        sigma = dp[2];                          //sigma is the standard mean deviation
669        length = dp[3];
670        radthick = dp[4];
671        facthick= dp[5];
672        rhoc = dp[6];
673        rhos = dp[7];
674        rhosolv = dp[8];
675        bkg = dp[9];
676       
677        lolim = exp(log(radius)-(4.*sigma));
678        if (lolim<0) {
679                lolim=0;                //to avoid numerical error when  va<0 (-ve r value)
680        }
681        uplim = exp(log(radius)+(4.*sigma));
682       
683        for(i=0;i<nord;i++) {
684                rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
685                AR=(1.0/(rad*sigma*sqrt(2.0*Pi)))*exp(-(0.5*((log(radius/rad))/sigma)*((log(radius/rad))/sigma)));
686                yyy = AR* Gauss20Wt[i] * CSCylIntegration(q,rad,radthick,facthick,rhoc,rhos,rhosolv,length);
687                Rsqryyy= Gauss20Wt[i] * AR * (rad+radthick)*(rad+radthick);             //SRK normalize to total dimensions
688                summ += yyy;
689                Rsqrsumm += Rsqryyy;
690        }
691       
692        answer = (uplim-lolim)/2.0*summ;
693        Rsqr = (uplim-lolim)/2.0*Rsqrsumm;
694        //normalize by average cylinder volume
695        Vpoly = Pi*Rsqr*(length+2*facthick);
696        answer /= Vpoly;
697        //convert to [cm-1]
698        answer *= 1.0e8;
699        //Scale
700        answer *= scale;
701        // add in the background
702        answer += bkg;
703       
704        return answer;
705}
706
707/*      OblateFormX  :  calculates the form factor of a core-shell Oblate ellipsoid at the given x-value p->x
708the ellipsoid has a core-shell structure
709
710*/
711double
712OblateForm(double dp[], double q)
713{
714        int i;
715        double scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg;
716        int nord=76;                    //order of integration
717        double uplim,lolim;             //upper and lower integration limits
718        double summ,zi,yyy,answer,oblatevol;                    //running tally of integration
719        double Pi;
720       
721        Pi = 4.0*atan(1.0);
722       
723        lolim = 0.0;
724        uplim = 1.0;
725       
726        summ = 0.0;                     //initialize intergral
727       
728       
729        scale = dp[0];          //make local copies in case memory moves
730        crmaj = dp[1];
731        crmin = dp[2];
732        trmaj = dp[3];
733        trmin = dp[4];
734        delpc = dp[5];
735        delps = dp[6];
736        bkg = dp[7];                       
737       
738        for(i=0;i<nord;i++) {
739                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
740                yyy = Gauss76Wt[i] * gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,q);
741                summ += yyy;
742        }
743       
744        answer = (uplim-lolim)/2.0*summ;
745        // normalize by particle volume
746        oblatevol = 4*Pi/3*trmaj*trmaj*trmin;
747        answer /= oblatevol;
748       
749        //convert to [cm-1]
750        answer *= 1.0e8;
751        //Scale
752        answer *= scale;
753        // add in the background
754        answer += bkg;
755       
756        return answer;
757}
758
759/*      ProlateFormX  :  calculates the form factor of a core-shell Prolate ellipsoid at the given x-value p->x
760the ellipsoid has a core-shell structure
761
762*/
763double
764ProlateForm(double dp[], double q)
765{
766        int i;
767        double scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg;
768        int nord=76;                    //order of integration
769        double uplim,lolim;             //upper and lower integration limits
770        double summ,zi,yyy,answer,prolatevol;                   //running tally of integration
771        double Pi;
772       
773        Pi = 4.0*atan(1.0);
774       
775        lolim = 0.0;
776        uplim = 1.0;
777       
778        summ = 0.0;                     //initialize intergral
779       
780        scale = dp[0];          //make local copies in case memory moves
781        crmaj = dp[1];
782        crmin = dp[2];
783        trmaj = dp[3];
784        trmin = dp[4];
785        delpc = dp[5];
786        delps = dp[6];
787        bkg = dp[7];                       
788       
789        for(i=0;i<nord;i++) {
790                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
791                yyy = Gauss76Wt[i] * gfn2(zi,crmaj,crmin,trmaj,trmin,delpc,delps,q);
792                summ += yyy;
793        }
794       
795        answer = (uplim-lolim)/2.0*summ;
796        // normalize by particle volume
797        prolatevol = 4*Pi/3*trmaj*trmin*trmin;
798        answer /= prolatevol;
799       
800        //convert to [cm-1]
801        answer *= 1.0e8;
802        //Scale
803        answer *= scale;
804        // add in the background
805        answer += bkg;
806       
807        return answer;
808}
809
810
811/*      StackedDiscsX  :  calculates the form factor of a stacked "tactoid" of core shell disks
812like clay platelets that are not exfoliated
813
814*/
815double
816StackedDiscs(double dp[], double q)
817{
818        int i;
819        double scale,length,bkg,rcore,thick,rhoc,rhol,rhosolv,N,gsd;            //local variables of coefficient wave
820        double va,vb,vcyl,summ,yyy,zi,halfheight,d,answer;
821        int nord=76;                    //order of integration
822        double Pi;
823       
824       
825        Pi = 4.0*atan(1.0);
826       
827        va = 0.0;
828        vb = Pi/2.0;
829       
830        summ = 0.0;                     //initialize intergral
831       
832        scale = dp[0];
833        rcore = dp[1];
834        length = dp[2];
835        thick = dp[3];
836        rhoc = dp[4];
837        rhol = dp[5];
838        rhosolv = dp[6];
839        N = dp[7];
840        gsd = dp[8];
841        bkg = dp[9];
842       
843        d=2.0*thick+length;
844        halfheight = length/2.0;
845       
846        for(i=0;i<nord;i++) {
847                zi = ( Gauss76Z[i]*(vb-va) + vb + va )/2.0;
848                yyy = Gauss76Wt[i] * Stackdisc_kern(q, rcore, rhoc,rhol,rhosolv, halfheight,thick,zi,gsd,d,N);
849                summ += yyy;
850        }
851       
852        answer = (vb-va)/2.0*summ;
853        // length is the total core length
854        vcyl=Pi*rcore*rcore*(2.0*thick+length)*N;
855        answer /= vcyl;
856        //Convert to [cm-1]
857        answer *= 1.0e8;
858        //Scale
859        answer *= scale;
860        // add in the background
861        answer += bkg;
862       
863        return answer;
864}
865
866
867/*      LamellarFFX  :  calculates the form factor of a lamellar structure - no S(q) effects included
868
869*/
870double
871LamellarFF(double dp[], double q)
872{
873        double scale,del,sig,contr,bkg;         //local variables of coefficient wave
874        double inten, qval,Pq;
875        double Pi;
876       
877       
878        Pi = 4.0*atan(1.0);
879        scale = dp[0];
880        del = dp[1];
881        sig = dp[2]*del;
882        contr = dp[3];
883        bkg = dp[4];
884       
885       
886        Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig));
887       
888        inten = 2.0*Pi*scale*Pq/(qval*qval);            //this is now dimensionless...
889       
890        inten /= del;                   //normalize by the thickness (in A)
891       
892        inten *= 1.0e8;         // 1/A to 1/cm
893       
894        return(inten+bkg);
895}
896
897/*      LamellarPSX  :  calculates the form factor of a lamellar structure - with S(q) effects included
898-------
899------- resolution effects ARE included, but only a CONSTANT default value, not the real q-dependent resolution!!
900
901        */
902double
903LamellarPS(double dp[], double q)
904{
905        double scale,dd,del,sig,contr,NN,Cp,bkg;                //local variables of coefficient wave
906        double inten, qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ;
907        double Pi,Euler,dQDefault,fii;
908        int ii,NNint;
909       
910        Euler = 0.5772156649;           // Euler's constant
911        dQDefault = 0.0025;             //[=] 1/A, q-resolution, default value
912        dQ = dQDefault;
913       
914        Pi = 4.0*atan(1.0);
915        qval = q;
916       
917        scale = dp[0];
918        dd = dp[1];
919        del = dp[2];
920        sig = dp[3]*del;
921        contr = dp[4];
922        NN = trunc(dp[5]);              //be sure that NN is an integer
923        Cp = dp[6];
924        bkg = dp[7];
925       
926        Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig));
927       
928        NNint = (int)NN;                //cast to an integer for the loop
929        ii=0;
930        Sq = 0.0;
931        for(ii=1;ii<(NNint-1);ii+=1) {
932       
933                fii = (double)ii;               //do I really need to do this?
934               
935                temp = 0.0;
936                alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler);
937                t1 = 2.0*dQ*dQ*dd*dd*alpha;
938                t2 = 2.0*qval*qval*dd*dd*alpha;
939                t3 = dQ*dQ*dd*dd*ii*ii;
940               
941                temp = 1.0-ii/NN;
942                temp *= cos(dd*qval*ii/(1.0+t1));
943                temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
944                temp /= sqrt(1.0+t1);
945               
946                Sq += temp;
947        }
948       
949        Sq *= 2.0;
950        Sq += 1.0;
951       
952        inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval);
953       
954        inten *= 1.0e8;         // 1/A to 1/cm
955       
956    return(inten+bkg);
957}
958
959
960/*      LamellarPS_HGX  :  calculates the form factor of a lamellar structure - with S(q) effects included
961-------
962------- resolution effects ARE included, but only a CONSTANT default value, not the real q-dependent resolution!!
963
964        */
965double
966LamellarPS_HG(double dp[], double q)
967{
968        double scale,dd,delT,delH,SLD_T,SLD_H,SLD_S,NN,Cp,bkg;          //local variables of coefficient wave
969        double inten,qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ,drh,drt;
970        double Pi,Euler,dQDefault,fii;
971        int ii,NNint;
972       
973       
974        Euler = 0.5772156649;           // Euler's constant
975        dQDefault = 0.0025;             //[=] 1/A, q-resolution, default value
976        dQ = dQDefault;
977       
978        Pi = 4.0*atan(1.0);
979        qval= q;
980       
981        scale = dp[0];
982        dd = dp[1];
983        delT = dp[2];
984        delH = dp[3];
985        SLD_T = dp[4];
986        SLD_H = dp[5];
987        SLD_S = dp[6];
988        NN = trunc(dp[7]);              //be sure that NN is an integer
989        Cp = dp[8];
990        bkg = dp[9];
991       
992       
993        drh = SLD_H - SLD_S;
994        drt = SLD_T - SLD_S;    //correction 13FEB06 by L.Porcar
995       
996        Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT);
997        Pq *= Pq;
998        Pq *= 4.0/(qval*qval);
999       
1000        NNint = (int)NN;                //cast to an integer for the loop
1001        ii=0;
1002        Sq = 0.0;
1003        for(ii=1;ii<(NNint-1);ii+=1) {
1004       
1005                fii = (double)ii;               //do I really need to do this?
1006               
1007                temp = 0.0;
1008                alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler);
1009                t1 = 2.0*dQ*dQ*dd*dd*alpha;
1010                t2 = 2.0*qval*qval*dd*dd*alpha;
1011                t3 = dQ*dQ*dd*dd*ii*ii;
1012               
1013                temp = 1.0-ii/NN;
1014                temp *= cos(dd*qval*ii/(1.0+t1));
1015                temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
1016                temp /= sqrt(1.0+t1);
1017               
1018                Sq += temp;
1019        }
1020       
1021        Sq *= 2.0;
1022        Sq += 1.0;
1023       
1024        inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval);
1025       
1026        inten *= 1.0e8;         // 1/A to 1/cm
1027       
1028        return(inten+bkg);
1029}
1030
1031/*      LamellarFF_HGX  :  calculates the form factor of a lamellar structure - no S(q) effects included
1032but extra SLD for head groups is included
1033
1034*/
1035double
1036LamellarFF_HG(double dp[], double q)
1037{
1038        double scale,delT,delH,slds,sldh,sldt,bkg;              //local variables of coefficient wave
1039        double inten, qval,Pq,drh,drt;
1040        double Pi;
1041       
1042       
1043        Pi = 4.0*atan(1.0);
1044        qval= q;
1045        scale = dp[0];
1046        delT = dp[1];
1047        delH = dp[2];
1048        sldt = dp[3];
1049        sldh = dp[4];
1050        slds = dp[5];
1051        bkg = dp[6];
1052       
1053       
1054        drh = sldh - slds;
1055        drt = sldt - slds;              //correction 13FEB06 by L.Porcar
1056       
1057        Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT);
1058        Pq *= Pq;
1059        Pq *= 4.0/(qval*qval);
1060       
1061        inten = 2.0*Pi*scale*Pq/(qval*qval);            //dimensionless...
1062       
1063        inten /= 2.0*(delT+delH);                       //normalize by the bilayer thickness
1064       
1065        inten *= 1.0e8;         // 1/A to 1/cm
1066       
1067        return(inten+bkg);
1068}
1069
1070/*      FlexExclVolCylX  :  calculates the form factor of a flexible cylinder with a circular cross section
1071-- incorporates Wei-Ren Chen's fixes - 2006
1072
1073        */
1074double
1075FlexExclVolCyl(double dp[], double q)
1076{
1077        double scale,L,B,bkg,rad,qr,cont;
1078        double Pi,flex,crossSect,answer;
1079       
1080       
1081        Pi = 4.0*atan(1.0);
1082       
1083        scale = dp[0];          //make local copies in case memory moves
1084        L = dp[1];
1085        B = dp[2];
1086        rad = dp[3];
1087        cont = dp[4];
1088        bkg = dp[5];
1089       
1090   
1091        qr = q*rad;
1092       
1093        flex = Sk_WR(q,L,B);
1094       
1095        crossSect = (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
1096        flex *= crossSect;
1097        flex *= Pi*rad*rad*L;
1098        flex *= cont*cont;
1099        flex *= 1.0e8;
1100        answer = scale*flex + bkg;
1101       
1102        return answer;
1103}
1104
1105/*      FlexCyl_EllipX  :  calculates the form factor of a flexible cylinder with an elliptical cross section
1106-- incorporates Wei-Ren Chen's fixes - 2006
1107
1108        */
1109double
1110FlexCyl_Ellip(double dp[], double q)
1111{
1112        double scale,L,B,bkg,rad,qr,cont,ellRatio;
1113        double Pi,flex,crossSect,answer;
1114       
1115       
1116        Pi = 4.0*atan(1.0);
1117        scale = dp[0];          //make local copies in case memory moves
1118        L = dp[1];
1119        B = dp[2];
1120        rad = dp[3];
1121        ellRatio = dp[4];
1122        cont = dp[5];
1123        bkg = dp[6];
1124       
1125        qr = q*rad;
1126       
1127        flex = Sk_WR(q,L,B);
1128       
1129        crossSect = EllipticalCross_fn(q,rad,(rad*ellRatio));
1130        flex *= crossSect;
1131        flex *= Pi*rad*rad*ellRatio*L;
1132        flex *= cont*cont;
1133        flex *= 1.0e8;
1134        answer = scale*flex + bkg;
1135       
1136        return answer;
1137}
1138
1139double
1140EllipticalCross_fn(double qq, double a, double b)
1141{
1142    double uplim,lolim,Pi,summ,arg,zi,yyy,answer;
1143    int i,nord=76;
1144   
1145    Pi = 4.0*atan(1.0);
1146    lolim=0.0;
1147    uplim=Pi/2.0;
1148    summ=0.0;
1149   
1150    for(i=0;i<nord;i++) {
1151                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1152                arg = qq*sqrt(a*a*sin(zi)*sin(zi)+b*b*cos(zi)*cos(zi));
1153                yyy = pow((2.0 * NR_BessJ1(arg) / arg),2);
1154                yyy *= Gauss76Wt[i];
1155                summ += yyy;
1156    }
1157    answer = (uplim-lolim)/2.0*summ;
1158    answer *= 2.0/Pi;
1159    return(answer);
1160       
1161}
1162/*      FlexCyl_PolyLenX  :  calculates the form factor of a flecible cylinder at the given x-value p->x
1163the cylinder has a polydisperse Length
1164
1165*/
1166double
1167FlexCyl_PolyLen(double dp[], double q)
1168{
1169        int i;
1170        double scale,radius,length,pd,delrho,bkg,lb;            //local variables of coefficient wave
1171        int nord=20;                    //order of integration
1172        double uplim,lolim;             //upper and lower integration limits
1173        double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
1174        double range,zz,Pi;
1175       
1176        Pi = 4.0*atan(1.0);
1177        range = 3.4;
1178       
1179        summ = 0.0;                     //initialize intergral
1180        scale = dp[0];                  //make local copies in case memory moves
1181        length = dp[1];                 //radius
1182        pd = dp[2];                     // average length
1183        lb = dp[3];
1184        radius = dp[4];
1185        delrho = dp[5];
1186        bkg = dp[6];
1187       
1188        zz = (1.0/pd)*(1.0/pd) - 1.0;
1189       
1190        lolim = length*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
1191        if(lolim<0) {
1192                lolim = 0;
1193        }
1194        if(pd>0.3) {
1195                range = 3.4 + (pd-0.3)*18.0;
1196        }
1197        uplim = length*(1.0+range*pd);
1198       
1199        for(i=0;i<nord;i++) {
1200                zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1201                yyy = Gauss20Wt[i] * FlePolyLen_kernel(q,radius,length,lb,zz,delrho,zi);
1202                summ += yyy;
1203        }
1204       
1205        answer = (uplim-lolim)/2.0*summ;
1206        //normalize by average cylinder volume (first moment), using the average length
1207        Vpoly=Pi*radius*radius*length;
1208        answer /= Vpoly;
1209       
1210        answer *=delrho*delrho;
1211       
1212        //convert to [cm-1]
1213        answer *= 1.0e8;
1214        //Scale
1215        answer *= scale;
1216        // add in the background
1217        answer += bkg;
1218       
1219        return answer;
1220}
1221
1222/*      FlexCyl_PolyLenX  :  calculates the form factor of a flexible cylinder at the given x-value p->x
1223the cylinder has a polydisperse cross sectional radius
1224
1225*/
1226double
1227FlexCyl_PolyRad(double dp[], double q)
1228{
1229        int i;
1230        double scale,radius,length,pd,delrho,bkg,lb;            //local variables of coefficient wave
1231        int nord=76;                    //order of integration
1232        double uplim,lolim;             //upper and lower integration limits
1233        double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
1234        double range,zz,Pi;
1235       
1236       
1237        Pi = 4.0*atan(1.0);
1238        range = 3.4;
1239       
1240        summ = 0.0;                     //initialize intergral
1241       
1242        scale = dp[0];                  //make local copies in case memory moves
1243        length = dp[1];                 //radius
1244        lb = dp[2];                     // average length
1245        radius = dp[3];
1246        pd = dp[4];
1247        delrho = dp[5];
1248        bkg = dp[6];
1249       
1250        zz = (1.0/pd)*(1.0/pd) - 1.0;
1251       
1252        lolim = radius*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
1253        if(lolim<0) {
1254                lolim = 0;
1255        }
1256        if(pd>0.3) {
1257                range = 3.4 + (pd-0.3)*18.0;
1258        }
1259        uplim = radius*(1.0+range*pd);
1260       
1261        for(i=0;i<nord;i++) {
1262                //zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1263                //yyy = Gauss20Wt[i] * FlePolyRad_kernel(q,radius,length,lb,zz,delrho,zi);
1264                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1265                yyy = Gauss76Wt[i] * FlePolyRad_kernel(q,radius,length,lb,zz,delrho,zi);
1266                summ += yyy;
1267        }
1268       
1269        answer = (uplim-lolim)/2.0*summ;
1270        //normalize by average cylinder volume (second moment), using the average radius
1271        Vpoly = Pi*radius*radius*length*(zz+2.0)/(zz+1.0);
1272        answer /= Vpoly;
1273       
1274        answer *=delrho*delrho;
1275       
1276        //convert to [cm-1]
1277        answer *= 1.0e8;
1278        //Scale
1279        answer *= scale;
1280        // add in the background
1281        answer += bkg;
1282       
1283        return answer;
1284}
1285
1286/////////functions for WRC implementation of flexible cylinders
1287static double
1288Sk_WR(double q, double L, double b)
1289{
1290        //
1291        double p1,p2,p1short,p2short,q0,qconnect;
1292        double C,epsilon,ans,q0short,Sexvmodify,pi;
1293       
1294        pi = 4.0*atan(1.0);
1295       
1296        p1 = 4.12;
1297        p2 = 4.42;
1298        p1short = 5.36;
1299        p2short = 5.62;
1300        q0 = 3.1;
1301        qconnect = q0/b;
1302        //     
1303        q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0);
1304       
1305        //
1306        if(L/b > 10.0) {
1307                C = 3.06/pow((L/b),0.44);
1308                epsilon = 0.176;
1309        } else {
1310                C = 1.0;
1311                epsilon = 0.170;
1312        }
1313        //
1314       
1315        if( L > 4*b ) { // Longer Chains
1316                if (q*b <= 3.1) {               //Modified by Yun on Oct. 15,
1317                        Sexvmodify = Sexvnew(q, L, b);
1318                        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);
1319                } else { //q(i)*b > 3.1
1320                        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);
1321                } 
1322        } else { //L <= 4*b Shorter Chains
1323                if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) {
1324                        if (q*b<=0.01) {
1325                                ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0;
1326                        } else {
1327                                ans = Sdebye1(q,L,b);
1328                        }
1329                } else {        //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3)
1330                        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);
1331                }
1332        }
1333       
1334        return(ans);
1335        //return(a2long(q, L, b, p1, p2, q0));
1336}
1337
1338//WR named this w (too generic)
1339static double
1340w_WR(double x)
1341{
1342    double yy;
1343    yy = 0.5*(1 + tanh((x - 1.523)/0.1477));
1344       
1345    return (yy);
1346}
1347
1348//
1349static double
1350u1(double q, double L, double b)
1351{
1352    double yy;
1353       
1354    yy = Rgsquareshort(q,L,b)*q*q;
1355   
1356    return (yy);
1357}
1358
1359// was named u
1360static double
1361u_WR(double q, double L, double b)
1362{
1363    double yy;
1364    yy = Rgsquare(q,L,b)*q*q;
1365    return (yy);
1366}
1367
1368
1369
1370//
1371static double
1372Rgsquarezero(double q, double L, double b)
1373{   
1374    double yy;
1375    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))));
1376   
1377    return (yy);
1378}
1379
1380//
1381static double
1382Rgsquareshort(double q, double L, double b)
1383{   
1384    double yy;
1385    yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b);
1386   
1387    return (yy);
1388}
1389
1390//
1391static double
1392Rgsquare(double q, double L, double b)
1393{
1394    double yy;
1395    yy = AlphaSquare(L/b)*L*b/6.0;
1396   
1397    return (yy);
1398}
1399
1400//
1401static double
1402AlphaSquare(double x)
1403{   
1404    double yy;
1405    yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) );
1406       
1407    return (yy);
1408}
1409
1410// ?? funciton is not used - but should the log actually be log10???
1411static double
1412miu(double x)
1413{   
1414    double yy;
1415    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)));
1416   
1417    return (yy);
1418}
1419
1420//
1421static double
1422Sdebye(double q, double L, double b)
1423{   
1424    double yy;
1425    yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2));
1426       
1427    return (yy);
1428}
1429
1430//
1431static double
1432Sdebye1(double q, double L, double b)
1433{   
1434    double yy;
1435    yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) );
1436   
1437    return (yy);
1438}
1439
1440//
1441static double
1442Sexv(double q, double L, double b)
1443{   
1444    double yy,C1,C2,C3,miu,Rg2;
1445    C1=1.22;
1446    C2=0.4288;
1447    C3=-1.651;
1448    miu = 0.585;
1449       
1450    Rg2 = Rgsquare(q,L,b);
1451   
1452    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)));
1453   
1454    return (yy);
1455}
1456
1457// this must be WR modified version
1458static double
1459Sexvnew(double q, double L, double b)
1460{   
1461    double yy,C1,C2,C3,miu;
1462    double del=1.05,C_star2,Rg2;
1463   
1464    C1=1.22;
1465    C2=0.4288;
1466    C3=-1.651;
1467    miu = 0.585;
1468       
1469    //calculating the derivative to decide on the corection (cutoff) term?
1470    // I have modified this from WRs original code
1471   
1472    if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) {
1473        C_star2 = 0.0;
1474    } else {
1475        C_star2 = 1.0;
1476    }
1477       
1478    Rg2 = Rgsquare(q,L,b);
1479   
1480        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)));
1481       
1482    return (yy);
1483}
1484
1485// these are the messy ones
1486static double
1487a2short(double q, double L, double b, double p1short, double p2short, double q0)
1488{   
1489    double yy,Rg2_sh;
1490    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p;
1491    double pi;
1492       
1493    E = 2.718281828459045091;
1494        pi = 4.0*atan(1.0);
1495    Rg2_sh = Rgsquareshort(q,L,b);
1496    Rg2_sh2 = Rg2_sh*Rg2_sh;
1497    t1 = ((q0*q0*Rg2_sh)/(b*b));
1498    Et1 = pow(E,t1);
1499    Emt1 =pow(E,-t1); 
1500    q02 = q0*q0;
1501    q0p = pow(q0,(-4.0 + p2short) );
1502   
1503    //E is the number e
1504    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)))))));
1505       
1506    return (yy);
1507}
1508
1509//
1510static double
1511a1short(double q, double L, double b, double p1short, double p2short, double q0)
1512{   
1513    double yy,Rg2_sh;
1514    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3;
1515        double pi;
1516   
1517    E = 2.718281828459045091;
1518        pi = 4.0*atan(1.0);
1519    Rg2_sh = Rgsquareshort(q,L,b);
1520    Rg2_sh2 = Rg2_sh*Rg2_sh;
1521    b3 = b*b*b;
1522    t1 = ((q0*q0*Rg2_sh)/(b*b));
1523    Et1 = pow(E,t1);
1524    Emt1 =pow(E,-t1); 
1525    q02 = q0*q0;
1526    q0p = pow(q0,(-4.0 + p1short) );
1527   
1528    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)))))); 
1529       
1530    return(yy);
1531}
1532
1533// this one will be lots of trouble
1534static double
1535a2long(double q, double L, double b, double p1, double p2, double q0)
1536{
1537    double yy,C1,C2,C3,C4,C5,miu,C,Rg2;
1538    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi;
1539    double E,b2,b3,b4,q02,q03,q04,q05,Rg22;
1540   
1541    pi = 4.0*atan(1.0);
1542    E = 2.718281828459045091;
1543    if( L/b > 10.0) {
1544        C = 3.06/pow((L/b),0.44);
1545    } else {
1546        C = 1.0;
1547    }
1548       
1549    C1 = 1.22;
1550    C2 = 0.4288;
1551    C3 = -1.651;
1552    C4 = 1.523;
1553    C5 = 0.1477;
1554    miu = 0.585;
1555       
1556    Rg2 = Rgsquare(q,L,b);
1557    Rg22 = Rg2*Rg2;
1558    b2 = b*b;
1559    b3 = b*b*b;
1560    b4 = b3*b;
1561    q02 = q0*q0;
1562    q03 = q0*q0*q0;
1563    q04 = q03*q0;
1564    q05 = q04*q0;
1565   
1566    t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) ));
1567   
1568    t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + (14*b3*pow(E,(-((q02*Rg2)/b2))))/(15*q03*Rg2) + (2*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7*b2)/(15*q02*Rg2)))*Rg2)/b)))/L;
1569   
1570    t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(2*C5);
1571   
1572    t4 = (b4*sqrt(Rg2)*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22);
1573   
1574    t5 = (2*b4*(((2*q0*Rg2)/b - (2*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);
1575   
1576    t6 = (8*b4*b*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22);
1577   
1578    t7 = (((-((3*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 3/miu)))/miu)) - (2*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 2/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 1/miu)))/miu));
1579   
1580    t8 = ((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
1581   
1582    t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7*b2)/(15*q02*Rg2))) + (7*b2)/(15*q02*Rg2))))/L;
1583   
1584    t10 = (2*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);
1585       
1586   
1587    yy = ((-1*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1)/miu))))*((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))));
1588       
1589    return (yy);
1590}
1591
1592//need to define this on my own
1593static double
1594sech_WR(double x)
1595{       
1596        return(1/cosh(x));
1597}
1598
1599//
1600static double
1601a1long(double q, double L, double b, double p1, double p2, double q0)
1602{   
1603    double yy,C,C1,C2,C3,C4,C5,miu,Rg2;
1604    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15;
1605    double E,pi;
1606    double b2,b3,b4,q02,q03,q04,q05,Rg22;
1607       
1608    pi = 4.0*atan(1.0);
1609    E = 2.718281828459045091;
1610   
1611    if( L/b > 10.0) {
1612        C = 3.06/pow((L/b),0.44);
1613    } else {
1614        C = 1.0;
1615    }
1616       
1617    C1 = 1.22;
1618    C2 = 0.4288;
1619    C3 = -1.651;
1620    C4 = 1.523;
1621    C5 = 0.1477;
1622    miu = 0.585;
1623       
1624    Rg2 = Rgsquare(q,L,b);
1625    Rg22 = Rg2*Rg2;
1626    b2 = b*b;
1627    b3 = b*b*b;
1628    b4 = b3*b;
1629    q02 = q0*q0;
1630    q03 = q0*q0*q0;
1631    q04 = q03*q0;
1632    q05 = q04*q0;
1633   
1634    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))));
1635   
1636    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))))));
1637   
1638    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))));
1639   
1640    t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
1641       
1642    t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2))));
1643   
1644    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)));
1645   
1646    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));
1647   
1648    t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));
1649   
1650    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))))));
1651   
1652    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))))));
1653   
1654    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));
1655   
1656    t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
1657   
1658    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))));
1659   
1660    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))))));
1661   
1662    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))));
1663       
1664   
1665    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)))))))))));
1666       
1667    return (yy);
1668}
1669
1670
1671
1672///////////////
1673
1674//
1675//     FUNCTION gfn2:    CONTAINS F(Q,A,B,mu)**2  AS GIVEN
1676//                       BY (53) AND (56,57) IN CHEN AND
1677//                       KOTLARCHYK REFERENCE
1678//
1679//     <PROLATE ELLIPSOIDS>
1680//
1681double
1682gfn2(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq)
1683{
1684        // local variables
1685        double aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,gfn2,pi43,gfn,Pi;
1686       
1687        Pi = 4.0*atan(1.0);
1688       
1689        pi43=4.0/3.0*Pi;
1690        aa = crmaj;
1691        bb = crmin;
1692        u2 = (aa*aa*xx*xx + bb*bb*(1.0-xx*xx));
1693        ut2 = (trmaj*trmaj*xx*xx + trmin*trmin*(1.0-xx*xx));
1694        uq = sqrt(u2)*qq;
1695        ut= sqrt(ut2)*qq;
1696        vc = pi43*aa*bb*bb;
1697        vt = pi43*trmaj*trmin*trmin;
1698        gfnc = 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq*vc*delpc;
1699        gfnt = 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut*vt*delps;
1700        gfn = gfnc+gfnt;
1701        gfn2 = gfn*gfn;
1702       
1703        return (gfn2);
1704}
1705
1706//
1707//     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN
1708//                       BY (53) & (58-59) IN CHEN AND
1709//                       KOTLARCHYK REFERENCE
1710//
1711//       <OBLATE ELLIPSOID>
1712// function gfn4 for oblate ellipsoids
1713double
1714gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq)
1715{
1716        // local variables
1717        double aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn4,pi43,Pi;
1718       
1719        Pi = 4.0*atan(1.0);
1720        pi43=4.0/3.0*Pi;
1721        aa = crmaj;
1722        bb = crmin;
1723        u2 = (bb*bb*xx*xx + aa*aa*(1.0-xx*xx));
1724        ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx));
1725        uq = sqrt(u2)*qq;
1726        ut= sqrt(ut2)*qq;
1727        vc = pi43*aa*aa*bb;
1728        vt = pi43*trmaj*trmaj*trmin;
1729        gfnc = 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq*vc*delpc;
1730        gfnt = 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut*vt*delps;
1731        tgfn = gfnc+gfnt;
1732        gfn4 = tgfn*tgfn;
1733       
1734        return (gfn4);
1735}
1736
1737double
1738FlePolyLen_kernel(double q, double radius, double length, double lb, double zz, double delrho, double zi)
1739{
1740    double Pq,vcyl,dl;
1741    double Pi,qr;
1742   
1743    Pi = 4.0*atan(1.0);
1744    qr = q*radius;
1745   
1746    Pq = Sk_WR(q,zi,lb);                //does not have cross section term
1747    Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
1748   
1749    vcyl=Pi*radius*radius*zi;
1750    Pq *= vcyl*vcyl;
1751   
1752    dl = SchulzPoint_cpr(zi,length,zz);
1753    return (Pq*dl);     
1754       
1755}
1756
1757double
1758FlePolyRad_kernel(double q, double ravg, double Lc, double Lb, double zz, double delrho, double zi)
1759{
1760    double Pq,vcyl,dr;
1761    double Pi,qr;
1762   
1763    Pi = 4.0*atan(1.0);
1764    qr = q*zi;
1765   
1766    Pq = Sk_WR(q,Lc,Lb);                //does not have cross section term
1767    Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
1768   
1769    vcyl=Pi*zi*zi*Lc;
1770    Pq *= vcyl*vcyl;
1771   
1772    dr = SchulzPoint_cpr(zi,ravg,zz);
1773    return (Pq*dr);     
1774       
1775}
1776
1777double
1778CSCylIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length)
1779{
1780        double answer,halfheight,Pi;
1781        double lolim,uplim,summ,yyy,zi;
1782        int nord,i;
1783       
1784        // set up the integration end points
1785        Pi = 4.0*atan(1.0);
1786        nord = 76;
1787        lolim = 0;
1788        uplim = Pi/2;
1789        halfheight = length/2.0;
1790       
1791        summ = 0.0;                             // initialize integral
1792        i=0;
1793        for(i=0;i<nord;i++) {
1794                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1795                yyy = Gauss76Wt[i] * CScyl(qq, rad, radthick, facthick, rhoc,rhos,rhosolv, halfheight, zi);
1796                summ += yyy;
1797        }
1798       
1799        // calculate value of integral to return
1800        answer = (uplim-lolim)/2.0*summ;
1801        return (answer);
1802}
1803
1804double
1805CScyl(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length, double dum)
1806{       
1807        // qq is the q-value for the calculation (1/A)
1808        // radius is the core radius of the cylinder (A)
1809        //  radthick and facthick are the radial and face layer thicknesses
1810        // rho(n) are the respective SLD's
1811        // length is the *Half* CORE-LENGTH of the cylinder
1812        // dum is the dummy variable for the integration (theta)
1813       
1814        double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval;
1815        double Pi;
1816       
1817        Pi = 4.0*atan(1.0); 
1818       
1819        dr1 = rhoc-rhos;
1820        dr2 = rhos-rhosolv;
1821        vol1 = Pi*rad*rad*(2.0*length);
1822        vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick);
1823       
1824        besarg1 = qq*rad*sin(dum);
1825        besarg2 = qq*(rad+radthick)*sin(dum);
1826        sinarg1 = qq*length*cos(dum);
1827        sinarg2 = qq*(length+facthick)*cos(dum);
1828       
1829        t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1;
1830        t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2;
1831       
1832        retval = ((t1+t2)*(t1+t2))*sin(dum);
1833        return (retval);
1834   
1835}
1836
1837
1838double
1839CoreShellCylKernel(double qq, double rcore, double thick, double rhoc, double rhos, double rhosolv, double length, double dum)
1840{
1841       
1842    double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval;
1843    double Pi;
1844   
1845    Pi = 4.0*atan(1.0);
1846   
1847    dr1 = rhoc-rhos;
1848    dr2 = rhos-rhosolv;
1849    vol1 = Pi*rcore*rcore*(2.0*length);
1850    vol2 = Pi*(rcore+thick)*(rcore+thick)*(2.0*length+2.0*thick);
1851   
1852    besarg1 = qq*rcore*sin(dum);
1853    besarg2 = qq*(rcore+thick)*sin(dum);
1854    sinarg1 = qq*length*cos(dum);
1855    sinarg2 = qq*(length+thick)*cos(dum);
1856   
1857    t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1;
1858    t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2;
1859   
1860    retval = ((t1+t2)*(t1+t2))*sin(dum);
1861       
1862    return (retval);
1863}
1864
1865double
1866Cyl_PolyLenKernel(double q, double radius, double len_avg, double zz, double delrho, double dumLen)
1867{
1868       
1869    double halfheight,uplim,lolim,zi,summ,yyy,Pi;
1870    double answer,dr,Vcyl;
1871    int i,nord;
1872   
1873    Pi = 4.0*atan(1.0);
1874    lolim = 0;
1875    uplim = Pi/2.0;
1876    halfheight = dumLen/2.0;
1877    nord=20;
1878    summ = 0.0;
1879   
1880    //do the cylinder orientational average
1881    for(i=0;i<nord;i++) {
1882                zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1883                yyy = Gauss20Wt[i] * CylKernel(q, radius, halfheight, zi);
1884                summ += yyy;
1885    }
1886    answer = (uplim-lolim)/2.0*summ;
1887    // Multiply by contrast^2
1888    answer *= delrho*delrho;
1889    // don't do the normal scaling to volume here
1890    // instead, multiply by VCyl^2 to get rid of the normalization for this radius of cylinder
1891    Vcyl = Pi*radius*radius*dumLen;
1892    answer *= Vcyl*Vcyl;
1893   
1894    dr = SchulzPoint_cpr(dumLen,len_avg,zz);
1895    return(dr*answer);
1896}
1897
1898
1899double
1900Stackdisc_kern(double qq, double rcore, double rhoc, double rhol, double rhosolv, double length, double thick, double dum, double gsd, double d, double N)
1901{               
1902        // qq is the q-value for the calculation (1/A)
1903        // rcore is the core radius of the cylinder (A)
1904        // rho(n) are the respective SLD's
1905        // length is the *Half* CORE-LENGTH of the cylinder = L (A)
1906        // dum is the dummy variable for the integration (x in Feigin's notation)
1907       
1908        //Local variables
1909        double totald,dr1,dr2,besarg1,besarg2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt;
1910        double Pi;
1911        int kk;
1912       
1913        Pi = 4.0*atan(1.0);
1914       
1915        dr1 = rhoc-rhosolv;
1916        dr2 = rhol-rhosolv;
1917        area = Pi*rcore*rcore;
1918        totald=2.0*(thick+length);
1919       
1920        besarg1 = qq*rcore*sin(dum);
1921        besarg2 = qq*rcore*sin(dum);
1922       
1923        sinarg1 = qq*length*cos(dum);
1924        sinarg2 = qq*(length+thick)*cos(dum);
1925       
1926        t1 = 2*area*(2*length)*dr1*(sin(sinarg1)/sinarg1)*(NR_BessJ1(besarg1)/besarg1);
1927        t2 = 2*area*dr2*(totald*sin(sinarg2)/sinarg2-2*length*sin(sinarg1)/sinarg1)*(NR_BessJ1(besarg2)/besarg2);
1928       
1929        retval =((t1+t2)*(t1+t2))*sin(dum);
1930       
1931        // loop for the structure facture S(q)
1932        sqq=0.0;
1933        for(kk=1;kk<N;kk+=1) {
1934                dexpt=qq*cos(dum)*qq*cos(dum)*d*d*gsd*gsd*kk/2.0;
1935                sqq=sqq+(N-kk)*cos(qq*cos(dum)*d*kk)*exp(-1.*dexpt);
1936        }                       
1937       
1938        // end of loop for S(q)
1939        sqq=1.0+2.0*sqq/N;
1940        retval *= sqq;
1941   
1942        return(retval);
1943}
1944
1945
1946double
1947Cyl_PolyRadKernel(double q, double radius, double length, double zz, double delrho, double dumRad)
1948{
1949       
1950    double halfheight,uplim,lolim,zi,summ,yyy,Pi;
1951    double answer,dr,Vcyl;
1952    int i,nord;
1953   
1954    Pi = 4.0*atan(1.0);
1955    lolim = 0;
1956    uplim = Pi/2.0;
1957    halfheight = length/2.0;
1958        //    nord=20;
1959    nord=76;
1960    summ = 0.0;
1961   
1962    //do the cylinder orientational average
1963        //    for(i=0;i<nord;i++) {
1964        //            zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1965        //            yyy = Gauss20Wt[i] * CylKernel(q, dumRad, halfheight, zi);
1966        //            summ += yyy;
1967        //    }
1968    for(i=0;i<nord;i++) {
1969                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
1970                yyy = Gauss76Wt[i] * CylKernel(q, dumRad, halfheight, zi);
1971                summ += yyy;
1972    }
1973    answer = (uplim-lolim)/2.0*summ;
1974    // Multiply by contrast^2
1975    answer *= delrho*delrho;
1976    // don't do the normal scaling to volume here
1977    // instead, multiply by VCyl^2 to get rid of the normalization for this radius of cylinder
1978    Vcyl = Pi*dumRad*dumRad*length;
1979    answer *= Vcyl*Vcyl;
1980   
1981    dr = SchulzPoint_cpr(dumRad,radius,zz);
1982    return(dr*answer);
1983}
1984
1985double
1986SchulzPoint_cpr(double dumRad, double radius, double zz)
1987{
1988    double dr;
1989   
1990    dr = zz*log(dumRad) - gammaln(zz+1.0) + (zz+1.0)*log((zz+1.0)/radius)-(dumRad/radius*(zz+1.0));
1991    return(exp(dr));
1992}
1993
1994static double
1995gammaln(double xx)
1996{
1997        double x,y,tmp,ser;
1998        static double cof[6]={76.18009172947146,-86.50532032941677,
1999                24.01409824083091,-1.231739572450155,
2000                0.1208650973866179e-2,-0.5395239384953e-5};
2001        int j;
2002       
2003        y=x=xx;
2004        tmp=x+5.5;
2005        tmp -= (x+0.5)*log(tmp);
2006        ser=1.000000000190015;
2007        for (j=0;j<=5;j++) ser += cof[j]/++y;
2008        return -tmp+log(2.5066282746310005*ser/x);
2009}
2010
2011
2012double
2013EllipsoidKernel(double qq, double a, double nua, double dum)
2014{
2015    double arg,nu,retval;               //local variables
2016       
2017    nu = nua/a;
2018    arg = qq*a*sqrt(1+dum*dum*(nu*nu-1));
2019   
2020    retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg);
2021    retval *= retval;
2022    retval *= 9;
2023   
2024    return(retval);
2025}//Function EllipsoidKernel()
2026
2027double
2028HolCylKernel(double qq, double rcore, double rshell, double length, double dum)
2029{
2030    double gamma,arg1,arg2,lam1,lam2,psi,sinarg,t2,retval;              //local variables
2031   
2032    gamma = rcore/rshell;
2033    arg1 = qq*rshell*sqrt(1-dum*dum);           //1=shell (outer radius)
2034    arg2 = qq*rcore*sqrt(1-dum*dum);                    //2=core (inner radius)
2035    lam1 = 2*NR_BessJ1(arg1)/arg1;
2036    lam2 = 2*NR_BessJ1(arg2)/arg2;
2037    psi = 1/(1-gamma*gamma)*(lam1 -  gamma*gamma*lam2);         //SRK 10/19/00
2038   
2039    sinarg = qq*length*dum/2;
2040    t2 = sin(sinarg)/sinarg;
2041   
2042    retval = psi*psi*t2*t2;
2043   
2044    return(retval);
2045}//Function HolCylKernel()
2046
2047double
2048PPKernel(double aa, double mu, double uu)
2049{
2050    // mu passed in is really mu*sqrt(1-sig^2)
2051    double arg1,arg2,Pi,tmp1,tmp2;                      //local variables
2052       
2053    Pi = 4.0*atan(1.0);
2054   
2055    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
2056    arg1 = (mu/2)*cos(Pi*uu/2);
2057    arg2 = (mu*aa/2)*sin(Pi*uu/2);
2058    if(arg1==0) {
2059                tmp1 = 1;
2060        } else {
2061                tmp1 = sin(arg1)*sin(arg1)/arg1/arg1;
2062    }
2063   
2064    if (arg2==0) {
2065                tmp2 = 1;
2066        } else {
2067                tmp2 = sin(arg2)*sin(arg2)/arg2/arg2;
2068    }
2069       
2070    return (tmp1*tmp2);
2071       
2072}//Function PPKernel()
2073
2074
2075double
2076TriaxialKernel(double q, double aa, double bb, double cc, double dx, double dy)
2077{
2078       
2079    double arg,val,pi;                  //local variables
2080       
2081    pi = 4.0*atan(1.0);
2082   
2083    arg = aa*aa*cos(pi*dx/2)*cos(pi*dx/2);
2084    arg += bb*bb*sin(pi*dx/2)*sin(pi*dx/2)*(1-dy*dy);
2085    arg += cc*cc*dy*dy;
2086    arg = q*sqrt(arg);
2087    val = 9 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) );
2088   
2089    return (val);
2090       
2091}//Function TriaxialKernel()
2092
2093
2094double
2095CylKernel(double qq, double rr,double h, double theta)
2096{
2097       
2098        // qq is the q-value for the calculation (1/A)
2099        // rr is the radius of the cylinder (A)
2100        // h is the HALF-LENGTH of the cylinder = L/2 (A)
2101       
2102    double besarg,bj,retval,d1,t1,b1,t2,b2;              //Local variables
2103       
2104       
2105    besarg = qq*rr*sin(theta);
2106       
2107    bj =NR_BessJ1(besarg);
2108       
2109        //* Computing 2nd power */
2110    d1 = sin(qq * h * cos(theta));
2111    t1 = d1 * d1;
2112        //* Computing 2nd power */
2113    d1 = bj;
2114    t2 = d1 * d1 * 4.0 * sin(theta);
2115        //* Computing 2nd power */
2116    d1 = qq * h * cos(theta);
2117    b1 = d1 * d1;
2118        //* Computing 2nd power */
2119    d1 = qq * rr * sin(theta);
2120    b2 = d1 * d1;
2121    retval = t1 * t2 / b1 / b2;
2122       
2123    return (retval);
2124       
2125}//Function CylKernel()
2126
2127double
2128EllipCylKernel(double qq, double ra,double nu, double theta)
2129{
2130        //this is the function LAMBDA1^2 in Feigin's notation   
2131        // qq is the q-value for the calculation (1/A)
2132        // ra is the transformed radius"a" in Feigin's notation
2133        // nu is the ratio (major radius)/(minor radius) of the Ellipsoid [=] ---
2134        // theta is the dummy variable of the integration
2135       
2136        double retval,arg;               //Local variables
2137       
2138        arg = qq*ra*sqrt((1+nu*nu)/2+(1-nu*nu)*cos(theta)/2);
2139       
2140        retval = 2*NR_BessJ1(arg)/arg;
2141       
2142        //square it
2143        retval *= retval;
2144       
2145    return(retval);
2146       
2147}//Function EllipCylKernel()
2148
2149double NR_BessJ1(double x)
2150{
2151        double ax,z;
2152        double xx,y,ans,ans1,ans2;
2153       
2154        if ((ax=fabs(x)) < 8.0) {
2155                y=x*x;
2156                ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
2157                                                                                                  +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
2158                ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
2159                                                                                           +y*(99447.43394+y*(376.9991397+y*1.0))));
2160                ans=ans1/ans2;
2161        } else {
2162                z=8.0/ax;
2163                y=z*z;
2164                xx=ax-2.356194491;
2165                ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
2166                                                                   +y*(0.2457520174e-5+y*(-0.240337019e-6))));
2167                ans2=0.04687499995+y*(-0.2002690873e-3
2168                                                          +y*(0.8449199096e-5+y*(-0.88228987e-6
2169                                                                                                         +y*0.105787412e-6)));
2170                ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
2171                if (x < 0.0) ans = -ans;
2172        }
2173       
2174        return(ans);
2175}
Note: See TracBrowser for help on using the repository browser.