Changeset e7678b2 in sasmodels


Ignore:
Timestamp:
Feb 29, 2016 8:21:55 AM (7 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

Location:
sasmodels
Files:
1 added
11 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/kernel_template.c

    rdeac08c re7678b2  
    8888#  define M_PI_4 0.7853981633974483 
    8989#endif 
     90#ifndef M_E 
     91#  define M_E 2.718281828459045091 
     92#endif 
    9093 
    9194// Non-standard function library 
  • sasmodels/models/core_shell_bicelle.c

    r8007311 re7678b2  
    2929} 
    3030 
     31inline double sinc(double x) {return x==0 ? 1.0 : sin(x)/x;} 
     32 
    3133static double 
    3234bicelle_kernel(double qq, 
     
    4143              double dum) 
    4244{ 
    43         double dr1,dr2,dr3; 
    44         double besarg1,besarg2; 
    45         double vol1,vol2,vol3; 
    46         double sinarg1,sinarg2; 
    47         double t1,t2,t3; 
    48         double retval,si1,si2,be1,be2; 
     45    double si1,si2,be1,be2; 
    4946 
    50         dr1 = rhoc-rhoh; 
    51         dr2 = rhor-rhosolv; 
    52         dr3=  rhoh-rhor; 
    53         vol1 = M_PI*rad*rad*(2.0*length); 
    54         vol2 = M_PI*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick); 
    55         vol3= M_PI*(rad)*(rad)*(2.0*length+2.0*facthick); 
    56         besarg1 = qq*rad*sin(dum); 
    57         besarg2 = qq*(rad+radthick)*sin(dum); 
    58         sinarg1 = qq*length*cos(dum); 
    59         sinarg2 = qq*(length+facthick)*cos(dum); 
     47    const double dr1 = rhoc-rhoh; 
     48    const double dr2 = rhor-rhosolv; 
     49    const double dr3 = rhoh-rhor; 
     50    const double vol1 = M_PI*rad*rad*(2.0*length); 
     51    const double vol2 = M_PI*(rad+radthick)*(rad+radthick)*2.0*(length+facthick); 
     52    const double vol3 = M_PI*rad*rad*2.0*(length+facthick); 
     53    double sn,cn; 
     54    SINCOS(dum, sn, cn); 
     55    double besarg1 = qq*rad*sn; 
     56    double besarg2 = qq*(rad+radthick)*sn; 
     57    double sinarg1 = qq*length*cn; 
     58    double sinarg2 = qq*(length+facthick)*cn; 
    6059 
    61         if(besarg1 == 0) { 
    62                 be1 = 0.5; 
    63         } else { 
    64                 be1 = J1(besarg1)/besarg1; 
    65         } 
    66         if(besarg2 == 0) { 
    67                 be2 = 0.5; 
    68         } else { 
    69                 be2 = J1(besarg2)/besarg2; 
    70         } 
    71         if(sinarg1 == 0) { 
    72                 si1 = 1.0; 
    73         } else { 
    74                 si1 = sin(sinarg1)/sinarg1; 
    75         } 
    76         if(sinarg2 == 0) { 
    77                 si2 = 1.0; 
    78         } else { 
    79                 si2 = sin(sinarg2)/sinarg2; 
    80         } 
    81         t1 = 2.0*vol1*dr1*si1*be1; 
    82         t2 = 2.0*vol2*dr2*si2*be2; 
    83         t3 = 2.0*vol3*dr3*si2*be1; 
     60    be1 = J1c(besarg1); 
     61    be2 = J1c(besarg2); 
     62    si1 = sinc(sinarg1); 
     63    si2 = sinc(sinarg2); 
    8464 
    85         retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum); 
    86         return(retval); 
     65    const double t = vol1*dr1*si1*be1 + 
     66                     vol2*dr2*si2*be2 + 
     67                     vol3*dr3*si2*be1; 
     68 
     69    const double retval = t*t*sn; 
     70 
     71    return(retval); 
    8772 
    8873} 
     
    9984                   double rhosolv) 
    10085{ 
     86    // set up the integration end points 
     87    const double uplim = M_PI/4; 
     88    const double halfheight = length/2.0; 
    10189 
     90    double summ = 0.0; 
     91    for(int i=0;i<N_POINTS_76;i++) { 
     92        double zi = (Gauss76Z[i] + 1.0)*uplim; 
     93        double yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick, 
     94                             halfheight, rhoc, rhoh, rhor,rhosolv, zi); 
     95        summ += yyy; 
     96    } 
    10297 
    103         double answer,halfheight; 
    104         double lolim,uplim,summ,yyy,zi; 
    105         int nord,i; 
    106  
    107         // set up the integration end points 
    108         nord = 76; 
    109         lolim = 0.0; 
    110         uplim = M_PI/2; 
    111         halfheight = length/2.0; 
    112  
    113         summ = 0.0; 
    114         i=0; 
    115         for(i=0;i<nord;i++) { 
    116                 zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 
    117                 yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick, 
    118                                      halfheight, rhoc, rhoh, rhor,rhosolv, zi); 
    119                 summ += yyy; 
    120         } 
    121  
    122         // calculate value of integral to return 
    123         answer = (uplim-lolim)/2.0*summ; 
    124         return(answer); 
     98    // calculate value of integral to return 
     99    double answer = uplim*summ; 
     100    return(answer); 
    125101} 
    126102 
     
    138114          double phi) 
    139115{ 
    140     double cyl_x, cyl_y; 
    141     double alpha, cos_val; 
    142     double answer; 
    143  
    144116    //convert angle degree to radian 
    145     theta *= M_PI/180.0; 
    146     phi *= M_PI/180.0; 
     117    theta *= M_PI_180; 
     118    phi *= M_PI_180; 
    147119 
    148120    // Cylinder orientation 
    149     cyl_x = cos(theta) * cos(phi); 
    150     cyl_y = sin(theta); 
     121    const double cyl_x = cos(theta) * cos(phi); 
     122    const double cyl_y = sin(theta); 
    151123 
    152124    // Compute the angle btw vector q and the axis of the cylinder 
    153     cos_val = cyl_x*q_x + cyl_y*q_y; 
    154     alpha = acos( cos_val ); 
     125    const double cos_val = cyl_x*q_x + cyl_y*q_y; 
     126    const double alpha = acos( cos_val ); 
    155127 
    156128    // Get the kernel 
    157     answer = bicelle_kernel(q, radius, rim_thickness, face_thickness, 
     129    double answer = bicelle_kernel(q, radius, rim_thickness, face_thickness, 
    158130                           length/2.0, core_sld, face_sld, rim_sld, 
    159131                           solvent_sld, alpha) / fabs(sin(alpha)); 
  • sasmodels/models/core_shell_bicelle.py

    rfa8011eb re7678b2  
    9191# pylint: enable=bad-whitespace, line-too-long 
    9292 
    93 source = ["lib/J1.c", "lib/gauss76.c", "core_shell_bicelle.c"] 
     93source = ["lib/Si.c", "lib/J1.c", "lib/J1c.c", "lib/gauss76.c", "core_shell_bicelle.c"] 
    9494 
    9595demo = dict(scale=1, background=0, 
  • sasmodels/models/core_shell_ellipsoid.c

    r81dd619 re7678b2  
    4444          double solvent_sld) 
    4545{ 
    46         double delpc,delps; 
    47         double uplim,lolim;             //upper and lower integration limits 
    48         double summ,zi,yyy,answer; //running tally of integration 
    4946 
    50         lolim = 0.0; 
    51         uplim = 1.0; 
     47    //upper and lower integration limits 
     48    const double lolim = 0.0; 
     49    const double uplim = 1.0; 
    5250 
    53         summ = 0.0;      //initialize intergral 
     51    double summ = 0.0;   //initialize intergral 
    5452 
    55         delpc = core_sld - shell_sld; //core - shell 
    56         delps = shell_sld - solvent_sld; //shell - solvent 
     53    const double delpc = core_sld - shell_sld;    //core - shell 
     54    const double delps = shell_sld - solvent_sld; //shell - solvent 
    5755 
    58         for(int i=0;i<76;i++) { 
    59                 zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 
    60                 yyy = Gauss76Wt[i] * gfn4(zi, 
    61                                           equat_core, 
    62                                           polar_core, 
    63                                           equat_shell, 
    64                                           polar_shell, 
    65                                           delpc, 
    66                                           delps, 
    67                                           q); 
    68                 summ += yyy; 
    69         } 
     56    for(int i=0;i<N_POINTS_76;i++) { 
     57        double zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 
     58        double yyy = Gauss76Wt[i] * gfn4(zi, 
     59                                  equat_core, 
     60                                  polar_core, 
     61                                  equat_shell, 
     62                                  polar_shell, 
     63                                  delpc, 
     64                                  delps, 
     65                                  q); 
     66        summ += yyy; 
     67    } 
    7068 
    71         answer = (uplim-lolim)/2.0*summ; 
     69    double answer = (uplim-lolim)/2.0*summ; 
    7270 
    73         //convert to [cm-1] 
    74         answer *= 1.0e-4; 
     71    //convert to [cm-1] 
     72    answer *= 1.0e-4; 
    7573 
    76         return answer; 
     74    return answer; 
    7775} 
    7876 
     
    8987          double phi) 
    9088{ 
    91     double cyl_x, cyl_y; 
    92     double cos_val; 
    93     double answer; 
    94     double sldcs,sldss; 
    95  
    9689    //convert angle degree to radian 
    97     theta = theta * M_PI/180.0; 
    98     phi = phi * M_PI/180.0; 
     90    theta = theta * M_PI_180; 
     91    phi = phi * M_PI_180; 
    9992 
    10093 
    10194    // ellipsoid orientation, the axis of the rotation is consistent with the ploar axis. 
    102     cyl_x = cos(theta) * cos(phi); 
    103     cyl_y = sin(theta); 
     95    const double cyl_x = cos(theta) * cos(phi); 
     96    const double cyl_y = sin(theta); 
    10497 
    105     sldcs = core_sld - shell_sld; 
    106     sldss = shell_sld- solvent_sld; 
     98    const double sldcs = core_sld - shell_sld; 
     99    const double sldss = shell_sld- solvent_sld; 
    107100 
    108101    // Compute the angle btw vector q and the 
    109102    // axis of the cylinder 
    110     cos_val = cyl_x*q_x + cyl_y*q_y; 
     103    const double cos_val = cyl_x*q_x + cyl_y*q_y; 
    111104 
    112105    // Call the IGOR library function to get the kernel: MUST use gfn4 not gf2 because of the def of params. 
    113     answer = gfn4(cos_val, 
     106    double answer = gfn4(cos_val, 
    114107                  equat_core, 
    115108                  polar_core, 
  • sasmodels/models/core_shell_ellipsoid_xt.c

    r81bb668 re7678b2  
    3030                   double x_polar_shell) 
    3131{ 
    32         double equat_shell, polar_shell; 
    33     equat_shell = equat_core + t_shell; 
    34     polar_shell = equat_core*x_core + t_shell*x_polar_shell; 
     32    const double equat_shell = equat_core + t_shell; 
     33    const double polar_shell = equat_core*x_core + t_shell*x_polar_shell; 
    3534    double vol = 4.0*M_PI/3.0*equat_shell*equat_shell*polar_shell; 
    3635    return vol; 
     
    4746          double solvent_sld) 
    4847{ 
    49         double delpc,delps; 
    50         double uplim,lolim;             //upper and lower integration limits 
    51         double summ,zi,yyy,answer; //running tally of integration 
    52         double polar_core, equat_shell, polar_shell; 
     48    const double lolim = 0.0; 
     49    const double uplim = 1.0; 
    5350 
    54         lolim = 0.0; 
    55         uplim = 1.0; 
     51    double summ = 0.0;   //initialize intergral 
    5652 
    57         summ = 0.0;      //initialize intergral 
    58  
    59         delpc = core_sld - shell_sld; //core - shell 
    60         delps = shell_sld - solvent_sld; //shell - solvent 
     53    const double delpc = core_sld - shell_sld; //core - shell 
     54    const double delps = shell_sld - solvent_sld; //shell - solvent 
    6155 
    6256 
    63     polar_core = equat_core*x_core; 
    64     equat_shell = equat_core + t_shell; 
    65     polar_shell = equat_core*x_core + t_shell*x_polar_shell; 
     57    const double polar_core = equat_core*x_core; 
     58    const double equat_shell = equat_core + t_shell; 
     59    const double polar_shell = equat_core*x_core + t_shell*x_polar_shell; 
    6660 
    67         for(int i=0;i<76;i++) { 
    68                 zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 
    69                 yyy = Gauss76Wt[i] * gfn4(zi, 
    70                                           equat_core, 
     61    for(int i=0;i<N_POINTS_76;i++) { 
     62        double zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 
     63        double yyy = Gauss76Wt[i] * gfn4(zi, 
     64                                  equat_core, 
    7165                                  polar_core, 
    7266                                  equat_shell, 
    7367                                  polar_shell, 
    74                                           delpc, 
    75                                           delps, 
    76                                           q); 
    77                 summ += yyy; 
    78         } 
     68                                  delpc, 
     69                                  delps, 
     70                                  q); 
     71        summ += yyy; 
     72    } 
    7973 
    80         answer = (uplim-lolim)/2.0*summ; 
     74    double answer = (uplim-lolim)/2.0*summ; 
     75    //convert to [cm-1] 
     76    answer *= 1.0e-4; 
    8177 
    82         //convert to [cm-1] 
    83         answer *= 1.0e-4; 
    84  
    85         return answer; 
     78    return answer; 
    8679} 
    8780 
     
    9891          double phi) 
    9992{ 
    100     double cyl_x, cyl_y; 
    101     double cos_val; 
    102     double answer; 
    103     double sldcs,sldss; 
    104         double polar_core, equat_shell, polar_shell; 
    105  
    10693    //convert angle degree to radian 
    107     theta = theta * M_PI/180.0; 
    108     phi = phi * M_PI/180.0; 
    109  
     94    theta = theta * M_PI_180; 
     95    phi = phi * M_PI_180; 
    11096 
    11197    // ellipsoid orientation, the axis of the rotation is consistent with the ploar axis. 
    112     cyl_x = cos(theta) * cos(phi); 
    113     cyl_y = sin(theta); 
     98    const double cyl_x = cos(theta) * cos(phi); 
     99    const double cyl_y = sin(theta); 
    114100 
    115     sldcs = core_sld - shell_sld; 
    116     sldss = shell_sld- solvent_sld; 
     101    const double sldcs = core_sld - shell_sld; 
     102    const double sldss = shell_sld- solvent_sld; 
    117103 
    118104    // Compute the angle btw vector q and the 
    119105    // axis of the cylinder 
    120     cos_val = cyl_x*q_x + cyl_y*q_y; 
     106    const double cos_val = cyl_x*q_x + cyl_y*q_y; 
    121107 
    122     polar_core = equat_core*x_core; 
    123     equat_shell = equat_core + t_shell; 
    124     polar_shell = equat_core*x_core + t_shell*x_polar_shell; 
     108    const double polar_core = equat_core*x_core; 
     109    const double equat_shell = equat_core + t_shell; 
     110    const double polar_shell = equat_core*x_core + t_shell*x_polar_shell; 
    125111 
    126112    // Call the IGOR library function to get the kernel: 
    127113    // MUST use gfn4 not gf2 because of the def of params. 
    128     answer = gfn4(cos_val, 
     114    double answer = gfn4(cos_val, 
    129115                  equat_core, 
    130116                  polar_core, 
  • sasmodels/models/flexible_cylinder.c

    rf94d8a2 re7678b2  
    2222{ 
    2323 
    24         double Pi = 4.0*atan(1.0); 
    25  
    26         double cont = sld-solvent_sld; 
    27         double qr = q*radius; 
    28         double flex = Sk_WR(q,length,kuhn_length); 
    29     double crossSect = (2.0*J1(qr)/qr)*(2.0*J1(qr)/qr); 
    30  
    31         flex *= crossSect; 
    32         flex *= Pi*radius*radius*length; 
    33         flex *= cont*cont; 
    34         flex *= 1.0e-4; 
    35  
    36         return flex; 
     24    const double cont = sld-solvent_sld; 
     25    const double qr = q*radius; 
     26    //const double crossSect = (2.0*J1(qr)/qr)*(2.0*J1(qr)/qr); 
     27    const double crossSect = J1c(qr); 
     28    double flex = Sk_WR(q,length,kuhn_length); 
     29    flex *= crossSect*crossSect; 
     30    flex *= M_PI*radius*radius*length; 
     31    flex *= cont*cont; 
     32    flex *= 1.0e-4; 
     33    return flex; 
    3734} 
    3835 
     
    4542{ 
    4643 
    47         double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld); 
    48  
    49         return result; 
     44    double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld); 
     45    return result; 
    5046} 
    5147 
     
    5753            double solvent_sld) 
    5854{ 
    59         double q; 
    60         q = sqrt(qx*qx+qy*qy); 
    61         double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld); 
     55    double q; 
     56    q = sqrt(qx*qx+qy*qy); 
     57    double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld); 
    6258 
    63         return result; 
     59    return result; 
    6460} 
  • sasmodels/models/flexible_cylinder.py

    r13ed84c re7678b2  
    9090    ] 
    9191# pylint: enable=bad-whitespace, line-too-long 
    92 source = ["lib/J1.c", "lib/wrc_cyl.c", "flexible_cylinder.c"] 
     92source = ["lib/J1.c", "lib/J1c.c", "lib/wrc_cyl.c", "flexible_cylinder.c"] 
    9393 
    9494demo = dict(scale=1.0, background=0.0001, 
  • sasmodels/models/flexible_cylinder_ex.c

    r504abee re7678b2  
    1717elliptical_crosssection(double q, double a, double b) 
    1818{ 
    19     double uplim,lolim,Pi,summ,arg,zi,yyy,answer; 
    20     int i,nord=76; 
     19    double summ=0.0; 
    2120 
    22     Pi = 4.0*atan(1.0); 
    23     lolim=0.0; 
    24     uplim=Pi/2.0; 
    25     summ=0.0; 
     21    for(int i=0;i<N_POINTS_76;i++) { 
     22        double zi = ( Gauss76Z[i] + 1.0 )*M_PI/4.0; 
     23        double sn, cn; 
     24        SINCOS(zi, sn, cn); 
     25        double arg = q*sqrt(a*a*sn*sn+b*b*cn*cn); 
     26        double yyy = pow((double)J1c(arg),2); 
     27        yyy *= Gauss76Wt[i]; 
     28        summ += yyy; 
     29    } 
    2630 
    27     for(i=0;i<nord;i++) { 
    28                 zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 
    29                 arg = q*sqrt(a*a*sin(zi)*sin(zi)+b*b*cos(zi)*cos(zi)); 
    30                 yyy = pow((2.0 * J1(arg) / arg),2); 
    31                 yyy *= Gauss76Wt[i]; 
    32                 summ += yyy; 
    33     } 
    34     answer = (uplim-lolim)/2.0*summ; 
    35     answer *= 2.0/Pi; 
    36     return(answer); 
     31    summ /= 2.0; 
     32    return(summ); 
    3733 
    3834} 
     
    4743{ 
    4844 
    49         double Pi,flex,crossSect, cont; 
     45    double flex,crossSect, cont; 
    5046 
    51         Pi = 4.0*atan(1.0); 
    52         cont = sld - solvent_sld; 
    53         crossSect = elliptical_crosssection(q,radius,(radius*axis_ratio)); 
     47    cont = sld - solvent_sld; 
     48    crossSect = elliptical_crosssection(q,radius,(radius*axis_ratio)); 
    5449 
    55         flex = Sk_WR(q,length,kuhn_length); 
    56         flex *= crossSect; 
    57         flex *= Pi*radius*radius*axis_ratio*axis_ratio*length; 
    58         flex *= cont*cont; 
    59         flex *= 1.0e-4; 
     50    flex = Sk_WR(q,length,kuhn_length); 
     51    flex *= crossSect; 
     52    flex *= M_PI*radius*radius*axis_ratio*axis_ratio*length; 
     53    flex *= cont*cont; 
     54    flex *= 1.0e-4; 
    6055 
    61         return flex; 
     56    return flex; 
    6257} 
    6358 
     
    7166{ 
    7267 
    73         double result = flexible_cylinder_ex_kernel(q, 
    74                         length, 
    75                         kuhn_length, 
    76                         radius, 
    77                         axis_ratio, 
    78                         sld, 
    79                         solvent_sld); 
     68    double result = flexible_cylinder_ex_kernel(q, 
     69                    length, 
     70                    kuhn_length, 
     71                    radius, 
     72                    axis_ratio, 
     73                    sld, 
     74                    solvent_sld); 
    8075 
    81         return result; 
     76    return result; 
    8277} 
    8378 
     
    9085            double solvent_sld) 
    9186{ 
    92         double q; 
    93         q = sqrt(qx*qx+qy*qy); 
    94         double result = flexible_cylinder_ex_kernel(q, 
    95                         length, 
    96                         kuhn_length, 
    97                         radius, 
    98                         axis_ratio, 
    99                         sld, 
    100                         solvent_sld); 
     87    double q; 
     88    q = sqrt(qx*qx+qy*qy); 
     89    double result = flexible_cylinder_ex_kernel(q, 
     90                    length, 
     91                    kuhn_length, 
     92                    radius, 
     93                    axis_ratio, 
     94                    sld, 
     95                    solvent_sld); 
    10196 
    102         return result; 
     97    return result; 
    10398} 
  • sasmodels/models/flexible_cylinder_ex.py

    r13ed84c re7678b2  
    113113# pylint: enable=bad-whitespace, line-too-long 
    114114 
    115 source = ["lib/J1.c", "lib/gauss76.c", "lib/wrc_cyl.c", "flexible_cylinder_ex.c"] 
     115source = ["lib/J1.c", "lib/J1c.c", "lib/gauss76.c", "lib/wrc_cyl.c", "flexible_cylinder_ex.c"] 
    116116 
    117117demo = dict(scale=1.0, background=0.0001, 
  • sasmodels/models/lib/Si.c

    rf12357f re7678b2  
    44double Si(double x) 
    55{ 
    6         double out; 
    7         double pi = 4.0*atan(1.0); 
     6    double out; 
    87 
    9         if (x >= pi*6.2/4.0){ 
    10                 double out_sin = 0.0; 
    11                 double out_cos = 0.0; 
    12                 out = pi/2.0; 
     8    if (x >= M_PI*6.2/4.0){ 
     9        double out_sin = 0.0; 
     10        double out_cos = 0.0; 
     11        out = M_PI/2.0; 
    1312 
    14                 // Explicitly writing factorial values triples the speed of the calculation 
    15                 out_cos = 1./x - 2./pow(x,3) + 24./pow(x,5) - 720./pow(x,7); 
    16                 out_sin = 1./pow(x,2) - 6./pow(x,4) + 120./pow(x,6) - 5040./pow(x,8); 
     13        // Explicitly writing factorial values triples the speed of the calculation 
     14        out_cos = 1./x - 2./pow(x,3) + 24./pow(x,5) - 720./pow(x,7); 
     15        out_sin = 1./pow(x,2) - 6./pow(x,4) + 120./pow(x,6) - 5040./pow(x,8); 
    1716 
    18                 out -= cos(x) * out_cos; 
    19                 out -= sin(x) * out_sin; 
    20                 return out; 
    21         } 
     17        out -= cos(x) * out_cos; 
     18        out -= sin(x) * out_sin; 
     19        return out; 
     20    } 
    2221 
    23         // Explicitly writing factorial values triples the speed of the calculation 
    24         out = x - pow(x, 3)/18. + pow(x,5)/600. - pow(x,7)/35280. + pow(x,9)/3265920. - pow(x,11)/439084800.; 
     22    // Explicitly writing factorial values triples the speed of the calculation 
     23    out = x - pow(x, 3)/18. + pow(x,5)/600. - pow(x,7)/35280. + pow(x,9)/3265920. - pow(x,11)/439084800.; 
    2524 
    26         return out; 
     25    return out; 
    2726} 
  • 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.