Changes in / [bc1dac6:1cdfd47] in sasmodels


Ignore:
Files:
1 added
17 edited

Legend:

Unmodified
Added
Removed

  • .travis.yml

    r5c36bf1 r335271e  
    66    env: 
    77    - PY=2.7 
    8     #- DEPLOY=True 
     8    - DEPLOY=True 
    99  - os: linux 
    1010    env: 

  • doc/guide/magnetism/magnetism.rst

    rbefe905 r4f5afc9  
    3939 
    4040.. math:: 
    41     -- &= (1-u_i)(1-u_f) \\ 
    42     -+ &= (1-u_i)(u_f) \\ 
    43     +- &= (u_i)(1-u_f) \\ 
    44     ++ &= (u_i)(u_f) 
     41    -- &= ((1-u_i)(1-u_f))^{1/4} \\ 
     42    -+ &= ((1-u_i)(u_f))^{1/4} \\ 
     43    +- &= ((u_i)(1-u_f))^{1/4} \\ 
     44    ++ &= ((u_i)(u_f))^{1/4} 
    4545 
    4646Ideally the experiment would measure the pure spin states independently and 
    … …  
    104104| 2015-05-02 Steve King 
    105105| 2017-11-15 Paul Kienzle 
    106 | 2018-06-02 Adam Washington 

  • doc/guide/plugin.rst

    rf796469 r7e6bc45e  
    822822 
    823823        :code:`source = ["lib/Si.c", ...]` 
    824         (`Si.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_Si.c>`_) 
     824        (`Si.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/Si.c>`_) 
    825825 
    826826    sas_3j1x_x(x): 

  • sasmodels/compare.py

    r1e7b202a rd86f0fc  
    107107    -title="note" adds note to the plot title, after the model name 
    108108    -weights shows weights plots for the polydisperse parameters 
    109     -profile shows the sld profile if the model has a plottable sld profile 
    110109 
    111110    === output options === 
    … …  
    646645 
    647646def make_data(opts): 
    648     # type: (Dict[str, Any], float) -> Tuple[Data, np.ndarray] 
     647    # type: (Dict[str, Any]) -> Tuple[Data, np.ndarray] 
    649648    """ 
    650649    Generate an empty dataset, used with the model to set Q points 
    … …  
    668667        if opts['zero']: 
    669668            q = np.hstack((0, q)) 
    670         # TODO: provide command line control of lambda and Delta lambda/lambda 
    671         #L, dLoL = 5, 0.14/np.sqrt(6)  # wavelength and 14% triangular FWHM 
    672         L, dLoL = 0, 0 
    673         data = empty_data1D(q, resolution=res, L=L, dL=L*dLoL) 
     669        data = empty_data1D(q, resolution=res) 
    674670        index = slice(None, None) 
    675671    return data, index 
    … …  
    756752    """ 
    757753    for k in range(opts['sets']): 
    758         if k > 0: 
     754        if k > 1: 
    759755            # print a separate seed for each dataset for better reproducibility 
    760756            new_seed = np.random.randint(1000000) 
    761             print("=== Set %d uses -random=%i ==="%(k+1, new_seed)) 
     757            print("Set %d uses -random=%i"%(k+1, new_seed)) 
    762758            np.random.seed(new_seed) 
    763759        opts['pars'] = parse_pars(opts, maxdim=maxdim) 
    … …  
    766762        result = run_models(opts, verbose=True) 
    767763        if opts['plot']: 
    768             if opts['is2d'] and k > 0: 
    769                 import matplotlib.pyplot as plt 
    770                 plt.figure() 
    771764            limits = plot_models(opts, result, limits=limits, setnum=k) 
    772765        if opts['show_weights']: 
    … …  
    776769            dim = base._kernel.dim 
    777770            plot_weights(model_info, get_mesh(model_info, base_pars, dim=dim)) 
    778         if opts['show_profile']: 
    779             import pylab 
    780             base, comp = opts['engines'] 
    781             base_pars, comp_pars = opts['pars'] 
    782             have_base = base._kernel.info.profile is not None 
    783             have_comp = ( 
    784                 comp is not None 
    785                 and comp._kernel.info.profile is not None 
    786                 and base_pars != comp_pars 
    787             ) 
    788             if have_base or have_comp: 
    789                 pylab.figure() 
    790             if have_base: 
    791                 plot_profile(base._kernel.info, **base_pars) 
    792             if have_comp: 
    793                 plot_profile(comp._kernel.info, label='comp', **comp_pars) 
    794                 pylab.legend() 
    795771    if opts['plot']: 
    796772        import matplotlib.pyplot as plt 
    … …  
    798774    return limits 
    799775 
    800 def plot_profile(model_info, label='base', **args): 
    801     # type: (ModelInfo, List[Tuple[float, np.ndarray, np.ndarray]]) -> None 
    802     """ 
    803     Plot the profile returned by the model profile method. 
    804  
    805     *model_info* defines model parameters, etc. 
    806  
    807     *mesh* is a list of tuples containing (*value*, *dispersity*, *weights*) 
    808     for each parameter, where (*dispersity*, *weights*) pairs are the 
    809     distributions to be plotted. 
    810     """ 
    811     import pylab 
    812  
    813     args = dict((k, v) for k, v in args.items() 
    814                 if "_pd" not in k 
    815                    and ":" not in k 
    816                    and k not in ("background", "scale", "theta", "phi", "psi")) 
    817     args = args.copy() 
    818  
    819     args.pop('scale', 1.) 
    820     args.pop('background', 0.) 
    821     z, rho = model_info.profile(**args) 
    822     #pylab.interactive(True) 
    823     pylab.plot(z, rho, '-', label=label) 
    824     pylab.grid(True) 
    825     #pylab.show() 
    826  
    827  
    828  
    829776def run_models(opts, verbose=False): 
    830777    # type: (Dict[str, Any]) -> Dict[str, Any] 
    … …  
    836783    base_n, comp_n = opts['count'] 
    837784    base_pars, comp_pars = opts['pars'] 
    838     base_data, comp_data = opts['data'] 
     785    data = opts['data'] 
    839786 
    840787    comparison = comp is not None 
    … …  
    850797            print("%s t=%.2f ms, intensity=%.0f" 
    851798                  % (base.engine, base_time, base_value.sum())) 
    852         _show_invalid(base_data, base_value) 
     799        _show_invalid(data, base_value) 
    853800    except ImportError: 
    854801        traceback.print_exc() 
    … …  
    862809                print("%s t=%.2f ms, intensity=%.0f" 
    863810                      % (comp.engine, comp_time, comp_value.sum())) 
    864             _show_invalid(base_data, comp_value) 
     811            _show_invalid(data, comp_value) 
    865812        except ImportError: 
    866813            traceback.print_exc() 
    … …  
    916863    have_base, have_comp = (base_value is not None), (comp_value is not None) 
    917864    base, comp = opts['engines'] 
    918     base_data, comp_data = opts['data'] 
     865    data = opts['data'] 
    919866    use_data = (opts['datafile'] is not None) and (have_base ^ have_comp) 
    920867 
    921868    # Plot if requested 
    922869    view = opts['view'] 
    923     #view = 'log' 
    924870    if limits is None: 
    925871        vmin, vmax = np.inf, -np.inf 
    … …  
    935881        if have_comp: 
    936882            plt.subplot(131) 
    937         plot_theory(base_data, base_value, view=view, use_data=use_data, limits=limits) 
     883        plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
    938884        plt.title("%s t=%.2f ms"%(base.engine, base_time)) 
    939885        #cbar_title = "log I" 
    … …  
    942888            plt.subplot(132) 
    943889        if not opts['is2d'] and have_base: 
    944             plot_theory(comp_data, base_value, view=view, use_data=use_data, limits=limits) 
    945         plot_theory(comp_data, comp_value, view=view, use_data=use_data, limits=limits) 
     890            plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
     891        plot_theory(data, comp_value, view=view, use_data=use_data, limits=limits) 
    946892        plt.title("%s t=%.2f ms"%(comp.engine, comp_time)) 
    947893        #cbar_title = "log I" 
    … …  
    959905            err[err > cutoff] = cutoff 
    960906        #err,errstr = base/comp,"ratio" 
    961         # Note: base_data only since base and comp have same q values (though 
    962         # perhaps different resolution), and we are plotting the difference 
    963         # at each q 
    964         plot_theory(base_data, None, resid=err, view=errview, use_data=use_data) 
     907        plot_theory(data, None, resid=err, view=errview, use_data=use_data) 
    965908        plt.xscale('log' if view == 'log' and not opts['is2d'] else 'linear') 
    966909        plt.legend(['P%d'%(k+1) for k in range(setnum+1)], loc='best') 
    … …  
    996939OPTIONS = [ 
    997940    # Plotting 
    998     'plot', 'noplot', 
    999     'weights', 'profile', 
     941    'plot', 'noplot', 'weights', 
    1000942    'linear', 'log', 'q4', 
    1001943    'rel', 'abs', 
    … …  
    11301072        'qmax'      : 0.05, 
    11311073        'nq'        : 128, 
    1132         'res'       : '0.0', 
     1074        'res'       : 0.0, 
    11331075        'noise'     : 0.0, 
    11341076        'accuracy'  : 'Low', 
    … …  
    11511093        'count'     : '1', 
    11521094        'show_weights' : False, 
    1153         'show_profile' : False, 
    11541095        'sphere'    : 0, 
    11551096        'ngauss'    : '0', 
    … …  
    11711112        elif arg.startswith('-q='): 
    11721113            opts['qmin'], opts['qmax'] = [float(v) for v in arg[3:].split(':')] 
    1173         elif arg.startswith('-res='):      opts['res'] = arg[5:] 
     1114        elif arg.startswith('-res='):      opts['res'] = float(arg[5:]) 
    11741115        elif arg.startswith('-noise='):    opts['noise'] = float(arg[7:]) 
    11751116        elif arg.startswith('-sets='):     opts['sets'] = int(arg[6:]) 
    … …  
    12121153        elif arg == '-default': opts['use_demo'] = False 
    12131154        elif arg == '-weights': opts['show_weights'] = True 
    1214         elif arg == '-profile': opts['show_profile'] = True 
    12151155        elif arg == '-html':    opts['html'] = True 
    12161156        elif arg == '-help':    opts['html'] = True 
    … …  
    12301170    if opts['qmin'] is None: 
    12311171        opts['qmin'] = 0.001*opts['qmax'] 
     1172    if opts['datafile'] is not None: 
     1173        data = load_data(os.path.expanduser(opts['datafile'])) 
     1174    else: 
     1175        data, _ = make_data(opts) 
    12321176 
    12331177    comparison = any(PAR_SPLIT in v for v in values) 
    … …  
    12691213        opts['cutoff'] = [float(opts['cutoff'])]*2 
    12701214 
    1271     if PAR_SPLIT in opts['res']: 
    1272         opts['res'] = [float(k) for k in opts['res'].split(PAR_SPLIT, 2)] 
    1273         comparison = True 
    1274     else: 
    1275         opts['res'] = [float(opts['res'])]*2 
    1276  
    1277     if opts['datafile'] is not None: 
    1278         data = load_data(os.path.expanduser(opts['datafile'])) 
    1279     else: 
    1280         # Hack around the fact that make_data doesn't take a pair of resolutions 
    1281         res = opts['res'] 
    1282         opts['res'] = res[0] 
    1283         data0, _ = make_data(opts) 
    1284         if res[0] != res[1]: 
    1285             opts['res'] = res[1] 
    1286             data1, _ = make_data(opts) 
    1287         else: 
    1288             data1 = data0 
    1289         opts['res'] = res 
    1290         data = data0, data1 
    1291  
    1292     base = make_engine(model_info[0], data[0], opts['engine'][0], 
     1215    base = make_engine(model_info[0], data, opts['engine'][0], 
    12931216                       opts['cutoff'][0], opts['ngauss'][0]) 
    12941217    if comparison: 
    1295         comp = make_engine(model_info[1], data[1], opts['engine'][1], 
     1218        comp = make_engine(model_info[1], data, opts['engine'][1], 
    12961219                           opts['cutoff'][1], opts['ngauss'][1]) 
    12971220    else: 
    … …  
    14061329 
    14071330    # Evaluate preset parameter expressions 
    1408     # Note: need to replace ':' with '_' in parameter names and expressions 
    1409     # in order to support math on magnetic parameters. 
    14101331    context = MATH.copy() 
    14111332    context['np'] = np 
    1412     context.update((k.replace(':', '_'), v) for k, v in pars.items()) 
     1333    context.update(pars) 
    14131334    context.update((k, v) for k, v in presets.items() if isinstance(v, float)) 
    1414     #for k,v in sorted(context.items()): print(k, v) 
    14151335    for k, v in presets.items(): 
    14161336        if not isinstance(v, float) and not k.endswith('_type'): 
    1417             presets[k] = eval(v.replace(':', '_'), context) 
     1337            presets[k] = eval(v, context) 
    14181338    context.update(presets) 
    1419     context.update((k.replace(':', '_'), v) for k, v in presets2.items() if isinstance(v, float)) 
     1339    context.update((k, v) for k, v in presets2.items() if isinstance(v, float)) 
    14201340    for k, v in presets2.items(): 
    14211341        if not isinstance(v, float) and not k.endswith('_type'): 
    1422             presets2[k] = eval(v.replace(':', '_'), context) 
     1342            presets2[k] = eval(v, context) 
    14231343 
    14241344    # update parameters with presets 
    … …  
    14501370    path = os.path.dirname(info.filename) 
    14511371    url = "file://" + path.replace("\\", "/")[2:] + "/" 
    1452     rst2html.view_html_wxapp(html, url) 
     1372    rst2html.view_html_qtapp(html, url) 
    14531373 
    14541374def explore(opts): 

  • sasmodels/data.py

    r65fbf7c rd86f0fc  
    3636 
    3737import numpy as np  # type: ignore 
    38 from numpy import sqrt, sin, cos, pi 
    3938 
    4039# pylint: disable=unused-import 
    … …  
    302301 
    303302 
    304 def empty_data1D(q, resolution=0.0, L=0., dL=0.): 
     303def empty_data1D(q, resolution=0.0): 
    305304    # type: (np.ndarray, float) -> Data1D 
    306     r""" 
     305    """ 
    307306    Create empty 1D data using the given *q* as the x value. 
    308307 
    309     rms *resolution* $\Delta q/q$ defaults to 0%.  If wavelength *L* and rms 
    310     wavelength divergence *dL* are defined, then *resolution* defines 
    311     rms $\Delta \theta/\theta$ for the lowest *q*, with $\theta$ derived from 
    312     $q = 4\pi/\lambda \sin(\theta)$. 
     308    *resolution* dq/q defaults to 5%. 
    313309    """ 
    314310 
    … …  
    317313    Iq, dIq = None, None 
    318314    q = np.asarray(q) 
    319     if L != 0 and resolution != 0: 
    320         theta = np.arcsin(q*L/(4*pi)) 
    321         dtheta = theta[0]*resolution 
    322         ## Solving Gaussian error propagation from 
    323         ##   Dq^2 = (dq/dL)^2 DL^2 + (dq/dtheta)^2 Dtheta^2 
    324         ## gives 
    325         ##   (Dq/q)^2 = (DL/L)**2 + (Dtheta/tan(theta))**2 
    326         ## Take the square root and multiply by q, giving 
    327         ##   Dq = (4*pi/L) * sqrt((sin(theta)*dL/L)**2 + (cos(theta)*dtheta)**2) 
    328         dq = (4*pi/L) * sqrt((sin(theta)*dL/L)**2 + (cos(theta)*dtheta)**2) 
    329     else: 
    330         dq = resolution * q 
    331  
    332     data = Data1D(q, Iq, dx=dq, dy=dIq) 
     315    data = Data1D(q, Iq, dx=resolution * q, dy=dIq) 
    333316    data.filename = "fake data" 
    334317    return data 

  • sasmodels/details.py

    r885753a r108e70e  
    243243    offset = np.cumsum(np.hstack((0, length))) 
    244244    call_details = make_details(kernel.info, length, offset[:-1], offset[-1]) 
    245     # Pad value array to a 32 value boundary 
     245    # Pad value array to a 32 value boundaryd 
    246246    data_len = nvalues + 2*sum(len(v) for v in dispersity) 
    247247    extra = (32 - data_len%32)%32 
    … …  
    250250    is_magnetic = convert_magnetism(kernel.info.parameters, data) 
    251251    #call_details.show() 
    252     #print("data", data) 
    253252    return call_details, data, is_magnetic 
    254253 
    … …  
    297296    mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues] 
    298297    mag = mag.reshape(-1, 3) 
    299     if np.any(mag[:, 0] != 0.0): 
    300         M0 = mag[:, 0].copy() 
     298    scale = mag[:, 0] 
     299    if np.any(scale): 
    301300        theta, phi = radians(mag[:, 1]), radians(mag[:, 2]) 
    302         mag[:, 0] = +M0*cos(theta)*cos(phi)  # mx 
    303         mag[:, 1] = +M0*sin(theta) # my 
    304         mag[:, 2] = -M0*cos(theta)*sin(phi)  # mz 
     301        cos_theta = cos(theta) 
     302        mag[:, 0] = scale*cos_theta*cos(phi)  # mx 
     303        mag[:, 1] = scale*sin(theta)  # my 
     304        mag[:, 2] = -scale*cos_theta*sin(phi)  # mz 
    305305        return True 
    306306    else: 

  • sasmodels/jitter.py

    r1198f90 rb3703f5  
    774774        # set small jitter as 0 if multiple pd dims 
    775775        dims = sum(v > 0 for v in jitter) 
    776         limit = [0, 0.5, 5][dims] 
     776        limit = [0, 0, 0.5, 5][dims] 
    777777        jitter = [0 if v < limit else v for v in jitter] 
    778778        axes.cla() 

  • sasmodels/kernel_iq.c

    r7c35fda rd86f0fc  
    7979//     du * (m_sigma_y + 1j*m_sigma_z); 
    8080// weights for spin crosssections: dd du real, ud real, uu, du imag, ud imag 
    81 static void set_spin_weights(double in_spin, double out_spin, double weight[6]) 
     81static void set_spin_weights(double in_spin, double out_spin, double spins[4]) 
    8282{ 
    8383  in_spin = clip(in_spin, 0.0, 1.0); 
    8484  out_spin = clip(out_spin, 0.0, 1.0); 
    85   // Previous version of this function took the square root of the weights, 
    86   // under the assumption that  
    87   // 
    88   //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...) 
    89   // 
    90   // However, since the weights are applied to the final intensity and 
    91   // are not interned inside the I(q) function, we want the full 
    92   // weight and not the square root.  Any function using 
    93   // set_spin_weights as part of calculating an amplitude will need to 
    94   // manually take that square root, but there is currently no such 
    95   // function. 
    96   weight[0] = (1.0-in_spin) * (1.0-out_spin); // dd 
    97   weight[1] = (1.0-in_spin) * out_spin;       // du 
    98   weight[2] = in_spin * (1.0-out_spin);       // ud 
    99   weight[3] = in_spin * out_spin;             // uu 
    100   weight[4] = weight[1]; // du.imag 
    101   weight[5] = weight[2]; // ud.imag 
     85  spins[0] = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); // dd 
     86  spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du real 
     87  spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud real 
     88  spins[3] = sqrt(sqrt(in_spin * out_spin));             // uu 
     89  spins[4] = spins[1]; // du imag 
     90  spins[5] = spins[2]; // ud imag 
    10291} 
    10392 
    10493// Compute the magnetic sld 
    10594static double mag_sld( 
    106   const unsigned int xs, // 0=dd, 1=du.real, 2=ud.real, 3=uu, 4=du.imag, 5=ud.imag 
     95  const unsigned int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag 
    10796  const double qx, const double qy, 
    10897  const double px, const double py, 
    … …  
    117106      case 0: // uu => sld - D M_perpx 
    118107          return sld - px*perp; 
    119       case 1: // ud.real => -D M_perpy 
     108      case 1: // ud real => -D M_perpy 
    120109          return py*perp; 
    121       case 2: // du.real => -D M_perpy 
     110      case 2: // du real => -D M_perpy 
    122111          return py*perp; 
    123       case 3: // dd => sld + D M_perpx 
     112      case 3: // dd real => sld + D M_perpx 
    124113          return sld + px*perp; 
    125114    } 
    126115  } else { 
    127116    if (xs== 4) { 
    128       return -mz;  // du.imag => +D M_perpz 
     117      return -mz;  // ud imag => -D M_perpz 
    129118    } else { // index == 5 
    130       return +mz;  // ud.imag => -D M_perpz 
     119      return mz;   // du imag => D M_perpz 
    131120    } 
    132121  } 
    … …  
    323312  //     up_angle = values[NUM_PARS+4]; 
    324313  // TODO: could precompute more magnetism parameters before calling the kernel. 
    325   double xs_weights[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill 
     314  double spins[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill 
    326315  double cos_mspin, sin_mspin; 
    327   set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], xs_weights); 
     316  set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], spins); 
    328317  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin); 
    329318#endif // MAGNETIC 
    … …  
    672661            // loop over uu, ud real, du real, dd, ud imag, du imag 
    673662            for (unsigned int xs=0; xs<6; xs++) { 
    674               const double xs_weight = xs_weights[xs]; 
     663              const double xs_weight = spins[xs]; 
    675664              if (xs_weight > 1.e-8) { 
    676665                // Since the cross section weight is significant, set the slds 
    … …  
    685674                  local_values.vector[sld_index] = 
    686675                    mag_sld(xs, qx, qy, px, py, values[sld_index+2], mx, my, mz); 
    687 //if (q_index==0) printf("%d: (qx,qy)=(%g,%g) xs=%d sld%d=%g p=(%g,%g) m=(%g,%g,%g)\n", 
    688 //  q_index, qx, qy, xs, sk, local_values.vector[sld_index], px, py, mx, my, mz); 
    689676                } 
    690677                scattering += xs_weight * CALL_KERNEL(); 

  • sasmodels/kerneldll.py

    r33969b6 rbf94e6e  
    123123    # add openmp support if not running on a mac 
    124124    if sys.platform != "darwin": 
    125         # OpenMP seems to be broken on gcc 5.4.0 (ubuntu 16.04.9) 
    126         # Shut it off for all unix until we can investigate. 
    127         #CC.append("-fopenmp") 
    128         pass 
     125        CC.append("-fopenmp") 
    129126    def compile_command(source, output): 
    130127        """unix compiler command""" 

  • sasmodels/models/core_shell_parallelepiped.c

    rdbf1a60 re077231  
    5959 
    6060    // outer integral (with gauss points), integration limits = 0, 1 
    61     // substitute d_cos_alpha for sin_alpha d_alpha 
    6261    double outer_sum = 0; //initialize integral 
    6362    for( int i=0; i<GAUSS_N; i++) { 
    6463        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
    6564        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 
     65 
     66        // inner integral (with gauss points), integration limits = 0, pi/2 
    6667        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 
    6768        const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q); 
    68  
    69         // inner integral (with gauss points), integration limits = 0, 1 
    70         // substitute beta = PI/2 u (so 2/PI * d_(PI/2 * beta) = d_beta) 
    7169        double inner_sum = 0.0; 
    7270        for(int j=0; j<GAUSS_N; j++) { 
    73             const double u = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
     71            const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
    7472            double sin_beta, cos_beta; 
    75             SINCOS(M_PI_2*u, sin_beta, cos_beta); 
     73            SINCOS(M_PI_2*beta, sin_beta, cos_beta); 
    7674            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 
    7775            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); 
    … …  
    9391            inner_sum += GAUSS_W[j] * f * f; 
    9492        } 
    95         // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    9693        inner_sum *= 0.5; 
    9794        // now sum up the outer integral 
    9895        outer_sum += GAUSS_W[i] * inner_sum; 
    9996    } 
    100     // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    10197    outer_sum *= 0.5; 
    10298 

  • sasmodels/models/core_shell_parallelepiped.py

    rf89ec96 r97be877  
    44 
    55Calculates the form factor for a rectangular solid with a core-shell structure. 
    6 The thickness and the scattering length density of the shell or "rim" can be 
    7 different on each (pair) of faces. The three dimensions of the core of the 
    8 parallelepiped (strictly here a cuboid) may be given in *any* size order as 
    9 long as the particles are randomly oriented (i.e. take on all possible 
    10 orientations see notes on 2D below). To avoid multiple fit solutions, 
    11 especially with Monte-Carlo fit methods, it may be advisable to restrict their 
    12 ranges. There may be a number of closely similar "best fits", so some trial and 
    13 error, or fixing of some dimensions at expected values, may help. 
     6The thickness and the scattering length density of the shell or 
     7"rim" can be different on each (pair) of faces. 
    148 
    159The form factor is normalized by the particle volume $V$ such that 
    … …  
    1711.. math:: 
    1812 
    19     I(q) = \frac{\text{scale}}{V} \langle P(q,\alpha,\beta) \rangle 
    20     + \text{background} 
     13    I(q) = \text{scale}\frac{\langle f^2 \rangle}{V} + \text{background} 
    2114 
    2215where $\langle \ldots \rangle$ is an average over all possible orientations 
    23 of the rectangular solid, and the usual $\Delta \rho^2 \ V^2$ term cannot be 
    24 pulled out of the form factor term due to the multiple slds in the model. 
    25  
    26 The core of the solid is defined by the dimensions $A$, $B$, $C$ here shown 
    27 such that $A < B < C$. 
    28  
    29 .. figure:: img/parallelepiped_geometry.jpg 
    30  
    31    Core of the core shell parallelepiped with the corresponding definition 
    32    of sides. 
    33  
     16of the rectangular solid. 
     17 
     18The function calculated is the form factor of the rectangular solid below. 
     19The core of the solid is defined by the dimensions $A$, $B$, $C$ such that 
     20$A < B < C$. 
     21 
     22.. image:: img/core_shell_parallelepiped_geometry.jpg 
    3423 
    3524There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension 
    3625(on the $BC$ faces). There are similar slabs on the $AC$ $(=t_B)$ and $AB$ 
    37 $(=t_C)$ faces. The projection in the $AB$ plane is 
    38  
    39 .. figure:: img/core_shell_parallelepiped_projection.jpg 
    40  
    41    AB cut through the core-shell parallelipiped showing the cross secion of 
    42    four of the six shell slabs. As can be seen, this model leaves **"gaps"** 
    43    at the corners of the solid. 
    44  
    45  
    46 The total volume of the solid is thus given as 
     26$(=t_C)$ faces. The projection in the $AB$ plane is then 
     27 
     28.. image:: img/core_shell_parallelepiped_projection.jpg 
     29 
     30The volume of the solid is 
    4731 
    4832.. math:: 
    4933 
    5034    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
     35 
     36**meaning that there are "gaps" at the corners of the solid.** 
    5137 
    5238The intensity calculated follows the :ref:`parallelepiped` model, with the 
    5339core-shell intensity being calculated as the square of the sum of the 
    54 amplitudes of the core and the slabs on the edges. The scattering amplitude is 
    55 computed for a particular orientation of the core-shell parallelepiped with 
    56 respect to the scattering vector and then averaged over all possible 
    57 orientations, where $\alpha$ is the angle between the $z$ axis and the $C$ axis 
    58 of the parallelepiped, and $\beta$ is the angle between the projection of the 
    59 particle in the $xy$ detector plane and the $y$ axis. 
    60  
    61 .. math:: 
    62  
    63     P(q)=\frac {\int_{0}^{\pi/2}\int_{0}^{\pi/2}F^2(q,\alpha,\beta) \ sin\alpha 
    64     \ d\alpha \ d\beta} {\int_{0}^{\pi/2} \ sin\alpha \ d\alpha \ d\beta} 
    65  
    66 and 
    67  
    68 .. math:: 
    69  
    70     F(q,\alpha,\beta) 
     40amplitudes of the core and the slabs on the edges. 
     41 
     42the scattering amplitude is computed for a particular orientation of the 
     43core-shell parallelepiped with respect to the scattering vector and then 
     44averaged over all possible orientations, where $\alpha$ is the angle between 
     45the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is 
     46the angle between projection of the particle in the $xy$ detector plane 
     47and the $y$ axis. 
     48 
     49.. math:: 
     50 
     51    F(Q) 
    7152    &= (\rho_\text{core}-\rho_\text{solvent}) 
    7253       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 
    7354    &+ (\rho_\text{A}-\rho_\text{solvent}) 
    74         \left[S(Q_A, A+2t_A) - S(Q_A, A)\right] S(Q_B, B) S(Q_C, C) \\ 
     55        \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\ 
    7556    &+ (\rho_\text{B}-\rho_\text{solvent}) 
    7657        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ 
    … …  
    8263.. math:: 
    8364 
    84     S(Q_X, L) = L \frac{\sin (\tfrac{1}{2} Q_X L)}{\tfrac{1}{2} Q_X L} 
     65    S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} 
    8566 
    8667and 
    … …  
    8869.. math:: 
    8970 
    90     Q_A &= q \sin\alpha \sin\beta \\ 
    91     Q_B &= q \sin\alpha \cos\beta \\ 
    92     Q_C &= q \cos\alpha 
     71    Q_A &= \sin\alpha \sin\beta \\ 
     72    Q_B &= \sin\alpha \cos\beta \\ 
     73    Q_C &= \cos\alpha 
    9374 
    9475 
    9576where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ 
    96 are the scattering lengths of the parallelepiped core, and the rectangular 
     77are the scattering length of the parallelepiped core, and the rectangular 
    9778slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 
    9879is the scattering length of the solvent. 
    9980 
    100 .. note::  
    101  
    102    the code actually implements two substitutions: $d(cos\alpha)$ is 
    103    substituted for -$sin\alpha \ d\alpha$ (note that in the 
    104    :ref:`parallelepiped` code this is explicitly implemented with 
    105    $\sigma = cos\alpha$), and $\beta$ is set to $\beta = u \pi/2$ so that 
    106    $du = \pi/2 \ d\beta$.  Thus both integrals go from 0 to 1 rather than 0 
    107    to $\pi/2$. 
    108  
    10981FITTING NOTES 
    11082~~~~~~~~~~~~~ 
    11183 
    112 #. There are many parameters in this model. Hold as many fixed as possible with 
    113    known values, or you will certainly end up at a solution that is unphysical. 
    114  
    115 #. The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
    116    based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    117    and length $(C+2t_C)$ values, after appropriately sorting the three 
    118    dimensions to give an oblate or prolate particle, to give an effective radius 
    119    for $S(q)$ when $P(q) * S(q)$ is applied. 
    120  
    121 #. For 2d data the orientation of the particle is required, described using 
    122    angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, where $\theta$ 
    123    and $\phi$ define the orientation of the director in the laboratry reference 
    124    frame of the beam direction ($z$) and detector plane ($x-y$ plane), while 
    125    the angle $\Psi$ is effectively the rotational angle around the particle 
    126    $C$ axis. For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the 
    127    $B$ axis oriented parallel to the y-axis of the detector with $A$ along 
    128    the x-axis. For other $\theta$, $\phi$ values, the order of rotations 
    129    matters. In particular, the parallelepiped must first be rotated $\theta$ 
    130    degrees in the $x-z$ plane before rotating $\phi$ degrees around the $z$ 
    131    axis (in the $x-y$ plane). Applying orientational distribution to the 
    132    particle orientation (i.e  `jitter` to one or more of these angles) can get 
    133    more confusing as `jitter` is defined **NOT** with respect to the laboratory 
    134    frame but the particle reference frame. It is thus highly recmmended to 
    135    read :ref:`orientation` for further details of the calculation and angular 
    136    dispersions. 
    137  
    138 .. note:: For 2d, constraints must be applied during fitting to ensure that the 
    139    order of sides chosen is not altered, and hence that the correct definition 
    140    of angles is preserved. For the default choice shown here, that means 
    141    ensuring that the inequality $A < B < C$ is not violated,  The calculation 
    142    will not report an error, but the results may be not correct. 
     84If the scale is set equal to the particle volume fraction, $\phi$, the returned 
     85value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However, 
     86**no interparticle interference effects are included in this calculation.** 
     87 
     88There are many parameters in this model. Hold as many fixed as possible with 
     89known values, or you will certainly end up at a solution that is unphysical. 
     90 
     91The returned value is in units of |cm^-1|, on absolute scale. 
     92 
     93NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
     94based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
     95and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 
     96to give an oblate or prolate particle, to give an effective radius, 
     97for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     98 
     99For 2d data the orientation of the particle is required, described using 
     100angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further 
     101details of the calculation and angular dispersions see :ref:`orientation`. 
     102The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, 
     103$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 
     104 
     105For 2d, constraints must be applied during fitting to ensure that the 
     106inequality $A < B < C$ is not violated, and hence the correct definition 
     107of angles is preserved. The calculation will not report an error, 
     108but the results may be not correct. 
    143109 
    144110.. figure:: img/parallelepiped_angle_definition.png 
    145111 
    146112    Definition of the angles for oriented core-shell parallelepipeds. 
    147     Note that rotation $\theta$, initially in the $x-z$ plane, is carried 
     113    Note that rotation $\theta$, initially in the $xz$ plane, is carried 
    148114    out first, then rotation $\phi$ about the $z$ axis, finally rotation 
    149     $\Psi$ is now around the $C$ axis of the particle. The neutron or X-ray 
    150     beam is along the $z$ axis and the detecotr defines the $x-y$ plane. 
     115    $\Psi$ is now around the axis of the cylinder. The neutron or X-ray 
     116    beam is along the $z$ axis. 
    151117 
    152118.. figure:: img/parallelepiped_angle_projection.png 
    … …  
    154120    Examples of the angles for oriented core-shell parallelepipeds against the 
    155121    detector plane. 
    156  
    157  
    158 Validation 
    159 ---------- 
    160  
    161 Cross-checked against hollow rectangular prism and rectangular prism for equal 
    162 thickness overlapping sides, and by Monte Carlo sampling of points within the 
    163 shape for non-uniform, non-overlapping sides. 
    164  
    165122 
    166123References 
    … …  
    178135 
    179136* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    180 * **Converted to sasmodels by:** Miguel Gonzalez **Date:** February 26, 2016 
     137* **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016 
    181138* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017 
    182 * **Last Reviewed by:** Paul Butler **Date:** May 24, 2018 - documentation 
    183   updated 
     139* Cross-checked against hollow rectangular prism and rectangular prism for 
     140  equal thickness overlapping sides, and by Monte Carlo sampling of points 
     141  within the shape for non-uniform, non-overlapping sides. 
    184142""" 
    185143 

  • sasmodels/models/mass_surface_fractal.py

    r7994359 r2d81cfe  
    3939    The surface ( $D_s$ ) and mass ( $D_m$ ) fractal dimensions are only 
    4040    valid if $0 < surface\_dim < 6$ , $0 < mass\_dim < 6$ , and 
    41     $(surface\_dim + mass\_dim ) < 6$ .  
    42     Older versions of sasview may have the default primary particle radius 
    43     larger than the cluster radius, this was an error, also present in the  
    44     Schmidt review paper below. The primary particle should be the smaller  
    45     as described in the original Hurd et.al. who also point out that  
    46     polydispersity in the primary particle sizes may affect their  
    47     apparent surface fractal dimension. 
    48      
     41    $(surface\_dim + mass\_dim ) < 6$ . 
     42 
    4943 
    5044References 
    5145---------- 
    5246 
    53 .. [#] P Schmidt, *J Appl. Cryst.*, 24 (1991) 414-435 Equation(19) 
    54 .. [#] A J Hurd, D W Schaefer, J E Martin, *Phys. Rev. A*, 
    55    35 (1987) 2361-2364 Equation(2) 
     47P Schmidt, *J Appl. Cryst.*, 24 (1991) 414-435 Equation(19) 
    5648 
    57 Authorship and Verification 
    58 ---------------------------- 
    59  
    60 * **Converted to sasmodels by:** Piotr Rozyczko **Date:** Jan 20, 2016 
    61 * **Last Reviewed by:** Richard Heenan **Date:** May 30, 2018 
     49A J Hurd, D W Schaefer, J E Martin, *Phys. Rev. A*, 
     5035 (1987) 2361-2364 Equation(2) 
    6251""" 
    6352 
    … …  
    7867        rg_primary    =  rg 
    7968        background   =  background 
     69        Ref: Schmidt, J Appl Cryst, eq(19), (1991), 24, 414-435 
    8070        Hurd, Schaefer, Martin, Phys Rev A, eq(2),(1987),35, 2361-2364 
    8171        Note that 0 < Ds< 6 and 0 < Dm < 6. 
    … …  
    8878    ["fractal_dim_mass", "",      1.8, [0.0, 6.0], "", "Mass fractal dimension"], 
    8979    ["fractal_dim_surf", "",      2.3, [0.0, 6.0], "", "Surface fractal dimension"], 
    90     ["rg_cluster",       "Ang", 4000., [0.0, inf], "", "Cluster radius of gyration"], 
    91     ["rg_primary",       "Ang",  86.7, [0.0, inf], "", "Primary particle radius of gyration"], 
     80    ["rg_cluster",       "Ang",  86.7, [0.0, inf], "", "Cluster radius of gyration"], 
     81    ["rg_primary",       "Ang", 4000., [0.0, inf], "", "Primary particle radius of gyration"], 
    9282] 
    9383# pylint: enable=bad-whitespace, line-too-long 
    … …  
    117107            fractal_dim_mass=1.8, 
    118108            fractal_dim_surf=2.3, 
    119             rg_cluster=4000.0, 
    120             rg_primary=86.7) 
     109            rg_cluster=86.7, 
     110            rg_primary=4000.0) 
    121111 
    122112tests = [ 
    123113 
    124     # Accuracy tests based on content in test/utest_other_models.py  All except first, changed so rg_cluster is the larger, RKH 30 May 2018 
    125     [{'fractal_dim_mass':   1.8, 
     114    # Accuracy tests based on content in test/utest_other_models.py 
     115    [{'fractal_dim_mass':      1.8, 
    126116      'fractal_dim_surf':   2.3, 
    127117      'rg_cluster':   86.7, 
    … …  
    133123    [{'fractal_dim_mass':      3.3, 
    134124      'fractal_dim_surf':   1.0, 
    135       'rg_cluster': 4000.0, 
    136       'rg_primary':   90.0, 
    137      }, 0.001, 0.0932516614456], 
     125      'rg_cluster':   90.0, 
     126      'rg_primary': 4000.0, 
     127     }, 0.001, 0.18562699016], 
    138128 
    139129    [{'fractal_dim_mass':      1.3, 
    140       'fractal_dim_surf':   2.0, 
    141       'rg_cluster': 2000.0, 
    142       'rg_primary':   90.0, 
     130      'fractal_dim_surf':   1.0, 
     131      'rg_cluster':   90.0, 
     132      'rg_primary': 2000.0, 
    143133      'background':    0.8, 
    144      }, 0.001, 1.28296431786], 
     134     }, 0.001, 1.16539753641], 
    145135 
    146136    [{'fractal_dim_mass':      2.3, 
    147       'fractal_dim_surf':   3.1, 
    148       'rg_cluster':  1000.0, 
    149       'rg_primary':  30.0, 
     137      'fractal_dim_surf':   1.0, 
     138      'rg_cluster':   90.0, 
     139      'rg_primary': 1000.0, 
    150140      'scale':        10.0, 
    151141      'background':    0.0, 
    152      }, 0.051, 0.00333804044899], 
     142     }, 0.051, 0.000169548800377], 
    153143    ] 

  • sasmodels/models/parallelepiped.c

    rdbf1a60 r108e70e  
    3838            inner_total += GAUSS_W[j] * square(si1 * si2); 
    3939        } 
    40         // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    4140        inner_total *= 0.5; 
    4241 
    … …  
    4443        outer_total += GAUSS_W[i] * inner_total * si * si; 
    4544    } 
    46     // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    4745    outer_total *= 0.5; 
    4846 

  • sasmodels/models/parallelepiped.py

    rf89ec96 ref07e95  
    22# Note: model title and parameter table are inserted automatically 
    33r""" 
     4The form factor is normalized by the particle volume. 
     5For information about polarised and magnetic scattering, see 
     6the :ref:`magnetism` documentation. 
     7 
    48Definition 
    59---------- 
    610 
    7 This model calculates the scattering from a rectangular solid 
    8 (:numref:`parallelepiped-image`). 
    9 If you need to apply polydispersity, see also :ref:`rectangular-prism`. For 
    10 information about polarised and magnetic scattering, see 
    11 the :ref:`magnetism` documentation. 
     11 This model calculates the scattering from a rectangular parallelepiped 
     12 (\:numref:`parallelepiped-image`\). 
     13 If you need to apply polydispersity, see also :ref:`rectangular-prism`. 
    1214 
    1315.. _parallelepiped-image: 
    … …  
    1921 
    2022The three dimensions of the parallelepiped (strictly here a cuboid) may be 
    21 given in *any* size order as long as the particles are randomly oriented (i.e. 
    22 take on all possible orientations see notes on 2D below). To avoid multiple fit 
    23 solutions, especially with Monte-Carlo fit methods, it may be advisable to 
    24 restrict their ranges. There may be a number of closely similar "best fits", so 
    25 some trial and error, or fixing of some dimensions at expected values, may 
    26 help. 
    27  
    28 The form factor is normalized by the particle volume and the 1D scattering 
    29 intensity $I(q)$ is then calculated as: 
     23given in *any* size order. To avoid multiple fit solutions, especially 
     24with Monte-Carlo fit methods, it may be advisable to restrict their ranges. 
     25There may be a number of closely similar "best fits", so some trial and 
     26error, or fixing of some dimensions at expected values, may help. 
     27 
     28The 1D scattering intensity $I(q)$ is calculated as: 
    3029 
    3130.. Comment by Miguel Gonzalez: 
    … …  
    4039 
    4140    I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2 
    42            \left< P(q, \alpha, \beta) \right> + \text{background} 
     41           \left< P(q, \alpha) \right> + \text{background} 
    4342 
    4443where the volume $V = A B C$, the contrast is defined as 
    45 $\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, $P(q, \alpha, \beta)$ 
    46 is the form factor corresponding to a parallelepiped oriented 
    47 at an angle $\alpha$ (angle between the long axis C and $\vec q$), and $\beta$ 
    48 (the angle between the projection of the particle in the $xy$ detector plane 
    49 and the $y$ axis) and the averaging $\left<\ldots\right>$ is applied over all 
    50 orientations. 
     44$\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, 
     45$P(q, \alpha)$ is the form factor corresponding to a parallelepiped oriented 
     46at an angle $\alpha$ (angle between the long axis C and $\vec q$), 
     47and the averaging $\left<\ldots\right>$ is applied over all orientations. 
    5148 
    5249Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the 
    53 form factor is given by (Mittelbach and Porod, 1961 [#Mittelbach]_) 
     50form factor is given by (Mittelbach and Porod, 1961) 
    5451 
    5552.. math:: 
    … …  
    6966    \mu &= qB 
    7067 
    71 where substitution of $\sigma = cos\alpha$ and $\beta = \pi/2 \ u$ have been 
    72 applied. 
    73  
    74 For **oriented** particles, the 2D scattering intensity, $I(q_x, q_y)$, is 
    75 given as: 
    76  
    77 .. math:: 
    78  
    79     I(q_x, q_y) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2 P(q_x, q_y) 
     68The scattering intensity per unit volume is returned in units of |cm^-1|. 
     69 
     70NB: The 2nd virial coefficient of the parallelepiped is calculated based on 
     71the averaged effective radius, after appropriately sorting the three 
     72dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and 
     73length $(= C)$ values, and used as the effective radius for 
     74$S(q)$ when $P(q) \cdot S(q)$ is applied. 
     75 
     76For 2d data the orientation of the particle is required, described using 
     77angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
     78of the calculation and angular dispersions see :ref:`orientation` . 
     79 
     80.. Comment by Miguel Gonzalez: 
     81   The following text has been commented because I think there are two 
     82   mistakes. Psi is the rotational angle around C (but I cannot understand 
     83   what it means against the q plane) and psi=0 corresponds to a||x and b||y. 
     84 
     85   The angle $\Psi$ is the rotational angle around the $C$ axis against 
     86   the $q$ plane. For example, $\Psi = 0$ when the $B$ axis is parallel 
     87   to the $x$-axis of the detector. 
     88 
     89The angle $\Psi$ is the rotational angle around the $C$ axis. 
     90For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the $B$ axis 
     91oriented parallel to the y-axis of the detector with $A$ along the x-axis. 
     92For other $\theta$, $\phi$ values, the parallelepiped has to be first rotated 
     93$\theta$ degrees in the $z-x$ plane and then $\phi$ degrees around the $z$ axis, 
     94before doing a final rotation of $\Psi$ degrees around the resulting $C$ axis 
     95of the particle to obtain the final orientation of the parallelepiped. 
     96 
     97.. _parallelepiped-orientation: 
     98 
     99.. figure:: img/parallelepiped_angle_definition.png 
     100 
     101    Definition of the angles for oriented parallelepiped, shown with $A<B<C$. 
     102 
     103.. figure:: img/parallelepiped_angle_projection.png 
     104 
     105    Examples of the angles for an oriented parallelepiped against the 
     106    detector plane. 
     107 
     108On introducing "Orientational Distribution" in the angles, "distribution of 
     109theta" and "distribution of phi" parameters will appear. These are actually 
     110rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped, 
     111perpendicular to the $a$ x $c$ and $b$ x $c$ faces. (When $\theta = \phi = 0$ 
     112these are parallel to the $Y$ and $X$ axes of the instrument.) The third 
     113orientation distribution, in $\psi$, is about the $c$ axis of the particle, 
     114perpendicular to the $a$ x $b$ face. Some experimentation may be required to 
     115understand the 2d patterns fully as discussed in :ref:`orientation` . 
     116 
     117For a given orientation of the parallelepiped, the 2D form factor is 
     118calculated as 
     119 
     120.. math:: 
     121 
     122    P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}{2}qA\cos\alpha)}\right]^2 
     123                  \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}{2}qB\cos\beta)}\right]^2 
     124                  \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}{2}qC\cos\gamma)}\right]^2 
     125 
     126with 
     127 
     128.. math:: 
     129 
     130    \cos\alpha &= \hat A \cdot \hat q, \\ 
     131    \cos\beta  &= \hat B \cdot \hat q, \\ 
     132    \cos\gamma &= \hat C \cdot \hat q 
     133 
     134and the scattering intensity as: 
     135 
     136.. math:: 
     137 
     138    I(q_x, q_y) = \frac{\text{scale}}{V} V^2 \Delta\rho^2 P(q_x, q_y) 
    80139            + \text{background} 
    81140 
    … …  
    89148   with scale being the volume fraction. 
    90149 
    91 Where $P(q_x, q_y)$ for a given orientation of the form factor is calculated as 
    92  
    93 .. math:: 
    94  
    95     P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1} 
    96                    {2}qA\cos\alpha)}\right]^2 
    97                   \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1} 
    98                    {2}qB\cos\beta)}\right]^2 
    99                   \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1} 
    100                    {2}qC\cos\gamma)}\right]^2 
    101  
    102 with 
    103  
    104 .. math:: 
    105  
    106     \cos\alpha &= \hat A \cdot \hat q, \\ 
    107     \cos\beta  &= \hat B \cdot \hat q, \\ 
    108     \cos\gamma &= \hat C \cdot \hat q 
    109  
    110  
    111 FITTING NOTES 
    112 ~~~~~~~~~~~~~ 
    113  
    114 #. The 2nd virial coefficient of the parallelepiped is calculated based on 
    115    the averaged effective radius, after appropriately sorting the three 
    116    dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and 
    117    length $(= C)$ values, and used as the effective radius for 
    118    $S(q)$ when $P(q) \cdot S(q)$ is applied. 
    119  
    120 #. For 2d data the orientation of the particle is required, described using 
    121    angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, where $\theta$ 
    122    and $\phi$ define the orientation of the director in the laboratry reference 
    123    frame of the beam direction ($z$) and detector plane ($x-y$ plane), while 
    124    the angle $\Psi$ is effectively the rotational angle around the particle 
    125    $C$ axis. For $\theta = 0$ and $\phi = 0$, $\Psi = 0$ corresponds to the 
    126    $B$ axis oriented parallel to the y-axis of the detector with $A$ along 
    127    the x-axis. For other $\theta$, $\phi$ values, the order of rotations 
    128    matters. In particular, the parallelepiped must first be rotated $\theta$ 
    129    degrees in the $x-z$ plane before rotating $\phi$ degrees around the $z$ 
    130    axis (in the $x-y$ plane). Applying orientational distribution to the 
    131    particle orientation (i.e  `jitter` to one or more of these angles) can get 
    132    more confusing as `jitter` is defined **NOT** with respect to the laboratory 
    133    frame but the particle reference frame. It is thus highly recmmended to 
    134    read :ref:`orientation` for further details of the calculation and angular 
    135    dispersions. 
    136  
    137 .. note:: For 2d, constraints must be applied during fitting to ensure that the 
    138    order of sides chosen is not altered, and hence that the correct definition 
    139    of angles is preserved. For the default choice shown here, that means 
    140    ensuring that the inequality $A < B < C$ is not violated,  The calculation 
    141    will not report an error, but the results may be not correct. 
    142     
    143 .. _parallelepiped-orientation: 
    144  
    145 .. figure:: img/parallelepiped_angle_definition.png 
    146  
    147     Definition of the angles for oriented parallelepiped, shown with $A<B<C$. 
    148  
    149 .. figure:: img/parallelepiped_angle_projection.png 
    150  
    151     Examples of the angles for an oriented parallelepiped against the 
    152     detector plane. 
    153  
    154 .. Comment by Paul Butler 
    155    I am commenting this section out as we are trying to minimize the amount of 
    156    oritentational detail here and encourage the user to go to the full 
    157    orientation documentation so that changes can be made in just one place. 
    158    below is the commented paragrah: 
    159    On introducing "Orientational Distribution" in the angles, "distribution of 
    160    theta" and "distribution of phi" parameters will appear. These are actually 
    161    rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped, 
    162    perpendicular to the $a$ x $c$ and $b$ x $c$ faces. (When $\theta = \phi = 0$ 
    163    these are parallel to the $Y$ and $X$ axes of the instrument.) The third 
    164    orientation distribution, in $\psi$, is about the $c$ axis of the particle, 
    165    perpendicular to the $a$ x $b$ face. Some experimentation may be required to 
    166    understand the 2d patterns fully as discussed in :ref:`orientation` . 
    167  
    168150 
    169151Validation 
    … …  
    174156angles. 
    175157 
     158 
    176159References 
    177160---------- 
    178161 
    179 .. [#Mittelbach] P Mittelbach and G Porod, *Acta Physica Austriaca*, 
    180    14 (1961) 185-211 
    181 .. [#] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
     162P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211 
     163 
     164R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
    182165 
    183166Authorship and Verification 
    … …  
    186169* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    187170* **Last Modified by:**  Paul Kienzle **Date:** April 05, 2017 
    188 * **Last Reviewed by:**  Miguel Gonzales and Paul Butler **Date:** May 24, 
    189   2018 - documentation updated 
     171* **Last Reviewed by:**  Richard Heenan **Date:** April 06, 2017 
    190172""" 
    191173 

  • sasmodels/resolution.py

    r0b9c6df r2d81cfe  
    2020MINIMUM_RESOLUTION = 1e-8 
    2121MINIMUM_ABSOLUTE_Q = 0.02  # relative to the minimum q in the data 
    22 PINHOLE_N_SIGMA = 2.5 # From: Barker & Pedersen 1995 JAC 
    2322 
    2423class Resolution(object): 
    … …  
    6665    *q_calc* is the list of points to calculate, or None if this should 
    6766    be estimated from the *q* and *q_width*. 
    68  
    69     *nsigma* is the width of the resolution function.  Should be 2.5. 
    70     See :func:`pinhole_resolution` for details. 
    71     """ 
    72     def __init__(self, q, q_width, q_calc=None, nsigma=PINHOLE_N_SIGMA): 
     67    """ 
     68    def __init__(self, q, q_width, q_calc=None, nsigma=3): 
    7369        #*min_step* is the minimum point spacing to use when computing the 
    7470        #underlying model.  It should be on the order of 
    … …  
    8682 
    8783        # Protect against models which are not defined for very low q.  Limit 
    88         # the smallest q value evaluated to 0.02*min.  Note that negative q 
    89         # values are trimmed even for broad resolution.  Although not possible 
    90         # from the geometry, they may appear since we are using a truncated 
    91         # gaussian to represent resolution rather than a skew distribution. 
     84        # the smallest q value evaluated (in absolute) to 0.02*min 
    9285        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q) 
    93         self.q_calc = self.q_calc[self.q_calc >= cutoff] 
     86        self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff] 
    9487 
    9588        # Build weight matrix from calculated q values 
    9689        self.weight_matrix = pinhole_resolution( 
    97             self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION), 
    98             nsigma=nsigma) 
     90            self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION)) 
     91        self.q_calc = abs(self.q_calc) 
    9992 
    10093    def apply(self, theory): 
    … …  
    108101    *q* points at which the data is measured. 
    109102 
    110     *qx_width* slit width in qx 
    111  
    112     *qy_width* slit height in qy 
     103    *dqx* slit width in qx 
     104 
     105    *dqy* slit height in qy 
    113106 
    114107    *q_calc* is the list of points to calculate, or None if this should 
    … …  
    161154 
    162155 
    163 def pinhole_resolution(q_calc, q, q_width, nsigma=PINHOLE_N_SIGMA): 
    164     r""" 
     156def pinhole_resolution(q_calc, q, q_width): 
     157    """ 
    165158    Compute the convolution matrix *W* for pinhole resolution 1-D data. 
    166159 
    … …  
    169162    *W*, the resolution smearing can be computed using *dot(W,q)*. 
    170163 
    171     Note that resolution is limited to $\pm 2.5 \sigma$.[1]  The true resolution 
    172     function is a broadened triangle, and does not extend over the entire 
    173     range $(-\infty, +\infty)$.  It is important to impose this limitation 
    174     since some models fall so steeply that the weighted value in gaussian 
    175     tails would otherwise dominate the integral. 
    176  
    177164    *q_calc* must be increasing.  *q_width* must be greater than zero. 
    178  
    179     [1] Barker, J. G., and J. S. Pedersen. 1995. Instrumental Smearing Effects 
    180     in Radially Symmetric Small-Angle Neutron Scattering by Numerical and 
    181     Analytical Methods. Journal of Applied Crystallography 28 (2): 105--14. 
    182     https://doi.org/10.1107/S0021889894010095. 
    183165    """ 
    184166    # The current algorithm is a midpoint rectangle rule.  In the test case, 
    … …  
    188170    cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :]) 
    189171    weights = cdf[1:] - cdf[:-1] 
    190     # Limit q range to +/- 2.5 sigma 
    191     qhigh = q + nsigma*q_width 
    192     #qlow = q - nsigma*q_width  # linear limits 
    193     qlow = q*q/qhigh  # log limits 
    194     weights[q_calc[:, None] < qlow[None, :]] = 0. 
    195     weights[q_calc[:, None] > qhigh[None, :]] = 0. 
    196172    weights /= np.sum(weights, axis=0)[None, :] 
    197173    return weights 
    … …  
    518494 
    519495 
    520 def gaussian(q, q0, dq, nsigma=2.5): 
    521     """ 
    522     Return the truncated Gaussian resolution function. 
     496def gaussian(q, q0, dq): 
     497    """ 
     498    Return the Gaussian resolution function. 
    523499 
    524500    *q0* is the center, *dq* is the width and *q* are the points to evaluate. 
    525501    """ 
    526     # Calculate the density of the tails; the resulting gaussian needs to be 
    527     # scaled by this amount in ordere to integrate to 1.0 
    528     two_tail_density = 2 * (1 + erf(-nsigma/sqrt(2)))/2 
    529     return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq)/(1-two_tail_density) 
     502    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq) 
    530503 
    531504 
    … …  
    585558 
    586559 
    587 def romberg_pinhole_1d(q, q_width, form, pars, nsigma=2.5): 
     560def romberg_pinhole_1d(q, q_width, form, pars, nsigma=5): 
    588561    """ 
    589562    Romberg integration for pinhole resolution. 
    … …  
    705678        np.testing.assert_equal(output, self.y) 
    706679 
    707     # TODO: turn pinhole/slit demos into tests 
    708  
    709     @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    710680    def test_pinhole(self): 
    711681        """ 
    … …  
    716686        theory = 12.0-1000.0*resolution.q_calc 
    717687        output = resolution.apply(theory) 
    718         # Note: answer came from output of previous run.  Non-integer 
    719         # values at ends come from the fact that q_calc does not 
    720         # extend beyond q, and so the weights don't balance. 
    721688        answer = [ 
    722             10.47037734, 9.86925860, 
    723             9., 8., 7., 6., 5., 4., 
    724             3.13074140, 2.52962266, 
     689            10.44785079, 9.84991299, 8.98101708, 
     690            7.99906585, 6.99998311, 6.00001689, 
     691            5.00093415, 4.01898292, 3.15008701, 2.55214921, 
    725692            ] 
    726693        np.testing.assert_allclose(output, answer, atol=1e-8) 
    … …  
    765732        self._compare(q, output, answer, 1e-6) 
    766733 
    767     @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    768734    def test_pinhole(self): 
    769735        """ 
    … …  
    780746        self._compare(q, output, answer, 3e-4) 
    781747 
    782     @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    783748    def test_pinhole_romberg(self): 
    784749        """ 
    … …  
    796761        #                     2*np.pi/pars['radius']/200) 
    797762        #tol = 0.001 
    798         ## The default 2.5 sigma and no extra points gets 1% 
     763        ## The default 3 sigma and no extra points gets 1% 
    799764        q_calc = None  # type: np.ndarray 
    800765        tol = 0.01 
    … …  
    11151080 
    11161081    if isinstance(resolution, Slit1D): 
    1117         width, height = resolution.qx_width, resolution.qy_width 
     1082        width, height = resolution.dqx, resolution.dqy 
    11181083        Iq_romb = romberg_slit_1d(resolution.q, width, height, model, pars) 
    11191084    else: 

  • sasmodels/rst2html.py

    r1fbadb2 r2d81cfe  
    5656    # others don't work properly with math_output! 
    5757    if math_output == "mathjax": 
    58         # TODO: this is copied from docs/conf.py; there should be only one 
    59         mathjax_path = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-MML-AM_CHTML" 
    60         settings = {"math_output": math_output + " " + mathjax_path} 
     58        settings = {"math_output": math_output} 
    6159    else: 
    6260        settings = {"math-output": math_output} 

  • sasmodels/sasview_model.py

    rd533590 r05df1de  
    686686        self._intermediate_results = getattr(calculator, 'results', None) 
    687687        calculator.release() 
    688         #self._model.release() 
     688        self._model.release() 
    689689        return result 
    690690 
    … …  
    719719 
    720720    def set_dispersion(self, parameter, dispersion): 
    721         # type: (str, weights.Dispersion) -> None 
     721        # type: (str, weights.Dispersion) -> Dict[str, Any] 
    722722        """ 
    723723        Set the dispersion object for a model parameter 
    … …  
    742742            self.dispersion[parameter] = dispersion.get_pars() 
    743743        else: 
    744             raise ValueError("%r is not a dispersity or orientation parameter" 
    745                              % parameter) 
     744            raise ValueError("%r is not a dispersity or orientation parameter") 
    746745 
    747746    def _dispersion_mesh(self): 
Note: See TracChangeset for help on using the changeset viewer.