source: sasview/sansmodels/src/libigor/libCylinder.c @ 8e91f01

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 8e91f01 was 8e91f01, checked in by Jae Cho <jhjcho@…>, 15 years ago

set smear = 0 in LamellaPS instead of using arbitrary constant 0.025

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