Changeset f4cf580 in sasmodels


Ignore:
Timestamp:
Sep 2, 2014 1:15:40 PM (11 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
87985ca
Parents:
5d4777d
Message:

resolve remaining differences between sasview and sasmodels

Files:
7 edited

Legend:

Unmodified
Added
Removed

  • compare-new.py

    r5d4777d rf4cf580  
    33 
    44import sys 
     5import math 
    56 
    67import numpy as np 
    … …  
    1516    Load a sasview model given the model name. 
    1617    """ 
    17     modelname = modelname 
     18    # convert model parameters from sasmodel form to sasview form 
     19    #print "old",sorted(pars.items()) 
     20    modelname, pars = revert_model(modelname, pars) 
     21    #print "new",sorted(pars.items()) 
    1822    sans = __import__('sans.models.'+modelname) 
    1923    ModelClass = getattr(getattr(sans.models,modelname,None),modelname,None) 
    … …  
    4145    return kernel 
    4246 
    43 def load_dll(modelname, dtype='single'): 
     47def load_ctypes(modelname, dtype='single'): 
    4448    sasmodels = __import__('sasmodels.models.'+modelname) 
    4549    module = getattr(sasmodels.models, modelname, None) 
    … …  
    4953def randomize(p, v): 
    5054    """ 
    51     Randomizing pars using sasview names. 
    52     Guessing at angles and slds. 
    53     """ 
    54     if any(p.endswith(s) for s in ('_pd_n','_pd_nsigmas','_pd_type')): 
     55    Randomizing parameter. 
     56 
     57    Guess the parameter type from name. 
     58    """ 
     59    if any(p.endswith(s) for s in ('_pd_n','_pd_nsigma','_pd_type')): 
    5560        return v 
    56     if any(s in p for s in ('theta','phi','psi')): 
    57         return np.random.randint(0,90) 
    58     elif any(s in p for s in ('sld','Sld','SLD')): 
    59         return np.random.rand()*1e-5 
     61    elif any(s in p for s in ('theta','phi','psi')): 
     62        # orientation in [-180,180], orientation pd in [0,45] 
     63        if p.endswith('_pd'): 
     64            return 45*np.random.rand() 
     65        else: 
     66            return 360*np.random.rand() - 180 
     67    elif 'sld' in p: 
     68        # sld in in [-0.5,10] 
     69        return 10.5*np.random.rand() - 0.5 
     70    elif p.endswith('_pd'): 
     71        # length pd in [0,1] 
     72        return np.random.rand() 
    6073    else: 
    61         return np.random.rand() 
     74        # length, scale, background in [0,200] 
     75        return 200*np.random.rand() 
    6276 
    6377def parlist(pars): 
    … …  
    7488        if p.endswith("_pd"): pars[p] = 0 
    7589 
    76 def compare(gpuname, gpupars, Ncpu, Ngpu, opts): 
    77     #from sasmodels.core import load_data 
    78     #data = load_data('December/DEC07098.DAT') 
     90def compare(name, pars, Ncpu, Ngpu, opts, set_pars): 
     91 
     92    # Sort out data 
    7993    qmax = 1.0 if '-highq' in opts else (0.2 if '-midq' in opts else 0.05) 
    8094    if "-1d" in opts: 
    8195        from sasmodels.bumps_model import empty_data1D 
    82         qmax = np.log10(qmax) 
     96        qmax = math.log10(qmax) 
    8397        data = empty_data1D(np.logspace(qmax-3, qmax, 128)) 
     98        index = slice(None, None) 
    8499    else: 
    85100        from sasmodels.bumps_model import empty_data2D, set_beam_stop 
    86101        data = empty_data2D(np.linspace(-qmax, qmax, 128)) 
    87102        set_beam_stop(data, 0.004) 
     103        index = ~data.mask 
    88104    is2D = hasattr(data, 'qx_data') 
     105 
     106    # modelling accuracy is determined by dtype and cutoff 
    89107    dtype = 'double' if '-double' in opts else 'single' 
    90108    cutoff_opts = [s for s in opts if s.startswith('-cutoff')] 
    91109    cutoff = float(cutoff_opts[0].split('=')[1]) if cutoff_opts else 1e-5 
    92110 
    93  
    94     if '-random' in opts: 
    95         gpupars = dict((p,randomize(p,v)) for p,v in gpupars.items()) 
    96     if '-pars' in opts: print "pars",parlist(gpupars) 
    97     if '-mono' in opts: suppress_pd(gpupars) 
    98  
    99     cpuname, cpupars = revert_model(gpuname, gpupars) 
    100  
    101     try: 
    102         gpumodel = load_opencl(gpuname, dtype=dtype) 
    103     except Exception,exc: 
    104         print exc 
    105         print "... trying again with single precision" 
    106         gpumodel = load_opencl(gpuname, dtype='single') 
    107     model = BumpsModel(data, gpumodel, cutoff=cutoff, **gpupars) 
     111    # randomize parameters 
     112    random_opts = [s for s in opts if s.startswith('-random')] 
     113    if random_opts: 
     114        if '=' in random_opts[0]: 
     115            seed = int(random_opts[0].split('=')[1]) 
     116        else: 
     117            seed = int(np.random.rand()*10000) 
     118        np.random.seed(seed) 
     119        print "Randomize using -random=%i"%seed 
     120        # Note: the sort guarantees order of calls to random number generator 
     121        pars = dict((p,randomize(p,v)) for p,v in sorted(pars.items())) 
     122        # The capped cylinder model has a constraint on its parameters 
     123        if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']: 
     124            pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius'] 
     125    pars.update(set_pars) 
     126 
     127    # parameter selection 
     128    if '-mono' in opts: 
     129        suppress_pd(pars) 
     130    if '-pars' in opts: 
     131        print "pars",parlist(pars) 
     132 
     133    # OpenCl calculation 
    108134    if Ngpu > 0: 
     135        try: 
     136            model = load_opencl(name, dtype=dtype) 
     137        except Exception,exc: 
     138            print exc 
     139            print "... trying again with single precision" 
     140            model = load_opencl(name, dtype='single') 
     141        problem = BumpsModel(data, model, cutoff=cutoff, **pars) 
    109142        toc = tic() 
    110         for i in range(Ngpu): 
     143        for _ in range(Ngpu): 
    111144            #pars['scale'] = np.random.rand() 
    112             model.update() 
    113             gpu = model.theory() 
     145            problem.update() 
     146            gpu = problem.theory() 
    114147        gpu_time = toc()*1000./Ngpu 
    115         print "ocl t=%.1f ms, intensity=%f"%(gpu_time, sum(gpu[~np.isnan(gpu)])) 
     148        print "opencl t=%.1f ms, intensity=%.0f"%(gpu_time, sum(gpu[index])) 
    116149        #print max(gpu), min(gpu) 
    117150 
    118     comp = None 
    119     if Ncpu > 0 and "-dll" in opts: 
    120         dllmodel = load_dll(gpuname, dtype=dtype) 
    121         model = BumpsModel(data, dllmodel, cutoff=cutoff, **gpupars) 
     151    # ctypes/sasview calculation 
     152    if Ncpu > 0 and "-ctypes" in opts: 
     153        model = load_ctypes(name, dtype=dtype) 
     154        problem = BumpsModel(data, model, cutoff=cutoff, **pars) 
    122155        toc = tic() 
    123         for i in range(Ncpu): 
    124             model.update() 
    125             cpu = model.theory() 
     156        for _ in range(Ncpu): 
     157            problem.update() 
     158            cpu = problem.theory() 
    126159        cpu_time = toc()*1000./Ncpu 
    127         comp = "dll" 
    128  
    129     elif 0: # Hack to check new vs old for GpuCylinder 
    130         from Models.code_cylinder_f import GpuCylinder as oldgpu 
    131         from sasmodel import SasModel 
    132         oldmodel = SasModel(data, oldgpu, dtype=dtype, **cpupars) 
     160        comp = "ctypes" 
     161        print "ctypes t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index])) 
     162    elif Ncpu > 0: 
     163        model = sasview_model(name, **pars) 
    133164        toc = tic() 
    134         for i in range(Ngpu): 
    135             oldmodel.update() 
    136             cpu = oldmodel.theory() 
    137         cpu_time = toc()*1000./Ngpu 
    138         comp = "old" 
    139  
    140     elif Ncpu > 0: 
    141         cpumodel = sasview_model(cpuname, **cpupars) 
    142         toc = tic() 
    143         if is2D: 
    144             for i in range(Ncpu): 
    145                 cpu = cpumodel.evalDistribution([data.qx_data, data.qy_data]) 
    146         else: 
    147             for i in range(Ncpu): 
    148                 cpu = cpumodel.evalDistribution(data.x) 
     165        for _ in range(Ncpu): 
     166            if is2D: 
     167                cpu = model.evalDistribution([data.qx_data, data.qy_data]) 
     168            else: 
     169                cpu = model.evalDistribution(data.x) 
    149170        cpu_time = toc()*1000./Ncpu 
    150         comp = 'sasview' 
    151         #print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[model.index])) 
    152  
    153     if comp: 
    154         print "%s t=%.1f ms, intensity=%f"%(comp, cpu_time, sum(cpu[model.index])) 
     171        comp = "sasview" 
     172        print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index])) 
     173 
     174    # Compare, but only if computing both forms 
    155175    if Ngpu > 0 and Ncpu > 0: 
    156176        #print "speedup %.2g"%(cpu_time/gpu_time) 
    … …  
    158178        #cpu *= max(gpu/cpu) 
    159179        resid, relerr = np.zeros_like(gpu), np.zeros_like(gpu) 
    160         resid[model.index] = (gpu - cpu)[model.index] 
    161         relerr[model.index] = resid[model.index]/cpu[model.index] 
    162         print "max(|ocl-%s|)"%comp, max(abs(resid[model.index])) 
    163         print "max(|(ocl-%s)/ocl|)"%comp, max(abs(relerr[model.index])) 
    164  
     180        resid[index] = (gpu - cpu)[index] 
     181        relerr[index] = resid[index]/cpu[index] 
     182        print "max(|ocl-%s|)"%comp, max(abs(resid[index])) 
     183        print "max(|(ocl-%s)/ocl|)"%comp, max(abs(relerr[index])) 
     184 
     185    # Plot if requested 
    165186    if '-noplot' in opts: return 
    166  
    167187    import matplotlib.pyplot as plt 
    168188    if Ncpu > 0: 
    … …  
    173193        if Ncpu > 0: plt.subplot(132) 
    174194        plot_data(data, gpu, scale='log') 
    175         plt.title("ocl t=%.1f ms"%gpu_time) 
     195        plt.title("opencl t=%.1f ms"%gpu_time) 
    176196    if Ncpu > 0 and Ngpu > 0: 
    177197        plt.subplot(133) 
    178198        err = resid if '-abs' in opts else relerr 
     199        errstr = "abs err" if '-abs' in opts else "rel err" 
     200        #err,errstr = gpu/cpu,"ratio" 
    179201        plot_data(data, err, scale='linear') 
    180         plt.title("max rel err = %.3g"%max(abs(err[model.index]))) 
     202        plt.title("max %s = %.3g"%(errstr, max(abs(err[index])))) 
    181203        if is2D: plt.colorbar() 
    182204    plt.show() 
    … …  
    200222    -lowq/-midq/-highq use q values up to 0.05, 0.2 or 1.0 
    201223    -2d/-1d uses 1d or 2d random data 
    202     -preset/-random randomizes the parameters 
     224    -preset/-random[=seed] preset or random parameters 
    203225    -poly/-mono force monodisperse/polydisperse 
    204     -sasview/-dll whether cpu is tested using sasview or dll 
     226    -sasview/-ctypes whether cpu is tested using sasview or ctypes 
    205227    -cutoff=1e-5/value cutoff for including a point in polydispersity 
    206228    -nopars/-pars prints the parameter set or not 
    … …  
    214236    %s 
    215237""" 
     238 
     239VALID_OPTIONS = [ 
     240    'plot','noplot', 
     241    'single','double', 
     242    'lowq','midq','highq', 
     243    '2d','1d', 
     244    'preset','random', 
     245    'poly','mono', 
     246    'sasview','ctypes', 
     247    'nopars','pars', 
     248    'rel','abs', 
     249    ] 
    216250 
    217251def main(): 
    … …  
    227261        sys.exit(1) 
    228262 
     263    valid_opts = set(VALID_OPTIONS) 
     264    invalid = [o[1:] for o in opts 
     265               if o[1:] not in valid_opts 
     266                  and not o.startswith('-cutoff=') 
     267                  and not o.startswith('-random=')] 
     268    if invalid: 
     269        print "Invalid options: %s"%(", ".join(invalid)) 
     270        sys.exit(1) 
     271 
    229272    name, pars = MODELS[args[0]]() 
    230273    Nopencl = int(args[1]) if len(args) > 1 else 5 
    … …  
    238281 
    239282    # Fill in parameters given on the command line 
     283    set_pars = {} 
    240284    for arg in args[3:]: 
    241285        k,v = arg.split('=') 
    242286        if k not in pars: 
    243287            # extract base name without distribution 
    244             # style may be either a.d or a_pd_d 
    245             s = set((p.split('_pd')[0]).split('.')[0] for p in pars) 
     288            s = set(p.split('_pd')[0] for p in pars) 
    246289            print "%r invalid; parameters are: %s"%(k,", ".join(sorted(s))) 
    247290            sys.exit(1) 
    248         pars[k] = float(v) if not v.endswith('type') else v 
    249  
    250     compare(name, pars, Nsasview, Nopencl, opts) 
     291        set_pars[k] = float(v) if not v.endswith('type') else v 
     292 
     293    compare(name, pars, Nsasview, Nopencl, opts, set_pars) 
    251294 
    252295# =========================================================================== 
    … …  
    265308    pars = dict( 
    266309        scale=1, background=0, 
    267         sld=6e-6, solvent_sld=1e-6, 
     310        sld=6, solvent_sld=1, 
    268311        #radius=5, length=20, 
    269312        radius=260, length=290, 
    270         theta=30, phi=15, 
     313        theta=30, phi=0, 
    271314        radius_pd=.2, radius_pd_n=1, 
    272315        length_pd=.2,length_pd_n=1, 
    273         theta_pd=15, theta_pd_n=25, 
     316        theta_pd=15, theta_pd_n=45, 
    274317        phi_pd=15, phi_pd_n=1, 
    275318        ) 
    … …  
    280323    pars = dict( 
    281324        scale=1, background=0, 
    282         sld=6e-6, solvent_sld=1e-6, 
    283         radius=260, cap_radius=80000, length=290, 
     325        sld=6, solvent_sld=1, 
     326        radius=260, cap_radius=290, length=290, 
    284327        theta=30, phi=15, 
    285328        radius_pd=.2, radius_pd_n=1, 
    286329        cap_radius_pd=.2, cap_radius_pd_n=1, 
    287330        length_pd=.2, length_pd_n=1, 
    288         theta_pd=15, theta_pd_n=25, 
     331        theta_pd=15, theta_pd_n=45, 
    289332        phi_pd=15, phi_pd_n=1, 
    290333        ) 
    … …  
    296339    pars = dict( 
    297340        scale=1, background=0, 
    298         core_sld=6e-6, shell_sld=8e-6, solvent_sld=1e-6, 
    299         radius=325, thickness=25, length=34.2709, 
    300         theta=90, phi=0, 
    301         radius_pd=0.1, radius_pd_n=10, 
    302         length_pd=0.1, length_pd_n=5, 
    303         thickness_pd=0.1, thickness_pd_n=5, 
    304         theta_pd=15.8, theta_pd_n=5, 
    305         phi_pd=0.0008748, phi_pd_n=5, 
     341        core_sld=6, shell_sld=8, solvent_sld=1, 
     342        radius=45, thickness=25, length=340, 
     343        theta=30, phi=15, 
     344        radius_pd=.2, radius_pd_n=1, 
     345        length_pd=.2, length_pd_n=1, 
     346        thickness_pd=.2, thickness_pd_n=1, 
     347        theta_pd=15, theta_pd_n=45, 
     348        phi_pd=15, phi_pd_n=1, 
    306349        ) 
    307350    return 'core_shell_cylinder', pars 
    … …  
    312355    pars = dict( 
    313356        scale=1, background=0, 
    314         sld=6e-6, solvent_sld=1e-6, 
     357        sld=6, solvent_sld=1, 
    315358        rpolar=50, requatorial=30, 
    316         theta=0, phi=0, 
    317         rpolar_pd=0.3, rpolar_pd_n=10, 
    318         requatorial_pd=0, requatorial_pd_n=10, 
    319         theta_pd=0, theta_pd_n=45, 
    320         phi_pd=0, phi_pd_n=45, 
     359        theta=30, phi=15, 
     360        rpolar_pd=.2, rpolar_pd_n=1, 
     361        requatorial_pd=.2, requatorial_pd_n=1, 
     362        theta_pd=15, theta_pd_n=45, 
     363        phi_pd=15, phi_pd_n=1, 
    321364        ) 
    322365    return 'ellipsoid', pars 
    … …  
    327370    pars = dict( 
    328371        scale=1, background=0, 
    329         sld=6e-6, solvent_sld=1e-6, 
    330         theta=0, phi=0, psi=0, 
     372        sld=6, solvent_sld=1, 
     373        theta=30, phi=15, psi=5, 
    331374        req_minor=25, req_major=36, rpolar=50, 
    332         theta_pd=0, theta_pd_n=5, 
    333         phi_pd=0, phi_pd_n=5, 
    334         psi_pd=0, psi_pd_n=5, 
    335         req_minor_pd=0, req_minor_pd_n=5, 
    336         req_major_pd=0, req_major_pd_n=5, 
    337         rpolar_pd=.3, rpolar_pd_n=25, 
     375        req_minor_pd=0, req_minor_pd_n=1, 
     376        req_major_pd=0, req_major_pd_n=1, 
     377        rpolar_pd=.2, rpolar_pd_n=1, 
     378        theta_pd=15, theta_pd_n=45, 
     379        phi_pd=15, phi_pd_n=1, 
     380        psi_pd=15, psi_pd_n=1, 
    338381        ) 
    339382    return 'triaxial_ellipsoid', pars 
    … …  
    343386    pars = dict( 
    344387        scale=1, background=0, 
    345         sld=6e-6, solvent_sld=1e-6, 
     388        sld=6, solvent_sld=1, 
    346389        radius=120, 
    347         radius_pd=.3, radius_pd_n=5, 
     390        radius_pd=.2, radius_pd_n=45, 
    348391        ) 
    349392    return 'sphere', pars 
    … …  
    353396    pars = dict( 
    354397        scale=1, background=0, 
    355         sld=6e-6, solvent_sld=1e-6, 
     398        sld=6, solvent_sld=1, 
    356399        thickness=40, 
    357         thickness_pd= 0.3, thickness_pd_n=40, 
     400        thickness_pd= 0.2, thickness_pd_n=40, 
    358401        ) 
    359402    return 'lamellar', pars 

  • sasmodels/convert.py

    r5d4777d rf4cf580  
    3737        name='sphere', 
    3838        sldSph='sld', sldSolv='solvent_sld', 
    39         #radius='radius',  # DO NOT LIST IDENTICAL PARAMETERS! 
     39        radius='radius',  # listing identical parameters is optional 
    4040        ), 
    4141    } 
    … …  
    5959    newpars = pars.copy() 
    6060    for old,new in mapping.items(): 
     61        if old == new: continue 
    6162        # Bumps style parameter names 
    6263        for variant in ("", "_pd", "_pd_n", "_pd_nsigma", "_pd_type"): 

  • sasmodels/gen.py

    r5d4777d rf4cf580  
    352352  const real I = %(fn)s(%(qcall)s, %(pcall)s); 
    353353  if (I>=REAL(0.0)) { // scattering cannot be negative 
    354     ret += weight*I; 
     354    ret += weight*I%(sasview_spherical)s; 
    355355    norm += weight; 
    356356    %(volume_norm)s 
    … …  
    364364SPHERICAL_CORRECTION="""\ 
    365365// Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta 
    366 real spherical_correction = (Ntheta>1 ? fabs(cos(M_PI_180*phi)) : REAL(1.0));\ 
     366real spherical_correction = (Ntheta>1 ? fabs(sin(M_PI_180*theta)) : REAL(1.0));\ 
     367""" 
     368# Use this to reproduce sasview behaviour 
     369SASVIEW_SPHERICAL_CORRECTION="""\ 
     370// Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta 
     371real spherical_correction = (Ntheta>1 ? fabs(cos(M_PI_180*theta))*M_PI_2 : REAL(1.0));\ 
    367372""" 
    368373 
    … …  
    472477    fq_pars = [p[0] for p in info['parameters'][len(COMMON_PARAMETERS):] 
    473478               if p[0] in set(fixed_pars+pd_pars)] 
    474     if False and "phi" in pd_pars: 
     479    if False and "theta" in pd_pars: 
    475480        spherical_correction = [indent(SPHERICAL_CORRECTION, depth)] 
    476481        weights = [p+"_w" for p in pd_pars]+['spherical_correction'] 
     482        sasview_spherical = "" 
     483    elif "theta" in pd_pars: 
     484        spherical_correction = [indent(SASVIEW_SPHERICAL_CORRECTION,depth)] 
     485        weights = [p+"_w" for p in pd_pars] 
     486        sasview_spherical = "*spherical_correction" 
    477487    else: 
    478488        spherical_correction = [] 
    479489        weights = [p+"_w" for p in pd_pars] 
     490        sasview_spherical = "" 
    480491    subst = { 
    481492        'weight_product': "*".join(weights), 
    … …  
    484495        'qcall': q_pars['qcall'], 
    485496        'pcall': ", ".join(fq_pars), # skip scale and background 
     497        'sasview_spherical': sasview_spherical, 
    486498        } 
    487499    loop_body = [indent(LOOP_BODY%subst, depth)] 

  • sasmodels/gpu.py

    r5d4777d rf4cf580  
    2222devices, where it can be combined with other structure factors and form 
    2323factors and have instrumental resolution effects applied. 
    24  
    25  
    2624""" 
    2725import numpy as np 
    … …  
    5957    source, info = gen.make(kernel_module) 
    6058    ## for debugging, save source to a .cl file, edit it, and reload as model 
    61     open(info['name']+'.cl','w').write(source) 
     59    #open(info['name']+'.cl','w').write(source) 
    6260    #source = open(info['name']+'.cl','r').read() 
    6361    return GpuModel(source, info, dtype) 

  • sasmodels/models/capped_cylinder.c

    r5d4777d rf4cf580  
    2222    const real upper = REAL(1.0); 
    2323    const real lower = h/cap_radius; // integral lower bound 
    24     const real m = q*cos_alpha*cap_radius; // cos argument slope 
    25     const real b = q*cos_alpha*(REAL(0.5)*length-h); // cos argument intercept 
    26     const real qrst = q*sin_alpha*cap_radius; // Q*R*sin(theta) 
     24    // cos term in integral is: 
     25    //    cos (q (R t - h + L/2) cos(alpha)) 
     26    // so turn it into: 
     27    //    cos (m t + b) 
     28    // where: 
     29    //    m = q R cos(alpha) 
     30    //    b = q(L/2-h) cos(alpha) 
     31    const real m = q*cap_radius*cos_alpha; // cos argument slope 
     32    const real b = q*(REAL(0.5)*length-h)*cos_alpha; // cos argument intercept 
     33    const real qrst = q*cap_radius*sin_alpha; // Q*R*sin(theta) 
    2734    real total = REAL(0.0); 
    2835    for (int i=0; i<76 ;i++) { 
    … …  
    3138        const real t = REAL(0.5)*(Gauss76Z[i]*(upper-lower)+upper+lower); 
    3239        const real radical = REAL(1.0) - t*t; 
    33         const real caparg = qrst*sqrt(radical); // cap bessel function arg 
    34         const real be = (caparg == REAL(0.0) ? REAL(0.5) : J1(caparg)/caparg); 
     40        const real arg = qrst*sqrt(radical); // cap bessel function arg 
     41        const real be = (arg == REAL(0.0) ? REAL(0.5) : J1(arg)/arg); 
    3542        const real Fq = cos(m*t + b) * radical * be; 
    3643        total += Gauss76Wt[i] * Fq; 
    … …  
    8491        // faster since we don't multiply and divide sin(alpha). 
    8592        const real besarg = q*radius*sn; 
    86         const real siarg = REAL(0.5)*q*length*cn; 
     93        const real siarg = q*REAL(0.5)*length*cn; 
    8794        // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1 
    8895        const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg); 
    … …  
    131138    const real cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn); 
    132139 
    133     // The following is CylKernel() / sin(alpha), but we are doing it in place 
    134     // to avoid sin(alpha)/sin(alpha) for alpha = 0.  It is also a teensy bit 
    135     // faster since we don't multiply and divide sin(alpha). 
    136140    const real besarg = q*radius*sn; 
    137     const real siarg = REAL(0.5)*q*length*cn; 
     141    const real siarg = q*REAL(0.5)*length*cn; 
    138142    // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1 
    139143    const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg); 

  • sasmodels/models/capped_cylinder.py

    r5d4777d rf4cf580  
    11r""" 
    22Calculates the scattering from a cylinder with spherical section end-caps. 
    3 This model simply becomes the `convex_lens`_ when the length of the cylinder 
     3This model simply becomes the a convex lens when the length of the cylinder 
    44$L=0$, that is, a sphereocylinder with end caps that have a radius larger 
    55than that of the cylinder and the center of the end cap radius lies within 
    … …  
    128128    [ "radius", "Ang",  20, [0, inf], "volume", 
    129129      "Cylinder radius" ], 
     130    # TODO: use an expression for cap radius with fixed bounds. 
     131    # The current form requires cap radius R bigger than cylinder radius r. 
     132    # Could instead use R/r in [1,inf], r/R in [0,1], or the angle between 
     133    # cylinder and cap in [0,90].  The problem is similar for the barbell 
     134    # model.  Propose r/R in [0,1] in both cases, with the model specifying 
     135    # cylinder radius in the capped cylinder model and sphere radius in the 
     136    # barbell model.  This leads to the natural value of zero for no cap 
     137    # in the capped cylinder, and zero for no bar in the barbell model.  In 
     138    # both models, one would be a pill. 
    130139    [ "cap_radius", "Ang", 20, [0, inf], "volume", 
    131140      "Cap radius" ], 

  • sasmodels/models/sphere.py

    r5d4777d rf4cf580  
    3030 
    3131Iq = """ 
     32    const real qr = q*radius; 
    3233    real sn, cn; 
    33     const real qr = q*radius; 
    3434    SINCOS(qr, sn, cn); 
    35     const real bes = (qr==REAL(0.0) ? REAL(1.0) : (REAL(3.0)*(sn-qr*cn))/(qr*qr*qr)); 
    36     const real f = bes * (sld - solvent_sld) * form_volume(radius); 
    37     return REAL(1.0e-4)*f*f; 
     35    const real bes = qr==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(sn-qr*cn)/(qr*qr*qr); 
     36    const real fq = bes * (sld - solvent_sld) * form_volume(radius); 
     37    return REAL(1.0e-4)*fq*fq; 
    3838    """ 
    3939 
Note: See TracChangeset for help on using the changeset viewer.