Changeset da63656 in sasmodels for sasmodels/models


Ignore:
Timestamp:
May 4, 2016 11:37:42 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:
24d5b30
Parents:
13b99fd (diff), 47e498b (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:

sasmodels/generate.py
sasmodels/kerneldll.py

Location:
sasmodels/models
Files:
2 added
23 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/lib/sas_gamma.c

    ra062d18 r6ce9048  
    77*/ 
    88 
     9#if defined(NEED_TGAMMA) 
     10static double cephes_stirf(double x) 
     11{ 
     12        const double MAXSTIR=143.01608; 
     13        const double SQTPI=2.50662827463100050242E0; 
     14        double y, w, v; 
     15 
     16        w = 1.0 / x; 
     17 
     18        w = (((( 
     19                7.87311395793093628397E-4*w 
     20                - 2.29549961613378126380E-4)*w 
     21                - 2.68132617805781232825E-3)*w 
     22                + 3.47222221605458667310E-3)*w 
     23                + 8.33333333333482257126E-2)*w 
     24                + 1.0; 
     25        y = exp(x); 
     26        if (x > MAXSTIR) 
     27        { /* Avoid overflow in pow() */ 
     28                v = pow(x, 0.5 * x - 0.25); 
     29                y = v * (v / y); 
     30        } 
     31        else 
     32        { 
     33                y = pow(x, x - 0.5) / y; 
     34        } 
     35        y = SQTPI * y * w; 
     36        return(y); 
     37} 
     38 
     39static double tgamma(double x) { 
     40        double p, q, z; 
     41        int sgngam; 
     42        int i; 
     43 
     44        sgngam = 1; 
     45        if (isnan(x)) 
     46                return(x); 
     47        q = fabs(x); 
     48 
     49        if (q > 33.0) 
     50        { 
     51                if (x < 0.0) 
     52                { 
     53                        p = floor(q); 
     54                        if (p == q) 
     55                        { 
     56                                return (NAN); 
     57                        } 
     58                        i = p; 
     59                        if ((i & 1) == 0) 
     60                                sgngam = -1; 
     61                        z = q - p; 
     62                        if (z > 0.5) 
     63                        { 
     64                                p += 1.0; 
     65                                z = q - p; 
     66                        } 
     67                        z = q * sin(M_PI * z); 
     68                        if (z == 0.0) 
     69                        { 
     70                                return(NAN); 
     71                        } 
     72                        z = fabs(z); 
     73                        z = M_PI / (z * cephes_stirf(q)); 
     74                } 
     75                else 
     76                { 
     77                        z = cephes_stirf(x); 
     78                } 
     79                return(sgngam * z); 
     80        } 
     81 
     82        z = 1.0; 
     83        while (x >= 3.0) 
     84        { 
     85                x -= 1.0; 
     86                z *= x; 
     87        } 
     88 
     89        while (x < 0.0) 
     90        { 
     91                if (x > -1.E-9) 
     92                        goto small; 
     93                z /= x; 
     94                x += 1.0; 
     95        } 
     96 
     97        while (x < 2.0) 
     98        { 
     99                if (x < 1.e-9) 
     100                        goto small; 
     101                z /= x; 
     102                x += 1.0; 
     103        } 
     104 
     105        if (x == 2.0) 
     106                return(z); 
     107 
     108        x -= 2.0; 
     109        p = ((((( 
     110                +1.60119522476751861407E-4*x 
     111                + 1.19135147006586384913E-3)*x 
     112                + 1.04213797561761569935E-2)*x 
     113                + 4.76367800457137231464E-2)*x 
     114                + 2.07448227648435975150E-1)*x 
     115                + 4.94214826801497100753E-1)*x 
     116                + 9.99999999999999996796E-1; 
     117        q = (((((( 
     118                -2.31581873324120129819E-5*x 
     119                + 5.39605580493303397842E-4)*x 
     120                - 4.45641913851797240494E-3)*x 
     121                + 1.18139785222060435552E-2)*x 
     122                + 3.58236398605498653373E-2)*x 
     123                - 2.34591795718243348568E-1)*x 
     124                + 7.14304917030273074085E-2)*x 
     125                + 1.00000000000000000320E0; 
     126        return(z * p / q); 
     127 
     128small: 
     129        if (x == 0.0) 
     130        { 
     131                return (NAN); 
     132        } 
     133        else 
     134                return(z / ((1.0 + 0.5772156649015329 * x) * x)); 
     135} 
     136#endif 
     137 
     138 
    9139inline double sas_gamma( double x) { return tgamma(x+1)/x; } 

  • 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.py

    rec45c4f rd2bb604  
    8282    """ 
    8383 
    84 Iqxy = """ 
    85     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    86     """ 
    87  
    8884# ER defaults to 0.0 
    8985# 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/onion.py

    r416609b rfa5fd8d  
    293293 
    294294#             ["name", "units", default, [lower, upper], "type","description"], 
    295 parameters = [["core_sld", "1e-6/Ang^2", 1.0, [-inf, inf], "", 
     295parameters = [["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "", 
    296296               "Core scattering length density"], 
    297297              ["core_radius", "Ang", 200., [0, inf], "volume", 
    298298               "Radius of the core"], 
    299               ["solvent_sld", "1e-6/Ang^2", 6.4, [-inf, inf], "", 
     299              ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "", 
    300300               "Solvent scattering length density"], 
    301               ["n", "", 1, [0, 10], "volume", 
     301              ["n_shells", "", 1, [0, 10], "volume", 
    302302               "number of shells"], 
    303               ["in_sld[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "", 
     303              ["sld_in[n_shells]", "1e-6/Ang^2", 1.7, [-inf, inf], "", 
    304304               "scattering length density at the inner radius of shell k"], 
    305               ["out_sld[n]", "1e-6/Ang^2", 2.0, [-inf, inf], "", 
     305              ["sld_out[n_shells]", "1e-6/Ang^2", 2.0, [-inf, inf], "", 
    306306               "scattering length density at the outer radius of shell k"], 
    307               ["thickness[n]", "Ang", 40., [0, inf], "volume", 
     307              ["thickness[n_shells]", "Ang", 40., [0, inf], "volume", 
    308308               "Thickness of shell k"], 
    309               ["A[n]", "", 1.0, [-inf, inf], "", 
     309              ["A[n_shells]", "", 1.0, [-inf, inf], "", 
    310310               "Decay rate of shell k"], 
    311311              ] 
    312312 
    313 #source = ["lib/sph_j1c.c", "onion.c"] 
    314  
    315 def Iq(q, *args, **kw): 
    316     return q 
    317  
    318 def Iqxy(qx, *args, **kw): 
    319     return qx 
    320  
    321  
    322 def shape(core_sld, core_radius, solvent_sld, n, in_sld, out_sld, thickness, A): 
     313source = ["lib/sph_j1c.c", "onion.c"] 
     314 
     315#def Iq(q, *args, **kw): 
     316#    return q 
     317 
     318profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 
     319def profile(core_sld, core_radius, solvent_sld, n_shells, 
     320            in_sld, out_sld, thickness, A): 
    323321    """ 
    324     SLD profile 
     322    Returns shape profile with x=radius, y=SLD. 
    325323    """ 
    326324 
    327     total_radius = 1.25*(sum(thickness[:n]) + core_radius + 1) 
     325    total_radius = 1.25*(sum(thickness[:n_shells]) + core_radius + 1) 
    328326    dr = total_radius/400  # 400 points for a smooth plot 
    329327 
    … …  
    338336 
    339337    # add in the shells 
    340     for k in range(n): 
     338    for k in range(n_shells): 
    341339        # Left side of each shells 
    342340        r0 = r[-1] 
    … …  
    365363    beta.append(solvent_sld) 
    366364 
    367     return np.asarray(r), np.asarray(beta) 
     365    return np.asarray(r), np.asarray(beta)*1e-6 
    368366 
    369367def ER(core_radius, n, thickness): 
    … …  
    374372 
    375373demo = { 
    376     "solvent_sld": 2.2, 
    377     "core_sld": 1.0, 
     374    "sld_solvent": 2.2, 
     375    "sld_core": 1.0, 
    378376    "core_radius": 100, 
    379377    "n": 4, 
    380     "in_sld": [0.5, 1.5, 0.9, 2.0], 
    381     "out_sld": [nan, 0.9, 1.2, 1.6], 
     378    "sld_in": [0.5, 1.5, 0.9, 2.0], 
     379    "sld_out": [nan, 0.9, 1.2, 1.6], 
    382380    "thickness": [50, 75, 150, 75], 
    383381    "A": [0, -1, 1e-4, 1], 

  • 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/spherical_sld.py

    rec45c4f rdd71228  
    11r""" 
    2 This model calculates an empirical functional form for SAS data using SpericalSLD profile 
    3  
    4 Similarly to the OnionExpShellModel, this model provides the form factor, P(q), for a multi-shell sphere, 
    5 where the interface between the each neighboring shells can be described by one of a number of functions 
    6 including error, power-law, and exponential functions. This model is to calculate the scattering intensity 
    7 by building a continuous custom SLD profile against the radius of the particle. The SLD profile is composed 
    8 of a flat core, a flat solvent, a number (up to 9 ) flat shells, and the interfacial layers between 
    9 the adjacent flat shells (or core, and solvent) (see below). 
     2This model calculates an empirical functional form for SAS data using 
     3SpericalSLD profile 
     4 
     5Similarly to the OnionExpShellModel, this model provides the form factor, 
     6P(q), for a multi-shell sphere, where the interface between the each neighboring 
     7shells can be described by one of a number of functions including error, 
     8power-law, and exponential functions. 
     9This model is to calculate the scattering intensity by building a continuous 
     10custom SLD profile against the radius of the particle. 
     11The SLD profile is composed of a flat core, a flat solvent, a number (up to 9 ) 
     12flat shells, and the interfacial layers between the adjacent flat shells 
     13(or core, and solvent) (see below). 
    1014 
    1115.. figure:: img/spherical_sld_profile.gif 
    … …  
    1317    Exemplary SLD profile 
    1418 
    15 Unlike the <onion> model (using an analytical integration), 
    16 the interfacial layers here are sub-divided and numerically integrated assuming each of the sub-layers are described 
    17 by a line function. The number of the sub-layer can be given by users by setting the integer values of npts_inter. 
    18 The form factor is normalized by the total volume of the sphere. 
     19Unlike the <onion> model (using an analytical integration), the interfacial 
     20layers here are sub-divided and numerically integrated assuming each of the 
     21sub-layers are described by a line function. 
     22The number of the sub-layer can be given by users by setting the integer values 
     23of npts_inter. The form factor is normalized by the total volume of the sphere. 
    1924 
    2025Definition 
    … …  
    2934    \sum_{\text{flat}_i=0}^N f_{\text{flat}_i} +f_\text{solvent} 
    3035 
    31 For a spherically symmetric particle with a particle density $\rho_x(r)$ the sld function can be defined as: 
     36For a spherically symmetric particle with a particle density $\rho_x(r)$ 
     37the sld function can be defined as: 
    3238 
    3339.. math:: 
    … …  
    3945 
    4046.. math:: 
    41     f_\text{core} = 4 \pi \int_{0}^{r_\text{core}} \rho_\text{core} \frac{\sin(qr)} {qr} r^2 dr = 
     47    f_\text{core} = 4 \pi \int_{0}^{r_\text{core}} \rho_\text{core} 
     48    \frac{\sin(qr)} {qr} r^2 dr = 
    4249    3 \rho_\text{core} V(r_\text{core}) 
    43     \Big[ \frac{\sin(qr_\text{core}) - qr_\text{core} \cos(qr_\text{core})} {qr_\text{core}^3} \Big] 
    44  
    45     f_{\text{inter}_i} = 4 \pi \int_{\Delta t_{ \text{inter}_i } } \rho_{ \text{inter}_i } \frac{\sin(qr)} {qr} r^2 dr 
    46  
    47     f_{\text{shell}_i} = 4 \pi \int_{\Delta t_{ \text{inter}_i } } \rho_{ \text{flat}_i } \frac{\sin(qr)} {qr} r^2 dr = 
    48     3 \rho_{ \text{flat}_i } V ( r_{ \text{inter}_i } + \Delta t_{ \text{inter}_i } ) 
    49     \Big[ \frac{\sin(qr_{\text{inter}_i} + \Delta t_{ \text{inter}_i } ) - q (r_{\text{inter}_i} + 
    50     \Delta t_{ \text{inter}_i }) \cos(q( r_{\text{inter}_i} + \Delta t_{ \text{inter}_i } ) ) } 
     50    \Big[ \frac{\sin(qr_\text{core}) - qr_\text{core} \cos(qr_\text{core})} 
     51    {qr_\text{core}^3} \Big] 
     52 
     53    f_{\text{inter}_i} = 4 \pi \int_{\Delta t_{ \text{inter}_i } } 
     54    \rho_{ \text{inter}_i } \frac{\sin(qr)} {qr} r^2 dr 
     55 
     56    f_{\text{shell}_i} = 4 \pi \int_{\Delta t_{ \text{inter}_i } } 
     57    \rho_{ \text{flat}_i } \frac{\sin(qr)} {qr} r^2 dr = 
     58    3 \rho_{ \text{flat}_i } V ( r_{ \text{inter}_i } + 
     59    \Delta t_{ \text{inter}_i } ) 
     60    \Big[ \frac{\sin(qr_{\text{inter}_i} + \Delta t_{ \text{inter}_i } ) 
     61    - q (r_{\text{inter}_i} + \Delta t_{ \text{inter}_i }) 
     62    \cos(q( r_{\text{inter}_i} + \Delta t_{ \text{inter}_i } ) ) } 
    5163    {q ( r_{\text{inter}_i} + \Delta t_{ \text{inter}_i } )^3 }  \Big] 
    5264    -3 \rho_{ \text{flat}_i } V(r_{ \text{inter}_i }) 
    53     \Big[ \frac{\sin(qr_{\text{inter}_i}) - qr_{\text{flat}_i} \cos(qr_{\text{inter}_i}) } {qr_{\text{inter}_i}^3} \Big] 
    54  
    55     f_\text{solvent} = 4 \pi \int_{r_N}^{\infty} \rho_\text{solvent} \frac{\sin(qr)} {qr} r^2 dr = 
     65    \Big[ \frac{\sin(qr_{\text{inter}_i}) - qr_{\text{flat}_i} 
     66    \cos(qr_{\text{inter}_i}) } {qr_{\text{inter}_i}^3} \Big] 
     67 
     68    f_\text{solvent} = 4 \pi \int_{r_N}^{\infty} \rho_\text{solvent} 
     69    \frac{\sin(qr)} {qr} r^2 dr = 
    5670    3 \rho_\text{solvent} V(r_N) 
    5771    \Big[ \frac{\sin(qr_N) - qr_N \cos(qr_N)} {qr_N^3} \Big] 
    … …  
    6680.. math:: 
    6781    \rho_{{inter}_i} (r) = \begin{cases} 
    68     B \exp\Big( \frac {\pm A(r - r_{\text{flat}_i})} {\Delta t_{ \text{inter}_i }} \Big) +C  & \text{for} A \neq 0 \\ 
    69     B \Big( \frac {(r - r_{\text{flat}_i})} {\Delta t_{ \text{inter}_i }} \Big) +C  & \text{for} A = 0 \\ 
     82    B \exp\Big( \frac {\pm A(r - r_{\text{flat}_i})} 
     83    {\Delta t_{ \text{inter}_i }} \Big) +C  & \text{for} A \neq 0 \\ 
     84    B \Big( \frac {(r - r_{\text{flat}_i})} 
     85    {\Delta t_{ \text{inter}_i }} \Big) +C  & \text{for} A = 0 \\ 
    7086    \end{cases} 
    7187 
    … …  
    7490.. math:: 
    7591    \rho_{{inter}_i} (r) = \begin{cases} 
    76     \pm B \Big( \frac {(r - r_{\text{flat}_i} )} {\Delta t_{ \text{inter}_i }} \Big) ^A  +C  & \text{for} A \neq 0 \\ 
     92    \pm B \Big( \frac {(r - r_{\text{flat}_i} )} {\Delta t_{ \text{inter}_i }} 
     93    \Big) ^A  +C  & \text{for} A \neq 0 \\ 
    7794    \rho_{\text{flat}_{i+1}}  & \text{for} A = 0 \\ 
    7895    \end{cases} 
    … …  
    8299.. math:: 
    83100    \rho_{{inter}_i} (r) = \begin{cases} 
    84     B \text{erf} \Big( \frac { A(r - r_{\text{flat}_i})} {\sqrt{2} \Delta t_{ \text{inter}_i }} \Big) +C  & \text{for} A \neq 0 \\ 
    85     B \Big( \frac {(r - r_{\text{flat}_i} )} {\Delta t_{ \text{inter}_i }} \Big)  +C  & \text{for} A = 0 \\ 
     101    B \text{erf} \Big( \frac { A(r - r_{\text{flat}_i})} 
     102    {\sqrt{2} \Delta t_{ \text{inter}_i }} \Big) +C  & \text{for} A \neq 0 \\ 
     103    B \Big( \frac {(r - r_{\text{flat}_i} )} {\Delta t_{ \text{inter}_i }} 
     104    \Big)  +C  & \text{for} A = 0 \\ 
    86105    \end{cases} 
    87106 
    88 The functions are normalized so that they vary between 0 and 1, and they are constrained such that the SLD 
    89 is continuous at the boundaries of the interface as well as each sub-layers. Thus B and C are determined. 
    90  
    91 Once $\rho_{\text{inter}_i}$ is found at the boundary of the sub-layer of the interface, we can find its contribution 
    92 to the form factor $P(q)$ 
    93  
    94 .. math:: 
    95     f_{\text{inter}_i} = 4 \pi \int_{\Delta t_{ \text{inter}_i } } \rho_{ \text{inter}_i } \frac{\sin(qr)} {qr} r^2 dr = 
     107The functions are normalized so that they vary between 0 and 1, and they are 
     108constrained such that the SLD is continuous at the boundaries of the interface 
     109as well as each sub-layers. Thus B and C are determined. 
     110 
     111Once $\rho_{\text{inter}_i}$ is found at the boundary of the sub-layer of the 
     112interface, we can find its contribution to the form factor $P(q)$ 
     113 
     114.. math:: 
     115    f_{\text{inter}_i} = 4 \pi \int_{\Delta t_{ \text{inter}_i } } 
     116    \rho_{ \text{inter}_i } \frac{\sin(qr)} {qr} r^2 dr = 
    96117    4 \pi \sum_{j=0}^{npts_{\text{inter}_i} -1 } 
    97     \int_{r_j}^{r_{j+1}} \rho_{ \text{inter}_i } (r_j) \frac{\sin(qr)} {qr} r^2 dr \approx 
     118    \int_{r_j}^{r_{j+1}} \rho_{ \text{inter}_i } (r_j) 
     119    \frac{\sin(qr)} {qr} r^2 dr \approx 
    98120 
    99121    4 \pi \sum_{j=0}^{npts_{\text{inter}_i} -1 } \Big[ 
    100     3 ( \rho_{ \text{inter}_i } ( r_{j+1} ) - \rho_{ \text{inter}_i } ( r_{j} ) V ( r_{ \text{sublayer}_j } ) 
    101     \Big[ \frac {r_j^2 \beta_\text{out}^2 \sin(\beta_\text{out}) - (\beta_\text{out}^2-2) \cos(\beta_\text{out}) } 
     122    3 ( \rho_{ \text{inter}_i } ( r_{j+1} ) - \rho_{ \text{inter}_i } 
     123    ( r_{j} ) V ( r_{ \text{sublayer}_j } ) 
     124    \Big[ \frac {r_j^2 \beta_\text{out}^2 \sin(\beta_\text{out}) 
     125    - (\beta_\text{out}^2-2) \cos(\beta_\text{out}) } 
    102126    {\beta_\text{out}^4 } \Big] 
    103127 
    104     - 3 ( \rho_{ \text{inter}_i } ( r_{j+1} ) - \rho_{ \text{inter}_i } ( r_{j} ) V ( r_{ \text{sublayer}_j-1 } ) 
    105     \Big[ \frac {r_{j-1}^2 \sin(\beta_\text{in}) - (\beta_\text{in}^2-2) \cos(\beta_\text{in}) } 
     128    - 3 ( \rho_{ \text{inter}_i } ( r_{j+1} ) - \rho_{ \text{inter}_i } 
     129    ( r_{j} ) V ( r_{ \text{sublayer}_j-1 } ) 
     130    \Big[ \frac {r_{j-1}^2 \sin(\beta_\text{in}) 
     131    - (\beta_\text{in}^2-2) \cos(\beta_\text{in}) } 
    106132    {\beta_\text{in}^4 } \Big] 
    107133 
    … …  
    120146    V(a) = \frac {4\pi}{3}a^3 
    121147 
    122     a_\text{in} ~ \frac{r_j}{r_{j+1} -r_j} \text{, } a_\text{out} ~ \frac{r_{j+1}}{r_{j+1} -r_j} 
     148    a_\text{in} ~ \frac{r_j}{r_{j+1} -r_j} \text{, } a_\text{out} 
     149    ~ \frac{r_{j+1}}{r_{j+1} -r_j} 
    123150 
    124151    \beta_\text{in} = qr_j \text{, } \beta_\text{out} = qr_{j+1} 
    125152 
    126153 
    127 We assume the $\rho_{\text{inter}_i} (r)$ can be approximately linear within a sub-layer $j$ 
     154We assume the $\rho_{\text{inter}_i} (r)$ can be approximately linear 
     155within a sub-layer $j$ 
    128156 
    129157Finally form factor can be calculated by 
    … …  
    131159.. math:: 
    132160 
    133     P(q) = \frac{[f]^2} {V_\text{particle}} \text{where} V_\text{particle} = V(r_{\text{shell}_N}) 
     161    P(q) = \frac{[f]^2} {V_\text{particle}} \text{where} V_\text{particle} 
     162    = V(r_{\text{shell}_N}) 
    134163 
    135164For 2D data the scattering intensity is calculated in the same way as 1D, 
    … …  
    150179 
    151180.. note:: 
    152     The outer most radius is used as the effective radius for S(Q) when $P(Q) * S(Q)$ is applied. 
     181    The outer most radius is used as the effective radius for S(Q) 
     182    when $P(Q) * S(Q)$ is applied. 
    153183 
    154184References 
    155185---------- 
    156 L A Feigin and D I Svergun, Structure Analysis by Small-Angle X-Ray and Neutron Scattering, Plenum Press, New York, (1987) 
     186L A Feigin and D I Svergun, Structure Analysis by Small-Angle X-Ray 
     187and Neutron Scattering, Plenum Press, New York, (1987) 
    157188 
    158189""" 
    … …  
    170201# pylint: disable=bad-whitespace, line-too-long 
    171202#            ["name", "units", default, [lower, upper], "type", "description"], 
    172 parameters = [["n_shells",         "",               1,      [0, 9],         "", "number of shells"], 
    173               ["radius_core",      "Ang",            50.0,   [0, inf],       "", "intern layer thickness"], 
    174               ["sld_core",         "1e-6/Ang^2",     2.07,   [-inf, inf],    "", "sld function flat"], 
    175               ["sld_flat[n]",      "1e-6/Ang^2",     4.06,   [-inf, inf],    "", "sld function flat"], 
    176               ["thick_flat[n]",    "Ang",            100.0,  [0, inf],       "", "flat layer_thickness"], 
    177               ["func_inter[n]",    "",               0,      [0, 4],         "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4"], 
    178               ["thick_inter[n]",   "Ang",            50.0,   [0, inf],       "", "intern layer thickness"], 
    179               ["inter_nu[n]",      "",               2.5,    [-inf, inf],    "", "steepness parameter"], 
    180               ["npts_inter",       "",               35,     [0, 35],        "", "number of points in each sublayer Must be odd number"], 
    181               ["sld_solvent",      "1e-6/Ang^2",     1.0,    [-inf, inf],    "", "sld function solvent"], 
     203parameters = [["n_shells",          "",               1,      [0, 9],         "volume", "number of shells"], 
     204              ["npts_inter",        "",               35,     [0, inf],        "", "number of points in each sublayer Must be odd number"], 
     205              ["radius_core",       "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness"], 
     206              ["sld_core",          "1e-6/Ang^2",     2.07,   [-inf, inf],    "", "sld function flat"], 
     207              ["sld_solvent",       "1e-6/Ang^2",     1.0,    [-inf, inf],    "", "sld function solvent"], 
     208              ["func_inter0",       "",               0,      [0, 4],         "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4"], 
     209              ["thick_inter0",      "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness for core layer"], 
     210              ["nu_inter0",         "",               2.5,    [-inf, inf],    "", "steepness parameter for core layer"], 
     211              ["sld_flat[n_shells]",      "1e-6/Ang^2",     4.06,   [-inf, inf],    "", "sld function flat"], 
     212              ["thick_flat[n_shells]",    "Ang",            100.0,  [0, inf],       "volume", "flat layer_thickness"], 
     213              ["func_inter[n_shells]",    "",               0,      [0, 4],         "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4"], 
     214              ["thick_inter[n_shells]",   "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness"], 
     215              ["nu_inter[n_shells]",      "",               2.5,    [-inf, inf],    "", "steepness parameter"], 
    182216              ] 
    183217# pylint: enable=bad-whitespace, line-too-long 
    184 #source = ["lib/librefl.c",  "lib/sph_j1c.c", "spherical_sld.c"] 
    185  
    186 def Iq(q, *args, **kw): 
    187     return q 
    188  
    189 def Iqxy(qx, *args, **kw): 
    190     return qx 
    191  
    192  
    193 demo = dict( 
    194     n_shells=4, 
    195     scale=1.0, 
    196     solvent_sld=1.0, 
    197     background=0.0, 
    198     npts_inter=35.0, 
    199     ) 
     218source = ["lib/librefl.c",  "lib/sph_j1c.c", "spherical_sld.c"] 
     219 
     220profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 
     221def profile(n_shells, radius_core,  sld_core,  sld_solvent, sld_flat, 
     222            thick_flat, func_inter, thick_inter, nu_inter, npts_inter): 
     223    """ 
     224    Returns shape profile with x=radius, y=SLD. 
     225    """ 
     226 
     227    z = [] 
     228    beta = [] 
     229    z0 = 0 
     230    # two sld points for core 
     231    z.append(0) 
     232    beta.append(sld_core) 
     233    z.append(radius_core) 
     234    beta.append(sld_core) 
     235    z0 += radius_core 
     236 
     237    for i in range(1, n_shells+2): 
     238        dz = thick_inter[i-1]/npts_inter 
     239        # j=0 for interface, j=1 for flat layer 
     240        for j in range(0, 2): 
     241            # interation for sub-layers 
     242            for n_s in range(0, npts_inter+1): 
     243                if j == 1: 
     244                    if i == n_shells+1: 
     245                        break 
     246                    # shift half sub thickness for the first point 
     247                    z0 -= dz#/2.0 
     248                    z.append(z0) 
     249                    #z0 -= dz/2.0 
     250                    z0 += thick_flat[i] 
     251                    sld_i = sld_flat[i] 
     252                    beta.append(sld_flat[i]) 
     253                    dz = 0 
     254                else: 
     255                    nu = nu_inter[i-1] 
     256                    # decide which sld is which, sld_r or sld_l 
     257                    if i == 1: 
     258                        sld_l = sld_core 
     259                    else: 
     260                        sld_l = sld_flat[i-1] 
     261                    if i == n_shells+1: 
     262                        sld_r = sld_solvent 
     263                    else: 
     264                        sld_r = sld_flat[i] 
     265                    # get function type 
     266                    func_idx = func_inter[i-1] 
     267                    # calculate the sld 
     268                    sld_i = intersldfunc(func_idx, npts_inter, n_s, nu, 
     269                                            sld_l, sld_r) 
     270                # append to the list 
     271                z.append(z0) 
     272                beta.append(sld_i) 
     273                z0 += dz 
     274                if j == 1: 
     275                    break 
     276    z.append(z0) 
     277    beta.append(sld_solvent) 
     278    z_ext = z0/5.0 
     279    z.append(z0+z_ext) 
     280    beta.append(sld_solvent) 
     281    # return sld profile (r, beta) 
     282    return np.asarray(z), np.asarray(beta)*1e-6 
     283 
     284def ER(n_shells, radius_core, thick_inter0, thick_inter, thick_flat): 
     285    total_thickness = thick_inter0 
     286    total_thickness += np.sum(thick_inter[:n_shells], axis=0) 
     287    total_thickness += np.sum(thick_flat[:n_shells], axis=0) 
     288    return total_thickness + radius_core 
     289 
     290 
     291demo = { 
     292    "n_shells": 4, 
     293    "npts_inter": 35.0, 
     294    "radius_core": 50.0, 
     295    "sld_core": 2.07, 
     296    "sld_solvent": 1.0, 
     297    "thick_inter0": 50.0, 
     298    "func_inter0": 0, 
     299    "nu_inter0": 2.5, 
     300    "sld_flat":[4.0,3.5,4.0,3.5], 
     301    "thick_flat":[100.0,100.0,100.0,100.0], 
     302    "func_inter":[0,0,0,0], 
     303    "thick_inter":[50.0,50.0,50.0,50.0], 
     304    "nu_inter":[2.5,2.5,2.5,2.5], 
     305    } 
    200306 
    201307#TODO: Not working yet 
    202308tests = [ 
    203309    # Accuracy tests based on content in test/utest_extra_models.py 
    204     [{'npts_iter':35, 
    205         'sld_solv':1, 
    206         'radius_core':50.0, 
    207         'sld_core':2.07, 
    208         'func_inter2':0.0, 
    209         'thick_inter2':50, 
    210         'nu_inter2':2.5, 
    211         'sld_flat2':4, 
    212         'thick_flat2':100, 
    213         'func_inter1':0.0, 
    214         'thick_inter1':50, 
    215         'nu_inter1':2.5, 
    216         'background': 0.0, 
     310    [{"n_shells":4, 
     311        'npts_inter':35, 
     312        "radius_core":50.0, 
     313        "sld_core":2.07, 
     314        "sld_solvent": 1.0, 
     315        "sld_flat":[4.0,3.5,4.0,3.5], 
     316        "thick_flat":[100.0,100.0,100.0,100.0], 
     317        "func_inter":[0,0,0,0], 
     318        "thick_inter":[50.0,50.0,50.0,50.0], 
     319        "nu_inter":[2.5,2.5,2.5,2.5] 
    217320    }, 0.001, 0.001], 
    218321] 

  • sasmodels/models/squarewell.py

    rec45c4f rd2bb604  
    123123""" 
    124124 
    125 Iqxy = """ 
    126     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    127     """ 
    128  
    129125# ER defaults to 0.0 
    130126# VR defaults to 1.0 

  • 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.