Changeset f4cf580 in sasmodels


Ignore:
Timestamp:
Sep 2, 2014 1:15:40 PM (6 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.