Changeset 50e1e40 in sasmodels


Ignore:
Timestamp:
Mar 1, 2016 7:49:00 PM (9 years ago)
Author:
Paul Kienzle <pkienzle@…>
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:
ad90df9
Parents:
a4a7308
Message:

use lib functions to simplify barbell, capped_cylinder, cylinder, ellipsoid, triaxial_ellipsoid; fix barbell calculation

Location:
sasmodels/models
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/barbell.c

    r3c97ff0 r50e1e40  
    77 
    88//barbell kernel - same as dumbell 
    9 double _bell_kernel(double q, double h, double bell_radius, 
    10         double length, double sin_alpha, double cos_alpha); 
    11 double _bell_kernel(double q, double h, double bell_radius, 
    12         double length, double sin_alpha, double cos_alpha) 
     9static double 
     10_bell_kernel(double q, double h, double bell_radius, 
     11             double half_length, double sin_alpha, double cos_alpha) 
    1312{ 
     13    // translate a point in [-1,1] to a point in [lower,upper] 
    1414    const double upper = 1.0; 
    15     const double lower = -1.0*h/bell_radius; 
     15    const double lower = h/bell_radius; 
     16    const double zm = 0.5*(upper-lower); 
     17    const double zb = 0.5*(upper+lower); 
    1618 
     19    // cos term in integral is: 
     20    //    cos (q (R t - h + L/2) cos(alpha)) 
     21    // so turn it into: 
     22    //    cos (m t + b) 
     23    // where: 
     24    //    m = q R cos(alpha) 
     25    //    b = q(L/2-h) cos(alpha) 
     26    const double m = q*bell_radius*cos_alpha; // cos argument slope 
     27    const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 
     28    const double qrst = q*bell_radius*sin_alpha; // Q*R*sin(theta) 
    1729    double total = 0.0; 
    1830    for (int i = 0; i < 76; i++){ 
    19         const double t = 0.5*(Gauss76Z[i]*(upper-lower)+upper+lower); 
    20         const double arg1 = q*cos_alpha*(bell_radius*t+h+length*0.5); 
    21         const double arg2 = q*bell_radius*sin_alpha*sqrt(1.0-t*t); 
    22         const double be = (arg2 == 0.0 ? 0.5 :J1(arg2)/arg2); 
    23         const double Fq = cos(arg1)*(1.0-t*t)*be; 
     31        const double t = Gauss76Z[i]*zm + zb; 
     32        const double radical = 1.0 - t*t; 
     33        const double bj = J1c(qrst*sqrt(radical)); 
     34        const double Fq = cos(m*t + b) * radical * bj; 
    2435        total += Gauss76Wt[i] * Fq; 
    2536    } 
    26     const double integral = 0.5*(upper-lower)*total; 
    27     return 4.0*M_PI*bell_radius*bell_radius*bell_radius*integral; 
     37    // translate dx in [-1,1] to dx in [lower,upper] 
     38    const double integral = total*zm; 
     39    const double bell_Fq = 2*M_PI*cube(bell_radius)*integral; 
     40    return bell_Fq; 
    2841} 
    2942 
    … …  
    3447 
    3548    // bell radius should never be less than radius when this is called 
    36     const double hdist = sqrt(bell_radius*bell_radius - radius*radius); 
    37     const double p1 = 2.0*bell_radius*bell_radius*bell_radius/3.0; 
    38     const double p2 = bell_radius*bell_radius*hdist; 
    39     const double p3 = hdist*hdist*hdist/3.0; 
     49    const double hdist = sqrt(square(bell_radius) - square(radius)); 
     50    const double p1 = 2.0/3.0*cube(bell_radius); 
     51    const double p2 = square(bell_radius)*hdist; 
     52    const double p3 = cube(hdist)/3.0; 
    4053 
    41     return M_PI*radius*radius*length + 2.0*M_PI*(p1+p2-p3); 
     54    return M_PI*square(radius)*length + 2.0*M_PI*(p1+p2-p3); 
    4255} 
    4356 
    44 double Iq(double q, double sld, 
    45         double solvent_sld, 
    46         double bell_radius, 
    47         double radius, 
    48         double length) 
     57double Iq(double q, double sld, double solvent_sld, 
     58          double bell_radius, double radius, double length) 
    4959{ 
    50     double sn, cn; // slots to hold sincos function output 
     60    // Exclude invalid inputs. 
     61    if (bell_radius < radius) return NAN; 
     62    const double h = -sqrt(bell_radius*bell_radius - radius*radius); 
     63    const double half_length = 0.5*length; 
    5164 
    52     if (bell_radius < radius) return NAN; 
    53  
    54     const double lower = 0.0; 
    55     const double upper = M_PI_2; 
    56     const double h = sqrt(bell_radius*bell_radius-radius*radius); 
     65    // translate a point in [-1,1] to a point in [0, pi/2] 
     66    const double zm = M_PI_4; 
     67    const double zb = M_PI_4; 
    5768    double total = 0.0; 
    5869    for (int i = 0; i < 76; i++){ 
    59         const double alpha= 0.5*(Gauss76Z[i]*(upper-lower) + upper + lower); 
    60         SINCOS(alpha, sn, cn); 
     70        const double alpha= Gauss76Z[i]*zm + zb; 
     71        double sin_alpha, cos_alpha; // slots to hold sincos function output 
     72        SINCOS(alpha, sin_alpha, cos_alpha); 
    6173 
    62         const double bell_Fq = _bell_kernel(q, h, bell_radius, length, sn, cn); 
     74        const double bell_Fq = _bell_kernel(q, h, bell_radius, half_length, sin_alpha, cos_alpha); 
     75        const double bj = J1c(q*radius*sin_alpha); 
     76        const double si = sinc(q*half_length*cos_alpha); 
     77        const double cyl_Fq = M_PI*radius*radius*length*bj*si; 
     78        const double Aq = bell_Fq + cyl_Fq; 
     79        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
     80    } 
     81    // translate dx in [-1,1] to dx in [lower,upper] 
     82    const double form = total*zm; 
    6383 
    64         const double arg1 = q*length*0.5*cn; 
    65         const double arg2 = q*radius*sn; 
    66         // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1 
    67         const double be = (arg2 == 0.0 ? 0.5 :J1(arg2)/arg2); 
    68         const double si = (arg1 == 0.0 ? 1.0 :sin(arg1)/arg1); 
    69         const double cyl_Fq = M_PI*radius*radius*length*si*2.0*be; 
    70  
    71         const double Aq = cyl_Fq + bell_Fq; 
    72         total += Gauss76Wt[i] * Aq * Aq * sn; 
    73     } 
    74  
    75     const double form = total*(upper-lower)*0.5; 
    76  
    77     //Contrast and volume normalization 
     84    //Contrast 
    7885    const double s = (sld - solvent_sld); 
    79     return form*1.0e-4*s*s; //form_volume(bell_radius,radius,length); 
     86    return 1.0e-4 * s * s * form; 
    8087} 
    8188 
    8289 
    83  
    8490double Iqxy(double qx, double qy, 
    85         double sld, 
    86         double solvent_sld, 
    87         double bell_radius, 
    88         double radius, 
    89         double length, 
    90         double theta, 
    91         double phi) 
     91        double sld, double solvent_sld, 
     92        double bell_radius, double radius, double length, 
     93        double theta, double phi) 
    9294{ 
    93      double sn, cn; // slots to hold sincos function output 
     95    // Compute angle alpha between q and the cylinder axis 
     96    double sn, cn; // slots to hold sincos function output 
     97    SINCOS(theta*M_PI_180, sn, cn); 
     98    const double q = sqrt(qx*qx+qy*qy); 
     99    const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
     100    const double alpha = acos(cos_val); // rod angle relative to q 
    94101 
    95102    // Exclude invalid inputs. 
    96103    if (bell_radius < radius) return NAN; 
     104    const double h = -sqrt(square(bell_radius) - square(radius)); 
     105    const double half_length = 0.5*length; 
    97106 
    98     // Compute angle alpha between q and the cylinder axis 
    99     SINCOS(theta*M_PI_180, sn, cn); 
    100     // # The following correction factor exists in sasview, but it can't be 
    101     // # right, so we are leaving it out for now. 
    102     const double q = sqrt(qx*qx+qy*qy); 
    103     const double cos_val = cn*cos(phi*M_PI_180)*qx + sn*qy; 
    104     const double alpha = acos(cos_val); // rod angle relative to q 
    105     SINCOS(alpha, sn, cn); 
    106  
    107     const double h = sqrt(bell_radius*bell_radius - radius*radius); // negative h 
    108     const double bell_Fq = _bell_kernel(q, h, bell_radius, length, sn, cn)/sn; 
    109  
    110     const double besarg = q*radius*sn; 
    111     const double siarg = q*0.5*length*cn; 
    112     // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1 
    113     const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg); 
    114     const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 
    115     const double cyl_Fq = M_PI*radius*radius*length*2.0*bj*si; 
    116  
    117     // Volume weighted average F(q) 
     107    double sin_alpha, cos_alpha; // slots to hold sincos function output 
     108    SINCOS(alpha, sin_alpha, cos_alpha); 
     109    const double bell_Fq = _bell_kernel(q, h, bell_radius, half_length, sin_alpha, cos_alpha); 
     110    const double bj = J1c(q*radius*sin_alpha); 
     111    const double si = sinc(q*half_length*cos_alpha); 
     112    const double cyl_Fq = M_PI*radius*radius*length*bj*si; 
    118113    const double Aq = cyl_Fq + bell_Fq; 
    119114 
    120     // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1 
     115    // Multiply by contrast^2 and convert to cm-1 
    121116    const double s = (sld - solvent_sld); 
    122     return 1.0e-4 * Aq * Aq * s * s; // form_volume(radius, cap_radius, length); 
     117    return 1.0e-4 * square(s * Aq); 
    123118} 

  • sasmodels/models/barbell.py

    rdcdf29d r50e1e40  
    121121# pylint: enable=bad-whitespace, line-too-long 
    122122 
    123 source = ["lib/J1.c", "lib/gauss76.c", "barbell.c"] 
     123source = ["lib/J1c.c", "lib/gauss76.c", "barbell.c"] 
    124124 
    125125# parameters for demo 

  • sasmodels/models/capped_cylinder.c

    r139c528 r50e1e40  
    1313//   length is the cylinder length, or the separation between the lens halves 
    1414//   alpha is the angle of the cylinder wrt q. 
    15 double _cap_kernel(double q, double h, double cap_radius, double length, 
    16                  double sin_alpha, double cos_alpha); 
    17 double _cap_kernel(double q, double h, double cap_radius, double length, 
    18                  double sin_alpha, double cos_alpha) 
     15static double 
     16_cap_kernel(double q, double h, double cap_radius, 
     17                      double half_length, double sin_alpha, double cos_alpha) 
    1918{ 
    20     // For speed, we are pre-calculating terms which are constant over the 
    21     // kernel. 
     19    // translate a point in [-1,1] to a point in [lower,upper] 
    2220    const double upper = 1.0; 
    2321    const double lower = h/cap_radius; // integral lower bound 
    2422    const double zm = 0.5*(upper-lower); 
    2523    const double zb = 0.5*(upper+lower); 
     24 
    2625    // cos term in integral is: 
    2726    //    cos (q (R t - h + L/2) cos(alpha)) 
    … …  
    3231    //    b = q(L/2-h) cos(alpha) 
    3332    const double m = q*cap_radius*cos_alpha; // cos argument slope 
    34     const double b = q*(0.5*length-h)*cos_alpha; // cos argument intercept 
     33    const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 
    3534    const double qrst = q*cap_radius*sin_alpha; // Q*R*sin(theta) 
    3635    double total = 0.0; 
    3736    for (int i=0; i<76 ;i++) { 
    38         // translate a point in [-1,1] to a point in [lower,upper] 
    39         //const double t = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    4037        const double t = Gauss76Z[i]*zm + zb; 
    4138        const double radical = 1.0 - t*t; 
    42         const double arg = qrst*sqrt(radical); // cap bessel function arg 
    43         const double be = (arg == 0.0 ? 0.5 : J1(arg)/arg); 
    44         const double Fq = cos(m*t + b) * radical * be; 
     39        const double bj = J1c(qrst*sqrt(radical)); 
     40        const double Fq = cos(m*t + b) * radical * bj; 
    4541        total += Gauss76Wt[i] * Fq; 
    4642    } 
    4743    // translate dx in [-1,1] to dx in [lower,upper] 
    48     //const double form = (upper-lower)/2.0*total; 
    49     const double integral = 0.5*(upper-lower)*total; 
    50     return 4.0*M_PI*cap_radius*cap_radius*cap_radius*integral; 
     44    const double integral = total*zm; 
     45    const double cap_Fq = 2*M_PI*cube(cap_radius)*integral; 
     46    return cap_Fq; 
    5147} 
    5248 
    … …  
    6662    // The first part is clearly V_cyl.  The second part requires some work: 
    6763    //    (R-h)^2 => h_c^2 
    68     //    (2R+h) => 2R+ h_c-h_c + h => 2R + (R-h)-hc + h => 3R-h_c 
     64    //    (2R+h) => 2R+ h_c-h_c + h => 2R + (R-h)-h_c + h => 3R-h_c 
    6965    // And so: 
    7066    //    2/3 pi (R-h)^2 (2R + h) => 2/3 pi h_c^2 (3 r_cap - h_c) 
    … …  
    7773    //        = pi (r^2 (L+hc) + hc^3/3) 
    7874    const double hc = cap_radius - sqrt(cap_radius*cap_radius - radius*radius); 
    79     return M_PI*(radius*radius*(length+hc) + 0.333333333333333*hc*hc*hc); 
     75    return M_PI*(radius*radius*(length+hc) + hc*hc*hc/3.0); 
    8076} 
    8177 
    82 double Iq(double q, 
    83     double sld, 
    84     double solvent_sld, 
    85     double radius, 
    86     double cap_radius, 
    87     double length) 
     78double Iq(double q, double sld, double solvent_sld, 
     79          double radius, double cap_radius, double length) 
    8880{ 
    89     double sn, cn; // slots to hold sincos function output 
    90  
    9181    // Exclude invalid inputs. 
    9282    if (cap_radius < radius) return NAN; 
     83    const double h = sqrt(cap_radius*cap_radius - radius*radius); 
     84    const double half_length = 0.5*length; 
    9385 
    94     const double lower = 0.0; 
    95     const double upper = M_PI_2; 
    96     const double h = sqrt(cap_radius*cap_radius - radius*radius); // negative h 
     86    // translate a point in [-1,1] to a point in [0, pi/2] 
     87    const double zm = M_PI_4; 
     88    const double zb = M_PI_4; 
    9789    double total = 0.0; 
    9890    for (int i=0; i<76 ;i++) { 
    99         // translate a point in [-1,1] to a point in [lower,upper] 
    100         const double alpha= 0.5*(Gauss76Z[i]*(upper-lower) + upper + lower); 
    101         SINCOS(alpha, sn, cn); 
     91        const double alpha= Gauss76Z[i]*zm + zb; 
     92        double sin_alpha, cos_alpha; // slots to hold sincos function output 
     93        SINCOS(alpha, sin_alpha, cos_alpha); 
    10294 
    103         const double cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn); 
    104  
    105         // The following is CylKernel() / sin(alpha), but we are doing it in place 
    106         // to avoid sin(alpha)/sin(alpha) for alpha = 0.  It is also a teensy bit 
    107         // faster since we don't multiply and divide sin(alpha). 
    108         const double besarg = q*radius*sn; 
    109         const double siarg = q*0.5*length*cn; 
    110         // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1 
    111         const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg); 
    112         const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 
    113         const double cyl_Fq = M_PI*radius*radius*length*2.0*bj*si; 
    114  
    115         // Volume weighted average F(q) 
    116         const double Aq = cyl_Fq + cap_Fq; 
    117         total += Gauss76Wt[i] * Aq * Aq * sn; // sn for spherical coord integration 
     95        const double cap_Fq = _cap_kernel(q, h, cap_radius, half_length, sin_alpha, cos_alpha); 
     96        const double bj = J1c(q*radius*sin_alpha); 
     97        const double si = sinc(q*half_length*cos_alpha); 
     98        const double cyl_Fq = M_PI*radius*radius*length*bj*si; 
     99        const double Aq = cap_Fq + cyl_Fq; 
     100        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; // sin_alpha for spherical coord integration 
    118101    } 
    119102    // translate dx in [-1,1] to dx in [lower,upper] 
    120     const double form = total * 0.5*(upper-lower); 
     103    const double form = total * zm; 
    121104 
    122     // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1 
    123     // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 
    124     // The additional volume factor is for polydisperse volume normalization. 
     105    // Contrast 
    125106    const double s = (sld - solvent_sld); 
    126     return 1.0e-4 * form * s * s; // form_volume(radius, cap_radius, length); 
     107    return 1.0e-4 * s * s * form; 
    127108} 
    128109 
    129110 
    130111double Iqxy(double qx, double qy, 
    131     double sld, 
    132     double solvent_sld, 
    133     double radius, 
    134     double cap_radius, 
    135     double length, 
    136     double theta, 
    137     double phi) 
     112    double sld, double solvent_sld, double radius, 
     113    double cap_radius, double length, 
     114    double theta, double phi) 
    138115{ 
    139     double sn, cn; // slots to hold sincos function output 
     116    // Compute angle alpha between q and the cylinder axis 
     117    double sn, cn; 
     118    SINCOS(theta*M_PI_180, sn, cn); 
     119    const double q = sqrt(qx*qx+qy*qy); 
     120    const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
     121    const double alpha = acos(cos_val); // rod angle relative to q 
    140122 
    141123    // Exclude invalid inputs. 
    142124    if (cap_radius < radius) return NAN; 
     125    const double h = sqrt(cap_radius*cap_radius - radius*radius); 
     126    const double half_length = 0.5*length; 
    143127 
    144     // Compute angle alpha between q and the cylinder axis 
    145     SINCOS(theta*M_PI_180, sn, cn); 
    146     // # The following correction factor exists in sasview, but it can't be 
    147     // # right, so we are leaving it out for now. 
    148     const double q = sqrt(qx*qx+qy*qy); 
    149     const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
    150     const double alpha = acos(cos_val); // rod angle relative to q 
    151     SINCOS(alpha, sn, cn); 
     128    double sin_alpha, cos_alpha; // slots to hold sincos function output 
     129    SINCOS(alpha, sin_alpha, cos_alpha); 
     130    const double cap_Fq = _cap_kernel(q, h, cap_radius, half_length, sin_alpha, cos_alpha); 
     131    const double bj = J1c(q*radius*sin_alpha); 
     132    const double si = sinc(q*half_length*cos_alpha); 
     133    const double cyl_Fq = M_PI*radius*radius*length*bj*si; 
     134    const double Aq = cap_Fq + cyl_Fq; 
    152135 
    153     const double h = sqrt(cap_radius*cap_radius - radius*radius); // negative h 
    154     const double cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn); 
    155  
    156     const double besarg = q*radius*sn; 
    157     const double siarg = q*0.5*length*cn; 
    158     // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1 
    159     const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg); 
    160     const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 
    161     const double cyl_Fq = M_PI*radius*radius*length*2.0*bj*si; 
    162  
    163     // Volume weighted average F(q) 
    164     const double Aq = cyl_Fq + cap_Fq; 
    165  
    166     // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1 
     136    // Multiply by contrast^2 and convert to cm-1 
    167137    const double s = (sld - solvent_sld); 
    168     return 1.0e-4 * Aq * Aq * s * s; // form_volume(radius, cap_radius, length); 
     138    return 1.0e-4 * square(s * Aq); 
    169139} 

  • sasmodels/models/capped_cylinder.py

    rdcdf29d r50e1e40  
    147147# pylint: enable=bad-whitespace, line-too-long 
    148148 
    149 source = ["lib/J1.c", "lib/gauss76.c", "capped_cylinder.c"] 
     149source = ["lib/J1c.c", "lib/gauss76.c", "capped_cylinder.c"] 
    150150 
    151151demo = dict(scale=1, background=0, 

  • sasmodels/models/cylinder.c

    r3c97ff0 r50e1e40  
    33double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    44    double radius, double length, double theta, double phi); 
    5  
    6 // twovd = 2 * volume * delta_rho 
    7 // besarg = q * R * sin(alpha) 
    8 // siarg = q * L/2 * cos(alpha) 
    9 double _cyl(double besarg, double siarg); 
    10 double _cyl(double besarg, double siarg) 
    11 { 
    12     const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg); 
    13     const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 
    14     return si*bj; 
    15 } 
    165 
    176double form_volume(double radius, double length) 
    … …  
    2615    double length) 
    2716{ 
     17    // TODO: return NaN if radius<0 or length<0? 
     18    // precompute qr and qh to save time in the loop 
    2819    const double qr = q*radius; 
    2920    const double qh = q*0.5*length; 
     21 
     22    // translate a point in [-1,1] to a point in [0, pi/2] 
     23    const double zm = M_PI_4; 
     24    const double zb = M_PI_4; 
     25 
    3026    double total = 0.0; 
    31     // double lower=0, upper=M_PI_2; 
    3227    for (int i=0; i<76 ;i++) { 
    33         // translate a point in [-1,1] to a point in [lower,upper] 
    34         //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    35         const double alpha = M_PI_4*(Gauss76Z[i] + 1.0); 
     28        const double alpha = Gauss76Z[i]*zm + zb; 
    3629        double sn, cn; 
    3730        SINCOS(alpha, sn, cn); 
    38         // For a bit of efficiency, we are moving the 2 V delta rho constant 
    39         // factor, 2Vdrho, out of the loop, so this is fq/2Vdrho rather than fq. 
    40         const double fq = _cyl(qr*sn, qh*cn); 
    41         total += Gauss76Wt[i] * fq * fq * sn; 
     31        const double fq = sinc(qh*cn) * J1c(qr*sn); 
     32        total += Gauss76Wt[i] * fq*fq * sn; 
    4233    } 
    4334    // translate dx in [-1,1] to dx in [lower,upper] 
    44     //const double form = (upper-lower)/2.0*total; 
    45     const double twoVdrho = 2.0*(sld-solvent_sld)*form_volume(radius, length); 
    46     return 1.0e-4 * twoVdrho * twoVdrho * total * M_PI_4; 
     35    const double form = total*zm; 
     36    const double s = (sld - solvent_sld) * form_volume(radius, length); 
     37    return 1.0e-4 * s * s * form; 
    4738} 
    4839 
    … …  
    5647    double phi) 
    5748{ 
    58     // TODO: check that radius<0 and length<0 give zero scattering. 
    59     // This should be the case since the polydispersity weight vector should 
    60     // be zero length, and this function never called. 
     49    // TODO: return NaN if radius<0 or length<0? 
    6150    double sn, cn; // slots to hold sincos function output 
    6251 
    6352    // Compute angle alpha between q and the cylinder axis 
    6453    SINCOS(theta*M_PI_180, sn, cn); 
    65     const double q = sqrt(qx*qx+qy*qy); 
     54    const double q = sqrt(qx*qx + qy*qy); 
    6655    const double cos_val = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q); 
    6756    const double alpha = acos(cos_val); 
     57 
    6858    SINCOS(alpha, sn, cn); 
    69     //sn = sqrt(1.0 - cos_val*cos_val); 
    70     //sn = 1.0 - 0.5*cos_val*cos_val;  // if cos_val is very small 
    71     //cn = cos_val; 
    72  
    73     const double twovd = 2.0*(sld-solvent_sld)*form_volume(radius, length); 
    74     const double fq = twovd * _cyl(q*radius*sn, q*0.5*length*cn); 
    75     return 1.0e-4 * fq * fq; 
     59    const double fq = sinc(q*0.5*length*cn) * J1c(q*radius*sn); 
     60    const double s = (sld-solvent_sld) * form_volume(radius, length); 
     61    return 1.0e-4 * square(s * fq); 
    7662} 

  • sasmodels/models/cylinder.py

    r0c3a226 r50e1e40  
    141141             ] 
    142142 
    143 source = ["lib/J1.c", "lib/gauss76.c", "cylinder.c"] 
     143source = ["lib/J1c.c", "lib/gauss76.c", "cylinder.c"] 
    144144 
    145145def ER(radius, length): 

  • sasmodels/models/ellipsoid.c

    r9c461c7 r50e1e40  
    1717double form_volume(double rpolar, double requatorial) 
    1818{ 
    19     return 1.333333333333333*M_PI*rpolar*requatorial*requatorial; 
     19    return M_4PI_3*rpolar*requatorial*requatorial; 
    2020} 
    2121 
    … …  
    2626    double requatorial) 
    2727{ 
    28     //const double lower = 0.0; 
    29     //const double upper = 1.0; 
     28    // translate a point in [-1,1] to a point in [0, 1] 
     29    const double zm = 0.5; 
     30    const double zb = 0.5; 
    3031    double total = 0.0; 
    3132    for (int i=0;i<76;i++) { 
    3233        //const double sin_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 
    33         const double sin_alpha = 0.5*(Gauss76Z[i] + 1.0); 
     34        const double sin_alpha = Gauss76Z[i]*zm + zb; 
    3435        total += Gauss76Wt[i] * _ellipsoid_kernel(q, rpolar, requatorial, sin_alpha); 
    3536    } 
    36     //const double form = (upper-lower)/2*total; 
    37     const double form = 0.5*total; 
     37    // translate dx in [-1,1] to dx in [lower,upper] 
     38    const double form = total*zm; 
    3839    const double s = (sld - solvent_sld) * form_volume(rpolar, requatorial); 
    39     return 1.0e-4 * form * s * s; 
     40    return 1.0e-4 * s * s * form; 
    4041} 
    4142 

  • sasmodels/models/sphere.py

    rc691551 r50e1e40  
    8383# This should perhaps be volume normalized? 
    8484form_volume = """ 
    85     return 1.333333333333333*M_PI*radius*radius*radius; 
     85    return M_4PI_3*cube(radius); 
    8686    """ 
    8787 

  • sasmodels/models/triaxial_ellipsoid.c

    r9c461c7 r50e1e40  
    2020 
    2121    double sn, cn; 
    22     //double st, ct; 
    23     //const double lower = 0.0; 
    24     //const double upper = 1.0; 
     22    // translate a point in [-1,1] to a point in [0, 1] 
     23    const double zm = 0.5; 
     24    const double zb = 0.5; 
    2525    double outer = 0.0; 
    2626    for (int i=0;i<76;i++) { 
    … …  
    3434        double inner = 0.0; 
    3535        for (int j=0;j<76;j++) { 
    36             const double y = 0.5*(Gauss76Z[j] + 1.0); 
    37             const double t = q*sqrt(acosx2 + bsinx2*(1.0-y*y) + c2*y*y); 
     36            const double ysq = square(Gauss76Z[j]*zm + zb); 
     37            const double t = q*sqrt(acosx2 + bsinx2*(1.0-ysq) + c2*ysq); 
    3838            const double fq = sph_j1c(t); 
    3939            inner += Gauss76Wt[j] * fq * fq ; 
    … …  
    4141        outer += Gauss76Wt[i] * 0.5 * inner; 
    4242    } 
    43     //const double fq2 = (upper-lower)/2*outer; 
    44     const double fq2 = 0.5*outer; 
     43    // translate dx in [-1,1] to dx in [lower,upper] 
     44    const double fqsq = outer*zm; 
    4545    const double s = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar); 
    46     return 1.0e-4 * fq2 * s * s; 
     46    return 1.0e-4 * s * s * fqsq; 
    4747} 
    4848 
    … …  
    6262    double sphi, cphi; 
    6363    double spsi, cpsi; 
    64     double st, ct; 
    6564 
    6665    const double q = sqrt(qx*qx + qy*qy); 
    … …  
    7675                          + req_major*req_major*cmu*cmu 
    7776                          + rpolar*rpolar*calpha*calpha); 
    78     SINCOS(t, st, ct); 
    79     const double fq = ( t==0.0 ? 1.0 : 3.0*(st-t*ct)/(t*t*t) ); 
     77    const double fq = sph_j1c(t); 
    8078    const double s = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar); 
    8179 
    82     return 1.0e-4 * fq * fq * s * s; 
     80    return 1.0e-4 * square(s * fq); 
    8381} 
    8482 
Note: See TracChangeset for help on using the changeset viewer.