Changeset 994d77f in sasmodels


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

Location:
sasmodels
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/gen.py

    r78356b31 r994d77f  
    3737order so that functions are defined before they are used. 
    3838 
    39 Floating point values should be declared as *real*.  Depending on how the 
    40 function is called, a macro will replace *real* with *float* or *double*. 
    41 Unfortunately, MacOSX is picky about floating point constants, which 
    42 should be defined with value + 'f' if they are of type *float* or just 
    43 a bare value if they are of type *double*.  The solution is a macro 
    44 *REAL(value)* which adds the 'f' if compiling for single precision 
    45 floating point.  [Note: we could write a clever regular expression 
    46 which automatically detects real-valued constants.  If we wanted to 
    47 make our code more C-like, we could define variables with double but 
    48 replace them with float before compiling for single precision.] 
     39Floating point values should be declared as *double*.  For single precision 
     40calculations, *double* will be replaced by *float*.  The single precision 
     41conversion will also tag floating point constants with "f" to make them 
     42single precision constants.  When using integral values in floating point 
     43expressions, they should be expressed as floating point values by including 
     44a decimal point.  This includes 0., 1. and 2. 
    4945 
    5046OpenCL has a *sincos* function which can improve performance when both 
    … …  
    5248this function does not exist in C99, all use of *sincos* should be 
    5349replaced by the macro *SINCOS(value,sn,cn)* where *sn* and *cn* are 
    54 previously declared *real* values.  *value* may be an expression.  When 
    55 compiled for systems without OpenCL, *SINCOS* will be replaced by 
    56 *sin* and *cos* calls. 
     50previously declared *double* variables.  When compiled for systems without 
     51OpenCL, *SINCOS* will be replaced by *sin* and *cos* calls.   If *value* is 
     52an expression, it will appear twice in this case; whether or not it will be 
     53evaluated twice depends on the quality of the compiler. 
    5754 
    5855If the input parameters are invalid, the scattering calculator should 
    … …  
    172169# TODO: identify model files which have changed since loading and reload them. 
    173170 
    174 __all__ = ["make, doc", "sources"] 
     171__all__ = ["make, doc", "sources", "use_single"] 
    175172 
    176173import sys 
    177174import os 
    178175import os.path 
     176import re 
    179177 
    180178import numpy as np 
    … …  
    229227#ifndef USE_OPENCL 
    230228#  include <math.h> 
    231 #  define REAL(x) (x) 
    232 #  ifndef real 
    233 #      define real double 
    234 #  endif 
    235229#  define global 
    236230#  define local 
    … …  
    253247// is not enabled on the card, which is why these constants may be missing 
    254248#ifndef M_PI 
    255 #  define M_PI REAL(3.141592653589793) 
     249#  define M_PI 3.141592653589793 
    256250#endif 
    257251#ifndef M_PI_2 
    258 #  define M_PI_2 REAL(1.570796326794897) 
     252#  define M_PI_2 1.570796326794897 
    259253#endif 
    260254#ifndef M_PI_4 
    261 #  define M_PI_4 REAL(0.7853981633974483) 
     255#  define M_PI_4 0.7853981633974483 
    262256#endif 
    263257 
    264258// Non-standard pi/180, used for converting between degrees and radians 
    265259#ifndef M_PI_180 
    266 #  define M_PI_180 REAL(0.017453292519943295) 
     260#  define M_PI_180 0.017453292519943295 
    267261#endif 
    268262""" 
    … …  
    274268KERNEL_1D = { 
    275269    'fn': "Iq", 
    276     'q_par_decl': "global const real *q,", 
    277     'qinit': "const real qi = q[i];", 
     270    'q_par_decl': "global const double *q,", 
     271    'qinit': "const double qi = q[i];", 
    278272    'qcall': "qi", 
    279273    'qwork': ["q"], 
    … …  
    282276KERNEL_2D = { 
    283277    'fn': "Iqxy", 
    284     'q_par_decl': "global const real *qx,\n    global const real *qy,", 
    285     'qinit': "const real qxi = qx[i];\n    const real qyi = qy[i];", 
     278    'q_par_decl': "global const double *qx,\n    global const double *qy,", 
     279    'qinit': "const double qxi = qx[i];\n    const double qyi = qy[i];", 
    286280    'qcall': "qxi, qyi", 
    287281    'qwork': ["qx", "qy"], 
    … …  
    296290kernel void %(name)s( 
    297291    %(q_par_decl)s 
    298     global real *result, 
     292    global double *result, 
    299293#ifdef USE_OPENCL 
    300     global real *loops_g, 
     294    global double *loops_g, 
    301295#else 
    302296    const int Nq, 
    303297#endif 
    304     local real *loops, 
    305     const real cutoff, 
     298    local double *loops, 
     299    const double cutoff, 
    306300    %(par_decl)s 
    307301    ) 
    … …  
    324318  { 
    325319    %(qinit)s 
    326     real ret=REAL(0.0), norm=REAL(0.0); 
    327     real vol=REAL(0.0), norm_vol=REAL(0.0); 
     320    double ret=0.0, norm=0.0; 
     321    double vol=0.0, norm_vol=0.0; 
    328322%(loops)s 
    329     if (vol*norm_vol != REAL(0.0)) { 
     323    if (vol*norm_vol != 0.0) { 
    330324      ret *= norm_vol/vol; 
    331325    } 
    … …  
    340334LOOP_OPEN="""\ 
    341335for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) { 
    342   const real %(name)s = loops[2*(%(name)s_i%(offset)s)]; 
    343   const real %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\ 
     336  const double %(name)s = loops[2*(%(name)s_i%(offset)s)]; 
     337  const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\ 
    344338""" 
    345339 
    … …  
    349343# it will also be added here. 
    350344LOOP_BODY="""\ 
    351 const real weight = %(weight_product)s; 
     345const double weight = %(weight_product)s; 
    352346if (weight > cutoff) { 
    353   const real I = %(fn)s(%(qcall)s, %(pcall)s); 
    354   if (I>=REAL(0.0)) { // scattering cannot be negative 
     347  const double I = %(fn)s(%(qcall)s, %(pcall)s); 
     348  if (I>=0.0) { // scattering cannot be negative 
    355349    ret += weight*I%(sasview_spherical)s; 
    356350    norm += weight; 
    … …  
    365359SPHERICAL_CORRECTION="""\ 
    366360// Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta 
    367 real spherical_correction = (Ntheta>1 ? fabs(sin(M_PI_180*theta)) : REAL(1.0));\ 
     361double spherical_correction = (Ntheta>1 ? fabs(sin(M_PI_180*theta)) : 1.0);\ 
    368362""" 
    369363# Use this to reproduce sasview behaviour 
    370364SASVIEW_SPHERICAL_CORRECTION="""\ 
    371365// Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta 
    372 real spherical_correction = (Ntheta>1 ? fabs(cos(M_PI_180*theta))*M_PI_2 : REAL(1.0));\ 
     366double spherical_correction = (Ntheta>1 ? fabs(cos(M_PI_180*theta))*M_PI_2 : 1.0);\ 
    373367""" 
    374368 
    … …  
    377371# to call the form_volume function from the user supplied kernel, and accumulate 
    378372# a normalized weight. 
    379 VOLUME_NORM="""const real vol_weight = %(weight)s; 
     373VOLUME_NORM="""const double vol_weight = %(weight)s; 
    380374    vol += vol_weight*form_volume(%(pars)s); 
    381375    norm_vol += vol_weight;\ 
    … …  
    384378# functions defined as strings in the .py module 
    385379WORK_FUNCTION="""\ 
    386 real %(name)s(%(pars)s); 
    387 real %(name)s(%(pars)s) 
     380double %(name)s(%(pars)s); 
     381double %(name)s(%(pars)s) 
    388382{ 
    389383%(body)s 
    … …  
    422416    """ 
    423417    return info['name'] + "_" + ("Iqxy" if is_2D else "Iq") 
     418 
     419 
     420def use_single(source): 
     421    """ 
     422    Convert code from double precision to single precision. 
     423    """ 
     424    source = re.sub(r'(^|[^a-zA-Z0-9_])double($|[^a-zA-Z0-9_])', 
     425                    r'\1float\2', source) 
     426    source = re.sub(r'[^a-zA-Z_](\d*[.]\d+|\d+[.]\d*)([eE][+-]?\d+)?', 
     427                    r'\g<0>f', source) 
     428    return source 
    424429 
    425430 
    … …  
    503508    # declarations for non-pd followed by pd pars 
    504509    # e.g., 
    505     #     const real sld, 
     510    #     const double sld, 
    506511    #     const int Nradius 
    507     fixed_par_decl = ",\n    ".join("const real %s"%p for p in fixed_pars) 
     512    fixed_par_decl = ",\n    ".join("const double %s"%p for p in fixed_pars) 
    508513    pd_par_decl = ",\n    ".join("const int N%s"%p for p in pd_pars) 
    509514    if fixed_par_decl and pd_par_decl: 
    … …  
    518523        # kernel name is, e.g., cylinder_Iq 
    519524        'name': kernel_name(info, is_2D), 
    520         # to declare, e.g., global real q[], 
     525        # to declare, e.g., global double q[], 
    521526        'q_par_decl': q_pars['q_par_decl'], 
    522         # to declare, e.g., real sld, int Nradius, int Nlength 
     527        # to declare, e.g., double sld, int Nradius, int Nlength 
    523528        'par_decl': par_decl, 
    524529        # to copy global to local pd pars we need, e.g., Nradius+Nlength 
    525530        'pd_length': "+".join('N'+p for p in pd_pars), 
    526         # the q initializers, e.g., real qi = q[i]; 
     531        # the q initializers, e.g., double qi = q[i]; 
    527532        'qinit': q_pars['qinit'], 
    528533        # the actual polydispersity loop 
    … …  
    537542        subst = { 
    538543            'name': fn, 
    539             'pars': ", ".join("real "+p for p in q_pars['qwork']+fq_pars), 
     544            'pars': ", ".join("double "+p for p in q_pars['qwork']+fq_pars), 
    540545            'body': info[fn], 
    541546            } 
    … …  
    607612        subst = { 
    608613            'name': "form_volume", 
    609             'pars': ", ".join("real "+p for p in info['partype']['volume']), 
     614            'pars': ", ".join("double "+p for p in info['partype']['volume']), 
    610615            'body': info['form_volume'], 
    611616            } 

  • sasmodels/gpu.py

    rd087487b r994d77f  
    2929from . import gen 
    3030 
    31 F32_DEFS = """\ 
    32 #define REAL(x) (x##f) 
    33 #define real float 
    34 """ 
    35  
    36  
    3731F64_DEFS = """\ 
    3832#ifdef cl_khr_fp64 
    3933#  pragma OPENCL EXTENSION cl_khr_fp64: enable 
    4034#endif 
    41 #define REAL(x) (x) 
    42 #define real double 
    4335""" 
    4436 
    … …  
    116108        raise RuntimeError("Double precision not supported for devices") 
    117109 
    118     header = F64_DEFS if dtype == gen.F64 else F32_DEFS 
     110    header = F64_DEFS if dtype == gen.F64 else "" 
     111    if dtype == gen.F32: 
     112        source = gen.use_single(source) 
    119113    # Note: USE_SINCOS makes the intel cpu slower under opencl 
    120114    if context.devices[0].type == cl.device_type.GPU: 

  • 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} 

  • sasmodels/models/core_shell_cylinder.c

    r5d4777d r994d77f  
    1 real form_volume(real radius, real thickness, real length); 
    2 real Iq(real q, real core_sld, real shell_sld, real solvent_sld, 
    3     real radius, real thickness, real length); 
    4 real Iqxy(real qx, real qy, real core_sld, real shell_sld, real solvent_sld, 
    5     real radius, real thickness, real length, real theta, real phi); 
     1double form_volume(double radius, double thickness, double length); 
     2double Iq(double q, double core_sld, double shell_sld, double solvent_sld, 
     3    double radius, double thickness, double length); 
     4double Iqxy(double qx, double qy, double core_sld, double shell_sld, double solvent_sld, 
     5    double radius, double thickness, double length, double theta, double phi); 
    66 
    77// twovd = 2 * volume * delta_rho 
    88// besarg = q * R * sin(alpha) 
    99// siarg = q * L/2 * cos(alpha) 
    10 real _cyl(real twovd, real besarg, real siarg); 
    11 real _cyl(real twovd, real besarg, real siarg) 
     10double _cyl(double twovd, double besarg, double siarg); 
     11double _cyl(double twovd, double besarg, double siarg) 
    1212{ 
    13     const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg); 
    14     const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg); 
     13    const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg); 
     14    const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 
    1515    return twovd*si*bj; 
    1616} 
    1717 
    18 real form_volume(real radius, real thickness, real length) 
     18double form_volume(double radius, double thickness, double length) 
    1919{ 
    2020    return M_PI*(radius+thickness)*(radius+thickness)*(length+2*thickness); 
    2121} 
    2222 
    23 real Iq(real q, 
    24     real core_sld, 
    25     real shell_sld, 
    26     real solvent_sld, 
    27     real radius, 
    28     real thickness, 
    29     real length) 
     23double Iq(double q, 
     24    double core_sld, 
     25    double shell_sld, 
     26    double solvent_sld, 
     27    double radius, 
     28    double thickness, 
     29    double length) 
    3030{ 
    3131    // precalculate constants 
    32     const real core_qr = q*radius; 
    33     const real core_qh = q*REAL(0.5)*length; 
    34     const real core_twovd = REAL(2.0) * form_volume(radius,0,length) 
     32    const double core_qr = q*radius; 
     33    const double core_qh = q*0.5*length; 
     34    const double core_twovd = 2.0 * form_volume(radius,0,length) 
    3535                            * (core_sld-shell_sld); 
    36     const real shell_qr = q*(radius + thickness); 
    37     const real shell_qh = q*(REAL(0.5)*length + thickness); 
    38     const real shell_twovd = REAL(2.0) * form_volume(radius,thickness,length) 
     36    const double shell_qr = q*(radius + thickness); 
     37    const double shell_qh = q*(0.5*length + thickness); 
     38    const double shell_twovd = 2.0 * form_volume(radius,thickness,length) 
    3939                             * (shell_sld-solvent_sld); 
    40     real total = REAL(0.0); 
    41     // real lower=0, upper=M_PI_2; 
     40    double total = 0.0; 
     41    // double lower=0, upper=M_PI_2; 
    4242    for (int i=0; i<76 ;i++) { 
    4343        // translate a point in [-1,1] to a point in [lower,upper] 
    44         //const real alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    45         real sn, cn; 
    46         const real alpha = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
     44        //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
     45        double sn, cn; 
     46        const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
    4747        SINCOS(alpha, sn, cn); 
    48         const real fq = _cyl(core_twovd, core_qr*sn, core_qh*cn) 
     48        const double fq = _cyl(core_twovd, core_qr*sn, core_qh*cn) 
    4949            + _cyl(shell_twovd, shell_qr*sn, shell_qh*cn); 
    5050        total += Gauss76Wt[i] * fq * fq * sn; 
    5151    } 
    5252    // translate dx in [-1,1] to dx in [lower,upper] 
    53     //const real form = (upper-lower)/2.0*total; 
    54     return REAL(1.0e-4) * total * M_PI_4; 
     53    //const double form = (upper-lower)/2.0*total; 
     54    return 1.0e-4 * total * M_PI_4; 
    5555} 
    5656 
    5757 
    58 real Iqxy(real qx, real qy, 
    59     real core_sld, 
    60     real shell_sld, 
    61     real solvent_sld, 
    62     real radius, 
    63     real thickness, 
    64     real length, 
    65     real theta, 
    66     real phi) 
     58double Iqxy(double qx, double qy, 
     59    double core_sld, 
     60    double shell_sld, 
     61    double solvent_sld, 
     62    double radius, 
     63    double thickness, 
     64    double length, 
     65    double theta, 
     66    double phi) 
    6767{ 
    68     real sn, cn; // slots to hold sincos function output 
     68    double sn, cn; // slots to hold sincos function output 
    6969 
    7070    // Compute angle alpha between q and the cylinder axis 
    … …  
    7272    // # The following correction factor exists in sasview, but it can't be 
    7373    // # right, so we are leaving it out for now. 
    74     // const real correction = fabs(cn)*M_PI_2; 
    75     const real q = sqrt(qx*qx+qy*qy); 
    76     const real cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
    77     const real alpha = acos(cos_val); 
     74    // const double correction = fabs(cn)*M_PI_2; 
     75    const double q = sqrt(qx*qx+qy*qy); 
     76    const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
     77    const double alpha = acos(cos_val); 
    7878 
    79     const real core_qr = q*radius; 
    80     const real core_qh = q*REAL(0.5)*length; 
    81     const real core_twovd = REAL(2.0) * form_volume(radius,0,length) 
     79    const double core_qr = q*radius; 
     80    const double core_qh = q*0.5*length; 
     81    const double core_twovd = 2.0 * form_volume(radius,0,length) 
    8282                            * (core_sld-shell_sld); 
    83     const real shell_qr = q*(radius + thickness); 
    84     const real shell_qh = q*(REAL(0.5)*length + thickness); 
    85     const real shell_twovd = REAL(2.0) * form_volume(radius,thickness,length) 
     83    const double shell_qr = q*(radius + thickness); 
     84    const double shell_qh = q*(0.5*length + thickness); 
     85    const double shell_twovd = 2.0 * form_volume(radius,thickness,length) 
    8686                             * (shell_sld-solvent_sld); 
    8787 
    8888    SINCOS(alpha, sn, cn); 
    89     const real fq = _cyl(core_twovd, core_qr*sn, core_qh*cn) 
     89    const double fq = _cyl(core_twovd, core_qr*sn, core_qh*cn) 
    9090        + _cyl(shell_twovd, shell_qr*sn, shell_qh*cn); 
    91     return REAL(1.0e-4) * fq * fq; 
     91    return 1.0e-4 * fq * fq; 
    9292} 

  • sasmodels/models/cylinder.c

    r0496031 r994d77f  
    1 real form_volume(real radius, real length); 
    2 real Iq(real q, real sld, real solvent_sld, real radius, real length); 
    3 real Iqxy(real qx, real qy, real sld, real solvent_sld, 
    4     real radius, real length, real theta, real phi); 
     1double form_volume(double radius, double length); 
     2double Iq(double q, double sld, double solvent_sld, double radius, double length); 
     3double Iqxy(double qx, double qy, double sld, double solvent_sld, 
     4    double radius, double length, double theta, double phi); 
    55 
    66// twovd = 2 * volume * delta_rho 
    77// besarg = q * R * sin(alpha) 
    88// siarg = q * L/2 * cos(alpha) 
    9 real _cyl(real twovd, real besarg, real siarg); 
    10 real _cyl(real twovd, real besarg, real siarg) 
     9double _cyl(double twovd, double besarg, double siarg); 
     10double _cyl(double twovd, double besarg, double siarg) 
    1111{ 
    12     const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg); 
    13     const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg); 
     12    const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg); 
     13    const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 
    1414    return twovd*si*bj; 
    1515} 
    1616 
    17 real form_volume(real radius, real length) 
     17double form_volume(double radius, double length) 
    1818{ 
    1919    return M_PI*radius*radius*length; 
    2020} 
    2121 
    22 real Iq(real q, 
    23     real sld, 
    24     real solvent_sld, 
    25     real radius, 
    26     real length) 
     22double Iq(double q, 
     23    double sld, 
     24    double solvent_sld, 
     25    double radius, 
     26    double length) 
    2727{ 
    28     const real qr = q*radius; 
    29     const real qh = q*REAL(0.5)*length; 
    30     const real twovd = REAL(2.0)*(sld-solvent_sld)*form_volume(radius, length); 
    31     real total = REAL(0.0); 
    32     // real lower=0, upper=M_PI_2; 
     28    const double qr = q*radius; 
     29    const double qh = q*0.5*length; 
     30    const double twovd = 2.0*(sld-solvent_sld)*form_volume(radius, length); 
     31    double total = 0.0; 
     32    // double lower=0, upper=M_PI_2; 
    3333    for (int i=0; i<76 ;i++) { 
    3434        // translate a point in [-1,1] to a point in [lower,upper] 
    35         //const real alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    36         const real alpha = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
    37         real sn, cn; 
     35        //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
     36        const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
     37        double sn, cn; 
    3838        SINCOS(alpha, sn, cn); 
    39         const real fq = _cyl(twovd, qr*sn, qh*cn); 
     39        const double fq = _cyl(twovd, qr*sn, qh*cn); 
    4040        total += Gauss76Wt[i] * fq * fq * sn; 
    4141    } 
    4242    // translate dx in [-1,1] to dx in [lower,upper] 
    43     //const real form = (upper-lower)/2.0*total; 
    44     return REAL(1.0e-4) * total * M_PI_4; 
     43    //const double form = (upper-lower)/2.0*total; 
     44    return 1.0e-4 * total * M_PI_4; 
    4545} 
    4646 
    4747 
    48 real Iqxy(real qx, real qy, 
    49     real sld, 
    50     real solvent_sld, 
    51     real radius, 
    52     real length, 
    53     real theta, 
    54     real phi) 
     48double Iqxy(double qx, double qy, 
     49    double sld, 
     50    double solvent_sld, 
     51    double radius, 
     52    double length, 
     53    double theta, 
     54    double phi) 
    5555{ 
    5656    // TODO: check that radius<0 and length<0 give zero scattering. 
    5757    // This should be the case since the polydispersity weight vector should 
    5858    // be zero length, and this function never called. 
    59     real sn, cn; // slots to hold sincos function output 
     59    double sn, cn; // slots to hold sincos function output 
    6060 
    6161    // Compute angle alpha between q and the cylinder axis 
    6262    SINCOS(theta*M_PI_180, sn, cn); 
    63     const real q = sqrt(qx*qx+qy*qy); 
    64     const real cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
    65     const real alpha = acos(cos_val); 
     63    const double q = sqrt(qx*qx+qy*qy); 
     64    const double cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
     65    const double alpha = acos(cos_val); 
    6666 
    67     const real twovd = REAL(2.0)*(sld-solvent_sld)*form_volume(radius, length); 
     67    const double twovd = 2.0*(sld-solvent_sld)*form_volume(radius, length); 
    6868    SINCOS(alpha, sn, cn); 
    69     const real fq = _cyl(twovd, q*radius*sn, q*REAL(0.5)*length*cn); 
    70     return REAL(1.0e-4) * fq * fq; 
     69    const double fq = _cyl(twovd, q*radius*sn, q*0.5*length*cn); 
     70    return 1.0e-4 * fq * fq; 
    7171} 

  • sasmodels/models/cylinder_clone.c

    r5d4777d r994d77f  
    1 real form_volume(real radius, real length); 
    2 real Iq(real q, real sld, real solvent_sld, real radius, real length); 
    3 real Iqxy(real qx, real qy, real sld, real solvent_sld, real radius, real length, real theta, real phi); 
     1double form_volume(double radius, double length); 
     2double Iq(double q, double sld, double solvent_sld, double radius, double length); 
     3double Iqxy(double qx, double qy, double sld, double solvent_sld, double radius, double length, double theta, double phi); 
    44 
    55 
    … …  
    77// besarg = q * R * sin(alpha) 
    88// siarg = q * L/2 * cos(alpha) 
    9 real _cyl(real twovd, real besarg, real siarg, real alpha); 
    10 real _cyl(real twovd, real besarg, real siarg, real alpha) 
     9double _cyl(double twovd, double besarg, double siarg, double alpha); 
     10double _cyl(double twovd, double besarg, double siarg, double alpha) 
    1111{ 
    12     const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg); 
    13     const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg); 
     12    const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg); 
     13    const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 
    1414    return twovd*si*bj; 
    1515} 
    1616 
    17 real form_volume(real radius, real length) 
     17double form_volume(double radius, double length) 
    1818{ 
    1919    return M_PI*radius*radius*length; 
    2020} 
    21 real Iq(real q, 
    22     real sldCyl, 
    23     real sldSolv, 
    24     real radius, 
    25     real length) 
     21double Iq(double q, 
     22    double sldCyl, 
     23    double sldSolv, 
     24    double radius, 
     25    double length) 
    2626{ 
    27     const real qr = q*radius; 
    28     const real qh = q*REAL(0.5)*length; 
    29     const real twovd = REAL(2.0)*(sldCyl-sldSolv)*form_volume(radius, length); 
    30     real total = REAL(0.0); 
    31     // real lower=0, upper=M_PI_2; 
     27    const double qr = q*radius; 
     28    const double qh = q*0.5*length; 
     29    const double twovd = 2.0*(sldCyl-sldSolv)*form_volume(radius, length); 
     30    double total = 0.0; 
     31    // double lower=0, upper=M_PI_2; 
    3232    for (int i=0; i<76 ;i++) { 
    3333        // translate a point in [-1,1] to a point in [lower,upper] 
    34         //const real alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    35         const real alpha = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
    36         real sn, cn; 
     34        //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
     35        const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
     36        double sn, cn; 
    3737        SINCOS(alpha, sn, cn); 
    38         const real fq = _cyl(twovd, qr*sn, qh*cn, alpha); 
     38        const double fq = _cyl(twovd, qr*sn, qh*cn, alpha); 
    3939        total += Gauss76Wt[i] * fq * fq * sn; 
    4040    } 
    4141    // translate dx in [-1,1] to dx in [lower,upper] 
    42     //const real form = (upper-lower)/2.0*total; 
    43     return REAL(1.0e8) * total * M_PI_4; 
     42    //const double form = (upper-lower)/2.0*total; 
     43    return 1.0e8 * total * M_PI_4; 
    4444} 
    4545 
    46 real Iqxy(real qx, real qy, 
    47     real sldCyl, 
    48     real sldSolv, 
    49     real radius, 
    50     real length, 
    51     real cyl_theta, 
    52     real cyl_phi) 
     46double Iqxy(double qx, double qy, 
     47    double sldCyl, 
     48    double sldSolv, 
     49    double radius, 
     50    double length, 
     51    double cyl_theta, 
     52    double cyl_phi) 
    5353{ 
    54     real sn, cn; // slots to hold sincos function output 
     54    double sn, cn; // slots to hold sincos function output 
    5555 
    5656    // Compute angle alpha between q and the cylinder axis 
    … …  
    5858    // # The following correction factor exists in sasview, but it can't be 
    5959    // # right, so we are leaving it out for now. 
    60     const real spherical_integration = fabs(cn)*M_PI_2; 
    61     const real q = sqrt(qx*qx+qy*qy); 
    62     const real cos_val = cn*cos(cyl_phi*M_PI_180)*(qx/q) + sn*(qy/q); 
    63     const real alpha = acos(cos_val); 
     60    const double spherical_integration = fabs(cn)*M_PI_2; 
     61    const double q = sqrt(qx*qx+qy*qy); 
     62    const double cos_val = cn*cos(cyl_phi*M_PI_180)*(qx/q) + sn*(qy/q); 
     63    const double alpha = acos(cos_val); 
    6464 
    65     const real qr = q*radius; 
    66     const real qh = q*REAL(0.5)*length; 
    67     const real twovd = REAL(2.0)*(sldCyl-sldSolv)*form_volume(radius, length); 
     65    const double qr = q*radius; 
     66    const double qh = q*0.5*length; 
     67    const double twovd = 2.0*(sldCyl-sldSolv)*form_volume(radius, length); 
    6868    SINCOS(alpha, sn, cn); 
    69     const real fq = _cyl(twovd, qr*sn, qh*cn, alpha); 
    70     return REAL(1.0e8) * fq * fq * spherical_integration; 
     69    const double fq = _cyl(twovd, qr*sn, qh*cn, alpha); 
     70    return 1.0e8 * fq * fq * spherical_integration; 
    7171} 

  • sasmodels/models/ellipsoid.c

    r5d4777d r994d77f  
    1 real form_volume(real rpolar, real requatorial); 
    2 real Iq(real q, real sld, real solvent_sld, real rpolar, real requatorial); 
    3 real Iqxy(real qx, real qy, real sld, real solvent_sld, 
    4     real rpolar, real requatorial, real theta, real phi); 
     1double form_volume(double rpolar, double requatorial); 
     2double Iq(double q, double sld, double solvent_sld, double rpolar, double requatorial); 
     3double Iqxy(double qx, double qy, double sld, double solvent_sld, 
     4    double rpolar, double requatorial, double theta, double phi); 
    55 
    6 real _ellipsoid_kernel(real q, real rpolar, real requatorial, real cos_alpha); 
    7 real _ellipsoid_kernel(real q, real rpolar, real requatorial, real cos_alpha) 
     6double _ellipsoid_kernel(double q, double rpolar, double requatorial, double cos_alpha); 
     7double _ellipsoid_kernel(double q, double rpolar, double requatorial, double cos_alpha) 
    88{ 
    9     real sn, cn; 
    10     real ratio = rpolar/requatorial; 
    11     const real u = q*requatorial*sqrt(REAL(1.0) 
    12                    + cos_alpha*cos_alpha*(ratio*ratio - REAL(1.0))); 
     9    double sn, cn; 
     10    double ratio = rpolar/requatorial; 
     11    const double u = q*requatorial*sqrt(1.0 
     12                   + cos_alpha*cos_alpha*(ratio*ratio - 1.0)); 
    1313    SINCOS(u, sn, cn); 
    14     const real f = ( u==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(sn-u*cn)/(u*u*u) ); 
     14    const double f = ( u==0.0 ? 1.0 : 3.0*(sn-u*cn)/(u*u*u) ); 
    1515    return f*f; 
    1616} 
    1717 
    18 real form_volume(real rpolar, real requatorial) 
     18double form_volume(double rpolar, double requatorial) 
    1919{ 
    20     return REAL(1.333333333333333)*M_PI*rpolar*requatorial*requatorial; 
     20    return 1.333333333333333*M_PI*rpolar*requatorial*requatorial; 
    2121} 
    2222 
    23 real Iq(real q, 
    24     real sld, 
    25     real solvent_sld, 
    26     real rpolar, 
    27     real requatorial) 
     23double Iq(double q, 
     24    double sld, 
     25    double solvent_sld, 
     26    double rpolar, 
     27    double requatorial) 
    2828{ 
    29     //const real lower = REAL(0.0); 
    30     //const real upper = REAL(1.0); 
    31     real total = REAL(0.0); 
     29    //const double lower = 0.0; 
     30    //const double upper = 1.0; 
     31    double total = 0.0; 
    3232    for (int i=0;i<76;i++) { 
    33         //const real cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 
    34         const real cos_alpha = REAL(0.5)*(Gauss76Z[i] + REAL(1.0)); 
     33        //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 
     34        const double cos_alpha = 0.5*(Gauss76Z[i] + 1.0); 
    3535        total += Gauss76Wt[i] * _ellipsoid_kernel(q, rpolar, requatorial, cos_alpha); 
    3636    } 
    37     //const real form = (upper-lower)/2*total; 
    38     const real form = REAL(0.5)*total; 
    39     const real s = (sld - solvent_sld) * form_volume(rpolar, requatorial); 
    40     return REAL(1.0e-4) * form * s * s; 
     37    //const double form = (upper-lower)/2*total; 
     38    const double form = 0.5*total; 
     39    const double s = (sld - solvent_sld) * form_volume(rpolar, requatorial); 
     40    return 1.0e-4 * form * s * s; 
    4141} 
    4242 
    43 real Iqxy(real qx, real qy, 
    44     real sld, 
    45     real solvent_sld, 
    46     real rpolar, 
    47     real requatorial, 
    48     real theta, 
    49     real phi) 
     43double Iqxy(double qx, double qy, 
     44    double sld, 
     45    double solvent_sld, 
     46    double rpolar, 
     47    double requatorial, 
     48    double theta, 
     49    double phi) 
    5050{ 
    51     real sn, cn; 
     51    double sn, cn; 
    5252 
    53     const real q = sqrt(qx*qx + qy*qy); 
     53    const double q = sqrt(qx*qx + qy*qy); 
    5454    SINCOS(theta*M_PI_180, sn, cn); 
    55     const real cos_alpha = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
    56     const real form = _ellipsoid_kernel(q, rpolar, requatorial, cos_alpha); 
    57     const real s = (sld - solvent_sld) * form_volume(rpolar, requatorial); 
     55    const double cos_alpha = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
     56    const double form = _ellipsoid_kernel(q, rpolar, requatorial, cos_alpha); 
     57    const double s = (sld - solvent_sld) * form_volume(rpolar, requatorial); 
    5858 
    59     return REAL(1.0e-4) * form * s * s; 
     59    return 1.0e-4 * form * s * s; 
    6060} 
    6161 

  • sasmodels/models/lamellar.py

    r19dcb933 r994d77f  
    7777# This should perhaps be volume normalized? 
    7878form_volume = """ 
    79     return REAL(1.0); 
     79    return 1.0; 
    8080    """ 
    8181 
    8282Iq = """ 
    83     const real sub = sld - solvent_sld; 
    84     const real qsq = q*q; 
    85     return REAL(4.0e-4)*M_PI*sub*sub/qsq * (REAL(1.0)-cos(q*thickness)) 
     83    const double sub = sld - solvent_sld; 
     84    const double qsq = q*q; 
     85    return 4.0e-4*M_PI*sub*sub/qsq * (1.0-cos(q*thickness)) 
    8686        / (thickness*qsq); 
    8787    """ 
    … …  
    8989Iqxy = """ 
    9090    // never called since no orientation or magnetic parameters. 
    91     return REAL(-1.0); 
     91    return -1.0; 
    9292    """ 
    9393 

  • sasmodels/models/lib/J1.c

    r14de349 r994d77f  
    1 real J1(real x); 
    2 real J1(real x) 
     1double J1(double x); 
     2double J1(double x) 
    33{ 
    4   const real ax = fabs(x); 
    5   if (ax < REAL(8.0)) { 
    6     const real y = x*x; 
    7     const real ans1 = x*(REAL(72362614232.0) 
    8               + y*(REAL(-7895059235.0) 
    9               + y*(REAL(242396853.1) 
    10               + y*(REAL(-2972611.439) 
    11               + y*(REAL(15704.48260) 
    12               + y*(REAL(-30.16036606))))))); 
    13     const real ans2 = REAL(144725228442.0) 
    14               + y*(REAL(2300535178.0) 
    15               + y*(REAL(18583304.74) 
    16               + y*(REAL(99447.43394) 
    17               + y*(REAL(376.9991397) 
     4  const double ax = fabs(x); 
     5  if (ax < 8.0) { 
     6    const double y = x*x; 
     7    const double ans1 = x*(72362614232.0 
     8              + y*(-7895059235.0 
     9              + y*(242396853.1 
     10              + y*(-2972611.439 
     11              + y*(15704.48260 
     12              + y*(-30.16036606)))))); 
     13    const double ans2 = 144725228442.0 
     14              + y*(2300535178.0 
     15              + y*(18583304.74 
     16              + y*(99447.43394 
     17              + y*(376.9991397 
    1818              + y)))); 
    1919    return ans1/ans2; 
    2020  } else { 
    21     const real y = REAL(64.0)/(ax*ax); 
    22     const real xx = ax - REAL(2.356194491); 
    23     const real ans1 = REAL(1.0) 
    24               + y*(REAL(0.183105e-2) 
    25               + y*(REAL(-0.3516396496e-4) 
    26               + y*(REAL(0.2457520174e-5) 
    27               + y*REAL(-0.240337019e-6)))); 
    28     const real ans2 = REAL(0.04687499995) 
    29               + y*(REAL(-0.2002690873e-3) 
    30               + y*(REAL(0.8449199096e-5) 
    31               + y*(REAL(-0.88228987e-6) 
    32               + y*REAL(0.105787412e-6)))); 
    33     real sn,cn; 
     21    const double y = 64.0/(ax*ax); 
     22    const double xx = ax - 2.356194491; 
     23    const double ans1 = 1.0 
     24              + y*(0.183105e-2 
     25              + y*(-0.3516396496e-4 
     26              + y*(0.2457520174e-5 
     27              + y*-0.240337019e-6))); 
     28    const double ans2 = 0.04687499995 
     29              + y*(-0.2002690873e-3 
     30              + y*(0.8449199096e-5 
     31              + y*(-0.88228987e-6 
     32              + y*0.105787412e-6))); 
     33    double sn,cn; 
    3434    SINCOS(xx, sn, cn); 
    35     const real ans = sqrt(REAL(0.636619772)/ax) * (cn*ans1 - (REAL(8.0)/ax)*sn*ans2); 
    36     return (x < REAL(0.0)) ? -ans : ans; 
     35    const double ans = sqrt(0.636619772/ax) * (cn*ans1 - (8.0/ax)*sn*ans2); 
     36    return (x < 0.0) ? -ans : ans; 
    3737  } 
    3838} 

  • sasmodels/models/lib/gauss150.c

    r14de349 r994d77f  
    99 
    1010// Gaussians 
    11 constant real Gauss150Z[150]={ 
    12         REAL(-0.9998723404457334), 
    13         REAL(-0.9993274305065947), 
    14         REAL(-0.9983473449340834), 
    15         REAL(-0.9969322929775997), 
    16         REAL(-0.9950828645255290), 
    17         REAL(-0.9927998590434373), 
    18         REAL(-0.9900842691660192), 
    19         REAL(-0.9869372772712794), 
    20         REAL(-0.9833602541697529), 
    21         REAL(-0.9793547582425894), 
    22         REAL(-0.9749225346595943), 
    23         REAL(-0.9700655145738374), 
    24         REAL(-0.9647858142586956), 
    25         REAL(-0.9590857341746905), 
    26         REAL(-0.9529677579610971), 
    27         REAL(-0.9464345513503147), 
    28         REAL(-0.9394889610042837), 
    29         REAL(-0.9321340132728527), 
    30         REAL(-0.9243729128743136), 
    31         REAL(-0.9162090414984952), 
    32         REAL(-0.9076459563329236), 
    33         REAL(-0.8986873885126239), 
    34         REAL(-0.8893372414942055), 
    35         REAL(-0.8795995893549102), 
    36         REAL(-0.8694786750173527), 
    37         REAL(-0.8589789084007133), 
    38         REAL(-0.8481048644991847), 
    39         REAL(-0.8368612813885015), 
    40         REAL(-0.8252530581614230), 
    41         REAL(-0.8132852527930605), 
    42         REAL(-0.8009630799369827), 
    43         REAL(-0.7882919086530552), 
    44         REAL(-0.7752772600680049), 
    45         REAL(-0.7619248049697269), 
    46         REAL(-0.7482403613363824), 
    47         REAL(-0.7342298918013638), 
    48         REAL(-0.7198995010552305), 
    49         REAL(-0.7052554331857488), 
    50         REAL(-0.6903040689571928), 
    51         REAL(-0.6750519230300931), 
    52         REAL(-0.6595056411226444), 
    53         REAL(-0.6436719971150083), 
    54         REAL(-0.6275578900977726), 
    55         REAL(-0.6111703413658551), 
    56         REAL(-0.5945164913591590), 
    57         REAL(-0.5776035965513142), 
    58         REAL(-0.5604390262878617), 
    59         REAL(-0.5430302595752546), 
    60         REAL(-0.5253848818220803), 
    61         REAL(-0.5075105815339176), 
    62         REAL(-0.4894151469632753), 
    63         REAL(-0.4711064627160663), 
    64         REAL(-0.4525925063160997), 
    65         REAL(-0.4338813447290861), 
    66         REAL(-0.4149811308476706), 
    67         REAL(-0.3959000999390257), 
    68         REAL(-0.3766465660565522), 
    69         REAL(-0.3572289184172501), 
    70         REAL(-0.3376556177463400), 
    71         REAL(-0.3179351925907259), 
    72         REAL(-0.2980762356029071), 
    73         REAL(-0.2780873997969574), 
    74         REAL(-0.2579773947782034), 
    75         REAL(-0.2377549829482451), 
    76         REAL(-0.2174289756869712), 
    77         REAL(-0.1970082295132342), 
    78         REAL(-0.1765016422258567), 
    79         REAL(-0.1559181490266516), 
    80         REAL(-0.1352667186271445), 
    81         REAL(-0.1145563493406956), 
    82         REAL(-0.0937960651617229), 
    83         REAL(-0.0729949118337358), 
    84         REAL(-0.0521619529078925), 
    85         REAL(-0.0313062657937972), 
    86         REAL(-0.0104369378042598), 
    87         REAL(0.0104369378042598), 
    88         REAL(0.0313062657937972), 
    89         REAL(0.0521619529078925), 
    90         REAL(0.0729949118337358), 
    91         REAL(0.0937960651617229), 
    92         REAL(0.1145563493406956), 
    93         REAL(0.1352667186271445), 
    94         REAL(0.1559181490266516), 
    95         REAL(0.1765016422258567), 
    96         REAL(0.1970082295132342), 
    97         REAL(0.2174289756869712), 
    98         REAL(0.2377549829482451), 
    99         REAL(0.2579773947782034), 
    100         REAL(0.2780873997969574), 
    101         REAL(0.2980762356029071), 
    102         REAL(0.3179351925907259), 
    103         REAL(0.3376556177463400), 
    104         REAL(0.3572289184172501), 
    105         REAL(0.3766465660565522), 
    106         REAL(0.3959000999390257), 
    107         REAL(0.4149811308476706), 
    108         REAL(0.4338813447290861), 
    109         REAL(0.4525925063160997), 
    110         REAL(0.4711064627160663), 
    111         REAL(0.4894151469632753), 
    112         REAL(0.5075105815339176), 
    113         REAL(0.5253848818220803), 
    114         REAL(0.5430302595752546), 
    115         REAL(0.5604390262878617), 
    116         REAL(0.5776035965513142), 
    117         REAL(0.5945164913591590), 
    118         REAL(0.6111703413658551), 
    119         REAL(0.6275578900977726), 
    120         REAL(0.6436719971150083), 
    121         REAL(0.6595056411226444), 
    122         REAL(0.6750519230300931), 
    123         REAL(0.6903040689571928), 
    124         REAL(0.7052554331857488), 
    125         REAL(0.7198995010552305), 
    126         REAL(0.7342298918013638), 
    127         REAL(0.7482403613363824), 
    128         REAL(0.7619248049697269), 
    129         REAL(0.7752772600680049), 
    130         REAL(0.7882919086530552), 
    131         REAL(0.8009630799369827), 
    132         REAL(0.8132852527930605), 
    133         REAL(0.8252530581614230), 
    134         REAL(0.8368612813885015), 
    135         REAL(0.8481048644991847), 
    136         REAL(0.8589789084007133), 
    137         REAL(0.8694786750173527), 
    138         REAL(0.8795995893549102), 
    139         REAL(0.8893372414942055), 
    140         REAL(0.8986873885126239), 
    141         REAL(0.9076459563329236), 
    142         REAL(0.9162090414984952), 
    143         REAL(0.9243729128743136), 
    144         REAL(0.9321340132728527), 
    145         REAL(0.9394889610042837), 
    146         REAL(0.9464345513503147), 
    147         REAL(0.9529677579610971), 
    148         REAL(0.9590857341746905), 
    149         REAL(0.9647858142586956), 
    150         REAL(0.9700655145738374), 
    151         REAL(0.9749225346595943), 
    152         REAL(0.9793547582425894), 
    153         REAL(0.9833602541697529), 
    154         REAL(0.9869372772712794), 
    155         REAL(0.9900842691660192), 
    156         REAL(0.9927998590434373), 
    157         REAL(0.9950828645255290), 
    158         REAL(0.9969322929775997), 
    159         REAL(0.9983473449340834), 
    160         REAL(0.9993274305065947), 
    161         REAL(0.9998723404457334) 
     11constant double Gauss150Z[150]={ 
     12        -0.9998723404457334, 
     13        -0.9993274305065947, 
     14        -0.9983473449340834, 
     15        -0.9969322929775997, 
     16        -0.9950828645255290, 
     17        -0.9927998590434373, 
     18        -0.9900842691660192, 
     19        -0.9869372772712794, 
     20        -0.9833602541697529, 
     21        -0.9793547582425894, 
     22        -0.9749225346595943, 
     23        -0.9700655145738374, 
     24        -0.9647858142586956, 
     25        -0.9590857341746905, 
     26        -0.9529677579610971, 
     27        -0.9464345513503147, 
     28        -0.9394889610042837, 
     29        -0.9321340132728527, 
     30        -0.9243729128743136, 
     31        -0.9162090414984952, 
     32        -0.9076459563329236, 
     33        -0.8986873885126239, 
     34        -0.8893372414942055, 
     35        -0.8795995893549102, 
     36        -0.8694786750173527, 
     37        -0.8589789084007133, 
     38        -0.8481048644991847, 
     39        -0.8368612813885015, 
     40        -0.8252530581614230, 
     41        -0.8132852527930605, 
     42        -0.8009630799369827, 
     43        -0.7882919086530552, 
     44        -0.7752772600680049, 
     45        -0.7619248049697269, 
     46        -0.7482403613363824, 
     47        -0.7342298918013638, 
     48        -0.7198995010552305, 
     49        -0.7052554331857488, 
     50        -0.6903040689571928, 
     51        -0.6750519230300931, 
     52        -0.6595056411226444, 
     53        -0.6436719971150083, 
     54        -0.6275578900977726, 
     55        -0.6111703413658551, 
     56        -0.5945164913591590, 
     57        -0.5776035965513142, 
     58        -0.5604390262878617, 
     59        -0.5430302595752546, 
     60        -0.5253848818220803, 
     61        -0.5075105815339176, 
     62        -0.4894151469632753, 
     63        -0.4711064627160663, 
     64        -0.4525925063160997, 
     65        -0.4338813447290861, 
     66        -0.4149811308476706, 
     67        -0.3959000999390257, 
     68        -0.3766465660565522, 
     69        -0.3572289184172501, 
     70        -0.3376556177463400, 
     71        -0.3179351925907259, 
     72        -0.2980762356029071, 
     73        -0.2780873997969574, 
     74        -0.2579773947782034, 
     75        -0.2377549829482451, 
     76        -0.2174289756869712, 
     77        -0.1970082295132342, 
     78        -0.1765016422258567, 
     79        -0.1559181490266516, 
     80        -0.1352667186271445, 
     81        -0.1145563493406956, 
     82        -0.0937960651617229, 
     83        -0.0729949118337358, 
     84        -0.0521619529078925, 
     85        -0.0313062657937972, 
     86        -0.0104369378042598, 
     87        0.0104369378042598, 
     88        0.0313062657937972, 
     89        0.0521619529078925, 
     90        0.0729949118337358, 
     91        0.0937960651617229, 
     92        0.1145563493406956, 
     93        0.1352667186271445, 
     94        0.1559181490266516, 
     95        0.1765016422258567, 
     96        0.1970082295132342, 
     97        0.2174289756869712, 
     98        0.2377549829482451, 
     99        0.2579773947782034, 
     100        0.2780873997969574, 
     101        0.2980762356029071, 
     102        0.3179351925907259, 
     103        0.3376556177463400, 
     104        0.3572289184172501, 
     105        0.3766465660565522, 
     106        0.3959000999390257, 
     107        0.4149811308476706, 
     108        0.4338813447290861, 
     109        0.4525925063160997, 
     110        0.4711064627160663, 
     111        0.4894151469632753, 
     112        0.5075105815339176, 
     113        0.5253848818220803, 
     114        0.5430302595752546, 
     115        0.5604390262878617, 
     116        0.5776035965513142, 
     117        0.5945164913591590, 
     118        0.6111703413658551, 
     119        0.6275578900977726, 
     120        0.6436719971150083, 
     121        0.6595056411226444, 
     122        0.6750519230300931, 
     123        0.6903040689571928, 
     124        0.7052554331857488, 
     125        0.7198995010552305, 
     126        0.7342298918013638, 
     127        0.7482403613363824, 
     128        0.7619248049697269, 
     129        0.7752772600680049, 
     130        0.7882919086530552, 
     131        0.8009630799369827, 
     132        0.8132852527930605, 
     133        0.8252530581614230, 
     134        0.8368612813885015, 
     135        0.8481048644991847, 
     136        0.8589789084007133, 
     137        0.8694786750173527, 
     138        0.8795995893549102, 
     139        0.8893372414942055, 
     140        0.8986873885126239, 
     141        0.9076459563329236, 
     142        0.9162090414984952, 
     143        0.9243729128743136, 
     144        0.9321340132728527, 
     145        0.9394889610042837, 
     146        0.9464345513503147, 
     147        0.9529677579610971, 
     148        0.9590857341746905, 
     149        0.9647858142586956, 
     150        0.9700655145738374, 
     151        0.9749225346595943, 
     152        0.9793547582425894, 
     153        0.9833602541697529, 
     154        0.9869372772712794, 
     155        0.9900842691660192, 
     156        0.9927998590434373, 
     157        0.9950828645255290, 
     158        0.9969322929775997, 
     159        0.9983473449340834, 
     160        0.9993274305065947, 
     161        0.9998723404457334 
    162162}; 
    163163 
    164 constant real Gauss150Wt[150]={ 
    165         REAL(0.0003276086705538), 
    166         REAL(0.0007624720924706), 
    167         REAL(0.0011976474864367), 
    168         REAL(0.0016323569986067), 
    169         REAL(0.0020663664924131), 
    170         REAL(0.0024994789888943), 
    171         REAL(0.0029315036836558), 
    172         REAL(0.0033622516236779), 
    173         REAL(0.0037915348363451), 
    174         REAL(0.0042191661429919), 
    175         REAL(0.0046449591497966), 
    176         REAL(0.0050687282939456), 
    177         REAL(0.0054902889094487), 
    178         REAL(0.0059094573005900), 
    179         REAL(0.0063260508184704), 
    180         REAL(0.0067398879387430), 
    181         REAL(0.0071507883396855), 
    182         REAL(0.0075585729801782), 
    183         REAL(0.0079630641773633), 
    184         REAL(0.0083640856838475), 
    185         REAL(0.0087614627643580), 
    186         REAL(0.0091550222717888), 
    187         REAL(0.0095445927225849), 
    188         REAL(0.0099300043714212), 
    189         REAL(0.0103110892851360), 
    190         REAL(0.0106876814158841), 
    191         REAL(0.0110596166734735), 
    192         REAL(0.0114267329968529), 
    193         REAL(0.0117888704247183), 
    194         REAL(0.0121458711652067), 
    195         REAL(0.0124975796646449), 
    196         REAL(0.0128438426753249), 
    197         REAL(0.0131845093222756), 
    198         REAL(0.0135194311690004), 
    199         REAL(0.0138484622795371), 
    200         REAL(0.0141714592928592), 
    201         REAL(0.0144882814685445), 
    202         REAL(0.0147987907597169), 
    203         REAL(0.0151028518701744), 
    204         REAL(0.0154003323133401), 
    205         REAL(0.0156911024699895), 
    206         REAL(0.0159750356447283), 
    207         REAL(0.0162520081211971), 
    208         REAL(0.0165218992159766), 
    209         REAL(0.0167845913311726), 
    210         REAL(0.0170399700056559), 
    211         REAL(0.0172879239649355), 
    212         REAL(0.0175283451696437), 
    213         REAL(0.0177611288626114), 
    214         REAL(0.0179861736145128), 
    215         REAL(0.0182033813680609), 
    216         REAL(0.0184126574807331), 
    217         REAL(0.0186139107660094), 
    218         REAL(0.0188070535331042), 
    219         REAL(0.0189920016251754), 
    220         REAL(0.0191686744559934), 
    221         REAL(0.0193369950450545), 
    222         REAL(0.0194968900511231), 
    223         REAL(0.0196482898041878), 
    224         REAL(0.0197911283358190), 
    225         REAL(0.0199253434079123), 
    226         REAL(0.0200508765398072), 
    227         REAL(0.0201676730337687), 
    228         REAL(0.0202756819988200), 
    229         REAL(0.0203748563729175), 
    230         REAL(0.0204651529434560), 
    231         REAL(0.0205465323660984), 
    232         REAL(0.0206189591819181), 
    233         REAL(0.0206824018328499), 
    234         REAL(0.0207368326754401), 
    235         REAL(0.0207822279928917), 
    236         REAL(0.0208185680053983), 
    237         REAL(0.0208458368787627), 
    238         REAL(0.0208640227312962), 
    239         REAL(0.0208731176389954), 
    240         REAL(0.0208731176389954), 
    241         REAL(0.0208640227312962), 
    242         REAL(0.0208458368787627), 
    243         REAL(0.0208185680053983), 
    244         REAL(0.0207822279928917), 
    245         REAL(0.0207368326754401), 
    246         REAL(0.0206824018328499), 
    247         REAL(0.0206189591819181), 
    248         REAL(0.0205465323660984), 
    249         REAL(0.0204651529434560), 
    250         REAL(0.0203748563729175), 
    251         REAL(0.0202756819988200), 
    252         REAL(0.0201676730337687), 
    253         REAL(0.0200508765398072), 
    254         REAL(0.0199253434079123), 
    255         REAL(0.0197911283358190), 
    256         REAL(0.0196482898041878), 
    257         REAL(0.0194968900511231), 
    258         REAL(0.0193369950450545), 
    259         REAL(0.0191686744559934), 
    260         REAL(0.0189920016251754), 
    261         REAL(0.0188070535331042), 
    262         REAL(0.0186139107660094), 
    263         REAL(0.0184126574807331), 
    264         REAL(0.0182033813680609), 
    265         REAL(0.0179861736145128), 
    266         REAL(0.0177611288626114), 
    267         REAL(0.0175283451696437), 
    268         REAL(0.0172879239649355), 
    269         REAL(0.0170399700056559), 
    270         REAL(0.0167845913311726), 
    271         REAL(0.0165218992159766), 
    272         REAL(0.0162520081211971), 
    273         REAL(0.0159750356447283), 
    274         REAL(0.0156911024699895), 
    275         REAL(0.0154003323133401), 
    276         REAL(0.0151028518701744), 
    277         REAL(0.0147987907597169), 
    278         REAL(0.0144882814685445), 
    279         REAL(0.0141714592928592), 
    280         REAL(0.0138484622795371), 
    281         REAL(0.0135194311690004), 
    282         REAL(0.0131845093222756), 
    283         REAL(0.0128438426753249), 
    284         REAL(0.0124975796646449), 
    285         REAL(0.0121458711652067), 
    286         REAL(0.0117888704247183), 
    287         REAL(0.0114267329968529), 
    288         REAL(0.0110596166734735), 
    289         REAL(0.0106876814158841), 
    290         REAL(0.0103110892851360), 
    291         REAL(0.0099300043714212), 
    292         REAL(0.0095445927225849), 
    293         REAL(0.0091550222717888), 
    294         REAL(0.0087614627643580), 
    295         REAL(0.0083640856838475), 
    296         REAL(0.0079630641773633), 
    297         REAL(0.0075585729801782), 
    298         REAL(0.0071507883396855), 
    299         REAL(0.0067398879387430), 
    300         REAL(0.0063260508184704), 
    301         REAL(0.0059094573005900), 
    302         REAL(0.0054902889094487), 
    303         REAL(0.0050687282939456), 
    304         REAL(0.0046449591497966), 
    305         REAL(0.0042191661429919), 
    306         REAL(0.0037915348363451), 
    307         REAL(0.0033622516236779), 
    308         REAL(0.0029315036836558), 
    309         REAL(0.0024994789888943), 
    310         REAL(0.0020663664924131), 
    311         REAL(0.0016323569986067), 
    312         REAL(0.0011976474864367), 
    313         REAL(0.0007624720924706), 
    314         REAL(0.0003276086705538) 
     164constant double Gauss150Wt[150]={ 
     165        0.0003276086705538, 
     166        0.0007624720924706, 
     167        0.0011976474864367, 
     168        0.0016323569986067, 
     169        0.0020663664924131, 
     170        0.0024994789888943, 
     171        0.0029315036836558, 
     172        0.0033622516236779, 
     173        0.0037915348363451, 
     174        0.0042191661429919, 
     175        0.0046449591497966, 
     176        0.0050687282939456, 
     177        0.0054902889094487, 
     178        0.0059094573005900, 
     179        0.0063260508184704, 
     180        0.0067398879387430, 
     181        0.0071507883396855, 
     182        0.0075585729801782, 
     183        0.0079630641773633, 
     184        0.0083640856838475, 
     185        0.0087614627643580, 
     186        0.0091550222717888, 
     187        0.0095445927225849, 
     188        0.0099300043714212, 
     189        0.0103110892851360, 
     190        0.0106876814158841, 
     191        0.0110596166734735, 
     192        0.0114267329968529, 
     193        0.0117888704247183, 
     194        0.0121458711652067, 
     195        0.0124975796646449, 
     196        0.0128438426753249, 
     197        0.0131845093222756, 
     198        0.0135194311690004, 
     199        0.0138484622795371, 
     200        0.0141714592928592, 
     201        0.0144882814685445, 
     202        0.0147987907597169, 
     203        0.0151028518701744, 
     204        0.0154003323133401, 
     205        0.0156911024699895, 
     206        0.0159750356447283, 
     207        0.0162520081211971, 
     208        0.0165218992159766, 
     209        0.0167845913311726, 
     210        0.0170399700056559, 
     211        0.0172879239649355, 
     212        0.0175283451696437, 
     213        0.0177611288626114, 
     214        0.0179861736145128, 
     215        0.0182033813680609, 
     216        0.0184126574807331, 
     217        0.0186139107660094, 
     218        0.0188070535331042, 
     219        0.0189920016251754, 
     220        0.0191686744559934, 
     221        0.0193369950450545, 
     222        0.0194968900511231, 
     223        0.0196482898041878, 
     224        0.0197911283358190, 
     225        0.0199253434079123, 
     226        0.0200508765398072, 
     227        0.0201676730337687, 
     228        0.0202756819988200, 
     229        0.0203748563729175, 
     230        0.0204651529434560, 
     231        0.0205465323660984, 
     232        0.0206189591819181, 
     233        0.0206824018328499, 
     234        0.0207368326754401, 
     235        0.0207822279928917, 
     236        0.0208185680053983, 
     237        0.0208458368787627, 
     238        0.0208640227312962, 
     239        0.0208731176389954, 
     240        0.0208731176389954, 
     241        0.0208640227312962, 
     242        0.0208458368787627, 
     243        0.0208185680053983, 
     244        0.0207822279928917, 
     245        0.0207368326754401, 
     246        0.0206824018328499, 
     247        0.0206189591819181, 
     248        0.0205465323660984, 
     249        0.0204651529434560, 
     250        0.0203748563729175, 
     251        0.0202756819988200, 
     252        0.0201676730337687, 
     253        0.0200508765398072, 
     254        0.0199253434079123, 
     255        0.0197911283358190, 
     256        0.0196482898041878, 
     257        0.0194968900511231, 
     258        0.0193369950450545, 
     259        0.0191686744559934, 
     260        0.0189920016251754, 
     261        0.0188070535331042, 
     262        0.0186139107660094, 
     263        0.0184126574807331, 
     264        0.0182033813680609, 
     265        0.0179861736145128, 
     266        0.0177611288626114, 
     267        0.0175283451696437, 
     268        0.0172879239649355, 
     269        0.0170399700056559, 
     270        0.0167845913311726, 
     271        0.0165218992159766, 
     272        0.0162520081211971, 
     273        0.0159750356447283, 
     274        0.0156911024699895, 
     275        0.0154003323133401, 
     276        0.0151028518701744, 
     277        0.0147987907597169, 
     278        0.0144882814685445, 
     279        0.0141714592928592, 
     280        0.0138484622795371, 
     281        0.0135194311690004, 
     282        0.0131845093222756, 
     283        0.0128438426753249, 
     284        0.0124975796646449, 
     285        0.0121458711652067, 
     286        0.0117888704247183, 
     287        0.0114267329968529, 
     288        0.0110596166734735, 
     289        0.0106876814158841, 
     290        0.0103110892851360, 
     291        0.0099300043714212, 
     292        0.0095445927225849, 
     293        0.0091550222717888, 
     294        0.0087614627643580, 
     295        0.0083640856838475, 
     296        0.0079630641773633, 
     297        0.0075585729801782, 
     298        0.0071507883396855, 
     299        0.0067398879387430, 
     300        0.0063260508184704, 
     301        0.0059094573005900, 
     302        0.0054902889094487, 
     303        0.0050687282939456, 
     304        0.0046449591497966, 
     305        0.0042191661429919, 
     306        0.0037915348363451, 
     307        0.0033622516236779, 
     308        0.0029315036836558, 
     309        0.0024994789888943, 
     310        0.0020663664924131, 
     311        0.0016323569986067, 
     312        0.0011976474864367, 
     313        0.0007624720924706, 
     314        0.0003276086705538 
    315315}; 

  • sasmodels/models/lib/gauss20.c

    r14de349 r994d77f  
    99 
    1010// Gaussians 
    11 constant real Gauss20Wt[20]={ 
    12         REAL(.0176140071391521), 
    13         REAL(.0406014298003869), 
    14         REAL(.0626720483341091), 
    15         REAL(.0832767415767047), 
    16         REAL(.10193011981724), 
    17         REAL(.118194531961518), 
    18         REAL(.131688638449177), 
    19         REAL(.142096109318382), 
    20         REAL(.149172986472604), 
    21         REAL(.152753387130726), 
    22         REAL(.152753387130726), 
    23         REAL(.149172986472604), 
    24         REAL(.142096109318382), 
    25         REAL(.131688638449177), 
    26         REAL(.118194531961518), 
    27         REAL(.10193011981724), 
    28         REAL(.0832767415767047), 
    29         REAL(.0626720483341091), 
    30         REAL(.0406014298003869), 
    31         REAL(.0176140071391521) 
     11constant double Gauss20Wt[20]={ 
     12        .0176140071391521, 
     13        .0406014298003869, 
     14        .0626720483341091, 
     15        .0832767415767047, 
     16        .10193011981724, 
     17        .118194531961518, 
     18        .131688638449177, 
     19        .142096109318382, 
     20        .149172986472604, 
     21        .152753387130726, 
     22        .152753387130726, 
     23        .149172986472604, 
     24        .142096109318382, 
     25        .131688638449177, 
     26        .118194531961518, 
     27        .10193011981724, 
     28        .0832767415767047, 
     29        .0626720483341091, 
     30        .0406014298003869, 
     31        .0176140071391521 
    3232}; 
    3333 
    34 constant real Gauss20Z[20]={ 
    35         REAL(-.993128599185095), 
    36         REAL(-.963971927277914), 
    37         REAL(-.912234428251326), 
    38         REAL(-.839116971822219), 
    39         REAL(-.746331906460151), 
    40         REAL(-.636053680726515), 
    41         REAL(-.510867001950827), 
    42         REAL(-.37370608871542), 
    43         REAL(-.227785851141645), 
    44         REAL(-.076526521133497), 
    45         REAL(.0765265211334973), 
    46         REAL(.227785851141645), 
    47         REAL(.37370608871542), 
    48         REAL(.510867001950827), 
    49         REAL(.636053680726515), 
    50         REAL(.746331906460151), 
    51         REAL(.839116971822219), 
    52         REAL(.912234428251326), 
    53         REAL(.963971927277914), 
    54         REAL(.993128599185095) 
     34constant double Gauss20Z[20]={ 
     35        -.993128599185095, 
     36        -.963971927277914, 
     37        -.912234428251326, 
     38        -.839116971822219, 
     39        -.746331906460151, 
     40        -.636053680726515, 
     41        -.510867001950827, 
     42        -.37370608871542, 
     43        -.227785851141645, 
     44        -.076526521133497, 
     45        .0765265211334973, 
     46        .227785851141645, 
     47        .37370608871542, 
     48        .510867001950827, 
     49        .636053680726515, 
     50        .746331906460151, 
     51        .839116971822219, 
     52        .912234428251326, 
     53        .963971927277914, 
     54        .993128599185095 
    5555}; 

  • sasmodels/models/lib/gauss76.c

    r14de349 r994d77f  
    99 
    1010// Gaussians 
    11 constant real Gauss76Wt[76]={ 
    12         REAL(.00126779163408536),               //0 
    13         REAL(.00294910295364247), 
    14         REAL(.00462793522803742), 
    15         REAL(.00629918049732845), 
    16         REAL(.00795984747723973), 
    17         REAL(.00960710541471375), 
    18         REAL(.0112381685696677), 
    19         REAL(.0128502838475101), 
    20         REAL(.0144407317482767), 
    21         REAL(.0160068299122486), 
    22         REAL(.0175459372914742),                //10 
    23         REAL(.0190554584671906), 
    24         REAL(.020532847967908), 
    25         REAL(.0219756145344162), 
    26         REAL(.0233813253070112), 
    27         REAL(.0247476099206597), 
    28         REAL(.026072164497986), 
    29         REAL(.0273527555318275), 
    30         REAL(.028587223650054), 
    31         REAL(.029773487255905), 
    32         REAL(.0309095460374916),                //20 
    33         REAL(.0319934843404216), 
    34         REAL(.0330234743977917), 
    35         REAL(.0339977794120564), 
    36         REAL(.0349147564835508), 
    37         REAL(.0357728593807139), 
    38         REAL(.0365706411473296), 
    39         REAL(.0373067565423816), 
    40         REAL(.0379799643084053), 
    41         REAL(.0385891292645067), 
    42         REAL(.0391332242205184),                //30 
    43         REAL(.0396113317090621), 
    44         REAL(.0400226455325968), 
    45         REAL(.040366472122844), 
    46         REAL(.0406422317102947), 
    47         REAL(.0408494593018285), 
    48         REAL(.040987805464794), 
    49         REAL(.0410570369162294), 
    50         REAL(.0410570369162294), 
    51         REAL(.040987805464794), 
    52         REAL(.0408494593018285),                //40 
    53         REAL(.0406422317102947), 
    54         REAL(.040366472122844), 
    55         REAL(.0400226455325968), 
    56         REAL(.0396113317090621), 
    57         REAL(.0391332242205184), 
    58         REAL(.0385891292645067), 
    59         REAL(.0379799643084053), 
    60         REAL(.0373067565423816), 
    61         REAL(.0365706411473296), 
    62         REAL(.0357728593807139),                //50 
    63         REAL(.0349147564835508), 
    64         REAL(.0339977794120564), 
    65         REAL(.0330234743977917), 
    66         REAL(.0319934843404216), 
    67         REAL(.0309095460374916), 
    68         REAL(.029773487255905), 
    69         REAL(.028587223650054), 
    70         REAL(.0273527555318275), 
    71         REAL(.026072164497986), 
    72         REAL(.0247476099206597),                //60 
    73         REAL(.0233813253070112), 
    74         REAL(.0219756145344162), 
    75         REAL(.020532847967908), 
    76         REAL(.0190554584671906), 
    77         REAL(.0175459372914742), 
    78         REAL(.0160068299122486), 
    79         REAL(.0144407317482767), 
    80         REAL(.0128502838475101), 
    81         REAL(.0112381685696677), 
    82         REAL(.00960710541471375),               //70 
    83         REAL(.00795984747723973), 
    84         REAL(.00629918049732845), 
    85         REAL(.00462793522803742), 
    86         REAL(.00294910295364247), 
    87         REAL(.00126779163408536)                //75 (indexed from 0) 
     11constant double Gauss76Wt[76]={ 
     12        .00126779163408536,             //0 
     13        .00294910295364247, 
     14        .00462793522803742, 
     15        .00629918049732845, 
     16        .00795984747723973, 
     17        .00960710541471375, 
     18        .0112381685696677, 
     19        .0128502838475101, 
     20        .0144407317482767, 
     21        .0160068299122486, 
     22        .0175459372914742,              //10 
     23        .0190554584671906, 
     24        .020532847967908, 
     25        .0219756145344162, 
     26        .0233813253070112, 
     27        .0247476099206597, 
     28        .026072164497986, 
     29        .0273527555318275, 
     30        .028587223650054, 
     31        .029773487255905, 
     32        .0309095460374916,              //20 
     33        .0319934843404216, 
     34        .0330234743977917, 
     35        .0339977794120564, 
     36        .0349147564835508, 
     37        .0357728593807139, 
     38        .0365706411473296, 
     39        .0373067565423816, 
     40        .0379799643084053, 
     41        .0385891292645067, 
     42        .0391332242205184,              //30 
     43        .0396113317090621, 
     44        .0400226455325968, 
     45        .040366472122844, 
     46        .0406422317102947, 
     47        .0408494593018285, 
     48        .040987805464794, 
     49        .0410570369162294, 
     50        .0410570369162294, 
     51        .040987805464794, 
     52        .0408494593018285,              //40 
     53        .0406422317102947, 
     54        .040366472122844, 
     55        .0400226455325968, 
     56        .0396113317090621, 
     57        .0391332242205184, 
     58        .0385891292645067, 
     59        .0379799643084053, 
     60        .0373067565423816, 
     61        .0365706411473296, 
     62        .0357728593807139,              //50 
     63        .0349147564835508, 
     64        .0339977794120564, 
     65        .0330234743977917, 
     66        .0319934843404216, 
     67        .0309095460374916, 
     68        .029773487255905, 
     69        .028587223650054, 
     70        .0273527555318275, 
     71        .026072164497986, 
     72        .0247476099206597,              //60 
     73        .0233813253070112, 
     74        .0219756145344162, 
     75        .020532847967908, 
     76        .0190554584671906, 
     77        .0175459372914742, 
     78        .0160068299122486, 
     79        .0144407317482767, 
     80        .0128502838475101, 
     81        .0112381685696677, 
     82        .00960710541471375,             //70 
     83        .00795984747723973, 
     84        .00629918049732845, 
     85        .00462793522803742, 
     86        .00294910295364247, 
     87        .00126779163408536              //75 (indexed from 0) 
    8888}; 
    8989 
    90 constant real Gauss76Z[76]={ 
    91         REAL(-.999505948362153),                //0 
    92         REAL(-.997397786355355), 
    93         REAL(-.993608772723527), 
    94         REAL(-.988144453359837), 
    95         REAL(-.981013938975656), 
    96         REAL(-.972229228520377), 
    97         REAL(-.961805126758768), 
    98         REAL(-.949759207710896), 
    99         REAL(-.936111781934811), 
    100         REAL(-.92088586125215), 
    101         REAL(-.904107119545567),                //10 
    102         REAL(-.885803849292083), 
    103         REAL(-.866006913771982), 
    104         REAL(-.844749694983342), 
    105         REAL(-.822068037328975), 
    106         REAL(-.7980001871612), 
    107         REAL(-.77258672828181), 
    108         REAL(-.74587051350361), 
    109         REAL(-.717896592387704), 
    110         REAL(-.688712135277641), 
    111         REAL(-.658366353758143),                //20 
    112         REAL(-.626910417672267), 
    113         REAL(-.594397368836793), 
    114         REAL(-.560882031601237), 
    115         REAL(-.526420920401243), 
    116         REAL(-.491072144462194), 
    117         REAL(-.454895309813726), 
    118         REAL(-.417951418780327), 
    119         REAL(-.380302767117504), 
    120         REAL(-.342012838966962), 
    121         REAL(-.303146199807908),                //30 
    122         REAL(-.263768387584994), 
    123         REAL(-.223945802196474), 
    124         REAL(-.183745593528914), 
    125         REAL(-.143235548227268), 
    126         REAL(-.102483975391227), 
    127         REAL(-.0615595913906112), 
    128         REAL(-.0205314039939986), 
    129         REAL(.0205314039939986), 
    130         REAL(.0615595913906112), 
    131         REAL(.102483975391227),                 //40 
    132         REAL(.143235548227268), 
    133         REAL(.183745593528914), 
    134         REAL(.223945802196474), 
    135         REAL(.263768387584994), 
    136         REAL(.303146199807908), 
    137         REAL(.342012838966962), 
    138         REAL(.380302767117504), 
    139         REAL(.417951418780327), 
    140         REAL(.454895309813726), 
    141         REAL(.491072144462194),         //50 
    142         REAL(.526420920401243), 
    143         REAL(.560882031601237), 
    144         REAL(.594397368836793), 
    145         REAL(.626910417672267), 
    146         REAL(.658366353758143), 
    147         REAL(.688712135277641), 
    148         REAL(.717896592387704), 
    149         REAL(.74587051350361), 
    150         REAL(.77258672828181), 
    151         REAL(.7980001871612),   //60 
    152         REAL(.822068037328975), 
    153         REAL(.844749694983342), 
    154         REAL(.866006913771982), 
    155         REAL(.885803849292083), 
    156         REAL(.904107119545567), 
    157         REAL(.92088586125215), 
    158         REAL(.936111781934811), 
    159         REAL(.949759207710896), 
    160         REAL(.961805126758768), 
    161         REAL(.972229228520377),         //70 
    162         REAL(.981013938975656), 
    163         REAL(.988144453359837), 
    164         REAL(.993608772723527), 
    165         REAL(.997397786355355), 
    166         REAL(.999505948362153)          //75 
     90constant double Gauss76Z[76]={ 
     91        -.999505948362153,              //0 
     92        -.997397786355355, 
     93        -.993608772723527, 
     94        -.988144453359837, 
     95        -.981013938975656, 
     96        -.972229228520377, 
     97        -.961805126758768, 
     98        -.949759207710896, 
     99        -.936111781934811, 
     100        -.92088586125215, 
     101        -.904107119545567,              //10 
     102        -.885803849292083, 
     103        -.866006913771982, 
     104        -.844749694983342, 
     105        -.822068037328975, 
     106        -.7980001871612, 
     107        -.77258672828181, 
     108        -.74587051350361, 
     109        -.717896592387704, 
     110        -.688712135277641, 
     111        -.658366353758143,              //20 
     112        -.626910417672267, 
     113        -.594397368836793, 
     114        -.560882031601237, 
     115        -.526420920401243, 
     116        -.491072144462194, 
     117        -.454895309813726, 
     118        -.417951418780327, 
     119        -.380302767117504, 
     120        -.342012838966962, 
     121        -.303146199807908,              //30 
     122        -.263768387584994, 
     123        -.223945802196474, 
     124        -.183745593528914, 
     125        -.143235548227268, 
     126        -.102483975391227, 
     127        -.0615595913906112, 
     128        -.0205314039939986, 
     129        .0205314039939986, 
     130        .0615595913906112, 
     131        .102483975391227,                       //40 
     132        .143235548227268, 
     133        .183745593528914, 
     134        .223945802196474, 
     135        .263768387584994, 
     136        .303146199807908, 
     137        .342012838966962, 
     138        .380302767117504, 
     139        .417951418780327, 
     140        .454895309813726, 
     141        .491072144462194,               //50 
     142        .526420920401243, 
     143        .560882031601237, 
     144        .594397368836793, 
     145        .626910417672267, 
     146        .658366353758143, 
     147        .688712135277641, 
     148        .717896592387704, 
     149        .74587051350361, 
     150        .77258672828181, 
     151        .7980001871612, //60 
     152        .822068037328975, 
     153        .844749694983342, 
     154        .866006913771982, 
     155        .885803849292083, 
     156        .904107119545567, 
     157        .92088586125215, 
     158        .936111781934811, 
     159        .949759207710896, 
     160        .961805126758768, 
     161        .972229228520377,               //70 
     162        .981013938975656, 
     163        .988144453359837, 
     164        .993608772723527, 
     165        .997397786355355, 
     166        .999505948362153                //75 
    167167}; 

  • sasmodels/models/sphere.py

    r19dcb933 r994d77f  
    8686# This should perhaps be volume normalized? 
    8787form_volume = """ 
    88     return REAL(1.333333333333333)*M_PI*radius*radius*radius; 
     88    return 1.333333333333333*M_PI*radius*radius*radius; 
    8989    """ 
    9090 
    9191Iq = """ 
    92     const real qr = q*radius; 
    93     real sn, cn; 
     92    const double qr = q*radius; 
     93    double sn, cn; 
    9494    SINCOS(qr, sn, cn); 
    95     const real bes = qr==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(sn-qr*cn)/(qr*qr*qr); 
    96     const real fq = bes * (sld - solvent_sld) * form_volume(radius); 
    97     return REAL(1.0e-4)*fq*fq; 
     95    const double bes = qr==0.0 ? 1.0 : 3.0*(sn-qr*cn)/(qr*qr*qr); 
     96    const double fq = bes * (sld - solvent_sld) * form_volume(radius); 
     97    return 1.0e-4*fq*fq; 
    9898    """ 
    9999 
    … …  
    101101Iqxy = """ 
    102102    // never called since no orientation or magnetic parameters. 
    103     return REAL(-1.0); 
     103    //return -1.0; 
     104    return Iq(sqrt(qx*qx + qy*qy), sld, solvent_sld, radius); 
    104105    """ 
    105106 

  • sasmodels/models/triaxial_ellipsoid.c

    r5d4777d r994d77f  
    1 real form_volume(real req_minor, real req_major, real rpolar); 
    2 real Iq(real q, real sld, real solvent_sld, 
    3     real req_minor, real req_major, real rpolar); 
    4 real Iqxy(real qx, real qy, real sld, real solvent_sld, 
    5     real req_minor, real req_major, real rpolar, real theta, real phi, real psi); 
     1double form_volume(double req_minor, double req_major, double rpolar); 
     2double Iq(double q, double sld, double solvent_sld, 
     3    double req_minor, double req_major, double rpolar); 
     4double Iqxy(double qx, double qy, double sld, double solvent_sld, 
     5    double req_minor, double req_major, double rpolar, double theta, double phi, double psi); 
    66 
    7 real form_volume(real req_minor, real req_major, real rpolar) 
     7double form_volume(double req_minor, double req_major, double rpolar) 
    88{ 
    9     return REAL(1.333333333333333)*M_PI*req_minor*req_major*rpolar; 
     9    return 1.333333333333333*M_PI*req_minor*req_major*rpolar; 
    1010} 
    1111 
    12 real Iq(real q, 
    13     real sld, 
    14     real solvent_sld, 
    15     real req_minor, 
    16     real req_major, 
    17     real rpolar) 
     12double Iq(double q, 
     13    double sld, 
     14    double solvent_sld, 
     15    double req_minor, 
     16    double req_major, 
     17    double rpolar) 
    1818{ 
    19     // if (req_minor > req_major || req_major > rpolar) return REAL(-1.0);  // Exclude invalid region 
     19    // if (req_minor > req_major || req_major > rpolar) return -1.0;  // Exclude invalid region 
    2020 
    21     real sn, cn; 
    22     real st, ct; 
    23     //const real lower = REAL(0.0); 
    24     //const real upper = REAL(1.0); 
    25     real outer = REAL(0.0); 
     21    double sn, cn; 
     22    double st, ct; 
     23    //const double lower = 0.0; 
     24    //const double upper = 1.0; 
     25    double outer = 0.0; 
    2626    for (int i=0;i<76;i++) { 
    27         //const real cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 
    28         const real x = REAL(0.5)*(Gauss76Z[i] + REAL(1.0)); 
     27        //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 
     28        const double x = 0.5*(Gauss76Z[i] + 1.0); 
    2929        SINCOS(M_PI_2*x, sn, cn); 
    30         const real acosx2 = req_minor*req_minor*cn*cn; 
    31         const real bsinx2 = req_major*req_major*sn*sn; 
    32         const real c2 = rpolar*rpolar; 
     30        const double acosx2 = req_minor*req_minor*cn*cn; 
     31        const double bsinx2 = req_major*req_major*sn*sn; 
     32        const double c2 = rpolar*rpolar; 
    3333 
    34         real inner = REAL(0.0); 
     34        double inner = 0.0; 
    3535        for (int j=0;j<76;j++) { 
    36             const real y = REAL(0.5)*(Gauss76Z[j] + REAL(1.0)); 
    37             const real t = q*sqrt(acosx2 + bsinx2*(REAL(1.0)-y*y) + c2*y*y); 
     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); 
    3838            SINCOS(t, st, ct); 
    39             const real fq = ( t==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(st-t*ct)/(t*t*t) ); 
     39            const double fq = ( t==0.0 ? 1.0 : 3.0*(st-t*ct)/(t*t*t) ); 
    4040            inner += Gauss76Wt[j] * fq * fq ; 
    4141        } 
    42         outer += Gauss76Wt[i] * REAL(0.5) * inner; 
     42        outer += Gauss76Wt[i] * 0.5 * inner; 
    4343    } 
    44     //const real fq2 = (upper-lower)/2*outer; 
    45     const real fq2 = REAL(0.5)*outer; 
    46     const real s = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar); 
    47     return REAL(1.0e-4) * fq2 * s * s; 
     44    //const double fq2 = (upper-lower)/2*outer; 
     45    const double fq2 = 0.5*outer; 
     46    const double s = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar); 
     47    return 1.0e-4 * fq2 * s * s; 
    4848} 
    4949 
    50 real Iqxy(real qx, real qy, 
    51     real sld, 
    52     real solvent_sld, 
    53     real req_minor, 
    54     real req_major, 
    55     real rpolar, 
    56     real theta, 
    57     real phi, 
    58     real psi) 
     50double Iqxy(double qx, double qy, 
     51    double sld, 
     52    double solvent_sld, 
     53    double req_minor, 
     54    double req_major, 
     55    double rpolar, 
     56    double theta, 
     57    double phi, 
     58    double psi) 
    5959{ 
    60     // if (req_minor > req_major || req_major > rpolar) return REAL(-1.0);  // Exclude invalid region 
     60    // if (req_minor > req_major || req_major > rpolar) return -1.0;  // Exclude invalid region 
    6161 
    62     real stheta, ctheta; 
    63     real sphi, cphi; 
    64     real spsi, cpsi; 
    65     real st, ct; 
     62    double stheta, ctheta; 
     63    double sphi, cphi; 
     64    double spsi, cpsi; 
     65    double st, ct; 
    6666 
    67     const real q = sqrt(qx*qx + qy*qy); 
    68     const real qxhat = qx/q; 
    69     const real qyhat = qy/q; 
     67    const double q = sqrt(qx*qx + qy*qy); 
     68    const double qxhat = qx/q; 
     69    const double qyhat = qy/q; 
    7070    SINCOS(theta*M_PI_180, stheta, ctheta); 
    7171    SINCOS(phi*M_PI_180, sphi, cphi); 
    7272    SINCOS(psi*M_PI_180, spsi, cpsi); 
    73     const real calpha = ctheta*cphi*qxhat + stheta*qyhat; 
    74     const real cnu = (-cphi*spsi*stheta + sphi*cpsi)*qxhat + spsi*ctheta*qyhat; 
    75     const real cmu = (-stheta*cpsi*cphi - spsi*sphi)*qxhat + ctheta*cpsi*qyhat; 
    76     const real t = q*sqrt(req_minor*req_minor*cnu*cnu 
     73    const double calpha = ctheta*cphi*qxhat + stheta*qyhat; 
     74    const double cnu = (-cphi*spsi*stheta + sphi*cpsi)*qxhat + spsi*ctheta*qyhat; 
     75    const double cmu = (-stheta*cpsi*cphi - spsi*sphi)*qxhat + ctheta*cpsi*qyhat; 
     76    const double t = q*sqrt(req_minor*req_minor*cnu*cnu 
    7777                          + req_major*req_major*cmu*cmu 
    7878                          + rpolar*rpolar*calpha*calpha); 
    7979    SINCOS(t, st, ct); 
    80     const real fq = ( t==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(st-t*ct)/(t*t*t) ); 
    81     const real s = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar); 
     80    const double fq = ( t==0.0 ? 1.0 : 3.0*(st-t*ct)/(t*t*t) ); 
     81    const double s = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar); 
    8282 
    83     return REAL(1.0e-4) * fq * fq * s * s; 
     83    return 1.0e-4 * fq * fq * s * s; 
    8484} 
    8585 
Note: See TracChangeset for help on using the changeset viewer.