Changeset 50e1e40 in sasmodels for sasmodels/models/capped_cylinder.c


Ignore:
Timestamp:
Mar 1, 2016 5:49:00 PM (8 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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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} 
Note: See TracChangeset for help on using the changeset viewer.