Changes in / [a0d75ce:ff1fff5] in sasmodels


Ignore:
Files:
56 edited

Legend:

Unmodified
Added
Removed

  • README.rst

    r84bc3c1 re30d645  
    4646lamellar.py is an example of a single file model with embedded C code. 
    4747 
    48 Magnetism hasn't been implemented yet.  We may want a separate Imagnetic 
    49 calculator with the extra parameters and calculations.   We should 
    50 return all desired spin states together so we can share the work of 
    51 computing the form factors for the different magnetic contrasts.  This 
    52 will mean extending the data handler to support multiple cross sections 
    53 in the same data set. 
    54  
    5548|TravisStatus|_ 
    5649 

  • sasmodels/compare.py

    ra0d75ce r6831fa0  
    3333import datetime 
    3434import traceback 
     35import re 
    3536 
    3637import numpy as np  # type: ignore 
    … …  
    3839from . import core 
    3940from . import kerneldll 
    40 from . import weights 
     41from . import exception 
    4142from .data import plot_theory, empty_data1D, empty_data2D 
    4243from .direct_model import DirectModel 
    4344from .convert import revert_name, revert_pars, constrain_new_to_old 
     45from .generate import FLOAT_RE 
    4446 
    4547try: 
    4648    from typing import Optional, Dict, Any, Callable, Tuple 
    47 except: 
     49except Exception: 
    4850    pass 
    4951else: 
    … …  
    8183    -default/-demo* use demo vs default parameters 
    8284    -html shows the model docs instead of running the model 
    83     -title="graph title" adds a title to the plot 
     85    -title="note" adds note to the plot title, after the model name 
    8486 
    8587Any two calculation engines can be selected for comparison: 
    … …  
    466468            if k.endswith('.type'): 
    467469                par = k[:-5] 
     470                if v == 'gaussian': continue 
    468471                cls = dispersers[v if v != 'rectangle' else 'rectangula'] 
    469472                handle = cls() 
    470473                model[0].disperser_handles[par] = handle 
    471                 model[0].set_dispersion(par, handle) 
     474                try: 
     475                    model[0].set_dispersion(par, handle) 
     476                except Exception: 
     477                    exception.annotate_exception("while setting %s to %r" 
     478                                                 %(par, v)) 
     479                    raise 
     480 
    472481 
    473482        #print("sasview pars",oldpars) 
    … …  
    605614    parameters. 
    606615    """ 
    607     n_base, n_comp = opts['n1'], opts['n2'] 
    608     pars = opts['pars'] 
     616    n_base, n_comp = opts['count'] 
     617    pars, pars2 = opts['pars'] 
    609618    data = opts['data'] 
    610619 
    … …  
    630639    if n_comp > 0: 
    631640        try: 
    632             comp_raw, comp_time = time_calculation(comp, pars, n_comp) 
     641            comp_raw, comp_time = time_calculation(comp, pars2, n_comp) 
    633642            comp_value = np.ma.masked_invalid(comp_raw) 
    634643            print("%s t=%.2f ms, intensity=%.0f" 
    … …  
    680689        else: 
    681690            err, errstr, errview = abs(relerr), "rel err", "log" 
     691        #sorted = np.sort(err.flatten()) 
     692        #cutoff = sorted[int(sorted.size*0.95)] 
     693        #err[err>cutoff] = cutoff 
    682694        #err,errstr = base/comp,"ratio" 
    683695        plot_theory(data, None, resid=err, view=errview, use_data=False) 
    … …  
    691703    fig = plt.gcf() 
    692704    extra_title = ' '+opts['title'] if opts['title'] else '' 
    693     fig.suptitle(opts['name'] + extra_title) 
    694  
    695     if n_comp > 0 and n_base > 0 and '-hist' in opts: 
     705    fig.suptitle(":".join(opts['name']) + extra_title) 
     706 
     707    if n_comp > 0 and n_base > 0 and opts['show_hist']: 
    696708        plt.figure() 
    697709        v = relerr 
    … …  
    791803    return pars 
    792804 
     805INTEGER_RE = re.compile("^[+-]?[1-9][0-9]*$") 
     806def isnumber(str): 
     807    match = FLOAT_RE.match(str) 
     808    isfloat = (match and not str[match.end():]) 
     809    return isfloat or INTEGER_RE.match(str) 
    793810 
    794811def parse_opts(argv): 
    … …  
    813830        print("expected parameters: model N1 N2") 
    814831 
    815     name = positional_args[0] 
    816     try: 
    817         model_info = core.load_model_info(name) 
    818     except ImportError as exc: 
    819         print(str(exc)) 
    820         print("Could not find model; use one of:\n    " + models) 
    821         return None 
    822  
    823832    invalid = [o[1:] for o in flags 
    824833               if o[1:] not in NAME_OPTIONS 
    … …  
    828837        return None 
    829838 
     839    name = positional_args[0] 
     840    n1 = int(positional_args[1]) if len(positional_args) > 1 else 1 
     841    n2 = int(positional_args[2]) if len(positional_args) > 2 else 1 
    830842 
    831843    # pylint: disable=bad-whitespace 
    … …  
    895907    # pylint: enable=bad-whitespace 
    896908 
     909    if ':' in name: 
     910        name, name2 = name.split(':',2) 
     911    else: 
     912        name2 = name 
     913    try: 
     914        model_info = core.load_model_info(name) 
     915        model_info2 = core.load_model_info(name2) if name2 != name else model_info 
     916    except ImportError as exc: 
     917        print(str(exc)) 
     918        print("Could not find model; use one of:\n    " + models) 
     919        return None 
     920 
     921    # Get demo parameters from model definition, or use default parameters 
     922    # if model does not define demo parameters 
     923    pars = get_pars(model_info, opts['use_demo']) 
     924    pars2 = get_pars(model_info2, opts['use_demo']) 
     925    # randomize parameters 
     926    #pars.update(set_pars)  # set value before random to control range 
     927    if opts['seed'] > -1: 
     928        pars = randomize_pars(model_info, pars, seed=opts['seed']) 
     929        if model_info != model_info2: 
     930            pars2 = randomize_pars(model_info2, pars2, seed=opts['seed']) 
     931        else: 
     932            pars2 = pars.copy() 
     933        print("Randomize using -random=%i"%opts['seed']) 
     934    if opts['mono']: 
     935        pars = suppress_pd(pars) 
     936        pars2 = suppress_pd(pars2) 
     937 
     938    # Fill in parameters given on the command line 
     939    presets = {} 
     940    presets2 = {} 
     941    for arg in values: 
     942        k, v = arg.split('=', 1) 
     943        if k not in pars and k not in pars2: 
     944            # extract base name without polydispersity info 
     945            s = set(p.split('_pd')[0] for p in pars) 
     946            print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s)))) 
     947            return None 
     948        v1, v2 = v.split(':',2) if ':' in v else (v,v) 
     949        if v1 and k in pars: 
     950            presets[k] = float(v1) if isnumber(v1) else v1 
     951        if v2 and k in pars2: 
     952            presets2[k] = float(v2) if isnumber(v2) else v2 
     953 
     954    # Evaluate preset parameter expressions 
     955    context = pars.copy() 
     956    context.update((k,v) for k,v in presets.items() if isinstance(v, float)) 
     957    for k, v in presets.items(): 
     958        if not isinstance(v, float) and not k.endswith('_type'): 
     959            presets[k] = eval(v, context) 
     960    context = pars.copy() 
     961    context.update(presets) 
     962    context.update((k,v) for k,v in presets2.items() if isinstance(v, float)) 
     963    for k, v in presets2.items(): 
     964        if not isinstance(v, float) and not k.endswith('_type'): 
     965            presets2[k] = eval(v, context) 
     966 
     967    # update parameters with presets 
     968    pars.update(presets)  # set value after random to control value 
     969    pars2.update(presets2)  # set value after random to control value 
     970    #import pprint; pprint.pprint(model_info) 
     971    constrain_pars(model_info, pars) 
     972    constrain_pars(model_info2, pars2) 
     973 
     974    same_model = name == name2 and pars == pars 
    897975    if len(engines) == 0: 
    898         engines.extend(['single', 'double']) 
     976        if same_model: 
     977            engines.extend(['single', 'double']) 
     978        else: 
     979            engines.extend(['single', 'single']) 
    899980    elif len(engines) == 1: 
    900         if engines[0][0] == 'double': 
     981        if not same_model: 
     982            engines.append(engines[0]) 
     983        elif engines[0] == 'double': 
    901984            engines.append('single') 
    902985        else: 
    … …  
    905988        del engines[2:] 
    906989 
    907     n1 = int(positional_args[1]) if len(positional_args) > 1 else 1 
    908     n2 = int(positional_args[2]) if len(positional_args) > 2 else 1 
    909990    use_sasview = any(engine == 'sasview' and count > 0 
    910991                      for engine, count in zip(engines, [n1, n2])) 
    911  
    912     # Get demo parameters from model definition, or use default parameters 
    913     # if model does not define demo parameters 
    914     pars = get_pars(model_info, opts['use_demo']) 
    915  
    916  
    917     # Fill in parameters given on the command line 
    918     presets = {} 
    919     for arg in values: 
    920         k, v = arg.split('=', 1) 
    921         if k not in pars: 
    922             # extract base name without polydispersity info 
    923             s = set(p.split('_pd')[0] for p in pars) 
    924             print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s)))) 
    925             return None 
    926         presets[k] = float(v) if not k.endswith('type') else v 
    927  
    928     # randomize parameters 
    929     #pars.update(set_pars)  # set value before random to control range 
    930     if opts['seed'] > -1: 
    931         pars = randomize_pars(model_info, pars, seed=opts['seed']) 
    932         print("Randomize using -random=%i"%opts['seed']) 
    933     if opts['mono']: 
    934         pars = suppress_pd(pars) 
    935     pars.update(presets)  # set value after random to control value 
    936     #import pprint; pprint.pprint(model_info) 
    937     constrain_pars(model_info, pars) 
    938992    if use_sasview: 
    939993        constrain_new_to_old(model_info, pars) 
     994        constrain_new_to_old(model_info2, pars2) 
     995 
    940996    if opts['show_pars']: 
    941997        print(str(parlist(model_info, pars, opts['is2d']))) 
    … …  
    9481004        base = None 
    9491005    if n2: 
    950         comp = make_engine(model_info, data, engines[1], opts['cutoff']) 
     1006        comp = make_engine(model_info2, data, engines[1], opts['cutoff']) 
    9511007    else: 
    9521008        comp = None 
    … …  
    9551011    # Remember it all 
    9561012    opts.update({ 
    957         'name'      : name, 
    958         'def'       : model_info, 
    959         'n1'        : n1, 
    960         'n2'        : n2, 
    961         'presets'   : presets, 
    962         'pars'      : pars, 
    9631013        'data'      : data, 
     1014        'name'      : [name, name2], 
     1015        'def'       : [model_info, model_info2], 
     1016        'count'     : [n1, n2], 
     1017        'presets'   : [presets, presets2], 
     1018        'pars'      : [pars, pars2], 
    9641019        'engines'   : [base, comp], 
    9651020    }) 
    … …  
    9761031    from .generate import view_html_from_info 
    9771032    app = wx.App() if wx.GetApp() is None else None 
    978     view_html_from_info(opts['def']) 
     1033    view_html_from_info(opts['def'][0]) 
    9791034    if app: app.MainLoop() 
    9801035 
    … …  
    10151070        config_matplotlib() 
    10161071        self.opts = opts 
    1017         model_info = opts['def'] 
    1018         pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars']) 
     1072        model_info = opts['def'][0] 
     1073        pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars'][0]) 
    10191074        # Initialize parameter ranges, fixing the 2D parameters for 1D data. 
    10201075        if not opts['is2d']: 
    … …  
    10621117        pars = dict((k, v.value) for k, v in self.pars.items()) 
    10631118        pars.update(self.pd_types) 
    1064         self.opts['pars'] = pars 
     1119        self.opts['pars'][0] = pars 
     1120        self.opts['pars'][1] = pars 
    10651121        limits = compare(self.opts, limits=self.limits) 
    10661122        if self.limits is None: 

  • sasmodels/conversion_table.py

    r5f1acda r68425bf  
    173173        "CSParallelepipedModel", 
    174174        { 
     175            "sld_core": "sld_pcore", 
     176            "sld_a": "sld_rimA", 
     177            "sld_b": "sld_rimB", 
     178            "sld_c": "sld_rimC", 
     179            "sld_solvent": "sld_solv", 
     180            "length_a": "shortA", 
     181            "length_b": "midB", 
     182            "length_c": "longC", 
     183            "thick_rim_a": "rimA", 
     184            "thick_rim_c": "rimC", 
     185            "thick_rim_b": "rimB", 
     186            "theta": "parallel_theta", 
    175187            "phi": "parallel_phi", 
    176188            "psi": "parallel_psi", 
    177             "sld_core": "sld_pcore", 
    178             "sld_c": "sld_rimC", 
    179             "sld_b": "sld_rimB", 
    180             "sld_solvent": "sld_solv", 
    181             "length_a": "shortA", 
    182             "sld_a": "sld_rimA", 
    183             "length_b": "midB", 
    184             "thick_rimc": "rimC", 
    185             "theta": "parallel_theta", 
    186             "thick_rim_a": "rimA", 
    187             "length_c": "longC", 
    188             "thick_rim_b": "rimB" 
    189189        } 
    190190    ], 
    … …  
    257257        "EllipticalCylinderModel", 
    258258        { 
     259            "axis_ratio": "r_ratio", 
     260            "radius_minor": "r_minor", 
     261            "sld": "sldCyl", 
     262            "sld_solvent": "sldSolv", 
     263            "theta": "cyl_theta", 
    259264            "phi": "cyl_phi", 
    260265            "psi": "cyl_psi", 
    261             "theta": "cyl_theta", 
    262             "sld": "sldCyl", 
    263             "axis_ratio": "r_ratio", 
    264             "sld_solvent": "sldSolv" 
    265266        } 
    266267    ], 
    … …  
    729730            "theta": "axis_theta", 
    730731            "sld_solvent": "solvent_sld", 
    731             "n_stacking": "n_stacking" 
     732            "n_stacking": "n_stacking", 
     733            "thick_layer": "layer_thick", 
     734            "thick_core": "core_thick", 
    732735        } 
    733736    ], 

  • sasmodels/kernel_header.c

    re1d6983 r11ca2ab  
    148148inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; } 
    149149 
     150 #define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 
     151     SINCOS(phi*M_PI_180, sn, cn); \ 
     152     q = sqrt(qx*qx + qy*qy); \ 
     153     cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * cos(theta*M_PI_180));  \ 
     154     sn = sqrt(1 - cn*cn); \ 
     155 } while (0) 
     156 
     157 #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ 
     158    q = sqrt(qx*qx + qy*qy); \ 
     159    const double qxhat = qx/q; \ 
     160    const double qyhat = qy/q; \ 
     161    double sin_theta, cos_theta; \ 
     162    double sin_phi, cos_phi; \ 
     163    double sin_psi, cos_psi; \ 
     164    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 
     165    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 
     166    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 
     167    cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat; \ 
     168    cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat; \ 
     169    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ 
     170 } while (0) 

  • sasmodels/models/barbell.c

    r0d6e865 r3a48772  
    3939    // translate dx in [-1,1] to dx in [lower,upper] 
    4040    const double integral = total*zm; 
    41     const double bell_Fq = 2*M_PI*cube(radius_bell)*integral; 
    42     return bell_Fq; 
     41    const double bell_fq = 2.0*M_PI*cube(radius_bell)*integral; 
     42    return bell_fq; 
    4343} 
     44 
     45static double 
     46_fq(double q, double h, 
     47    double radius_bell, double radius, double half_length, 
     48    double sin_alpha, double cos_alpha) 
     49{ 
     50    const double bell_fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha); 
     51    const double bj = sas_J1c(q*radius*sin_alpha); 
     52    const double si = sinc(q*half_length*cos_alpha); 
     53    const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
     54    const double Aq = bell_fq + cyl_fq; 
     55    return Aq; 
     56} 
     57 
    4458 
    4559double form_volume(double radius_bell, 
    … …  
    4761        double length) 
    4862{ 
    49  
    5063    // bell radius should never be less than radius when this is called 
    5164    const double hdist = sqrt(square(radius_bell) - square(radius)); 
    … …  
    7184        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    7285        SINCOS(alpha, sin_alpha, cos_alpha); 
    73  
    74         const double bell_Fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha); 
    75         const double bj = sas_J1c(q*radius*sin_alpha); 
    76         const double si = sinc(q*half_length*cos_alpha); 
    77         const double cyl_Fq = M_PI*radius*radius*length*bj*si; 
    78         const double Aq = bell_Fq + cyl_Fq; 
     86        const double Aq = _fq(q, h, radius_bell, radius, half_length, sin_alpha, cos_alpha); 
    7987        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
    8088    } 
    … …  
    93101        double theta, double phi) 
    94102{ 
    95     // Compute angle alpha between q and the cylinder axis 
    96     double sn, cn; // slots to hold sincos function output 
    97     SINCOS(phi*M_PI_180, sn, cn); 
    98     const double q = sqrt(qx*qx + qy*qy); 
    99     const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 
    100     const double alpha = acos(cos_val); // rod angle relative to q 
     103    double q, sin_alpha, cos_alpha; 
     104    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    101105 
    102106    const double h = -sqrt(square(radius_bell) - square(radius)); 
    103     const double half_length = 0.5*length; 
    104  
    105     double sin_alpha, cos_alpha; // slots to hold sincos function output 
    106     SINCOS(alpha, sin_alpha, cos_alpha); 
    107     const double bell_Fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha); 
    108     const double bj = sas_J1c(q*radius*sin_alpha); 
    109     const double si = sinc(q*half_length*cos_alpha); 
    110     const double cyl_Fq = M_PI*radius*radius*length*bj*si; 
    111     const double Aq = cyl_Fq + bell_Fq; 
     107    const double Aq = _fq(q, h, radius_bell, radius, 0.5*length, sin_alpha, cos_alpha); 
    112108 
    113109    // Multiply by contrast^2 and convert to cm-1 

  • sasmodels/models/bcc_paracrystal.c

    r0bef47b r4962519  
    3131        const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta; 
    3232        const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta; 
     33 
    3334        const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9)); 
    34         result = pow(1.0-(temp10*temp10),3)*temp6 
     35        result = cube(1.0 - (temp10*temp10))*temp6 
    3536            / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    3637              * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
    … …  
    8182 
    8283    return answer; 
    83  
    84  
    8584} 
    8685 
    8786 
    88 double Iqxy(double qx, double qy, double dnn, 
    89     double d_factor, double radius,double sld, double solvent_sld, 
    90     double theta, double phi, double psi){ 
     87double Iqxy(double qx, double qy, 
     88    double dnn, double d_factor, double radius, 
     89    double sld, double solvent_sld, 
     90    double theta, double phi, double psi) 
     91{ 
     92    double q, cos_a1, cos_a2, cos_a3; 
     93    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
    9194 
    92   double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z, 
    93   //double q_z; 
    94   double cos_val_b3, cos_val_b2, cos_val_b1; 
    95   double a1_dot_q, a2_dot_q,a3_dot_q; 
    96   double answer; 
    97   double Zq, Fkq, Fkq_2; 
     95    const double a1 = +cos_a3 - cos_a1 + cos_a2; 
     96    const double a2 = +cos_a3 + cos_a1 - cos_a2; 
     97    const double a3 = -cos_a3 + cos_a1 + cos_a2; 
    9898 
    99   //convert to q and make scaled values 
    100   double q = sqrt(qx*qx+qy*qy); 
    101   double q_x = qx/q; 
    102   double q_y = qy/q; 
     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); 
    103106 
    104   //convert angle degree to radian 
    105   theta = theta * M_PI_180; 
    106   phi = phi * M_PI_180; 
    107   psi = psi * M_PI_180; 
    108  
    109   const double Da = d_factor*dnn; 
    110   const double s1 = dnn/sqrt(0.75); 
    111  
    112  
    113   //the occupied volume of the lattice 
    114   const double latticescale = 2.0*sphere_volume(radius/s1); 
    115   // q vector 
    116   //q_z = 0.0; // for SANS; assuming qz is negligible 
    117   /// Angles here are respect to detector coordinate 
    118   ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854)) 
    119     // b3 axis orientation 
    120     b3_x = cos(theta) * cos(phi); 
    121     b3_y = sin(theta); 
    122     //b3_z = -cos(theta) * sin(phi); 
    123     cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z; 
    124  
    125     //alpha = acos(cos_val_b3); 
    126     // b1 axis orientation 
    127     b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi); 
    128     b1_y = sin(psi)*cos(theta); 
    129     cos_val_b1 = b1_x*q_x + b1_y*q_y; 
    130     // b2 axis orientation 
    131     b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi); 
    132         b2_y = cos(theta)*cos(psi); 
    133     cos_val_b2 = b2_x*q_x + b2_y*q_y; 
    134  
    135     // The following test should always pass 
    136     if (fabs(cos_val_b3)>1.0) { 
    137       //printf("bcc_ana_2D: Unexpected error: cos()>1\n"); 
    138       cos_val_b3 = 1.0; 
    139     } 
    140     if (fabs(cos_val_b2)>1.0) { 
    141       //printf("bcc_ana_2D: Unexpected error: cos()>1\n"); 
    142       cos_val_b2 = 1.0; 
    143     } 
    144     if (fabs(cos_val_b1)>1.0) { 
    145       //printf("bcc_ana_2D: Unexpected error: cos()>1\n"); 
    146       cos_val_b1 = 1.0; 
    147     } 
    148     // Compute the angle btw vector q and the a3 axis 
    149     a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3); 
    150  
    151     // a1 axis 
    152     a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1); 
    153  
    154     // a2 axis 
    155     a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2); 
    156  
    157  
    158     // Get Fkq and Fkq_2 
    159     Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q)); 
    160     Fkq_2 = Fkq*Fkq; 
    161     // Call Zq=Z1*Z2*Z3 
    162     Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2); 
    163     Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2); 
    164     Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2); 
    165  
    166   // Use SphereForm directly from libigor 
    167   answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale; 
    168  
    169   return answer; 
    170  } 
     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; 
     111} 

  • sasmodels/models/binary_hard_sphere.c

    re481a39 r4f79d94  
    77    ); 
    88     
    9 double Iqxy(double qx, double qy, 
    10     double lg_radius, double sm_radius, 
    11     double lg_vol_frac, double sm_vol_frac, 
    12     double lg_sld, double sm_sld, double solvent_sld 
    13     ); 
    14  
    159void calculate_psfs(double qval, 
    1610    double r2, double nf2, 
    … …  
    5549    // /* do form factor calculations  */ 
    5650     
    57     v1 = 4.0*M_PI/3.0*r1*r1*r1; 
    58     v2 = 4.0*M_PI/3.0*r2*r2*r2; 
     51    v1 = M_4PI_3*r1*r1*r1; 
     52    v2 = M_4PI_3*r2*r2*r2; 
    5953     
    6054    n1 = phi1/v1; 
    … …  
    6458    qr2 = r2*q; 
    6559 
    66     //if (qr1 == 0){ 
    67         //sc1 = 1.0/3.0; 
    68     //}else{ 
    69         //sc1 = (sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1; 
    70     //} 
    71     //if (qr2 == 0){ 
    72         //sc2 = 1.0/3.0; 
    73     //}else{ 
    74         //sc2 = (sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2; 
    75     //} 
    7660    sc1 = sph_j1c(qr1); 
    7761    sc2 = sph_j1c(qr2); 
    78     b1 = r1*r1*r1*(rho1-rhos)*4.0/3.0*M_PI*sc1; 
    79     b2 = r2*r2*r2*(rho2-rhos)*4.0/3.0*M_PI*sc2; 
     62    b1 = r1*r1*r1*(rho1-rhos)*M_4PI_3*sc1; 
     63    b2 = r2*r2*r2*(rho2-rhos)*M_4PI_3*sc2; 
    8064    inten = n1*b1*b1*psf11; 
    8165    inten += sqrt(n1*n2)*2.0*b1*b2*psf12; 
    … …  
    8872} 
    8973 
    90  
    91 double Iqxy(double qx, double qy, 
    92     double lg_radius, double sm_radius, 
    93     double lg_vol_frac, double sm_vol_frac, 
    94     double lg_sld, double sm_sld, double solvent_sld) 
    95      
    96 { 
    97     double q = sqrt(qx*qx + qy*qy); 
    98     return Iq(q, 
    99         lg_radius, sm_radius, 
    100         lg_vol_frac, sm_vol_frac, 
    101         lg_sld, sm_sld, solvent_sld); 
    102 } 
    10374 
    10475void calculate_psfs(double qval, 

  • sasmodels/models/capped_cylinder.c

    r0d6e865 r3a48772  
    4545    // translate dx in [-1,1] to dx in [lower,upper] 
    4646    const double integral = total*zm; 
    47     const double cap_Fq = 2*M_PI*cube(radius_cap)*integral; 
     47    const double cap_Fq = 2.0*M_PI*cube(radius_cap)*integral; 
    4848    return cap_Fq; 
     49} 
     50 
     51static double 
     52_fq(double q, double h, double radius_cap, double radius, double half_length, 
     53    double sin_alpha, double cos_alpha) 
     54{ 
     55    const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha); 
     56    const double bj = sas_J1c(q*radius*sin_alpha); 
     57    const double si = sinc(q*half_length*cos_alpha); 
     58    const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
     59    const double Aq = cap_Fq + cyl_Fq; 
     60    return Aq; 
    4961} 
    5062 
    … …  
    93105        SINCOS(alpha, sin_alpha, cos_alpha); 
    94106 
    95         const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha); 
    96         const double bj = sas_J1c(q*radius*sin_alpha); 
    97         const double si = sinc(q*half_length*cos_alpha); 
    98         const double cyl_Fq = M_PI*radius*radius*length*bj*si; 
    99         const double Aq = cap_Fq + cyl_Fq; 
    100         total += Gauss76Wt[i] * Aq * Aq * sin_alpha; // sin_alpha for spherical coord integration 
     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; 
    101110    } 
    102111    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    114123    double theta, double phi) 
    115124{ 
    116     // Compute angle alpha between q and the cylinder axis 
    117     double sn, cn; 
    118     SINCOS(phi*M_PI_180, sn, cn); 
    119     const double q = sqrt(qx*qx+qy*qy); 
    120     const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 
    121     const double alpha = acos(cos_val); // rod angle relative to q 
     125    double q, sin_alpha, cos_alpha; 
     126    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    122127 
    123128    const double h = sqrt(radius_cap*radius_cap - radius*radius); 
    124     const double half_length = 0.5*length; 
    125  
    126     double sin_alpha, cos_alpha; // slots to hold sincos function output 
    127     SINCOS(alpha, sin_alpha, cos_alpha); 
    128     const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha); 
    129     const double bj = sas_J1c(q*radius*sin_alpha); 
    130     const double si = sinc(q*half_length*cos_alpha); 
    131     const double cyl_Fq = M_PI*radius*radius*length*bj*si; 
    132     const double Aq = cap_Fq + cyl_Fq; 
     129    const double Aq = _fq(q, h, radius_cap, radius, 0.5*length, sin_alpha, cos_alpha); 
    133130 
    134131    // Multiply by contrast^2 and convert to cm-1 

  • sasmodels/models/core_shell_bicelle.c

    r0d6e865 r5bddd89  
    2626double form_volume(double radius, double thick_rim, double thick_face, double length) 
    2727{ 
    28     return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2*thick_face); 
     28    return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2.0*thick_face); 
    2929} 
    3030 
    … …  
    3939              double rhor, 
    4040              double rhosolv, 
    41               double dum) 
     41              double sin_alpha, 
     42              double cos_alpha) 
    4243{ 
    4344    double si1,si2,be1,be2; 
    … …  
    4950    const double vol2 = M_PI*(rad+radthick)*(rad+radthick)*2.0*(length+facthick); 
    5051    const double vol3 = M_PI*rad*rad*2.0*(length+facthick); 
    51     double sn,cn; 
    52     SINCOS(dum, sn, cn); 
    53     double besarg1 = qq*rad*sn; 
    54     double besarg2 = qq*(rad+radthick)*sn; 
    55     double sinarg1 = qq*length*cn; 
    56     double sinarg2 = qq*(length+facthick)*cn; 
     52    double besarg1 = qq*rad*sin_alpha; 
     53    double besarg2 = qq*(rad+radthick)*sin_alpha; 
     54    double sinarg1 = qq*length*cos_alpha; 
     55    double sinarg2 = qq*(length+facthick)*cos_alpha; 
    5756 
    5857    be1 = sas_J1c(besarg1); 
    … …  
    6564                     vol3*dr3*si2*be1; 
    6665 
    67     const double retval = t*t*sn; 
     66    const double retval = t*t*sin_alpha; 
    6867 
    69     return(retval); 
     68    return retval; 
    7069 
    7170} 
    … …  
    8382{ 
    8483    // set up the integration end points 
    85     const double uplim = M_PI/4; 
    86     const double halfheight = length/2.0; 
     84    const double uplim = M_PI_4; 
     85    const double halfheight = 0.5*length; 
    8786 
    8887    double summ = 0.0; 
    8988    for(int i=0;i<N_POINTS_76;i++) { 
    90         double zi = (Gauss76Z[i] + 1.0)*uplim; 
     89        double alpha = (Gauss76Z[i] + 1.0)*uplim; 
     90        double sin_alpha, cos_alpha; // slots to hold sincos function output 
     91        SINCOS(alpha, sin_alpha, cos_alpha); 
    9192        double yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick, 
    92                              halfheight, rhoc, rhoh, rhor,rhosolv, zi); 
     93                             halfheight, rhoc, rhoh, rhor, rhosolv, 
     94                             sin_alpha, cos_alpha); 
    9395        summ += yyy; 
    9496    } 
    … …  
    9698    // calculate value of integral to return 
    9799    double answer = uplim*summ; 
    98     return(answer); 
     100    return answer; 
    99101} 
    100102 
    101103static double 
    102 bicelle_kernel_2d(double q, double q_x, double q_y, 
     104bicelle_kernel_2d(double qx, double qy, 
    103105          double radius, 
    104106          double thick_rim, 
    … …  
    112114          double phi) 
    113115{ 
    114     //convert angle degree to radian 
    115     theta *= M_PI_180; 
    116     phi *= M_PI_180; 
     116    double q, sin_alpha, cos_alpha; 
     117    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    117118 
    118     // Cylinder orientation 
    119     const double cyl_x = sin(theta) * cos(phi); 
    120     const double cyl_y = sin(theta) * sin(phi); 
    121  
    122     // Compute the angle btw vector q and the axis of the cylinder 
    123     const double cos_val = cyl_x*q_x + cyl_y*q_y; 
    124     const double alpha = acos( cos_val ); 
    125  
    126     // Get the kernel 
    127119    double answer = bicelle_kernel(q, radius, thick_rim, thick_face, 
    128                            length/2.0, core_sld, face_sld, rim_sld, 
    129                            solvent_sld, alpha) / fabs(sin(alpha)); 
     120                           0.5*length, core_sld, face_sld, rim_sld, 
     121                           solvent_sld, sin_alpha, cos_alpha) / fabs(sin_alpha); 
    130122 
    131123    answer *= 1.0e-4; 
    … …  
    162154          double phi) 
    163155{ 
    164     double q; 
    165     q = sqrt(qx*qx+qy*qy); 
    166     double intensity = bicelle_kernel_2d(q, qx/q, qy/q, 
     156    double intensity = bicelle_kernel_2d(qx, qy, 
    167157                      radius, 
    168158                      thick_rim, 

  • sasmodels/models/core_shell_cylinder.c

    r0d6e865 r5bddd89  
    1111double _cyl(double twovd, double besarg, double siarg) 
    1212{ 
    13     const double bj = (besarg == 0.0 ? 0.5 : 0.5*sas_J1c(besarg)); 
    14     const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 
    15     return twovd*si*bj; 
     13    return twovd * sinc(siarg) * sas_J1c(besarg); 
    1614} 
    1715 
    1816double form_volume(double radius, double thickness, double length) 
    1917{ 
    20     return M_PI*(radius+thickness)*(radius+thickness)*(length+2*thickness); 
     18    return M_PI*(radius+thickness)*(radius+thickness)*(length+2.0*thickness); 
    2119} 
    2220 
    … …  
    6664    double phi) 
    6765{ 
    68     double sn, cn; // slots to hold sincos function output 
    69  
    70     // Compute angle alpha between q and the cylinder axis 
    71     SINCOS(phi*M_PI_180, sn, cn); 
    72     // # The following correction factor exists in sasview, but it can't be 
    73     // # right, so we are leaving it out for now. 
    74     // const double correction = fabs(cn)*M_PI_2; 
    75     const double q = sqrt(qx*qx+qy*qy); 
    76     const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 
    77     const double alpha = acos(cos_val); 
     66    double q, sin_alpha, cos_alpha; 
     67    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    7868 
    7969    const double core_qr = q*radius; 
    … …  
    8676                             * (shell_sld-solvent_sld); 
    8777 
    88     SINCOS(alpha, sn, cn); 
    89     const double fq = _cyl(core_twovd, core_qr*sn, core_qh*cn) 
    90         + _cyl(shell_twovd, shell_qr*sn, shell_qh*cn); 
     78    const double fq = _cyl(core_twovd, core_qr*sin_alpha, core_qh*cos_alpha) 
     79        + _cyl(shell_twovd, shell_qr*sin_alpha, shell_qh*cos_alpha); 
    9180    return 1.0e-4 * fq * fq; 
    9281} 

  • sasmodels/models/core_shell_ellipsoid.c

    r0d6e865 r5bddd89  
    3232    const double equat_shell = radius_equat_core + thick_shell; 
    3333    const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 
    34     double vol = 4.0*M_PI/3.0*equat_shell*equat_shell*polar_shell; 
     34    double vol = M_4PI_3*equat_shell*equat_shell*polar_shell; 
    3535    return vol; 
    3636} 
    … …  
    6060 
    6161    for(int i=0;i<N_POINTS_76;i++) { 
    62         double zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 
     62        double zi = 0.5*( Gauss76Z[i]*(uplim-lolim) + uplim + lolim ); 
    6363        double yyy = Gauss76Wt[i] * gfn4(zi, 
    6464                                  radius_equat_core, 
    … …  
    7272    } 
    7373 
    74     double answer = (uplim-lolim)/2.0*summ; 
     74    double answer = 0.5*(uplim-lolim)*summ; 
    7575    //convert to [cm-1] 
    7676    answer *= 1.0e-4; 
    … …  
    8080 
    8181static double 
    82 core_shell_ellipsoid_xt_kernel_2d(double q, double q_x, double q_y, 
     82core_shell_ellipsoid_xt_kernel_2d(double qx, double qy, 
    8383          double radius_equat_core, 
    8484          double x_core, 
    … …  
    9191          double phi) 
    9292{ 
    93     //convert angle degree to radian 
    94     theta = theta * M_PI_180; 
    95     phi = phi * M_PI_180; 
    96  
    97     // ellipsoid orientation, the axis of the rotation is consistent with the ploar axis. 
    98     const double cyl_x = sin(theta) * cos(phi); 
    99     const double cyl_y = sin(theta) * sin(phi); 
     93    double q, sin_alpha, cos_alpha; 
     94    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    10095 
    10196    const double sldcs = core_sld - shell_sld; 
    10297    const double sldss = shell_sld- solvent_sld; 
    103  
    104     // Compute the angle btw vector q and the 
    105     // axis of the cylinder 
    106     const double cos_val = cyl_x*q_x + cyl_y*q_y; 
    10798 
    10899    const double polar_core = radius_equat_core*x_core; 
    … …  
    112103    // Call the IGOR library function to get the kernel: 
    113104    // MUST use gfn4 not gf2 because of the def of params. 
    114     double answer = gfn4(cos_val, 
     105    double answer = gfn4(cos_alpha, 
    115106                  radius_equat_core, 
    116107                  polar_core, 
    … …  
    160151          double phi) 
    161152{ 
    162     double q; 
    163     q = sqrt(qx*qx+qy*qy); 
    164     double intensity = core_shell_ellipsoid_xt_kernel_2d(q, qx/q, qy/q, 
     153    double intensity = core_shell_ellipsoid_xt_kernel_2d(qx, qy, 
    165154                       radius_equat_core, 
    166155                       x_core, 

  • sasmodels/models/core_shell_parallelepiped.c

    r2222134 r14838a3  
    3535    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    3636     
    37     double t1, t2, t3, t4, tmp, answer;    
    38     double mu = q * length_b; 
     37    const double mu = 0.5 * q * length_b; 
    3938     
    4039    //calculate volume before rescaling (in original code, but not used) 
    41     //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);             
     40    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);         
    4241    //double vol = length_a * length_b * length_c +  
    4342    //       2.0 * thick_rim_a * length_b * length_c +  
    … …  
    4645     
    4746    // Scale sides by B 
    48     double a_scaled = length_a / length_b; 
    49     double c_scaled = length_c / length_b; 
     47    const double a_scaled = length_a / length_b; 
     48    const double c_scaled = length_c / length_b; 
    5049 
    51     // DelRho values (note that drC is not used later)        
    52         double dr0 = core_sld-solvent_sld; 
    53         double drA = arim_sld-solvent_sld; 
    54         double drB = brim_sld-solvent_sld; 
    55         //double drC = crim_sld-solvent_sld; 
     50    // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG) 
     51    // the a_scaled in the definition of tb was present in CSPPKernel in libCylinder.c, 
     52    // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions 
     53    // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!) 
     54    double ta = (a_scaled + 2.0*thick_rim_a)/length_b; 
     55    double tb = (a_scaled + 2.0*thick_rim_b)/length_b; 
    5656 
    57     //Order of integration 
    58     int nordi=76;                                
    59     int nordj=76; 
    60      
    61     // outer integral (with nordi gauss points), integration limits = 0, 1 
    62     double summ = 0; //initialize integral 
     57    double Vin = length_a * length_b * length_c; 
     58    //double Vot = (length_a * length_b * length_c + 
     59    //            2.0 * thick_rim_a * length_b * length_c + 
     60    //            2.0 * length_a * thick_rim_b * length_c + 
     61    //            2.0 * length_a * length_b * thick_rim_c); 
     62    double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 
     63    double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
    6364 
    64     for( int i=0; i<nordi; i++) { 
    65                  
    66         // inner integral (with nordj gauss points), integration limits = 0, 1 
    67          
    68         double summj = 0.0; 
    69             double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );          
    70                  
    71             for(int j=0; j<nordj; j++) { 
     65    // Scale factors (note that drC is not used later) 
     66    const double drho0 = (core_sld-solvent_sld); 
     67    const double drhoA = (arim_sld-solvent_sld); 
     68    const double drhoB = (brim_sld-solvent_sld); 
     69    //const double drC_Vot = (crim_sld-solvent_sld)*Vot; 
    7270 
    73             double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    74             double mudum = mu * sqrt(1.0-sigma*sigma); 
     71    // Precompute scale factors for combining cross terms from the shape 
     72    const double scale23 = drhoA*V1; 
     73    const double scale14 = drhoB*V2; 
     74    const double scale12 = drho0*Vin - scale23 - scale14; 
    7575 
    76                 double Vin = length_a * length_b * length_c; 
    77                 //double Vot = (length_a * length_b * length_c + 
    78             //            2.0 * thick_rim_a * length_b * length_c + 
    79             //            2.0 * length_a * thick_rim_b * length_c + 
    80             //            2.0 * length_a * length_b * thick_rim_c); 
    81                 double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 
    82                 double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
    83          
    84             // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG) 
    85             // the a_scaled in the definition of tb was present in CSPPKernel in libCylinder.c, 
    86             // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions 
    87             // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!)             
    88             double ta = (a_scaled+2.0*thick_rim_a)/length_b;  
    89             double tb = (a_scaled+2.0*thick_rim_b)/length_b; 
    90      
    91                 double arg1 = (0.5*mudum*a_scaled) * sin(0.5*M_PI*uu); 
    92                 double arg2 = (0.5*mudum) * cos(0.5*M_PI*uu); 
    93                 double arg3=  (0.5*mudum*ta) * sin(0.5*M_PI*uu); 
    94                 double arg4=  (0.5*mudum*tb) * cos(0.5*M_PI*uu); 
     76    // outer integral (with gauss points), integration limits = 0, 1 
     77    double outer_total = 0; //initialize integral 
    9578 
    96                 if(arg1==0.0){ 
    97                         t1 = 1.0; 
    98                 } else { 
    99                         t1 = (sin(arg1)/arg1);                //defn for CSPP model sin(arg1)/arg1    test:  (sin(arg1)/arg1)*(sin(arg1)/arg1)    
    100                 } 
    101                 if(arg2==0.0){ 
    102                         t2 = 1.0; 
    103                 } else { 
    104                         t2 = (sin(arg2)/arg2);           //defn for CSPP model sin(arg2)/arg2   test: (sin(arg2)/arg2)*(sin(arg2)/arg2)     
    105                 }        
    106                 if(arg3==0.0){ 
    107                         t3 = 1.0; 
    108                 } else { 
    109                         t3 = sin(arg3)/arg3; 
    110                 } 
    111                 if(arg4==0.0){ 
    112                         t4 = 1.0; 
    113                 } else { 
    114                         t4 = sin(arg4)/arg4; 
    115                 } 
    116              
     79    for( int i=0; i<76; i++) { 
     80        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
     81        double mu_proj = mu * sqrt(1.0-sigma*sigma); 
     82 
     83        // inner integral (with gauss points), integration limits = 0, 1 
     84        double inner_total = 0.0; 
     85        for(int j=0; j<76; j++) { 
     86            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
     87            double sin_uu, cos_uu; 
     88            SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
     89            const double si1 = sinc(mu_proj * sin_uu * a_scaled); 
     90            const double si2 = sinc(mu_proj * cos_uu); 
     91            const double si3 = sinc(mu_proj * sin_uu * ta); 
     92            const double si4 = sinc(mu_proj * cos_uu * tb); 
     93 
    11794            // Expression in libCylinder.c (neither drC nor Vot are used) 
    118                 tmp =( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 )*( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 );   //  correct FF : square of sum of phase factors             
    119              
     95            const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4; 
     96 
    12097            // To note also that in csparallelepiped.cpp, there is a function called 
    12198            // cspkernel, which seems to make more sense and has the following comment: 
    … …  
    126103            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c         
    127104             
    128             summj += Gauss76Wt[j] * tmp; 
     105            //  correct FF : sum of square of phase factors 
     106            inner_total += Gauss76Wt[j] * form * form; 
     107        } 
     108        inner_total *= 0.5; 
    129109 
    130         } 
    131                  
    132         // value of the inner integral 
    133         answer = 0.5 * summj; 
     110        // now sum up the outer integral 
     111        const double si = sinc(mu * c_scaled * sigma); 
     112        outer_total += Gauss76Wt[i] * inner_total * si * si; 
     113    } 
     114    outer_total *= 0.5; 
    134115 
    135                 // finish the outer integral  
    136                 double arg = 0.5 * mu* c_scaled *sigma; 
    137                 if ( arg == 0.0 ) { 
    138                         answer *= 1.0; 
    139                 } else { 
    140                 answer *= sin(arg)*sin(arg)/arg/arg; 
    141                 } 
    142                  
    143                 // now sum up the outer integral 
    144                 summ += Gauss76Wt[i] * answer; 
    145  
    146     } 
    147      
    148         answer = 0.5 * summ; 
    149  
    150         //convert from [1e-12 A-1] to [cm-1] 
    151         answer *= 1.0e-4; 
    152          
    153         return answer; 
     116    //convert from [1e-12 A-1] to [cm-1] 
     117    return 1.0e-4 * outer_total; 
    154118} 
    155119 
    … …  
    170134    double psi) 
    171135{ 
    172     double tmp1, tmp2, tmp3, tmpt1, tmpt2, tmpt3;    
     136    double q, cos_val_a, cos_val_b, cos_val_c; 
     137    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a); 
    173138 
    174     double q = sqrt(qx*qx+qy*qy); 
    175     double qx_scaled = qx/q; 
    176     double qy_scaled = qy/q; 
     139    // cspkernel in csparallelepiped recoded here 
     140    const double dr0 = core_sld-solvent_sld; 
     141    const double drA = arim_sld-solvent_sld; 
     142    const double drB = brim_sld-solvent_sld; 
     143    const double drC = crim_sld-solvent_sld; 
    177144 
    178     // Convert angles given in degrees to radians 
    179     theta *= M_PI_180; 
    180     phi   *= M_PI_180; 
    181     psi   *= M_PI_180; 
    182      
    183     // Parallelepiped c axis orientation 
    184     double cparallel_x = cos(theta) * cos(phi); 
    185     double cparallel_y = sin(theta); 
    186      
    187     // Compute angle between q and parallelepiped axis 
    188     double cos_val_c = cparallel_x*qx_scaled + cparallel_y*qy_scaled;// + cparallel_z*qz; 
    189  
    190     // Parallelepiped a axis orientation 
    191     double parallel_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi); 
    192     double parallel_y = sin(psi)*cos(theta); 
    193     double cos_val_a = parallel_x*qx_scaled + parallel_y*qy_scaled; 
    194  
    195     // Parallelepiped b axis orientation 
    196     double bparallel_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi); 
    197     double bparallel_y = cos(theta)*cos(psi); 
    198     double cos_val_b = bparallel_x*qx_scaled + bparallel_y*qy_scaled; 
    199  
    200     // The following tests should always pass 
    201     if (fabs(cos_val_c)>1.0) { 
    202       //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    203       cos_val_c = 1.0; 
    204     } 
    205     if (fabs(cos_val_a)>1.0) { 
    206       //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    207       cos_val_a = 1.0; 
    208     } 
    209     if (fabs(cos_val_b)>1.0) { 
    210       //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    211       cos_val_b = 1.0; 
    212     } 
    213      
    214     // cspkernel in csparallelepiped recoded here  
    215     double dr0 = core_sld-solvent_sld; 
    216         double drA = arim_sld-solvent_sld; 
    217         double drB = brim_sld-solvent_sld; 
    218         double drC = crim_sld-solvent_sld; 
    219         double Vin = length_a * length_b * length_c; 
     145    double Vin = length_a * length_b * length_c; 
     146    double V1 = 2.0 * thick_rim_a * length_b * length_c;    // incorrect V1(aa*bb*cc+2*ta*bb*cc) 
     147    double V2 = 2.0 * length_a * thick_rim_b * length_c;    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
     148    double V3 = 2.0 * length_a * length_b * thick_rim_c; 
    220149    // As for the 1D case, Vot is not used 
    221         //double Vot = (length_a * length_b * length_c + 
     150    //double Vot = (length_a * length_b * length_c + 
    222151    //              2.0 * thick_rim_a * length_b * length_c + 
    223152    //              2.0 * length_a * thick_rim_b * length_c + 
    224153    //              2.0 * length_a * length_b * thick_rim_c); 
    225         double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 
    226         double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
    227     double V3 = (2.0 * length_a * length_b * thick_rim_c); 
     154 
    228155    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    229156    // the scaling by B. The use of length_a for the 3 of them seems clearly a mistake to me, 
    230157    // but for the moment I let it like this until understanding better the code. 
    231         double ta = length_a + 2.0*thick_rim_a; 
     158    double ta = length_a + 2.0*thick_rim_a; 
    232159    double tb = length_a + 2.0*thick_rim_b; 
    233160    double tc = length_a + 2.0*thick_rim_c; 
    234161    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    235     double argA = 0.5*q*length_a*cos_val_a; 
    236     double argB = 0.5*q*length_b*cos_val_b; 
    237     double argC = 0.5*q*length_c*cos_val_c; 
    238     double argtA = 0.5*q*ta*cos_val_a; 
    239     double argtB = 0.5*q*tb*cos_val_b; 
    240     double argtC = 0.5*q*tc*cos_val_c; 
     162    double siA = sinc(0.5*q*length_a*cos_val_a); 
     163    double siB = sinc(0.5*q*length_b*cos_val_b); 
     164    double siC = sinc(0.5*q*length_c*cos_val_c); 
     165    double siAt = sinc(0.5*q*ta*cos_val_a); 
     166    double siBt = sinc(0.5*q*tb*cos_val_b); 
     167    double siCt = sinc(0.5*q*tc*cos_val_c); 
    241168     
    242     if(argA==0.0) { 
    243         tmp1 = 1.0; 
    244     } else { 
    245         tmp1 = sin(argA)/argA; 
    246     } 
    247     if (argB==0.0) { 
    248         tmp2 = 1.0; 
    249     } else { 
    250         tmp2 = sin(argB)/argB; 
    251     } 
    252     if (argC==0.0) { 
    253         tmp3 = 1.0; 
    254     } else { 
    255         tmp3 = sin(argC)/argC; 
    256     } 
    257     if(argtA==0.0) { 
    258         tmpt1 = 1.0; 
    259     } else { 
    260         tmpt1 = sin(argtA)/argtA; 
    261     } 
    262     if (argtB==0.0) { 
    263         tmpt2 = 1.0; 
    264     } else { 
    265         tmpt2 = sin(argtB)/argtB; 
    266     } 
    267     if (argtC==0.0) { 
    268         tmpt3 = 1.0; 
    269     } else { 
    270         tmpt3 = sin(argtC)*sin(argtC)/argtC/argtC; 
    271     } 
    272169 
    273170    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
    274171    // in the 1D code, but should be checked! 
    275     double f = ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1 +  
    276                drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3); 
     172    double f = ( dr0*siA*siB*siC*Vin 
     173               + drA*(siAt-siA)*siB*siC*V1 
     174               + drB*siA*(siBt-siB)*siC*V2 
     175               + drC*siA*siB*(siCt*siCt-siC)*V3); 
    277176    
    278177    return 1.0e-4 * f * f; 

  • sasmodels/models/core_shell_parallelepiped.py

    r2222134 r14838a3  
    164164# parameters for demo 
    165165demo = dict(scale=1, background=0.0, 
    166             sld_core=1e-6, sld_a=2e-6, sld_b=4e-6, 
    167             sld_c=2e-6, sld_solvent=6e-6, 
     166            sld_core=1, sld_a=2, sld_b=4, sld_c=2, sld_solvent=6, 
    168167            length_a=35, length_b=75, length_c=400, 
    169168            thick_rim_a=10, thick_rim_b=10, thick_rim_c=10, 
    … …  
    177176            theta_pd=10, theta_pd_n=1, 
    178177            phi_pd=10, phi_pd_n=1, 
    179             psi_pd=10, psi_pd_n=10) 
     178            psi_pd=10, psi_pd_n=1) 
    180179 
    181180qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 

  • sasmodels/models/core_shell_sphere.c

    r2c74c11 r3a48772  
    1616double form_volume(double radius, double thickness) 
    1717{ 
    18     return 4.0 * M_PI / 3.0 * pow((radius + thickness), 3); 
     18    return M_4PI_3 * cube(radius + thickness); 
    1919} 

  • sasmodels/models/core_shell_sphere.py

    r9a4811a r4962519  
    9595    """ 
    9696    return (1, 1) 
    97     whole = 4.0 * pi / 3.0 * pow((radius + thickness), 3) 
    98     core = 4.0 * pi / 3.0 * radius * radius * radius 
     97    whole = 4.0/3.0 * pi * (radius + thickness)**3 
     98    core = 4.0/3.0 * pi * radius**3 
    9999    return whole, whole - core 
    100100 

  • sasmodels/models/cylinder.c

    r0d6e865 r11ca2ab  
    3333        // alpha(theta,phi) the projection of the cylinder on the detector plane 
    3434        SINCOS(alpha, sn, cn); 
    35         total += Gauss76Wt[i] * fq(q, sn, cn, radius, length)* fq(q, sn, cn, radius, length) * sn; 
     35        total += Gauss76Wt[i] * square(fq(q, sn, cn, radius, length)) * sn; 
    3636    } 
    3737    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    5858    double phi) 
    5959{ 
    60     double sn, cn; // slots to hold sincos function output 
    61  
    62     // Compute angle alpha between q and the cylinder axis 
    63     SINCOS(phi*M_PI_180, sn, cn); 
    64     const double q = sqrt(qx*qx + qy*qy); 
    65     const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 
    66  
    67     const double alpha = acos(cos_val); 
    68  
    69     SINCOS(alpha, sn, cn); 
     60    double q, sin_alpha, cos_alpha; 
     61    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    7062    const double s = (sld-solvent_sld) * form_volume(radius, length); 
    71     return 1.0e-4 * square(s * fq(q, sn, cn, radius, length)); 
     63    return 1.0e-4 * square(s * fq(q, sin_alpha, cos_alpha, radius, length)); 
    7264} 

  • sasmodels/models/dab.py

    ra807206 r4962519  
    5959 
    6060Iq = """ 
    61     double numerator   = pow(cor_length, 3); 
    62     double denominator = pow(1 + pow(q*cor_length,2), 2); 
     61    double numerator   = cube(cor_length); 
     62    double denominator = square(1 + square(q*cor_length)); 
    6363     
    6464    return numerator / denominator ; 

  • sasmodels/models/ellipsoid.c

    r0d6e865 r5bddd89  
    4949    double phi) 
    5050{ 
    51     double sn, cn; 
    52  
    53     const double q = sqrt(qx*qx + qy*qy); 
    54     SINCOS(phi*M_PI_180, sn, cn); 
    55     const double cos_alpha = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 
    56     const double alpha = acos(cos_alpha); 
    57     SINCOS(alpha, sn, cn); 
    58     const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, sn); 
     51    double q, sin_alpha, cos_alpha; 
     52    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     53    const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, sin_alpha); 
    5954    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    6055 

  • sasmodels/models/elliptical_cylinder.c

    ra807206 r68425bf  
    66 
    77 
    8 double _elliptical_cylinder_kernel(double q, double radius_minor, double r_ratio, double theta); 
    98 
    10 double _elliptical_cylinder_kernel(double q, double radius_minor, double r_ratio, double theta) 
     9static double 
     10_kernel(double q, double radius_minor, double r_ratio, double theta) 
    1111{ 
    1212    // This is the function LAMBDA1^2 in Feigin's notation 
    … …  
    3030} 
    3131 
    32 double Iq(double q, double radius_minor, double r_ratio, double length, 
    33           double sld, double solvent_sld) { 
     32double 
     33Iq(double q, double radius_minor, double r_ratio, double length, 
     34   double sld, double solvent_sld) 
     35{ 
     36    // orientational average limits 
     37    const double va = 0.0; 
     38    const double vb = 1.0; 
     39    // inner integral limits 
     40    const double vaj=0.0; 
     41    const double vbj=M_PI; 
    3442 
    35     const int nordi=76; //order of integration 
    36     const int nordj=20; 
    37     double va,vb;       //upper and lower integration limits 
    38     double summ,zi,yyy,answer;         //running tally of integration 
    39     double summj,vaj,vbj,zij,arg,si;            //for the inner integration 
    40  
    41     // orientational average limits 
    42     va = 0.0; 
    43     vb = 1.0; 
    44     // inner integral limits 
    45     vaj=0.0; 
    46     vbj=M_PI; 
     43    const double radius_major = r_ratio * radius_minor; 
     44    const double rA = 0.5*(square(radius_major) + square(radius_minor)); 
     45    const double rB = 0.5*(square(radius_major) - square(radius_minor)); 
    4746 
    4847    //initialize integral 
    49     summ = 0.0; 
    50  
    51     const double delrho = sld - solvent_sld; 
    52  
    53     for(int i=0;i<nordi;i++) { 
     48    double outer_sum = 0.0; 
     49    for(int i=0;i<76;i++) { 
    5450        //setup inner integral over the ellipsoidal cross-section 
    55         summj=0; 
    56         zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;     //the "x" dummy 
    57         arg = radius_minor*sqrt(1.0-zi*zi); 
    58         for(int j=0;j<nordj;j++) { 
     51        const double cos_val = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
     52        const double sin_val = sqrt(1.0 - cos_val*cos_val); 
     53        //const double arg = radius_minor*sin_val; 
     54        double inner_sum=0; 
     55        for(int j=0;j<20;j++) { 
    5956            //20 gauss points for the inner integral 
    60             zij = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0;        //the "y" dummy 
    61             yyy = Gauss20Wt[j] * _elliptical_cylinder_kernel(q, arg, r_ratio, zij); 
    62             summj += yyy; 
     57            const double theta = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     58            const double r = sin_val*sqrt(rA - rB*cos(theta)); 
     59            const double be = sas_J1c(q*r); 
     60            inner_sum += Gauss20Wt[j] * be * be; 
    6361        } 
    6462        //now calculate the value of the inner integral 
    65         answer = (vbj-vaj)/2.0*summj; 
     63        inner_sum *= 0.5*(vbj-vaj); 
    6664 
    6765        //now calculate outer integral 
    68         arg = q*length*zi/2.0; 
    69         si = square(sinc(arg)); 
    70         yyy = Gauss76Wt[i] * answer * si; 
    71         summ += yyy; 
     66        const double si = sinc(q*0.5*length*cos_val); 
     67        outer_sum += Gauss76Wt[i] * inner_sum * si * si; 
    7268    } 
     69    outer_sum *= 0.5*(vb-va); 
    7370 
    7471    //divide integral by Pi 
    75     answer = (vb-va)/2.0*summ/M_PI; 
    76     // Multiply by contrast^2 
    77     answer *= delrho*delrho; 
     72    const double form = outer_sum/M_PI; 
    7873 
     74    // scale by contrast and volume, and convert to to 1/cm units 
    7975    const double vol = form_volume(radius_minor, r_ratio, length); 
    80     return answer*vol*vol*1.0e-4; 
     76    const double delrho = sld - solvent_sld; 
     77    return 1.0e-4*square(delrho*vol)*form; 
    8178} 
    8279 
    8380 
    84 double Iqxy(double qx, double qy, double radius_minor, double r_ratio, double length, 
    85             double sld, double solvent_sld, double theta, double phi, double psi) { 
    86     const double _theta = theta * M_PI / 180.0; 
    87     const double _phi = phi * M_PI / 180.0; 
    88     const double _psi = psi * M_PI / 180.0; 
    89     const double q = sqrt(qx*qx+qy*qy); 
    90     const double q_x = qx/q; 
    91     const double q_y = qy/q; 
     81double 
     82Iqxy(double qx, double qy, 
     83     double radius_minor, double r_ratio, double length, 
     84     double sld, double solvent_sld, 
     85     double theta, double phi, double psi) 
     86{ 
     87    double q, cos_val, cos_mu, cos_nu; 
     88    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val, cos_mu, cos_nu); 
    9289 
    93     //Cylinder orientation 
    94     double cyl_x = cos(_theta) * cos(_phi); 
    95     double cyl_y = sin(_theta); 
    96  
    97     //cyl_z = -cos(_theta) * sin(_phi); 
    98  
    99     // q vector 
    100     //q_z = 0; 
    101  
    102     // Note: cos(alpha) = 0 and 1 will get an 
    103     // undefined value from CylKernel 
    104     //alpha = acos( cos_val ); 
    105  
    106     //ellipse orientation: 
    107     // the elliptical corss section was transformed and projected 
    108     // into the detector plane already through sin(alpha)and furthermore psi remains as same 
    109     // on the detector plane. 
    110     // So, all we need is to calculate the angle (nu) of the minor axis of the ellipse wrt 
    111     // the wave vector q. 
    112  
    113     //x- y- component on the detector plane. 
    114     const double ella_x =  -cos(_phi)*sin(_psi) * sin(_theta)+sin(_phi)*cos(_psi); 
    115     const double ella_y =  sin(_psi)*cos(_theta); 
    116     const double ellb_x =  -sin(_theta)*cos(_psi)*cos(_phi)-sin(_psi)*sin(_phi); 
    117     const double ellb_y =  cos(_theta)*cos(_psi); 
    118  
    119     // Compute the angle btw vector q and the 
    120     // axis of the cylinder 
    121     double cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z; 
    122  
    123     // calculate the axis of the ellipse wrt q-coord. 
    124     double cos_nu = ella_x*q_x + ella_y*q_y; 
    125     double cos_mu = ellb_x*q_x + ellb_y*q_y; 
    126  
    127     // The following test should always pass 
    128     if (fabs(cos_val)>1.0) { 
    129       //printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    130       cos_val = 1.0; 
    131     } 
    132     if (fabs(cos_nu)>1.0) { 
    133       //printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n"); 
    134       cos_nu = 1.0; 
    135     } 
    136     if (fabs(cos_mu)>1.0) { 
    137       //printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n"); 
    138       cos_mu = 1.0; 
    139     } 
    140  
    141     const double r_major = r_ratio * radius_minor; 
    142     const double qr = q*sqrt( r_major*r_major*cos_nu*cos_nu + radius_minor*radius_minor*cos_mu*cos_mu ); 
    143     const double qL = q*length*cos_val/2.0; 
    144  
    145     double Be; 
    146     if (qr==0){ 
    147       Be = 0.5; 
    148     }else{ 
    149       //Be = NR_BessJ1(qr)/qr; 
    150       Be = 0.5*sas_J1c(qr); 
    151     } 
    152  
    153     double Si; 
    154     if (qL==0){ 
    155       Si = 1.0; 
    156     }else{ 
    157       Si = sin(qL)/qL; 
    158     } 
    159  
    160     const double k = 2.0 * Be * Si; 
     90    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 
     91    // Given:    radius_major = r_ratio * radius_minor 
     92    const double r = radius_minor*sqrt(square(r_ratio*cos_nu) + cos_mu*cos_mu); 
     93    const double be = sas_J1c(q*r); 
     94    const double si = sinc(q*0.5*length*cos_val); 
     95    const double Aq = be * si; 
     96    const double delrho = sld - solvent_sld; 
    16197    const double vol = form_volume(radius_minor, r_ratio, length); 
    162     return (sld - solvent_sld) * (sld - solvent_sld) * k * k *vol*vol*1.0e-4; 
     98    return 1.0e-4 * square(delrho * vol * Aq); 
    16399} 

  • sasmodels/models/fcc_paracrystal.c

    r0bef47b r4962519  
    3232 
    3333        const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 
    34         result = pow((1.0-(temp10*temp10)),3)*temp6 
     34        result = cube(1.0-(temp10*temp10))*temp6 
    3535            / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    3636              * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
    … …  
    8585} 
    8686 
     87double Iqxy(double qx, double qy, 
     88    double dnn, double d_factor, double radius, 
     89    double sld, double solvent_sld, 
     90    double theta, double phi, double psi) 
     91{ 
     92    double q, cos_a1, cos_a2, cos_a3; 
     93    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
    8794 
    88 double Iqxy(double qx, double qy, double dnn, 
    89     double d_factor, double radius,double sld, double solvent_sld, 
    90     double theta, double phi, double psi){ 
     95    const double a1 = cos_a2 + cos_a3; 
     96    const double a2 = cos_a3 + cos_a1; 
     97    const double a3 = cos_a2 + cos_a1; 
     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); 
    91105 
    92   double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z, 
    93   // double q_z; 
    94   double cos_val_b3, cos_val_b2, cos_val_b1; 
    95   double a1_dot_q, a2_dot_q,a3_dot_q; 
    96   double answer; 
    97   double Zq, Fkq, Fkq_2; 
     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); 
    98107 
    99   //convert to q and make scaled values 
    100   double q = sqrt(qx*qx+qy*qy); 
    101   double q_x = qx/q; 
    102   double q_y = qy/q; 
    103  
    104   //convert angle degree to radian 
    105   theta = theta * M_PI_180; 
    106   phi = phi * M_PI_180; 
    107   psi = psi * M_PI_180; 
    108  
    109   const double Da = d_factor*dnn; 
    110   const double s1 = dnn/sqrt(0.75); 
    111  
    112  
    113   //the occupied volume of the lattice 
    114   const double latticescale = 2.0*sphere_volume(radius/s1); 
    115   // q vector 
    116   // q_z = 0.0; // for SANS; assuming qz is negligible 
    117   /// Angles here are respect to detector coordinate 
    118   ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854)) 
    119     // b3 axis orientation 
    120     b3_x = cos(theta) * cos(phi); 
    121     b3_y = sin(theta); 
    122     //b3_z = -cos(theta) * sin(phi); 
    123     cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z; 
    124  
    125     //alpha = acos(cos_val_b3); 
    126     // b1 axis orientation 
    127     b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi); 
    128     b1_y = sin(psi)*cos(theta); 
    129     cos_val_b1 = b1_x*q_x + b1_y*q_y; 
    130     // b2 axis orientation 
    131     b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi); 
    132         b2_y = cos(theta)*cos(psi); 
    133     cos_val_b2 = b2_x*q_x + b2_y*q_y; 
    134  
    135     // The following test should always pass 
    136     if (fabs(cos_val_b3)>1.0) { 
    137       //printf("FCC_ana_2D: Unexpected error: cos()>1\n"); 
    138       cos_val_b3 = 1.0; 
    139     } 
    140     if (fabs(cos_val_b2)>1.0) { 
    141       //printf("FCC_ana_2D: Unexpected error: cos()>1\n"); 
    142       cos_val_b2 = 1.0; 
    143     } 
    144     if (fabs(cos_val_b1)>1.0) { 
    145       //printf("FCC_ana_2D: Unexpected error: cos()>1\n"); 
    146       cos_val_b1 = 1.0; 
    147     } 
    148     // Compute the angle btw vector q and the a3 axis 
    149     a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3); 
    150  
    151     // a1 axis 
    152     a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1); 
    153  
    154     // a2 axis 
    155     a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2); 
    156  
    157  
    158     // Get Fkq and Fkq_2 
    159     Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q)); 
    160     Fkq_2 = Fkq*Fkq; 
    161     // Call Zq=Z1*Z2*Z3 
    162     Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2); 
    163     Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2); 
    164     Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2); 
    165  
    166   // Use SphereForm directly from libigor 
    167   answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale; 
    168  
    169   return answer; 
    170  } 
     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; 
     112} 

  • sasmodels/models/flexible_cylinder_elliptical.c

    rb66d38e r4962519  
    22double Iq(double q, double length, double kuhn_length, double radius, 
    33          double axis_ratio, double sld, double solvent_sld); 
    4 double Iqxy(double qx, double qy, double length, double kuhn_length, 
    5             double radius, double axis_ratio, double sld, double solvent_sld); 
    64double flexible_cylinder_ex_kernel(double q, double length, double kuhn_length, 
    75                                double radius, double axis_ratio, double sld, 
    … …  
    2018 
    2119    for(int i=0;i<N_POINTS_76;i++) { 
    22         double zi = ( Gauss76Z[i] + 1.0 )*M_PI/4.0; 
     20        double zi = ( Gauss76Z[i] + 1.0 )*M_PI_4; 
    2321        double sn, cn; 
    2422        SINCOS(zi, sn, cn); 
    25         double arg = q*sqrt(a*a*sn*sn+b*b*cn*cn); 
    26         double yyy = pow((double)sas_J1c(arg),2); 
    27         yyy *= Gauss76Wt[i]; 
    28         summ += yyy; 
     23        double arg = q*sqrt(a*a*sn*sn + b*b*cn*cn); 
     24        double yyy = sas_J1c(arg); 
     25        summ += Gauss76Wt[i] * yyy * yyy; 
    2926    } 
    3027 
    … …  
    7774} 
    7875 
    79 double Iqxy(double qx, double qy, 
    80             double length, 
    81             double kuhn_length, 
    82             double radius, 
    83             double axis_ratio, 
    84             double sld, 
    85             double solvent_sld) 
    86 { 
    87     double q; 
    88     q = sqrt(qx*qx+qy*qy); 
    89     double result = flexible_cylinder_ex_kernel(q, 
    90                     length, 
    91                     kuhn_length, 
    92                     radius, 
    93                     axis_ratio, 
    94                     sld, 
    95                     solvent_sld); 
    96  
    97     return result; 
    98 } 

  • sasmodels/models/fractal_core_shell.c

    ra807206 r3a48772  
    1313double form_volume(double radius, double thickness) 
    1414{ 
    15     return 4.0 * M_PI / 3.0 * pow((radius + thickness), 3); 
     15    return M_4PI_3 * cube(radius + thickness); 
    1616} 
    1717 

  • sasmodels/models/fractal_core_shell.py

    ra807206 r4962519  
    100100        @param thickness: shell thickness 
    101101    """ 
    102     whole = 4.0 * pi / 3.0 * pow((radius + thickness), 3) 
    103     core = 4.0 * pi / 3.0 * radius * radius * radius 
     102    whole = 4.0/3.0 * pi * (radius + thickness)**3 
     103    core = 4.0/3.0 * pi * radius**3 
    104104    return whole, whole-core 
    105105 

  • sasmodels/models/fuzzy_sphere.py

    r9a4811a r3a48772  
    8686# This should perhaps be volume normalized? 
    8787form_volume = """ 
    88     return 1.333333333333333*M_PI*radius*radius*radius; 
     88    return M_4PI_3*cube(radius); 
    8989    """ 
    9090 

  • sasmodels/models/hayter_msa.c

    r0bef47b r4962519  
    2323        double SIdiam, diam, Kappa, cs, IonSt; 
    2424        double  Perm, Beta; 
    25         double pi, charge; 
     25        double charge; 
    2626        int ierr; 
    2727         
    28         pi = M_PI; 
    29  
    3028        diam=2*radius_effective;                //in A 
    3129 
    … …  
    3836        charge=zz*Elcharge;             //in Coulomb (C) 
    3937        SIdiam = diam*1.0E-10;          //in m 
    40         Vp=4.0*pi/3.0*(SIdiam/2.0)*(SIdiam/2.0)*(SIdiam/2.0);   //in m^3 
     38        Vp=M_4PI_3*cube(SIdiam/2.0);    //in m^3 
    4139        cs=csalt*6.022E23*1.0E3;        //# salt molecules/m^3 
    4240         
    … …  
    5048        Kappa=sqrt(2*Beta*IonSt/Perm);     //Kappa calc from Ionic strength 
    5149                                                                           //   Kappa=2/SIdiam                                  // Use to compare with HP paper 
    52         gMSAWave[5]=Beta*charge*charge/(pi*Perm*SIdiam*pow((2.0+Kappa*SIdiam),2)); 
     50        gMSAWave[5]=Beta*charge*charge/(M_PI*Perm*SIdiam*square(2.0+Kappa*SIdiam)); 
    5351         
    5452        //         Finally set up dimensionless parameters  

  • sasmodels/models/hollow_cylinder.c

    r0d6e865 r5bddd89  
    11double form_volume(double radius, double thickness, double length); 
    2  
    32double Iq(double q, double radius, double thickness, double length, double sld, 
    43        double solvent_sld); 
    … …  
    98 
    109// From Igor library 
    11 static double hollow_cylinder_scaling( 
    12     double integrand, double delrho, double volume) 
     10static double 
     11_hollow_cylinder_scaling(double integrand, double delrho, double volume) 
    1312{ 
    14     double answer; 
    15     // Multiply by contrast^2 
    16     answer = integrand*delrho*delrho; 
    17  
    18     //normalize by cylinder volume 
    19     answer *= volume*volume; 
    20  
    21     //convert to [cm-1] 
    22     answer *= 1.0e-4; 
    23  
    24     return answer; 
     13    return 1.0e-4 * square(volume * delrho * integrand); 
    2514} 
    2615 
    2716 
    28 static double _hollow_cylinder_kernel( 
    29     double q, double radius, double thickness, double length, double dum) 
     17static double 
     18_hollow_cylinder_kernel(double q, 
     19    double radius, double thickness, double length, double sin_val, double cos_val) 
    3020{ 
    31     const double qs = q*sqrt(1.0-dum*dum); 
     21    const double qs = q*sin_val; 
    3222    const double lam1 = sas_J1c((radius+thickness)*qs); 
    3323    const double lam2 = sas_J1c(radius*qs); 
    … …  
    3525    //Note: lim_{r -> r_c} psi = J0(radius_core*qs) 
    3626    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); //SRK 10/19/00 
    37     const double t2 = sinc(q*length*dum/2.0); 
    38     return square(psi*t2); 
     27    const double t2 = sinc(0.5*q*length*cos_val); 
     28    return psi*t2; 
    3929} 
    4030 
    41  
    42 static double hollow_cylinder_analytical_2D_scaled( 
    43     double q, double q_x, double q_y, double radius, double thickness, 
    44     double length, double sld, double solvent_sld, double theta, double phi) 
    45 { 
    46     double cyl_x, cyl_y; //, cyl_z 
    47     //double q_z; 
    48     double vol, cos_val, delrho; 
    49     double answer; 
    50     //convert angle degree to radian 
    51     theta = theta * M_PI_180; 
    52     phi = phi * M_PI_180; 
    53     delrho = solvent_sld - sld; 
    54  
    55     // Cylinder orientation 
    56     cyl_x = sin(theta) * cos(phi); 
    57     cyl_y = sin(theta) * sin(phi); 
    58     //cyl_z = -cos(theta) * sin(phi); 
    59  
    60     // q vector 
    61     //q_z = 0; 
    62  
    63     // Compute the angle btw vector q and the 
    64     // axis of the cylinder 
    65     cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z; 
    66  
    67     answer = _hollow_cylinder_kernel(q, radius, thickness, length, cos_val); 
    68  
    69     vol = form_volume(radius, thickness, length); 
    70     answer = hollow_cylinder_scaling(answer, delrho, vol); 
    71  
    72     return answer; 
    73 } 
    74  
    75  
    76 double form_volume(double radius, double thickness, double length) 
     31double 
     32form_volume(double radius, double thickness, double length) 
    7733{ 
    7834    double v_shell = M_PI*length*((radius+thickness)*(radius+thickness)-radius*radius); 
    … …  
    8137 
    8238 
    83 double Iq(double q, double radius, double thickness, double length, 
     39double 
     40Iq(double q, double radius, double thickness, double length, 
    8441    double sld, double solvent_sld) 
    8542{ 
    86     int i; 
    87     double lower,upper,zi, inter;               //upper and lower integration limits 
    88     double summ,answer,delrho;                  //running tally of integration 
    89     double norm,volume; //final calculation variables 
     43    const double lower = 0.0; 
     44    const double upper = 1.0;           //limits of numerical integral 
    9045 
    91     lower = 0.0; 
    92     upper = 1.0;                //limits of numerical integral 
    93  
    94     summ = 0.0;                 //initialize intergral 
    95     for (i=0;i<76;i++) { 
    96         zi = ( Gauss76Z[i] * (upper-lower) + lower + upper )/2.0; 
    97         inter = Gauss76Wt[i] * _hollow_cylinder_kernel(q, radius, thickness, length, zi); 
    98         summ += inter; 
     46    double summ = 0.0;                  //initialize intergral 
     47    for (int i=0;i<76;i++) { 
     48        const double cos_val = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
     49        const double sin_val = sqrt(1.0 - cos_val*cos_val); 
     50        const double inter = _hollow_cylinder_kernel(q, radius, thickness, length, 
     51                                                     sin_val, cos_val); 
     52        summ += Gauss76Wt[i] * inter; 
    9953    } 
    10054 
    101     norm = summ*(upper-lower)/2.0; 
    102     volume = form_volume(radius, thickness, length); 
    103     delrho = solvent_sld - sld; 
    104     answer = hollow_cylinder_scaling(norm, delrho, volume); 
     55    const double Aq = 0.5*summ*(upper-lower); 
     56    const double volume = form_volume(radius, thickness, length); 
     57    return _hollow_cylinder_scaling(Aq, solvent_sld - sld, volume); 
     58} 
    10559 
    106     return(answer); 
     60double 
     61Iqxy(double qx, double qy, 
     62    double radius, double thickness, double length, 
     63    double sld, double solvent_sld, double theta, double phi) 
     64{ 
     65    double q, sin_alpha, cos_alpha; 
     66    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     67    const double Aq = _hollow_cylinder_kernel(q, radius, thickness, length, 
     68        sin_alpha, cos_alpha); 
     69 
     70    const double vol = form_volume(radius, thickness, length); 
     71    return _hollow_cylinder_scaling(Aq, solvent_sld-sld, vol); 
    10772} 
    10873 
    10974 
    110 double Iqxy(double qx, double qy, double radius, double thickness, 
    111     double length, double sld, double solvent_sld, double theta, double phi) 
    112 { 
    113     const double q = sqrt(qx*qx+qy*qy); 
    114     return hollow_cylinder_analytical_2D_scaled(q, qx/q, qy/q, radius, thickness, length, sld, solvent_sld, theta, phi); 
    115 } 

  • sasmodels/models/hollow_rectangular_prism.c

    ra807206 rab2aea8  
    22double Iq(double q, double sld, double solvent_sld, double length_a,  
    33          double b2a_ratio, double c2a_ratio, double thickness); 
    4 double Iqxy(double qx, double qy, double sld, double solvent_sld,  
    5             double length_a, double b2a_ratio, double c2a_ratio, double thickness); 
    64 
    75double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
    86{ 
    9     double b_side = length_a * b2a_ratio; 
    10     double c_side = length_a * c2a_ratio; 
     7    double length_b = length_a * b2a_ratio; 
     8    double length_c = length_a * c2a_ratio; 
    119    double a_core = length_a - 2.0*thickness; 
    12     double b_core = b_side - 2.0*thickness; 
    13     double c_core = c_side - 2.0*thickness; 
     10    double b_core = length_b - 2.0*thickness; 
     11    double c_core = length_c - 2.0*thickness; 
    1412    double vol_core = a_core * b_core * c_core; 
    15     double vol_total = length_a * b_side * c_side; 
     13    double vol_total = length_a * length_b * length_c; 
    1614    double vol_shell = vol_total - vol_core; 
    1715    return vol_shell; 
    … …  
    2826    double termA1, termA2, termB1, termB2, termC1, termC2; 
    2927     
    30     double b_side = length_a * b2a_ratio; 
    31     double c_side = length_a * c2a_ratio; 
     28    double length_b = length_a * b2a_ratio; 
     29    double length_c = length_a * c2a_ratio; 
    3230    double a_half = 0.5 * length_a; 
    33     double b_half = 0.5 * b_side; 
    34     double c_half = 0.5 * c_side; 
     31    double b_half = 0.5 * length_b; 
     32    double c_half = 0.5 * length_c; 
     33    double vol_total = length_a * length_b * length_c; 
     34    double vol_core = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness); 
    3535 
    36    //Integration limits to use in Gaussian quadrature 
     36    //Integration limits to use in Gaussian quadrature 
    3737    double v1a = 0.0; 
    38     double v1b = 0.5 * M_PI;  //theta integration limits 
     38    double v1b = M_PI_2;  //theta integration limits 
    3939    double v2a = 0.0; 
    40     double v2b = 0.5 * M_PI;  //phi integration limits 
     40    double v2b = M_PI_2;  //phi integration limits 
    4141     
    42     //Order of integration 
    43     int nordi=76;                                
    44     int nordj=76; 
     42    double outer_sum = 0.0; 
     43     
     44    for(int i=0; i<76; i++) { 
    4545 
    46     double sumi = 0.0; 
    47      
    48     for(int i=0; i<nordi; i++) { 
     46        double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );     
    4947 
    50             double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );  
     48        double termC1 = sinc(q * c_half * cos(theta)); 
     49        double termC2 = sinc(q * (c_half-thickness)*cos(theta)); 
    5150 
    52             double arg = q * c_half * cos(theta); 
    53             if (fabs(arg) > 1.e-16) {termC1 = sin(arg)/arg;} else {termC1 = 1.0;} 
    54             arg = q * (c_half-thickness)*cos(theta); 
    55             if (fabs(arg) > 1.e-16) {termC2 = sin(arg)/arg;} else {termC2 = 1.0;} 
    56  
    57             double sumj = 0.0; 
     51        double inner_sum = 0.0; 
    5852         
    59             for(int j=0; j<nordj; j++) { 
     53        for(int j=0; j<76; j++) { 
    6054 
    6155            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b );  
    … …  
    6357            // Amplitude AP from eqn. (13), rewritten to avoid round-off effects when arg=0 
    6458 
    65                 arg = q * a_half * sin(theta) * sin(phi); 
    66                 if (fabs(arg) > 1.e-16) {termA1 = sin(arg)/arg;} else {termA1 = 1.0;} 
    67                 arg = q * (a_half-thickness) * sin(theta) * sin(phi); 
    68                 if (fabs(arg) > 1.e-16) {termA2 = sin(arg)/arg;} else {termA2 = 1.0;} 
     59            termA1 = sinc(q * a_half * sin(theta) * sin(phi)); 
     60            termA2 = sinc(q * (a_half-thickness) * sin(theta) * sin(phi)); 
    6961 
    70                 arg = q * b_half * sin(theta) * cos(phi); 
    71                 if (fabs(arg) > 1.e-16) {termB1 = sin(arg)/arg;} else {termB1 = 1.0;} 
    72                 arg = q * (b_half-thickness) * sin(theta) * cos(phi); 
    73                 if (fabs(arg) > 1.e-16) {termB2 = sin(arg)/arg;} else {termB2 = 1.0;} 
     62            termB1 = sinc(q * b_half * sin(theta) * cos(phi)); 
     63            termB2 = sinc(q * (b_half-thickness) * sin(theta) * cos(phi)); 
    7464 
    75             double AP1 = (length_a*b_side*c_side) * termA1 * termB1 * termC1; 
    76             double AP2 = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness) * termA2 * termB2 * termC2; 
    77             double AP = AP1 - AP2; 
     65            double AP1 = vol_total * termA1 * termB1 * termC1; 
     66            double AP2 = vol_core * termA2 * termB2 * termC2; 
    7867 
    79                 sumj += Gauss76Wt[j] * (AP*AP); 
     68            inner_sum += Gauss76Wt[j] * square(AP1-AP2); 
    8069 
    81             } 
     70        } 
    8271 
    83             sumj = 0.5 * (v2b-v2a) * sumj; 
    84             sumi += Gauss76Wt[i] * sumj * sin(theta); 
     72        inner_sum = 0.5 * (v2b-v2a) * inner_sum; 
     73        outer_sum += Gauss76Wt[i] * inner_sum * sin(theta); 
    8574 
    8675    } 
    8776 
    88     double answer = 0.5*(v1b-v1a)*sumi; 
     77    double answer = 0.5*(v1b-v1a)*outer_sum; 
    8978 
    9079    // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization) 
    9180    // The factor 2 is due to the different theta integration limit (pi/2 instead of pi) 
    92     answer *= (2.0/M_PI); 
     81    answer /= M_PI_2; 
    9382 
    9483    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 
    95     answer *= (sld-solvent_sld)*(sld-solvent_sld); 
     84    answer *= square(sld-solvent_sld); 
    9685 
    9786    // Convert from [1e-12 A-1] to [cm-1] 
    … …  
    10190     
    10291} 
    103  
    104 double Iqxy(double qx, double qy, 
    105     double sld, 
    106     double solvent_sld, 
    107     double length_a, 
    108     double b2a_ratio, 
    109     double c2a_ratio, 
    110     double thickness) 
    111 { 
    112     double q = sqrt(qx*qx + qy*qy); 
    113     double intensity = Iq(q, sld, solvent_sld, length_a, b2a_ratio, c2a_ratio, thickness);  
    114     return intensity;     
    115 } 

  • sasmodels/models/hollow_rectangular_prism.py

    ra807206 rab2aea8  
    1313the difference of the amplitudes of two massive parallelepipeds 
    1414differing in their outermost dimensions in each direction by the 
    15 same length increment :math:`2\Delta` (Nayuk, 2012). 
     15same length increment $2\Delta$ (Nayuk, 2012). 
    1616 
    1717As in the case of the massive parallelepiped model (:ref:`rectangular-prism`), 
    … …  
    5858 
    5959.. math:: 
    60   I(q) = \text{scale} \times V \times (\rho_{\text{p}} - 
    61   \rho_{\text{solvent}})^2 \times P(q) + \text{background} 
     60  I(q) = \text{scale} \times V \times (\rho_\text{p} - 
     61  \rho_\text{solvent})^2 \times P(q) + \text{background} 
    6262 
    63 where $\rho_{\text{p}}$ is the scattering length of the parallelepiped, 
    64 $\rho_{\text{solvent}}$ is the scattering length of the solvent, 
     63where $\rho_\text{p}$ is the scattering length of the parallelepiped, 
     64$\rho_\text{solvent}$ is the scattering length of the solvent, 
    6565and (if the data are in absolute units) *scale* represents the volume fraction 
    6666(which is unitless). 
    … …  
    148148# parameters for demo 
    149149demo = dict(scale=1, background=0, 
    150             sld=6.3e-6, sld_solvent=1.0e-6, 
     150            sld=6.3, sld_solvent=1.0, 
    151151            length_a=35, b2a_ratio=1, c2a_ratio=1, thickness=1, 
    152152            length_a_pd=0.1, length_a_pd_n=10, 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    ra807206 rab2aea8  
    22double Iq(double q, double sld, double solvent_sld, double length_a,  
    33          double b2a_ratio, double c2a_ratio); 
    4 double Iqxy(double qx, double qy, double sld, double solvent_sld,  
    5             double length_a, double b2a_ratio, double c2a_ratio); 
    64 
    75double form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
    86{ 
    9     double b_side = length_a * b2a_ratio; 
    10     double c_side = length_a * c2a_ratio; 
    11     double vol_shell = 2.0 * (length_a*b_side + length_a*c_side + b_side*c_side); 
     7    double length_b = length_a * b2a_ratio; 
     8    double length_c = length_a * c2a_ratio; 
     9    double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c); 
    1210    return vol_shell; 
    1311} 
    … …  
    2018    double c2a_ratio) 
    2119{ 
    22     double b_side = length_a * b2a_ratio; 
    23     double c_side = length_a * c2a_ratio; 
    24     double a_half = 0.5 * length_a; 
    25     double b_half = 0.5 * b_side; 
    26     double c_half = 0.5 * c_side; 
     20    const double length_b = length_a * b2a_ratio; 
     21    const double length_c = length_a * c2a_ratio; 
     22    const double a_half = 0.5 * length_a; 
     23    const double b_half = 0.5 * length_b; 
     24    const double c_half = 0.5 * length_c; 
    2725 
    2826   //Integration limits to use in Gaussian quadrature 
    29     double v1a = 0.0; 
    30     double v1b = 0.5 * M_PI;  //theta integration limits 
    31     double v2a = 0.0; 
    32     double v2b = 0.5 * M_PI;  //phi integration limits 
     27    const double v1a = 0.0; 
     28    const double v1b = M_PI_2;  //theta integration limits 
     29    const double v2a = 0.0; 
     30    const double v2b = M_PI_2;  //phi integration limits 
    3331     
    34     //Order of integration 
    35     int nordi=76;                                
    36     int nordj=76; 
     32    double outer_sum = 0.0; 
     33    for(int i=0; i<76; i++) { 
     34        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
    3735 
    38     double sumi = 0.0; 
    39      
    40     for(int i=0; i<nordi; i++) { 
     36        double sin_theta, cos_theta; 
     37        double sin_c, cos_c; 
     38        SINCOS(theta, sin_theta, cos_theta); 
     39        SINCOS(q*c_half*cos_theta, sin_c, cos_c); 
    4140 
    42             double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );  
    43          
    4441        // To check potential problems if denominator goes to zero here !!! 
    45         double termAL_theta = 8.0*cos(q*c_half*cos(theta)) / (q*q*sin(theta)*sin(theta)); 
    46         double termAT_theta = 8.0*sin(q*c_half*cos(theta)) / (q*q*sin(theta)*cos(theta)); 
     42        const double termAL_theta = 8.0 * cos_c / (q*q*sin_theta*sin_theta); 
     43        const double termAT_theta = 8.0 * sin_c / (q*q*sin_theta*cos_theta); 
    4744 
    48             double sumj = 0.0; 
    49          
    50             for(int j=0; j<nordj; j++) { 
     45        double inner_sum = 0.0; 
     46        for(int j=0; j<76; j++) { 
     47            const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
    5148 
    52             double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b );  
    53              
     49            double sin_phi, cos_phi; 
     50            double sin_a, cos_a; 
     51            double sin_b, cos_b; 
     52            SINCOS(phi, sin_phi, cos_phi); 
     53            SINCOS(q*a_half*sin_theta*sin_phi, sin_a, cos_a); 
     54            SINCOS(q*b_half*sin_theta*cos_phi, sin_b, cos_b); 
     55 
    5456            // Amplitude AL from eqn. (7c) 
    55             double AL = termAL_theta * sin(q*a_half*sin(theta)*sin(phi)) *  
    56                 sin(q*b_half*sin(theta)*cos(phi)) / (sin(phi)*cos(phi)); 
     57            const double AL = termAL_theta 
     58                * sin_a*sin_b / (sin_phi*cos_phi); 
    5759 
    5860            // Amplitude AT from eqn. (9) 
    59             double AT = termAT_theta * (  (cos(q*a_half*sin(theta)*sin(phi))*sin(q*b_half*sin(theta)*cos(phi))/cos(phi))  
    60                 + (cos(q*b_half*sin(theta)*cos(phi))*sin(q*a_half*sin(theta)*sin(phi))/sin(phi)) ); 
     61            const double AT = termAT_theta 
     62                * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi ); 
    6163 
    62             sumj += Gauss76Wt[j] * (AL+AT)*(AL+AT); 
     64            inner_sum += Gauss76Wt[j] * square(AL+AT); 
     65        } 
    6366 
    64             } 
    65  
    66             sumj = 0.5 * (v2b-v2a) * sumj; 
    67             sumi += Gauss76Wt[i] * sumj * sin(theta); 
    68  
     67        inner_sum *= 0.5 * (v2b-v2a); 
     68        outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 
    6969    } 
    7070 
    71     double answer = 0.5*(v1b-v1a)*sumi; 
     71    outer_sum *= 0.5*(v1b-v1a); 
    7272 
    7373    // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization) 
    7474    // The factor 2 is due to the different theta integration limit (pi/2 instead of pi) 
    75     answer *= (2.0/M_PI); 
     75    double answer = outer_sum/M_PI_2; 
    7676 
    7777    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 
    78     answer *= (sld-solvent_sld)*(sld-solvent_sld); 
     78    answer *= square(sld-solvent_sld); 
    7979 
    8080    // Convert from [1e-12 A-1] to [cm-1] 
    … …  
    8282 
    8383    return answer; 
    84      
    8584} 
    86  
    87 double Iqxy(double qx, double qy, 
    88     double sld, 
    89     double solvent_sld, 
    90     double length_a, 
    91     double b2a_ratio, 
    92     double c2a_ratio) 
    93 { 
    94     double q = sqrt(qx*qx + qy*qy); 
    95     double intensity = Iq(q, sld, solvent_sld, length_a, b2a_ratio, c2a_ratio);  
    96     return intensity;     
    97 } 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.py

    ra807206 rab2aea8  
    33r""" 
    44 
    5 This model provides the form factor, *P(q)*, for a hollow rectangular 
     5This model provides the form factor, $P(q)$, for a hollow rectangular 
    66prism with infinitely thin walls. It computes only the 1D scattering, not the 2D. 
    77 
    … …  
    1414 
    1515Assuming a hollow parallelepiped with infinitely thin walls, edge lengths 
    16 :math:`A \le B \le C` and presenting an orientation with respect to the 
    17 scattering vector given by |theta| and |phi|, where |theta| is the angle 
    18 between the *z* axis and the longest axis of the parallelepiped *C*, and 
    19 |phi| is the angle between the scattering vector (lying in the *xy* plane) 
    20 and the *y* axis, the form factor is given by 
     16$A \le B \le C$ and presenting an orientation with respect to the 
     17scattering vector given by $\theta$ and $\phi$, where $\theta$ is the angle 
     18between the $z$ axis and the longest axis of the parallelepiped $C$, and 
     19$\phi$ is the angle between the scattering vector (lying in the $xy$ plane) 
     20and the $y$ axis, the form factor is given by 
    2121 
    2222.. math:: 
    23   P(q) =  \frac{1}{V^2} \frac{2}{\pi} \int_0^{\frac{\pi}{2}} 
    24   \int_0^{\frac{\pi}{2}} [A_L(q)+A_T(q)]^2 \sin\theta d\theta d\phi 
     23 
     24    P(q) = \frac{1}{V^2} \frac{2}{\pi} \int_0^{\frac{\pi}{2}} 
     25           \int_0^{\frac{\pi}{2}} [A_L(q)+A_T(q)]^2 \sin\theta\,d\theta\,d\phi 
    2526 
    2627where 
    2728 
    2829.. math:: 
    29   V = 2AB + 2AC + 2BC 
    3030 
    31 .. math:: 
    32   A_L(q) =  8 \times \frac{ \sin \bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr) 
    33                               \sin \bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr) 
    34                               \cos \bigl( q \frac{C}{2} \cos\theta \bigr) } 
    35                             {q^2 \, \sin^2\theta \, \sin\phi \cos\phi} 
    36  
    37 .. math:: 
    38   A_T(q) =  A_F(q) \times \frac{2 \, \sin \bigl( q \frac{C}{2} \cos\theta \bigr)}{q \, \cos\theta} 
     31    V &= 2AB + 2AC + 2BC \\ 
     32    A_L(q) &=  8 \times \frac{ 
     33            \sin \left( \tfrac{1}{2} q A \sin\phi \sin\theta \right) 
     34            \sin \left( \tfrac{1}{2} q B \cos\phi \sin\theta \right) 
     35            \cos \left( \tfrac{1}{2} q C \cos\theta \right) 
     36        }{q^2 \, \sin^2\theta \, \sin\phi \cos\phi} \\ 
     37    A_T(q) &=  A_F(q) \times 
     38      \frac{2\,\sin \left( \tfrac{1}{2} q C \cos\theta \right)}{q\,\cos\theta} 
    3939 
    4040and 
    4141 
    4242.. math:: 
    43   A_F(q) =  4 \frac{ \cos \bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr) 
    44                        \sin \bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr) } 
     43 
     44  A_F(q) =  4 \frac{ \cos \left( \tfrac{1}{2} q A \sin\phi \sin\theta \right) 
     45                       \sin \left( \tfrac{1}{2} q B \cos\phi \sin\theta \right) } 
    4546                     {q \, \cos\phi \, \sin\theta} + 
    46               4 \frac{ \sin \bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr) 
    47                        \cos \bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr) } 
     47              4 \frac{ \sin \left( \tfrac{1}{2} q A \sin\phi \sin\theta \right) 
     48                       \cos \left( \tfrac{1}{2} q B \cos\phi \sin\theta \right) } 
    4849                     {q \, \sin\phi \, \sin\theta} 
    4950 
    … …  
    5152 
    5253.. math:: 
    53   I(q) = \mbox{scale} \times V \times (\rho_{\mbox{p}} - \rho_{\mbox{solvent}})^2 \times P(q) 
    5454 
    55 where *V* is the volume of the rectangular prism, :math:`\rho_{\mbox{p}}` 
    56 is the scattering length of the parallelepiped, :math:`\rho_{\mbox{solvent}}` 
     55  I(q) = \text{scale} \times V \times (\rho_\text{p} - \rho_\text{solvent})^2 \times P(q) 
     56 
     57where $V$ is the volume of the rectangular prism, $\rho_\text{p}$ 
     58is the scattering length of the parallelepiped, $\rho_\text{solvent}$ 
    5759is the scattering length of the solvent, and (if the data are in absolute 
    5860units) *scale* represents the volume fraction (which is unitless). 
    … …  
    127129# parameters for demo 
    128130demo = dict(scale=1, background=0, 
    129             sld=6.3e-6, sld_solvent=1.0e-6, 
     131            sld=6.3, sld_solvent=1.0, 
    130132            length_a=35, b2a_ratio=1, c2a_ratio=1, 
    131133            length_a_pd=0.1, length_a_pd_n=10, 

  • sasmodels/models/lamellar_stack_paracrystal.c

    r0bef47b r4962519  
    6767        double Snq; 
    6868 
    69         Snq = an/( (double)Nlayers*pow((1.0+ww*ww-2.0*ww*cos(qval*davg)),2) ); 
     69        Snq = an/( (double)Nlayers*square(1.0+ww*ww-2.0*ww*cos(qval*davg)) ); 
    7070         
    7171        return(Snq); 

  • sasmodels/models/lib/Si.c

    r0278e3f ra5b6997  
    33double Si(double x) 
    44{ 
    5     double out; 
     5    if (x >= M_PI*6.2/4.0){ 
     6        const double z = 1./(x*x); 
     7        // Explicitly writing factorial values triples the speed of the calculation 
     8        const double out_cos = (((-720.*z + 24.)*z - 2.)*z + 1.)/x; 
     9        const double out_sin = (((-5040.*z + 120.)*z - 6.)*z + 1)*z; 
    610 
    7     if (x >= M_PI*6.2/4.0){ 
    8         double out_sin = 0.0; 
    9         double out_cos = 0.0; 
    10         out = M_PI/2.0; 
    11  
     11        double cos_x, sin_x; 
     12        SINCOS(x, cos_x, sin_x); 
     13        return M_PI_2 - cos_x*out_cos - sin_x*out_sin; 
     14    } else { 
     15        const double z = x*x; 
    1216        // Explicitly writing factorial values triples the speed of the calculation 
    13         out_cos = 1./x - 2./pow(x,3) + 24./pow(x,5) - 720./pow(x,7); 
    14         out_sin = 1./pow(x,2) - 6./pow(x,4) + 120./pow(x,6) - 5040./pow(x,8); 
    15  
    16         out -= cos(x) * out_cos; 
    17         out -= sin(x) * out_sin; 
    18         return out; 
     17        return (((((-1./439084800.*z 
     18            + 1./3265920.)*z 
     19            - 1./35280.)*z 
     20            + 1./600.)*z 
     21            - 1./18.)*z 
     22            + 1.)*x; 
    1923    } 
    20  
    21     // Explicitly writing factorial values triples the speed of the calculation 
    22     out = x - pow(x, 3)/18. + pow(x,5)/600. - pow(x,7)/35280. + pow(x,9)/3265920. - pow(x,11)/439084800.; 
    23  
    24     return out; 
    2524} 

  • sasmodels/models/lib/core_shell.c

    r7d4b2ae r3a48772  
    1717    const double core_contrast = core_sld - shell_sld; 
    1818    const double core_bes = sph_j1c(core_qr); 
    19     const double core_volume = 4.0 * M_PI / 3.0 * radius * radius * radius; 
     19    const double core_volume = M_4PI_3 * cube(radius); 
    2020    double f = core_volume * core_bes * core_contrast; 
    2121 
    … …  
    2424    const double shell_contrast = shell_sld - solvent_sld; 
    2525    const double shell_bes = sph_j1c(shell_qr); 
    26     const double shell_volume = 4.0 * M_PI / 3.0 * pow((radius + thickness), 3); 
     26    const double shell_volume = M_4PI_3 * cube(radius + thickness); 
    2727    f += shell_volume * shell_bes * shell_contrast; 
    2828    return f * f * 1.0e-4; 

  • sasmodels/models/lib/gfn.c

    r177c1a1 r3a48772  
    1212{ 
    1313    // local variables 
    14     const double pi43=4.0/3.0*M_PI; 
    1514    const double aa = crmaj; 
    1615    const double bb = crmin; 
    … …  
    2019    //const double siq = (uq == 0.0 ? 1.0 : 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq); 
    2120    const double siq = sph_j1c(uq); 
    22     const double vc = pi43*aa*aa*bb; 
     21    const double vc = M_4PI_3*aa*aa*bb; 
    2322    const double gfnc = siq*vc*delpc; 
    2423 
    2524    const double ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx)); 
    2625    const double ut= sqrt(ut2)*qq; 
    27     const double vt = pi43*trmaj*trmaj*trmin; 
     26    const double vt = M_4PI_3*trmaj*trmaj*trmin; 
    2827    //const double sit = (ut == 0.0 ? 1.0 : 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut); 
    2928    const double sit = sph_j1c(ut); 
    … …  
    3332    const double result = tgfn*tgfn; 
    3433 
    35     return (result); 
     34    return result; 
    3635} 

  • sasmodels/models/lib/wrc_cyl.c

    rba32cdd r4962519  
    1414    //return t; 
    1515 
    16     return pow( (1.0 + (x/3.12)*(x/3.12) + 
    17          (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); 
     16    return pow(1.0+square(x/3.12)+cube(x/8.67), 0.176/3.0); 
    1817} 
    1918 
    … …  
    2322{ 
    2423    const double r = b/L; 
    25     return (L*b/6.0) * 
    26            (1.0 - r*1.5  + 1.5*r*r - 0.75*r*r*r*(1.0 - exp(-2.0/r))); 
     24    return (L*b/6.0) * (1.0 - r*(1.5 + r*(1.5 + r*0.75*expm1(-2.0/r)))); 
    2725} 
    2826 
    … …  
    4139} 
    4240 
    43 static inline double 
     41static double 
    4442sech_WR(double x) 
    4543{ 
    … …  
    5149{ 
    5250    double C; 
    53     const double onehalf = 1.0/2.0; 
    5451 
    5552    if( L/b > 10.0) { 
    … …  
    8683 
    8784    const double t2 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    88          Rg02/b2))*((1.0 + onehalf*(((-1.0) - 
     85         Rg02/b2))*((1.0 + 0.5*(((-1.0) - 
    8986         tanh((-C4 + Rgb/C5))))))); 
    9087 
    … …  
    112109 
    113110    const double t9 = (2.0*b4*(((2.0*q0*Rg2)/b - 
    114          (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))*((1.0 + onehalf*(((-1.0) - 
     111         (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))*((1.0 + 0.5*(((-1.0) - 
    115112         tanh(((-C4) + Rgb)/C5)))))); 
    116113 
    117114    const double t10 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    118          Rg02/b2))*((1.0 + onehalf*(((-1.0) - tanh(((-C4) + 
     115         Rg02/b2))*((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 
    119116         Rgb)/C5)))))); 
    120117 
    … …  
    131128 
    132129    const double t14 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    133           Rg02/b2))*((1.0 + onehalf*(((-1.0) - tanh(((-C4) + 
     130          Rg02/b2))*((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 
    134131          Rgb)/C5)))))); 
    135132 
    … …  
    140137 
    141138    double yy = (pow(q0,p1)*(((-((b*M_PI)/(L*q0))) +t1/L +t2/(q04*Rg22) + 
    142         onehalf*t3*t4)) + (t5*((pow(q0,(p1 - p2))* 
     139        0.5*t3*t4)) + (t5*((pow(q0,(p1 - p2))* 
    143140        (((-pow(q0,(-p1)))*(((b2*M_PI)/(L*q02) +t6/L +t7/(2.0*C5) - 
    144         t8/(C5*q04*Rg22) + t9/(q04*Rg22) -t10/(q05*Rg22) + onehalf*t11*t12)) - 
     141        t8/(C5*q04*Rg22) + t9/(q04*Rg22) -t10/(q05*Rg22) + 0.5*t11*t12)) - 
    145142        b*p1*pow(q0,((-1.0) - p1))*(((-((b*M_PI)/(L*q0))) + t13/L + 
    146         t14/(q04*Rg22) + onehalf*t15*((1.0 + tanh(((-C4) + 
     143        t14/(q04*Rg22) + 0.5*t15*((1.0 + tanh(((-C4) + 
    147144        Rgb)/C5))))))))))); 
    148145 
    … …  
    154151{ 
    155152    double C; 
    156     const double onehalf = 1.0/2.0; 
    157153 
    158154    if( L/b > 10.0) { 
    … …  
    201197    const double t5 = (2.0*b4*(((2.0*q0*Rg2)/b - 
    202198         (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))* 
    203          ((1.0 + onehalf*(((-1.0) - tanh(((-C4) + 
     199         ((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 
    204200         Rgb)/C5))))))/(q04*Rg22); 
    205201 
    206202    const double t6 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    207          Rg02/b2))*((1.0 + onehalf*(((-1) - tanh(((-C4) + 
     203         Rg02/b2))*((1.0 + 0.5*(((-1) - tanh(((-C4) + 
    208204         Rgb)/C5))))))/(q05*Rg22); 
    209205 
    … …  
    219215 
    220216    const double t10 = (2.0*b4*(((-1) + pow((double)M_E,(-(Rg02/b2))) + 
    221           Rg02/b2))*((1.0 + onehalf*(((-1) - tanh(((-C4) + 
     217          Rg02/b2))*((1.0 + 0.5*(((-1) - tanh(((-C4) + 
    222218          Rgb)/C5))))))/(q04*Rg22); 
    223219 
    224220    const double yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*M_PI)/(L*q02) + 
    225          t2 + t3 - t4 + t5 - t6 + onehalf*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))* 
     221         t2 + t3 - t4 + t5 - t6 + 0.5*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))* 
    226222         (((-((b*M_PI)/(L*q0))) + t9 + t10 + 
    227          onehalf*((C3*pow(((Rgb)),((-3.0)/miu)) + 
     223         0.5*((C3*pow(((Rgb)),((-3.0)/miu)) + 
    228224         C2*pow(((Rgb)),((-2.0)/miu)) + 
    229225         C1*pow(((Rgb)),((-1.0)/miu))))* 

  • sasmodels/models/linear_pearls.c

    r2c74c11 r4962519  
    1919{ 
    2020    // Pearl volume 
    21     double pearl_vol = 4.0 /3.0 * M_PI * pow(radius, 3.0); 
     21    double pearl_vol = M_4PI_3 * cube(radius); 
    2222    // Return total volume 
    2323    return num_pearls * pearl_vol;; 
    … …  
    3535    double contrast_pearl = pearl_sld - solvent_sld; 
    3636    //each volume 
    37     double pearl_vol = 4.0 /3.0 * M_PI * pow(radius, 3.0); 
     37    double pearl_vol = M_4PI_3 * cube(radius); 
    3838    //total volume 
    3939    double tot_vol = num_pearls * pearl_vol; 
    … …  
    4343    double separation = edge_sep + 2.0 * radius; 
    4444 
    45     double x=q*radius; 
    46  
    47     // Try Taylor on x*xos(x) 
    48         // double out_cos = x - pow(x,3)/2 + pow(x,5)/24 - pow(x,7)/720 + pow(x,9)/40320; 
    49     // psi -= x*out_cos; 
    50  
    5145    //sine functions of a pearl 
    52     double psi = sin(q * radius); 
    53     psi -= x * cos(x); 
    54     psi /= pow((q * radius), 3.0); 
     46    double psi = sph_j1c(q * radius); 
    5547 
    5648    // N pearls contribution 
    … …  
    6153    } 
    6254    // form factor for num_pearls 
    63     double form_factor = n_contrib; 
    64     form_factor *= pow((m_s*psi*3.0), 2.0); 
    65     form_factor /= (tot_vol * 1.0e4); 
     55    double form_factor = 1.0e-4 * n_contrib * square(m_s*psi) / tot_vol; 
    6656 
    6757    return form_factor; 

  • sasmodels/models/linear_pearls.py

    r40a87fa r4962519  
    6363single = False 
    6464 
    65 source = ["linear_pearls.c"] 
     65source = ["lib/sph_j1c.c", "linear_pearls.c"] 
    6666 
    6767demo = dict(scale=1.0, background=0.0, 

  • sasmodels/models/mass_fractal.c

    ra807206 r3a48772  
    3434double form_volume(double radius){ 
    3535 
    36     return 1.333333333333333*M_PI*radius*radius*radius; 
     36    return M_4PI_3*radius*radius*radius; 
    3737} 
    3838 

  • sasmodels/models/mass_surface_fractal.c

    r30fbe2e r3a48772  
    3737double form_volume(double radius){ 
    3838 
    39     return 1.333333333333333*M_PI*radius*radius*radius; 
     39    return M_4PI_3*radius*radius*radius; 
    4040} 
    4141 

  • sasmodels/models/multilayer_vesicle.c

    r2c74c11 r3a48772  
    1919 
    2020        // layer 1 
    21         voli = 4.0*M_PI/3.0*ri*ri*ri; 
     21        voli = M_4PI_3*ri*ri*ri; 
    2222        fval += voli*sldi*sph_j1c(ri*q); 
    2323 
    … …  
    2525 
    2626        // layer 2 
    27         voli = 4.0*M_PI/3.0*ri*ri*ri; 
     27        voli = M_4PI_3*ri*ri*ri; 
    2828        fval -= voli*sldi*sph_j1c(ri*q); 
    2929 

  • sasmodels/models/parallelepiped.c

    ra807206 r14838a3  
    11double form_volume(double length_a, double length_b, double length_c); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a, double length_b, double length_c); 
     2double Iq(double q, double sld, double solvent_sld, 
     3    double length_a, double length_b, double length_c); 
    34double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    4     double length_a, double length_b, double length_c, double theta, double phi, double psi); 
    5  
    6 // From Igor library 
    7 double _pkernel(double a, double b,double c, double ala, double alb, double alc); 
    8 double _pkernel(double a, double b,double c, double ala, double alb, double alc){ 
    9     double argA,argB,argC,tmp1,tmp2,tmp3; 
    10     //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    11     argA = 0.5*a*ala; 
    12     argB = 0.5*b*alb; 
    13     argC = 0.5*c*alc; 
    14     if(argA==0.0) { 
    15         tmp1 = 1.0; 
    16     } else { 
    17         tmp1 = sin(argA)*sin(argA)/argA/argA; 
    18     } 
    19     if (argB==0.0) { 
    20         tmp2 = 1.0; 
    21     } else { 
    22         tmp2 = sin(argB)*sin(argB)/argB/argB; 
    23     } 
    24     if (argC==0.0) { 
    25         tmp3 = 1.0; 
    26     } else { 
    27         tmp3 = sin(argC)*sin(argC)/argC/argC; 
    28     } 
    29     return (tmp1*tmp2*tmp3); 
    30  
    31 } 
    32  
     5    double length_a, double length_b, double length_c, 
     6    double theta, double phi, double psi); 
    337 
    348double form_volume(double length_a, double length_b, double length_c) 
    … …  
    4519    double length_c) 
    4620{ 
    47     double tmp1, tmp2; 
    48      
    49     double mu = q * length_b; 
     21    const double mu = 0.5 * q * length_b; 
    5022     
    5123    // Scale sides by B 
    52     double a_scaled = length_a / length_b; 
    53     double c_scaled = length_c / length_b; 
     24    const double a_scaled = length_a / length_b; 
     25    const double c_scaled = length_c / length_b; 
    5426         
    55     //Order of integration 
    56     int nordi=76;                                
    57     int nordj=76; 
     27    // outer integral (with gauss points), integration limits = 0, 1 
     28    double outer_total = 0; //initialize integral 
    5829 
    59     // outer integral (with nordi gauss points), integration limits = 0, 1 
    60     double summ = 0; //initialize integral 
     30    for( int i=0; i<76; i++) { 
     31        const double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
     32        const double mu_proj = mu * sqrt(1.0-sigma*sigma); 
    6133 
    62     for( int i=0; i<nordi; i++) { 
    63                  
    64         // inner integral (with nordj gauss points), integration limits = 0, 1 
    65          
    66         double summj = 0.0; 
    67             double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );          
    68                  
    69             for(int j=0; j<nordj; j++) { 
     34        // inner integral (with gauss points), integration limits = 0, 1 
     35        // corresponding to angles from 0 to pi/2. 
     36        double inner_total = 0.0; 
     37        for(int j=0; j<76; j++) { 
     38            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
     39            double sin_uu, cos_uu; 
     40            SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
     41            const double si1 = sinc(mu_proj * sin_uu * a_scaled); 
     42            const double si2 = sinc(mu_proj * cos_uu); 
     43            inner_total += Gauss76Wt[j] * square(si1 * si2); 
     44        } 
     45        inner_total *= 0.5; 
    7046 
    71             double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    72             double mudum = mu * sqrt(1.0-sigma*sigma); 
    73                 double arg1 = 0.5 * mudum * cos(0.5*M_PI*uu); 
    74                 double arg2 = 0.5 * mudum * a_scaled * sin(0.5*M_PI*uu); 
    75             if(arg1==0.0) { 
    76                 tmp1 = 1.0; 
    77             } else { 
    78                 tmp1 = sin(arg1)*sin(arg1)/arg1/arg1; 
    79             } 
    80             if (arg2==0.0) { 
    81                 tmp2 = 1.0; 
    82             } else { 
    83                 tmp2 = sin(arg2)*sin(arg2)/arg2/arg2; 
    84             } 
     47        const double si = sinc(mu * c_scaled * sigma); 
     48        outer_total += Gauss76Wt[i] * inner_total * si * si; 
     49    } 
     50    outer_total *= 0.5; 
    8551 
    86             summj += Gauss76Wt[j] * tmp1 * tmp2; 
    87         } 
    88                  
    89         // value of the inner integral 
    90         double answer = 0.5 * summj; 
    91  
    92         double arg = 0.5 * mu * c_scaled * sigma; 
    93         if ( arg == 0.0 ) { 
    94             answer *= 1.0; 
    95         } else { 
    96             answer *= sin(arg)*sin(arg)/arg/arg; 
    97         } 
    98                  
    99             // sum of outer integral 
    100         summ += Gauss76Wt[i] * answer; 
    101          
    102     }    
    103     
    104     const double vd = (sld-solvent_sld) * form_volume(length_a, length_b, length_c); 
    105      
    106     // convert from [1e-12 A-1] to [cm-1] and 0.5 factor for outer integral 
    107     return 1.0e-4 * 0.5 * vd * vd * summ; 
    108      
     52    // Multiply by contrast^2 and convert from [1e-12 A-1] to [cm-1] 
     53    const double V = form_volume(length_a, length_b, length_c); 
     54    const double drho = (sld-solvent_sld); 
     55    return 1.0e-4 * square(drho * V) * outer_total; 
    10956} 
    11057 
    … …  
    12067    double psi) 
    12168{ 
    122     double q = sqrt(qx*qx+qy*qy); 
    123     double qx_scaled = qx/q; 
    124     double qy_scaled = qy/q; 
     69    double q, cos_val_a, cos_val_b, cos_val_c; 
     70    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a); 
    12571 
    126     // Convert angles given in degrees to radians 
    127     theta *= M_PI_180; 
    128     phi   *= M_PI_180; 
    129     psi   *= M_PI_180; 
    130      
    131     // Parallelepiped c axis orientation 
    132     double cparallel_x = cos(theta) * cos(phi); 
    133     double cparallel_y = sin(theta); 
    134      
    135     // Compute angle between q and parallelepiped axis 
    136     double cos_val_c = cparallel_x*qx_scaled + cparallel_y*qy_scaled;// + cparallel_z*qz; 
    137  
    138     // Parallelepiped a axis orientation 
    139     double parallel_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi); 
    140     double parallel_y = sin(psi)*cos(theta); 
    141     double cos_val_a = parallel_x*qx_scaled + parallel_y*qy_scaled; 
    142  
    143     // Parallelepiped b axis orientation 
    144     double bparallel_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi); 
    145     double bparallel_y = cos(theta)*cos(psi); 
    146     double cos_val_b = bparallel_x*qx_scaled + bparallel_y*qy_scaled; 
    147  
    148     // The following tests should always pass 
    149     if (fabs(cos_val_c)>1.0) { 
    150       //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    151       cos_val_c = 1.0; 
    152     } 
    153     if (fabs(cos_val_a)>1.0) { 
    154       //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    155       cos_val_a = 1.0; 
    156     } 
    157     if (fabs(cos_val_b)>1.0) { 
    158       //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    159       cos_val_b = 1.0; 
    160     } 
    161      
    162     // Call the IGOR library function to get the kernel 
    163     double form = _pkernel( q*length_a, q*length_b, q*length_c, cos_val_a, cos_val_b, cos_val_c); 
    164    
    165     // Multiply by contrast^2 
    166     const double vd = (sld - solvent_sld) * form_volume(length_a, length_b, length_c); 
    167     return 1.0e-4 * vd * vd * form; 
     72    const double siA = sinc(0.5*q*length_a*cos_val_a); 
     73    const double siB = sinc(0.5*q*length_b*cos_val_b); 
     74    const double siC = sinc(0.5*q*length_c*cos_val_c); 
     75    const double V = form_volume(length_a, length_b, length_c); 
     76    const double drho = (sld - solvent_sld); 
     77    const double form = V * drho * siA * siB * siC; 
     78    // Square and convert from [1e-12 A-1] to [cm-1] 
     79    return 1.0e-4 * form * form; 
    16880} 

  • sasmodels/models/parallelepiped.py

    r416f5c7 red0827a  
    221221# parameters for demo 
    222222demo = dict(scale=1, background=0, 
    223             sld=6.3e-6, sld_solvent=1.0e-6, 
     223            sld=6.3, sld_solvent=1.0, 
    224224            length_a=35, length_b=75, length_c=400, 
    225225            theta=45, phi=30, psi=15, 

  • sasmodels/models/pearl_necklace.c

    ra807206 r3a48772  
    2323        double num_strings = num_pearls - 1.0; 
    2424         
    25         //Pi 
    26         double pi = 4.0*atan(1.0); 
    27          
    2825        // center to center distance between the neighboring pearls 
    2926        double A_s = edge_sep + 2.0 * radius; 
    … …  
    3633         
    3734        // each volume 
    38         double string_vol = edge_sep * pi * thick_string * thick_string / 4.0; 
    39         double pearl_vol = 4.0 / 3.0 * pi * radius * radius * radius; 
     35        double string_vol = edge_sep * M_PI * thick_string * thick_string / 4.0; 
     36        double pearl_vol = M_4PI_3 * radius * radius * radius; 
    4037 
    4138        //total volume 
    … …  
    110107        double total_vol; 
    111108 
    112         double pi = 4.0*atan(1.0); 
    113109        double number_of_strings = num_pearls - 1.0; 
    114110         
    115         double string_vol = edge_sep * pi * thick_string * thick_string / 4.0; 
    116         double pearl_vol = 4.0 / 3.0 * pi * radius * radius * radius; 
     111        double string_vol = edge_sep * M_PI * thick_string * thick_string / 4.0; 
     112        double pearl_vol = M_4PI_3 * radius * radius * radius; 
    117113 
    118114        total_vol = number_of_strings * string_vol; 

  • sasmodels/models/pearl_necklace.py

    ra807206 r4962519  
    112112    """ 
    113113    tot_vol = volume(radius, edge_sep, thick_string, num_pearls) 
    114     rad_out = pow((3.0*tot_vol/4.0/pi), 0.33333) 
     114    rad_out = (tot_vol/(4.0/3.0*pi)) ** (1./3.) 
    115115    return rad_out 
    116116 

  • sasmodels/models/porod.py

    r40a87fa r4962519  
    4141    """ 
    4242    with errstate(divide='ignore'): 
    43         return power(q, -4) 
     43        return q**-4 
    4444 
    4545Iq.vectorized = True  # Iq accepts an array of q values 

  • sasmodels/models/rectangular_prism.c

    ra807206 rab2aea8  
    22double Iq(double q, double sld, double solvent_sld, double length_a,  
    33          double b2a_ratio, double c2a_ratio); 
    4 double Iqxy(double qx, double qy, double sld, double solvent_sld,  
    5             double length_a, double b2a_ratio, double c2a_ratio); 
    64 
    75double form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
    … …  
    1715    double c2a_ratio) 
    1816{ 
    19     double termA, termB, termC; 
    20      
    21     double b_side = length_a * b2a_ratio; 
    22     double c_side = length_a * c2a_ratio; 
    23     double volume = length_a * b_side * c_side; 
    24     double a_half = 0.5 * length_a; 
    25     double b_half = 0.5 * b_side; 
    26     double c_half = 0.5 * c_side; 
     17    const double length_b = length_a * b2a_ratio; 
     18    const double length_c = length_a * c2a_ratio; 
     19    const double a_half = 0.5 * length_a; 
     20    const double b_half = 0.5 * length_b; 
     21    const double c_half = 0.5 * length_c; 
    2722 
    2823   //Integration limits to use in Gaussian quadrature 
    29     double v1a = 0.0; 
    30     double v1b = 0.5 * M_PI;  //theta integration limits 
    31     double v2a = 0.0; 
    32     double v2b = 0.5 * M_PI;  //phi integration limits 
     24    const double v1a = 0.0; 
     25    const double v1b = M_PI_2;  //theta integration limits 
     26    const double v2a = 0.0; 
     27    const double v2b = M_PI_2;  //phi integration limits 
    3328     
    34     //Order of integration 
    35     int nordi=76;                                
    36     int nordj=76; 
     29    double outer_sum = 0.0; 
     30    for(int i=0; i<76; i++) { 
     31        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
     32        double sin_theta, cos_theta; 
     33        SINCOS(theta, sin_theta, cos_theta); 
    3734 
    38     double sumi = 0.0; 
    39      
    40     for(int i=0; i<nordi; i++) { 
     35        const double termC = sinc(q * c_half * cos_theta); 
    4136 
    42             double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );  
     37        double inner_sum = 0.0; 
     38        for(int j=0; j<76; j++) { 
     39            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
     40            double sin_phi, cos_phi; 
     41            SINCOS(phi, sin_phi, cos_phi); 
    4342 
    44             double arg = q * c_half * cos(theta); 
    45             if (fabs(arg) > 1.e-16) {termC = sin(arg)/arg;} else {termC = 1.0;}   
    46  
    47             double sumj = 0.0; 
    48          
    49             for(int j=0; j<nordj; j++) { 
    50  
    51             double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b );  
    52  
    53                 // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0 
    54  
    55                 arg = q * a_half * sin(theta) * sin(phi);  
    56                 if (fabs(arg) > 1.e-16) {termA = sin(arg)/arg;} else {termA = 1.0;} 
    57                 
    58                 arg = q * b_half * sin(theta) * cos(phi);  
    59                 if (fabs(arg) > 1.e-16) {termB = sin(arg)/arg;} else {termB = 1.0;}        
    60                 
    61                 double AP = termA * termB * termC;   
    62  
    63                 sumj += Gauss76Wt[j] * (AP*AP); 
    64  
    65             } 
    66  
    67             sumj = 0.5 * (v2b-v2a) * sumj; 
    68             sumi += Gauss76Wt[i] * sumj * sin(theta); 
    69  
     43            // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0 
     44            const double termA = sinc(q * a_half * sin_theta * sin_phi); 
     45            const double termB = sinc(q * b_half * sin_theta * cos_phi); 
     46            const double AP = termA * termB * termC; 
     47            inner_sum += Gauss76Wt[j] * AP * AP; 
     48        } 
     49        inner_sum = 0.5 * (v2b-v2a) * inner_sum; 
     50        outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 
    7051    } 
    7152 
    72     double answer = 0.5*(v1b-v1a)*sumi; 
     53    double answer = 0.5*(v1b-v1a)*outer_sum; 
    7354 
    7455    // Normalize by Pi (Eqn. 16).  
    … …  
    7758    // The factor 2 appears because the theta integral has been defined between  
    7859    // 0 and pi/2, instead of 0 to pi. 
    79     answer *= (2.0/M_PI); //Form factor P(q) 
     60    answer /= M_PI_2; //Form factor P(q) 
    8061 
    8162    // Multiply by contrast^2 and volume^2 
    82     answer *= (sld-solvent_sld)*(sld-solvent_sld)*volume*volume; 
     63    const double volume = length_a * length_b * length_c; 
     64    answer *= square((sld-solvent_sld)*volume); 
    8365 
    8466    // Convert from [1e-12 A-1] to [cm-1]  
    … …  
    8668 
    8769    return answer; 
    88      
    8970} 
    90  
    91 double Iqxy(double qx, double qy, 
    92     double sld, 
    93     double solvent_sld, 
    94     double length_a, 
    95     double b2a_ratio, 
    96     double c2a_ratio) 
    97 { 
    98     double q = sqrt(qx*qx + qy*qy); 
    99     double intensity = Iq(q, sld, solvent_sld, length_a, b2a_ratio, c2a_ratio);  
    100     return intensity;     
    101 } 

  • sasmodels/models/rectangular_prism.py

    ra807206 rab2aea8  
    33r""" 
    44 
    5 This model provides the form factor, *P(q)*, for a rectangular prism. 
     5This model provides the form factor, $P(q)$, for a rectangular prism. 
    66 
    77Note that this model is almost totally equivalent to the existing 
    88:ref:`parallelepiped` model. 
    99The only difference is that the way the relevant 
    10 parameters are defined here (*a*, *b/a*, *c/a* instead of *a*, *b*, *c*) 
     10parameters are defined here ($a$, $b/a$, $c/a$ instead of $a$, $b$, $c$) 
    1111which allows use of polydispersity with this model while keeping the shape of 
    12 the prism (e.g. setting *b/a* = 1 and *c/a* = 1 and applying polydispersity 
     12the prism (e.g. setting $b/a = 1$ and $c/a = 1$ and applying polydispersity 
    1313to *a* will generate a distribution of cubes of different sizes). 
    1414Note also that, contrary to :ref:`parallelepiped`, it does not compute 
    … …  
    2424Note also that the angle definitions used in the code and the present 
    2525documentation correspond to those used in (Nayuk, 2012) (see Fig. 1 of 
    26 that reference), with |theta| corresponding to |alpha| in that paper, 
     26that reference), with $\theta$ corresponding to $\alpha$ in that paper, 
    2727and not to the usual convention used for example in the 
    2828:ref:`parallelepiped` model. As the present model does not compute 
    … …  
    3030 
    3131In this model the scattering from a massive parallelepiped with an 
    32 orientation with respect to the scattering vector given by |theta| 
    33 and |phi| 
     32orientation with respect to the scattering vector given by $\theta$ 
     33and $\phi$ 
    3434 
    3535.. math:: 
    36   A_P\,(q) =  \frac{\sin \bigl( q \frac{C}{2} \cos\theta \bigr)}{\left( q \frac{C}{2} 
    37   \cos\theta \right)} \, \times \, \frac{\sin \bigl( q \frac{A}{2} \sin\theta \sin\phi 
    38   \bigr)}{\left( q \frac{A}{2} \sin\theta \sin\phi \right)} \, \times \, \frac{\sin \bigl( 
    39   q \frac{B}{2} \sin\theta \cos\phi \bigr)}{\left( q \frac{B}{2} \sin\theta \cos\phi \right)} 
    4036 
    41 where *A*, *B* and *C* are the sides of the parallelepiped and must fulfill 
    42 :math:`A \le B \le C`, |theta| is the angle between the *z* axis and the 
    43 longest axis of the parallelepiped *C*, and |phi| is the angle between the 
    44 scattering vector (lying in the *xy* plane) and the *y* axis. 
     37  A_P\,(q) = 
     38      \frac{\sin \left( \tfrac{1}{2}qC \cos\theta \right) }{\tfrac{1}{2} qC \cos\theta} 
     39      \,\times\, 
     40      \frac{\sin \left( \tfrac{1}{2}qA \cos\theta \right) }{\tfrac{1}{2} qA \cos\theta} 
     41      \,\times\ , 
     42      \frac{\sin \left( \tfrac{1}{2}qB \cos\theta \right) }{\tfrac{1}{2} qB \cos\theta} 
     43 
     44where $A$, $B$ and $C$ are the sides of the parallelepiped and must fulfill 
     45$A \le B \le C$, $\theta$ is the angle between the $z$ axis and the 
     46longest axis of the parallelepiped $C$, and $\phi$ is the angle between the 
     47scattering vector (lying in the $xy$ plane) and the $y$ axis. 
    4548 
    4649The normalized form factor in 1D is obtained averaging over all possible 
    … …  
    4851 
    4952.. math:: 
    50   P(q) =  \frac{2}{\pi} \times \, \int_0^{\frac{\pi}{2}} \, 
     53  P(q) =  \frac{2}{\pi} \int_0^{\frac{\pi}{2}} \, 
    5154  \int_0^{\frac{\pi}{2}} A_P^2(q) \, \sin\theta \, d\theta \, d\phi 
    5255 
    … …  
    5457 
    5558.. math:: 
    56   I(q) = \mbox{scale} \times V \times (\rho_{\mbox{p}} - 
    57   \rho_{\mbox{solvent}})^2 \times P(q) 
     59  I(q) = \text{scale} \times V \times (\rho_\text{p} - 
     60  \rho_\text{solvent})^2 \times P(q) 
    5861 
    59 where *V* is the volume of the rectangular prism, :math:`\rho_{\mbox{p}}` 
    60 is the scattering length of the parallelepiped, :math:`\rho_{\mbox{solvent}}` 
     62where $V$ is the volume of the rectangular prism, $\rho_\text{p}$ 
     63is the scattering length of the parallelepiped, $\rho_\text{solvent}$ 
    6164is the scattering length of the solvent, and (if the data are in absolute 
    6265units) *scale* represents the volume fraction (which is unitless). 
    … …  
    125128# parameters for demo 
    126129demo = dict(scale=1, background=0, 
    127             sld=6.3e-6, sld_solvent=1.0e-6, 
     130            sld=6.3, sld_solvent=1.0, 
    128131            length_a=35, b2a_ratio=1, c2a_ratio=1, 
    129132            length_a_pd=0.1, length_a_pd_n=10, 

  • sasmodels/models/sc_paracrystal.c

    r0bef47b r4962519  
    4949        double da = d_factor*dnn; 
    5050        double temp1 = qq*qq*da*da; 
    51         double temp2 = pow( 1.0-exp(-1.0*temp1) ,3); 
     51        double temp2 = cube(-expm1(-temp1)); 
    5252        double temp3 = qq*dnn; 
    5353        double temp4 = 2.0*exp(-0.5*temp1); 
    5454        double temp5 = exp(-1.0*temp1); 
    5555 
    56         double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5); 
    57         integrand *= 2.0/M_PI; 
     56        double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5)/M_PI_2; 
    5857 
    5958        return(integrand); 
    6059} 
    6160 
    62 static 
    63 double sc_crystal_kernel(double q, 
     61double Iq(double q, 
    6462          double dnn, 
    6563          double d_factor, 
    … …  
    6967{ 
    7068        const double va = 0.0; 
    71         const double vb = M_PI/2.0; //orientation average, outer integral 
     69        const double vb = M_PI_2; //orientation average, outer integral 
    7270 
    7371    double summ=0.0; 
    … …  
    103101} 
    104102 
    105 static 
    106 double sc_crystal_kernel_2d(double q, double q_x, double q_y, 
     103double Iqxy(double qx, double qy, 
    107104          double dnn, 
    108105          double d_factor, 
    … …  
    114111          double psi) 
    115112{ 
    116     //convert angle degree to radian 
    117     theta = theta * M_PI_180; 
    118     phi = phi * M_PI_180; 
    119     psi = psi * M_PI_180; 
     113    double q, cos_a1, cos_a2, cos_a3; 
     114    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
    120115 
    121     const double qda_2 = pow(q*d_factor*dnn,2.0); 
     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*cos_a1)/cosh_qd) 
     121                    * tanh_qd/(1. - cos(qd*cos_a2)/cosh_qd) 
     122                    * tanh_qd/(1. - cos(qd*cos_a3)/cosh_qd); 
    122123 
    123     double snt, cnt; 
    124     SINCOS(theta, snt, cnt); 
    125  
    126     double snp, cnp; 
    127     SINCOS(phi, snp, cnp); 
    128  
    129     double sns, cns; 
    130     SINCOS(psi, sns, cns); 
    131  
    132     /// Angles here are respect to detector coordinate instead of against 
    133     //  q coordinate(PRB 36, 3, 1754) 
    134     // a3 axis orientation 
    135  
    136     const double a3_x = cnt * cnp; 
    137     const double a3_y = snt; 
    138  
    139     // Compute the angle btw vector q and the a3 axis 
    140     double cos_val_a3 = a3_x*q_x + a3_y*q_y; 
    141  
    142     // a1 axis orientation 
    143     const double a1_x = -cnp*sns * snt+snp*cns; 
    144     const double a1_y = sns*cnt; 
    145  
    146     double cos_val_a1 = a1_x*q_x + a1_y*q_y; 
    147  
    148     // a2 axis orientation 
    149     const double a2_x = -snt*cns*cnp-sns*snp; 
    150     const double a2_y = cnt*cns; 
    151  
    152     // a2 axis 
    153     const double cos_val_a2 =  a2_x*q_x + a2_y*q_y; 
    154  
    155     // The following test should always pass 
    156     if (fabs(cos_val_a3)>1.0) { 
    157         //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    158         cos_val_a3 = 1.0; 
    159     } 
    160     if (fabs(cos_val_a1)>1.0) { 
    161         //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    162         cos_val_a1 = 1.0; 
    163     } 
    164     if (fabs(cos_val_a2)>1.0) { 
    165         //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    166         cos_val_a3 = 1.0; 
    167     } 
    168  
    169     const double a3_dot_q = dnn*q*cos_val_a3; 
    170     const double a1_dot_q = dnn*q*cos_val_a1; 
    171     const double a2_dot_q = dnn*q*cos_val_a2; 
    172  
    173     // Call Zq=Z1*Z2*Z3 
    174     double Zq = (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a1_dot_q)+exp(-qda_2)); 
    175     Zq *= (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a2_dot_q)+exp(-qda_2)); 
    176     Zq *= (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a3_dot_q)+exp(-qda_2)); 
    177  
    178     // Use SphereForm directly from libigor 
    179     double answer = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; 
    180  
    181     //consider scales 
    182     const double latticeScale = sphere_volume(radius/dnn); 
    183     answer *= latticeScale; 
    184  
    185     return answer; 
     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; 
    186128} 
    187  
    188 double Iq(double q, 
    189           double dnn, 
    190           double d_factor, 
    191           double radius, 
    192           double sphere_sld, 
    193           double solvent_sld) 
    194 { 
    195     return sc_crystal_kernel(q, 
    196               dnn, 
    197               d_factor, 
    198               radius, 
    199               sphere_sld, 
    200               solvent_sld); 
    201 } 
    202  
    203 // Iqxy is never called since no orientation or magnetic parameters. 
    204 double Iqxy(double qx, double qy, 
    205             double dnn, 
    206             double d_factor, 
    207             double radius, 
    208             double sphere_sld, 
    209             double solvent_sld, 
    210             double theta, 
    211             double phi, 
    212             double psi) 
    213 { 
    214     double q = sqrt(qx*qx + qy*qy); 
    215  
    216  
    217     return sc_crystal_kernel_2d(q, qx/q, qy/q, 
    218                   dnn, 
    219                   d_factor, 
    220                   radius, 
    221                   sphere_sld, 
    222                   solvent_sld, 
    223                   theta, 
    224                   phi, 
    225                   psi); 
    226  
    227 } 
    228  

  • sasmodels/models/stacked_disks.c

    r0d6e865 r6831fa0  
    1414          double solvent_sld); 
    1515 
     16double Iqxy(double qx, double qy, 
     17          double thick_core, 
     18          double thick_layer, 
     19          double radius, 
     20          double n_stacking, 
     21          double sigma_dnn, 
     22          double core_sld, 
     23          double layer_sld, 
     24          double solvent_sld, 
     25          double theta, 
     26          double phi); 
     27 
    1628static 
    1729double _kernel(double qq, 
    … …  
    2234               double halfheight, 
    2335               double thick_layer, 
    24                double zi, 
     36               double sin_alpha, 
     37               double cos_alpha, 
    2538               double sigma_dnn, 
    2639               double d, 
    … …  
    3447        // zi is the dummy variable for the integration (x in Feigin's notation) 
    3548 
    36         const double besarg1 = qq*radius*sin(zi); 
    37         const double besarg2 = qq*radius*sin(zi); 
     49        const double besarg1 = qq*radius*sin_alpha; 
     50        //const double besarg2 = qq*radius*sin_alpha; 
    3851 
    39         const double sinarg1 = qq*halfheight*cos(zi); 
    40         const double sinarg2 = qq*(halfheight+thick_layer)*cos(zi); 
     52        const double sinarg1 = qq*halfheight*cos_alpha; 
     53        const double sinarg2 = qq*(halfheight+thick_layer)*cos_alpha; 
    4154 
    4255        const double be1 = sas_J1c(besarg1); 
    43         const double be2 = sas_J1c(besarg2); 
    44         const double si1 = sin(sinarg1)/sinarg1; 
    45         const double si2 = sin(sinarg2)/sinarg2; 
     56        //const double be2 = sas_J1c(besarg2); 
     57        const double be2 = be1; 
     58        const double si1 = sinc(sinarg1); 
     59        const double si2 = sinc(sinarg2); 
    4660 
    4761        const double dr1 = (core_sld-solvent_sld); 
    4862        const double dr2 = (layer_sld-solvent_sld); 
    4963        const double area = M_PI*radius*radius; 
    50         const double totald=2.0*(thick_layer+halfheight); 
     64        const double totald = 2.0*(thick_layer+halfheight); 
    5165 
    5266        const double t1 = area*(2.0*halfheight)*dr1*(si1)*(be1); 
    … …  
    5468 
    5569 
    56         double retval =((t1+t2)*(t1+t2))*sin(zi); 
     70        double retval =((t1+t2)*(t1+t2)); 
    5771 
    5872        // loop for the structure facture S(q) 
    5973        double sqq=0.0; 
    6074        for(int kk=1;kk<n_stacking;kk+=1) { 
    61                 double dexpt=qq*cos(zi)*qq*cos(zi)*d*d*sigma_dnn*sigma_dnn*kk/2.0; 
    62                 sqq=sqq+(n_stacking-kk)*cos(qq*cos(zi)*d*kk)*exp(-1.*dexpt); 
     75                double qd_cos_alpha = qq*d*cos_alpha; 
     76                double dexpt=square(qd_cos_alpha*sigma_dnn)*kk/2.0; 
     77                sqq += (n_stacking-kk)*cos(qd_cos_alpha*kk)*exp(-1.*dexpt); 
    6378        } 
    64  
    6579        // end of loop for S(q) 
    6680        sqq=1.0+2.0*sqq/n_stacking; 
    6781 
    68         retval *= sqq; 
    69  
    70         return(retval); 
     82        return retval * sqq; 
    7183} 
    7284 
    … …  
    88100        double summ = 0.0;      //initialize integral 
    89101 
    90         double d=2.0*thick_layer+thick_core; 
    91         double halfheight = thick_core/2.0; 
     102        double d = 2.0*thick_layer+thick_core; 
     103        double halfheight = 0.5*thick_core; 
    92104 
    93         for(int i=0;i<N_POINTS_76;i++) { 
    94                 double zi = (Gauss76Z[i] + 1.0)*M_PI/4.0; 
    95                 double yyy = Gauss76Wt[i] * 
    96                     _kernel(q, 
     105        for(int i=0; i<N_POINTS_76; i++) { 
     106                double zi = (Gauss76Z[i] + 1.0)*M_PI_4; 
     107                double sin_alpha, cos_alpha; // slots to hold sincos function output 
     108                SINCOS(zi, sin_alpha, cos_alpha); 
     109                double yyy = _kernel(q, 
    97110                                   radius, 
    98111                                   core_sld, 
    … …  
    101114                                   halfheight, 
    102115                                   thick_layer, 
    103                                    zi, 
     116                                   sin_alpha, 
     117                                   cos_alpha, 
    104118                                   sigma_dnn, 
    105119                                   d, 
    106120                                   n_stacking); 
    107                 summ += yyy; 
     121                summ += Gauss76Wt[i] * yyy * sin_alpha; 
    108122        } 
    109123 
    110         double answer = M_PI/4.0*summ; 
     124        double answer = M_PI_4*summ; 
    111125 
    112126        //Convert to [cm-1] 
    … …  
    116130} 
    117131 
    118 static double stacked_disks_kernel_2d(double q, double q_x, double q_y, 
    119                             double thick_core, 
    120                             double thick_layer, 
    121                             double radius, 
    122                             double n_stacking, 
    123                             double sigma_dnn, 
    124                             double core_sld, 
    125                             double layer_sld, 
    126                             double solvent_sld, 
    127                             double theta, 
    128                             double phi) 
    129 { 
    130  
    131     double ct, st, cp, sp; 
    132  
    133     //convert angle degree to radian 
    134     theta = theta * M_PI/180.0; 
    135     phi = phi * M_PI/180.0; 
    136  
    137     SINCOS(theta, st, ct); 
    138     SINCOS(phi, sp, cp); 
    139  
    140     // silence compiler warnings about unused variable 
    141     (void) sp; 
    142  
    143     // parallelepiped orientation 
    144     const double cyl_x = st * cp; 
    145     const double cyl_y = st * sp; 
    146  
    147     // Compute the angle btw vector q and the 
    148     // axis of the parallelepiped 
    149     const double cos_val = cyl_x*q_x + cyl_y*q_y; 
    150  
    151     // Note: cos(alpha) = 0 and 1 will get an 
    152     // undefined value from Stackdisc_kern 
    153     double alpha = acos( cos_val ); 
    154  
    155     // Call the IGOR library function to get the kernel 
    156     double d = 2 * thick_layer + thick_core; 
    157     double halfheight = thick_core/2.0; 
    158     double answer = _kernel(q, 
    159                      radius, 
    160                      core_sld, 
    161                      layer_sld, 
    162                      solvent_sld, 
    163                      halfheight, 
    164                      thick_layer, 
    165                      alpha, 
    166                      sigma_dnn, 
    167                      d, 
    168                      n_stacking); 
    169  
    170     answer /= sin(alpha); 
    171     //convert to [cm-1] 
    172     answer *= 1.0e-4; 
    173  
    174     return answer; 
    175 } 
    176  
    177132double form_volume(double thick_core, 
    178133                   double thick_layer, 
    179134                   double radius, 
    180135                   double n_stacking){ 
    181     double d = 2 * thick_layer + thick_core; 
    182     return acos(-1.0) * radius * radius * d * n_stacking; 
     136    double d = 2.0 * thick_layer + thick_core; 
     137    return M_PI * radius * radius * d * n_stacking; 
    183138} 
    184139 
    … …  
    203158                    solvent_sld); 
    204159} 
     160 
     161 
     162double 
     163Iqxy(double qx, double qy, 
     164     double thick_core, 
     165     double thick_layer, 
     166     double radius, 
     167     double n_stacking, 
     168     double sigma_dnn, 
     169     double core_sld, 
     170     double layer_sld, 
     171     double solvent_sld, 
     172     double theta, 
     173     double phi) 
     174{ 
     175    double q, sin_alpha, cos_alpha; 
     176    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     177 
     178    double d = 2.0 * thick_layer + thick_core; 
     179    double halfheight = 0.5*thick_core; 
     180    double answer = _kernel(q, 
     181                     radius, 
     182                     core_sld, 
     183                     layer_sld, 
     184                     solvent_sld, 
     185                     halfheight, 
     186                     thick_layer, 
     187                     sin_alpha, 
     188                     cos_alpha, 
     189                     sigma_dnn, 
     190                     d, 
     191                     n_stacking); 
     192 
     193    //convert to [cm-1] 
     194    answer *= 1.0e-4; 
     195 
     196    return answer; 
     197} 
     198 

  • sasmodels/models/star_polymer.c

    r2c74c11 r3a48772  
    66{ 
    77 
    8     double u_2 = radius2 * pow(q,2); 
     8    double u_2 = radius2 * q * q; 
    99    double v = u_2 * arms / (3.0 * arms - 2.0); 
    1010 
    11     double term1 = v - 1.0 + exp(-v); 
    12     double term2 = ((arms - 1.0)/2.0)* pow((1.0 - exp(-v)),2.0); 
     11    double term1 = v + expm1(-v); 
     12    double term2 = ((arms - 1.0)/2.0) * square(expm1(-v)); 
    1313 
    14     return (2.0 * (term1 + term2)) / (arms * pow(v,2.0)); 
     14    return (2.0 * (term1 + term2)) / (arms * v * v); 
    1515 
    1616} 

  • sasmodels/models/surface_fractal.c

    ra807206 rb716cc6  
     1// Don't need invalid test since fractal_dim_surf is not polydisperse 
     2// #define INVALID(v) (v.fractal_dim_surf <= 1.0 || v.fractal_dim_surf >= 3.0) 
     3 
    14double form_volume(double radius); 
    25 
    … …  
    1114    double cutoff_length) 
    1215{ 
    13     double pq, sq, mmo, result; 
     16    // calculate P(q) 
     17    const double pq = square(sph_j1c(q*radius)); 
    1418 
    15     //Replaced the original formula with Taylor expansion near zero. 
    16     //pq = pow((3.0*(sin(q*radius) - q*radius*cos(q*radius))/pow((q*radius),3)),2); 
     19    // calculate S(q) 
     20    // Note: lim q->0 S(q) = -gamma(mmo) cutoff_length^mmo (mmo cutoff_length) 
     21    // however, the surface fractal formula is invalid outside the range 
     22    const double mmo = 5.0 - fractal_dim_surf; 
     23    const double sq = sas_gamma(mmo) * pow(cutoff_length, mmo) 
     24           * pow(1.0 + square(q*cutoff_length), -0.5*mmo) 
     25           * sin(-mmo * atan(q*cutoff_length)) / q; 
    1726 
    18     pq = sph_j1c(q*radius); 
    19     pq = pq*pq; 
     27    // Empirically determined that the results are valid within this range. 
     28    // Above 1/r, the form starts to oscillate;  below 
     29    //const double result = (q > 5./(3-fractal_dim_surf)/cutoff_length) && q < 1./radius 
     30    //                      ? pq * sq : 0.); 
    2031 
    21     //calculate S(q) 
    22     mmo = 5.0 - fractal_dim_surf; 
    23     sq  = sas_gamma(mmo)*sin(-(mmo)*atan(q*cutoff_length)); 
    24     sq *= pow(cutoff_length, mmo); 
    25     sq /= pow((1.0 + (q*cutoff_length)*(q*cutoff_length)),(mmo/2.0)); 
    26     sq /= q; 
     32    double result = pq * sq; 
    2733 
    28     //combine and return 
    29     result = pq * sq; 
    30  
    31     return result; 
     34    // exclude negative results 
     35    return result > 0. ? result : 0.; 
    3236} 
    33 double form_volume(double radius){ 
    34  
    35     return 1.333333333333333*M_PI*radius*radius*radius; 
     37double form_volume(double radius) 
     38{ 
     39    return M_4PI_3*cube(radius); 
    3640} 
    3741 

  • sasmodels/models/surface_fractal.py

    ra807206 rb716cc6  
    1010.. math:: 
    1111 
    12     I(q) = scale \times P(q)S(q) + background 
    13  
    14 .. math:: 
    15  
    16     P(q) = F(qR)^2 
    17  
    18 .. math:: 
    19  
    20     F(x) = \frac{3\left[sin(x)-xcos(x)\right]}{x^3} 
    21  
    22 .. math:: 
    23  
    24     S(q) = \frac{\Gamma(5-D_S)\zeta^{5-D_S}}{\left[1+(q\zeta)^2 
    25     \right]^{(5-D_S)/2}} 
    26     \frac{sin\left[(D_S - 5) tan^{-1}(q\zeta) \right]}{q} 
    27  
    28 .. math:: 
    29  
    30     scale = scale\_factor \times NV^2(\rho_{particle} - \rho_{solvent})^2 
    31  
    32 .. math:: 
    33  
    34     V = \frac{4}{3}\pi R^3 
     12    I(q) &= \text{scale} \times P(q)S(q) + \text{background} \\ 
     13    P(q) &= F(qR)^2 \\ 
     14    F(x) &= \frac{3\left[\sin(x)-x\cos(x)\right]}{x^3} \\ 
     15    S(q) &= \Gamma(5-D_S)\xi^{\,5-D_S}\left[1+(q\xi)^2 \right]^{-(5-D_S)/2} 
     16            \sin\left[-(5-D_S) \tan^{-1}(q\xi) \right] q^{-1} \\ 
     17    \text{scale} &= \phi N V^2(\rho_\text{particle} - \rho_\text{solvent})^2 \\ 
     18    V &= \frac{4}{3}\pi R^3 
    3519 
    3620where $R$ is the radius of the building block, $D_S$ is the **surface** fractal 
    37 dimension,| \zeta\|  is the cut-off length, $\rho_{solvent}$ is the scattering 
    38 length density of the solvent, 
    39 and $\rho_{particle}$ is the scattering length density of particles. 
     21dimension, $\xi$ is the cut-off length, $\rho_\text{solvent}$ is the scattering 
     22length density of the solvent, $\rho_\text{particle}$ is the scattering 
     23length density of particles and $\phi$ is the volume fraction of the particles. 
    4024 
    4125.. note:: 
    42     The surface fractal dimension $D_s$ is only valid if $1<surface\_dim<3$. 
    43     It is also only valid over a limited $q$ range (see the reference for 
    44     details) 
     26 
     27    The surface fractal dimension is only valid if $1<D_S<3$. The result is 
     28    only valid over a limited $q$ range, $\tfrac{5}{3-D_S}\xi^{\,-1} < q < R^{-1}$. 
     29    See the reference for details. 
    4530 
    4631 
    … …  
    8974source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "surface_fractal.c"] 
    9075 
    91 demo = dict(scale=1, background=0, 
     76demo = dict(scale=1, background=1e-5, 
    9277            radius=10, fractal_dim_surf=2.0, cutoff_length=500) 
    9378 

  • sasmodels/models/teubner_strey.py

    rb3f2a24 r4962519  
    7171 
    7272import numpy as np 
    73 from numpy import inf,power,pi 
     73from numpy import inf, pi 
    7474 
    7575name = "teubner_strey" 
    … …  
    9393    drho2 = (sld-sld_solvent)*(sld-sld_solvent) 
    9494    k = 2.0*pi*xi/d 
    95     a2 = power(1.0+power(k,2.0),2.0) 
    96     c1 = -2.0*xi*xi*power(k,2.0)+2*xi*xi 
    97     c2 = power(xi,4.0) 
     95    a2 = (1.0 + k**2)**2 
     96    c1 = 2.0*xi**2 * (1.0 - k**2) 
     97    c2 = xi**4 
    9898    prefactor = 8.0*pi*volfraction*(1.0-volfraction)*drho2*c2/xi 
    9999    #k2 = (2.0*pi/d)*(2.0*pi/d) 

  • sasmodels/models/triaxial_ellipsoid.c

    ra807206 r3a48772  
    1010double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
    1111{ 
    12     return 1.333333333333333*M_PI*radius_equat_minor*radius_equat_major*radius_polar; 
     12    return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar; 
    1313} 
    1414 
    … …  
    5858    double psi) 
    5959{ 
    60     double stheta, ctheta; 
    61     double sphi, cphi; 
    62     double spsi, cpsi; 
     60    double q, calpha, cmu, cnu; 
     61    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, calpha, cmu, cnu); 
    6362 
    64     const double q = sqrt(qx*qx + qy*qy); 
    65     const double qxhat = qx/q; 
    66     const double qyhat = qy/q; 
    67     SINCOS(theta*M_PI_180, stheta, ctheta); 
    68     SINCOS(phi*M_PI_180, sphi, cphi); 
    69     SINCOS(psi*M_PI_180, spsi, cpsi); 
    70     const double calpha = ctheta*cphi*qxhat + stheta*qyhat; 
    71     const double cnu = (-cphi*spsi*stheta + sphi*cpsi)*qxhat + spsi*ctheta*qyhat; 
    72     const double cmu = (-stheta*cpsi*cphi - spsi*sphi)*qxhat + ctheta*cpsi*qyhat; 
    7363    const double t = q*sqrt(radius_equat_minor*radius_equat_minor*cnu*cnu 
    7464                          + radius_equat_major*radius_equat_major*cmu*cmu 

  • sasmodels/models/vesicle.c

    r2c74c11 r3a48772  
    88{ 
    99    //note that for the vesicle model, the volume is ONLY the shell volume 
    10     double volume; 
    11     volume =4.*M_PI*(radius+thickness)*(radius+thickness)*(radius+thickness)/3; 
    12     volume -=4.*M_PI*radius*radius*radius/3.; 
    13     return volume; 
     10    return M_4PI_3*(cube(radius+thickness) - cube(radius)); 
    1411} 
    1512 
    … …  
    3229    // core first, then add in shell 
    3330    contrast = sld_solvent-sld; 
    34     vol = 4.0*M_PI/3.0*radius*radius*radius; 
    35     f = vol*sph_j1c(q*radius)*contrast; 
     31    vol = M_4PI_3*cube(radius); 
     32    f = vol * sph_j1c(q*radius) * contrast; 
    3633  
    3734    //now the shell. No volume normalization as this is done by the caller 
    3835    contrast = sld-sld_solvent; 
    39     vol = 4.0*M_PI/3.0*(radius+thickness)*(radius+thickness)*(radius+thickness); 
    40     f += vol*sph_j1c(q*(radius+thickness))*contrast; 
     36    vol = M_4PI_3*cube(radius+thickness); 
     37    f += vol * sph_j1c(q*(radius+thickness)) * contrast; 
    4138 
    4239    //rescale to [cm-1].  
    43     f2 = volfraction*f*f*1.0e-4; 
     40    f2 = volfraction * f*f*1.0e-4; 
    4441     
    45     return(f2); 
     42    return f2; 
    4643} 

  • sasmodels/models/vesicle.py

    re77872e r3a48772  
    116116    ''' 
    117117 
    118     whole = 4. * pi * (radius + thickness) ** 3. / 3. 
    119     core = 4. * pi * radius ** 3. / 3. 
     118    whole = 4./3. * pi * (radius + thickness)**3 
     119    core = 4./3. * pi * radius**3 
    120120    return whole, whole - core 
    121121 
Note: See TracChangeset for help on using the changeset viewer.