source: sasview/sansmodels/src/sans/models/libigor/libCylinder.h @ e161fac

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since e161fac was 6a2c931, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Re #5 fixing samsmodels compilation on MSVC
  • Property mode set to 100644
File size: 15.3 KB
Line 
1/*
2        libCylinderFit.h -- equates for CylinderFit XOP
3*/
4#if defined(_MSC_VER)
5#include "winFuncs.h"
6#endif
7
8/* Prototypes */
9/* IGOR Fit Functions */
10double CylinderForm(double dp[], double q);
11double EllipCyl76(double dp[], double q);
12double EllipCyl20(double dp[], double q);
13double TriaxialEllipsoid(double dp[], double q);
14double Parallelepiped(double dp[], double q);
15double HollowCylinder(double dp[], double q);
16double EllipsoidForm(double dp[], double q);
17double Cyl_PolyRadius(double dp[], double q);
18double Cyl_PolyLength(double dp[], double q);
19double CoreShellCylinder(double dp[], double q);
20double OblateForm(double dp[], double q);
21double ProlateForm(double dp[], double q);
22double FlexExclVolCyl(double dp[], double q);
23double FlexCyl_PolyLen(double dp[], double q);
24double FlexCyl_PolyRad(double dp[], double q);
25double FlexCyl_Ellip(double dp[], double q);
26double PolyCoShCylinder(double dp[], double q);
27double StackedDiscs(double dp[], double q);
28double LamellarFF(double dp[], double q);
29double LamellarFF_HG(double dp[], double q);
30double LamellarPS(double dp[], double q);
31double LamellarPS_HG(double dp[], double q);
32double Lamellar_ParaCrystal(double dp[], double q);
33double Spherocylinder(double dp[], double q);
34double ConvexLens(double dp[], double q);
35double Dumbbell(double dp[], double q);
36double CappedCylinder(double dp[], double q);
37double Barbell(double dp[], double q);
38double PolyCoreBicelle(double dp[], double q);
39double CSParallelepiped(double dp[], double q);
40
41/* internal functions */
42double CylKernel(double qq, double rr,double h, double theta);
43double NR_BessJ1(double x);
44double EllipCylKernel(double qq, double ra,double nu, double theta);
45double TriaxialKernel(double q, double aa, double bb, double cc, double dx, double dy);
46double PPKernel(double aa, double mu, double uu);
47double HolCylKernel(double qq, double rcore, double rshell, double length, double dum);
48double EllipsoidKernel(double qq, double a, double va, double dum);
49double Cyl_PolyRadKernel(double q, double radius, double length, double zz, double delrho, double dumRad);
50double SchulzPoint_cpr(double dumRad, double radius, double zz);
51double Cyl_PolyLenKernel(double q, double radius, double len_avg, double zz, double delrho, double dumLen);
52double CoreShellCylKernel(double qq, double rcore, double thick, double rhoc, double rhos, double rhosolv, double length, double dum);
53double gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq);
54double gfn2(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq);
55double FlePolyLen_kernel(double q, double radius, double length, double lb, double zz, double delrho, double zi);
56double FlePolyRad_kernel(double q, double ravg, double Lc, double Lb, double zz, double delrho, double zi);
57double EllipticalCross_fn(double qq, double a, double b);
58double CScyl(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length, double dum);
59double CSCylIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length);
60double Stackdisc_kern(double qq, double rcore, double rhoc, double rhol, double rhosolv, double length, double thick, double dum, double gsd, double d, double N);
61double paraCryst_sn(double ww, double qval, double davg, long nl, double an);
62double paraCryst_an(double ww, double qval, double davg, long nl);
63double SphCyl_kernel(double w[], double x, double tt, double Theta);
64double ConvLens_kernel(double w[], double x, double tt, double theta);
65double Dumb_kernel(double w[], double x, double tt, double theta);
66double BicelleKernel(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length, double dum);
67double BicelleIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length);
68double CSPPKernel(double dp[], double mu, double uu);
69
70/////////functions for WRC implementation of flexible cylinders
71static double
72gammaln(double xx)
73{
74        double x,y,tmp,ser;
75        static double cof[6]={76.18009172947146,-86.50532032941677,
76                24.01409824083091,-1.231739572450155,
77                0.1208650973866179e-2,-0.5395239384953e-5};
78        int j;
79
80        y=x=xx;
81        tmp=x+5.5;
82        tmp -= (x+0.5)*log(tmp);
83        ser=1.000000000190015;
84        for (j=0;j<=5;j++) ser += cof[j]/++y;
85        return -tmp+log(2.5066282746310005*ser/x);
86}
87
88//
89static double
90AlphaSquare(double x)
91{
92    double yy;
93    yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) );
94
95    return (yy);
96}
97
98//
99static double
100Rgsquarezero(double q, double L, double b)
101{
102    double yy;
103    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))));
104
105    return (yy);
106}
107
108//
109static double
110Rgsquareshort(double q, double L, double b)
111{
112    double yy;
113    yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b);
114
115    return (yy);
116}
117
118//
119static double
120Rgsquare(double q, double L, double b)
121{
122    double yy;
123    yy = AlphaSquare(L/b)*L*b/6.0;
124
125    return (yy);
126}
127
128// ?? funciton is not used - but should the log actually be log10???
129static double
130miu(double x)
131{
132    double yy;
133    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)));
134
135    return (yy);
136}
137
138//WR named this w (too generic)
139static double
140w_WR(double x)
141{
142    double yy;
143    yy = 0.5*(1 + tanh((x - 1.523)/0.1477));
144
145    return (yy);
146}
147
148//
149static double
150u1(double q, double L, double b)
151{
152    double yy;
153
154    yy = Rgsquareshort(q,L,b)*q*q;
155
156    return (yy);
157}
158
159// was named u
160static double
161u_WR(double q, double L, double b)
162{
163    double yy;
164    yy = Rgsquare(q,L,b)*q*q;
165    return (yy);
166}
167
168
169//
170static double
171Sdebye1(double q, double L, double b)
172{
173    double yy;
174    yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) );
175
176    return (yy);
177}
178
179
180//
181static double
182Sdebye(double q, double L, double b)
183{
184    double yy;
185    yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2));
186
187    return (yy);
188}
189
190//
191static double
192Sexv(double q, double L, double b)
193{
194    double yy,C1,C2,C3,miu,Rg2;
195    C1=1.22;
196    C2=0.4288;
197    C3=-1.651;
198    miu = 0.585;
199
200    Rg2 = Rgsquare(q,L,b);
201
202    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)));
203
204    return (yy);
205}
206
207// this must be WR modified version
208static double
209Sexvnew(double q, double L, double b)
210{
211    double yy,C1,C2,C3,miu;
212    double del=1.05,C_star2,Rg2;
213
214    C1=1.22;
215    C2=0.4288;
216    C3=-1.651;
217    miu = 0.585;
218
219    //calculating the derivative to decide on the corection (cutoff) term?
220    // I have modified this from WRs original code
221
222    if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) {
223        C_star2 = 0.0;
224    } else {
225        C_star2 = 1.0;
226    }
227
228    Rg2 = Rgsquare(q,L,b);
229
230        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)));
231
232    return (yy);
233}
234
235// these are the messy ones
236static double
237a2short(double q, double L, double b, double p1short, double p2short, double q0)
238{
239    double yy,Rg2_sh;
240    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p;
241    double pi;
242
243    E = 2.718281828459045091;
244        pi = 4.0*atan(1.0);
245    Rg2_sh = Rgsquareshort(q,L,b);
246    Rg2_sh2 = Rg2_sh*Rg2_sh;
247    t1 = ((q0*q0*Rg2_sh)/(b*b));
248    Et1 = pow(E,t1);
249    Emt1 =pow(E,-t1);
250    q02 = q0*q0;
251    q0p = pow(q0,(-4.0 + p2short) );
252
253    //E is the number e
254    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)))))));
255
256    return (yy);
257}
258
259//
260static double
261a1short(double q, double L, double b, double p1short, double p2short, double q0)
262{
263    double yy,Rg2_sh;
264    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3;
265        double pi;
266
267    E = 2.718281828459045091;
268        pi = 4.0*atan(1.0);
269    Rg2_sh = Rgsquareshort(q,L,b);
270    Rg2_sh2 = Rg2_sh*Rg2_sh;
271    b3 = b*b*b;
272    t1 = ((q0*q0*Rg2_sh)/(b*b));
273    Et1 = pow(E,t1);
274    Emt1 =pow(E,-t1);
275    q02 = q0*q0;
276    q0p = pow(q0,(-4.0 + p1short) );
277
278    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))))));
279
280    return(yy);
281}
282
283
284//need to define this on my own
285static double
286sech_WR(double x)
287{
288        return(1/cosh(x));
289}
290
291// this one will be lots of trouble
292static double
293a2long(double q, double L, double b, double p1, double p2, double q0)
294{
295    double yy,C1,C2,C3,C4,C5,miu,C,Rg2;
296    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi;
297    double E,b2,b3,b4,q02,q03,q04,q05,Rg22;
298
299    pi = 4.0*atan(1.0);
300    E = 2.718281828459045091;
301    if( L/b > 10.0) {
302        C = 3.06/pow((L/b),0.44);
303    } else {
304        C = 1.0;
305    }
306
307    C1 = 1.22;
308    C2 = 0.4288;
309    C3 = -1.651;
310    C4 = 1.523;
311    C5 = 0.1477;
312    miu = 0.585;
313
314    Rg2 = Rgsquare(q,L,b);
315    Rg22 = Rg2*Rg2;
316    b2 = b*b;
317    b3 = b*b*b;
318    b4 = b3*b;
319    q02 = q0*q0;
320    q03 = q0*q0*q0;
321    q04 = q03*q0;
322    q05 = q04*q0;
323
324    t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) ));
325
326    t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7*b2)/(15.0*q02*Rg2)))*Rg2)/b)))/L;
327
328    t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2.0))/(2.0*C5);
329
330    t4 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22);
331
332    t5 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);
333
334    t6 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22);
335
336    t7 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu));
337
338    t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
339
340    t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))))/L;
341
342    t10 = (2.0*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);
343
344
345    yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))));
346
347    return (yy);
348}
349//
350static double
351a1long(double q, double L, double b, double p1, double p2, double q0)
352{
353    double yy,C,C1,C2,C3,C4,C5,miu,Rg2;
354    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15;
355    double E,pi;
356    double b2,b3,b4,q02,q03,q04,q05,Rg22;
357
358    pi = 4.0*atan(1.0);
359    E = 2.718281828459045091;
360
361    if( L/b > 10.0) {
362        C = 3.06/pow((L/b),0.44);
363    } else {
364        C = 1.0;
365    }
366
367    C1 = 1.22;
368    C2 = 0.4288;
369    C3 = -1.651;
370    C4 = 1.523;
371    C5 = 0.1477;
372    miu = 0.585;
373
374    Rg2 = Rgsquare(q,L,b);
375    Rg22 = Rg2*Rg2;
376    b2 = b*b;
377    b3 = b*b*b;
378    b4 = b3*b;
379    q02 = q0*q0;
380    q03 = q0*q0*q0;
381    q04 = q03*q0;
382    q05 = q04*q0;
383
384    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))));
385
386    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))))));
387
388    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))));
389
390    t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
391
392    t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2))));
393
394    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)));
395
396    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));
397
398    t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));
399
400    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))))));
401
402    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))))));
403
404    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));
405
406    t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
407
408    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))));
409
410    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))))));
411
412    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))));
413
414
415    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)))))))))));
416
417    return (yy);
418}
419
420
421static double
422Sk_WR(double q, double L, double b)
423{
424        //
425        double p1,p2,p1short,p2short,q0;
426        double C,ans,q0short,Sexvmodify,pi;
427
428        pi = 4.0*atan(1.0);
429
430        p1 = 4.12;
431        p2 = 4.42;
432        p1short = 5.36;
433        p2short = 5.62;
434        q0 = 3.1;
435        //
436        q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0);
437
438        //
439        if(L/b > 10.0) {
440                C = 3.06/pow((L/b),0.44);
441        } else {
442                C = 1.0;
443        }
444        //
445
446        if( L > 4*b ) { // Longer Chains
447                if (q*b <= 3.1) {               //Modified by Yun on Oct. 15,
448                        Sexvmodify = Sexvnew(q, L, b);
449                        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);
450                } else { //q(i)*b > 3.1
451                        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);
452                }
453        } else { //L <= 4*b Shorter Chains
454                if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) {
455                        if (q*b<=0.01) {
456                                ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0;
457                        } else {
458                                ans = Sdebye1(q,L,b);
459                        }
460                } else {        //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3)
461                        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);
462                }
463        }
464
465        return(ans);
466        //return(a2long(q, L, b, p1, p2, q0));
467}
468
Note: See TracBrowser for help on using the repository browser.