source: sans/Analysis/branches/ajj_23APR07/XOPs/SANSAnalysis/Cylinder.c @ 93

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

Add combined XOP code. Currently only contains XCode project file to build Universal binary suitable for Igor 6.

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