Changes in / [ba7302a:09141ff] in sasmodels


Ignore:
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • explore/precision.py

    r237c9cf r487e695  
    8181        elif xrange == "linear": 
    8282            lin_min, lin_max, lin_steps = 1, 1000, 2000 
    83             lin_min, lin_max, lin_steps = 0.001, 2, 2000 
    8483        elif xrange == "log": 
    8584            log_min, log_max, log_steps = -3, 5, 400 
    … …  
    322321    np_function=lambda x: np.fmod(x, 2*np.pi), 
    323322    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"), 
    371323) 
    372324 

  • sasmodels/compare.py

    rced5bd2 rb76191e  
    7373    -res=0 sets the resolution width dQ/Q if calculating with resolution 
    7474    -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 
    7675    -nq=128 sets the number of Q points in the data set 
    7776    -1d*/-2d computes 1d or 2d data 
    … …  
    769768    'accuracy', 'is2d' and 'view' parsed from the command line. 
    770769    """ 
    771     qmin, qmax, nq, res = opts['qmin'], opts['qmax'], opts['nq'], opts['res'] 
     770    qmax, nq, res = opts['qmax'], opts['nq'], opts['res'] 
    772771    if opts['is2d']: 
    773772        q = np.linspace(-qmax, qmax, nq)  # type: np.ndarray 
    … …  
    778777    else: 
    779778        if opts['view'] == 'log' and not opts['zero']: 
    780             q = np.logspace(math.log10(qmin), math.log10(qmax), nq) 
     779            qmax = math.log10(qmax) 
     780            q = np.logspace(qmax-3, qmax, nq) 
    781781        else: 
    782             q = np.linspace(qmin, qmax, nq) 
     782            q = np.linspace(0.001*qmax, qmax, nq) 
    783783        if opts['zero']: 
    784784            q = np.hstack((0, q)) 
    … …  
    955955        else: 
    956956            err, errstr, errview = abs(relerr), "rel err", "log" 
    957             if (err == 0.).all(): 
    958                 errview = 'linear' 
    959957        if 0:  # 95% cutoff 
    960958            sorted = np.sort(err.flatten()) 
    … …  
    962960            err[err > cutoff] = cutoff 
    963961        #err,errstr = base/comp,"ratio" 
    964         plot_theory(data, None, resid=err, view=view, use_data=use_data) 
    965         plt.yscale(errview) 
     962        plot_theory(data, None, resid=err, view=errview, use_data=use_data) 
     963        if view == 'linear': 
     964            plt.xscale('linear') 
    966965        plt.title("max %s = %.3g"%(errstr, abs(err).max())) 
    967966        #cbar_title = errstr if errview=="linear" else "log "+errstr 
    … …  
    10021001 
    10031002    # Data generation 
    1004     'data=', 'noise=', 'res=', 'nq=', 'q=', 
    1005     'lowq', 'midq', 'highq', 'exq', 'zero', 
     1003    'data=', 'noise=', 'res=', 
     1004    'nq=', 'lowq', 'midq', 'highq', 'exq', 'zero', 
    10061005    '2d', '1d', 
    10071006 
    … …  
    11231122        'view'      : 'log', 
    11241123        'is2d'      : False, 
    1125         'qmin'      : None, 
    11261124        'qmax'      : 0.05, 
    11271125        'nq'        : 128, 
    … …  
    11611159        elif arg == '-zero':    opts['zero'] = True 
    11621160        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(':')] 
    11651161        elif arg.startswith('-res='):      opts['res'] = float(arg[5:]) 
    11661162        elif arg.startswith('-noise='):    opts['noise'] = float(arg[7:]) 
    … …  
    12111207 
    12121208    # Create the computational engines 
    1213     if opts['qmin'] is None: 
    1214         opts['qmin'] = 0.001*opts['qmax'] 
    12151209    if opts['datafile'] is not None: 
    12161210        data = load_data(os.path.expanduser(opts['datafile'])) 

  • sasmodels/data.py

    r09141ff r09141ff  
    434434 
    435435    scale = data.x**4 if view == 'q4' else 1.0 
    436     xscale = yscale = 'linear' if view == 'linear' else 'log' 
    437436 
    438437    if use_data or use_theory: 
    … …  
    467466            plt.ylim(*limits) 
    468467 
    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) 
     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') 
    478475        plt.xlabel("$q$/A$^{-1}$") 
    479         plt.yscale(yscale) 
    480476        plt.ylabel('$I(q)$') 
    481477        title = ("data and model" if use_theory and use_data 
    … …  
    505501        plt.xlabel("$q$/A$^{-1}$") 
    506502        plt.ylabel('residuals') 
     503        plt.xscale('linear') 
    507504        plt.title('(model - Iq)/dIq') 
    508         plt.xscale(xscale) 
    509         plt.yscale('linear') 
    510505 
    511506 

  • sasmodels/models/flexible_cylinder.py

    r2573fa1 r31df0c9  
    9090    length = 10**np.random.uniform(2, 6) 
    9191    radius = 10**np.random.uniform(1, 3) 
    92     kuhn_length = 10**np.random.uniform(-2, 0)*length 
     92    kuhn_length = 10**np.random.uniform(-2, -0.7)*length  # at least 10 segments 
    9393    pars = dict( 
    9494        length=length, 
    … …  
    100100tests = [ 
    101101    # Accuracy tests based on content in test/utest_other_models.py 
    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], 
     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], 
    109110 
    110111    # Additional tests with larger range of parameters 
    111     [{'length':    1000.0,  # test T2 
     112    [{'length':    1000.0, 
    112113      'kuhn_length': 100.0, 
    113114      'radius':       20.0, 
    … …  
    116117      'background':    0.0001, 
    117118     }, 1.0, 0.000595345], 
    118     [{'length':        10.0,  # test T3 
     119    [{'length':        10.0, 
    119120      'kuhn_length': 800.0, 
    120121      'radius':        2.0, 
    … …  
    123124      'background':    0.001, 
    124125     }, 0.1, 1.55228], 
    125     [{'length':        100.0,  # test T4 
     126    [{'length':        100.0, 
    126127      'kuhn_length': 800.0, 
    127128      'radius':       50.0, 
    … …  
    132133    ] 
    133134 
    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. 

  • sasmodels/models/lib/wrc_cyl.c

    r18a2bfc r92ce163  
    22    Functions for WRC implementation of flexible cylinders 
    33*/ 
    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.); 
     4double Sk_WR(double q, double L, double b); 
     5 
     6 
     7static double 
     8AlphaSquare(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// 
     20static double 
     21Rgsquarezero(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// 
     29static double 
     30Rgsquareshort(double q, double L, double b) 
     31{ 
     32    return AlphaSquare(L/b) * Rgsquarezero(q,L,b); 
     33} 
     34 
     35// 
     36static double 
     37Rgsquare(double q, double L, double b) 
     38{ 
     39    return AlphaSquare(L/b)*L*b/6.0; 
     40} 
     41 
     42static double 
     43sech_WR(double x) 
     44{ 
     45    return(1/cosh(x)); 
     46} 
     47 
     48static double 
     49a1long(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); 
    5055    } else { 
    51         return 2.*(expm1(-qsq) + qsq)/(qsq*qsq); 
     56        C = 1.0; 
    5257    } 
    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) 
     58 
    6359    const double C1 = 1.22; 
    6460    const double C2 = 0.4288; 
    … …  
    6864    const double miu = 0.585; 
    6965 
    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; 
     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; 
    10972    const double b3 = b*b*b; 
    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  
     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 
     150static double 
     151a2long(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// 
     233static double 
     234a_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} 
     254static double 
     255a1short(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 
     262static double 
     263a2short(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) 
     270static double 
     271w_WR(double x) 
     272{ 
     273    return 0.5*(1 + tanh((x - 1.523)/0.1477)); 
     274} 
     275 
     276// 
     277static double 
     278u1(double q, double L, double b) 
     279{ 
     280    return Rgsquareshort(q,L,b)*q*q; 
     281} 
     282 
     283static double 
     284u_WR(double q, double L, double b) 
     285{ 
     286    return Rgsquare(q,L,b)*q*q; 
     287} 
     288 
     289static double 
     290Sdebye_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} 
     305static double 
     306Sdebye(double q, double L, double b) 
     307{ 
     308    double arg = u_WR(q,L,b); 
     309    return Sdebye_kernel(arg); 
     310} 
     311 
     312// 
     313static double 
     314Sdebye1(double q, double L, double b) 
     315{ 
     316    double arg = u1(q,L,b); 
     317    return Sdebye_kernel(arg); 
     318 
     319} 
     320 
     321// 
    142322static double 
    143323Sexv(double q, double L, double b) 
    144324{ 
    145     // Pedersen eq 13, corrected by Chen eq A.5, swapping w and 1-w 
     325 
    146326    const double C1=1.22; 
    147327    const double C2=0.4288; 
    148328    const double C3=-1.651; 
    149329    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 
     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 
     339static double 
     340Sexvnew(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; 
    172348    const double del=1.05; 
    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; 
     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 
     365double 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       } 
    183401    } 
    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 } 
     402 
     403  return(ans); 
     404} 
Note: See TracChangeset for help on using the changeset viewer.