source: sasview/sansmodels/src/libigor/libCylinder.c @ f9a1279

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since f9a1279 was 222d75c7, checked in by Jae Cho <jhjcho@…>, 15 years ago

reverted added lamellar_kernel

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