Changeset d5cc54c in sasmodels


Ignore:
Timestamp:
Jan 12, 2017 3:57:43 PM (8 years ago)
Author:
GitHub <noreply@…>
Parents:
15ec718 (diff), 7c2935c (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.
git-author:
jhbakker <j.h.bakker@…> (01/12/17 15:57:43)
git-committer:
GitHub <noreply@…> (01/12/17 15:57:43)
Message:

Merge 7c2935c1e6bf1f935aa7dd83519aeb98216c717a into 15ec718cb0314ef130e02c3ad1dee33479b15466

Files:
2 added
4 deleted
15 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/direct_model.py

    r0444c02 rb397165  
    3131from . import resolution2d 
    3232from .details import make_kernel_args, dispersion_mesh 
    33 from sas.sasgui.perspectives.fitting.fitpage import FitPage 
    3433 
    3534try: 
    … …  
    194193        # interpret data 
    195194        if hasattr(data, 'lam'): 
    196         #if not FitPage.no_transform.GetValue(): #if the no_transform radio button is not active DOES NOT WORK! not active before fitting 
    197195            self.data_type = 'sesans' 
    198196        elif hasattr(data, 'qx_data'): 

  • sasmodels/kernel_header.c

    r218cdbc rb7e8b94  
    148148inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; } 
    149149 
    150 #if 0 
     150#if 1 
     151//think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017 
    151152#define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 
    152153    SINCOS(phi*M_PI_180, sn, cn); \ 
    153154    q = sqrt(qx*qx + qy*qy); \ 
    154     cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * cos(theta*M_PI_180));  \ 
     155    cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \ 
    155156    sn = sqrt(1 - cn*cn); \ 
    156157    } while (0) 

  • sasmodels/kernelpy.py

    r777d369 rb397165  
    1515from .generate import F64 
    1616from .kernel import KernelModel, Kernel 
    17 import gc 
    1817 
    1918try: 
    … …  
    8079        """ 
    8180        self.q = None 
    82         gc.collect() 
    8381 
    8482class PyKernel(Kernel): 

  • sasmodels/models/barbell.py

    r0d6e865 rfcb33e4  
    2020.. math:: 
    2121 
    22     I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q)\right> 
     22    I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q,\alpha).sin(\alpha)\right> 
    2323 
    24 where the amplitude $A(q)$ is given as 
     24where the amplitude $A(q,\alpha)$ with the rod axis at angle $\alpha$ to $q$ is given as 
    2525 
    2626.. math:: 
    2727 
    2828    A(q) =&\ \pi r^2L 
    29         \frac{\sin\left(\tfrac12 qL\cos\theta\right)} 
    30              {\tfrac12 qL\cos\theta} 
    31         \frac{2 J_1(qr\sin\theta)}{qr\sin\theta} \\ 
     29        \frac{\sin\left(\tfrac12 qL\cos\alpha\right)} 
     30             {\tfrac12 qL\cos\alpha} 
     31        \frac{2 J_1(qr\sin\alpha)}{qr\sin\alpha} \\ 
    3232        &\ + 4 \pi R^3 \int_{-h/R}^1 dt 
    33         \cos\left[ q\cos\theta 
     33        \cos\left[ q\cos\alpha 
    3434            \left(Rt + h + {\tfrac12} L\right)\right] 
    3535        \times (1-t^2) 
    36         \frac{J_1\left[qR\sin\theta \left(1-t^2\right)^{1/2}\right]} 
    37              {qR\sin\theta \left(1-t^2\right)^{1/2}} 
     36        \frac{J_1\left[qR\sin\alpha \left(1-t^2\right)^{1/2}\right]} 
     37             {qR\sin\alpha \left(1-t^2\right)^{1/2}} 
    3838 
    3939The $\left<\ldots\right>$ brackets denote an average of the structure over 
    40 all orientations. $\left<A^2(q)\right>$ is then the form factor, $P(q)$. 
     40all orientations. $\left<A^2(q,\alpha)\right>$ is then the form factor, $P(q)$. 
    4141The scale factor is equivalent to the volume fraction of cylinders, each of 
    4242volume, $V$. Contrast $\Delta\rho$ is the difference of scattering length 
    … …  
    8585* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    8686* **Last Modified by:** Paul Butler **Date:** March 20, 2016 
    87 * **Last Reviewed by:** Paul Butler **Date:** March 20, 2016 
     87* **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 
    8888""" 
    8989from numpy import inf 

  • sasmodels/models/capped_cylinder.py

    r0d6e865 rfcb33e4  
    2121.. math:: 
    2222 
    23     I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q)\right> 
     23    I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q,\alpha).sin(\alpha)\right> 
    2424 
    25 where the amplitude $A(q)$ is given as 
     25where the amplitude $A(q,\alpha)$ with the rod axis at angle $\alpha$ to $q$ is given as 
    2626 
    2727.. math:: 
    2828 
    2929    A(q) =&\ \pi r^2L 
    30         \frac{\sin\left(\tfrac12 qL\cos\theta\right)} 
    31             {\tfrac12 qL\cos\theta} 
    32         \frac{2 J_1(qr\sin\theta)}{qr\sin\theta} \\ 
     30        \frac{\sin\left(\tfrac12 qL\cos\alpha\right)} 
     31            {\tfrac12 qL\cos\alpha} 
     32        \frac{2 J_1(qr\sin\alpha)}{qr\sin\alpha} \\ 
    3333        &\ + 4 \pi R^3 \int_{-h/R}^1 dt 
    34         \cos\left[ q\cos\theta 
     34        \cos\left[ q\cos\alpha 
    3535            \left(Rt + h + {\tfrac12} L\right)\right] 
    3636        \times (1-t^2) 
    37         \frac{J_1\left[qR\sin\theta \left(1-t^2\right)^{1/2}\right]} 
    38              {qR\sin\theta \left(1-t^2\right)^{1/2}} 
     37        \frac{J_1\left[qR\sin\alpha \left(1-t^2\right)^{1/2}\right]} 
     38             {qR\sin\alpha \left(1-t^2\right)^{1/2}} 
    3939 
    4040The $\left<\ldots\right>$ brackets denote an average of the structure over 
    … …  
    8888* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    8989* **Last Modified by:** Paul Butler **Date:** September 30, 2016 
    90 * **Last Reviewed by:** Richard Heenan **Date:** March 19, 2016 
     90* **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 
     91 
    9192""" 
    9293from numpy import inf 

  • sasmodels/models/core_shell_bicelle.py

    ra23639a rfcb33e4  
    4141 
    4242    I(Q,\alpha) = \frac{\text{scale}}{V_t} \cdot 
    43         F(Q,\alpha)^2 + \text{background} 
     43        F(Q,\alpha)^2.sin(\alpha) + \text{background} 
    4444 
    4545where 
    … …  
    8585* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    8686* **Last Modified by:** Paul Butler **Date:** September 30, 2016 
    87 * **Last Reviewed by:** Richard Heenan **Date:** October 5, 2016 
     87* **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 
    8888""" 
    8989 

  • sasmodels/models/core_shell_cylinder.py

    r755ecc2 rfcb33e4  
    99.. math:: 
    1010 
    11     I(q,\alpha) = \frac{\text{scale}}{V_s} F^2(q) + \text{background} 
     11    I(q,\alpha) = \frac{\text{scale}}{V_s} F^2(q,\alpha).sin(\alpha) + \text{background} 
    1212 
    1313where 
    … …  
    1515.. math:: 
    1616 
    17     F(q) = &\ (\rho_c - \rho_s) V_c 
     17    F(q,\alpha) = &\ (\rho_c - \rho_s) V_c 
    1818           \frac{\sin \left( q \tfrac12 L\cos\alpha \right)} 
    1919                {q \tfrac12 L\cos\alpha} 

  • sasmodels/models/core_shell_ellipsoid.py

    r73e08ae r7c2935c  
    162162 
    163163q = 0.1 
    164 phi = pi/6 
    165 qx = q*cos(phi) 
    166 qy = q*sin(phi) 
    167 # After redefinition of angles find new reasonable values for unit test 
    168 #tests = [ 
    169 #    # Accuracy tests based on content in test/utest_coreshellellipsoidXTmodel.py 
    170 #    [{'radius_equat_core': 200.0, 
    171 #      'x_core': 0.1, 
    172 #      'thick_shell': 50.0, 
    173 #      'x_polar_shell': 0.2, 
    174 #      'sld_core': 2.0, 
    175 #      'sld_shell': 1.0, 
    176 #      'sld_solvent': 6.3, 
    177 #      'background': 0.001, 
    178 #      'scale': 1.0, 
    179 #     }, 1.0, 0.00189402], 
     164# tests had in old coords theta=0, phi=0; new coords theta=90, phi=0 
     165qx = q*cos(pi/6.0) 
     166qy = q*sin(pi/6.0) 
     167# 11Jan2017 RKH sorted tests after redefinition of angles 
     168tests = [ 
     169     # Accuracy tests based on content in test/utest_coreshellellipsoidXTmodel.py 
     170    [{'radius_equat_core': 200.0, 
     171      'x_core': 0.1, 
     172      'thick_shell': 50.0, 
     173      'x_polar_shell': 0.2, 
     174      'sld_core': 2.0, 
     175      'sld_shell': 1.0, 
     176      'sld_solvent': 6.3, 
     177      'background': 0.001, 
     178      'scale': 1.0, 
     179     }, 1.0, 0.00189402], 
    180180 
    181181    # Additional tests with larger range of parameters 
    182 #    [{'background': 0.01}, 0.1, 11.6915], 
    183  
    184 #    [{'radius_equat_core': 20.0, 
    185 #      'x_core': 200.0, 
    186 #      'thick_shell': 54.0, 
    187 #      'x_polar_shell': 3.0, 
    188 #      'sld_core': 20.0, 
    189 #      'sld_shell': 10.0, 
    190 #      'sld_solvent': 6.0, 
    191 #      'background': 0.0, 
    192 #      'scale': 1.0, 
    193 #     }, 0.01, 8688.53], 
    194  
    195 #   [{'background': 0.001}, (0.4, 0.5), 0.00690673], 
    196  
    197 #   [{'radius_equat_core': 20.0, 
    198 #      'x_core': 200.0, 
    199 #      'thick_shell': 54.0, 
    200 #      'x_polar_shell': 3.0, 
    201 #      'sld_core': 20.0, 
    202 #      'sld_shell': 10.0, 
    203 #      'sld_solvent': 6.0, 
    204 #      'background': 0.01, 
    205 #      'scale': 0.01, 
    206 #     }, (qx, qy), 0.0100002], 
    207 #    ] 
     182    [{'background': 0.01}, 0.1, 11.6915], 
     183 
     184    [{'radius_equat_core': 20.0, 
     185      'x_core': 200.0, 
     186      'thick_shell': 54.0, 
     187      'x_polar_shell': 3.0, 
     188      'sld_core': 20.0, 
     189      'sld_shell': 10.0, 
     190      'sld_solvent': 6.0, 
     191      'background': 0.0, 
     192      'scale': 1.0, 
     193     }, 0.01, 8688.53], 
     194 
     195   # 2D tests 
     196   [{'background': 0.001, 
     197     'theta': 90.0, 
     198     'phi': 0.0, 
     199     }, (0.4, 0.5), 0.00690673], 
     200 
     201   [{'radius_equat_core': 20.0, 
     202      'x_core': 200.0, 
     203      'thick_shell': 54.0, 
     204      'x_polar_shell': 3.0, 
     205      'sld_core': 20.0, 
     206      'sld_shell': 10.0, 
     207      'sld_solvent': 6.0, 
     208      'background': 0.01, 
     209      'scale': 0.01, 
     210      'theta': 90.0, 
     211      'phi': 0.0, 
     212     }, (qx, qy), 0.01000025], 
     213    ] 

  • sasmodels/models/cylinder.py

    r4cdd0cc rb7e8b94  
    22# Note: model title and parameter table are inserted automatically 
    33r""" 
    4 The form factor is normalized by the particle volume V = \piR^2L. 
     4 
    55For information about polarised and magnetic scattering, see 
    66the :ref:`magnetism` documentation. 
    … …  
    1414.. math:: 
    1515 
    16     P(q,\alpha) = \frac{\text{scale}}{V} F^2(q,\alpha) + \text{background} 
     16    P(q,\alpha) = \frac{\text{scale}}{V} F^2(q,\alpha).sin(\alpha) + \text{background} 
    1717 
    1818where 
    … …  
    2525           \frac{J_1 \left(q R \sin \alpha\right)}{q R \sin \alpha} 
    2626 
    27 and $\alpha$ is the angle between the axis of the cylinder and $\vec q$, $V$ 
     27and $\alpha$ is the angle between the axis of the cylinder and $\vec q$, $V =\pi R^2L$ 
    2828is the volume of the cylinder, $L$ is the length of the cylinder, $R$ is the 
    2929radius of the cylinder, and $\Delta\rho$ (contrast) is the scattering length 
    … …  
    3535.. math:: 
    3636 
    37     F^2(q)=\int_{0}^{\pi/2}{F^2(q,\theta)\sin(\theta)d\theta} 
     37    F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha}=\int_{0}^{1}{F^2(q,u)du} 
    3838 
    3939 
    40 To provide easy access to the orientation of the cylinder, we define the 
    41 axis of the cylinder using two angles $\theta$ and $\phi$. Those angles 
     40Numerical integration is simplified by a change of variable to $u = cos(\alpha)$ with  
     41$sin(\alpha)=\sqrt{1-u^2}$.  
     42 
     43The output of the 1D scattering intensity function for randomly oriented 
     44cylinders is thus given by 
     45 
     46.. math:: 
     47 
     48    P(q) = \frac{\text{scale}}{V} 
     49        \int_0^{\pi/2} F^2(q,\alpha) \sin \alpha\ d\alpha + \text{background} 
     50 
     51 
     52NB: The 2nd virial coefficient of the cylinder is calculated based on the 
     53radius and length values, and used as the effective radius for $S(q)$ 
     54when $P(q) \cdot S(q)$ is applied. 
     55 
     56For oriented cylinders, we define the direction of the 
     57axis of the cylinder using two angles $\theta$ (note this is not the 
     58same as the scattering angle used in q) and $\phi$. Those angles 
    4259are defined in :numref:`cylinder-angle-definition` . 
    4360 
    … …  
    4865    Definition of the angles for oriented cylinders. 
    4966 
    50  
    51 NB: The 2nd virial coefficient of the cylinder is calculated based on the 
    52 radius and length values, and used as the effective radius for $S(q)$ 
    53 when $P(q) \cdot S(q)$ is applied. 
    54  
    55 The output of the 1D scattering intensity function for randomly oriented 
    56 cylinders is then given by 
    57  
    58 .. math:: 
    59  
    60     P(q) = \frac{\text{scale}}{V} 
    61         \int_0^{\pi/2} F^2(q,\alpha) \sin \alpha\ d\alpha + \text{background} 
    62  
    63 The $\theta$ and $\phi$ parameters are not used for the 1D output. 
     67The $\theta$ and $\phi$ parameters only appear in the model when fitting 2d data. 
    6468 
    6569Validation 
    … …  
    7478 
    7579    P(q) = \int_0^{\pi/2} d\phi 
    76         \int_0^\pi p(\alpha) P_0(q,\alpha) \sin \alpha\ d\alpha 
     80        \int_0^\pi p(\theta) P_0(q,\theta) \sin \theta\ d\theta 
    7781 
    7882 
    7983where $p(\theta,\phi) = 1$ is the probability distribution for the orientation 
    80 and $P_0(q,\alpha)$ is the scattering intensity for the fully oriented 
     84and $P_0(q,\theta)$ is the scattering intensity for the fully oriented 
    8185system, and then comparing to the 1D result. 
    8286 
    … …  
    145149 
    146150qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
    147 # After redefinition of angles, find new tests values  
    148 #tests = [[{}, 0.2, 0.042761386790780453], 
    149 #         [{}, [0.2], [0.042761386790780453]], 
    150 #         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852], 
    151 #         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]], 
    152 #        ] 
     151# After redefinition of angles, find new tests values.  Was 10 10 in old coords 
     152tests = [[{}, 0.2, 0.042761386790780453], 
     153        [{}, [0.2], [0.042761386790780453]], 
     154#  new coords     
     155        [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852], 
     156        [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]], 
     157# old coords   [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852], 
     158#              [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]], 
     159        ] 
    153160del qx, qy  # not necessary to delete, but cleaner 
    154161# ADDED by:  RKH  ON: 18Mar2016 renamed sld's etc 

  • sasmodels/models/elliptical_cylinder.py

    ra807206 rfcb33e4  
    2020 
    2121    I(\vec q)=\frac{1}{V_\text{cyl}}\int{d\psi}\int{d\phi}\int{ 
    22         p(\theta,\phi,\psi)F^2(\vec q,\alpha,\psi)\sin(\theta)d\theta} 
     22        p(\theta,\phi,\psi)F^2(\vec q,\alpha,\psi)\sin(\alpha)d\alpha} 
    2323 
    2424with the functions 
    … …  
    2626.. math:: 
    2727 
    28     F(\vec q,\alpha,\psi) = 2\frac{J_1(a)\sin(b)}{ab} 
     28    F(q,\alpha,\psi) = 2\frac{J_1(a)\sin(b)}{ab} 
    2929 
    3030where 
    … …  
    3232.. math:: 
    3333 
    34     a &= \vec q\sin(\alpha)\left[ 
    35         r^2_\text{major}\sin^2(\psi)+r^2_\text{minor}\cos(\psi) \right]^{1/2} 
     34    a = qr'\sin(\alpha) 
     35     
     36    b = q\frac{L}{2}\cos(\alpha) 
     37     
     38    r'=\frac{r_{minor}}{\sqrt{2}}\sqrt{(1+\nu^{2}) + (1-\nu^{2})cos(\psi)} 
    3639 
    37     b &= \vec q\frac{L}{2}\cos(\alpha) 
    3840 
    39 and the angle $\Psi$ is defined as the orientation of the major axis of the 
     41and the angle $\psi$ is defined as the orientation of the major axis of the 
    4042ellipse with respect to the vector $\vec q$. The angle $\alpha$ is the angle 
    4143between the axis of the cylinder and $\vec q$. 
    … …  
    9597 
    9698L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and 
    97 Neutron Scattering*, Plenum, New York, (1987) 
     99Neutron Scattering*, Plenum, New York, (1987) [see table 3.4] 
     100 
     101Authorship and Verification 
     102---------------------------- 
     103 
     104* **Author:** 
     105* **Last Modified by:**  
     106* **Last Reviewed by:**  Richard Heenan - corrected equation in docs **Date:** December 21, 2016 
     107 
    98108""" 
    99109 

  • sasmodels/models/stacked_disks.py

    r07300ea rb7e8b94  
    148148            sld_layer=0.0, 
    149149            sld_solvent=5.0, 
    150             theta=0, 
     150            theta=90, 
    151151            phi=0) 
    152 #After redefinition to spherical coordinates find new reasonable test values 
    153 #tests = [ 
    154 #    # Accuracy tests based on content in test/utest_extra_models.py. 
    155 #    # Added 2 tests with n_stacked = 5 using SasView 3.1.2 - PDB 
    156 #    [{'thick_core': 10.0, 
    157 #      'thick_layer': 15.0, 
    158 #      'radius': 3000.0, 
    159 #      'n_stacking': 1.0, 
    160 #      'sigma_d': 0.0, 
    161 #      'sld_core': 4.0, 
    162 #      'sld_layer': -0.4, 
    163 #      'solvent_sd': 5.0, 
    164 #      'theta': 0.0, 
    165 #      'phi': 0.0, 
    166 #      'scale': 0.01, 
    167 #      'background': 0.001, 
    168 #     }, 0.001, 5075.12], 
    169  
    170 #    [{'thick_core': 10.0, 
    171 #      'thick_layer': 15.0, 
    172 #      'radius': 3000.0, 
    173 #      'n_stacking': 5.0, 
    174 #      'sigma_d': 0.0, 
    175 #      'sld_core': 4.0, 
    176 #      'sld_layer': -0.4, 
    177 #      'solvent_sd': 5.0, 
    178 #      'theta': 0.0, 
    179 #      'phi': 0.0, 
    180 #      'scale': 0.01, 
    181 #      'background': 0.001, 
    182 #     }, 0.001, 5065.12793824], 
    183  
    184 #    [{'thick_core': 10.0, 
    185 #      'thick_layer': 15.0, 
    186 #      'radius': 3000.0, 
    187 #      'n_stacking': 5.0, 
    188 #      'sigma_d': 0.0, 
    189 #      'sld_core': 4.0, 
    190 #      'sld_layer': -0.4, 
    191 #      'solvent_sd': 5.0, 
    192 #      'theta': 0.0, 
    193 #      'phi': 0.0, 
    194 #      'scale': 0.01, 
    195 #      'background': 0.001, 
    196 #     }, 0.164, 0.0127673597265], 
    197  
    198 #    [{'thick_core': 10.0, 
    199 #      'thick_layer': 15.0, 
    200 #      'radius': 3000.0, 
    201 #      'n_stacking': 1.0, 
    202 #      'sigma_d': 0.0, 
    203 #      'sld_core': 4.0, 
    204 #      'sld_layer': -0.4, 
    205 #      'solvent_sd': 5.0, 
    206 #      'theta': 0.0, 
    207 #      'phi': 0.0, 
    208 #      'scale': 0.01, 
    209 #      'background': 0.001, 
    210 #     }, [0.001, 90.0], [5075.12, 0.001]], 
    211  
    212 #    [{'thick_core': 10.0, 
    213 #      'thick_layer': 15.0, 
    214 #      'radius': 3000.0, 
    215 #      'n_stacking': 1.0, 
    216 #      'sigma_d': 0.0, 
    217 #      'sld_core': 4.0, 
    218 #      'sld_layer': -0.4, 
    219 #      'solvent_sd': 5.0, 
    220 #      'theta': 0.0, 
    221 #      'phi': 0.0, 
    222 #      'scale': 0.01, 
    223 #      'background': 0.001, 
    224 #     }, ([0.4, 0.5]), [0.00105074, 0.00121761]], 
    225  
    226 #    [{'thick_core': 10.0, 
    227 #      'thick_layer': 15.0, 
    228 #      'radius': 3000.0, 
    229 #      'n_stacking': 1.0, 
    230 #      'sigma_d': 0.0, 
    231 #      'sld_core': 4.0, 
    232 #      'sld_layer': -0.4, 
    233 #     'solvent_sd': 5.0, 
    234 #      'theta': 0.0, 
    235 #      'phi': 0.0, 
    236 #      'scale': 0.01, 
    237 #      'background': 0.001, 
    238 #     }, ([1.3, 1.57]), [0.0010039, 0.0010038]], 
    239 #    ] 
    240 # 21Mar2016   RKH notes that unit tests all have n_stacking=1, ought to test other values 
    241  
     152# After redefinition of spherical coordinates - 
     153# tests had in old coords theta=0, phi=0; new coords theta=90, phi=0 
     154# but should not matter here as so far all the tests are 1D not 2D 
     155tests = [ 
     156# Accuracy tests based on content in test/utest_extra_models.py. 
     157# Added 2 tests with n_stacked = 5 using SasView 3.1.2 - PDB; for which alas q=0.001 values seem closer to n_stacked=1 not 5, changed assuming my 4.1 code OK, RKH 
     158    [{'thick_core': 10.0, 
     159      'thick_layer': 15.0, 
     160      'radius': 3000.0, 
     161      'n_stacking': 1.0, 
     162      'sigma_d': 0.0, 
     163      'sld_core': 4.0, 
     164      'sld_layer': -0.4, 
     165      'solvent_sd': 5.0, 
     166      'theta': 90.0, 
     167      'phi': 0.0, 
     168      'scale': 0.01, 
     169      'background': 0.001, 
     170     }, 0.001, 5075.12], 
     171    [{'thick_core': 10.0, 
     172      'thick_layer': 15.0, 
     173      'radius': 3000.0, 
     174      'n_stacking': 5, 
     175      'sigma_d': 0.0, 
     176      'sld_core': 4.0, 
     177      'sld_layer': -0.4, 
     178      'solvent_sd': 5.0, 
     179      'theta': 90.0, 
     180      'phi': 0.0, 
     181      'scale': 0.01, 
     182      'background': 0.001, 
     183#     }, 0.001, 5065.12793824],    n_stacking=1 not 5 ? slight change in value here 11jan2017, check other cpu types 
     184#     }, 0.001, 5075.11570493], 
     185     }, 0.001, 25325.635693], 
     186    [{'thick_core': 10.0, 
     187      'thick_layer': 15.0, 
     188      'radius': 3000.0, 
     189      'n_stacking': 5, 
     190      'sigma_d': 0.0, 
     191      'sld_core': 4.0, 
     192      'sld_layer': -0.4, 
     193      'solvent_sd': 5.0, 
     194      'theta': 90.0, 
     195      'phi': 0.0, 
     196      'scale': 0.01, 
     197      'background': 0.001, 
     198#     }, 0.164, 0.0127673597265],    n_stacking=1 not 5 ?  slight change in value here 11jan2017, check other cpu types 
     199#     }, 0.164, 0.01480812968], 
     200     }, 0.164, 0.0598367986509], 
     201 
     202    [{'thick_core': 10.0, 
     203      'thick_layer': 15.0, 
     204      'radius': 3000.0, 
     205      'n_stacking': 1.0, 
     206      'sigma_d': 0.0, 
     207      'sld_core': 4.0, 
     208      'sld_layer': -0.4, 
     209      'solvent_sd': 5.0, 
     210      'theta': 90.0, 
     211      'phi': 0.0, 
     212      'scale': 0.01, 
     213      'background': 0.001, 
     214# second test here was at q=90, changed it to q=5, note I(q) is then just value of flat background 
     215     }, [0.001, 5.0], [5075.12, 0.001]], 
     216 
     217    [{'thick_core': 10.0, 
     218      'thick_layer': 15.0, 
     219      'radius': 3000.0, 
     220      'n_stacking': 1.0, 
     221      'sigma_d': 0.0, 
     222      'sld_core': 4.0, 
     223      'sld_layer': -0.4, 
     224      'solvent_sd': 5.0, 
     225      'theta': 90.0, 
     226      'phi': 0.0, 
     227      'scale': 0.01, 
     228      'background': 0.001, 
     229     }, ([0.4, 0.5]), [0.00105074, 0.00121761]], 
     230 
     231    [{'thick_core': 10.0, 
     232      'thick_layer': 15.0, 
     233      'radius': 3000.0, 
     234      'n_stacking': 1.0, 
     235      'sigma_d': 0.0, 
     236      'sld_core': 4.0, 
     237      'sld_layer': -0.4, 
     238     'solvent_sd': 5.0, 
     239      'theta': 90.0, 
     240      'phi': 0.0, 
     241      'scale': 0.01, 
     242      'background': 0.001, 
     243     }, ([1.3, 1.57]), [0.0010039, 0.0010038]], 
     244    ] 
     245# 11Jan2017   RKH checking unit test agai, note they are all 1D, no 2D 
     246 

  • sasmodels/resolution.py

    r26b848d rb397165  
    99from numpy import sqrt, log, log10, exp, pi  # type: ignore 
    1010import numpy as np  # type: ignore 
    11  
    12 from sasmodels import sesans 
    13 from sasmodels.sesans import SesansTransform as SesansTransform 
    1411 
    1512__all__ = ["Resolution", "Perfect1D", "Pinhole1D", "Slit1D", 
    … …  
    4643        raise NotImplementedError("Subclass does not define the apply function") 
    4744 
     45 
    4846class Perfect1D(Resolution): 
    4947    """ 
    … …  
    5755    def apply(self, theory): 
    5856        return theory 
     57 
    5958 
    6059class Pinhole1D(Resolution): 

  • sasmodels/resolution2d.py

    r40a87fa r7e94989  
    6464        # just need q_calc and weights 
    6565        self.data = data 
    66         self.index = index 
    67  
    68         self.qx_data = data.qx_data[index] 
    69         self.qy_data = data.qy_data[index] 
    70         self.q_data = data.q_data[index] 
     66        self.index = index if index is not None else slice(None) 
     67 
     68        self.qx_data = data.qx_data[self.index] 
     69        self.qy_data = data.qy_data[self.index] 
     70        self.q_data = data.q_data[self.index] 
    7171 
    7272        dqx = getattr(data, 'dqx_data', None) 
    … …  
    7474        if dqx is not None and dqy is not None: 
    7575            # Here dqx and dqy mean dq_parr and dq_perp 
    76             self.dqx_data = dqx[index] 
    77             self.dqy_data = dqy[index] 
     76            self.dqx_data = dqx[self.index] 
     77            self.dqy_data = dqy[self.index] 
    7878            ## Remove singular points if exists 
    7979            self.dqx_data[self.dqx_data < SIGMA_ZERO] = SIGMA_ZERO 
    … …  
    126126        # The angle (phi) of the original q point 
    127127        q_phi = np.arctan(q_phi).repeat(nbins)\ 
    128             .reshape(nq, nbins).transpose().flatten() 
     128            .reshape([nq, nbins]).transpose().flatten() 
    129129        ## Find Gaussian weight for each dq bins: The weight depends only 
    130130        #  on r-direction (The integration may not need) 

  • sasmodels/sasview_model.py

    r8977226 r64614ad  
    583583                            % type(qdist)) 
    584584 
     585    def get_composition_models(self): 
     586        """ 
     587            Returns usable models that compose this model 
     588        """ 
     589        s_model = None 
     590        p_model = None 
     591        if hasattr(self._model_info, "composition") \ 
     592           and self._model_info.composition is not None: 
     593            p_model = _make_model_from_info(self._model_info.composition[1][0])() 
     594            s_model = _make_model_from_info(self._model_info.composition[1][1])() 
     595        return p_model, s_model 
     596 
    585597    def calculate_Iq(self, qx, qy=None): 
    586598        # type: (Sequence[float], Optional[Sequence[float]]) -> np.ndarray 

  • sasmodels/sesans.py

    r15ec718 rb397165  
    1414import numpy as np  # type: ignore 
    1515from numpy import pi, exp  # type: ignore 
    16 from scipy.special import j0 
    17 #from mpmath import j0 as j0 
     16from scipy.special import jv as besselj 
     17#import direct_model.DataMixin as model 
    1818         
    19 class SesansTransform(object): 
    20     #: Set of spin-echo lengths in the measured data 
    21     SE = None  # type: np.ndarray 
    22     #: Maximum acceptance of scattering vector in the spin echo encoding dimension (for ToF: Q of min(R) and max(lam)) 
    23     zaccept = None # type: float 
    24     #: Maximum size sensitivity; larger radius requires more computation 
    25     Rmax = None  # type: float 
    26     #: q values to calculate when computing transform 
    27     q = None  # type: np.ndarray 
    28  
    29     # transform arrays 
    30     _H = None  # type: np.ndarray 
    31     _H0 = None # type: np.ndarray 
    32  
    33     def set_transform(self, SE, zaccept, Rmax): 
    34         if self.SE is None or len(SE) != len(self.SE) or np.any(SE != self.SE) or zaccept != self.zaccept or Rmax != self.Rmax: 
    35             self.SE, self.zaccept, self.Rmax = SE, zaccept, Rmax 
    36             self._set_q() 
    37             self._set_hankel() 
    38  
    39     def apply(self, Iq): 
    40         G0 = np.dot(self._H0, Iq) 
    41         G = np.dot(self._H.T, Iq) 
    42         P = G - G0 
    43         return P 
    44  
    45     def _set_q(self): 
    46         #q_min = dq = 0.1 * 2*pi / self.Rmax 
    47  
    48         q_max = 2*pi / (self.SE[1]-self.SE[0]) 
    49         q_min = dq = 0.1 *2*pi / (np.size(self.SE) * self.SE[-1]) 
    50  
    51         #q_min = dq = q_max / 100000 
    52         q=np.arange(q_min, q_max, q_min) 
    53         self.q = q 
    54         self.dq = dq 
    55  
    56     def _set_hankel(self): 
    57         #Rmax = #value in text box somewhere in FitPage? 
    58         q = self.q 
    59         dq = self.dq 
    60         SElength = self.SE 
    61  
    62         H0 = dq / (2 * pi) * q 
    63         q=np.array(q,dtype='float32') 
    64         SElength=np.array(SElength,dtype='float32') 
    65  
    66         # Using numpy tile, dtype is conserved 
    67         repq=np.tile(q,(SElength.size,1)) 
    68         repSE=np.tile(SElength,(q.size,1)) 
    69         H = dq / (2 * pi) * j0(repSE*repq.T)*repq.T 
    70  
    71         # Using numpy meshgrid - meshgrid produces float64 from float32 inputs! Problem for 32-bit OS: Memerrors! 
    72         #H0 = dq / (2 * pi) * q 
    73         #repSE, repq = np.meshgrid(SElength, q) 
    74         #repq=np.array(repq,dtype='float32') 
    75         #repSE=np.array(repSE,dtype='float32') 
    76         #H = dq / (2 * pi) * j0(repSE*repq)*repq 
    77  
    78         self._H, self._H0 = H, H0 
    79  
    80 class SESANS1D(SesansTransform): 
    81     def __init__(self, data, _H0, _H, q_calc): 
    82         # x values of the data (Sasview requires these to be named "q") 
    83         self.q = data.x 
    84         self._H0 = _H0 
    85         self._H = _H 
    86         # Pysmear does some checks on the smearer object, these checks require the "data" object... 
    87         self.data=data 
    88         # q values of the SAS model 
    89         self.q_calc = q_calc # These are the MODEL's q values used by the smearer (in this case: the Hankel transform) 
    90     def apply(self, theory): 
    91         return SesansTransform.apply(self,theory) 
     19def make_q(q_max, Rmax): 
     20    r""" 
     21    Return a $q$ vector suitable for SESANS covering from $2\pi/ (10 R_{\max})$ 
     22    to $q_max$. This is the integration range of the Hankel transform; bigger range and  
     23    more points makes a better numerical integration. 
     24    Smaller q_min will increase reliable spin echo length range.  
     25    Rmax is the "radius" of the largest expected object and can be set elsewhere. 
     26    q_max is determined by the acceptance angle of the SESANS instrument. 
     27    """ 
     28    from sas.sascalc.data_util.nxsunit import Converter 
     29 
     30    q_min = dq = 0.1 * 2*pi / Rmax 
     31    return np.arange(q_min, 
     32                     Converter(q_max[1])(q_max[0], 
     33                                         units="1/A"), 
     34                     dq) 
     35     
     36def make_all_q(data): 
     37    """ 
     38    Return a $q$ vector suitable for calculating the total scattering cross section for 
     39    calculating the effect of finite acceptance angles on Time of Flight SESANS instruments. 
     40    If no acceptance is given, or unwanted (set "unwanted" flag in paramfile), no all_q vector is needed. 
     41    If the instrument has a rectangular acceptance, 2 all_q vectors are needed. 
     42    If the instrument has a circular acceptance, 1 all_q vector is needed 
     43     
     44    """ 
     45    if not data.has_no_finite_acceptance: 
     46        return [] 
     47    elif data.has_yz_acceptance(data): 
     48        # compute qx, qy 
     49        Qx, Qy = np.meshgrid(qx, qy) 
     50        return [Qx, Qy] 
     51    else: 
     52        # else only need q 
     53        # data.has_z_acceptance 
     54        return [q] 
     55 
     56def transform(data, q_calc, Iq_calc, qmono, Iq_mono): 
     57    """ 
     58    Decides which transform type is to be used, based on the experiment data file contents (header) 
     59    (2016-03-19: currently controlled from parameters script) 
     60    nqmono is the number of q vectors to be used for the detector integration 
     61    """ 
     62    nqmono = len(qmono) 
     63    if nqmono == 0: 
     64        result = call_hankel(data, q_calc, Iq_calc) 
     65    elif nqmono == 1: 
     66        q = qmono[0] 
     67        result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono) 
     68    else: 
     69        Qx, Qy = [qmono[0], qmono[1]] 
     70        Qx = np.reshape(Qx, nqx, nqy) 
     71        Qy = np.reshape(Qy, nqx, nqy) 
     72        Iq_mono = np.reshape(Iq_mono, nqx, nqy) 
     73        qx = Qx[0, :] 
     74        qy = Qy[:, 0] 
     75        result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono) 
     76 
     77    return result 
     78 
     79def call_hankel(data, q_calc, Iq_calc): 
     80    return hankel((data.x, data.x_unit), 
     81                  (data.lam, data.lam_unit), 
     82                  (data.sample.thickness, 
     83                   data.sample.thickness_unit), 
     84                  q_calc, Iq_calc) 
     85   
     86def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono): 
     87    return hankel(data.x, data.lam * 1e-9, 
     88                  data.sample.thickness / 10, 
     89                  q_calc, Iq_calc) 
     90                   
     91def call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono): 
     92    return hankel(data.x, data.y, data.lam * 1e-9, 
     93                  data.sample.thickness / 10, 
     94                  q_calc, Iq_calc) 
     95                         
     96def TotalScatter(model, parameters):  #Work in progress!! 
     97#    Calls a model with existing model parameters already in place, then integrate the product of q and I(q) from 0 to (4*pi/lambda) 
     98    allq = np.linspace(0,4*pi/wavelength,1000) 
     99    allIq = 1 
     100    integral = allq*allIq 
     101     
     102 
     103 
     104def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still 
     105#============================================================================== 
     106#     2D Cosine Transform if "wavelength" is a vector 
     107#============================================================================== 
     108#allq is the q-space needed to create the total scattering cross-section 
     109 
     110    Gprime = np.zeros_like(wavelength, 'd') 
     111    s = np.zeros_like(wavelength, 'd') 
     112    sd = np.zeros_like(wavelength, 'd') 
     113    Gprime = np.zeros_like(wavelength, 'd') 
     114    f = np.zeros_like(wavelength, 'd') 
     115    for i, wavelength_i in enumerate(wavelength): 
     116        z = magfield*wavelength_i 
     117        allq=np.linspace() #for calculating the Q-range of the  scattering power integral 
     118        allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow 
     119        alldq = (allq[1]-allq[0])*1e10 
     120        sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 
     121        s[i]=1-exp(-sigma) 
     122        for j, Iqy_j, qy_j in enumerate(qy): 
     123            for k, Iqz_k, qz_k in enumerate(qz): 
     124                Iq = np.sqrt(Iqy_j^2+Iqz_k^2) 
     125                q = np.sqrt(qy_j^2 + qz_k^2) 
     126                Gintegral = Iq*cos(z*Qz_k) 
     127                Gprime[i] += Gintegral 
     128#                sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i] 
     129#                s[i] += 1-exp(Totalscatter(modelname)*thickness) 
     130#                For now, work with standard 2-phase scatter 
     131 
     132 
     133                sd[i] += Iq 
     134        f[i] = 1-s[i]+sd[i] 
     135        P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i] 
     136 
     137 
     138 
     139 
     140def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname): 
     141#============================================================================== 
     142#     HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS 
     143#============================================================================== 
     144#acceptq is the q-space needed to create limited acceptance effect 
     145    SElength= wavelength*magfield 
     146    G = np.zeros_like(SElength, 'd') 
     147    threshold=2*pi*theta/wavelength 
     148    for i, SElength_i in enumerate(SElength): 
     149        allq=np.linspace() #for calculating the Q-range of the  scattering power integral 
     150        allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow 
     151        alldq = (allq[1]-allq[0])*1e10 
     152        sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 
     153        s[i]=1-exp(-sigma) 
     154 
     155        dq = (q[1]-q[0])*1e10 
     156        a = (x<threshold) 
     157        acceptq = a*q 
     158        acceptIq = a*Iq 
     159 
     160        G[i] = np.sum(besselj(0, acceptq*SElength_i)*acceptIq*acceptq*dq) 
     161 
     162#        G[i]=np.sum(integral) 
     163 
     164    G *= dq*1e10*2*pi 
     165 
     166    P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0])) 
     167     
     168def hankel(SElength, wavelength, thickness, q, Iq): 
     169    r""" 
     170    Compute the expected SESANS polarization for a given SANS pattern. 
     171 
     172    Uses the hankel transform followed by the exponential.  The values for *zz* 
     173    (or spin echo length, or delta), wavelength and sample thickness should 
     174    come from the dataset.  $q$ should be chosen such that the oscillations 
     175    in $I(q)$ are well sampled (e.g., $5 \cdot 2 \pi/d_{\max}$). 
     176 
     177    *SElength* [A] is the set of $z$ points at which to compute the 
     178    Hankel transform 
     179 
     180    *wavelength* [m]  is the wavelength of each individual point *zz* 
     181 
     182    *thickness* [cm] is the sample thickness. 
     183 
     184    *q* [A$^{-1}$] is the set of $q$ points at which the model has been 
     185    computed. These should be equally spaced. 
     186 
     187    *I* [cm$^{-1}$] is the value of the SANS model at *q* 
     188    """ 
     189 
     190    from sas.sascalc.data_util.nxsunit import Converter 
     191    wavelength = Converter(wavelength[1])(wavelength[0],"A") 
     192    thickness = Converter(thickness[1])(thickness[0],"A") 
     193    Iq = Converter("1/cm")(Iq,"1/A") # All models default to inverse centimeters 
     194    SElength = Converter(SElength[1])(SElength[0],"A") 
     195 
     196    G = np.zeros_like(SElength, 'd') 
     197#============================================================================== 
     198#     Hankel Transform method if "wavelength" is a scalar; mono-chromatic SESANS 
     199#============================================================================== 
     200    for i, SElength_i in enumerate(SElength): 
     201        integral = besselj(0, q*SElength_i)*Iq*q 
     202        G[i] = np.sum(integral) 
     203    G0 = np.sum(Iq*q) 
     204 
     205    # [m^-1] step size in q, needed for integration 
     206    dq = (q[1]-q[0]) 
     207 
     208    # integration step, convert q into [m**-1] and 2 pi circle integration 
     209    G *= dq*2*pi 
     210    G0 = np.sum(Iq*q)*dq*2*np.pi 
     211 
     212    P = exp(thickness*wavelength**2/(4*pi**2)*(G-G0)) 
     213 
     214    return P 
Note: See TracChangeset for help on using the changeset viewer.