Changeset c047acf in sasmodels for sasmodels/models/pringle.c


Ignore:
Timestamp:
Jul 27, 2016 3:55:54 AM (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:
1596de3
Parents:
3f8584a2
Message:

much faster pringle model which can run on single precision GPU

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/pringle.c

    r2c74c11 rc047acf  
    1 double form_volume(double radius, 
    2           double thickness); 
     1double form_volume(double radius, double thickness, double alpha, double beta); 
    32 
    43double Iq(double q, 
    … …  
    109          double sld_solvent); 
    1110 
     11 
    1212static 
    13 double pringleC(double radius, 
    14                 double alpha, 
    15                 double beta, 
    16                 double q, 
    17                 double phi, 
    18                 double n) { 
    19  
    20     double va, vb; 
    21     double bessargs, cosarg, bessargcb; 
    22     double r, retval, yyy; 
    23  
    24  
    25     va = 0; 
    26     vb = radius; 
     13void _integrate_bessel( 
     14    double radius, 
     15    double alpha, 
     16    double beta, 
     17    double q_sin_psi, 
     18    double q_cos_psi, 
     19    double n, 
     20    double *Sn, 
     21    double *Cn) 
     22{ 
     23    // translate gauss point z in [-1,1] to a point in [0, radius] 
     24    const double zm = 0.5*radius; 
     25    const double zb = 0.5*radius; 
    2726 
    2827    // evaluate at Gauss points 
    29     // remember to index from 0,size-1 
     28    double sumS = 0.0;          // initialize integral 
     29    double sumC = 0.0;          // initialize integral 
     30    double r; 
     31    for (int i=0; i < 76; i++) { 
     32        r = Gauss76Z[i]*zm + zb; 
    3033 
    31     double summ = 0.0;          // initialize integral 
    32     int ii = 0; 
    33     do { 
    34         // Using 76 Gauss points 
    35         r = (Gauss76Z[ii] * (vb - va) + vb + va) / 2.0; 
     34        const double qrs = r*q_sin_psi; 
     35        const double qrrc = r*r*q_cos_psi; 
    3636 
    37         bessargs = q*r*sin(phi); 
    38         cosarg = q*r*r*alpha*cos(phi); 
    39         bessargcb = q*r*r*beta*cos(phi); 
     37        double y = Gauss76Wt[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs); 
     38        double S, C; 
     39        SINCOS(alpha*qrrc, S, C); 
     40        sumS += y*S; 
     41        sumC += y*C; 
     42    } 
    4043 
    41         yyy = Gauss76Wt[ii]*r*cos(cosarg) 
    42                 *sas_JN(n, bessargcb) 
    43                 *sas_JN(2*n, bessargs); 
    44         summ += yyy; 
    45  
    46         ii += 1; 
    47     } while (ii < N_POINTS_76);                 // end of loop over quadrature points 
    48     // 
    49     // calculate value of integral to return 
    50  
    51     retval = (vb - va) / 2.0 * summ; 
    52     retval = retval / pow(r, 2.0); 
    53  
    54     return retval; 
     44    #if 1 
     45    // TODO: should the normalization be to R^2 or the last value, r^2 
     46    // the gaussian window does not go all the way from 0 to 1. 
     47    //radius = Gauss76Z[75] * zm + zb; 
     48    *Sn = zm*sumS / (r*r); 
     49    *Cn = zm*sumC / (r*r); 
     50    #else 
     51    *Sn = zm*sumS / (radius*radius); 
     52    *Cn = zm*sumC / (radius*radius); 
     53    #endif 
    5554} 
    5655 
    5756static 
    58 double pringleS(double radius, 
    59                 double alpha, 
    60                 double beta, 
    61                 double q, 
    62                 double phi, 
    63                 double n) { 
    64  
    65     double va, vb, summ; 
    66     double bessargs, sinarg, bessargcb; 
    67     double r, retval, yyy; 
    68     // set up the integration 
    69     // end points and weights 
    70  
    71     va = 0; 
    72     vb = radius; 
    73  
    74     // evaluate at Gauss points 
    75     // remember to index from 0,size-1 
    76  
    77     summ = 0.0;         // initialize integral 
    78     int ii = 0; 
    79     do { 
    80         // Using 76 Gauss points 
    81         r = (Gauss76Z[ii] * (vb - va) + vb + va) / 2.0; 
    82  
    83         bessargs = q*r*sin(phi); 
    84         sinarg = q*r*r*alpha*cos(phi); 
    85         bessargcb = q*r*r*beta*cos(phi); 
    86  
    87         yyy = Gauss76Wt[ii]*r*sin(sinarg) 
    88                     *sas_JN(n, bessargcb) 
    89                     *sas_JN(2*n, bessargs); 
    90         summ += yyy; 
    91  
    92         ii += 1; 
    93     } while (ii < N_POINTS_76); 
    94  
    95     // end of loop over quadrature points 
    96     // 
    97     // calculate value of integral to return 
    98  
    99     retval = (vb-va)/2.0*summ; 
    100     retval = retval/pow(r, 2.0); 
    101  
    102     return retval; 
     57double _sum_bessel_orders( 
     58    double radius, 
     59    double alpha, 
     60    double beta, 
     61    double q_sin_psi, 
     62    double q_cos_psi) 
     63{ 
     64    //calculate sum term from n = -3 to 3 
     65    //Note 1: 
     66    //    S_n(-x) = (-1)^S_n(x) 
     67    //    => S_n^2(-x) = S_n^2(x), 
     68    //    => sum_-k^k Sk = S_0^2 + 2*sum_1^kSk^2 
     69    //Note 2: 
     70    //    better precision to sum terms from smaller to larger 
     71    //    though it doesn't seem to make a difference in this case. 
     72    double Sn, Cn, sum; 
     73    sum = 0.0; 
     74    for (int n=3; n>0; n--) { 
     75      _integrate_bessel(radius, alpha, beta, q_sin_psi, q_cos_psi, n, &Sn, &Cn); 
     76      sum += 2.0*(Sn*Sn + Cn*Cn); 
     77    } 
     78    _integrate_bessel(radius, alpha, beta, q_sin_psi, q_cos_psi, 0, &Sn, &Cn); 
     79    sum += Sn*Sn+ Cn*Cn; 
     80    return sum; 
    10381} 
    10482 
    10583static 
    106 double _kernel(double thickness, 
    107                double radius, 
    108                double alpha, 
    109                double beta, 
    110                double q, 
    111                double phi) { 
     84double _integrate_psi( 
     85    double q, 
     86    double radius, 
     87    double thickness, 
     88    double alpha, 
     89    double beta) 
     90{ 
     91    // translate gauss point z in [-1,1] to a point in [0, pi/2] 
     92    const double zm = M_PI_4; 
     93    const double zb = M_PI_4; 
    11294 
    113     const double sincarg = q * thickness * cos(phi) / 2.0; 
    114     const double sincterm = pow(sin(sincarg) / sincarg, 2.0); 
     95    double sum = 0.0; 
     96    for (int i = 0; i < 76; i++) { 
     97        double psi = Gauss76Z[i]*zm + zb; 
     98        double sin_psi, cos_psi; 
     99        SINCOS(psi, sin_psi, cos_psi); 
     100        double bessel_term = _sum_bessel_orders(radius, alpha, beta, q*sin_psi, q*cos_psi); 
     101        double sinc_term = square(sinc(q * thickness * cos_psi / 2.0)); 
     102        double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term; 
     103        sum += Gauss76Wt[i] * pringle_kernel; 
     104    } 
    115105 
    116     //calculate sum term from n = -3 to 3 
    117     double sumterm = 0.0; 
    118     for (int nn = -3; nn <= 3; nn++) { 
    119         double powc = pringleC(radius, alpha, beta, q, phi, nn); 
    120         double pows = pringleS(radius, alpha, beta, q, phi, nn); 
    121         sumterm += pow(powc, 2.0) + pow(pows, 2.0); 
    122     } 
    123     double retval = 4.0 * sin(phi) * sumterm * sincterm; 
    124  
    125     return retval; 
    126  
     106    return zm * sum; 
    127107} 
    128108 
    129 static double pringles_kernel(double q, 
    130           double radius, 
    131           double thickness, 
    132           double alpha, 
    133           double beta, 
    134           double sld_pringle, 
    135           double sld_solvent) 
     109double form_volume(double radius, double thickness, double alpha, double beta) 
    136110{ 
    137  
    138     //upper and lower integration limits 
    139     const double lolim = 0.0; 
    140     const double uplim = M_PI / 2.0; 
    141  
    142     double summ = 0.0;                  //initialize integral 
    143  
    144     double delrho = sld_pringle - sld_solvent; //make contrast term 
    145  
    146     for (int i = 0; i < N_POINTS_76; i++) { 
    147         double phi = (Gauss76Z[i] * (uplim - lolim) + uplim + lolim) / 2.0; 
    148         summ += Gauss76Wt[i] * _kernel(thickness, radius, alpha, beta, q, phi); 
    149     } 
    150  
    151     double answer = (uplim - lolim) / 2.0 * summ; 
    152     answer *= delrho*delrho; 
    153  
    154     return answer; 
     111    // TODO: Normalize by form volume 
     112    //return M_PI*radius*radius*thickness; 
     113    return 1.0; 
    155114} 
    156115 
    157 double form_volume(double radius, 
    158         double thickness){ 
    159  
    160         return 1.0; 
     116double Iq( 
     117    double q, 
     118    double radius, 
     119    double thickness, 
     120    double alpha, 
     121    double beta, 
     122    double sld_pringle, 
     123    double sld_solvent) 
     124{ 
     125    double form = _integrate_psi(q, radius, thickness, alpha, beta); 
     126    double contrast = sld_pringle - sld_solvent; 
     127    double volume = M_PI*radius*radius*thickness; 
     128    // TODO: If normalize by form volume, need an extra volume here 
     129    //return 1.0e-4*form * square(contrast * volume); 
     130    return 1.0e-4*form * square(contrast) * volume; 
    161131} 
    162  
    163 double Iq(double q, 
    164           double radius, 
    165           double thickness, 
    166           double alpha, 
    167           double beta, 
    168           double sld_pringle, 
    169           double sld_solvent) 
    170 { 
    171     const double form = pringles_kernel(q, 
    172                   radius, 
    173                   thickness, 
    174                   alpha, 
    175                   beta, 
    176                   sld_pringle, 
    177                   sld_solvent); 
    178  
    179     return 1.0e-4*form*M_PI*radius*radius*thickness; 
    180 } 
Note: See TracChangeset for help on using the changeset viewer.