Changes in / [09141ff:ba7302a] in sasmodels
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified explore/precision.py ¶
r487e695 r237c9cf 81 81 elif xrange == "linear": 82 82 lin_min, lin_max, lin_steps = 1, 1000, 2000 83 lin_min, lin_max, lin_steps = 0.001, 2, 2000 83 84 elif xrange == "log": 84 85 log_min, log_max, log_steps = -3, 5, 400 … … 321 322 np_function=lambda x: np.fmod(x, 2*np.pi), 322 323 ocl_function=make_ocl("return fmod(q, 2*M_PI);", "sas_fmod"), 324 ) 325 add_function( 326 name="debye", 327 mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4, 328 np_function=lambda x: 2*(np.expm1(-x**2) + x**2)/x**4, 329 ocl_function=make_ocl(""" 330 const double qsq = q*q; 331 if (qsq < 1.0) { // Pade approximation 332 const double x = qsq; 333 if (0) { // 0.36 single 334 // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 4}] 335 return (x*x/180. + 1.)/((1./30.*x + 1./3.)*x + 1); 336 } else if (0) { // 1.0 for single 337 // padeapproximant[2*exp[-x^2] + x^2-1)/x^4, {x, 0, 6}] 338 const double A1=1./24., A2=1./84, A3=-1./3360; 339 const double B1=3./8., B2=3./56., B3=1./336.; 340 return (((A3*x + A2)*x + A1)*x + 1.)/(((B3*x + B2)*x + B1)*x + 1.); 341 } else if (1) { // 1.0 for single, 0.25 for double 342 // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 343 const double A1=1./15., A2=1./60, A3=0., A4=1./75600.; 344 const double B1=2./5., B2=1./15., B3=1./180., B4=1./5040.; 345 return ((((A4*x + A3)*x + A2)*x + A1)*x + 1.) 346 /((((B4*x + B3)*x + B2)*x + B1)*x + 1.); 347 } else { // 1.0 for single, 0.5 for double 348 // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 349 const double A1=1./12., A2=2./99., A3=1./2640., A4=1./23760., A5=-1./1995840.; 350 const double B1=5./12., B2=5./66., B3=1./132., B4=1./2376., B5=1./95040.; 351 return (((((A5*x + A4)*x + A3)*x + A2)*x + A1)*x + 1.) 352 /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.); 353 } 354 } else if (qsq < 1.) { // Taylor series; 0.9 for single, 0.25 for double 355 const double x = qsq; 356 const double C0 = +1.; 357 const double C1 = -1./3.; 358 const double C2 = +1./12.; 359 const double C3 = -1./60.; 360 const double C4 = +1./360.; 361 const double C5 = -1./2520.; 362 const double C6 = +1./20160.; 363 const double C7 = -1./181440.; 364 //return ((((C5*x + C4)*x + C3)*x + C2)*x + C1)*x + C0; 365 //return (((((C6*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0; 366 return ((((((C7*x + C6)*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0; 367 } else { 368 return 2.*(expm1(-qsq) + qsq)/(qsq*qsq); 369 } 370 """, "sas_debye"), 323 371 ) 324 372 -
TabularUnified sasmodels/compare.py ¶
rb76191e rced5bd2 73 73 -res=0 sets the resolution width dQ/Q if calculating with resolution 74 74 -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0 75 -q=min:max alternative specification of qrange 75 76 -nq=128 sets the number of Q points in the data set 76 77 -1d*/-2d computes 1d or 2d data … … 768 769 'accuracy', 'is2d' and 'view' parsed from the command line. 769 770 """ 770 qm ax, nq, res =opts['qmax'], opts['nq'], opts['res']771 qmin, qmax, nq, res = opts['qmin'], opts['qmax'], opts['nq'], opts['res'] 771 772 if opts['is2d']: 772 773 q = np.linspace(-qmax, qmax, nq) # type: np.ndarray … … 777 778 else: 778 779 if opts['view'] == 'log' and not opts['zero']: 779 qmax = math.log10(qmax) 780 q = np.logspace(qmax-3, qmax, nq) 780 q = np.logspace(math.log10(qmin), math.log10(qmax), nq) 781 781 else: 782 q = np.linspace( 0.001*qmax, qmax, nq)782 q = np.linspace(qmin, qmax, nq) 783 783 if opts['zero']: 784 784 q = np.hstack((0, q)) … … 955 955 else: 956 956 err, errstr, errview = abs(relerr), "rel err", "log" 957 if (err == 0.).all(): 958 errview = 'linear' 957 959 if 0: # 95% cutoff 958 960 sorted = np.sort(err.flatten()) … … 960 962 err[err > cutoff] = cutoff 961 963 #err,errstr = base/comp,"ratio" 962 plot_theory(data, None, resid=err, view=errview, use_data=use_data) 963 if view == 'linear': 964 plt.xscale('linear') 964 plot_theory(data, None, resid=err, view=view, use_data=use_data) 965 plt.yscale(errview) 965 966 plt.title("max %s = %.3g"%(errstr, abs(err).max())) 966 967 #cbar_title = errstr if errview=="linear" else "log "+errstr … … 1001 1002 1002 1003 # Data generation 1003 'data=', 'noise=', 'res=', 1004 ' nq=', 'lowq', 'midq', 'highq', 'exq', 'zero',1004 'data=', 'noise=', 'res=', 'nq=', 'q=', 1005 'lowq', 'midq', 'highq', 'exq', 'zero', 1005 1006 '2d', '1d', 1006 1007 … … 1122 1123 'view' : 'log', 1123 1124 'is2d' : False, 1125 'qmin' : None, 1124 1126 'qmax' : 0.05, 1125 1127 'nq' : 128, … … 1159 1161 elif arg == '-zero': opts['zero'] = True 1160 1162 elif arg.startswith('-nq='): opts['nq'] = int(arg[4:]) 1163 elif arg.startswith('-q='): 1164 opts['qmin'], opts['qmax'] = [float(v) for v in arg[3:].split(':')] 1161 1165 elif arg.startswith('-res='): opts['res'] = float(arg[5:]) 1162 1166 elif arg.startswith('-noise='): opts['noise'] = float(arg[7:]) … … 1207 1211 1208 1212 # Create the computational engines 1213 if opts['qmin'] is None: 1214 opts['qmin'] = 0.001*opts['qmax'] 1209 1215 if opts['datafile'] is not None: 1210 1216 data = load_data(os.path.expanduser(opts['datafile'])) -
TabularUnified sasmodels/data.py ¶
r09141ff r09141ff 434 434 435 435 scale = data.x**4 if view == 'q4' else 1.0 436 xscale = yscale = 'linear' if view == 'linear' else 'log' 436 437 437 438 if use_data or use_theory: … … 466 467 plt.ylim(*limits) 467 468 468 plt.xscale('linear' if not some_present or non_positive_x 469 else view if view is not None 470 else 'log') 471 plt.yscale('linear' 472 if view == 'q4' or not some_present or not all_positive 473 else view if view is not None 474 else 'log') 469 470 xscale = ('linear' if not some_present or non_positive_x 471 else view if view is not None 472 else 'log') 473 yscale = ('linear' 474 if view == 'q4' or not some_present or not all_positive 475 else view if view is not None 476 else 'log') 477 plt.xscale(xscale) 475 478 plt.xlabel("$q$/A$^{-1}$") 479 plt.yscale(yscale) 476 480 plt.ylabel('$I(q)$') 477 481 title = ("data and model" if use_theory and use_data … … 501 505 plt.xlabel("$q$/A$^{-1}$") 502 506 plt.ylabel('residuals') 503 plt.xscale('linear')504 507 plt.title('(model - Iq)/dIq') 508 plt.xscale(xscale) 509 plt.yscale('linear') 505 510 506 511 -
TabularUnified sasmodels/models/flexible_cylinder.py ¶
r31df0c9 r2573fa1 90 90 length = 10**np.random.uniform(2, 6) 91 91 radius = 10**np.random.uniform(1, 3) 92 kuhn_length = 10**np.random.uniform(-2, -0.7)*length # at least 10 segments92 kuhn_length = 10**np.random.uniform(-2, 0)*length 93 93 pars = dict( 94 94 length=length, … … 100 100 tests = [ 101 101 # Accuracy tests based on content in test/utest_other_models.py 102 # Currently fails in OCL 103 # [{'length': 1000.0, 104 # 'kuhn_length': 100.0, 105 # 'radius': 20.0, 106 # 'sld': 1.0, 107 # 'sld_solvent': 6.3, 108 # 'background': 0.0001, 109 # }, 0.001, 3509.2187], 102 [{'length': 1000.0, # test T1 103 'kuhn_length': 100.0, 104 'radius': 20.0, 105 'sld': 1.0, 106 'sld_solvent': 6.3, 107 'background': 0.0001, 108 }, 0.001, 3509.2187], 110 109 111 110 # Additional tests with larger range of parameters 112 [{'length': 1000.0, 111 [{'length': 1000.0, # test T2 113 112 'kuhn_length': 100.0, 114 113 'radius': 20.0, … … 117 116 'background': 0.0001, 118 117 }, 1.0, 0.000595345], 119 [{'length': 10.0, 118 [{'length': 10.0, # test T3 120 119 'kuhn_length': 800.0, 121 120 'radius': 2.0, … … 124 123 'background': 0.001, 125 124 }, 0.1, 1.55228], 126 [{'length': 100.0, 125 [{'length': 100.0, # test T4 127 126 'kuhn_length': 800.0, 128 127 'radius': 50.0, … … 133 132 ] 134 133 134 # There are a few branches in the code that ought to have test values: 135 # 136 # For length > 4 * kuhn_length 137 # if length > 10 * kuhn_length then C is scaled by 3.06 (L/b)^(-0.44) 138 # q*kuhn_length <= 3.1 => Sexv_new 139 # dS/dQ < 0 has different behaviour from dS/dQ >= 0 140 # T2 q*kuhn_length > 3.1 => a_long 141 # 142 # For length <= 4 * kuhn_length 143 # q*kuhn_length <= max(1.9/Rg_short, 3.0) => Sdebye((q*Rg)^2) 144 # q*Rg < 0.5 uses Pade approx, q*Rg > 1.0 uses math lib 145 # T3,T4 q*kuhn_length > max(1.9/Rg_short, 3.0) => a_short 146 # 147 # Note that the transitions between branches may be abrupt. You can see a 148 # several percent change around length=10*kuhn_length and length=4*kuhn_length 149 # using the following: 150 # 151 # sascomp flexible_cylinder -calc=double -sets=10 length=10*kuhn_length,10.000001*kuhn_length 152 # sascomp flexible_cylinder -calc=double -sets=10 length=4*kuhn_length,4.000001*kuhn_length 153 # 154 # The transition between low q and high q around q*kuhn_length = 3 seems 155 # to be good to 4 digits or better. This was tested by computing the value 156 # on each branches near the transition point and reporting the relative error 157 # for kuhn lengths of 10, 100 and 1000 and a variety of length:kuhn_length 158 # ratios. -
TabularUnified sasmodels/models/lib/wrc_cyl.c ¶
r92ce163 r18a2bfc 2 2 Functions for WRC implementation of flexible cylinders 3 3 */ 4 double Sk_WR(double q, double L, double b); 5 6 7 static double 8 AlphaSquare(double x) 9 { 10 // Potentially faster. Needs proper benchmarking. 11 // add native_powr to kernel_template 12 //double t = native_powr( (1.0 + (x/3.12)*(x/3.12) + 13 // (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); 14 //return t; 15 16 return pow(1.0+square(x/3.12)+cube(x/8.67), 0.176/3.0); 17 } 18 19 // 20 static double 21 Rgsquarezero(double q, double L, double b) 22 { 23 const double r = b/L; 24 return (L*b/6.0) * (1.0 + r*(-1.5 + r*(1.5 + r*0.75*expm1(-2.0/r)))); 25 26 } 27 28 // 29 static double 30 Rgsquareshort(double q, double L, double b) 31 { 32 return AlphaSquare(L/b) * Rgsquarezero(q,L,b); 33 } 34 35 // 36 static double 37 Rgsquare(double q, double L, double b) 38 { 39 return AlphaSquare(L/b)*L*b/6.0; 40 } 41 42 static double 43 sech_WR(double x) 44 { 45 return(1/cosh(x)); 46 } 47 48 static double 49 a1long(double q, double L, double b, double p1, double p2, double q0) 50 { 51 double C; 52 53 if( L/b > 10.0) { 54 C = 3.06/pow((L/b),0.44); 4 5 static double 6 Rgsquare(double L, double b) 7 { 8 const double x = L/b; 9 // Use Horner's method to evaluate Pedersen eq 15: 10 // alpha^2 = [1.0 + (x/3.12)^2 + (x/8.67)^3] ^ (0.176/3) 11 const double alphasq = 12 pow(1.0 + x*x*(1.027284681130835e-01 + 1.534414548417740e-03*x), 13 5.866666666666667e-02); 14 return alphasq*L*b/6.0; 15 } 16 17 static double 18 Rgsquareshort(double L, double b) 19 { 20 const double r = b/L; // = 1/n_b in Pedersen ref. 21 return Rgsquare(L, b) * (1.0 + r*(-1.5 + r*(1.5 + r*0.75*expm1(-2.0/r)))); 22 } 23 24 static double 25 w_WR(double x) 26 { 27 // Pedersen eq. 16: 28 // w = [1 + tanh((x-C4)/C5)]/2 29 const double C4 = 1.523; 30 const double C5 = 0.1477; 31 return 0.5 + 0.5*tanh((x - C4)/C5); 32 } 33 34 static double 35 Sdebye(double qsq) 36 { 37 #if FLOAT_SIZE>4 38 # define DEBYE_CUTOFF 0.25 // 1e-15 error 39 #else 40 # define DEBYE_CUTOFF 1.0 // 4e-7 error 41 #endif 42 43 if (qsq < DEBYE_CUTOFF) { 44 const double x = qsq; 45 // mathematica: PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 46 const double A1=1./15., A2=1./60, A3=0., A4=1./75600.; 47 const double B1=2./5., B2=1./15., B3=1./180., B4=1./5040.; 48 return ((((A4*x + A3)*x + A2)*x + A1)*x + 1.) 49 / ((((B4*x + B3)*x + B2)*x + B1)*x + 1.); 55 50 } else { 56 C = 1.0;51 return 2.*(expm1(-qsq) + qsq)/(qsq*qsq); 57 52 } 58 53 } 54 55 static double 56 a_long(double q, double L, double b/*, double p1, double p2, double q0*/) 57 { 58 const double p1 = 4.12; 59 const double p2 = 4.42; 60 const double q0 = 3.1; 61 62 // Constants C1, ..., C5 derived from least squares fit (Pedersen, eq 13,16) 59 63 const double C1 = 1.22; 60 64 const double C2 = 0.4288; … … 64 68 const double miu = 0.585; 65 69 66 const double Rg2 = Rgsquare(q,L,b); 67 const double Rg22 = Rg2*Rg2; 68 const double Rg = sqrt(Rg2); 69 const double Rgb = Rg*q0/b; 70 71 const double b2 = b*b; 70 const double C = (L/b>10.0 ? 3.06*pow(L/b, -0.44) : 1.0); 71 //printf("branch B-%d q=%g L=%g b=%g\n", C==1.0, q, L, b); 72 const double r2 = Rgsquare(L,b); 73 const double r = sqrt(r2); 74 const double qr_b = q0*r/b; 75 const double qr_b_sq = qr_b*qr_b; 76 const double qr_b_4 = qr_b_sq*qr_b_sq; 77 const double qr_b_miu = pow(qr_b, -1.0/miu); 78 const double em1_qr_b_sq = expm1(-qr_b_sq); 79 const double sech2 = 1.0/square(cosh((qr_b-C4)/C5)); 80 const double w = w_WR(qr_b); 81 82 const double t1 = pow(q0, 1.0 + p1 + p2)/(b*(p1-p2)); 83 const double t2 = C/(15.0*L) * ( 84 + 14.0*b*b*em1_qr_b_sq/(q0*qr_b_sq) 85 + 2.0*q0*r2*exp(-qr_b_sq)*(11.0 + 7.0/qr_b_sq)); 86 const double t11 = ((C3*qr_b_miu + C2)*qr_b_miu + C1)*qr_b_miu; 87 const double t3 = r*sech2/(2.*C5)*t11; 88 const double t4 = r*(em1_qr_b_sq + qr_b_sq)*sech2 / (C5*qr_b_4); 89 const double t5 = -4.0*r*qr_b*em1_qr_b_sq/qr_b_4 * (1.0 - w); 90 const double t10 = 2.0*(em1_qr_b_sq + qr_b_sq)/qr_b_4 * (1.0 - w); //=Sdebye*(1-w) 91 const double t6 = 4.0*b/q0 * t10; 92 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); 93 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; 94 const double t12 = b*b*M_PI/(L*q0*q0) + t2 + t3 - t4 + t5 - t6 + t7*w; 95 const double t13 = -b*M_PI/(L*q0) + t9 + t10 + t11*w; 96 97 const double a1 = pow(q0,p1)*t13 - t1*pow(q0,-p2)*(t12 + b*p1/q0*t13); 98 const double a2 = t1*pow(q0,-p1)*(t12 + b*p1/q0*t13); 99 100 const double ans = a1*pow(q*b, -p1) + a2*pow(q*b, -p2) + M_PI/(q*L); 101 return ans; 102 } 103 104 static double 105 _short(double r2, double exp_qr_b, double L, double b, double p1short, 106 double p2short, double q0) 107 { 108 const double qr2 = q0*q0 * r2; 72 109 const double b3 = b*b*b; 73 const double b4 = b3*b; 74 const double q02 = q0*q0; 75 const double q03 = q0*q0*q0; 76 const double q04 = q03*q0; 77 const double q05 = q04*q0; 78 79 const double Rg02 = Rg2*q02; 80 81 const double t1 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2))) * 82 ((11.0/15.0 + (7.0*b2)/(15.0*Rg02))) + 83 (7.0*b2)/(15.0*Rg02)))); 84 85 const double t2 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 86 Rg02/b2))*((1.0 + 0.5*(((-1.0) - 87 tanh((-C4 + Rgb/C5))))))); 88 89 const double t3 = ((C3*pow(Rgb,((-3.0)/miu)) + 90 C2*pow(Rgb,((-2.0)/miu)) + 91 C1*pow(Rgb,((-1.0)/miu)))); 92 93 const double t4 = ((1.0 + tanh(((-C4) + Rgb)/C5))); 94 95 const double t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - 96 b*p2*pow(q0,((-1.0) - p1 - p2)))); 97 98 const double t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) + 99 (14.0*b3*pow((double)M_E,(-(Rg02/b2))))/(15.0*q03*Rg2) + 100 (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*((11.0/15.0 + 101 (7.0*b2)/(15.0*Rg02)))*Rg2)/b))); 102 103 const double t7 = (Rg*((C3*pow(((Rgb)),((-3.0)/miu)) + 104 C2*pow(((Rgb)),((-2.0)/miu)) + 105 C1*pow(((Rgb)),((-1.0)/miu))))*pow(sech_WR(((-C4) + 106 Rgb)/C5),2)); 107 108 const double t8 = (b4*Rg*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 109 Rg02/b2))*pow(sech_WR(((-C4) + Rgb)/C5),2)); 110 111 const double t9 = (2.0*b4*(((2.0*q0*Rg2)/b - 112 (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))*((1.0 + 0.5*(((-1.0) - 113 tanh(((-C4) + Rgb)/C5)))))); 114 115 const double t10 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 116 Rg02/b2))*((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 117 Rgb)/C5)))))); 118 119 const double t11 = (((-((3.0*C3*Rg*pow(((Rgb)),((-1.0) - 120 3.0/miu)))/miu)) - (2.0*C2*Rg*pow(((Rgb)),((-1.0) - 121 2.0/miu)))/miu - (C1*Rg*pow(((Rgb)),((-1.0) - 122 1.0/miu)))/miu)); 123 124 const double t12 = ((1.0 + tanh(((-C4) + Rgb)/C5))); 125 126 const double t13 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2)))*((11.0/15.0 + 127 (7.0*b2)/(15.0*q02* Rg2))) + 128 (7.0*b2)/(15.0*Rg02)))); 129 130 const double t14 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 131 Rg02/b2))*((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 132 Rgb)/C5)))))); 133 134 const double t15 = ((C3*pow(((Rgb)),((-3.0)/miu)) + 135 C2*pow(((Rgb)),((-2.0)/miu)) + 136 C1*pow(((Rgb)),((-1.0)/miu)))); 137 138 139 double yy = (pow(q0,p1)*(((-((b*M_PI)/(L*q0))) +t1/L +t2/(q04*Rg22) + 140 0.5*t3*t4)) + (t5*((pow(q0,(p1 - p2))* 141 (((-pow(q0,(-p1)))*(((b2*M_PI)/(L*q02) +t6/L +t7/(2.0*C5) - 142 t8/(C5*q04*Rg22) + t9/(q04*Rg22) -t10/(q05*Rg22) + 0.5*t11*t12)) - 143 b*p1*pow(q0,((-1.0) - p1))*(((-((b*M_PI)/(L*q0))) + t13/L + 144 t14/(q04*Rg22) + 0.5*t15*((1.0 + tanh(((-C4) + 145 Rgb)/C5))))))))))); 146 147 return (yy); 148 } 149 150 static double 151 a2long(double q, double L, double b, double p1, double p2, double q0) 152 { 153 double C; 154 155 if( L/b > 10.0) { 156 C = 3.06/pow((L/b),0.44); 157 } else { 158 C = 1.0; 159 } 160 161 const double C1 = 1.22; 162 const double C2 = 0.4288; 163 const double C3 = -1.651; 164 const double C4 = 1.523; 165 const double C5 = 0.1477; 166 const double miu = 0.585; 167 168 const double Rg2 = Rgsquare(q,L,b); 169 const double Rg22 = Rg2*Rg2; 170 const double b2 = b*b; 171 const double b3 = b*b*b; 172 const double b4 = b3*b; 173 const double q02 = q0*q0; 174 const double q03 = q0*q0*q0; 175 const double q04 = q03*q0; 176 const double q05 = q04*q0; 177 const double Rg = sqrt(Rg2); 178 const double Rgb = Rg*q0/b; 179 const double Rg02 = Rg2*q02; 180 181 const double t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - 182 b*p2*pow(q0,((-1.0) - p1 - p2)) )); 183 184 const double t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + 185 (14.0*b3*pow((double)M_E,(-(Rg02/b2))))/(15.0*q03*Rg2) + 186 (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*((11.0/15.0 + 187 (7*b2)/(15.0*Rg02)))*Rg2)/b)))/L; 188 189 const double t3 = (Rg*((C3*pow(((Rgb)),((-3.0)/miu)) + 190 C2*pow(((Rgb)),((-2.0)/miu)) + 191 C1*pow(((Rgb)),((-1.0)/miu))))* 192 pow(sech_WR(((-C4) +Rgb)/C5),2.0))/(2.0*C5); 193 194 const double t4 = (b4*Rg*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 195 Rg02/b2))*pow(sech_WR(((-C4) + 196 Rgb)/C5),2))/(C5*q04*Rg22); 197 198 const double t5 = (2.0*b4*(((2.0*q0*Rg2)/b - 199 (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))* 200 ((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 201 Rgb)/C5))))))/(q04*Rg22); 202 203 const double t6 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 204 Rg02/b2))*((1.0 + 0.5*(((-1) - tanh(((-C4) + 205 Rgb)/C5))))))/(q05*Rg22); 206 207 const double t7 = (((-((3.0*C3*Rg*pow(((Rgb)),((-1.0) - 208 3.0/miu)))/miu)) - (2.0*C2*Rg*pow(((Rgb)), 209 ((-1.0) - 2.0/miu)))/miu - (C1*Rg*pow(((Rgb)), 210 ((-1.0) - 1.0/miu)))/miu)); 211 212 const double t8 = ((1.0 + tanh(((-C4) + Rgb)/C5))); 213 214 const double t9 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2)))*((11.0/15.0 + 215 (7.0*b2)/(15*Rg02))) + (7.0*b2)/(15.0*Rg02))))/L; 216 217 const double t10 = (2.0*b4*(((-1) + pow((double)M_E,(-(Rg02/b2))) + 218 Rg02/b2))*((1.0 + 0.5*(((-1) - tanh(((-C4) + 219 Rgb)/C5))))))/(q04*Rg22); 220 221 const double yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*M_PI)/(L*q02) + 222 t2 + t3 - t4 + t5 - t6 + 0.5*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))* 223 (((-((b*M_PI)/(L*q0))) + t9 + t10 + 224 0.5*((C3*pow(((Rgb)),((-3.0)/miu)) + 225 C2*pow(((Rgb)),((-2.0)/miu)) + 226 C1*pow(((Rgb)),((-1.0)/miu))))* 227 ((1.0 + tanh(((-C4) + Rgb)/C5)))))))))); 228 229 return (yy); 230 } 231 232 // 233 static double 234 a_short(double q, double L, double b, double p1short, 235 double p2short, double factor, double pdiff, double q0) 236 { 237 const double Rg2_sh = Rgsquareshort(q,L,b); 238 const double Rg2_sh2 = Rg2_sh*Rg2_sh; 239 const double b3 = b*b*b; 240 const double t1 = ((q0*q0*Rg2_sh)/(b*b)); 241 const double Et1 = exp(t1); 242 const double Emt1 = 1.0/Et1; 243 const double q02 = q0*q0; 244 const double q0p = pow(q0,(-4.0 + p1short) ); 245 246 double yy = ((factor/(L*(pdiff)*Rg2_sh2)* 247 ((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short + 248 2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 249 2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*M_PI*q02*q0*Rg2_sh2 + 250 Et1*p2short*M_PI*q02*q0*Rg2_sh2)))))); 251 252 return(yy); 253 } 254 static double 255 a1short(double q, double L, double b, double p1short, double p2short, double q0) 256 { 257 258 double factor = 1.0; 259 return a_short(q, L, b, p1short, p2short, factor, p1short - p2short, q0); 260 } 261 262 static double 263 a2short(double q, double L, double b, double p1short, double p2short, double q0) 264 { 265 double factor = -1.0; 266 return a_short(q, L, b, p2short, p1short, factor, p1short-p2short, q0); 267 } 268 269 //WR named this w (too generic) 270 static double 271 w_WR(double x) 272 { 273 return 0.5*(1 + tanh((x - 1.523)/0.1477)); 274 } 275 276 // 277 static double 278 u1(double q, double L, double b) 279 { 280 return Rgsquareshort(q,L,b)*q*q; 281 } 282 283 static double 284 u_WR(double q, double L, double b) 285 { 286 return Rgsquare(q,L,b)*q*q; 287 } 288 289 static double 290 Sdebye_kernel(double arg) 291 { 292 // ORIGINAL 293 double result = 2.0*(exp(-arg) + arg -1.0)/(arg*arg); 294 295 // CONVERSION 1 from http://herbie.uwplse.org/ 296 // 297 // exhibits discontinuity - needs more investigation 298 //double a1 = 1.0/6.0; 299 //double a2 = 1.0/72.0; 300 //double a3 = 1.0/24.0; 301 //double result = pow((1.0 - a1*arg - (a2+a3)*arg*arg), 2); 302 303 return result; 304 } 305 static double 306 Sdebye(double q, double L, double b) 307 { 308 double arg = u_WR(q,L,b); 309 return Sdebye_kernel(arg); 310 } 311 312 // 313 static double 314 Sdebye1(double q, double L, double b) 315 { 316 double arg = u1(q,L,b); 317 return Sdebye_kernel(arg); 318 319 } 320 321 // 110 const double q0p = pow(q0, -4.0 + p1short); 111 112 double yy = 1.0/(L*r2*r2) * b/exp_qr_b*q0p 113 * (8.0*b3*L 114 - 8.0*b3*exp_qr_b*L 115 + 2.0*b3*exp_qr_b*L*p2short 116 - 2.0*b*exp_qr_b*L*p2short*qr2 117 + 4.0*b*exp_qr_b*L*qr2 118 - 2.0*b3*L*p2short 119 + 4.0*b*L*qr2 120 - M_PI*exp_qr_b*qr2*q0*r2 121 + M_PI*exp_qr_b*p2short*qr2*q0*r2); 122 123 return yy; 124 } 125 static double 126 a_short(double qp, double L, double b 127 /*double p1short, double p2short*/, double q0) 128 { 129 const double p1short = 5.36; 130 const double p2short = 5.62; 131 132 const double r2 = Rgsquareshort(L,b); 133 const double exp_qr_b = exp(r2*square(q0/b)); 134 const double pdiff = p1short - p2short; 135 const double a1 = _short(r2,exp_qr_b,L,b,p1short,p2short,q0)/pdiff; 136 const double a2= -_short(r2,exp_qr_b,L,b,p2short,p1short,q0)/pdiff; 137 const double ans = a1*pow(qp*b, -p1short) + a2*pow(qp*b, -p2short) + M_PI/(qp*L); 138 return ans; 139 } 140 141 322 142 static double 323 143 Sexv(double q, double L, double b) 324 144 { 325 145 // Pedersen eq 13, corrected by Chen eq A.5, swapping w and 1-w 326 146 const double C1=1.22; 327 147 const double C2=0.4288; 328 148 const double C3=-1.651; 329 149 const double miu = 0.585; 330 const double qRg = q*sqrt(Rgsquare(q,L,b)); 331 const double x = pow(qRg, -1.0/miu); 332 333 double yy = (1.0 - w_WR(qRg))*Sdebye(q,L,b) + 334 w_WR(qRg)*x*(C1 + x*(C2 + x*C3)); 335 return (yy); 336 } 337 338 339 static double 340 Sexvnew(double q, double L, double b) 341 { 342 double yy; 343 344 const double C1 =1.22; 345 const double C2 =0.4288; 346 const double C3 =-1.651; 347 const double miu = 0.585; 150 const double qr = q*sqrt(Rgsquare(L,b)); 151 const double qr_miu = pow(qr, -1.0/miu); 152 const double w = w_WR(qr); 153 const double t10 = Sdebye(qr*qr)*(1.0 - w); 154 const double t11 = ((C3*qr_miu + C2)*qr_miu + C1)*qr_miu; 155 156 return t10 + w*t11; 157 } 158 159 // Modified by Yun on Oct. 15, 160 static double 161 Sexv_new(double q, double L, double b) 162 { 163 const double qr = q*sqrt(Rgsquare(L,b)); 164 const double qr2 = qr*qr; 165 const double C = (L/b > 10.0) ? 3.06*pow(L/b, -0.44) : 1.0; 166 const double t9 = C*b/L * (4.0 - exp(-qr2) * (11.0 + 7.0/qr2) + 7.0/qr2)/15.0; 167 168 const double Sexv_orig = Sexv(q, L, b); 169 170 // calculating the derivative to decide on the correction (cutoff) term? 171 // Note: this is modified from WRs original code 348 172 const double del=1.05; 349 const double qRg = q*sqrt(Rgsquare(q,L,b)); 350 const double x = pow(qRg, -1.0/miu); 351 352 353 //calculating the derivative to decide on the corection (cutoff) term? 354 // I have modified this from WRs original code 355 const double qdel = (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q); 356 const double C_star2 =(qdel >= 0.0) ? 0.0 : 1.0; 357 358 yy = (1.0 - w_WR(qRg))*Sdebye(q,L,b) + 359 C_star2*w_WR(qRg)* 360 x*(C1 + x*(C2 + x*C3)); 361 362 return (yy); 363 } 364 365 double Sk_WR(double q, double L, double b) 366 { 367 // 368 const double p1 = 4.12; 369 const double p2 = 4.42; 370 const double p1short = 5.36; 371 const double p2short = 5.62; 372 const double q0 = 3.1; 373 374 double q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0); 375 double Sexvmodify, ans; 376 377 const double C =(L/b > 10.0) ? 3.06/pow((L/b),0.44) : 1.0; 378 379 if( L > 4*b ) { // Longer Chains 380 if (q*b <= 3.1) { //Modified by Yun on Oct. 15, 381 Sexvmodify = Sexvnew(q, L, b); 382 ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) - 383 (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L); 384 385 } else { //q(i)*b > 3.1 386 ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) + 387 a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + M_PI/(q*L); 388 } 389 } else { //L <= 4*b Shorter Chains 390 if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) { 391 if (q*b<=0.01) { 392 ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0; 393 } else { 394 ans = Sdebye1(q,L,b); 395 } 396 } else { //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3) 397 ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) + 398 a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) + 399 M_PI/(q*L); 400 } 173 const double qdel = (Sexv(q*del,L,b) - Sexv_orig)/(q*(del - 1.0)); 174 175 if (qdel < 0) { 176 //printf("branch A1-%d q=%g L=%g b=%g\n", C==1.0, q, L, b); 177 return t9 + Sexv_orig; 178 } else { 179 //printf("branch A2-%d q=%g L=%g b=%g\n", C==1.0, q, L, b); 180 const double w = w_WR(qr); 181 const double t10 = Sdebye(qr*qr)*(1.0 - w); 182 return t9 + t10; 401 183 } 402 403 return(ans); 404 } 184 } 185 186 187 static double 188 Sk_WR(double q, double L, double b) 189 { 190 const double Rg_short = sqrt(Rgsquareshort(L, b)); 191 double q0short = fmax(1.9/Rg_short, 3.0); 192 double ans; 193 194 195 if( L > 4*b ) { // L > 4*b : Longer Chains 196 if (q*b <= 3.1) { 197 ans = Sexv_new(q, L, b); 198 } else { //q(i)*b > 3.1 199 ans = a_long(q, L, b /*, p1, p2, q0*/); 200 } 201 } else { // L <= 4*b : Shorter Chains 202 if (q*b <= q0short) { // q*b <= fmax(1.9/Rg_short, 3) 203 //printf("branch C-%d q=%g L=%g b=%g\n", square(q*Rg_short)<DEBYE_CUTOFF, q, L, b); 204 // Note that q0short is usually 3, but it will be greater than 3 205 // small enough b, depending on the L/b ratio: 206 // L/b == 1 => b < 2.37 207 // L/b == 2 => b < 1.36 208 // L/b == 3 => b < 1.00 209 // L/b == 4 => b < 0.816 210 // 2017-10-01 pkienzle: moved low q approximation into Sdebye() 211 ans = Sdebye(square(q*Rg_short)); 212 } else { // q*b > max(1.9/Rg_short, 3) 213 //printf("branch D q=%g L=%g b=%g\n", q, L, b); 214 ans = a_short(q, L, b /*, p1short, p2short*/, q0short); 215 } 216 } 217 218 return ans; 219 }
Note: See TracChangeset
for help on using the changeset viewer.