Changeset 73cbc5b in sasmodels


Ignore:
Timestamp:
Jun 7, 2017 11:36:23 AM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
4f0c993
Parents:
b34fc77 (diff), c1114bf (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into ticket-776-orientation

Files:
7 added
7 deleted
39 edited

Legend:

Unmodified
Added
Removed

  • doc/ref/gpu/gpu_computations.rst

    r3f5a566 r7e74ed5  
    3131from available OpenCL platforms. 
    3232 
     33OpenCL devices can be set from OpenCL options dialog in Fitting menu or as 
     34enviromental variables. 
     35 
     36**If you don't want to use OpenCL, you can select "No OpenCL" option from** 
     37**GUI dialog or set *SAS_OPENCL=None* in your environment settings** 
     38**This will only use normal programs.** 
     39 
    3340SasView prefers AMD or NVIDIA drivers for GPU, and prefers Intel or 
    3441Apple drivers for CPU. Both GPU and CPU are included on the assumption that CPU  
    … …  
    3946chose to run the model. 
    4047 
    41 **If you don't want to use OpenCL, you can set** *SAS_OPENCL=None* 
    42 **in your environment settings, and it will only use normal programs.** 
    43  
    44 If you want to use one of the other devices, you can run the following 
    45 from the python console in SasView:: 
     48If you want to use one of the other devices, you can select it from OpenCL 
     49options dialog (accessible from Fitting menu) or run the following from 
     50the python console in SasView:: 
    4651 
    4752    import pyopencl as cl 

  • sasmodels/generate.py

    rc4e3215 rbb4b509  
    492492 
    493493def test_tag_float(): 
    494  
    495     cases=""" 
     494    """check that floating point constants are properly identified and tagged with 'f'""" 
     495 
     496    cases = """ 
    496497ZP  : 0. 
    497498ZPF : 0.0,0.01,0.1 
    … …  
    519520""" 
    520521 
    521     output=""" 
     522    output = """ 
    522523ZP  : 0.f 
    523524ZPF : 0.0f,0.01f,0.1f 
    … …  
    611612    # type: (str, List[Parameter]) -> List[str] 
    612613    """ 
    613     Return a list of *prefix.parameter* from parameter items. 
     614    Return a list of *prefix+parameter* from parameter items. 
     615 
     616    *prefix* should be "v." if v is a struct. 
    614617    """ 
    615618    return [p.as_call_reference(prefix) for p in pars] 
    … …  
    733736        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt)) 
    734737 
    735     magpars = [k-2 for k,p in enumerate(partable.call_parameters) 
     738    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
    736739               if p.type == 'sld'] 
    737740 
    … …  
    742745    source.append("#define NUM_MAGNETIC %d" % partable.nmagnetic) 
    743746    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    744     for k,v in enumerate(magpars[:3]): 
     747    for k, v in enumerate(magpars[:3]): 
    745748        source.append("#define MAGNETIC_PAR%d %d"%(k+1, v)) 
    746749 
    … …  
    779782        "#undef CALL_IQ", 
    780783        "#undef KERNEL_NAME", 
    781          ] 
     784    ] 
    782785 
    783786    imagnetic = [ 
    … …  
    872875 
    873876# TODO: need a single source for rst_prolog; it is also in doc/rst_prolog 
    874 RST_PROLOG = """\ 
     877RST_PROLOG = r"""\ 
    875878.. |Ang| unicode:: U+212B 
    876879.. |Ang^-1| replace:: |Ang|\ :sup:`-1` 

  • sasmodels/kernel_header.c

    r9901384 r73cbc5b  
    1414#    define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0) 
    1515#  endif 
    16    // Intel CPU on Mac gives strange values for erf(); also on the tested 
     16   // Intel CPU on Mac gives strange values for erf(); on the verified 
    1717   // platforms (intel, nvidia, amd), the cephes erf() is significantly 
    1818   // faster than that available in the native OpenCL. 
    … …  
    5757         typedef int int32_t; 
    5858         #include <math.h> 
    59          // TODO: test isnan 
     59         // TODO: check isnan is correct 
    6060         inline double _isnan(double x) { return x != x; } // hope this doesn't optimize away! 
    6161         #undef isnan 

  • sasmodels/kernelcl.py

    rd2b939c rc1114bf  
    578578        # Free buffers 
    579579        for v in (details_b, values_b): 
    580             if v is not None: v.release() 
     580            if v is not None: 
     581                v.release() 
    581582 
    582583        pd_norm = self.result[self.q_input.nq] 
    583         scale = values[0]/(pd_norm if pd_norm!=0.0 else 1.0) 
     584        scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    584585        background = values[1] 
    585586        #print("scale",scale,values[0],self.result[self.q_input.nq],background) 

  • sasmodels/model_test.py

    reaa4458 r73cbc5b  
    8585        skip = [] 
    8686    for model_name in models: 
    87         if model_name in skip: continue 
     87        if model_name in skip: 
     88            continue 
    8889        model_info = load_model_info(model_name) 
    8990 
    … …  
    239240            multiple = [test for test in self.info.tests 
    240241                        if isinstance(test[2], list) 
    241                             and not all(result is None for result in test[2])] 
     242                        and not all(result is None for result in test[2])] 
    242243            tests_has_1D_multiple = any(isinstance(test[1][0], float) 
    243244                                        for test in multiple) 
    … …  
    262263            user_pars, x, y = test 
    263264            pars = expand_pars(self.info.parameters, user_pars) 
     265            invalid = invalid_pars(self.info.parameters, pars) 
     266            if invalid: 
     267                raise ValueError("Unknown parameters in test: " + ", ".join(invalid)) 
    264268 
    265269            if not isinstance(y, list): 
    … …  
    305309 
    306310    return ModelTestCase 
     311 
     312def invalid_pars(partable, pars): 
     313    # type: (ParameterTable, Dict[str, float]) 
     314    """ 
     315    Return a list of parameter names that are not part of the model. 
     316    """ 
     317    names = set(p.id for p in partable.call_parameters) 
     318    invalid = [] 
     319    for par in sorted(pars.keys()): 
     320        parts = par.split('_pd') 
     321        if len(parts) > 1 and parts[1] not in ("", "_n", "nsigma", "type"): 
     322            invalid.append(par) 
     323            continue 
     324        if parts[0] not in names: 
     325            invalid.append(par) 
     326    return invalid 
     327 
    307328 
    308329def is_near(target, actual, digits=5): 

  • sasmodels/modelinfo.py

    r9c44b7b rbb4b509  
    101101                limits = (float(low), float(high)) 
    102102            except Exception: 
    103                 raise ValueError("invalid limits for %s: %r"%(name,user_limits)) 
     103                raise ValueError("invalid limits for %s: %r"%(name, user_limits)) 
    104104            if low >= high: 
    105105                raise ValueError("require lower limit < upper limit") 
    … …  
    342342    def as_call_reference(self, prefix=""): 
    343343        # type: (str) -> str 
     344        """ 
     345        Return *prefix* + parameter name.  For struct references, use "v." 
     346        as the prefix. 
     347        """ 
    344348        # Note: if the parameter is a struct type, then we will need to use 
    345349        # &prefix+id.  For scalars and vectors we can just use prefix+id. 
    … …  
    420424        self.npars = sum(p.length for p in self.kernel_parameters) 
    421425        self.nmagnetic = sum(p.length for p in self.kernel_parameters 
    422                              if p.type=='sld') 
     426                             if p.type == 'sld') 
    423427        self.nvalues = 2 + self.npars 
    424428        if self.nmagnetic: 
    … …  
    457461        self.has_2d = any(p.type in ('orientation', 'magnetic') 
    458462                          for p in self.kernel_parameters) 
    459         self.magnetism_index = [k for k,p in enumerate(self.call_parameters) 
     463        self.magnetism_index = [k for k, p in enumerate(self.call_parameters) 
    460464                                if p.id.startswith('M0:')] 
    461465 
    … …  
    544548                              'magnetic', 'magnetic amplitude for '+p.description), 
    545549                    Parameter('mtheta:'+p.id, 'degrees', 0., [-90., 90.], 
    546                                'magnetic', 'magnetic latitude for '+p.description), 
     550                              'magnetic', 'magnetic latitude for '+p.description), 
    547551                    Parameter('mphi:'+p.id, 'degrees', 0., [-180., 180.], 
    548                                'magnetic', 'magnetic longitude for '+p.description), 
     552                              'magnetic', 'magnetic longitude for '+p.description), 
    549553                ]) 
    550554        #print("call parameters", full_list) 
    … …  
    683687 
    684688    if (model_info.Iq is None 
    685         and model_info.Iqxy is None 
    686         and model_info.Imagnetic is None 
    687         and model_info.form_volume is None): 
     689            and model_info.Iqxy is None 
     690            and model_info.Imagnetic is None 
     691            and model_info.form_volume is None): 
    688692        return 
    689693 
    … …  
    843847    #: it to False because they require double precision calculations. 
    844848    single = None           # type: bool 
     849    #: True if the model can be run as an opencl model.  If for some reason 
     850    #: the model cannot be run in opencl (e.g., because the model passes 
     851    #: functions by reference), then set this to false. 
     852    opencl = None           # type: bool 
    845853    #: True if the model is a structure factor used to model the interaction 
    846854    #: between form factor models.  This will default to False if it is not 

  • sasmodels/models/hollow_cylinder.py

    r9b79f29 rf102a96  
    8080    ["thickness",   "Ang",     10.0, [0, inf],    "volume",      "Cylinder wall thickness"], 
    8181    ["length",      "Ang",    400.0, [0, inf],    "volume",      "Cylinder total length"], 
    82     ["sld",         "1/Ang^2",  6.3, [-inf, inf], "sld",         "Cylinder sld"], 
    83     ["sld_solvent", "1/Ang^2",  1,   [-inf, inf], "sld",         "Solvent sld"], 
     82    ["sld",         "1e-6/Ang^2",  6.3, [-inf, inf], "sld",         "Cylinder sld"], 
     83    ["sld_solvent", "1e-6/Ang^2",  1,   [-inf, inf], "sld",         "Solvent sld"], 
    8484    ["theta",       "degrees", 90,   [-360, 360], "orientation", "Cylinder axis to beam angle"], 
    8585    ["phi",         "degrees",  0,   [-360, 360], "orientation", "Rotation about beam"], 

  • sasmodels/models/lib/sas_3j1x_x.c

    r473a9f1 reb2946f  
    4646double sas_3j1x_x(double q) 
    4747{ 
    48     if (q < SPH_J1C_CUTOFF) { 
     48    // 2017-05-18 PAK - support negative q 
     49    if (fabs(q) < SPH_J1C_CUTOFF) { 
    4950        const double q2 = q*q; 
    5051        return (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))));// + q2*(3./3991680.))))); 

  • sasmodels/models/lib/sas_J0.c

    rc8902ac reb2946f  
    236236        xx = x; 
    237237 
    238     if( x <= 2.0 ) { 
     238    // 2017-05-18 PAK - support negative x 
     239    if( xx <= 2.0 ) { 
    239240        z = xx * xx; 
    240         if( x < 1.0e-3 ) 
     241        if( xx < 1.0e-3 ) 
    241242            return( 1.0 - 0.25*z ); 
    242243 
    … …  
    245246    } 
    246247 
    247     q = 1.0/x; 
     248    q = 1.0/xx; 
    248249    w = sqrt(q); 
    249250 

  • sasmodels/models/lib/sas_J1.c

    r473a9f1 r5181ccc  
    109109{ 
    110110 
    111     double w, z, p, q, xn; 
     111    double w, z, p, q, abs_x, sign_x; 
    112112 
    113113    const double Z1 = 1.46819706421238932572E1; 
    114114    const double Z2 = 4.92184563216946036703E1; 
    115     const double THPIO4 =  2.35619449019234492885; 
    116     const double SQ2OPI = 0.79788456080286535588; 
    117  
    118     w = x; 
    119     if( x < 0 ) 
    120         w = -x; 
    121  
    122     if( w <= 5.0 ) { 
    123         z = x * x; 
     115 
     116    // 2017-05-18 PAK - mathematica and mpmath use J1(-x) = -J1(x) 
     117    if (x < 0) { 
     118        abs_x = -x; 
     119        sign_x = -1.0; 
     120    } else { 
     121        abs_x = x; 
     122        sign_x = 1.0; 
     123    } 
     124 
     125    if( abs_x <= 5.0 ) { 
     126        z = abs_x * abs_x; 
    124127        w = polevl( z, RPJ1, 3 ) / p1evl( z, RQJ1, 8 ); 
    125         w = w * x * (z - Z1) * (z - Z2); 
    126         return( w ); 
     128        w = w * abs_x * (z - Z1) * (z - Z2); 
     129        return( sign_x * w ); 
    127130    } 
    128131 
    129     w = 5.0/x; 
     132    w = 5.0/abs_x; 
    130133    z = w * w; 
    131  
    132134    p = polevl( z, PPJ1, 6)/polevl( z, PQJ1, 6 ); 
    133135    q = polevl( z, QPJ1, 7)/p1evl( z, QQJ1, 7 ); 
    134136 
    135     xn = x - THPIO4; 
    136  
    137     double sn, cn; 
    138     SINCOS(xn, sn, cn); 
    139     p = p * cn - w * q * sn; 
    140  
    141     return( p * SQ2OPI / sqrt(x) ); 
     137    // 2017-05-19 PAK improve accuracy using trig identies 
     138    // original: 
     139    //    const double THPIO4 =  2.35619449019234492885; 
     140    //    const double SQ2OPI = 0.79788456080286535588; 
     141    //    double sin_xn, cos_xn; 
     142    //    SINCOS(abs_x - THPIO4, sin_xn, cos_xn); 
     143    //    p = p * cos_xn - w * q * sin_xn; 
     144    //    return( sign_x * p * SQ2OPI / sqrt(abs_x) ); 
     145    // expanding p*cos(a - 3 pi/4) - wq sin(a - 3 pi/4) 
     146    //    [ p(sin(a) - cos(a)) + wq(sin(a) + cos(a)) / sqrt(2) 
     147    // note that sqrt(1/2) * sqrt(2/pi) = sqrt(1/pi) 
     148    const double SQRT1_PI = 0.56418958354775628; 
     149    double sin_x, cos_x; 
     150    SINCOS(abs_x, sin_x, cos_x); 
     151    p = p*(sin_x - cos_x) + w*q*(sin_x + cos_x); 
     152    return( sign_x * p * SQRT1_PI / sqrt(abs_x) ); 
    142153} 
    143154 
    … …  
    179190    }; 
    180191 
    181 float cephes_j1f(float x) 
     192float cephes_j1f(float xx) 
    182193{ 
    183194 
    184     float xx, w, z, p, q, xn; 
     195    float x, w, z, p, q, xn; 
    185196 
    186197    const float Z1 = 1.46819706421238932572E1; 
    187     const float THPIO4F =  2.35619449019234492885;    /* 3*pi/4 */ 
    188  
    189  
    190     xx = x; 
    191     if( xx < 0 ) 
    192         xx = -x; 
    193  
    194     if( xx <= 2.0 ) { 
    195         z = xx * xx; 
    196         p = (z-Z1) * xx * polevl( z, JPJ1, 4 ); 
    197         return( p ); 
     198 
     199 
     200    // 2017-05-18 PAK - mathematica and mpmath use J1(-x) = -J1(x) 
     201    x = xx; 
     202    if( x < 0 ) 
     203        x = -xx; 
     204 
     205    if( x <= 2.0 ) { 
     206        z = x * x; 
     207        p = (z-Z1) * x * polevl( z, JPJ1, 4 ); 
     208        return( xx < 0. ? -p : p ); 
    198209    } 
    199210 
    … …  
    203214    p = w * polevl( q, MO1J1, 7); 
    204215    w = q*q; 
    205     xn = q * polevl( w, PH1J1, 7) - THPIO4F; 
    206     p = p * cos(xn + xx); 
    207  
    208     return(p); 
     216    // 2017-05-19 PAK improve accuracy using trig identies 
     217    // original: 
     218    //    const float THPIO4F =  2.35619449019234492885;    /* 3*pi/4 */ 
     219    //    xn = q * polevl( w, PH1J1, 7) - THPIO4F; 
     220    //    p = p * cos(xn + x); 
     221    //    return( xx < 0. ? -p : p ); 
     222    // expanding cos(a + b - 3 pi/4) is 
     223    //    [sin(a)sin(b) + sin(a)cos(b) + cos(a)sin(b)-cos(a)cos(b)] / sqrt(2) 
     224    xn = q * polevl( w, PH1J1, 7); 
     225    float cos_xn, sin_xn; 
     226    float cos_x, sin_x; 
     227    SINCOS(xn, sin_xn, cos_xn);  // about xn and 1 
     228    SINCOS(x, sin_x, cos_x); 
     229    p *= M_SQRT1_2*(sin_xn*(sin_x+cos_x) + cos_xn*(sin_x-cos_x)); 
     230 
     231    return( xx < 0. ? -p : p ); 
    209232} 
    210233#endif 

  • sasmodels/models/lib/sas_erf.c

    rb3796fa reb2946f  
    165165        // the erf function instead and z < 1.0. 
    166166        //return (1.0 - cephes_erf(a)); 
    167         z = x * x; 
    168         y = x * polevl(z, TD, 4) / p1evl(z, UD, 5); 
    169  
    170         return y; 
     167        // 2017-05-18 PAK - use erf(a) rather than erf(|a|) 
     168        z = a * a; 
     169        y = a * polevl(z, TD, 4) / p1evl(z, UD, 5); 
     170 
     171        return 1.0 - y; 
    171172    } 
    172173 
    … …  
    279280        //is explicit here for the case < 1.0 
    280281        //return (1.0 - sas_erf(a)); 
    281         z = x * x; 
    282         y = x * polevl( z, TF, 6 ); 
    283  
    284         return y; 
     282        // 2017-05-18 PAK - use erf(a) rather than erf(|a|) 
     283        z = a * a; 
     284        y = a * polevl( z, TF, 6 ); 
     285 
     286        return 1.0 - y; 
    285287    } 
    286288 

  • sasmodels/models/stacked_disks.py

    r9802ab3 r9d50a1e  
    156156qx = q*cos(pi/6.0) 
    157157qy = q*sin(pi/6.0) 
    158 tests = [ 
    159158# Accuracy tests based on content in test/utest_extra_models.py. 
    160159# Added 2 tests with n_stacked = 5 using SasView 3.1.2 - PDB; 
    161160# for which alas q=0.001 values seem closer to n_stacked=1 not 5, 
    162161# changed assuming my 4.1 code OK, RKH 
    163     [{'thick_core': 10.0, 
    164       'thick_layer': 15.0, 
    165       'radius': 3000.0, 
    166       'n_stacking': 1.0, 
    167       'sigma_d': 0.0, 
    168       'sld_core': 4.0, 
    169       'sld_layer': -0.4, 
    170       'solvent_sd': 5.0, 
     162tests = [ 
     163    [{'thick_core': 10.0, 
     164      'thick_layer': 15.0, 
     165      'radius': 3000.0, 
     166      'n_stacking': 1.0, 
     167      'sigma_d': 0.0, 
     168      'sld_core': 4.0, 
     169      'sld_layer': -0.4, 
     170      'sld_solvent': 5.0, 
    171171      'theta': 90.0, 
    172172      'phi': 0.0, 
    … …  
    181181      'sld_core': 4.0, 
    182182      'sld_layer': -0.4, 
    183       'solvent_sd': 5.0, 
     183      'sld_solvent': 5.0, 
    184184      'theta': 90.0, 
    185185      'phi': 0.0, 
    … …  
    196196      'sld_core': 4.0, 
    197197      'sld_layer': -0.4, 
    198       'solvent_sd': 5.0, 
     198      'sld_solvent': 5.0, 
    199199      'theta': 90.0, 
    200200      'phi': 20.0, 
    201201      'scale': 0.01, 
    202       'background': 0.001}, 
    203       (qx, qy), 0.0491167089952  ], 
     202      'background': 0.001, 
     203     }, (qx, qy), 0.0491167089952], 
    204204    [{'thick_core': 10.0, 
    205205      'thick_layer': 15.0, 
    … …  
    209209      'sld_core': 4.0, 
    210210      'sld_layer': -0.4, 
    211       'solvent_sd': 5.0, 
     211      'sld_solvent': 5.0, 
    212212      'theta': 90.0, 
    213213      'phi': 0.0, 
    … …  
    225225      'sld_core': 4.0, 
    226226      'sld_layer': -0.4, 
    227       'solvent_sd': 5.0, 
     227      'sld_solvent': 5.0, 
    228228      'theta': 90.0, 
    229229      'phi': 0.0, 
    … …  
    240240      'sld_core': 4.0, 
    241241      'sld_layer': -0.4, 
    242       'solvent_sd': 5.0, 
     242      'sld_solvent': 5.0, 
    243243      'theta': 90.0, 
    244244      'phi': 0.0, 
    … …  
    246246      'background': 0.001, 
    247247     }, ([0.4, 0.5]), [0.00105074, 0.00121761]], 
    248     [{'thick_core': 10.0, 
    249       'thick_layer': 15.0, 
    250       'radius': 3000.0, 
    251       'n_stacking': 1.0, 
    252       'sigma_d': 0.0, 
    253       'sld_core': 4.0, 
    254       'sld_layer': -0.4, 
    255       'solvent_sd': 5.0, 
    256       'theta': 90.0, 
    257       'phi': 20.0, 
    258       'scale': 0.01, 
    259       'background': 0.001, 
    260      }, (qx, qy), 0.0341738733124 ], 
    261  
    262     [{'thick_core': 10.0, 
    263       'thick_layer': 15.0, 
    264       'radius': 3000.0, 
    265       'n_stacking': 1.0, 
    266       'sigma_d': 0.0, 
    267       'sld_core': 4.0, 
    268       'sld_layer': -0.4, 
    269      'solvent_sd': 5.0, 
     248#    [{'thick_core': 10.0, 
     249#      'thick_layer': 15.0, 
     250#      'radius': 3000.0, 
     251#      'n_stacking': 1.0, 
     252#      'sigma_d': 0.0, 
     253#      'sld_core': 4.0, 
     254#      'sld_layer': -0.4, 
     255#      'sld_solvent': 5.0, 
     256#      'theta': 90.0, 
     257#      'phi': 20.0, 
     258#      'scale': 0.01, 
     259#      'background': 0.001, 
     260# 2017-05-18 PAK temporarily suppress output of qx,qy test; j1 is not accurate for large qr 
     261#     }, (qx, qy), 0.0341738733124], 
     262#     }, (qx, qy), None], 
     263 
     264    [{'thick_core': 10.0, 
     265      'thick_layer': 15.0, 
     266      'radius': 3000.0, 
     267      'n_stacking': 1.0, 
     268      'sigma_d': 0.0, 
     269      'sld_core': 4.0, 
     270      'sld_layer': -0.4, 
     271      'sld_solvent': 5.0, 
    270272      'theta': 90.0, 
    271273      'phi': 0.0, 
    … …  
    274276     }, ([1.3, 1.57]), [0.0010039, 0.0010038]], 
    275277    ] 
    276 # 11Jan2017   RKH checking unit test agai, note they are all 1D, no 2D 
    277  
     278# 11Jan2017   RKH checking unit test again, note they are all 1D, no 2D 

  • sasmodels/models/two_power_law.py

    r40a87fa rbb4b509  
    120120     }, 0.150141, 0.125945], 
    121121 
    122     [{'coeffcent_1':    1.0, 
     122    [{'coefficent_1':    1.0, 
    123123      'crossover':  0.04, 
    124124      'power_1':    1.0, 
    … …  
    127127     }, 0.442528, 0.00166884], 
    128128 
    129     [{'coeffcent_1':    1.0, 
     129    [{'coefficent_1':    1.0, 
    130130      'crossover':  0.04, 
    131131      'power_1':    1.0, 

  • explore/jitter.py

    r85190c2 rd4c33d6  
    33dispersity and possible replacement algorithms. 
    44""" 
     5import sys 
     6 
    57import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot 
    68import matplotlib.pyplot as plt 
    79from matplotlib.widgets import Slider, CheckButtons 
    810from matplotlib import cm 
    9  
    1011import numpy as np 
    1112from numpy import pi, cos, sin, sqrt, exp, degrees, radians 
    1213 
    13 def draw_beam(ax): 
     14def draw_beam(ax, view=(0, 0)): 
    1415    #ax.plot([0,0],[0,0],[1,-1]) 
    1516    #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
    … …  
    2223    x = r*np.outer(np.cos(u), np.ones_like(v)) 
    2324    y = r*np.outer(np.sin(u), np.ones_like(v)) 
    24     z = np.outer(np.ones_like(u), v) 
     25    z = 1.3*np.outer(np.ones_like(u), v) 
     26 
     27    theta, phi = view 
     28    shape = x.shape 
     29    points = np.matrix([x.flatten(), y.flatten(), z.flatten()]) 
     30    points = Rz(phi)*Ry(theta)*points 
     31    x, y, z = [v.reshape(shape) for v in points] 
    2532 
    2633    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
    2734 
    28 def draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi): 
    29     size=[0.1, 0.4, 1.0] 
    30     view=[theta, phi, psi] 
    31     shimmy=[0,0,0] 
    32     #draw_shape = draw_parallelepiped 
    33     draw_shape = draw_ellipsoid 
     35def draw_jitter(ax, view, jitter): 
     36    size = [0.1, 0.4, 1.0] 
     37    draw_shape = draw_parallelepiped 
     38    #draw_shape = draw_ellipsoid 
    3439 
    3540    #np.random.seed(10) 
    … …  
    6469        [ 1,  1,  1], 
    6570    ] 
     71    dtheta, dphi, dpsi = jitter 
    6672    if dtheta == 0: 
    6773        cloud = [v for v in cloud if v[0] == 0] 
    … …  
    7076    if dpsi == 0: 
    7177        cloud = [v for v in cloud if v[2] == 0] 
    72     draw_shape(ax, size, view, shimmy, steps=100, alpha=0.8) 
     78    draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8) 
    7379    for point in cloud: 
    74         shimmy=[dtheta*point[0], dphi*point[1], dpsi*point[2]] 
    75         draw_shape(ax, size, view, shimmy, alpha=0.8) 
     80        delta = [dtheta*point[0], dphi*point[1], dpsi*point[2]] 
     81        draw_shape(ax, size, view, delta, alpha=0.8) 
    7682    for v in 'xyz': 
    7783        a, b, c = size 
    … …  
    8086        getattr(ax, v+'axis').label.set_text(v) 
    8187 
    82 def draw_ellipsoid(ax, size, view, shimmy, steps=25, alpha=1): 
     88def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): 
    8389    a,b,c = size 
    84     theta, phi, psi = view 
    85     dtheta, dphi, dpsi = shimmy 
    86  
    8790    u = np.linspace(0, 2 * np.pi, steps) 
    8891    v = np.linspace(0, np.pi, steps) 
    … …  
    9093    y = b*np.outer(np.sin(u), np.sin(v)) 
    9194    z = c*np.outer(np.ones_like(u), np.cos(v)) 
    92  
    93     shape = x.shape 
    94     points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
    95     points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points 
    96     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
    97     x,y,z = [v.reshape(shape) for v in points] 
     95    x, y, z = transform_xyz(view, jitter, x, y, z) 
    9896 
    9997    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
    10098 
    101 def draw_parallelepiped(ax, size, view, shimmy, alpha=1): 
     99    draw_labels(ax, view, jitter, [ 
     100         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     101         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     102         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     103         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     104         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     105         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     106    ]) 
     107 
     108def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): 
    102109    a,b,c = size 
    103     theta, phi, psi = view 
    104     dtheta, dphi, dpsi = shimmy 
    105  
    106110    x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    107111    y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
    … …  
    118122    ]) 
    119123 
    120     points = np.matrix([x,y,z]) 
    121     points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points 
    122     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
    123  
    124     x,y,z = [np.array(v).flatten() for v in points] 
     124    x, y, z = transform_xyz(view, jitter, x, y, z) 
    125125    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
    126126 
    127 def draw_sphere(ax, radius=10., steps=100): 
    128     u = np.linspace(0, 2 * np.pi, steps) 
    129     v = np.linspace(0, np.pi, steps) 
    130  
    131     x = radius * np.outer(np.cos(u), np.sin(v)) 
    132     y = radius * np.outer(np.sin(u), np.sin(v)) 
    133     z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 
    134     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    135  
    136 def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'): 
    137     theta_center = radians(theta) 
    138     phi_center = radians(phi) 
    139     flow_center = radians(flow) 
    140     dtheta = radians(dtheta) 
    141     dphi = radians(dphi) 
    142  
    143     # 10 point 3-sigma gaussian weights 
    144     t = np.linspace(-3., 3., 11) 
     127    draw_labels(ax, view, jitter, [ 
     128         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     129         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     130         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     131         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     132         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     133         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     134    ]) 
     135 
     136def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gauss'): 
     137    theta, phi, psi = view 
     138    dtheta, dphi, dpsi = jitter 
    145139    if dist == 'gauss': 
     140        t = np.linspace(-3, 3, n) 
    146141        weights = exp(-0.5*t**2) 
    147142    elif dist == 'rect': 
     143        t = np.linspace(0, 1, n) 
    148144        weights = np.ones_like(t) 
    149145    else: 
    150146        raise ValueError("expected dist to be 'gauss' or 'rect'") 
    151     theta = dtheta*t 
    152     phi = dphi*t 
    153  
    154     x = radius * np.outer(cos(phi), cos(theta)) 
    155     y = radius * np.outer(sin(phi), cos(theta)) 
    156     z = radius * np.outer(np.ones_like(phi), sin(theta)) 
    157     #w = np.outer(weights, weights*abs(cos(dtheta*t))) 
    158     w = np.outer(weights, weights*abs(cos(theta))) 
    159  
    160     x, y, z, w = [v.flatten() for v in (x,y,z,w)] 
    161     x, y, z = rotate(x, y, z, phi_center, theta_center, flow_center) 
    162  
    163     ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.) 
    164  
    165 def rotate(x, y, z, phi, theta, psi): 
    166     R = Rz(phi)*Ry(theta)*Rz(psi) 
    167     p = np.vstack([x,y,z]) 
    168     return R*p 
     147 
     148    # mesh in theta, phi formed by rotating z 
     149    z = np.matrix([[0], [0], [radius]]) 
     150    points = np.hstack([Rx(phi_i)*Ry(theta_i)*z 
     151                        for theta_i in dtheta*t 
     152                        for phi_i in dphi*t]) 
     153    # rotate relative to beam 
     154    points = orient_relative_to_beam(view, points) 
     155 
     156    w = np.outer(weights, weights) 
     157 
     158    x, y, z = [np.array(v).flatten() for v in points] 
     159    ax.scatter(x, y, z, c=w.flatten(), marker='o', vmin=0., vmax=1.) 
    169160 
    170161def Rx(angle): 
    … …  
    188179         [0., 0., 1.]] 
    189180    return np.matrix(R) 
     181 
     182def transform_xyz(view, jitter, x, y, z): 
     183    x, y, z = [np.asarray(v) for v in (x, y, z)] 
     184    shape = x.shape 
     185    points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
     186    points = apply_jitter(jitter, points) 
     187    points = orient_relative_to_beam(view, points) 
     188    x, y, z = [np.array(v).reshape(shape) for v in points] 
     189    return x, y, z 
     190 
     191def apply_jitter(jitter, points): 
     192    dtheta, dphi, dpsi = jitter 
     193    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 
     194    return points 
     195 
     196def orient_relative_to_beam(view, points): 
     197    theta, phi, psi = view 
     198    points = Rz(phi)*Ry(theta)*Rz(psi)*points 
     199    return points 
     200 
     201def draw_labels(ax, view, jitter, text): 
     202    labels, locations, orientations = zip(*text) 
     203    px, py, pz = zip(*locations) 
     204    dx, dy, dz = zip(*orientations) 
     205 
     206    px, py, pz = transform_xyz(view, jitter, px, py, pz) 
     207    dx, dy, dz = transform_xyz(view, jitter, dx, dy, dz) 
     208 
     209    for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)): 
     210        zdir = np.asarray(zdir).flatten() 
     211        ax.text(p[0], p[1], p[2], label, zdir=zdir) 
     212 
     213def draw_sphere(ax, radius=10., steps=100): 
     214    u = np.linspace(0, 2 * np.pi, steps) 
     215    v = np.linspace(0, np.pi, steps) 
     216 
     217    x = radius * np.outer(np.cos(u), np.sin(v)) 
     218    y = radius * np.outer(np.sin(u), np.sin(v)) 
     219    z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 
     220    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    190221 
    191222def main(): 
    … …  
    206237 
    207238    axcolor = 'lightgoldenrodyellow' 
     239 
    208240    axtheta  = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    209241    axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
    … …  
    212244    sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi) 
    213245    spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi) 
     246 
    214247    axdtheta  = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    215248    axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
    … …  
    217250    sdtheta = Slider(axdtheta, 'dTheta', 0, 30, valinit=dtheta) 
    218251    sdphi = Slider(axdphi, 'dPhi', 0, 30, valinit=dphi) 
    219     sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dphi) 
     252    sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dpsi) 
    220253 
    221254    def update(val, axis=None): 
    222         theta, phi, psi = stheta.val, sphi.val, spsi.val 
    223         dtheta, dphi, dpsi = sdtheta.val, sdphi.val, sdpsi.val 
     255        view = stheta.val, sphi.val, spsi.val 
     256        jitter = sdtheta.val, sdphi.val, sdpsi.val 
    224257        ax.cla() 
    225         draw_beam(ax) 
    226         draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi) 
    227         #if not axis.startswith('d'): 
    228         #    ax.view_init(elev=theta, azim=phi) 
     258        draw_beam(ax, (0, 0)) 
     259        draw_jitter(ax, view, jitter) 
     260        #draw_jitter(ax, view, (0,0,0)) 
     261        draw_mesh(ax, view, jitter) 
    229262        plt.gcf().canvas.draw() 
    230263 

  • sasmodels/details.py

    rccd5f01 rf39759c  
    1515 
    1616import numpy as np  # type: ignore 
    17 from numpy import pi, cos, sin 
     17from numpy import pi, cos, sin, radians 
    1818 
    1919try: 
    … …  
    2929 
    3030try: 
    31     from typing import List 
     31    from typing import List, Tuple, Sequence 
    3232except ImportError: 
    3333    pass 
    3434else: 
    3535    from .modelinfo import ModelInfo 
     36    from .modelinfo import ParameterTable 
    3637 
    3738 
    … …  
    5354    coordinates, the total circumference decreases as latitude varies from 
    5455    pi r^2 at the equator to 0 at the pole, and the weight associated 
    55     with a range of phi values needs to be scaled by this circumference. 
     56    with a range of latitude values needs to be scaled by this circumference. 
    5657    This scale factor needs to be updated each time the theta value 
    5758    changes.  *theta_par* indicates which of the values in the parameter 
    … …  
    231232    nvalues = kernel.info.parameters.nvalues 
    232233    scalars = [(v[0] if len(v) else np.NaN) for v, w in pairs] 
    233     values, weights = zip(*pairs[2:npars+2]) if npars else ((),()) 
     234    # skipping scale and background when building values and weights 
     235    values, weights = zip(*pairs[2:npars+2]) if npars else ((), ()) 
     236    weights = correct_theta_weights(kernel.info.parameters, values, weights) 
    234237    length = np.array([len(w) for w in weights]) 
    235238    offset = np.cumsum(np.hstack((0, length))) 
    … …  
    244247    return call_details, data, is_magnetic 
    245248 
     249def correct_theta_weights(parameters, values, weights): 
     250    # type: (ParameterTable, Sequence[np.ndarray], Sequence[np.ndarray]) -> Sequence[np.ndarray] 
     251    """ 
     252    If there is a theta parameter, update the weights of that parameter so that 
     253    the cosine weighting required for polar integration is preserved.  Avoid 
     254    evaluation strictly at the pole, which would otherwise send the weight to 
     255    zero. 
     256 
     257    Note: values and weights do not include scale and background 
     258    """ 
     259    # TODO: document code, explaining why scale and background are skipped 
     260    # given that we don't have scale and background in the list, we 
     261    # should be calling the variables something other than values and weights 
     262    # Apparently the parameters.theta_offset similarly skips scale and 
     263    # and background, so the indexing works out. 
     264    if parameters.theta_offset >= 0: 
     265        index = parameters.theta_offset 
     266        theta = values[index] 
     267        theta_weight = np.minimum(abs(cos(radians(theta))), 1e-6) 
     268        # copy the weights list so we can update it 
     269        weights = list(weights) 
     270        weights[index] = theta_weight*np.asarray(weights[index]) 
     271        weights = tuple(weights) 
     272    return weights 
     273 
    246274 
    247275def convert_magnetism(parameters, values): 
     276    # type: (ParameterTable, Sequence[np.ndarray]) -> bool 
    248277    """ 
    249278    Convert magnetism values from polar to rectangular coordinates. 
    … …  
    255284    scale = mag[:,0] 
    256285    if np.any(scale): 
    257         theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180. 
     286        theta, phi = radians(mag[:, 1]), radians(mag[:, 2]) 
    258287        cos_theta = cos(theta) 
    259288        mag[:, 0] = scale*cos_theta*cos(phi)  # mx 
    … …  
    269298    """ 
    270299    Create a mesh grid of dispersion parameters and weights. 
     300 
     301    pars is a list of pairs (values, weights), where the values are the 
     302    individual parameter values at which to evaluate the polydispersity 
     303    distribution and weights are the weights associated with each value. 
     304 
     305    Only the volume parameters should be included in this list.  Orientation 
     306    parameters do not affect the calculation of effective radius or volume 
     307    ratio. 
    271308 
    272309    Returns [p1,p2,...],w where pj is a vector of values for parameter j 

  • sasmodels/kernel_iq.c

    rbde38b5 rd4c33d6  
    2525    int32_t num_weights;        // total length of the weights vector 
    2626    int32_t num_active;         // number of non-trivial pd loops 
    27     int32_t theta_par;          // id of spherical correction variable 
     27    int32_t theta_par;          // id of spherical correction variable (not used) 
    2828} ProblemDetails; 
    2929 
    … …  
    173173 
    174174 
    175 #if MAX_PD>0 
    176   const int theta_par = details->theta_par; 
    177   const int fast_theta = (theta_par == p0); 
    178   const int slow_theta = (theta_par >= 0 && !fast_theta); 
    179   double spherical_correction = 1.0; 
    180 #else 
    181   // Note: if not polydisperse the weights cancel and we don't need the 
    182   // spherical correction. 
    183   const double spherical_correction = 1.0; 
    184 #endif 
    185  
    186175  int step = pd_start; 
    187176 
    … …  
    220209#endif 
    221210#if MAX_PD>0 
    222   if (slow_theta) { // Theta is not in inner loop 
    223     spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6); 
    224   } 
    225211  while(i0 < n0) { 
    226212    local_values.vector[p0] = v0[i0]; 
    227213    double weight0 = w0[i0] * weight1; 
    228214//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0); 
    229     if (fast_theta) { // Theta is in inner loop 
    230       spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6); 
    231     } 
    232215#else 
    233216    const double weight0 = 1.0; 
    … …  
    244227      // Note: weight==0 must always be excluded 
    245228      if (weight0 > cutoff) { 
    246         // spherical correction is set at a minimum of 1e-6, otherwise there 
    247         // would be problems looking at models with theta=90. 
    248         const double weight = weight0 * spherical_correction; 
    249         pd_norm += weight * CALL_VOLUME(local_values.table); 
     229        pd_norm += weight0 * CALL_VOLUME(local_values.table); 
    250230 
    251231        #ifdef USE_OPENMP 
    … …  
    304284#endif // !MAGNETIC 
    305285//printf("q_index:%d %g %g %g %g\n",q_index, scattering, weight, spherical_correction, weight0); 
    306           result[q_index] += weight * scattering; 
     286          result[q_index] += weight0 * scattering; 
    307287        } 
    308288      } 

  • sasmodels/kernel_iq.cl

    rbde38b5 rd4c33d6  
    2525    int32_t num_weights;        // total length of the weights vector 
    2626    int32_t num_active;         // number of non-trivial pd loops 
    27     int32_t theta_par;          // id of spherical correction variable 
     27    int32_t theta_par;          // id of spherical correction variable (not used) 
    2828} ProblemDetails; 
    2929 
    … …  
    169169 
    170170 
    171 #if MAX_PD>0 
    172   const int theta_par = details->theta_par; 
    173   const bool fast_theta = (theta_par == p0); 
    174   const bool slow_theta = (theta_par >= 0 && !fast_theta); 
    175   double spherical_correction = 1.0; 
    176 #else 
    177   // Note: if not polydisperse the weights cancel and we don't need the 
    178   // spherical correction. 
    179   const double spherical_correction = 1.0; 
    180 #endif 
    181  
    182171  int step = pd_start; 
    183172 
    … …  
    217206#endif 
    218207#if MAX_PD>0 
    219   if (slow_theta) { // Theta is not in inner loop 
    220     spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6); 
    221   } 
    222208  while(i0 < n0) { 
    223209    local_values.vector[p0] = v0[i0]; 
    224210    double weight0 = w0[i0] * weight1; 
    225211//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0); 
    226     if (fast_theta) { // Theta is in inner loop 
    227       spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6); 
    228     } 
    229212#else 
    230213    const double weight0 = 1.0; 
    … …  
    232215 
    233216//if (q_index == 0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, local_values.vector[i]); printf("\n"); } 
    234 //if (q_index == 0) printf("sphcor: %g\n", spherical_correction); 
    235217 
    236218    #ifdef INVALID 
    … …  
    241223      // Note: weight==0 must always be excluded 
    242224      if (weight0 > cutoff) { 
    243         // spherical correction is set at a minimum of 1e-6, otherwise there 
    244         // would be problems looking at models with theta=90. 
    245         const double weight = weight0 * spherical_correction; 
    246         pd_norm += weight * CALL_VOLUME(local_values.table); 
     225        pd_norm += weight0 * CALL_VOLUME(local_values.table); 
    247226 
    248227#if defined(MAGNETIC) && NUM_MAGNETIC > 0 
    … …  
    296275        const double scattering = CALL_IQ(q, q_index, local_values.table); 
    297276#endif // !MAGNETIC 
    298         this_result += weight * scattering; 
     277        this_result += weight0 * scattering; 
    299278      } 
    300279    } 

  • sasmodels/kernelpy.py

    rdbfd471 rd4c33d6  
    197197 
    198198    pd_norm = 0.0 
    199     spherical_correction = 1.0 
    200199    partial_weight = np.NaN 
    201200    weight = np.NaN 
    202201 
    203202    p0_par = call_details.pd_par[0] 
    204     p0_is_theta = (p0_par == call_details.theta_par) 
    205203    p0_length = call_details.pd_length[0] 
    206204    p0_index = p0_length 
    … …  
    219217            parameters[pd_par] = pd_value[pd_offset+pd_index] 
    220218            partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:]) 
    221             if call_details.theta_par >= 0: 
    222                 cor = sin(pi / 180 * parameters[call_details.theta_par]) 
    223                 spherical_correction = max(abs(cor), 1e-6) 
    224219            p0_index = loop_index%p0_length 
    225220 
    226221        weight = partial_weight * pd_weight[p0_offset + p0_index] 
    227222        parameters[p0_par] = pd_value[p0_offset + p0_index] 
    228         if p0_is_theta: 
    229             cor = cos(pi/180 * parameters[p0_par]) 
    230             spherical_correction = max(abs(cor), 1e-6) 
    231223        p0_index += 1 
    232224        if weight > cutoff: 
    … …  
    239231 
    240232            # update value and norm 
    241             weight *= spherical_correction 
    242233            total += weight * Iq 
    243234            pd_norm += weight * form_volume() 

  • sasmodels/models/barbell.c

    r592343f r2a0b2b1  
    1010//barbell kernel - same as dumbell 
    1111static double 
    12 _bell_kernel(double q, double h, double radius_bell, 
    13              double half_length, double sin_alpha, double cos_alpha) 
     12_bell_kernel(double qab, double qc, double h, double radius_bell, 
     13             double half_length) 
    1414{ 
    1515    // translate a point in [-1,1] to a point in [lower,upper] 
    … …  
    2626    //    m = q R cos(alpha) 
    2727    //    b = q(L/2-h) cos(alpha) 
    28     const double m = q*radius_bell*cos_alpha; // cos argument slope 
    29     const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 
    30     const double qrst = q*radius_bell*sin_alpha; // Q*R*sin(theta) 
     28    const double m = radius_bell*qc; // cos argument slope 
     29    const double b = (half_length-h)*qc; // cos argument intercept 
     30    const double qab_r = radius_bell*qab; // Q*R*sin(theta) 
    3131    double total = 0.0; 
    3232    for (int i = 0; i < 76; i++){ 
    3333        const double t = Gauss76Z[i]*zm + zb; 
    3434        const double radical = 1.0 - t*t; 
    35         const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
     35        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    3636        const double Fq = cos(m*t + b) * radical * bj; 
    3737        total += Gauss76Wt[i] * Fq; 
    … …  
    4444 
    4545static double 
    46 _fq(double q, double h, 
    47     double radius_bell, double radius, double half_length, 
    48     double sin_alpha, double cos_alpha) 
     46_fq(double qab, double qc, double h, 
     47    double radius_bell, double radius, double half_length) 
    4948{ 
    50     const double bell_fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha); 
    51     const double bj = sas_2J1x_x(q*radius*sin_alpha); 
    52     const double si = sas_sinx_x(q*half_length*cos_alpha); 
     49    const double bell_fq = _bell_kernel(qab, qc, h, radius_bell, half_length); 
     50    const double bj = sas_2J1x_x(radius*qab); 
     51    const double si = sas_sinx_x(half_length*qc); 
    5352    const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
    5453    const double Aq = bell_fq + cyl_fq; 
    … …  
    8483        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8584        SINCOS(alpha, sin_alpha, cos_alpha); 
    86         const double Aq = _fq(q, h, radius_bell, radius, half_length, sin_alpha, cos_alpha); 
     85        const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 
    8786        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
    8887    } 
    … …  
    103102    double q, sin_alpha, cos_alpha; 
    104103    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     104    const double qab = q*sin_alpha; 
     105    const double qc = q*cos_alpha; 
    105106 
    106107    const double h = -sqrt(square(radius_bell) - square(radius)); 
    107     const double Aq = _fq(q, h, radius_bell, radius, 0.5*length, sin_alpha, cos_alpha); 
     108    const double Aq = _fq(qab, qc, h, radius_bell, radius, 0.5*length); 
    108109 
    109110    // Multiply by contrast^2 and convert to cm-1 

  • sasmodels/models/bcc_paracrystal.c

    r50beefe r7e0b281  
    1 double form_volume(double radius); 
    2 double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
    3 double Iqxy(double qx, double qy, double dnn, 
    4     double d_factor, double radius,double sld, double solvent_sld, 
    5     double theta, double phi, double psi); 
     1static double 
     2bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4#if 0  // Equations as written in Matsuoka 
     5    const double a1 = (+qa + qb + qc)/2.0; 
     6    const double a2 = (-qa - qb + qc)/2.0; 
     7    const double a3 = (-qa + qb - qc)/2.0; 
     8#else 
     9    const double a1 = (+qa + qb - qc)/2.0; 
     10    const double a2 = (+qa - qb + qc)/2.0; 
     11    const double a3 = (-qa + qb + qc)/2.0; 
     12#endif 
    613 
    7 double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
    8 double _BCCeval(double Theta, double Phi, double temp1, double temp3); 
    9 double _sphereform(double q, double radius, double sld, double solvent_sld); 
     14#if 1 
     15    // Numerator: (1 - exp(a)^2)^3 
     16    //         => (-(exp(2a) - 1))^3 
     17    //         => -expm1(2a)^3 
     18    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a)) 
     19    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
     20    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1) 
     21    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     22    const double exp_arg = exp(arg); 
     23    const double Zq = -cube(expm1(2.0*arg)) 
     24        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
     25          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
     26          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
     27#else 
     28    // Alternate form, which perhaps is more approachable 
     29    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     30    const double sinh_qd = sinh(arg); 
     31    const double cosh_qd = cosh(arg); 
     32    const double Zq = sinh_qd/(cosh_qd - cos(dnn*a1)) 
     33                    * sinh_qd/(cosh_qd - cos(dnn*a2)) 
     34                    * sinh_qd/(cosh_qd - cos(dnn*a3)); 
     35#endif 
     36 
     37    return Zq; 
     38} 
    1039 
    1140 
    12 double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { 
    13  
    14         const double Da = d_factor*dnn; 
    15         const double temp1 = q*q*Da*Da; 
    16         const double temp3 = q*dnn; 
    17  
    18         double retVal = _BCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); 
    19         return(retVal); 
     41// occupied volume fraction calculated from lattice symmetry and sphere radius 
     42static double 
     43bcc_volume_fraction(double radius, double dnn) 
     44{ 
     45    return 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 
    2046} 
    2147 
    22 double _BCCeval(double Theta, double Phi, double temp1, double temp3) { 
    23  
    24         double result; 
    25         double sin_theta,cos_theta,sin_phi,cos_phi; 
    26         SINCOS(Theta, sin_theta, cos_theta); 
    27         SINCOS(Phi, sin_phi, cos_phi); 
    28  
    29         const double temp6 =  sin_theta; 
    30         const double temp7 =  sin_theta*cos_phi + sin_theta*sin_phi + cos_theta; 
    31         const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta; 
    32         const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta; 
    33  
    34         const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9)); 
    35         result = cube(1.0 - (temp10*temp10))*temp6 
    36             / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    37               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
    38               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10)); 
    39  
    40         return (result); 
    41 } 
    42  
    43 double form_volume(double radius){ 
     48static double 
     49form_volume(double radius) 
     50{ 
    4451    return sphere_volume(radius); 
    4552} 
    4653 
    4754 
    48 double Iq(double q, double dnn, 
     55static double Iq(double q, double dnn, 
    4956  double d_factor, double radius, 
    50   double sld, double solvent_sld){ 
     57  double sld, double solvent_sld) 
     58{ 
     59    // translate a point in [-1,1] to a point in [0, 2 pi] 
     60    const double phi_m = M_PI; 
     61    const double phi_b = M_PI; 
     62    // translate a point in [-1,1] to a point in [0, pi] 
     63    const double theta_m = M_PI_2; 
     64    const double theta_b = M_PI_2; 
    5165 
    52         //Volume fraction calculated from lattice symmetry and sphere radius 
    53         const double s1 = dnn/sqrt(0.75); 
    54         const double latticescale = 2.0*sphere_volume(radius/s1); 
    55  
    56     const double va = 0.0; 
    57     const double vb = 2.0*M_PI; 
    58     const double vaj = 0.0; 
    59     const double vbj = M_PI; 
    60  
    61     double summ = 0.0; 
    62     double answer = 0.0; 
    63         for(int i=0; i<150; i++) { 
    64                 //setup inner integral over the ellipsoidal cross-section 
    65                 double summj=0.0; 
    66                 const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;             //the outer dummy is phi 
    67                 for(int j=0;j<150;j++) { 
    68                         //20 gauss points for the inner integral 
    69                         double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;             //the inner dummy is theta 
    70                         double yyy = Gauss150Wt[j] * _BCC_Integrand(q,dnn,d_factor,ztheta,zphi); 
    71                         summj += yyy; 
    72                 } 
    73                 //now calculate the value of the inner integral 
    74                 double answer = (vbj-vaj)/2.0*summj; 
    75  
    76                 //now calculate outer integral 
    77                 summ = summ+(Gauss150Wt[i] * answer); 
    78         }               //final scaling is done at the end of the function, after the NT_FP64 case 
    79  
    80         answer = (vb-va)/2.0*summ; 
    81         answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale; 
    82  
    83     return answer; 
     66    double outer_sum = 0.0; 
     67    for(int i=0; i<150; i++) { 
     68        double inner_sum = 0.0; 
     69        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     70        double sin_theta, cos_theta; 
     71        SINCOS(theta, sin_theta, cos_theta); 
     72        const double qc = q*cos_theta; 
     73        const double qab = q*sin_theta; 
     74        for(int j=0;j<150;j++) { 
     75            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     76            double sin_phi, cos_phi; 
     77            SINCOS(phi, sin_phi, cos_phi); 
     78            const double qa = qab*cos_phi; 
     79            const double qb = qab*sin_phi; 
     80            const double form = bcc_Zq(qa, qb, qc, dnn, d_factor); 
     81            inner_sum += Gauss150Wt[j] * form; 
     82        } 
     83        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     84        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     85    } 
     86    outer_sum *= theta_m; 
     87    const double Zq = outer_sum/(4.0*M_PI); 
     88    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     89    return bcc_volume_fraction(radius, dnn) * Pq * Zq; 
    8490} 
    8591 
    8692 
    87 double Iqxy(double qx, double qy, 
     93static double Iqxy(double qx, double qy, 
    8894    double dnn, double d_factor, double radius, 
    8995    double sld, double solvent_sld, 
    … …  
    9298    double q, zhat, yhat, xhat; 
    9399    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     100    const double qa = q*xhat; 
     101    const double qb = q*yhat; 
     102    const double qc = q*zhat; 
    94103 
    95     const double a1 = +xhat - zhat + yhat; 
    96     const double a2 = +xhat + zhat - yhat; 
    97     const double a3 = -xhat + zhat + yhat; 
    98  
    99     const double qd = 0.5*q*dnn; 
    100     const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
    101     const double tanh_qd = tanh(arg); 
    102     const double cosh_qd = cosh(arg); 
    103     const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd) 
    104                     * tanh_qd/(1. - cos(qd*a2)/cosh_qd) 
    105                     * tanh_qd/(1. - cos(qd*a3)/cosh_qd); 
    106  
    107     const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; 
    108     //the occupied volume of the lattice 
    109     const double lattice_scale = 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 
    110     return lattice_scale * Fq; 
     104    q = sqrt(qa*qa + qb*qb + qc*qc); 
     105    const double Zq = bcc_Zq(qa, qb, qc, dnn, d_factor); 
     106    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     107    return bcc_volume_fraction(radius, dnn) * Pq * Zq; 
    111108} 

  • sasmodels/models/bcc_paracrystal.py

    r69e1afc r2a0b2b1  
    8181.. figure:: img/parallelepiped_angle_definition.png 
    8282 
    83     Orientation of the crystal with respect to the scattering plane, when  
     83    Orientation of the crystal with respect to the scattering plane, when 
    8484    $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis). 
    8585 

  • sasmodels/models/capped_cylinder.c

    r592343f r2a0b2b1  
    1414//   radius_cap is the radius of the lens 
    1515//   length is the cylinder length, or the separation between the lens halves 
    16 //   alpha is the angle of the cylinder wrt q. 
     16//   theta is the angle of the cylinder wrt q. 
    1717static double 
    18 _cap_kernel(double q, double h, double radius_cap, 
    19                       double half_length, double sin_alpha, double cos_alpha) 
     18_cap_kernel(double qab, double qc, double h, double radius_cap, 
     19            double half_length) 
    2020{ 
    2121    // translate a point in [-1,1] to a point in [lower,upper] 
    … …  
    2626 
    2727    // cos term in integral is: 
    28     //    cos (q (R t - h + L/2) cos(alpha)) 
     28    //    cos (q (R t - h + L/2) cos(theta)) 
    2929    // so turn it into: 
    3030    //    cos (m t + b) 
    3131    // where: 
    32     //    m = q R cos(alpha) 
    33     //    b = q(L/2-h) cos(alpha) 
    34     const double m = q*radius_cap*cos_alpha; // cos argument slope 
    35     const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 
    36     const double qrst = q*radius_cap*sin_alpha; // Q*R*sin(theta) 
     32    //    m = q R cos(theta) 
     33    //    b = q(L/2-h) cos(theta) 
     34    const double m = radius_cap*qc; // cos argument slope 
     35    const double b = (half_length-h)*qc; // cos argument intercept 
     36    const double qab_r = radius_cap*qab; // Q*R*sin(theta) 
    3737    double total = 0.0; 
    3838    for (int i=0; i<76 ;i++) { 
    3939        const double t = Gauss76Z[i]*zm + zb; 
    4040        const double radical = 1.0 - t*t; 
    41         const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
     41        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    4242        const double Fq = cos(m*t + b) * radical * bj; 
    4343        total += Gauss76Wt[i] * Fq; 
    … …  
    5050 
    5151static double 
    52 _fq(double q, double h, double radius_cap, double radius, double half_length, 
    53     double sin_alpha, double cos_alpha) 
     52_fq(double qab, double qc, double h, double radius_cap, double radius, double half_length) 
    5453{ 
    55     const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha); 
    56     const double bj = sas_2J1x_x(q*radius*sin_alpha); 
    57     const double si = sas_sinx_x(q*half_length*cos_alpha); 
     54    const double cap_Fq = _cap_kernel(qab, qc, h, radius_cap, half_length); 
     55    const double bj = sas_2J1x_x(radius*qab); 
     56    const double si = sas_sinx_x(half_length*qc); 
    5857    const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
    5958    const double Aq = cap_Fq + cyl_Fq; 
    … …  
    101100    double total = 0.0; 
    102101    for (int i=0; i<76 ;i++) { 
    103         const double alpha= Gauss76Z[i]*zm + zb; 
    104         double sin_alpha, cos_alpha; // slots to hold sincos function output 
    105         SINCOS(alpha, sin_alpha, cos_alpha); 
    106  
    107         const double Aq = _fq(q, h, radius_cap, radius, half_length, sin_alpha, cos_alpha); 
    108         // sin_alpha for spherical coord integration 
    109         total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
     102        const double theta = Gauss76Z[i]*zm + zb; 
     103        double sin_theta, cos_theta; // slots to hold sincos function output 
     104        SINCOS(theta, sin_theta, cos_theta); 
     105        const double qab = q*sin_theta; 
     106        const double qc = q*cos_theta; 
     107        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 
     108        // scale by sin_theta for spherical coord integration 
     109        total += Gauss76Wt[i] * Aq * Aq * sin_theta; 
    110110    } 
    111111    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    125125    double q, sin_alpha, cos_alpha; 
    126126    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     127    const double qab = q*sin_alpha; 
     128    const double qc = q*cos_alpha; 
    127129 
    128130    const double h = sqrt(radius_cap*radius_cap - radius*radius); 
    129     const double Aq = _fq(q, h, radius_cap, radius, 0.5*length, sin_alpha, cos_alpha); 
     131    const double Aq = _fq(qab, qc, h, radius_cap, radius, 0.5*length); 
    130132 
    131133    // Multiply by contrast^2 and convert to cm-1 

  • sasmodels/models/core_shell_bicelle.c

    rb260926 r2a0b2b1  
    1 double form_volume(double radius, double thick_rim, double thick_face, double length); 
    2 double Iq(double q, 
    3           double radius, 
    4           double thick_rim, 
    5           double thick_face, 
    6           double length, 
    7           double core_sld, 
    8           double face_sld, 
    9           double rim_sld, 
    10           double solvent_sld); 
    11  
    12  
    13 double Iqxy(double qx, double qy, 
    14           double radius, 
    15           double thick_rim, 
    16           double thick_face, 
    17           double length, 
    18           double core_sld, 
    19           double face_sld, 
    20           double rim_sld, 
    21           double solvent_sld, 
    22           double theta, 
    23           double phi); 
    24  
    25  
    26 double form_volume(double radius, double thick_rim, double thick_face, double length) 
     1static double 
     2form_volume(double radius, double thick_rim, double thick_face, double length) 
    273{ 
    28     return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2.0*thick_face); 
     4    return M_PI*square(radius+thick_rim)*(length+2.0*thick_face); 
    295} 
    306 
    317static double 
    32 bicelle_kernel(double q, 
    33               double rad, 
    34               double radthick, 
    35               double facthick, 
    36               double halflength, 
    37               double rhoc, 
    38               double rhoh, 
    39               double rhor, 
    40               double rhosolv, 
    41               double sin_alpha, 
    42               double cos_alpha) 
     8bicelle_kernel(double qab, 
     9    double qc, 
     10    double radius, 
     11    double thick_radius, 
     12    double thick_face, 
     13    double halflength, 
     14    double sld_core, 
     15    double sld_face, 
     16    double sld_rim, 
     17    double sld_solvent) 
    4318{ 
    44     const double dr1 = rhoc-rhoh; 
    45     const double dr2 = rhor-rhosolv; 
    46     const double dr3 = rhoh-rhor; 
    47     const double vol1 = M_PI*square(rad)*2.0*(halflength); 
    48     const double vol2 = M_PI*square(rad+radthick)*2.0*(halflength+facthick); 
    49     const double vol3 = M_PI*square(rad)*2.0*(halflength+facthick); 
     19    const double dr1 = sld_core-sld_face; 
     20    const double dr2 = sld_rim-sld_solvent; 
     21    const double dr3 = sld_face-sld_rim; 
     22    const double vol1 = M_PI*square(radius)*2.0*(halflength); 
     23    const double vol2 = M_PI*square(radius+thick_radius)*2.0*(halflength+thick_face); 
     24    const double vol3 = M_PI*square(radius)*2.0*(halflength+thick_face); 
    5025 
    51     const double be1 = sas_2J1x_x(q*(rad)*sin_alpha); 
    52     const double be2 = sas_2J1x_x(q*(rad+radthick)*sin_alpha); 
    53     const double si1 = sas_sinx_x(q*(halflength)*cos_alpha); 
    54     const double si2 = sas_sinx_x(q*(halflength+facthick)*cos_alpha); 
     26    const double be1 = sas_2J1x_x((radius)*qab); 
     27    const double be2 = sas_2J1x_x((radius+thick_radius)*qab); 
     28    const double si1 = sas_sinx_x((halflength)*qc); 
     29    const double si2 = sas_sinx_x((halflength+thick_face)*qc); 
    5530 
    5631    const double t = vol1*dr1*si1*be1 + 
    … …  
    5833                     vol3*dr3*si2*be1; 
    5934 
    60     const double retval = t*t; 
    61  
    62     return retval; 
    63  
     35    return t; 
    6436} 
    6537 
    6638static double 
    67 bicelle_integration(double q, 
    68                    double rad, 
    69                    double radthick, 
    70                    double facthick, 
    71                    double length, 
    72                    double rhoc, 
    73                    double rhoh, 
    74                    double rhor, 
    75                    double rhosolv) 
     39Iq(double q, 
     40    double radius, 
     41    double thick_radius, 
     42    double thick_face, 
     43    double length, 
     44    double sld_core, 
     45    double sld_face, 
     46    double sld_rim, 
     47    double sld_solvent) 
    7648{ 
    7749    // set up the integration end points 
    … …  
    7951    const double halflength = 0.5*length; 
    8052 
    81     double summ = 0.0; 
     53    double total = 0.0; 
    8254    for(int i=0;i<N_POINTS_76;i++) { 
    83         double alpha = (Gauss76Z[i] + 1.0)*uplim; 
    84         double sin_alpha, cos_alpha; // slots to hold sincos function output 
    85         SINCOS(alpha, sin_alpha, cos_alpha); 
    86         double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick, 
    87                              halflength, rhoc, rhoh, rhor, rhosolv, 
    88                              sin_alpha, cos_alpha); 
    89         summ += yyy*sin_alpha; 
     55        double theta = (Gauss76Z[i] + 1.0)*uplim; 
     56        double sin_theta, cos_theta; // slots to hold sincos function output 
     57        SINCOS(theta, sin_theta, cos_theta); 
     58        double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 
     59                                   halflength, sld_core, sld_face, sld_rim, sld_solvent); 
     60        total += Gauss76Wt[i]*fq*fq*sin_theta; 
    9061    } 
    9162 
    9263    // calculate value of integral to return 
    93     double answer = uplim*summ; 
    94     return answer; 
     64    double answer = total*uplim; 
     65    return 1.0e-4*answer; 
    9566} 
    9667 
    9768static double 
    98 bicelle_kernel_2d(double qx, double qy, 
    99           double radius, 
    100           double thick_rim, 
    101           double thick_face, 
    102           double length, 
    103           double core_sld, 
    104           double face_sld, 
    105           double rim_sld, 
    106           double solvent_sld, 
    107           double theta, 
    108           double phi) 
     69Iqxy(double qx, double qy, 
     70    double radius, 
     71    double thick_rim, 
     72    double thick_face, 
     73    double length, 
     74    double core_sld, 
     75    double face_sld, 
     76    double rim_sld, 
     77    double solvent_sld, 
     78    double theta, 
     79    double phi) 
    10980{ 
    11081    double q, sin_alpha, cos_alpha; 
    11182    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     83    const double qab = q*sin_alpha; 
     84    const double qc = q*cos_alpha; 
    11285 
    113     double answer = bicelle_kernel(q, radius, thick_rim, thick_face, 
     86    double fq = bicelle_kernel(qab, qc, radius, thick_rim, thick_face, 
    11487                           0.5*length, core_sld, face_sld, rim_sld, 
    115                            solvent_sld, sin_alpha, cos_alpha); 
    116     return 1.0e-4*answer; 
     88                           solvent_sld); 
     89    return 1.0e-4*fq*fq; 
    11790} 
    118  
    119 double Iq(double q, 
    120           double radius, 
    121           double thick_rim, 
    122           double thick_face, 
    123           double length, 
    124           double core_sld, 
    125           double face_sld, 
    126           double rim_sld, 
    127           double solvent_sld) 
    128 { 
    129     double intensity = bicelle_integration(q, radius, thick_rim, thick_face, 
    130                        length, core_sld, face_sld, rim_sld, solvent_sld); 
    131     return intensity*1.0e-4; 
    132 } 
    133  
    134  
    135 double Iqxy(double qx, double qy, 
    136           double radius, 
    137           double thick_rim, 
    138           double thick_face, 
    139           double length, 
    140           double core_sld, 
    141           double face_sld, 
    142           double rim_sld, 
    143           double solvent_sld, 
    144           double theta, 
    145           double phi) 
    146 { 
    147     double intensity = bicelle_kernel_2d(qx, qy, 
    148                       radius, 
    149                       thick_rim, 
    150                       thick_face, 
    151                       length, 
    152                       core_sld, 
    153                       face_sld, 
    154                       rim_sld, 
    155                       solvent_sld, 
    156                       theta, 
    157                       phi); 
    158  
    159     return intensity; 
    160 } 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    rdedcf34 r2a0b2b1  
    1717        double thick_face, 
    1818        double length, 
    19         double rhoc, 
    20         double rhoh, 
    21         double rhor, 
    22         double rhosolv) 
     19        double sld_core, 
     20        double sld_face, 
     21        double sld_rim, 
     22        double sld_solvent) 
    2323{ 
    24     double si1,si2,be1,be2; 
    2524     // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle 
    2625     // tested against limiting cases of cylinder, elliptical_cylinder, stacked_discs, and core_shell_bicelle 
    27      //    const double uplim = M_PI_4; 
    2826    const double halfheight = 0.5*length; 
    29     //const double va = 0.0; 
    30     //const double vb = 1.0; 
    31     // inner integral limits 
    32     //const double vaj=0.0; 
    33     //const double vbj=M_PI; 
    34  
    3527    const double r_major = r_minor * x_core; 
    3628    const double r2A = 0.5*(square(r_major) + square(r_minor)); 
    3729    const double r2B = 0.5*(square(r_major) - square(r_minor)); 
    38     const double dr1 = (rhoc-rhoh)   *M_PI*r_minor*r_major*(2.0*halfheight);; 
    39     const double dr2 = (rhor-rhosolv)*M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
    40     const double dr3 = (rhoh-rhor)   *M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
    41     //const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight); 
    42     //const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
    43     //const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
     30    const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight); 
     31    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
     32    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
     33    const double dr1 = vol1*(sld_core-sld_face); 
     34    const double dr2 = vol2*(sld_rim-sld_solvent); 
     35    const double dr3 = vol3*(sld_face-sld_rim); 
    4436 
    4537    //initialize integral 
    … …  
    4739    for(int i=0;i<76;i++) { 
    4840        //setup inner integral over the ellipsoidal cross-section 
    49         // since we generate these lots of times, why not store them somewhere? 
    50         //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
    51         const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0; 
    52         const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
    53         double inner_sum=0; 
    54         double sinarg1 = q*halfheight*cos_alpha; 
    55         double sinarg2 = q*(halfheight+thick_face)*cos_alpha; 
    56         si1 = sas_sinx_x(sinarg1); 
    57         si2 = sas_sinx_x(sinarg2); 
     41        //const double va = 0.0; 
     42        //const double vb = 1.0; 
     43        //const double cos_theta = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
     44        const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.0; 
     45        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     46        const double qab = q*sin_theta; 
     47        const double qc = q*cos_theta; 
     48        const double si1 = sas_sinx_x(halfheight*qc); 
     49        const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 
     50        double inner_sum=0.0; 
    5851        for(int j=0;j<76;j++) { 
    5952            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    60             //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    61             const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 
    62             const double rr = sqrt(r2A - r2B*cos(beta)); 
    63             double besarg1 = q*rr*sin_alpha; 
    64             double besarg2 = q*(rr+thick_rim)*sin_alpha; 
    65             be1 = sas_2J1x_x(besarg1); 
    66             be2 = sas_2J1x_x(besarg2); 
    67             inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 + 
    68                                               dr2*si2*be2 + 
    69                                               dr3*si2*be1); 
     53            // inner integral limits 
     54            //const double vaj=0.0; 
     55            //const double vbj=M_PI; 
     56            //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     57            const double phi = ( Gauss76Z[j] +1.0)*M_PI_2; 
     58            const double rr = sqrt(r2A - r2B*cos(phi)); 
     59            const double be1 = sas_2J1x_x(rr*qab); 
     60            const double be2 = sas_2J1x_x((rr+thick_rim)*qab); 
     61            const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 
     62 
     63            inner_sum += Gauss76Wt[j] * fq * fq; 
    7064        } 
    7165        //now calculate outer integral 
    … …  
    8377          double thick_face, 
    8478          double length, 
    85           double rhoc, 
    86           double rhoh, 
    87           double rhor, 
    88           double rhosolv, 
     79          double sld_core, 
     80          double sld_face, 
     81          double sld_rim, 
     82          double sld_solvent, 
    8983          double theta, 
    9084          double phi, 
    9185          double psi) 
    9286{ 
    93        // THIS NEEDS TESTING 
    9487    double q, xhat, yhat, zhat; 
    9588    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    96     const double dr1 = rhoc-rhoh; 
    97     const double dr2 = rhor-rhosolv; 
    98     const double dr3 = rhoh-rhor; 
     89    const double qa = q*xhat; 
     90    const double qb = q*yhat; 
     91    const double qc = q*zhat; 
     92 
     93    const double dr1 = sld_core-sld_face; 
     94    const double dr2 = sld_rim-sld_solvent; 
     95    const double dr3 = sld_face-sld_rim; 
    9996    const double r_major = r_minor*x_core; 
    10097    const double halfheight = 0.5*length; 
    … …  
    104101 
    105102    // Compute effective radius in rotated coordinates 
    106     const double r_hat = sqrt(square(r_major*xhat) + square(r_minor*yhat)); 
    107     const double rshell_hat = sqrt(square((r_major+thick_rim)*xhat) 
    108                                    + square((r_minor+thick_rim)*yhat)); 
    109     const double be1 = sas_2J1x_x( q*r_hat ); 
    110     const double be2 = sas_2J1x_x( q*rshell_hat ); 
    111     const double si1 = sas_sinx_x( q*halfheight*zhat ); 
    112     const double si2 = sas_sinx_x( q*(halfheight + thick_face)*zhat ); 
    113     const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1); 
    114     return 1.0e-4 * Aq; 
     103    const double qr_hat = sqrt(square(r_major*qa) + square(r_minor*qb)); 
     104    const double qrshell_hat = sqrt(square((r_major+thick_rim)*qa) 
     105                                   + square((r_minor+thick_rim)*qb)); 
     106    const double be1 = sas_2J1x_x( qr_hat ); 
     107    const double be2 = sas_2J1x_x( qrshell_hat ); 
     108    const double si1 = sas_sinx_x( halfheight*qc ); 
     109    const double si2 = sas_sinx_x( (halfheight + thick_face)*qc ); 
     110    const double fq = vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1; 
     111    return 1.0e-4 * fq*fq; 
    115112} 
    116  

  • sasmodels/models/core_shell_cylinder.c

    r592343f r2a0b2b1  
    1 double form_volume(double radius, double thickness, double length); 
    2 double Iq(double q, double core_sld, double shell_sld, double solvent_sld, 
    3     double radius, double thickness, double length); 
    4 double 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); 
    6  
    71// vd = volume * delta_rho 
    8 // besarg = q * R * sin(alpha) 
    9 // siarg = q * L/2 * cos(alpha) 
    10 double _cyl(double vd, double besarg, double siarg); 
    11 double _cyl(double vd, double besarg, double siarg) 
     2// besarg = q * R * sin(theta) 
     3// siarg = q * L/2 * cos(theta) 
     4static double _cyl(double vd, double besarg, double siarg) 
    125{ 
    136    return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg); 
    147} 
    158 
    16 double form_volume(double radius, double thickness, double length) 
     9static double 
     10form_volume(double radius, double thickness, double length) 
    1711{ 
    18     return M_PI*(radius+thickness)*(radius+thickness)*(length+2.0*thickness); 
     12    return M_PI*square(radius+thickness)*(length+2.0*thickness); 
    1913} 
    2014 
    21 double Iq(double q, 
     15static double 
     16Iq(double q, 
    2217    double core_sld, 
    2318    double shell_sld, 
    … …  
    2823{ 
    2924    // precalculate constants 
    30     const double core_qr = q*radius; 
    31     const double core_qh = q*0.5*length; 
     25    const double core_r = radius; 
     26    const double core_h = 0.5*length; 
    3227    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 
    33     const double shell_qr = q*(radius + thickness); 
    34     const double shell_qh = q*(0.5*length + thickness); 
     28    const double shell_r = (radius + thickness); 
     29    const double shell_h = (0.5*length + thickness); 
    3530    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    3631    double total = 0.0; 
    37     // double lower=0, upper=M_PI_2; 
    3832    for (int i=0; i<76 ;i++) { 
    39         // translate a point in [-1,1] to a point in [lower,upper] 
    40         //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    41         double sn, cn; 
    42         const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
    43         SINCOS(alpha, sn, cn); 
    44         const double fq = _cyl(core_vd, core_qr*sn, core_qh*cn) 
    45             + _cyl(shell_vd, shell_qr*sn, shell_qh*cn); 
    46         total += Gauss76Wt[i] * fq * fq * sn; 
     33        // translate a point in [-1,1] to a point in [0, pi/2] 
     34        //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
     35        double sin_theta, cos_theta; 
     36        const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4; 
     37        SINCOS(theta, sin_theta,  cos_theta); 
     38        const double qab = q*sin_theta; 
     39        const double qc = q*cos_theta; 
     40        const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
     41            + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
     42        total += Gauss76Wt[i] * fq * fq * sin_theta; 
    4743    } 
    4844    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    6460    double q, sin_alpha, cos_alpha; 
    6561    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     62    const double qab = q*sin_alpha; 
     63    const double qc = q*cos_alpha; 
    6664 
    67     const double core_qr = q*radius; 
    68     const double core_qh = q*0.5*length; 
     65    const double core_r = radius; 
     66    const double core_h = 0.5*length; 
    6967    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 
    70     const double shell_qr = q*(radius + thickness); 
    71     const double shell_qh = q*(0.5*length + thickness); 
     68    const double shell_r = (radius + thickness); 
     69    const double shell_h = (0.5*length + thickness); 
    7270    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    7371 
    74     const double fq = _cyl(core_vd, core_qr*sin_alpha, core_qh*cos_alpha) 
    75         + _cyl(shell_vd, shell_qr*sin_alpha, shell_qh*cos_alpha); 
     72    const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
     73        + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
    7674    return 1.0e-4 * fq * fq; 
    7775} 

  • sasmodels/models/core_shell_ellipsoid.c

    r0a3d9b2 r2a0b2b1  
    1 double form_volume(double radius_equat_core, 
    2                    double polar_core, 
    3                    double equat_shell, 
    4                    double polar_shell); 
    5 double Iq(double q, 
    6           double radius_equat_core, 
    7           double x_core, 
    8           double thick_shell, 
    9           double x_polar_shell, 
    10           double core_sld, 
    11           double shell_sld, 
    12           double solvent_sld); 
    131 
     2// Converted from Igor function gfn4, using the same pattern as ellipsoid 
     3// for evaluating the parts of the integral. 
     4//     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN 
     5//                       BY (53) & (58-59) IN CHEN AND 
     6//                       KOTLARCHYK REFERENCE 
     7// 
     8//       <OBLATE ELLIPSOID> 
     9static double 
     10_cs_ellipsoid_kernel(double qab, double qc, 
     11    double equat_core, double polar_core, 
     12    double equat_shell, double polar_shell, 
     13    double sld_core_shell, double sld_shell_solvent) 
     14{ 
     15    const double qr_core = sqrt(square(equat_core*qab) + square(polar_core*qc)); 
     16    const double si_core = sas_3j1x_x(qr_core); 
     17    const double volume_core = M_4PI_3*equat_core*equat_core*polar_core; 
     18    const double fq_core = si_core*volume_core*sld_core_shell; 
    1419 
    15 double Iqxy(double qx, double qy, 
    16           double radius_equat_core, 
    17           double x_core, 
    18           double thick_shell, 
    19           double x_polar_shell, 
    20           double core_sld, 
    21           double shell_sld, 
    22           double solvent_sld, 
    23           double theta, 
    24           double phi); 
     20    const double qr_shell = sqrt(square(equat_shell*qab) + square(polar_shell*qc)); 
     21    const double si_shell = sas_3j1x_x(qr_shell); 
     22    const double volume_shell = M_4PI_3*equat_shell*equat_shell*polar_shell; 
     23    const double fq_shell = si_shell*volume_shell*sld_shell_solvent; 
    2524 
     25    return fq_core + fq_shell; 
     26} 
    2627 
    27 double form_volume(double radius_equat_core, 
    28                    double x_core, 
    29                    double thick_shell, 
    30                    double x_polar_shell) 
     28static double 
     29form_volume(double radius_equat_core, 
     30    double x_core, 
     31    double thick_shell, 
     32    double x_polar_shell) 
    3133{ 
    3234    const double equat_shell = radius_equat_core + thick_shell; 
    … …  
    3739 
    3840static double 
    39 core_shell_ellipsoid_xt_kernel(double q, 
    40           double radius_equat_core, 
    41           double x_core, 
    42           double thick_shell, 
    43           double x_polar_shell, 
    44           double core_sld, 
    45           double shell_sld, 
    46           double solvent_sld) 
     41Iq(double q, 
     42    double radius_equat_core, 
     43    double x_core, 
     44    double thick_shell, 
     45    double x_polar_shell, 
     46    double core_sld, 
     47    double shell_sld, 
     48    double solvent_sld) 
    4749{ 
    48     const double lolim = 0.0; 
    49     const double uplim = 1.0; 
    50  
    51  
    52     const double delpc = core_sld - shell_sld; //core - shell 
    53     const double delps = shell_sld - solvent_sld; //shell - solvent 
    54  
     50    const double sld_core_shell = core_sld - shell_sld; 
     51    const double sld_shell_solvent = shell_sld - solvent_sld; 
    5552 
    5653    const double polar_core = radius_equat_core*x_core; 
    … …  
    5855    const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 
    5956 
    60     double summ = 0.0;   //initialize intergral 
     57    // translate from [-1, 1] => [0, 1] 
     58    const double m = 0.5; 
     59    const double b = 0.5; 
     60    double total = 0.0;     //initialize intergral 
    6161    for(int i=0;i<76;i++) { 
    62         double zi = 0.5*( Gauss76Z[i]*(uplim-lolim) + uplim + lolim ); 
    63         double yyy = gfn4(zi, radius_equat_core, polar_core, equat_shell, 
    64                           polar_shell, delpc, delps, q); 
    65         summ += Gauss76Wt[i] * yyy; 
     62        const double cos_theta = Gauss76Z[i]*m + b; 
     63        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     64        double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta, 
     65            radius_equat_core, polar_core, 
     66            equat_shell, polar_shell, 
     67            sld_core_shell, sld_shell_solvent); 
     68        total += Gauss76Wt[i] * fq * fq; 
    6669    } 
    67     summ *= 0.5*(uplim-lolim); 
     70    total *= m; 
    6871 
    6972    // convert to [cm-1] 
    70     return 1.0e-4 * summ; 
     73    return 1.0e-4 * total; 
    7174} 
    7275 
    7376static double 
    74 core_shell_ellipsoid_xt_kernel_2d(double qx, double qy, 
    75           double radius_equat_core, 
    76           double x_core, 
    77           double thick_shell, 
    78           double x_polar_shell, 
    79           double core_sld, 
    80           double shell_sld, 
    81           double solvent_sld, 
    82           double theta, 
    83           double phi) 
     77Iqxy(double qx, double qy, 
     78    double radius_equat_core, 
     79    double x_core, 
     80    double thick_shell, 
     81    double x_polar_shell, 
     82    double core_sld, 
     83    double shell_sld, 
     84    double solvent_sld, 
     85    double theta, 
     86    double phi) 
    8487{ 
    8588    double q, sin_alpha, cos_alpha; 
    8689    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     90    const double qab = q*sin_alpha; 
     91    const double qc = q*cos_alpha; 
    8792 
    88     const double sldcs = core_sld - shell_sld; 
    89     const double sldss = shell_sld- solvent_sld; 
     93    const double sld_core_shell = core_sld - shell_sld; 
     94    const double sld_shell_solvent = shell_sld - solvent_sld; 
    9095 
    9196    const double polar_core = radius_equat_core*x_core; 
    … …  
    9398    const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 
    9499 
    95     // Call the IGOR library function to get the kernel: 
    96     // MUST use gfn4 not gf2 because of the def of params. 
    97     double answer = gfn4(cos_alpha, 
    98                   radius_equat_core, 
    99                   polar_core, 
    100                   equat_shell, 
    101                   polar_shell, 
    102                   sldcs, 
    103                   sldss, 
    104                   q); 
     100    double fq = _cs_ellipsoid_kernel(qab, qc, 
     101                  radius_equat_core, polar_core, 
     102                  equat_shell, polar_shell, 
     103                  sld_core_shell, sld_shell_solvent); 
    105104 
    106105    //convert to [cm-1] 
    107     answer *= 1.0e-4; 
    108  
    109     return answer; 
     106    return 1.0e-4 * fq * fq; 
    110107} 
    111  
    112 double Iq(double q, 
    113           double radius_equat_core, 
    114           double x_core, 
    115           double thick_shell, 
    116           double x_polar_shell, 
    117           double core_sld, 
    118           double shell_sld, 
    119           double solvent_sld) 
    120 { 
    121     double intensity = core_shell_ellipsoid_xt_kernel(q, 
    122            radius_equat_core, 
    123            x_core, 
    124            thick_shell, 
    125            x_polar_shell, 
    126            core_sld, 
    127            shell_sld, 
    128            solvent_sld); 
    129  
    130     return intensity; 
    131 } 
    132  
    133  
    134 double Iqxy(double qx, double qy, 
    135           double radius_equat_core, 
    136           double x_core, 
    137           double thick_shell, 
    138           double x_polar_shell, 
    139           double core_sld, 
    140           double shell_sld, 
    141           double solvent_sld, 
    142           double theta, 
    143           double phi) 
    144 { 
    145     double intensity = core_shell_ellipsoid_xt_kernel_2d(qx, qy, 
    146                        radius_equat_core, 
    147                        x_core, 
    148                        thick_shell, 
    149                        x_polar_shell, 
    150                        core_sld, 
    151                        shell_sld, 
    152                        solvent_sld, 
    153                        theta, 
    154                        phi); 
    155  
    156     return intensity; 
    157 } 

  • sasmodels/models/core_shell_ellipsoid.py

    r9802ab3 r2a0b2b1  
    2525ellipsoid. This may have some undesirable effects if the aspect ratio of the 
    2626ellipsoid is large (ie, if $X << 1$ or $X >> 1$ ), when the $S(q)$ 
    27 - which assumes spheres - will not in any case be valid.  Generating a  
    28 custom product model will enable separate effective volume fraction and effective  
     27- which assumes spheres - will not in any case be valid.  Generating a 
     28custom product model will enable separate effective volume fraction and effective 
    2929radius in the $S(q)$. 
    3030 
    … …  
    4444 
    4545.. math:: 
    46     \begin{align}     
     46    \begin{align} 
    4747    F(q,\alpha) = &f(q,radius\_equat\_core,radius\_equat\_core.x\_core,\alpha) \\ 
    4848    &+ f(q,radius\_equat\_core + thick\_shell,radius\_equat\_core.x\_core + thick\_shell.x\_polar\_shell,\alpha) 
    49     \end{align}  
     49    \end{align} 
    5050 
    5151where 
    52   
     52 
    5353.. math:: 
    5454 
    … …  
    7777   F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha} 
    7878 
    79 For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,  
     79For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data, 
    8080see the :ref:`elliptical-cylinder` model for further information. 
    8181 
    … …  
    139139# pylint: enable=bad-whitespace, line-too-long 
    140140 
    141 source = ["lib/sas_3j1x_x.c", "lib/gfn.c", "lib/gauss76.c", 
    142           "core_shell_ellipsoid.c"] 
     141source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 
    143142 
    144143def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): 

  • sasmodels/models/core_shell_parallelepiped.c

    r92dfe0c r2a0b2b1  
    1 double form_volume(double length_a, double length_b, double length_c,  
     1double form_volume(double length_a, double length_b, double length_c, 
    22                   double thick_rim_a, double thick_rim_b, double thick_rim_c); 
    33double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld, 
    … …  
    99            double thick_rim_c, double theta, double phi, double psi); 
    1010 
    11 double form_volume(double length_a, double length_b, double length_c,  
     11double form_volume(double length_a, double length_b, double length_c, 
    1212                   double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    1313{ 
    1414    //return length_a * length_b * length_c; 
    15     return length_a * length_b * length_c +  
    16            2.0 * thick_rim_a * length_b * length_c +  
     15    return length_a * length_b * length_c + 
     16           2.0 * thick_rim_a * length_b * length_c + 
    1717           2.0 * thick_rim_b * length_a * length_c + 
    1818           2.0 * thick_rim_c * length_a * length_b; 
    … …  
    3434    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
    3535    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    36      
     36 
    3737    const double mu = 0.5 * q * length_b; 
    38      
     38 
    3939    //calculate volume before rescaling (in original code, but not used) 
    40     //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);         
    41     //double vol = length_a * length_b * length_c +  
    42     //       2.0 * thick_rim_a * length_b * length_c +  
     40    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     41    //double vol = length_a * length_b * length_c + 
     42    //       2.0 * thick_rim_a * length_b * length_c + 
    4343    //       2.0 * thick_rim_b * length_a * length_c + 
    4444    //       2.0 * thick_rim_c * length_a * length_b; 
    45      
     45 
    4646    // Scale sides by B 
    4747    const double a_scaled = length_a / length_b; 
    … …  
    101101            //   ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors 
    102102            // This is the function called by csparallelepiped_analytical_2D_scaled, 
    103             // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c         
    104              
     103            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c 
     104 
    105105            //  correct FF : sum of square of phase factors 
    106106            inner_total += Gauss76Wt[j] * form * form; 
    … …  
    136136    double q, zhat, yhat, xhat; 
    137137    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     138    const double qa = q*xhat; 
     139    const double qb = q*yhat; 
     140    const double qc = q*zhat; 
    138141 
    139142    // cspkernel in csparallelepiped recoded here 
    … …  
    160163    double tc = length_a + 2.0*thick_rim_c; 
    161164    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    162     double siA = sas_sinx_x(0.5*q*length_a*xhat); 
    163     double siB = sas_sinx_x(0.5*q*length_b*yhat); 
    164     double siC = sas_sinx_x(0.5*q*length_c*zhat); 
    165     double siAt = sas_sinx_x(0.5*q*ta*xhat); 
    166     double siBt = sas_sinx_x(0.5*q*tb*yhat); 
    167     double siCt = sas_sinx_x(0.5*q*tc*zhat); 
    168      
     165    double siA = sas_sinx_x(0.5*length_a*qa); 
     166    double siB = sas_sinx_x(0.5*length_b*qb); 
     167    double siC = sas_sinx_x(0.5*length_c*qc); 
     168    double siAt = sas_sinx_x(0.5*ta*qa); 
     169    double siBt = sas_sinx_x(0.5*tb*qb); 
     170    double siCt = sas_sinx_x(0.5*tc*qc); 
     171 
    169172 
    170173    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
    … …  
    174177               + drB*siA*(siBt-siB)*siC*V2 
    175178               + drC*siA*siB*(siCt*siCt-siC)*V3); 
    176     
     179 
    177180    return 1.0e-4 * f * f; 
    178181} 

  • sasmodels/models/cylinder.c

    r592343f rb34fc77  
    1 double form_volume(double radius, double length); 
    2 double fq(double q, double sn, double cn,double radius, double length); 
    3 double orient_avg_1D(double q, double radius, double length); 
    4 double Iq(double q, double sld, double solvent_sld, double radius, double length); 
    5 double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    6     double radius, double length, double theta, double phi); 
    7  
    81#define INVALID(v) (v.radius<0 || v.length<0) 
    92 
    10 double form_volume(double radius, double length) 
     3static double 
     4form_volume(double radius, double length) 
    115{ 
    126    return M_PI*radius*radius*length; 
    137} 
    148 
    15 double fq(double q, double sn, double cn, double radius, double length) 
     9static double 
     10fq(double qab, double qc, double radius, double length) 
    1611{ 
    17     // precompute qr and qh to save time in the loop 
    18     const double qr = q*radius; 
    19     const double qh = q*0.5*length;  
    20     return sas_2J1x_x(qr*sn) * sas_sinx_x(qh*cn); 
     12    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length); 
    2113} 
    2214 
    23 double orient_avg_1D(double q, double radius, double length) 
     15static double 
     16orient_avg_1D(double q, double radius, double length) 
    2417{ 
    2518    // translate a point in [-1,1] to a point in [0, pi/2] 
    2619    const double zm = M_PI_4; 
    27     const double zb = M_PI_4;  
     20    const double zb = M_PI_4; 
    2821 
    2922    double total = 0.0; 
    3023    for (int i=0; i<76 ;i++) { 
    31         const double alpha = Gauss76Z[i]*zm + zb; 
    32         double sn, cn; // slots to hold sincos function output 
    33         // alpha(theta,phi) the projection of the cylinder on the detector plane 
    34         SINCOS(alpha, sn, cn); 
    35         total += Gauss76Wt[i] * square( fq(q, sn, cn, radius, length) ) * sn; 
     24        const double theta = Gauss76Z[i]*zm + zb; 
     25        double sin_theta, cos_theta; // slots to hold sincos function output 
     26        // theta (theta,phi) the projection of the cylinder on the detector plane 
     27        SINCOS(theta , sin_theta, cos_theta); 
     28        const double form = fq(q*sin_theta, q*cos_theta, radius, length); 
     29        total += Gauss76Wt[i] * form * form * sin_theta; 
    3630    } 
    3731    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    3933} 
    4034 
    41 double Iq(double q, 
     35static double 
     36Iq(double q, 
    4237    double sld, 
    4338    double solvent_sld, 
    … …  
    4944} 
    5045 
    51  
    52 double Iqxy(double qx, double qy, 
     46static double 
     47Iqxy(double qx, double qy, 
    5348    double sld, 
    5449    double solvent_sld, 
    … …  
    6055    double q, sin_alpha, cos_alpha; 
    6156    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    62     //printf("sn: %g cn: %g\n", sin_alpha, cos_alpha); 
     57    const double qab = q*sin_alpha; 
     58    const double qc = q*cos_alpha; 
     59 
    6360    const double s = (sld-solvent_sld) * form_volume(radius, length); 
    64     const double form = fq(q, sin_alpha, cos_alpha, radius, length); 
     61    const double form = fq(qab, qc, radius, length); 
    6562    return 1.0e-4 * square(s * form); 
    6663} 

  • sasmodels/models/ellipsoid.c

    r3b571ae r2a0b2b1  
    1 double form_volume(double radius_polar, double radius_equatorial); 
    2 double Iq(double q, double sld, double sld_solvent, double radius_polar, double radius_equatorial); 
    3 double Iqxy(double qx, double qy, double sld, double sld_solvent, 
    4     double radius_polar, double radius_equatorial, double theta, double phi); 
    5  
    6 double form_volume(double radius_polar, double radius_equatorial) 
     1static double 
     2form_volume(double radius_polar, double radius_equatorial) 
    73{ 
    84    return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial; 
    95} 
    106 
    11 double Iq(double q, 
     7static  double 
     8Iq(double q, 
    129    double sld, 
    1310    double sld_solvent, 
    … …  
    4138} 
    4239 
    43 double Iqxy(double qx, double qy, 
     40static double 
     41Iqxy(double qx, double qy, 
    4442    double sld, 
    4543    double sld_solvent, 
    … …  
    5149    double q, sin_alpha, cos_alpha; 
    5250    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    53     const double r = sqrt(square(radius_equatorial*sin_alpha) 
    54                           + square(radius_polar*cos_alpha)); 
    55     const double f = sas_3j1x_x(q*r); 
     51    const double qab = q*sin_alpha; 
     52    const double qc = q*cos_alpha; 
     53 
     54    const double qr = sqrt(square(radius_equatorial*qab) + square(radius_polar*qc)); 
     55    const double f = sas_3j1x_x(qr); 
    5656    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    5757 
    5858    return 1.0e-4 * square(f * s); 
    5959} 
    60  

  • sasmodels/models/elliptical_cylinder.c

    r61104c8 r2a0b2b1  
    1 double form_volume(double radius_minor, double r_ratio, double length); 
    2 double Iq(double q, double radius_minor, double r_ratio, double length, 
    3           double sld, double solvent_sld); 
    4 double Iqxy(double qx, double qy, double radius_minor, double r_ratio, double length, 
    5             double sld, double solvent_sld, double theta, double phi, double psi); 
    6  
    7  
    8 double 
     1static double 
    92form_volume(double radius_minor, double r_ratio, double length) 
    103{ 
    … …  
    125} 
    136 
    14 double 
     7static double 
    158Iq(double q, double radius_minor, double r_ratio, double length, 
    169   double sld, double solvent_sld) 
    … …  
    6154 
    6255 
    63 double 
     56static double 
    6457Iqxy(double qx, double qy, 
    6558     double radius_minor, double r_ratio, double length, 
    … …  
    6962    double q, xhat, yhat, zhat; 
    7063    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     64    const double qa = q*xhat; 
     65    const double qb = q*yhat; 
     66    const double qc = q*zhat; 
    7167 
    7268    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 
    7369    // Given:    radius_major = r_ratio * radius_minor 
    74     const double r = radius_minor*sqrt(square(r_ratio*xhat) + square(yhat)); 
    75     const double be = sas_2J1x_x(q*r); 
    76     const double si = sas_sinx_x(q*zhat*0.5*length); 
     70    const double qr = radius_minor*sqrt(square(r_ratio*qa) + square(qb)); 
     71    const double be = sas_2J1x_x(qr); 
     72    const double si = sas_sinx_x(qc*0.5*length); 
    7773    const double Aq = be * si; 
    7874    const double delrho = sld - solvent_sld; 

  • sasmodels/models/fcc_paracrystal.c

    r50beefe r7e0b281  
    1 double form_volume(double radius); 
    2 double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
    3 double Iqxy(double qx, double qy, double dnn, 
    4     double d_factor, double radius,double sld, double solvent_sld, 
    5     double theta, double phi, double psi); 
     1static double 
     2fcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4#if 0  // Equations as written in Matsuoka 
     5    const double a1 = ( qa + qb)/2.0; 
     6    const double a2 = (-qa + qc)/2.0; 
     7    const double a3 = (-qa + qb)/2.0; 
     8#else 
     9    const double a1 = ( qa + qb)/2.0; 
     10    const double a2 = ( qa + qc)/2.0; 
     11    const double a3 = ( qb + qc)/2.0; 
     12#endif 
    613 
    7 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
    8 double _FCCeval(double Theta, double Phi, double temp1, double temp3); 
     14    // Numerator: (1 - exp(a)^2)^3 
     15    //         => (-(exp(2a) - 1))^3 
     16    //         => -expm1(2a)^3 
     17    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a)) 
     18    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
     19    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1) 
     20    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     21    const double exp_arg = exp(arg); 
     22    const double Zq = -cube(expm1(2.0*arg)) 
     23        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
     24          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
     25          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
     26 
     27    return Zq; 
     28} 
    929 
    1030 
    11 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { 
    12  
    13         const double Da = d_factor*dnn; 
    14         const double temp1 = q*q*Da*Da; 
    15         const double temp3 = q*dnn; 
    16  
    17         double retVal = _FCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); 
    18         return(retVal); 
     31// occupied volume fraction calculated from lattice symmetry and sphere radius 
     32static double 
     33fcc_volume_fraction(double radius, double dnn) 
     34{ 
     35    return 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 
    1936} 
    2037 
    21 double _FCCeval(double Theta, double Phi, double temp1, double temp3) { 
    22  
    23         double result; 
    24         double sin_theta,cos_theta,sin_phi,cos_phi; 
    25         SINCOS(Theta, sin_theta, cos_theta); 
    26         SINCOS(Phi, sin_phi, cos_phi); 
    27  
    28         const double temp6 =  sin_theta; 
    29         const double temp7 =  sin_theta*sin_phi + cos_theta; 
    30         const double temp8 = -sin_theta*cos_phi + cos_theta; 
    31         const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi; 
    32  
    33         const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 
    34         result = cube(1.0-(temp10*temp10))*temp6 
    35             / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    36               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
    37               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10)); 
    38  
    39         return (result); 
    40 } 
    41  
    42 double form_volume(double radius){ 
     38static double 
     39form_volume(double radius) 
     40{ 
    4341    return sphere_volume(radius); 
    4442} 
    4543 
    4644 
    47 double Iq(double q, double dnn, 
     45static double Iq(double q, double dnn, 
    4846  double d_factor, double radius, 
    49   double sld, double solvent_sld){ 
     47  double sld, double solvent_sld) 
     48{ 
     49    // translate a point in [-1,1] to a point in [0, 2 pi] 
     50    const double phi_m = M_PI; 
     51    const double phi_b = M_PI; 
     52    // translate a point in [-1,1] to a point in [0, pi] 
     53    const double theta_m = M_PI_2; 
     54    const double theta_b = M_PI_2; 
    5055 
    51         //Volume fraction calculated from lattice symmetry and sphere radius 
    52         const double s1 = dnn*sqrt(2.0); 
    53         const double latticescale = 4.0*sphere_volume(radius/s1); 
     56    double outer_sum = 0.0; 
     57    for(int i=0; i<150; i++) { 
     58        double inner_sum = 0.0; 
     59        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     60        double sin_theta, cos_theta; 
     61        SINCOS(theta, sin_theta, cos_theta); 
     62        const double qc = q*cos_theta; 
     63        const double qab = q*sin_theta; 
     64        for(int j=0;j<150;j++) { 
     65            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     66            double sin_phi, cos_phi; 
     67            SINCOS(phi, sin_phi, cos_phi); 
     68            const double qa = qab*cos_phi; 
     69            const double qb = qab*sin_phi; 
     70            const double form = fcc_Zq(qa, qb, qc, dnn, d_factor); 
     71            inner_sum += Gauss150Wt[j] * form; 
     72        } 
     73        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     74        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     75    } 
     76    outer_sum *= theta_m; 
     77    const double Zq = outer_sum/(4.0*M_PI); 
     78    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    5479 
    55     const double va = 0.0; 
    56     const double vb = 2.0*M_PI; 
    57     const double vaj = 0.0; 
    58     const double vbj = M_PI; 
    59  
    60     double summ = 0.0; 
    61     double answer = 0.0; 
    62         for(int i=0; i<150; i++) { 
    63                 //setup inner integral over the ellipsoidal cross-section 
    64                 double summj=0.0; 
    65                 const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;             //the outer dummy is phi 
    66                 for(int j=0;j<150;j++) { 
    67                         //20 gauss points for the inner integral 
    68                         double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;             //the inner dummy is theta 
    69                         double yyy = Gauss150Wt[j] * _FCC_Integrand(q,dnn,d_factor,ztheta,zphi); 
    70                         summj += yyy; 
    71                 } 
    72                 //now calculate the value of the inner integral 
    73                 double answer = (vbj-vaj)/2.0*summj; 
    74  
    75                 //now calculate outer integral 
    76                 summ = summ+(Gauss150Wt[i] * answer); 
    77         }               //final scaling is done at the end of the function, after the NT_FP64 case 
    78  
    79         answer = (vb-va)/2.0*summ; 
    80         answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale; 
    81  
    82     return answer; 
     80    return fcc_volume_fraction(radius, dnn) * Pq * Zq; 
     81} 
    8382 
    8483 
    85 } 
    86  
    87 double Iqxy(double qx, double qy, 
     84static double Iqxy(double qx, double qy, 
    8885    double dnn, double d_factor, double radius, 
    8986    double sld, double solvent_sld, 
    … …  
    9289    double q, zhat, yhat, xhat; 
    9390    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     91    const double qa = q*xhat; 
     92    const double qb = q*yhat; 
     93    const double qc = q*zhat; 
    9494 
    95     const double a1 = yhat + xhat; 
    96     const double a2 = xhat + zhat; 
    97     const double a3 = yhat + zhat; 
    98     const double qd = 0.5*q*dnn; 
    99     const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
    100     const double tanh_qd = tanh(arg); 
    101     const double cosh_qd = cosh(arg); 
    102     const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd) 
    103                     * tanh_qd/(1. - cos(qd*a2)/cosh_qd) 
    104                     * tanh_qd/(1. - cos(qd*a3)/cosh_qd); 
    105  
    106     //if (isnan(Zq)) printf("q:(%g,%g) qd: %g a1: %g a2: %g a3: %g arg: %g\n", qx, qy, qd, a1, a2, a3, arg); 
    107  
    108     const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; 
    109     //the occupied volume of the lattice 
    110     const double lattice_scale = 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 
    111     return lattice_scale * Fq; 
     95    q = sqrt(qa*qa + qb*qb + qc*qc); 
     96    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     97    const double Zq = fcc_Zq(qa, qb, qc, dnn, d_factor); 
     98    return fcc_volume_fraction(radius, dnn) * Pq * Zq; 
    11299} 

  • sasmodels/models/hollow_cylinder.c

    r592343f r2a0b2b1  
    11double form_volume(double radius, double thickness, double length); 
    22double Iq(double q, double radius, double thickness, double length, double sld, 
    3         double solvent_sld); 
     3    double solvent_sld); 
    44double Iqxy(double qx, double qy, double radius, double thickness, double length, double sld, 
    5         double solvent_sld, double theta, double phi); 
     5    double solvent_sld, double theta, double phi); 
    66 
    77//#define INVALID(v) (v.radius_core >= v.radius) 
    … …  
    1414} 
    1515 
    16  
    1716static double 
    18 _hollow_cylinder_kernel(double q, 
    19     double radius, double thickness, double length, double sin_val, double cos_val) 
     17_fq(double qab, double qc, 
     18    double radius, double thickness, double length) 
    2019{ 
    21     const double qs = q*sin_val; 
    22     const double lam1 = sas_2J1x_x((radius+thickness)*qs); 
    23     const double lam2 = sas_2J1x_x(radius*qs); 
     20    const double lam1 = sas_2J1x_x((radius+thickness)*qab); 
     21    const double lam2 = sas_2J1x_x(radius*qab); 
    2422    const double gamma_sq = square(radius/(radius+thickness)); 
    25     //Note: lim_{thickness -> 0} psi = sas_J0(radius*qs) 
    26     //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qs) 
    27     const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); //SRK 10/19/00 
    28     const double t2 = sas_sinx_x(0.5*q*length*cos_val); 
     23    //Note: lim_{thickness -> 0} psi = sas_J0(radius*qab) 
     24    //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qab) 
     25    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq);    //SRK 10/19/00 
     26    const double t2 = sas_sinx_x(0.5*length*qc); 
    2927    return psi*t2; 
    3028} 
    … …  
    4341{ 
    4442    const double lower = 0.0; 
    45     const double upper = 1.0;           //limits of numerical integral 
     43    const double upper = 1.0;        //limits of numerical integral 
    4644 
    47     double summ = 0.0;                  //initialize intergral 
     45    double summ = 0.0;            //initialize intergral 
    4846    for (int i=0;i<76;i++) { 
    49         const double cos_val = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
    50         const double sin_val = sqrt(1.0 - cos_val*cos_val); 
    51         const double inter = _hollow_cylinder_kernel(q, radius, thickness, length, 
    52                                                      sin_val, cos_val); 
    53         summ += Gauss76Wt[i] * inter * inter; 
     47        const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
     48        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     49        const double form = _fq(q*sin_theta, q*cos_theta, 
     50                                radius, thickness, length); 
     51        summ += Gauss76Wt[i] * form * form; 
    5452    } 
    5553 
    … …  
    6664    double q, sin_alpha, cos_alpha; 
    6765    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    68     const double Aq = _hollow_cylinder_kernel(q, radius, thickness, length, 
    69         sin_alpha, cos_alpha); 
     66    const double qab = q*sin_alpha; 
     67    const double qc = q*cos_alpha; 
     68 
     69    const double form = _fq(qab, qc, radius, thickness, length); 
    7070 
    7171    const double vol = form_volume(radius, thickness, length); 
    72     return _hollow_cylinder_scaling(Aq*Aq, solvent_sld-sld, vol); 
     72    return _hollow_cylinder_scaling(form*form, solvent_sld-sld, vol); 
    7373} 
    74  

  • sasmodels/models/parallelepiped.c

    rd605080 r2a0b2b1  
    2020{ 
    2121    const double mu = 0.5 * q * length_b; 
    22      
     22 
    2323    // Scale sides by B 
    2424    const double a_scaled = length_a / length_b; 
    2525    const double c_scaled = length_c / length_b; 
    26          
     26 
    2727    // outer integral (with gauss points), integration limits = 0, 1 
    2828    double outer_total = 0; //initialize integral 
    … …  
    6969    double q, xhat, yhat, zhat; 
    7070    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     71    const double qa = q*xhat; 
     72    const double qb = q*yhat; 
     73    const double qc = q*zhat; 
    7174 
    72     const double siA = sas_sinx_x(0.5*length_a*q*xhat); 
    73     const double siB = sas_sinx_x(0.5*length_b*q*yhat); 
    74     const double siC = sas_sinx_x(0.5*length_c*q*zhat); 
     75    const double siA = sas_sinx_x(0.5*length_a*qa); 
     76    const double siB = sas_sinx_x(0.5*length_b*qb); 
     77    const double siC = sas_sinx_x(0.5*length_c*qc); 
    7578    const double V = form_volume(length_a, length_b, length_c); 
    7679    const double drho = (sld - solvent_sld); 

  • sasmodels/models/sc_paracrystal.c

    r50beefe r7e0b281  
    1 double form_volume(double radius); 
     1static double 
     2sc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4    const double a1 = qa; 
     5    const double a2 = qb; 
     6    const double a3 = qc; 
    27 
    3 double Iq(double q, 
    4           double dnn, 
    5           double d_factor, 
    6           double radius, 
    7           double sphere_sld, 
    8           double solvent_sld); 
     8    // Numerator: (1 - exp(a)^2)^3 
     9    //         => (-(exp(2a) - 1))^3 
     10    //         => -expm1(2a)^3 
     11    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a)) 
     12    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
     13    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1) 
     14    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     15    const double exp_arg = exp(arg); 
     16    const double Zq = -cube(expm1(2.0*arg)) 
     17        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
     18          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
     19          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
    920 
    10 double Iqxy(double qx, double qy, 
    11             double dnn, 
    12             double d_factor, 
    13             double radius, 
    14             double sphere_sld, 
    15             double solvent_sld, 
    16             double theta, 
    17             double phi, 
    18             double psi); 
     21    return Zq; 
     22} 
    1923 
    20 double form_volume(double radius) 
     24// occupied volume fraction calculated from lattice symmetry and sphere radius 
     25static double 
     26sc_volume_fraction(double radius, double dnn) 
     27{ 
     28    return sphere_volume(radius/dnn); 
     29} 
     30 
     31static double 
     32form_volume(double radius) 
    2133{ 
    2234    return sphere_volume(radius); 
    2335} 
    2436 
    25 static double 
    26 sc_eval(double theta, double phi, double temp3, double temp4, double temp5) 
     37 
     38static double Iq(double q, double dnn, 
     39  double d_factor, double radius, 
     40  double sld, double solvent_sld) 
    2741{ 
    28     double cnt, snt; 
    29     SINCOS(theta, cnt, snt); 
     42    // translate a point in [-1,1] to a point in [0, 2 pi] 
     43    const double phi_m = M_PI_4; 
     44    const double phi_b = M_PI_4; 
     45    // translate a point in [-1,1] to a point in [0, pi] 
     46    const double theta_m = M_PI_4; 
     47    const double theta_b = M_PI_4; 
    3048 
    31     double cnp, snp; 
    32     SINCOS(phi, cnp, snp); 
    3349 
    34         double temp6 = snt; 
    35         double temp7 = -1.0*temp3*snt*cnp; 
    36         double temp8 = temp3*snt*snp; 
    37         double temp9 = temp3*cnt; 
    38         double result = temp6/((1.0-temp4*cos((temp7))+temp5)* 
    39                                (1.0-temp4*cos((temp8))+temp5)* 
    40                                (1.0-temp4*cos((temp9))+temp5)); 
    41         return (result); 
     50    double outer_sum = 0.0; 
     51    for(int i=0; i<150; i++) { 
     52        double inner_sum = 0.0; 
     53        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     54        double sin_theta, cos_theta; 
     55        SINCOS(theta, sin_theta, cos_theta); 
     56        const double qc = q*cos_theta; 
     57        const double qab = q*sin_theta; 
     58        for(int j=0;j<150;j++) { 
     59            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     60            double sin_phi, cos_phi; 
     61            SINCOS(phi, sin_phi, cos_phi); 
     62            const double qa = qab*cos_phi; 
     63            const double qb = qab*sin_phi; 
     64            const double form = sc_Zq(qa, qb, qc, dnn, d_factor); 
     65            inner_sum += Gauss150Wt[j] * form; 
     66        } 
     67        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     68        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     69    } 
     70    outer_sum *= theta_m; 
     71    const double Zq = outer_sum/M_PI_2; 
     72    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     73 
     74    return sc_volume_fraction(radius, dnn) * Pq * Zq; 
    4275} 
    4376 
    44 static double 
    45 sc_integrand(double dnn, double d_factor, double qq, double xx, double yy) 
    46 { 
    47     //Function to calculate integrand values for simple cubic structure 
    4877 
    49         double da = d_factor*dnn; 
    50         double temp1 = qq*qq*da*da; 
    51         double temp2 = cube(-expm1(-temp1)); 
    52         double temp3 = qq*dnn; 
    53         double temp4 = 2.0*exp(-0.5*temp1); 
    54         double temp5 = exp(-1.0*temp1); 
    55  
    56         double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5)/M_PI_2; 
    57  
    58         return(integrand); 
    59 } 
    60  
    61 double Iq(double q, 
    62           double dnn, 
    63           double d_factor, 
    64           double radius, 
    65           double sphere_sld, 
    66           double solvent_sld) 
    67 { 
    68         const double va = 0.0; 
    69         const double vb = M_PI_2; //orientation average, outer integral 
    70  
    71     double summ=0.0; 
    72     double answer=0.0; 
    73  
    74         for(int i=0;i<150;i++) { 
    75                 //setup inner integral over the ellipsoidal cross-section 
    76                 double summj=0.0; 
    77                 double zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0; 
    78                 for(int j=0;j<150;j++) { 
    79                         //150 gauss points for the inner integral 
    80                         double zij = ( Gauss150Z[j]*(vb-va) + va + vb )/2.0; 
    81                         double tmp = Gauss150Wt[j] * sc_integrand(dnn,d_factor,q,zi,zij); 
    82                         summj += tmp; 
    83                 } 
    84                 //now calculate the value of the inner integral 
    85                 answer = (vb-va)/2.0*summj; 
    86  
    87                 //now calculate outer integral 
    88                 double tmp = Gauss150Wt[i] * answer; 
    89                 summ += tmp; 
    90         }               //final scaling is done at the end of the function, after the NT_FP64 case 
    91  
    92         answer = (vb-va)/2.0*summ; 
    93  
    94         //Volume fraction calculated from lattice symmetry and sphere radius 
    95         // NB: 4/3 pi r^3 / dnn^3 = 4/3 pi(r/dnn)^3 
    96         const double latticeScale = sphere_volume(radius/dnn); 
    97  
    98         answer *= sphere_form(q, radius, sphere_sld, solvent_sld)*latticeScale; 
    99  
    100         return answer; 
    101 } 
    102  
    103 double Iqxy(double qx, double qy, 
    104           double dnn, 
    105           double d_factor, 
    106           double radius, 
    107           double sphere_sld, 
    108           double solvent_sld, 
    109           double theta, 
    110           double phi, 
    111           double psi) 
     78static double Iqxy(double qx, double qy, 
     79    double dnn, double d_factor, double radius, 
     80    double sld, double solvent_sld, 
     81    double theta, double phi, double psi) 
    11282{ 
    11383    double q, zhat, yhat, xhat; 
    11484    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     85    const double qa = q*xhat; 
     86    const double qb = q*yhat; 
     87    const double qc = q*zhat; 
    11588 
    116     const double qd = q*dnn; 
    117     const double arg = 0.5*square(qd*d_factor); 
    118     const double tanh_qd = tanh(arg); 
    119     const double cosh_qd = cosh(arg); 
    120     const double Zq = tanh_qd/(1. - cos(qd*zhat)/cosh_qd) 
    121                     * tanh_qd/(1. - cos(qd*yhat)/cosh_qd) 
    122                     * tanh_qd/(1. - cos(qd*xhat)/cosh_qd); 
    123  
    124     const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; 
    125     //the occupied volume of the lattice 
    126     const double lattice_scale = sphere_volume(radius/dnn); 
    127     return lattice_scale * Fq; 
     89    q = sqrt(qa*qa + qb*qb + qc*qc); 
     90    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     91    const double Zq = sc_Zq(qa, qb, qc, dnn, d_factor); 
     92    return sc_volume_fraction(radius, dnn) * Pq * Zq; 
    12893} 

  • sasmodels/models/sc_paracrystal.py

    r69e1afc r2a0b2b1  
    152152    [{}, 0.001, 10.3048], 
    153153    [{}, 0.215268, 0.00814889], 
    154     [{}, (0.414467), 0.001313289], 
     154    [{}, 0.414467, 0.001313289], 
    155155    [{'theta':10.0,'phi':20,'psi':30.0},(0.045,-0.035),18.0397138402 ], 
    156156    [{'theta':10.0,'phi':20,'psi':30.0},(0.023,0.045),0.0177333171285 ] 
    157157    ] 
    158  
    159  

  • sasmodels/models/stacked_disks.c

    r19f996b rb34fc77  
    11static double stacked_disks_kernel( 
    2     double q, 
     2    double qab, 
     3    double qc, 
    34    double halfheight, 
    45    double thick_layer, 
    … …  
    910    double layer_sld, 
    1011    double solvent_sld, 
    11     double sin_alpha, 
    12     double cos_alpha, 
    1312    double d) 
    1413 
    … …  
    2019    // zi is the dummy variable for the integration (x in Feigin's notation) 
    2120 
    22     const double besarg1 = q*radius*sin_alpha; 
    23     //const double besarg2 = q*radius*sin_alpha; 
     21    const double besarg1 = radius*qab; 
     22    //const double besarg2 = radius*qab; 
    2423 
    25     const double sinarg1 = q*halfheight*cos_alpha; 
    26     const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha; 
     24    const double sinarg1 = halfheight*qc; 
     25    const double sinarg2 = (halfheight+thick_layer)*qc; 
    2726 
    2827    const double be1 = sas_2J1x_x(besarg1); 
    … …  
    4342 
    4443    // loop for the structure factor S(q) 
    45     double qd_cos_alpha = q*d*cos_alpha; 
     44    double qd_cos_alpha = d*qc; 
    4645    //d*cos_alpha is the projection of d onto q (in other words the component 
    4746    //of d that is parallel to q. 
    … …  
    8483        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8584        SINCOS(zi, sin_alpha, cos_alpha); 
    86         double yyy = stacked_disks_kernel(q, 
     85        double yyy = stacked_disks_kernel(q*sin_alpha, q*cos_alpha, 
    8786                           halfheight, 
    8887                           thick_layer, 
    … …  
    9392                           layer_sld, 
    9493                           solvent_sld, 
    95                            sin_alpha, 
    96                            cos_alpha, 
    9794                           d); 
    9895        summ += Gauss76Wt[i] * yyy * sin_alpha; 
    … …  
    152149    double phi) 
    153150{ 
    154     int n_stacking = (int)(fp_n_stacking + 0.5); 
    155151    double q, sin_alpha, cos_alpha; 
    156152    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     153    const double qab = q*sin_alpha; 
     154    const double qc = q*cos_alpha; 
    157155 
     156    int n_stacking = (int)(fp_n_stacking + 0.5); 
    158157    double d = 2.0 * thick_layer + thick_core; 
    159158    double halfheight = 0.5*thick_core; 
    160     double answer = stacked_disks_kernel(q, 
     159    double answer = stacked_disks_kernel(qab, qc, 
    161160                     halfheight, 
    162161                     thick_layer, 
    … …  
    167166                     layer_sld, 
    168167                     solvent_sld, 
    169                      sin_alpha, 
    170                      cos_alpha, 
    171168                     d); 
    172169 

  • sasmodels/models/triaxial_ellipsoid.c

    r68dd6a9 r2a0b2b1  
    1 double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar); 
    2 double Iq(double q, double sld, double sld_solvent, 
    3     double radius_equat_minor, double radius_equat_major, double radius_polar); 
    4 double Iqxy(double qx, double qy, double sld, double sld_solvent, 
    5     double radius_equat_minor, double radius_equat_major, double radius_polar, double theta, double phi, double psi); 
    6  
    71//#define INVALID(v) (v.radius_equat_minor > v.radius_equat_major || v.radius_equat_major > v.radius_polar) 
    82 
    9  
    10 double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
     3static double 
     4form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
    115{ 
    126    return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar; 
    137} 
    148 
    15 double Iq(double q, 
     9static double 
     10Iq(double q, 
    1611    double sld, 
    1712    double sld_solvent, 
    … …  
    4540    // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi 
    4641    const double fqsq = outer/4.0;  // = outer*um*zm*8.0/(4.0*M_PI); 
    47     const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
    48     return 1.0e-4 * s * s * fqsq; 
     42    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
     43    const double drho = (sld - sld_solvent); 
     44    return 1.0e-4 * square(vol*drho) * fqsq; 
    4945} 
    5046 
    51 double Iqxy(double qx, double qy, 
     47static double 
     48Iqxy(double qx, double qy, 
    5249    double sld, 
    5350    double sld_solvent, 
    … …  
    6158    double q, xhat, yhat, zhat; 
    6259    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     60    const double qa = q*xhat; 
     61    const double qb = q*yhat; 
     62    const double qc = q*zhat; 
    6363 
    64     const double r = sqrt(square(radius_equat_minor*xhat) 
    65                           + square(radius_equat_major*yhat) 
    66                           + square(radius_polar*zhat)); 
    67     const double fq = sas_3j1x_x(q*r); 
    68     const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
     64    const double qr = sqrt(square(radius_equat_minor*qa) 
     65                           + square(radius_equat_major*qb) 
     66                           + square(radius_polar*qc)); 
     67    const double fq = sas_3j1x_x(qr); 
     68    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
     69    const double drho = (sld - sld_solvent); 
    6970 
    70     return 1.0e-4 * square(s * fq); 
     71    return 1.0e-4 * square(vol * drho * fq); 
    7172} 
    72  

  • sasmodels/weights.py

    r41e7f2e rd4c33d6  
    5555        """ 
    5656        sigma = self.width * center if relative else self.width 
     57        if not relative: 
     58            # For orientation, the jitter is relative to 0 not the angle 
     59            #center = 0 
     60            pass 
    5761        if sigma == 0 or self.npts < 2: 
    5862            if lb <= center <= ub: 
Note: See TracChangeset for help on using the changeset viewer.