source: sasmodels/sasmodels/models/lib/wrc_cyl.c @ db1c84b

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since db1c84b was db1c84b, checked in by smk78, 16 months ago

Added WRC reference to wrc_cyl.c

  • Property mode set to 100644
File size: 7.2 KB
Line 
1/*
2    Functions for WRC implementation of flexible cylinders. See
3    W R Chen, P D Butler and L J Magid,
4    Incorporating Intermicellar Interactions in the Fitting of
5    SANS Data from Cationic Wormlike Micelles.
6    Langmuir, 22(15) 2006 6539-6548
7*/
8
9static double
10Rgsquare(double L, double b)
11{
12    const double x = L/b;
13    // Use Horner's method to evaluate Pedersen eq 15:
14    //     alpha^2 = [1.0 + (x/3.12)^2 + (x/8.67)^3] ^ (0.176/3)
15    const double alphasq =
16        pow(1.0 + x*x*(1.027284681130835e-01 + 1.534414548417740e-03*x),
17            5.866666666666667e-02);
18    return alphasq*L*b/6.0;
19}
20
21static double
22Rgsquareshort(double L, double b)
23{
24    const double r = b/L;  // = 1/n_b in Pedersen ref.
25    return Rgsquare(L, b) * (1.0 + r*(-1.5 + r*(1.5 + r*0.75*expm1(-2.0/r))));
26}
27
28static double
29w_WR(double x)
30{
31    // Pedersen eq. 16:
32    //    w = [1 + tanh((x-C4)/C5)]/2
33    const double C4 = 1.523;
34    const double C5 = 0.1477;
35    return 0.5 + 0.5*tanh((x - C4)/C5);
36}
37
38static double
39Sdebye(double qsq)
40{
41#if FLOAT_SIZE>4
42#  define DEBYE_CUTOFF 0.25  // 1e-15 error
43#else
44#  define DEBYE_CUTOFF 1.0  // 4e-7 error
45#endif
46
47    if (qsq < DEBYE_CUTOFF) {
48        const double x = qsq;
49        // mathematica: PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}]
50        const double A1=1./15., A2=1./60, A3=0., A4=1./75600.;
51        const double B1=2./5., B2=1./15., B3=1./180., B4=1./5040.;
52        return ((((A4*x + A3)*x + A2)*x + A1)*x + 1.)
53             / ((((B4*x + B3)*x + B2)*x + B1)*x + 1.);
54    } else {
55        return 2.*(expm1(-qsq) + qsq)/(qsq*qsq);
56    }
57}
58
59static double
60a_long(double q, double L, double b/*, double p1, double p2, double q0*/)
61{
62    const double p1 = 4.12;
63    const double p2 = 4.42;
64    const double q0 = 3.1;
65
66    // Constants C1, ..., C5 derived from least squares fit (Pedersen, eq 13,16)
67    const double C1 = 1.22;
68    const double C2 = 0.4288;
69    const double C3 = -1.651;
70    const double C4 = 1.523;
71    const double C5 = 0.1477;
72    const double miu = 0.585;
73
74    const double C = (L/b>10.0 ? 3.06*pow(L/b, -0.44) : 1.0);
75    //printf("branch B-%d q=%g L=%g b=%g\n", C==1.0, q, L, b);
76    const double r2 = Rgsquare(L,b);
77    const double r = sqrt(r2);
78    const double qr_b = q0*r/b;
79    const double qr_b_sq = qr_b*qr_b;
80    const double qr_b_4 = qr_b_sq*qr_b_sq;
81    const double qr_b_miu = pow(qr_b, -1.0/miu);
82    const double em1_qr_b_sq = expm1(-qr_b_sq);
83    const double sech2 = 1.0/square(cosh((qr_b-C4)/C5));
84    const double w = w_WR(qr_b);
85
86    const double t1 = pow(q0, 1.0 + p1 + p2)/(b*(p1-p2));
87    const double t2 = C/(15.0*L) * (
88        + 14.0*b*b*em1_qr_b_sq/(q0*qr_b_sq)
89        + 2.0*q0*r2*exp(-qr_b_sq)*(11.0 + 7.0/qr_b_sq));
90    const double t11 = ((C3*qr_b_miu + C2)*qr_b_miu + C1)*qr_b_miu;
91    const double t3 = r*sech2/(2.*C5)*t11;
92    const double t4 = r*(em1_qr_b_sq + qr_b_sq)*sech2 / (C5*qr_b_4);
93    const double t5 = -4.0*r*qr_b*em1_qr_b_sq/qr_b_4 * (1.0 - w);
94    const double t10 = 2.0*(em1_qr_b_sq + qr_b_sq)/qr_b_4 * (1.0 - w); //=Sdebye*(1-w)
95    const double t6 = 4.0*b/q0 * t10;
96    const double t7 = r*((-3.0*C3*qr_b_miu -2.0*C2)*qr_b_miu -1.0*C1)*qr_b_miu/(miu*qr_b);
97    const double t9 = C*b/L * (4.0 - exp(-qr_b_sq) * (11.0 + 7.0/qr_b_sq) + 7.0/qr_b_sq)/15.0;
98    const double t12 = b*b*M_PI/(L*q0*q0) + t2 + t3 - t4 + t5 - t6 + t7*w;
99    const double t13 = -b*M_PI/(L*q0) + t9 + t10 + t11*w;
100
101    const double a1 = pow(q0,p1)*t13 - t1*pow(q0,-p2)*(t12 + b*p1/q0*t13);
102    const double a2 = t1*pow(q0,-p1)*(t12 + b*p1/q0*t13);
103
104    const double ans = a1*pow(q*b, -p1) + a2*pow(q*b, -p2) + M_PI/(q*L);
105    return ans;
106}
107
108static double
109_short(double r2, double exp_qr_b, double L, double b, double p1short,
110        double p2short, double q0)
111{
112    const double qr2 = q0*q0 * r2;
113    const double b3 = b*b*b;
114    const double q0p = pow(q0, -4.0 + p1short);
115
116    double yy = 1.0/(L*r2*r2) * b/exp_qr_b*q0p
117        * (8.0*b3*L
118           - 8.0*b3*exp_qr_b*L
119           + 2.0*b3*exp_qr_b*L*p2short
120           - 2.0*b*exp_qr_b*L*p2short*qr2
121           + 4.0*b*exp_qr_b*L*qr2
122           - 2.0*b3*L*p2short
123           + 4.0*b*L*qr2
124           - M_PI*exp_qr_b*qr2*q0*r2
125           + M_PI*exp_qr_b*p2short*qr2*q0*r2);
126
127    return yy;
128}
129static double
130a_short(double qp, double L, double b
131        /*double p1short, double p2short*/, double q0)
132{
133    const double p1short = 5.36;
134    const double p2short = 5.62;
135
136    const double r2 = Rgsquareshort(L,b);
137    const double exp_qr_b = exp(r2*square(q0/b));
138    const double pdiff = p1short - p2short;
139    const double a1 = _short(r2,exp_qr_b,L,b,p1short,p2short,q0)/pdiff;
140    const double a2= -_short(r2,exp_qr_b,L,b,p2short,p1short,q0)/pdiff;
141    const double ans = a1*pow(qp*b, -p1short) + a2*pow(qp*b, -p2short) + M_PI/(qp*L);
142    return ans;
143}
144
145
146static double
147Sexv(double q, double L, double b)
148{
149    // Pedersen eq 13, corrected by Chen eq A.5, swapping w and 1-w
150    const double C1=1.22;
151    const double C2=0.4288;
152    const double C3=-1.651;
153    const double miu = 0.585;
154    const double qr = q*sqrt(Rgsquare(L,b));
155    const double qr_miu = pow(qr, -1.0/miu);
156    const double w = w_WR(qr);
157    const double t10 = Sdebye(qr*qr)*(1.0 - w);
158    const double t11 = ((C3*qr_miu + C2)*qr_miu + C1)*qr_miu;
159
160    return t10 + w*t11;
161}
162
163// Modified by Yun on Oct. 15,
164static double
165Sexv_new(double q, double L, double b)
166{
167    const double qr = q*sqrt(Rgsquare(L,b));
168    const double qr2 = qr*qr;
169    const double C = (L/b > 10.0) ? 3.06*pow(L/b, -0.44) : 1.0;
170    const double t9 = C*b/L * (4.0 - exp(-qr2) * (11.0 + 7.0/qr2) + 7.0/qr2)/15.0;
171
172    const double Sexv_orig = Sexv(q, L, b);
173
174    // calculating the derivative to decide on the correction (cutoff) term?
175    // Note: this is modified from WRs original code
176    const double del=1.05;
177    const double qdel = (Sexv(q*del,L,b) - Sexv_orig)/(q*(del - 1.0));
178
179    if (qdel < 0) {
180        //printf("branch A1-%d q=%g L=%g b=%g\n", C==1.0, q, L, b);
181        return t9 + Sexv_orig;
182    } else {
183        //printf("branch A2-%d q=%g L=%g b=%g\n", C==1.0, q, L, b);
184        const double w = w_WR(qr);
185        const double t10 = Sdebye(qr*qr)*(1.0 - w);
186        return t9 + t10;
187    }
188}
189
190
191static double
192Sk_WR(double q, double L, double b)
193{
194    const double Rg_short = sqrt(Rgsquareshort(L, b));
195    double q0short = fmax(1.9/Rg_short, 3.0);
196    double ans;
197
198
199    if( L > 4*b ) { // L > 4*b : Longer Chains
200        if (q*b <= 3.1) {
201            ans = Sexv_new(q, L, b);
202        } else { //q(i)*b > 3.1
203            ans = a_long(q, L, b /*, p1, p2, q0*/);
204        }
205    } else { // L <= 4*b : Shorter Chains
206        if (q*b <= q0short) { // q*b <= fmax(1.9/Rg_short, 3)
207            //printf("branch C-%d q=%g L=%g b=%g\n", square(q*Rg_short)<DEBYE_CUTOFF, q, L, b);
208            // Note that q0short is usually 3, but it will be greater than 3
209            // small enough b, depending on the L/b ratio:
210            //     L/b == 1 => b < 2.37
211            //     L/b == 2 => b < 1.36
212            //     L/b == 3 => b < 1.00
213            //     L/b == 4 => b < 0.816
214            // 2017-10-01 pkienzle: moved low q approximation into Sdebye()
215            ans = Sdebye(square(q*Rg_short));
216        } else {  // q*b > max(1.9/Rg_short, 3)
217            //printf("branch D q=%g L=%g b=%g\n", q, L, b);
218            ans = a_short(q, L, b /*, p1short, p2short*/, q0short);
219        }
220    }
221
222    return ans;
223}
Note: See TracBrowser for help on using the repository browser.