Changeset 9a7ef88 in sasmodels


Ignore:
Timestamp:
Apr 17, 2018 8:16:45 AM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
5400ad4
Parents:
17a8c94 (diff), 7c3fb15 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into ticket-896

Files:
10 edited

Legend:

Unmodified
Added
Removed

  • doc/guide/pd/polydispersity.rst

    rf4ae8c4 r29afc50  
    2828sigmas $N_\sigma$ to include from the tails of the distribution, and the 
    2929number of points used to compute the average. The center of the distribution 
    30 is set by the value of the model parameter. 
    31  
    32 Volume parameters have polydispersity *PD* (not to be confused with a 
    33 molecular weight distributions in polymer science), but orientation parameters 
    34 use angular distributions of width $\sigma$. 
     30is set by the value of the model parameter. The meaning of a polydispersity  
     31parameter *PD* (not to be confused with a molecular weight distributions  
     32in polymer science) in a model depends on the type of parameter it is being  
     33applied too. 
     34 
     35The distribution width applied to *volume* (ie, shape-describing) parameters  
     36is relative to the center value such that $\sigma = \mathrm{PD} \cdot \bar x$.  
     37However, the distribution width applied to *orientation* (ie, angle-describing)  
     38parameters is just $\sigma = \mathrm{PD}$. 
    3539 
    3640$N_\sigma$ determines how far into the tails to evaluate the distribution, 
    … …  
    6973or angular orientations, use the Gaussian or Boltzmann distributions. 
    7074 
     75If applying polydispersion to parameters describing angles, use the Uniform  
     76distribution. Beware of using distributions that are always positive (eg, the  
     77Lognormal) because angles can be negative! 
     78 
    7179The array distribution allows a user-defined distribution to be applied. 
    7280 
    … …  
    215223The polydispersity in sasmodels is given by 
    216224 
    217 .. math:: \text{PD} = p = \sigma / x_\text{med} 
    218  
    219 The mean value of the distribution is given by $\bar x = \exp(\mu+ p^2/2)$ 
    220 and the peak value by $\max x = \exp(\mu - p^2)$. 
     225.. math:: \text{PD} = \sigma = p / x_\text{med} 
     226 
     227The mean value of the distribution is given by $\bar x = \exp(\mu+ \sigma^2/2)$ 
     228and the peak value by $\max x = \exp(\mu - \sigma^2)$. 
    221229 
    222230The variance (the square of the standard deviation) of the *lognormal* 
    … …  
    232240.. figure:: pd_lognormal.jpg 
    233241 
    234     Lognormal distribution. 
     242    Lognormal distribution for PD=0.1. 
    235243 
    236244For further information on the Lognormal distribution see: 
    … …  
    334342| 2017-05-08 Paul Kienzle 
    335343| 2018-03-20 Steve King 
     344| 2018-04-04 Steve King 

  • explore/asymint.py

    ra1c32c2 rc462169  
    3030import numpy as np 
    3131import mpmath as mp 
    32 from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10 
     32from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10, arccos 
    3333from numpy.polynomial.legendre import leggauss 
    3434from scipy.integrate import dblquad, simps, romb, romberg 
    … …  
    4141shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1] 
    4242Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2] 
     43 
     44DTYPE = 'd' 
    4345 
    4446class MPenv: 
    … …  
    7274    sas_sinx_x = staticmethod(sp.sas_sinx_x) 
    7375    pi = np.pi 
    74     mpf = staticmethod(float) 
     76    #mpf = staticmethod(float) 
     77    mpf = staticmethod(lambda x: np.array(x, DTYPE)) 
    7578 
    7679SLD = 3 
    … …  
    220223    #A, B, C = 4450, 14000, 47 
    221224    #A, B, C = 445, 140, 47  # integer for the sake of mpf 
    222     A, B, C = 6800, 114, 1380 
    223     DA, DB, DC = 2300, 21, 58 
     225    A, B, C = 114, 1380, 6800 
     226    DA, DB, DC = 21, 58, 2300 
    224227    SLDA, SLDB, SLDC = "5", "-0.3", "11.5" 
     228    ## default parameters from sasmodels 
     229    #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2 
     230    ## swap A-B-C to C-B-A 
     231    #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA 
    225232    #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3 
    226233    #SLD_SOLVENT,CONTRAST = 0, 4 
    … …  
    233240        DA, DC = DC, DA 
    234241        SLDA, SLDC = SLDC, SLDA 
     242    #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
    235243    NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
    236244    NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv) 
    … …  
    348356    return n**2, Iq 
    349357 
     358def gauss_quad_usub(q, n=150, dtype=DTYPE): 
     359    """ 
     360    Compute the integral using gaussian quadrature for n = 20, 76 or 150. 
     361 
     362    Use *u = sin theta* substitution, and restrict integration over a single 
     363    quadrant for shapes that are mirror symmetric about AB, AC and BC planes. 
     364 
     365    Note that this doesn't work for fcc/bcc paracrystals, which instead step 
     366    over the entire 4 pi surface uniformly in theta-phi. 
     367    """ 
     368    z, w = leggauss(n) 
     369    cos_theta = 0.5 * (z + 1) 
     370    theta = arccos(cos_theta) 
     371    phi = pi/2*(0.5 * (z + 1)) 
     372    Atheta, Aphi = np.meshgrid(theta, phi) 
     373    Aw = w[None, :] * w[:, None] 
     374    q, Atheta, Aphi, Aw = [np.asarray(v, dtype=dtype) for v in (q, Atheta, Aphi, Aw)] 
     375    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi) 
     376    Iq = np.sum(Zq*Aw)*0.25 
     377    return n**2, Iq 
     378 
    350379def gridded_2d(q, n=300): 
    351380    """ 
    … …  
    395424    print("gauss-1025", *gauss_quad_2d(Q, n=1025)) 
    396425    print("gauss-2049", *gauss_quad_2d(Q, n=2049)) 
     426    print("gauss-20 usub", *gauss_quad_usub(Q, n=20)) 
     427    print("gauss-76 usub", *gauss_quad_usub(Q, n=76)) 
     428    print("gauss-150 usub", *gauss_quad_usub(Q, n=150)) 
    397429    #gridded_2d(Q, n=2**8+1) 
    398430    gridded_2d(Q, n=2**10+1) 

  • sasmodels/compare.py

    rd86f0fc r5770493  
    752752    """ 
    753753    for k in range(opts['sets']): 
    754         if k > 1: 
     754        if k > 0: 
    755755            # print a separate seed for each dataset for better reproducibility 
    756756            new_seed = np.random.randint(1000000) 
    757             print("Set %d uses -random=%i"%(k+1, new_seed)) 
     757            print("=== Set %d uses -random=%i ==="%(k+1, new_seed)) 
    758758            np.random.seed(new_seed) 
    759759        opts['pars'] = parse_pars(opts, maxdim=maxdim) 
    … …  
    762762        result = run_models(opts, verbose=True) 
    763763        if opts['plot']: 
     764            if opts['is2d'] and k > 0: 
     765                import matplotlib.pyplot as plt 
     766                plt.figure() 
    764767            limits = plot_models(opts, result, limits=limits, setnum=k) 
    765768        if opts['show_weights']: 
    … …  
    13291332 
    13301333    # Evaluate preset parameter expressions 
     1334    # Note: need to replace ':' with '_' in parameter names and expressions 
     1335    # in order to support math on magnetic parameters. 
    13311336    context = MATH.copy() 
    13321337    context['np'] = np 
    1333     context.update(pars) 
     1338    context.update((k.replace(':', '_'), v) for k, v in pars.items()) 
    13341339    context.update((k, v) for k, v in presets.items() if isinstance(v, float)) 
     1340    #for k,v in sorted(context.items()): print(k, v) 
    13351341    for k, v in presets.items(): 
    13361342        if not isinstance(v, float) and not k.endswith('_type'): 
    1337             presets[k] = eval(v, context) 
     1343            presets[k] = eval(v.replace(':', '_'), context) 
    13381344    context.update(presets) 
    1339     context.update((k, v) for k, v in presets2.items() if isinstance(v, float)) 
     1345    context.update((k.replace(':', '_'), v) for k, v in presets2.items() if isinstance(v, float)) 
    13401346    for k, v in presets2.items(): 
    13411347        if not isinstance(v, float) and not k.endswith('_type'): 
    1342             presets2[k] = eval(v, context) 
     1348            presets2[k] = eval(v.replace(':', '_'), context) 
    13431349 
    13441350    # update parameters with presets 

  • sasmodels/details.py

    r108e70e r885753a  
    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 boundaryd 
     245    # Pad value array to a 32 value boundary 
    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) 
    252253    return call_details, data, is_magnetic 
    253254 
    … …  
    296297    mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues] 
    297298    mag = mag.reshape(-1, 3) 
    298     scale = mag[:, 0] 
    299     if np.any(scale): 
     299    if np.any(mag[:, 0] != 0.0): 
     300        M0 = mag[:, 0].copy() 
    300301        theta, phi = radians(mag[:, 1]), radians(mag[:, 2]) 
    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 
     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 
    305305        return True 
    306306    else: 

  • sasmodels/kernel_iq.c

    rd86f0fc rdc6f601  
    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 spins[4]) 
     81static void set_spin_weights(double in_spin, double out_spin, double weight[6]) 
    8282{ 
    8383  in_spin = clip(in_spin, 0.0, 1.0); 
    8484  out_spin = clip(out_spin, 0.0, 1.0); 
    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 
     85  // Note: sasview 3.1 scaled all slds by sqrt(weight) and assumed that 
     86  //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...) 
     87  // which is likely to be the case for simple models. 
     88  weight[0] = sqrt((1.0-in_spin) * (1.0-out_spin)); // dd 
     89  weight[1] = sqrt((1.0-in_spin) * out_spin);       // du.real 
     90  weight[2] = sqrt(in_spin * (1.0-out_spin));       // ud.real 
     91  weight[3] = sqrt(in_spin * out_spin);             // uu 
     92  weight[4] = weight[1]; // du.imag 
     93  weight[5] = weight[2]; // ud.imag 
    9194} 
    9295 
    9396// Compute the magnetic sld 
    9497static double mag_sld( 
    95   const unsigned int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag 
     98  const unsigned int xs, // 0=dd, 1=du.real, 2=ud.real, 3=uu, 4=du.imag, 5=ud.imag 
    9699  const double qx, const double qy, 
    97100  const double px, const double py, 
    … …  
    106109      case 0: // uu => sld - D M_perpx 
    107110          return sld - px*perp; 
    108       case 1: // ud real => -D M_perpy 
     111      case 1: // ud.real => -D M_perpy 
    109112          return py*perp; 
    110       case 2: // du real => -D M_perpy 
     113      case 2: // du.real => -D M_perpy 
    111114          return py*perp; 
    112       case 3: // dd real => sld + D M_perpx 
     115      case 3: // dd => sld + D M_perpx 
    113116          return sld + px*perp; 
    114117    } 
    115118  } else { 
    116119    if (xs== 4) { 
    117       return -mz;  // ud imag => -D M_perpz 
     120      return -mz;  // du.imag => +D M_perpz 
    118121    } else { // index == 5 
    119       return mz;   // du imag => D M_perpz 
     122      return +mz;  // ud.imag => -D M_perpz 
    120123    } 
    121124  } 
    … …  
    312315  //     up_angle = values[NUM_PARS+4]; 
    313316  // TODO: could precompute more magnetism parameters before calling the kernel. 
    314   double spins[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill 
     317  double xs_weights[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill 
    315318  double cos_mspin, sin_mspin; 
    316   set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], spins); 
     319  set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], xs_weights); 
    317320  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin); 
    318321#endif // MAGNETIC 
    … …  
    661664            // loop over uu, ud real, du real, dd, ud imag, du imag 
    662665            for (unsigned int xs=0; xs<6; xs++) { 
    663               const double xs_weight = spins[xs]; 
     666              const double xs_weight = xs_weights[xs]; 
    664667              if (xs_weight > 1.e-8) { 
    665668                // Since the cross section weight is significant, set the slds 
    … …  
    674677                  local_values.vector[sld_index] = 
    675678                    mag_sld(xs, qx, qy, px, py, values[sld_index+2], mx, my, mz); 
     679//if (q_index==0) printf("%d: (qx,qy)=(%g,%g) xs=%d sld%d=%g p=(%g,%g) m=(%g,%g,%g)\n", 
     680//  q_index, qx, qy, xs, sk, local_values.vector[sld_index], px, py, mx, my, mz); 
    676681                } 
    677682                scattering += xs_weight * CALL_KERNEL(); 

  • sasmodels/sasview_model.py

    r3221de0 r7c3fb15  
    593593            # Check whether we have a list of ndarrays [qx,qy] 
    594594            qx, qy = qdist 
    595             if not self._model_info.parameters.has_2d: 
    596                 return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2)) 
    597             else: 
    598                 return self.calculate_Iq(qx, qy) 
     595            return self.calculate_Iq(qx, qy) 
    599596 
    600597        elif isinstance(qdist, np.ndarray): 
    … …  
    677674        call_details, values, is_magnetic = make_kernel_args(calculator, pairs) 
    678675        #call_details.show() 
    679         #print("pairs", pairs) 
     676        #print("================ parameters ==================") 
     677        #for p, v in zip(parameters.call_parameters, pairs): print(p.name, v[0]) 
    680678        #for k, p in enumerate(self._model_info.parameters.call_parameters): 
    681679        #    print(k, p.name, *pairs[k]) 
    … …  
    721719 
    722720    def set_dispersion(self, parameter, dispersion): 
    723         # type: (str, weights.Dispersion) -> Dict[str, Any] 
     721        # type: (str, weights.Dispersion) -> None 
    724722        """ 
    725723        Set the dispersion object for a model parameter 
    … …  
    744742            self.dispersion[parameter] = dispersion.get_pars() 
    745743        else: 
    746             raise ValueError("%r is not a dispersity or orientation parameter") 
     744            raise ValueError("%r is not a dispersity or orientation parameter" 
     745                             % parameter) 
    747746 
    748747    def _dispersion_mesh(self): 
    … …  
    871870    CylinderModel().evalDistribution([0.1, 0.1]) 
    872871 
     872def magnetic_demo(): 
     873    Model = _make_standard_model('sphere') 
     874    model = Model() 
     875    model.setParam('M0:sld', 8) 
     876    q = np.linspace(-0.35, 0.35, 500) 
     877    qx, qy = np.meshgrid(q, q) 
     878    result = model.calculate_Iq(qx.flatten(), qy.flatten()) 
     879    result = result.reshape(qx.shape) 
     880 
     881    import pylab 
     882    pylab.imshow(np.log(result + 0.001)) 
     883    pylab.show() 
     884 
    873885if __name__ == "__main__": 
    874886    print("cylinder(0.1,0.1)=%g"%test_cylinder()) 
     887    #magnetic_demo() 
    875888    #test_product() 
    876889    #test_structure_factor() 

  • sasmodels/models/core_shell_parallelepiped.c

    re077231 rdbf1a60  
    5959 
    6060    // outer integral (with gauss points), integration limits = 0, 1 
     61    // substitute d_cos_alpha for sin_alpha d_alpha 
    6162    double outer_sum = 0; //initialize integral 
    6263    for( int i=0; i<GAUSS_N; i++) { 
    6364        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
    6465        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 
    65  
    66         // inner integral (with gauss points), integration limits = 0, pi/2 
    6766        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 
    6867        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) 
    6971        double inner_sum = 0.0; 
    7072        for(int j=0; j<GAUSS_N; j++) { 
    71             const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
     73            const double u = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
    7274            double sin_beta, cos_beta; 
    73             SINCOS(M_PI_2*beta, sin_beta, cos_beta); 
     75            SINCOS(M_PI_2*u, sin_beta, cos_beta); 
    7476            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 
    7577            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); 
    … …  
    9193            inner_sum += GAUSS_W[j] * f * f; 
    9294        } 
     95        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    9396        inner_sum *= 0.5; 
    9497        // now sum up the outer integral 
    9598        outer_sum += GAUSS_W[i] * inner_sum; 
    9699    } 
     100    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    97101    outer_sum *= 0.5; 
    98102 

  • sasmodels/models/core_shell_parallelepiped.py

    r97be877 r5bc6d21  
    1111.. math:: 
    1212 
    13     I(q) = \text{scale}\frac{\langle f^2 \rangle}{V} + \text{background} 
     13    I(q) = \frac{\text{scale}}{V} \langle P(q,\alpha,\beta) \rangle  
     14    + \text{background} 
    1415 
    1516where $\langle \ldots \rangle$ is an average over all possible orientations 
    16 of the rectangular solid. 
    17  
    18 The function calculated is the form factor of the rectangular solid below. 
     17of the rectangular solid, and the usual $\Delta \rho^2 \ V^2$ term cannot be 
     18pulled out of the form factor term due to the multiple slds in the model. 
     19 
    1920The core of the solid is defined by the dimensions $A$, $B$, $C$ such that 
    2021$A < B < C$. 
    2122 
    22 .. image:: img/core_shell_parallelepiped_geometry.jpg 
     23.. figure:: img/parallelepiped_geometry.jpg 
     24 
     25   Core of the core shell parallelepiped with the corresponding definition 
     26   of sides. 
     27 
    2328 
    2429There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension 
    2530(on the $BC$ faces). There are similar slabs on the $AC$ $(=t_B)$ and $AB$ 
    26 $(=t_C)$ faces. The projection in the $AB$ plane is then 
    27  
    28 .. image:: img/core_shell_parallelepiped_projection.jpg 
    29  
    30 The volume of the solid is 
     31$(=t_C)$ faces. The projection in the $AB$ plane is 
     32 
     33.. figure:: img/core_shell_parallelepiped_projection.jpg 
     34 
     35   AB cut through the core-shell parallelipiped showing the cross secion of 
     36   four of the six shell slabs. As can be seen, this model leaves **"gaps"** 
     37   at the corners of the solid. 
     38 
     39 
     40The total volume of the solid is thus given as 
    3141 
    3242.. math:: 
    3343 
    3444    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    35  
    36 **meaning that there are "gaps" at the corners of the solid.** 
    3745 
    3846The intensity calculated follows the :ref:`parallelepiped` model, with the 
    3947core-shell intensity being calculated as the square of the sum of the 
    40 amplitudes of the core and the slabs on the edges. 
    41  
    42 the scattering amplitude is computed for a particular orientation of the 
    43 core-shell parallelepiped with respect to the scattering vector and then 
    44 averaged over all possible orientations, where $\alpha$ is the angle between 
    45 the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is 
    46 the angle between projection of the particle in the $xy$ detector plane 
    47 and the $y$ axis. 
    48  
    49 .. math:: 
    50  
    51     F(Q) 
     48amplitudes of the core and the slabs on the edges. The scattering amplitude is 
     49computed for a particular orientation of the core-shell parallelepiped with 
     50respect to the scattering vector and then averaged over all possible 
     51orientations, where $\alpha$ is the angle between the $z$ axis and the $C$ axis 
     52of the parallelepiped, and $\beta$ is the angle between the projection of the 
     53particle in the $xy$ detector plane and the $y$ axis. 
     54 
     55.. math:: 
     56 
     57    P(q)=\frac {\int_{0}^{\pi/2}\int_{0}^{\pi/2}F^2(q,\alpha,\beta) \ sin\alpha 
     58    \ d\alpha \ d\beta} {\int_{0}^{\pi/2} \ sin\alpha \ d\alpha \ d\beta} 
     59 
     60and 
     61 
     62.. math:: 
     63 
     64    F(q,\alpha,\beta) 
    5265    &= (\rho_\text{core}-\rho_\text{solvent}) 
    5366       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 
    5467    &+ (\rho_\text{A}-\rho_\text{solvent}) 
    55         \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\ 
     68        \left[S(Q_A, A+2t_A) - S(Q_A, A)\right] S(Q_B, B) S(Q_C, C) \\ 
    5669    &+ (\rho_\text{B}-\rho_\text{solvent}) 
    5770        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ 
    … …  
    6376.. math:: 
    6477 
    65     S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} 
     78    S(Q_X, L) = L \frac{\sin (\tfrac{1}{2} Q_X L)}{\tfrac{1}{2} Q_X L} 
    6679 
    6780and 
    … …  
    6982.. math:: 
    7083 
    71     Q_A &= \sin\alpha \sin\beta \\ 
    72     Q_B &= \sin\alpha \cos\beta \\ 
    73     Q_C &= \cos\alpha 
     84    Q_A &= q \sin\alpha \sin\beta \\ 
     85    Q_B &= q \sin\alpha \cos\beta \\ 
     86    Q_C &= q \cos\alpha 
    7487 
    7588 
    7689where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ 
    77 are the scattering length of the parallelepiped core, and the rectangular 
     90are the scattering lengths of the parallelepiped core, and the rectangular 
    7891slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 
    7992is the scattering length of the solvent. 
     93 
     94.. note::  
     95 
     96   the code actually implements two substitutions: $d(cos\alpha)$ is 
     97   substituted for -$sin\alpha \ d\alpha$ (note that in the 
     98   :ref:`parallelepiped` code this is explicitly implemented with 
     99   $\sigma = cos\alpha$), and $\beta$ is set to $\beta = u \pi/2$ so that 
     100   $du = \pi/2 \ d\beta$.  Thus both integrals go from 0 to 1 rather than 0 
     101   to $\pi/2$. 
    80102 
    81103FITTING NOTES 
    … …  
    94116based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    95117and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 
    96 to give an oblate or prolate particle, to give an effective radius, 
    97 for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     118to give an oblate or prolate particle, to give an effective radius 
     119for $S(q)$ when $P(q) * S(q)$ is applied. 
    98120 
    99121For 2d data the orientation of the particle is required, described using 
    100 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further 
     122angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below. For further 
    101123details of the calculation and angular dispersions see :ref:`orientation`. 
    102124The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, 
    103125$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 
    104126 
    105 For 2d, constraints must be applied during fitting to ensure that the 
    106 inequality $A < B < C$ is not violated, and hence the correct definition 
    107 of angles is preserved. The calculation will not report an error, 
    108 but the results may be not correct. 
     127.. note:: For 2d, constraints must be applied during fitting to ensure that the 
     128   inequality $A < B < C$ is not violated, and hence the correct definition 
     129   of angles is preserved. The calculation will not report an error, 
     130   but the results may be not correct. 
    109131 
    110132.. figure:: img/parallelepiped_angle_definition.png 
    … …  
    113135    Note that rotation $\theta$, initially in the $xz$ plane, is carried 
    114136    out first, then rotation $\phi$ about the $z$ axis, finally rotation 
    115     $\Psi$ is now around the axis of the cylinder. The neutron or X-ray 
     137    $\Psi$ is now around the axis of the particle. The neutron or X-ray 
    116138    beam is along the $z$ axis. 
    117139 
    … …  
    120142    Examples of the angles for oriented core-shell parallelepipeds against the 
    121143    detector plane. 
     144 
     145 
     146Validation 
     147---------- 
     148 
     149Cross-checked against hollow rectangular prism and rectangular prism for equal 
     150thickness overlapping sides, and by Monte Carlo sampling of points within the 
     151shape for non-uniform, non-overlapping sides. 
     152 
    122153 
    123154References 
    … …  
    135166 
    136167* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    137 * **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016 
     168* **Converted to sasmodels by:** Miguel Gonzalez **Date:** February 26, 2016 
    138169* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017 
    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. 
    142170""" 
    143171 

  • sasmodels/models/parallelepiped.c

    r108e70e rdbf1a60  
    3838            inner_total += GAUSS_W[j] * square(si1 * si2); 
    3939        } 
     40        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    4041        inner_total *= 0.5; 
    4142 
    … …  
    4344        outer_total += GAUSS_W[i] * inner_total * si * si; 
    4445    } 
     46    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    4547    outer_total *= 0.5; 
    4648 

  • sasmodels/models/parallelepiped.py

    ref07e95 r5bc6d21  
    22# Note: model title and parameter table are inserted automatically 
    33r""" 
    4 The form factor is normalized by the particle volume. 
    5 For information about polarised and magnetic scattering, see 
    6 the :ref:`magnetism` documentation. 
    7  
    84Definition 
    95---------- 
    106 
    117 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`. 
     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 
     11the :ref:`magnetism` documentation. 
    1412 
    1513.. _parallelepiped-image: 
    … …  
    2624error, or fixing of some dimensions at expected values, may help. 
    2725 
    28 The 1D scattering intensity $I(q)$ is calculated as: 
     26The form factor is normalized by the particle volume and the 1D scattering 
     27intensity $I(q)$ is then calculated as: 
    2928 
    3029.. Comment by Miguel Gonzalez: 
    … …  
    3938 
    4039    I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2 
    41            \left< P(q, \alpha) \right> + \text{background} 
     40           \left< P(q, \alpha, \beta) \right> + \text{background} 
    4241 
    4342where the volume $V = A B C$, the contrast is defined as 
    44 $\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, 
    45 $P(q, \alpha)$ is the form factor corresponding to a parallelepiped oriented 
    46 at an angle $\alpha$ (angle between the long axis C and $\vec q$), 
    47 and the averaging $\left<\ldots\right>$ is applied over all orientations. 
     43$\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, $P(q, \alpha, \beta)$ 
     44is the form factor corresponding to a parallelepiped oriented 
     45at an angle $\alpha$ (angle between the long axis C and $\vec q$), and $\beta$ 
     46( the angle between the projection of the particle in the $xy$ detector plane 
     47and the $y$ axis) and the averaging $\left<\ldots\right>$ is applied over all 
     48orientations. 
    4849 
    4950Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the 
    50 form factor is given by (Mittelbach and Porod, 1961) 
     51form factor is given by (Mittelbach and Porod, 1961 [#Mittelbach]_) 
    5152 
    5253.. math:: 
    … …  
    6667    \mu &= qB 
    6768 
    68 The scattering intensity per unit volume is returned in units of |cm^-1|. 
     69where substitution of $\sigma = cos\alpha$ and $\beta = \pi/2 \ u$ have been 
     70applied. 
    6971 
    7072NB: The 2nd virial coefficient of the parallelepiped is calculated based on 
    … …  
    120122.. math:: 
    121123 
    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 
     124    P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1} 
     125                   {2}qA\cos\alpha)}\right]^2 
     126                  \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1} 
     127                   {2}qB\cos\beta)}\right]^2 
     128                  \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1} 
     129                   {2}qC\cos\gamma)}\right]^2 
    125130 
    126131with 
    … …  
    160165---------- 
    161166 
    162 P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211 
    163  
    164 R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
     167.. [#Mittelbach] P Mittelbach and G Porod, *Acta Physica Austriaca*, 
     168   14 (1961) 185-211 
     169.. [#] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
    165170 
    166171Authorship and Verification 
Note: See TracChangeset for help on using the changeset viewer.