Changeset e7678b2 in sasmodels for sasmodels/models/lib/wrc_cyl.c


Ignore:
Timestamp:
Feb 29, 2016 6:21:55 AM (8 years ago)
Author:
piotr
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
73860b6
Parents:
deac08c
Message:

Code review from PAK

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/lib/wrc_cyl.c

    r13ed84c re7678b2  
    55AlphaSquare(double x) 
    66{ 
     7    // Potentially faster. Needs proper benchmarking. 
     8    // add native_powr to kernel_template 
     9    //double t = native_powr( (1.0 + (x/3.12)*(x/3.12) + 
     10    //     (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); 
     11    //return t; 
     12 
    713    return pow( (1.0 + (x/3.12)*(x/3.12) + 
    814         (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); 
    … …  
    1319Rgsquarezero(double q, double L, double b) 
    1420{ 
    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)))); 
     21    const double r = b/L; 
     22    return (L*b/6.0) * 
     23           (1.0 - r*1.5  + 1.5*r*r - 0.75*r*r*r*(1.0 - exp(-2.0/r))); 
    1724} 
    1825 
    … …  
    3138} 
    3239 
    33 static double 
     40static inline double 
    3441sech_WR(double x) 
    3542{ 
    … …  
    4047a1long(double q, double L, double b, double p1, double p2, double q0) 
    4148{ 
    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    double C; 
     50    const double onehalf = 1.0/2.0; 
    4951 
    5052    if( L/b > 10.0) { 
    … …  
    5456    } 
    5557 
    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))))))))))); 
     58    const double C1 = 1.22; 
     59    const double C2 = 0.4288; 
     60    const double C3 = -1.651; 
     61    const double C4 = 1.523; 
     62    const double C5 = 0.1477; 
     63    const double miu = 0.585; 
     64 
     65    const double Rg2 = Rgsquare(q,L,b); 
     66    const double Rg22 = Rg2*Rg2; 
     67    const double Rg = sqrt(Rg2); 
     68    const double Rgb = Rg*q0/b; 
     69 
     70    const double b2 = b*b; 
     71    const double b3 = b*b*b; 
     72    const double b4 = b3*b; 
     73    const double q02 = q0*q0; 
     74    const double q03 = q0*q0*q0; 
     75    const double q04 = q03*q0; 
     76    const double q05 = q04*q0; 
     77 
     78    const double Rg02 = Rg2*q02; 
     79 
     80    const double t1 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2))) * 
     81         ((11.0/15.0 + (7.0*b2)/(15.0*Rg02))) + 
     82         (7.0*b2)/(15.0*Rg02)))); 
     83 
     84    const double t2 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
     85         Rg02/b2))*((1.0 + onehalf*(((-1.0) - 
     86         tanh((-C4 + Rgb/C5))))))); 
     87 
     88    const double t3 = ((C3*pow(Rgb,((-3.0)/miu)) + 
     89         C2*pow(Rgb,((-2.0)/miu)) + 
     90         C1*pow(Rgb,((-1.0)/miu)))); 
     91 
     92    const double t4 = ((1.0 + tanh(((-C4) + Rgb)/C5))); 
     93 
     94    const double t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - 
     95         b*p2*pow(q0,((-1.0) - p1 - p2)))); 
     96 
     97    const double t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) + 
     98         (14.0*b3*pow((double)M_E,(-(Rg02/b2))))/(15.0*q03*Rg2) + 
     99         (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*((11.0/15.0 + 
     100         (7.0*b2)/(15.0*Rg02)))*Rg2)/b))); 
     101 
     102    const double t7 = (Rg*((C3*pow(((Rgb)),((-3.0)/miu)) + 
     103         C2*pow(((Rgb)),((-2.0)/miu)) + 
     104         C1*pow(((Rgb)),((-1.0)/miu))))*pow(sech_WR(((-C4) + 
     105         Rgb)/C5),2)); 
     106 
     107    const double t8 = (b4*Rg*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
     108         Rg02/b2))*pow(sech_WR(((-C4) + Rgb)/C5),2)); 
     109 
     110    const double t9 = (2.0*b4*(((2.0*q0*Rg2)/b - 
     111         (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))*((1.0 + onehalf*(((-1.0) - 
     112         tanh(((-C4) + Rgb)/C5)))))); 
     113 
     114    const double t10 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
     115         Rg02/b2))*((1.0 + onehalf*(((-1.0) - tanh(((-C4) + 
     116         Rgb)/C5)))))); 
     117 
     118    const double t11 = (((-((3.0*C3*Rg*pow(((Rgb)),((-1.0) - 
     119          3.0/miu)))/miu)) - (2.0*C2*Rg*pow(((Rgb)),((-1.0) - 
     120          2.0/miu)))/miu - (C1*Rg*pow(((Rgb)),((-1.0) - 
     121          1.0/miu)))/miu)); 
     122 
     123    const double t12 = ((1.0 + tanh(((-C4) + Rgb)/C5))); 
     124 
     125    const double t13 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2)))*((11.0/15.0 + 
     126          (7.0*b2)/(15.0*q02* Rg2))) + 
     127          (7.0*b2)/(15.0*Rg02)))); 
     128 
     129    const double t14 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
     130          Rg02/b2))*((1.0 + onehalf*(((-1.0) - tanh(((-C4) + 
     131          Rgb)/C5)))))); 
     132 
     133    const double t15 = ((C3*pow(((Rgb)),((-3.0)/miu)) + 
     134        C2*pow(((Rgb)),((-2.0)/miu)) + 
     135        C1*pow(((Rgb)),((-1.0)/miu)))); 
     136 
     137 
     138    double yy = (pow(q0,p1)*(((-((b*M_PI)/(L*q0))) +t1/L +t2/(q04*Rg22) + 
     139        onehalf*t3*t4)) + (t5*((pow(q0,(p1 - p2))* 
     140        (((-pow(q0,(-p1)))*(((b2*M_PI)/(L*q02) +t6/L +t7/(2.0*C5) - 
     141        t8/(C5*q04*Rg22) + t9/(q04*Rg22) -t10/(q05*Rg22) + onehalf*t11*t12)) - 
     142        b*p1*pow(q0,((-1.0) - p1))*(((-((b*M_PI)/(L*q0))) + t13/L + 
     143        t14/(q04*Rg22) + onehalf*t15*((1.0 + tanh(((-C4) + 
     144        Rgb)/C5))))))))))); 
    137145 
    138146    return (yy); 
    … …  
    142150a2long(double q, double L, double b, double p1, double p2, double q0) 
    143151{ 
    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; 
     152    double C; 
     153    const double onehalf = 1.0/2.0; 
     154 
    150155    if( L/b > 10.0) { 
    151156        C = 3.06/pow((L/b),0.44); 
    … …  
    154159    } 
    155160 
    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)))))))))); 
     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 + onehalf*(((-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 + onehalf*(((-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 + onehalf*(((-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 + onehalf*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))* 
     223         (((-((b*M_PI)/(L*q0))) + t9 + t10 + 
     224         onehalf*((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)))))))))); 
    220228 
    221229    return (yy); 
    … …  
    224232// 
    225233static double 
    226 a1short(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)* 
     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)* 
    244247        ((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short + 
    245248        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)))))); 
     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)))))); 
    248251 
    249252    return(yy); 
    250253} 
     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} 
    251261 
    252262static double 
    253263a2short(double q, double L, double b, double p1short, double p2short, double q0) 
    254264{ 
    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); 
     265    double factor = -1.0; 
     266    return a_short(q, L, b, p2short, p1short, factor, p1short-p2short, q0); 
    277267} 
    278268 
    … …  
    301291{ 
    302292    // ORIGINAL 
    303     double result = 2.0*(exp(-arg) + arg -1.0)/(pow((arg),2)); 
     293    double result = 2.0*(exp(-arg) + arg -1.0)/(arg*arg); 
    304294 
    305295    // CONVERSION 1 from http://herbie.uwplse.org/ 
    … …  
    333323Sexv(double q, double L, double b) 
    334324{ 
    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  
     325 
     326    const double C1=1.22; 
     327    const double C2=0.4288; 
     328    const double C3=-1.651; 
     329    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)); 
    348335    return (yy); 
    349336} 
    … …  
    353340Sexvnew(double q, double L, double b) 
    354341{ 
    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; 
     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; 
     348    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 
    362352 
    363353    //calculating the derivative to decide on the corection (cutoff) term? 
    364354    // 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))); 
     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)); 
    377361 
    378362    return (yy); 
    … …  
    382366double Sk_WR(double q, double L, double b) 
    383367{ 
    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)) - 
     368    // 
     369    const double p1 = 4.12; 
     370    const double p2 = 4.42; 
     371    const double p1short = 5.36; 
     372    const double p2short = 5.62; 
     373    const double q0 = 3.1; 
     374 
     375    double q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0); 
     376    double Sexvmodify, ans; 
     377 
     378    const double C =(L/b > 10.0) ? 3.06/pow((L/b),0.44) : 1.0; 
     379 
     380    if( L > 4*b ) { // Longer Chains 
     381       if (q*b <= 3.1) {   //Modified by Yun on Oct. 15, 
     382         Sexvmodify = Sexvnew(q, L, b); 
     383         ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) - 
    413384            (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L); 
    414385 
    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); 
     386       } else { //q(i)*b > 3.1 
     387         ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) + 
     388               a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + M_PI/(q*L); 
     389       } 
     390    } else { //L <= 4*b Shorter Chains 
     391       if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) { 
     392         if (q*b<=0.01) { 
     393            ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0; 
     394         } else { 
     395            ans = Sdebye1(q,L,b); 
     396         } 
     397       } else {  //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3) 
     398         ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) + 
     399               a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) + 
     400               M_PI/(q*L); 
     401       } 
    418402    } 
    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   } 
    432403 
    433404  return(ans); 
Note: See TracChangeset for help on using the changeset viewer.