source: sasmodels/sasmodels/models/lib/wrc_cyl.c @ 13ed84c

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 13ed84c was 13ed84c, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

set single=False on all models that fail the single precision tests
  • Property mode set to 100644
File size: 11.4 KB
Line 
1/*
2    Functions for WRC implementation of flexible cylinders
3*/
4static double
5AlphaSquare(double x)
6{
7    return pow( (1.0 + (x/3.12)*(x/3.12) +
8         (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) );
9}
10
11//
12static double
13Rgsquarezero(double q, double L, double b)
14{
15    return (L*b/6.0) * (1.0 - 1.5*(b/L) + 1.5*pow((b/L),2) -
16         0.75*pow((b/L),3)*(1.0 - exp(-2.0*(L/b))));
17}
18
19//
20static double
21Rgsquareshort(double q, double L, double b)
22{
23    return AlphaSquare(L/b) * Rgsquarezero(q,L,b);
24}
25
26//
27static double
28Rgsquare(double q, double L, double b)
29{
30    return AlphaSquare(L/b)*L*b/6.0;
31}
32
33static double
34sech_WR(double x)
35{
36    return(1/cosh(x));
37}
38
39static double
40a1long(double q, double L, double b, double p1, double p2, double q0)
41{
42    double yy,C,C1,C2,C3,C4,C5,miu,Rg2;
43    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15;
44    double E,pi;
45    double b2,b3,b4,q02,q03,q04,q05,Rg22;
46
47    pi = 4.0*atan(1.0);
48    E = 2.718281828459045091;
49
50    if( L/b > 10.0) {
51        C = 3.06/pow((L/b),0.44);
52    } else {
53        C = 1.0;
54    }
55
56    C1 = 1.22;
57    C2 = 0.4288;
58    C3 = -1.651;
59    C4 = 1.523;
60    C5 = 0.1477;
61    miu = 0.585;
62
63    Rg2 = Rgsquare(q,L,b);
64    Rg22 = Rg2*Rg2;
65    b2 = b*b;
66    b3 = b*b*b;
67    b4 = b3*b;
68    q02 = q0*q0;
69    q03 = q0*q0*q0;
70    q04 = q03*q0;
71    q05 = q04*q0;
72
73    t1 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 +
74        (7.0*b2)/(15.0*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))));
75
76    t2 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
77        (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) -
78        tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
79
80    t3 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) +
81        C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) +
82        C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))));
83
84    t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
85
86    t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) -
87        b*p2*pow(q0,((-1.0) - p1 - p2))));
88
89    t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) +
90        (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) +
91        (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 +
92        (7.0*b2)/(15.0*q02*Rg2)))*Rg2)/b)));
93
94    t7 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) +
95        C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) +
96        C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) +
97        (sqrt(Rg2)*q0)/b)/C5),2));
98
99    t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
100        (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));
101
102    t9 = (2.0*b4*(((2.0*q0*Rg2)/b -
103        (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) -
104        tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
105
106    t10 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
107        (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) +
108        (sqrt(Rg2)*q0)/b)/C5))))));
109
110    t11 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) -
111        3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) -
112        2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) -
113        1.0/miu)))/miu));
114
115    t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
116
117    t13 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 +
118        (7.0*b2)/(15.0*q02* Rg2))) +
119        (7.0*b2)/(15.0*q02*Rg2))));
120
121    t14 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
122        (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) +
123        (sqrt(Rg2)*q0)/b)/C5))))));
124
125    t15 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) +
126        C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) +
127        C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))));
128
129
130    yy = (pow(q0,p1)*(((-((b*pi)/(L*q0))) +t1/L +t2/(q04*Rg22) +
131        1.0/2.0*t3*t4)) + (t5*((pow(q0,(p1 - p2))*
132        (((-pow(q0,(-p1)))*(((b2*pi)/(L*q02) +t6/L +t7/(2.0*C5) -
133        t8/(C5*q04*Rg22) + t9/(q04*Rg22) -t10/(q05*Rg22) + 1.0/2.0*t11*t12)) -
134        b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) + t13/L +
135        t14/(q04*Rg22) + 1.0/2.0*t15*((1.0 + tanh(((-C4) +
136        (sqrt(Rg2)*q0)/b)/C5)))))))))));
137
138    return (yy);
139}
140
141static double
142a2long(double q, double L, double b, double p1, double p2, double q0)
143{
144    double yy,C1,C2,C3,C4,C5,miu,C,Rg2;
145    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi;
146    double E,b2,b3,b4,q02,q03,q04,q05,Rg22;
147
148    pi = 4.0*atan(1.0);
149    E = 2.718281828459045091;
150    if( L/b > 10.0) {
151        C = 3.06/pow((L/b),0.44);
152    } else {
153        C = 1.0;
154    }
155
156    C1 = 1.22;
157    C2 = 0.4288;
158    C3 = -1.651;
159    C4 = 1.523;
160    C5 = 0.1477;
161    miu = 0.585;
162
163    Rg2 = Rgsquare(q,L,b);
164    Rg22 = Rg2*Rg2;
165    b2 = b*b;
166    b3 = b*b*b;
167    b4 = b3*b;
168    q02 = q0*q0;
169    q03 = q0*q0*q0;
170    q04 = q03*q0;
171    q05 = q04*q0;
172
173    t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) -
174        b*p2*pow(q0,((-1.0) - p1 - p2)) ));
175
176    t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) +
177        (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) +
178        (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 +
179        (7*b2)/(15.0*q02*Rg2)))*Rg2)/b)))/L;
180
181    t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) +
182        C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) +
183        C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*
184        pow(sech_WR(((-C4) +(sqrt(Rg2)*q0)/b)/C5),2.0))/(2.0*C5);
185
186    t4 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
187        (q02*Rg2)/b2))*pow(sech_WR(((-C4) +
188        (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22);
189
190    t5 = (2.0*b4*(((2.0*q0*Rg2)/b -
191        (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*
192        ((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) +
193        (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);
194
195    t6 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
196        (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) +
197        (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22);
198
199    t7 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) -
200        3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),
201        ((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),
202        ((-1.0) - 1.0/miu)))/miu));
203
204    t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
205
206    t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 +
207        (7.0    *b2)/(15*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))))/L;
208
209    t10 = (2.0*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) +
210        (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) +
211        (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);
212
213    yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) +
214        t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*
215        (((-((b*pi)/(L*q0))) + t9 + t10 +
216        1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) +
217        C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) +
218        C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*
219        ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))));
220
221    return (yy);
222}
223
224//
225static double
226a1short(double q, double L, double b, double p1short, double p2short, double q0)
227{
228    double yy,Rg2_sh;
229    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3;
230    double pi;
231
232    E = 2.718281828459045091;
233    pi = 4.0*atan(1.0);
234    Rg2_sh = Rgsquareshort(q,L,b);
235    Rg2_sh2 = Rg2_sh*Rg2_sh;
236    b3 = b*b*b;
237    t1 = ((q0*q0*Rg2_sh)/(b*b));
238    Et1 = pow(E,t1);
239    Emt1 =pow(E,-t1);
240    q02 = q0*q0;
241    q0p = pow(q0,(-4.0 + p1short) );
242
243    yy = ((1.0/(L*((p1short - p2short))*Rg2_sh2)*
244        ((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short +
245        2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh -
246        2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 +
247        Et1*p2short*pi*q02*q0*Rg2_sh2))))));
248
249    return(yy);
250}
251
252static double
253a2short(double q, double L, double b, double p1short, double p2short, double q0)
254{
255    double yy,Rg2_sh;
256    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p;
257    double pi;
258
259    E = 2.718281828459045091;
260    pi = 4.0*atan(1.0);
261    Rg2_sh = Rgsquareshort(q,L,b);
262    Rg2_sh2 = Rg2_sh*Rg2_sh;
263    t1 = ((q0*q0*Rg2_sh)/(b*b));
264    Et1 = pow(E,t1);
265    Emt1 =pow(E,-t1);
266    q02 = q0*q0;
267    q0p = pow(q0,(-4.0 + p2short) );
268
269    //E is the number e
270    yy = ((-(1.0/(L*((p1short - p2short))*Rg2_sh2)*
271        ((b*Emt1*q0p*((8.0*b*b*b*L - 8.0*b*b*b*Et1*L - 2.0*b*b*b*L*p1short +
272        2.0*b*b*b*Et1*L*p1short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh -
273        2.0*b*Et1*L*p1short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 +
274        Et1*p1short*pi*q02*q0*Rg2_sh2)))))));
275
276    return (yy);
277}
278
279//WR named this w (too generic)
280static double
281w_WR(double x)
282{
283    return 0.5*(1 + tanh((x - 1.523)/0.1477));
284}
285
286//
287static double
288u1(double q, double L, double b)
289{
290    return Rgsquareshort(q,L,b)*q*q;
291}
292
293static double
294u_WR(double q, double L, double b)
295{
296    return Rgsquare(q,L,b)*q*q;
297}
298
299static double
300Sdebye_kernel(double arg)
301{
302    // ORIGINAL
303    double result = 2.0*(exp(-arg) + arg -1.0)/(pow((arg),2));
304
305    // CONVERSION 1 from http://herbie.uwplse.org/
306    //
307    // exhibits discontinuity - needs more investigation
308    //double a1 = 1.0/6.0;
309    //double a2 = 1.0/72.0;
310    //double a3 = 1.0/24.0;
311    //double result = pow((1.0 - a1*arg - (a2+a3)*arg*arg), 2);
312
313    return result;
314}
315static double
316Sdebye(double q, double L, double b)
317{
318    double arg = u_WR(q,L,b);
319    return Sdebye_kernel(arg);
320}
321
322//
323static double
324Sdebye1(double q, double L, double b)
325{
326    double arg = u1(q,L,b);
327    return Sdebye_kernel(arg);
328
329}
330
331//
332static double
333Sexv(double q, double L, double b)
334{
335    double yy,C1,C2,C3,miu,Rg2;
336    C1=1.22;
337    C2=0.4288;
338    C3=-1.651;
339    miu = 0.585;
340
341    Rg2 = Rgsquare(q,L,b);
342
343    yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) +
344        w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) +
345        C2*pow((q*sqrt(Rg2)),(-2.0/miu)) +
346        C3*pow((q*sqrt(Rg2)),(-3.0/miu)));
347
348    return (yy);
349}
350
351
352static double
353Sexvnew(double q, double L, double b)
354{
355    double yy,C1,C2,C3,miu;
356    double del=1.05,C_star2,Rg2;
357
358    C1=1.22;
359    C2=0.4288;
360    C3=-1.651;
361    miu = 0.585;
362
363    //calculating the derivative to decide on the corection (cutoff) term?
364    // I have modified this from WRs original code
365
366    if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) {
367        C_star2 = 0.0;
368    } else {
369        C_star2 = 1.0;
370    }
371
372    Rg2 = Rgsquare(q,L,b);
373
374    yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) +
375        C_star2*w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) +
376        C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu)));
377
378    return (yy);
379}
380
381double Sk_WR(double q, double L, double b);
382double Sk_WR(double q, double L, double b)
383{
384  //
385  double p1,p2,p1short,p2short,q0;
386  double C,ans,q0short,Sexvmodify,pi;
387
388  pi = 4.0*atan(1.0);
389
390  p1 = 4.12;
391  p2 = 4.42;
392  p1short = 5.36;
393  p2short = 5.62;
394  q0 = 3.1;
395
396  q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0);
397
398
399  if(L/b > 10.0) {
400    C = 3.06/pow((L/b),0.44);
401  } else {
402    C = 1.0;
403  }
404  //
405
406  //
407  if( L > 4*b ) { // Longer Chains
408    if (q*b <= 3.1) {   //Modified by Yun on Oct. 15,
409
410      Sexvmodify = Sexvnew(q, L, b);
411
412      ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) -
413            (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L);
414
415    } else { //q(i)*b > 3.1
416      ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) +
417            a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + pi/(q*L);
418    }
419  } else { //L <= 4*b Shorter Chains
420    if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) {
421      if (q*b<=0.01) {
422        ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0;
423      } else {
424        ans = Sdebye1(q,L,b);
425      }
426    } else {  //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3)
427      ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) +
428            a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) +
429            pi/(q*L);
430    }
431  }
432
433  return(ans);
434}
Note: See TracBrowser for help on using the repository browser.