Changeset 56b2687 in sasmodels for sasmodels/models


Ignore:
Timestamp:
Jul 14, 2016 1:35:58 PM (8 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
98ba1fc
Parents:
61f8638 (diff), fa800e72 (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 polydisp

Conflicts:

README.rst
sasmodels/core.py
sasmodels/data.py
sasmodels/generate.py
sasmodels/kernelcl.py
sasmodels/kerneldll.py
sasmodels/sasview_model.py

Location:
sasmodels/models
Files:
2 added
23 edited
3 moved

Legend:

Unmodified
Added
Removed
  • sasmodels/models/core_shell_parallelepiped.py

    rec45c4f r500128b  
    66can be different on all three (pairs) of faces.** 
    77 
    8 The form factor is normalized by the particle volume *V* such that 
     8The form factor is normalized by the particle volume $V$ such that 
    99 
    10 *I(q)* = *scale* \* <*f*\ :sup:`2`> / *V* + *background* 
     10.. math:: 
    1111 
    12 where < > is an average over all possible orientations of the rectangular solid. 
     12    I(q) = \text{scale}\frac{\langle f^2 \rangle}{V} + \text{background} 
     13 
     14where $\langle \ldots \rangle$ is an average over all possible orientations 
     15of the rectangular solid. 
    1316 
    1417An instrument resolution smeared version of the model is also provided. 
     
    1922 
    2023The function calculated is the form factor of the rectangular solid below. 
    21 The core of the solid is defined by the dimensions *A*, *B*, *C* such that 
    22 *A* < *B* < *C*. 
     24The core of the solid is defined by the dimensions $A$, $B$, $C$ such that 
     25$A < B < C$. 
    2326 
    2427.. image:: img/core_shell_parallelepiped_geometry.jpg 
    2528 
    26 There are rectangular "slabs" of thickness $t_A$ that add to the *A* dimension 
    27 (on the *BC* faces). There are similar slabs on the *AC* $(=t_B)$ and *AB* 
    28 $(=t_C)$ faces. The projection in the *AB* plane is then 
     29There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension 
     30(on the $BC$ faces). There are similar slabs on the $AC$ $(=t_B)$ and $AB$ 
     31$(=t_C)$ faces. The projection in the $AB$ plane is then 
    2932 
    3033.. image:: img/core_shell_parallelepiped_projection.jpg 
     
    4346 
    4447**For the calculation of the form factor to be valid, the sides of the solid 
    45 MUST be chosen such that** *A* < *B* < *C*. 
     48MUST be chosen such that** $A < B < C$. 
    4649**If this inequality is not satisfied, the model will not report an error, 
    4750and the calculation will not be correct.** 
     
    4952FITTING NOTES 
    5053If the scale is set equal to the particle volume fraction, |phi|, the returned 
    51 value is the scattered intensity per unit volume; ie, *I(q)* = |phi| *P(q)*. 
     54value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 
    5255However, **no interparticle interference effects are included in this calculation.** 
    5356 
     
    5659 
    5760Constraints must be applied during fitting to ensure that the inequality 
    58 *A* < *B* < *C* is not violated. The calculation will not report an error, 
     61$A < B < C$ is not violated. The calculation will not report an error, 
    5962but the results will not be correct. 
    6063 
     
    6467based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    6568and length $(C+2t_C)$ values, and used as the effective radius 
    66 for *S(Q)* when *P(Q)* \* *S(Q)* is applied. 
     69for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
    6770 
    6871.. Comment by Miguel Gonzalez: 
     
    7174 
    7275To provide easy access to the orientation of the parallelepiped, we define the 
    73 axis of the cylinder using three angles |theta|, |phi| and |bigpsi|. 
     76axis of the cylinder using three angles $\theta$, $\phi$ and $\Psi$. 
    7477(see :ref:`cylinder orientation <cylinder-angle-definition>`). 
    75 The angle |bigpsi| is the rotational angle around the *long_c* axis against the 
    76 *q* plane. For example, |bigpsi| = 0 when the *short_b* axis is parallel to the 
     78The angle $\Psi$ is the rotational angle around the *long_c* axis against the 
     79$q$ plane. For example, $\Psi = 0$ when the *short_b* axis is parallel to the 
    7780*x*-axis of the detector. 
    7881 
  • sasmodels/models/hollow_rectangular_prism.py

    rec45c4f r117090a  
    33r""" 
    44 
    5 This model provides the form factor, *P(q)*, for a hollow rectangular 
    6 parallelepiped with a wall of thickness |bigdelta|. 
     5This model provides the form factor, $P(q)$, for a hollow rectangular 
     6parallelepiped with a wall of thickness $\Delta$. 
    77It computes only the 1D scattering, not the 2D. 
    88 
     
    2424  \int_0^{\frac{\pi}{2}} A_{P\Delta}^2(q) \, \sin\theta \, d\theta \, d\phi 
    2525 
    26 where |theta| is the angle between the *z* axis and the longest axis 
    27 of the parallelepiped, |phi| is the angle between the scattering vector 
    28 (lying in the *xy* plane) and the *y* axis, and 
     26where $\theta$ is the angle between the $z$ axis and the longest axis 
     27of the parallelepiped, $\phi$ is the angle between the scattering vector 
     28(lying in the $xy$ plane) and the $y$ axis, and 
    2929 
    3030.. math:: 
     
    4949  \end{align} 
    5050 
    51 where *A*, *B* and *C* are the external sides of the parallelepiped fulfilling 
    52 :math:`A \le B \le C`, and the volume *V* of the parallelepiped is 
     51where $A$, $B$ and $C$ are the external sides of the parallelepiped fulfilling 
     52$A \le B \le C$, and the volume $V$ of the parallelepiped is 
    5353 
    5454.. math:: 
     
    5858 
    5959.. math:: 
    60   I(q) = \mbox{scale} \times V \times (\rho_{\mbox{p}} - 
    61   \rho_{\mbox{solvent}})^2 \times P(q) 
     60  I(q) = \text{scale} \times V \times (\rho_{\text{p}} - 
     61  \rho_{\text{solvent}})^2 \times P(q) + \text{background} 
    6262 
    63 where :math:`\rho_{\mbox{p}}` is the scattering length of the parallelepiped, 
    64 :math:`\rho_{\mbox{solvent}}` is the scattering length of the solvent, 
     63where $\rho_{\text{p}}$ is the scattering length of the parallelepiped, 
     64$\rho_{\text{solvent}}$ is the scattering length of the solvent, 
    6565and (if the data are in absolute units) *scale* represents the volume fraction 
    6666(which is unitless). 
  • sasmodels/models/lamellar.py

    rd2bb604 r56b2687  
    99.. math:: 
    1010 
    11     I(q) = scale*\frac{2\pi P(q)}{q^2\delta } 
     11    I(q) = \text{scale}\frac{2\pi P(q)}{q^2\delta} + \text{background} 
    1212 
    1313 
     
    1616.. math:: 
    1717 
    18    P(q) = \frac{2\Delta\rho^2}{q^2}(1-cos(q\delta)) = \frac{4\Delta\rho^2}{q^2}sin^2(\frac{q\delta}{2}) 
     18   P(q) = \frac{2\Delta\rho^2}{q^2}(1-\cos(q\delta)) 
     19        = \frac{4\Delta\rho^2}{q^2}\sin^2\left(\frac{q\delta}{2}\right) 
    1920 
    2021where $\delta$ is the total layer thickness and $\Delta\rho$ is the scattering length density difference. 
  • sasmodels/models/sc_paracrystal.py

    rec45c4f r500128b  
    1313.. math:: 
    1414 
    15     I(q) = \frac{scale}{V_p}V_{lattice}P(q)Z(q) 
     15    I(q) = \text{scale}\frac{V_\text{lattice}P(q)Z(q)}{V_p} + \text{background} 
    1616 
    1717where scale is the volume fraction of spheres, $V_p$ is the volume of 
    18 the primary particle, $V_{lattice}$ is a volume correction for the crystal 
     18the primary particle, $V_\text{lattice}$ is a volume correction for the crystal 
    1919structure, $P(q)$ is the form factor of the sphere (normalized), and 
    2020$Z(q)$ is the paracrystalline structure factor for a simple cubic structure. 
     
    2828.. math:: 
    2929 
    30     V_{lattice}=\frac{4\pi}{3}\frac{R^3}{D^3} 
     30    V_\text{lattice}=\frac{4\pi}{3}\frac{R^3}{D^3} 
    3131 
    3232The distortion factor (one standard deviation) of the paracrystal is included 
  • sasmodels/models/squarewell.py

    rd2bb604 r56b2687  
    1414The interaction potential is: 
    1515 
    16   .. image:: img\squarewell.png 
     16  .. image:: img/squarewell.png 
    1717 
    1818.. math:: 
  • sasmodels/models/cylinder.c

    r26141cb re9b1663d  
    33double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    44    double radius, double length, double theta, double phi); 
     5 
     6#define INVALID(v) (v.radius<0 || v.length<0) 
    57 
    68double form_volume(double radius, double length) 
     
    1517    double length) 
    1618{ 
    17     // TODO: return NaN if radius<0 or length<0? 
    1819    // precompute qr and qh to save time in the loop 
    1920    const double qr = q*radius; 
     
    4748    double phi) 
    4849{ 
    49     // TODO: return NaN if radius<0 or length<0? 
    5050    double sn, cn; // slots to hold sincos function output 
    5151 
  • sasmodels/models/cylinder.py

    rf247314 r7ae2b7f  
    8282""" 
    8383 
    84 import numpy as np 
    85 from numpy import pi, inf 
     84import numpy as np  # type: ignore 
     85from numpy import pi, inf  # type: ignore 
    8686 
    8787name = "cylinder" 
  • sasmodels/models/flexible_cylinder.c

    re6408d0 r4937980  
    1 double form_volume(double length, double kuhn_length, double radius); 
    2 double Iq(double q, double length, double kuhn_length, double radius, 
    3           double sld, double solvent_sld); 
    4 double Iqxy(double qx, double qy, double length, double kuhn_length, 
    5             double radius, double sld, double solvent_sld); 
    6 double flexible_cylinder_kernel(double q, double length, double kuhn_length, 
    7                                 double radius, double sld, double solvent_sld); 
    8  
    9  
    10 double form_volume(double length, double kuhn_length, double radius) 
     1static double 
     2form_volume(length, kuhn_length, radius) 
    113{ 
    124    return 1.0; 
    135} 
    146 
    15 double flexible_cylinder_kernel(double q, 
    16           double length, 
    17           double kuhn_length, 
    18           double radius, 
    19           double sld, 
    20           double solvent_sld) 
     7static double 
     8Iq(double q, 
     9   double length, 
     10   double kuhn_length, 
     11   double radius, 
     12   double sld, 
     13   double solvent_sld) 
    2114{ 
    22  
    23     const double cont = sld-solvent_sld; 
    24     const double qr = q*radius; 
    25     //const double crossSect = (2.0*J1(qr)/qr)*(2.0*J1(qr)/qr); 
    26     const double crossSect = sas_J1c(qr); 
    27     double flex = Sk_WR(q,length,kuhn_length); 
    28     flex *= crossSect*crossSect; 
    29     flex *= M_PI*radius*radius*length; 
    30     flex *= cont*cont; 
    31     flex *= 1.0e-4; 
    32     return flex; 
     15    const double contrast = sld - solvent_sld; 
     16    const double cross_section = sas_J1c(q*radius); 
     17    const double volume = M_PI*radius*radius*length; 
     18    const double flex = Sk_WR(q, length, kuhn_length); 
     19    return 1.0e-4 * volume * square(contrast*cross_section) * flex; 
    3320} 
    34  
    35 double Iq(double q, 
    36           double length, 
    37           double kuhn_length, 
    38           double radius, 
    39           double sld, 
    40           double solvent_sld) 
    41 { 
    42  
    43     double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld); 
    44     return result; 
    45 } 
    46  
    47 double Iqxy(double qx, double qy, 
    48             double length, 
    49             double kuhn_length, 
    50             double radius, 
    51             double sld, 
    52             double solvent_sld) 
    53 { 
    54     double q; 
    55     q = sqrt(qx*qx+qy*qy); 
    56     double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld); 
    57  
    58     return result; 
    59 } 
  • sasmodels/models/gel_fit.c

    r30b4ddf r03cac08  
    1 double form_volume(void); 
    2  
    3 double Iq(double q, 
    4           double guinier_scale, 
    5           double lorentzian_scale, 
    6           double gyration_radius, 
    7           double fractal_exp, 
    8           double cor_length); 
    9  
    10 double Iqxy(double qx, double qy, 
    11           double guinier_scale, 
    12           double lorentzian_scale, 
    13           double gyration_radius, 
    14           double fractal_exp, 
    15           double cor_length); 
    16  
    17 static double _gel_fit_kernel(double q, 
     1static double Iq(double q, 
    182          double guinier_scale, 
    193          double lorentzian_scale, 
     
    248    // Lorentzian Term 
    259    ////////////////////////double a(x[i]*x[i]*zeta*zeta); 
    26     double lorentzian_term = q*q*cor_length*cor_length; 
     10    double lorentzian_term = square(q*cor_length); 
    2711    lorentzian_term = 1.0 + ((fractal_exp + 1.0)/3.0)*lorentzian_term; 
    2812    lorentzian_term = pow(lorentzian_term, (fractal_exp/2.0) ); 
     
    3014    // Exponential Term 
    3115    ////////////////////////double d(x[i]*x[i]*rg*rg); 
    32     double exp_term = q*q*gyration_radius*gyration_radius; 
     16    double exp_term = square(q*gyration_radius); 
    3317    exp_term = exp(-1.0*(exp_term/3.0)); 
    3418 
     
    3721    return result; 
    3822} 
    39 double form_volume(void){ 
    40     // Unused, so free to return garbage. 
    41     return NAN; 
    42 } 
    43  
    44 double Iq(double q, 
    45           double guinier_scale, 
    46           double lorentzian_scale, 
    47           double gyration_radius, 
    48           double fractal_exp, 
    49           double cor_length) 
    50 { 
    51     return _gel_fit_kernel(q, 
    52                           guinier_scale, 
    53                           lorentzian_scale, 
    54                           gyration_radius, 
    55                           fractal_exp, 
    56                           cor_length); 
    57 } 
    58  
    59 // Iqxy is never called since no orientation or magnetic parameters. 
    60 double Iqxy(double qx, double qy, 
    61           double guinier_scale, 
    62           double lorentzian_scale, 
    63           double gyration_radius, 
    64           double fractal_exp, 
    65           double cor_length) 
    66 { 
    67     double q = sqrt(qx*qx + qy*qy); 
    68     return _gel_fit_kernel(q, 
    69                           guinier_scale, 
    70                           lorentzian_scale, 
    71                           gyration_radius, 
    72                           fractal_exp, 
    73                           cor_length); 
    74 } 
    75  
  • sasmodels/models/hardsphere.py

    rec45c4f rd2bb604  
    149149   """ 
    150150 
    151 Iqxy = """ 
    152     // never called since no orientation or magnetic parameters. 
    153     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    154     """ 
    155  
    156151# ER defaults to 0.0 
    157152# VR defaults to 1.0 
  • sasmodels/models/hayter_msa.py

    rec45c4f rd2bb604  
    8787    return 1.0; 
    8888    """ 
    89 Iqxy = """ 
    90     // never called since no orientation or magnetic parameters. 
    91     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    92     """ 
    9389# ER defaults to 0.0 
    9490# VR defaults to 1.0 
  • sasmodels/models/lamellar_hg.py

    rec45c4f rd2bb604  
    101101    """ 
    102102 
    103 Iqxy = """ 
    104     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    105     """ 
    106  
    107103# ER defaults to 0.0 
    108104# VR defaults to 1.0 
  • sasmodels/models/lamellar_hg_stack_caille.py

    rec45c4f rd2bb604  
    120120    """ 
    121121 
    122 Iqxy = """ 
    123     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    124     """ 
    125  
    126122# ER defaults to 0.0 
    127123# VR defaults to 1.0 
  • sasmodels/models/lamellar_stack_caille.py

    rec45c4f rd2bb604  
    104104    """ 
    105105 
    106 Iqxy = """ 
    107     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    108     """ 
    109  
    110106# ER defaults to 0.0 
    111107# VR defaults to 1.0 
  • sasmodels/models/lamellar_stack_paracrystal.py

    rec45c4f rd2bb604  
    132132    """ 
    133133 
    134 Iqxy = """ 
    135     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    136     """ 
    137  
    138134# ER defaults to 0.0 
    139135# VR defaults to 1.0 
  • sasmodels/models/lib/sas_JN.c

    re6408d0 r4937980  
    4848*/ 
    4949 
    50 static double 
    51 sas_JN( int n, double x ) { 
     50double sas_JN( int n, double x ); 
     51 
     52double sas_JN( int n, double x ) { 
    5253 
    5354    const double MACHEP = 1.11022302462515654042E-16; 
  • sasmodels/models/lib/sph_j1c.c

    re6f1410 rba32cdd  
    77* using double precision that are the source. 
    88*/ 
     9double sph_j1c(double q); 
    910 
    1011// The choice of the number of terms in the series and the cutoff value for 
     
    4344#endif 
    4445 
    45 double sph_j1c(double q); 
    4646double sph_j1c(double q) 
    4747{ 
  • sasmodels/models/lib/sphere_form.c

    rad90df9 rba32cdd  
    1 inline double 
    2 sphere_volume(double radius) 
     1double sphere_volume(double radius); 
     2double sphere_form(double q, double radius, double sld, double solvent_sld); 
     3 
     4double sphere_volume(double radius) 
    35{ 
    46    return M_4PI_3*cube(radius); 
    57} 
    68 
    7 inline double 
    8 sphere_form(double q, double radius, double sld, double solvent_sld) 
     9double sphere_form(double q, double radius, double sld, double solvent_sld) 
    910{ 
    1011    const double fq = sphere_volume(radius) * sph_j1c(q*radius); 
  • sasmodels/models/lib/wrc_cyl.c

    re7678b2 rba32cdd  
    22    Functions for WRC implementation of flexible cylinders 
    33*/ 
     4double Sk_WR(double q, double L, double b); 
     5 
     6 
    47static double 
    58AlphaSquare(double x) 
     
    363366} 
    364367 
    365 double Sk_WR(double q, double L, double b); 
    366368double Sk_WR(double q, double L, double b) 
    367369{ 
  • sasmodels/models/onion.c

    rabdd01c r639c4e3  
    44    double thickness, double A) 
    55{ 
    6   const double vol = 4.0/3.0 * M_PI * r * r * r; 
     6  const double vol = M_4PI_3 * cube(r); 
    77  const double qr = q * r; 
    88  const double alpha = A * r/thickness; 
     
    1919    double sinqr, cosqr; 
    2020    SINCOS(qr, sinqr, cosqr); 
    21     fun = -3.0( 
     21    fun = -3.0*( 
    2222            ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq 
    2323                - (alpha*sinqr/qr - cosqr) 
     
    3232f_linear(double q, double r, double sld, double slope) 
    3333{ 
    34   const double vol = 4.0/3.0 * M_PI * r * r * r; 
     34  const double vol = M_4PI_3 * cube(r); 
    3535  const double qr = q * r; 
    3636  const double bes = sph_j1c(qr); 
     
    5252{ 
    5353  const double bes = sph_j1c(q * r); 
    54   const double vol = 4.0/3.0 * M_PI * r * r * r; 
     54  const double vol = M_4PI_3 * cube(r); 
    5555  return sld * vol * bes; 
    5656} 
     
    6464    r += thickness[i]; 
    6565  } 
    66   return 4.0/3.0 * M_PI * r * r * r; 
     66  return M_4PI_3*cube(r); 
    6767} 
    6868 
    6969static double 
    70 Iq(double q, double core_sld, double core_radius, double solvent_sld, 
    71     double n, double in_sld[], double out_sld[], double thickness[], 
     70Iq(double q, double sld_core, double core_radius, double sld_solvent, 
     71    double n, double sld_in[], double sld_out[], double thickness[], 
    7272    double A[]) 
    7373{ 
    7474  int i; 
    75   r = core_radius; 
    76   f = f_constant(q, r, core_sld); 
     75  double r = core_radius; 
     76  double f = f_constant(q, r, sld_core); 
    7777  for (i=0; i<n; i++){ 
    7878    const double r0 = r; 
     
    9292    } 
    9393  } 
    94   f -= f_constant(q, r, solvent_sld); 
    95   f2 = f * f * 1.0e-4; 
     94  f -= f_constant(q, r, sld_solvent); 
     95  const double f2 = f * f * 1.0e-4; 
    9696 
    9797  return f2; 
  • sasmodels/models/rpa.c

    rabdd01c r639c4e3  
    11double Iq(double q, double case_num, 
    2     double Na, double Phia, double va, double a_sld, double ba, 
    3     double Nb, double Phib, double vb, double b_sld, double bb, 
    4     double Nc, double Phic, double vc, double c_sld, double bc, 
    5     double Nd, double Phid, double vd, double d_sld, double bd, 
     2    double N[], double Phi[], double v[], double L[], double b[], 
    63    double Kab, double Kac, double Kad, 
    74    double Kbc, double Kbd, double Kcd 
    85    ); 
    96 
    10 double Iqxy(double qx, double qy, double case_num, 
    11     double Na, double Phia, double va, double a_sld, double ba, 
    12     double Nb, double Phib, double vb, double b_sld, double bb, 
    13     double Nc, double Phic, double vc, double c_sld, double bc, 
    14     double Nd, double Phid, double vd, double d_sld, double bd, 
    15     double Kab, double Kac, double Kad, 
    16     double Kbc, double Kbd, double Kcd 
    17     ); 
    18  
    19 double form_volume(void); 
    20  
    21 double form_volume(void) 
    22 { 
    23     return 1.0; 
    24 } 
    25  
    267double Iq(double q, double case_num, 
    27     double Na, double Phia, double va, double La, double ba, 
    28     double Nb, double Phib, double vb, double Lb, double bb, 
    29     double Nc, double Phic, double vc, double Lc, double bc, 
    30     double Nd, double Phid, double vd, double Ld, double bd, 
     8    double N[], double Phi[], double v[], double L[], double b[], 
    319    double Kab, double Kac, double Kad, 
    3210    double Kbc, double Kbd, double Kcd 
     
    3614#if 0  // Sasview defaults 
    3715  if (icase <= 1) { 
    38     Na=Nb=1000.0; 
    39     Phia=Phib=0.0000001; 
     16    N[0]=N[1]=1000.0; 
     17    Phi[0]=Phi[1]=0.0000001; 
    4018    Kab=Kac=Kad=Kbc=Kbd=-0.0004; 
    4119    La=Lb=1.0e-12; 
     
    4321    ba=bb=5.0; 
    4422  } else if (icase <= 4) { 
    45     Phia=0.0000001; 
     23    Phi[0]=0.0000001; 
    4624    Kab=Kac=Kad=-0.0004; 
    4725    La=1.0e-12; 
     
    5129#else 
    5230  if (icase <= 1) { 
    53     Na=Nb=0.0; 
    54     Phia=Phib=0.0; 
     31    N[0]=N[1]=0.0; 
     32    Phi[0]=Phi[1]=0.0; 
    5533    Kab=Kac=Kad=Kbc=Kbd=0.0; 
    56     La=Lb=Ld; 
    57     va=vb=vd; 
    58     ba=bb=0.0; 
     34    L[0]=L[1]=L[3]; 
     35    v[0]=v[1]=v[3]; 
     36    b[0]=b[1]=0.0; 
    5937  } else if (icase <= 4) { 
    60     Na = 0.0; 
    61     Phia=0.0; 
     38    N[0] = 0.0; 
     39    Phi[0]=0.0; 
    6240    Kab=Kac=Kad=0.0; 
    63     La=Ld; 
    64     va=vd; 
    65     ba=0.0; 
     41    L[0]=L[3]; 
     42    v[0]=v[3]; 
     43    b[0]=0.0; 
    6644  } 
    6745#endif 
    6846 
    69   const double Xa = q*q*ba*ba*Na/6.0; 
    70   const double Xb = q*q*bb*bb*Nb/6.0; 
    71   const double Xc = q*q*bc*bc*Nc/6.0; 
    72   const double Xd = q*q*bd*bd*Nd/6.0; 
     47  const double Xa = q*q*b[0]*b[0]*N[0]/6.0; 
     48  const double Xb = q*q*b[1]*b[1]*N[1]/6.0; 
     49  const double Xc = q*q*b[2]*b[2]*N[2]/6.0; 
     50  const double Xd = q*q*b[3]*b[3]*N[3]/6.0; 
    7351 
    7452  // limit as Xa goes to 0 is 1 
     
    9876#if 0 
    9977  const double S0aa = icase<5 
    100                       ? 1.0 : Na*Phia*va*Paa; 
     78                      ? 1.0 : N[0]*Phi[0]*v[0]*Paa; 
    10179  const double S0bb = icase<2 
    102                       ? 1.0 : Nb*Phib*vb*Pbb; 
    103   const double S0cc = Nc*Phic*vc*Pcc; 
    104   const double S0dd = Nd*Phid*vd*Pdd; 
     80                      ? 1.0 : N[1]*Phi[1]*v[1]*Pbb; 
     81  const double S0cc = N[2]*Phi[2]*v[2]*Pcc; 
     82  const double S0dd = N[3]*Phi[3]*v[3]*Pdd; 
    10583  const double S0ab = icase<8 
    106                       ? 0.0 : sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb; 
     84                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; 
    10785  const double S0ac = icase<9 
    108                       ? 0.0 : sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb); 
     86                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); 
    10987  const double S0ad = icase<9 
    110                       ? 0.0 : sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc); 
     88                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); 
    11189  const double S0bc = (icase!=4 && icase!=7 && icase!= 9) 
    112                       ? 0.0 : sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc; 
     90                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; 
    11391  const double S0bd = (icase!=4 && icase!=7 && icase!= 9) 
    114                       ? 0.0 : sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc); 
     92                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); 
    11593  const double S0cd = (icase==0 || icase==2 || icase==5) 
    116                       ? 0.0 : sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd; 
     94                      ? 0.0 : sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; 
    11795#else  // sasview equivalent 
    118 //printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,Nc,Phic,vc,Pcc); 
    119   double S0aa = Na*Phia*va*Paa; 
    120   double S0bb = Nb*Phib*vb*Pbb; 
    121   double S0cc = Nc*Phic*vc*Pcc; 
    122   double S0dd = Nd*Phid*vd*Pdd; 
    123   double S0ab = sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb; 
    124   double S0ac = sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb); 
    125   double S0ad = sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc); 
    126   double S0bc = sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc; 
    127   double S0bd = sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc); 
    128   double S0cd = sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd; 
     96//printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N[2],Phi[2],v[2],Pcc); 
     97  double S0aa = N[0]*Phi[0]*v[0]*Paa; 
     98  double S0bb = N[1]*Phi[1]*v[1]*Pbb; 
     99  double S0cc = N[2]*Phi[2]*v[2]*Pcc; 
     100  double S0dd = N[3]*Phi[3]*v[3]*Pdd; 
     101  double S0ab = sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; 
     102  double S0ac = sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); 
     103  double S0ad = sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); 
     104  double S0bc = sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; 
     105  double S0bd = sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); 
     106  double S0cd = sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; 
    129107switch(icase){ 
    130108  case 0: 
     
    311289  // Note: 1e-13 to convert from fm to cm for scattering length 
    312290  const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13; 
    313   const double Lad = icase<5 ? 0.0 : (La/va - Ld/vd)*sqrt_Nav; 
    314   const double Lbd = icase<2 ? 0.0 : (Lb/vb - Ld/vd)*sqrt_Nav; 
    315   const double Lcd = (Lc/vc - Ld/vd)*sqrt_Nav; 
     291  const double Lad = icase<5 ? 0.0 : (L[0]/v[0] - L[3]/v[3])*sqrt_Nav; 
     292  const double Lbd = icase<2 ? 0.0 : (L[1]/v[1] - L[3]/v[3])*sqrt_Nav; 
     293  const double Lcd = (L[2]/v[2] - L[3]/v[3])*sqrt_Nav; 
    316294 
    317295  const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33 
     
    321299 
    322300} 
    323  
    324 double Iqxy(double qx, double qy, 
    325     double case_num, 
    326     double Na, double Phia, double va, double a_sld, double ba, 
    327     double Nb, double Phib, double vb, double b_sld, double bb, 
    328     double Nc, double Phic, double vc, double c_sld, double bc, 
    329     double Nd, double Phid, double vd, double d_sld, double bd, 
    330     double Kab, double Kac, double Kad, 
    331     double Kbc, double Kbd, double Kcd 
    332     ) 
    333 { 
    334     double q = sqrt(qx*qx + qy*qy); 
    335     return Iq(q, 
    336         case_num, 
    337         Na, Phia, va, a_sld, ba, 
    338         Nb, Phib, vb, b_sld, bb, 
    339         Nc, Phic, vc, c_sld, bc, 
    340         Nd, Phid, vd, d_sld, bd, 
    341         Kab, Kac, Kad, 
    342         Kbc, Kbd, Kcd); 
    343 } 
  • sasmodels/models/rpa.py

    rec45c4f ra5b8477  
    8686#   ["name", "units", default, [lower, upper], "type","description"], 
    8787parameters = [ 
    88     ["case_num", CASES, 0, [0, 10], "", "Component organization"], 
     88    ["case_num", "", 1, [CASES], "", "Component organization"], 
    8989 
    90     ["Na", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
    91     ["Phia", "", 0.25, [0, 1], "", "volume fraction"], 
    92     ["va", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
    93     ["La", "fm", 10.0, [-inf, inf], "", "scattering length"], 
    94     ["ba", "Ang", 5.0, [0, inf], "", "segment length"], 
    95  
    96     ["Nb", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
    97     ["Phib", "", 0.25, [0, 1], "", "volume fraction"], 
    98     ["vb", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
    99     ["Lb", "fm", 10.0, [-inf, inf], "", "scattering length"], 
    100     ["bb", "Ang", 5.0, [0, inf], "", "segment length"], 
    101  
    102     ["Nc", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
    103     ["Phic", "", 0.25, [0, 1], "", "volume fraction"], 
    104     ["vc", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
    105     ["Lc", "fm", 10.0, [-inf, inf], "", "scattering length"], 
    106     ["bc", "Ang", 5.0, [0, inf], "", "segment length"], 
    107  
    108     ["Nd", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
    109     ["Phid", "", 0.25, [0, 1], "", "volume fraction"], 
    110     ["vd", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
    111     ["Ld", "fm", 10.0, [-inf, inf], "", "scattering length"], 
    112     ["bd", "Ang", 5.0, [0, inf], "", "segment length"], 
     90    ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
     91    ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"], 
     92    ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
     93    ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"], 
     94    ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"], 
    11395 
    11496    ["Kab", "", -0.0004, [-inf, inf], "", "Interaction parameter"], 
  • sasmodels/models/stickyhardsphere.py

    rec45c4f rd2bb604  
    171171""" 
    172172 
    173 Iqxy = """ 
    174     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    175     """ 
    176  
    177173# ER defaults to 0.0 
    178174# VR defaults to 1.0 
Note: See TracChangeset for help on using the changeset viewer.