Changeset 994d77f in sasmodels for sasmodels/models/capped_cylinder.c


Ignore:
Timestamp:
Oct 30, 2014 10:33:53 AM (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:
ef2861b
Parents:
d087487b
Message:

Convert double to float rather than using real

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/capped_cylinder.c

    rf4cf580 r994d77f  
    1 real form_volume(real radius, real cap_radius, real length); 
    2 real Iq(real q, real sld, real solvent_sld, 
    3     real radius, real cap_radius, real length); 
    4 real Iqxy(real qx, real qy, real sld, real solvent_sld, 
    5     real radius, real cap_radius, real length, real theta, real phi); 
     1double form_volume(double radius, double cap_radius, double length); 
     2double Iq(double q, double sld, double solvent_sld, 
     3    double radius, double cap_radius, double length); 
     4double Iqxy(double qx, double qy, double sld, double solvent_sld, 
     5    double radius, double cap_radius, double length, double theta, double phi); 
    66 
    77// Integral over a convex lens kernel for t in [h/R,1].  See the docs for 
    … …  
    1313//   length is the cylinder length, or the separation between the lens halves 
    1414//   alpha is the angle of the cylinder wrt q. 
    15 real _cap_kernel(real q, real h, real cap_radius, real length, 
    16                  real sin_alpha, real cos_alpha); 
    17 real _cap_kernel(real q, real h, real cap_radius, real length, 
    18                  real sin_alpha, real cos_alpha) 
     15double _cap_kernel(double q, double h, double cap_radius, double length, 
     16                 double sin_alpha, double cos_alpha); 
     17double _cap_kernel(double q, double h, double cap_radius, double length, 
     18                 double sin_alpha, double cos_alpha) 
    1919{ 
    2020    // For speed, we are pre-calculating terms which are constant over the 
    2121    // kernel. 
    22     const real upper = REAL(1.0); 
    23     const real lower = h/cap_radius; // integral lower bound 
     22    const double upper = 1.0; 
     23    const double lower = h/cap_radius; // integral lower bound 
    2424    // cos term in integral is: 
    2525    //    cos (q (R t - h + L/2) cos(alpha)) 
    … …  
    2929    //    m = q R cos(alpha) 
    3030    //    b = q(L/2-h) cos(alpha) 
    31     const real m = q*cap_radius*cos_alpha; // cos argument slope 
    32     const real b = q*(REAL(0.5)*length-h)*cos_alpha; // cos argument intercept 
    33     const real qrst = q*cap_radius*sin_alpha; // Q*R*sin(theta) 
    34     real total = REAL(0.0); 
     31    const double m = q*cap_radius*cos_alpha; // cos argument slope 
     32    const double b = q*(0.5*length-h)*cos_alpha; // cos argument intercept 
     33    const double qrst = q*cap_radius*sin_alpha; // Q*R*sin(theta) 
     34    double total = 0.0; 
    3535    for (int i=0; i<76 ;i++) { 
    3636        // translate a point in [-1,1] to a point in [lower,upper] 
    37         //const real t = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    38         const real t = REAL(0.5)*(Gauss76Z[i]*(upper-lower)+upper+lower); 
    39         const real radical = REAL(1.0) - t*t; 
    40         const real arg = qrst*sqrt(radical); // cap bessel function arg 
    41         const real be = (arg == REAL(0.0) ? REAL(0.5) : J1(arg)/arg); 
    42         const real Fq = cos(m*t + b) * radical * be; 
     37        //const double t = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
     38        const double t = 0.5*(Gauss76Z[i]*(upper-lower)+upper+lower); 
     39        const double radical = 1.0 - t*t; 
     40        const double arg = qrst*sqrt(radical); // cap bessel function arg 
     41        const double be = (arg == 0.0 ? 0.5 : J1(arg)/arg); 
     42        const double Fq = cos(m*t + b) * radical * be; 
    4343        total += Gauss76Wt[i] * Fq; 
    4444    } 
    4545    // translate dx in [-1,1] to dx in [lower,upper] 
    46     //const real form = (upper-lower)/2.0*total; 
    47     const real integral = REAL(0.5)*(upper-lower)*total; 
    48     return REAL(4.0)*M_PI*cap_radius*cap_radius*cap_radius*integral; 
     46    //const double form = (upper-lower)/2.0*total; 
     47    const double integral = 0.5*(upper-lower)*total; 
     48    return 4.0*M_PI*cap_radius*cap_radius*cap_radius*integral; 
    4949} 
    5050 
    51 real form_volume(real radius, real cap_radius, real length) 
     51double form_volume(double radius, double cap_radius, double length) 
    5252{ 
    5353    // cap radius should never be less than radius when this is called 
    … …  
    6060    //        = pi r^2 L + pi hc (r^2 + hc^2/3) 
    6161    //        = pi * (r^2 (L+hc) + hc^3/3) 
    62     const real hc = cap_radius - sqrt(cap_radius*cap_radius - radius*radius); 
    63     return M_PI*(radius*radius*(length+hc) + REAL(0.333333333333333)*hc*hc*hc); 
     62    const double hc = cap_radius - sqrt(cap_radius*cap_radius - radius*radius); 
     63    return M_PI*(radius*radius*(length+hc) + 0.333333333333333*hc*hc*hc); 
    6464} 
    6565 
    66 real Iq(real q, 
    67     real sld, 
    68     real solvent_sld, 
    69     real radius, 
    70     real cap_radius, 
    71     real length) 
     66double Iq(double q, 
     67    double sld, 
     68    double solvent_sld, 
     69    double radius, 
     70    double cap_radius, 
     71    double length) 
    7272{ 
    73     real sn, cn; // slots to hold sincos function output 
     73    double sn, cn; // slots to hold sincos function output 
    7474 
    7575    // Exclude invalid inputs. 
    76     if (cap_radius < radius) return REAL(-1.0); 
     76    if (cap_radius < radius) return -1.0; 
    7777 
    78     const real lower = REAL(0.0); 
    79     const real upper = M_PI_2; 
    80     const real h = sqrt(cap_radius*cap_radius - radius*radius); // negative h 
    81     real total = REAL(0.0); 
     78    const double lower = 0.0; 
     79    const double upper = M_PI_2; 
     80    const double h = sqrt(cap_radius*cap_radius - radius*radius); // negative h 
     81    double total = 0.0; 
    8282    for (int i=0; i<76 ;i++) { 
    8383        // translate a point in [-1,1] to a point in [lower,upper] 
    84         const real alpha= REAL(0.5)*(Gauss76Z[i]*(upper-lower) + upper + lower); 
     84        const double alpha= 0.5*(Gauss76Z[i]*(upper-lower) + upper + lower); 
    8585        SINCOS(alpha, sn, cn); 
    8686 
    87         const real cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn); 
     87        const double cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn); 
    8888 
    8989        // The following is CylKernel() / sin(alpha), but we are doing it in place 
    9090        // to avoid sin(alpha)/sin(alpha) for alpha = 0.  It is also a teensy bit 
    9191        // faster since we don't multiply and divide sin(alpha). 
    92         const real besarg = q*radius*sn; 
    93         const real siarg = q*REAL(0.5)*length*cn; 
     92        const double besarg = q*radius*sn; 
     93        const double siarg = q*0.5*length*cn; 
    9494        // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1 
    95         const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg); 
    96         const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg); 
    97         const real cyl_Fq = M_PI*radius*radius*length*REAL(2.0)*bj*si; 
     95        const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg); 
     96        const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 
     97        const double cyl_Fq = M_PI*radius*radius*length*2.0*bj*si; 
    9898 
    9999        // Volume weighted average F(q) 
    100         const real Aq = cyl_Fq + cap_Fq; 
     100        const double Aq = cyl_Fq + cap_Fq; 
    101101        total += Gauss76Wt[i] * Aq * Aq * sn; // sn for spherical coord integration 
    102102    } 
    103103    // translate dx in [-1,1] to dx in [lower,upper] 
    104     const real form = total * REAL(0.5)*(upper-lower); 
     104    const double form = total * 0.5*(upper-lower); 
    105105 
    106106    // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1 
    107107    // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 
    108108    // The additional volume factor is for polydisperse volume normalization. 
    109     const real s = (sld - solvent_sld); 
    110     return REAL(1.0e-4) * form * s * s; // form_volume(radius, cap_radius, length); 
     109    const double s = (sld - solvent_sld); 
     110    return 1.0e-4 * form * s * s; // form_volume(radius, cap_radius, length); 
    111111} 
    112112 
    113113 
    114 real Iqxy(real qx, real qy, 
    115     real sld, 
    116     real solvent_sld, 
    117     real radius, 
    118     real cap_radius, 
    119     real length, 
    120     real theta, 
    121     real phi) 
     114double Iqxy(double qx, double qy, 
     115    double sld, 
     116    double solvent_sld, 
     117    double radius, 
     118    double cap_radius, 
     119    double length, 
     120    double theta, 
     121    double phi) 
    122122{ 
    123     real sn, cn; // slots to hold sincos function output 
     123    double sn, cn; // slots to hold sincos function output 
    124124 
    125125    // Exclude invalid inputs. 
    126     if (cap_radius < radius) return REAL(-1.0); 
     126    if (cap_radius < radius) return -1.0; 
    127127 
    128128    // Compute angle alpha between q and the cylinder axis 
    … …  
    130130    // # The following correction factor exists in sasview, but it can't be 
    131131    // # right, so we are leaving it out for now. 
    132     const real q = sqrt(qx*qx+qy*qy); 
    133     const real cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
    134     const real alpha = acos(cos_val); // rod angle relative to q 
     132    const double q = sqrt(qx*qx+qy*qy); 
     133    const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
     134    const double alpha = acos(cos_val); // rod angle relative to q 
    135135    SINCOS(alpha, sn, cn); 
    136136 
    137     const real h = sqrt(cap_radius*cap_radius - radius*radius); // negative h 
    138     const real cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn); 
     137    const double h = sqrt(cap_radius*cap_radius - radius*radius); // negative h 
     138    const double cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn); 
    139139 
    140     const real besarg = q*radius*sn; 
    141     const real siarg = q*REAL(0.5)*length*cn; 
     140    const double besarg = q*radius*sn; 
     141    const double siarg = q*0.5*length*cn; 
    142142    // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1 
    143     const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg); 
    144     const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg); 
    145     const real cyl_Fq = M_PI*radius*radius*length*REAL(2.0)*bj*si; 
     143    const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg); 
     144    const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 
     145    const double cyl_Fq = M_PI*radius*radius*length*2.0*bj*si; 
    146146 
    147147    // Volume weighted average F(q) 
    148     const real Aq = cyl_Fq + cap_Fq; 
     148    const double Aq = cyl_Fq + cap_Fq; 
    149149 
    150150    // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1 
    151     const real s = (sld - solvent_sld); 
    152     return REAL(1.0e-4) * Aq * Aq * s * s; // form_volume(radius, cap_radius, length); 
     151    const double s = (sld - solvent_sld); 
     152    return 1.0e-4 * Aq * Aq * s * s; // form_volume(radius, cap_radius, length); 
    153153} 
Note: See TracChangeset for help on using the changeset viewer.