Changeset 7891a2a in sasmodels for sasmodels/models


Ignore:
Timestamp:
May 9, 2016 9:36:28 AM (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:
a326b8e
Parents:
24d5b30 (diff), 7bf4757 (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:
Paul Kienzle <pkienzle@…> (05/09/16 09:34:57)
git-committer:
Paul Kienzle <pkienzle@…> (05/09/16 09:36:28)
Message:

Merge branch 'master' into polydisp

Location:
sasmodels/models
Files:
2 added
23 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/__init__.py

    r9ccf580 r306e354  
    11""" 
    2     1D Modeling for SAS 
     21D Modeling for SAS 
    33""" 
    4 #from sas.models import * 
    54 
    6 import os 
    7 from distutils.filelist import findall 
    8  
    9 __version__ = "2.1.0" 
    10  
    11 def get_data_path(media): 
    12     """ 
    13     """ 
    14     # Check for data path in the package 
    15     path = os.path.join(os.path.dirname(__file__), media) 
    16     if os.path.isdir(path): 
    17         return path 
    18  
    19     # Check for data path next to exe/zip file. 
    20     # If we are inside a py2exe zip file, we need to go up 
    21     # to get to the directory containing  
    22     # the media for this module 
    23     path = os.path.dirname(__file__) 
    24     #Look for maximum n_dir up of the current dir to find media 
    25     n_dir = 12 
    26     for i in range(n_dir): 
    27         path, _ = os.path.split(path) 
    28         media_path = os.path.join(path, media) 
    29         if os.path.isdir(media_path): 
    30             module_media_path = os.path.join(media_path, 'models_media') 
    31             if os.path.isdir(module_media_path): 
    32                 return module_media_path 
    33             return media_path 
    34     
    35     raise RuntimeError('Could not find models media files') 
    36  
    37 def data_files(): 
    38     """ 
    39     Return the data files associated with media. 
    40      
    41     The format is a list of (directory, [files...]) pairs which can be 
    42     used directly in setup(...,data_files=...) for setup.py. 
    43  
    44     """ 
    45     data_files = [] 
    46     path_img = get_data_path(media=os.path.join("sasmodels","sasmodels","models","img")) 
    47     #path_img = get_data_path(media="img") 
    48     im_list = findall(path_img) 
    49     #for f in findall(path): 
    50     #    if os.path.isfile(f) and f not in im_list: 
    51     #        data_files.append(('media/models_media', [f])) 
    52      
    53     for f in im_list: 
    54         data_files.append(('media/models_media/img', [f])) 
    55  
    56     # Add precompiled models dlls to the top level directory "models" 
    57     path_img = get_data_path(media=os.path.join("..", "..", "models")) 
    58     is_list = findall(path_img) 
    59     for f in is_list: 
    60         if os.path.splitext(f)[1] == ".so": 
    61             data_files.append(('models', [f])) 
    62  
    63     return data_files 
    64  

  • 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 r24d5b30  
    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""" 
    159190 
     191import numpy as np 
    160192from numpy import inf 
    161193 
    … …  
    170202# pylint: disable=bad-whitespace, line-too-long 
    171203#            ["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"], 
     204parameters = [["n_shells",          "",               1,      [0, 9],         "volume", "number of shells"], 
     205              ["npts_inter",        "",               35,     [0, inf],        "", "number of points in each sublayer Must be odd number"], 
     206              ["radius_core",       "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness"], 
     207              ["sld_core",          "1e-6/Ang^2",     2.07,   [-inf, inf],    "", "sld function flat"], 
     208              ["sld_solvent",       "1e-6/Ang^2",     1.0,    [-inf, inf],    "", "sld function solvent"], 
     209              ["func_inter0",       "",               0,      [0, 4],         "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4"], 
     210              ["thick_inter0",      "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness for core layer"], 
     211              ["nu_inter0",         "",               2.5,    [-inf, inf],    "", "steepness parameter for core layer"], 
     212              ["sld_flat[n_shells]",      "1e-6/Ang^2",     4.06,   [-inf, inf],    "", "sld function flat"], 
     213              ["thick_flat[n_shells]",    "Ang",            100.0,  [0, inf],       "volume", "flat layer_thickness"], 
     214              ["func_inter[n_shells]",    "",               0,      [0, 4],         "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4"], 
     215              ["thick_inter[n_shells]",   "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness"], 
     216              ["nu_inter[n_shells]",      "",               2.5,    [-inf, inf],    "", "steepness parameter"], 
    182217              ] 
    183218# 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     ) 
     219source = ["lib/librefl.c",  "lib/sph_j1c.c", "spherical_sld.c"] 
     220 
     221profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 
     222def profile(n_shells, radius_core,  sld_core,  sld_solvent, sld_flat, 
     223            thick_flat, func_inter, thick_inter, nu_inter, npts_inter): 
     224    """ 
     225    Returns shape profile with x=radius, y=SLD. 
     226    """ 
     227 
     228    z = [] 
     229    beta = [] 
     230    z0 = 0 
     231    # two sld points for core 
     232    z.append(0) 
     233    beta.append(sld_core) 
     234    z.append(radius_core) 
     235    beta.append(sld_core) 
     236    z0 += radius_core 
     237 
     238    for i in range(1, n_shells+2): 
     239        dz = thick_inter[i-1]/npts_inter 
     240        # j=0 for interface, j=1 for flat layer 
     241        for j in range(0, 2): 
     242            # interation for sub-layers 
     243            for n_s in range(0, npts_inter+1): 
     244                if j == 1: 
     245                    if i == n_shells+1: 
     246                        break 
     247                    # shift half sub thickness for the first point 
     248                    z0 -= dz#/2.0 
     249                    z.append(z0) 
     250                    #z0 -= dz/2.0 
     251                    z0 += thick_flat[i] 
     252                    sld_i = sld_flat[i] 
     253                    beta.append(sld_flat[i]) 
     254                    dz = 0 
     255                else: 
     256                    nu = nu_inter[i-1] 
     257                    # decide which sld is which, sld_r or sld_l 
     258                    if i == 1: 
     259                        sld_l = sld_core 
     260                    else: 
     261                        sld_l = sld_flat[i-1] 
     262                    if i == n_shells+1: 
     263                        sld_r = sld_solvent 
     264                    else: 
     265                        sld_r = sld_flat[i] 
     266                    # get function type 
     267                    func_idx = func_inter[i-1] 
     268                    # calculate the sld 
     269                    sld_i = intersldfunc(func_idx, npts_inter, n_s, nu, 
     270                                            sld_l, sld_r) 
     271                # append to the list 
     272                z.append(z0) 
     273                beta.append(sld_i) 
     274                z0 += dz 
     275                if j == 1: 
     276                    break 
     277    z.append(z0) 
     278    beta.append(sld_solvent) 
     279    z_ext = z0/5.0 
     280    z.append(z0+z_ext) 
     281    beta.append(sld_solvent) 
     282    # return sld profile (r, beta) 
     283    return np.asarray(z), np.asarray(beta)*1e-6 
     284 
     285def ER(n_shells, radius_core, thick_inter0, thick_inter, thick_flat): 
     286    total_thickness = thick_inter0 
     287    total_thickness += np.sum(thick_inter[:n_shells], axis=0) 
     288    total_thickness += np.sum(thick_flat[:n_shells], axis=0) 
     289    return total_thickness + radius_core 
     290 
     291 
     292demo = { 
     293    "n_shells": 4, 
     294    "npts_inter": 35.0, 
     295    "radius_core": 50.0, 
     296    "sld_core": 2.07, 
     297    "sld_solvent": 1.0, 
     298    "thick_inter0": 50.0, 
     299    "func_inter0": 0, 
     300    "nu_inter0": 2.5, 
     301    "sld_flat":[4.0,3.5,4.0,3.5], 
     302    "thick_flat":[100.0,100.0,100.0,100.0], 
     303    "func_inter":[0,0,0,0], 
     304    "thick_inter":[50.0,50.0,50.0,50.0], 
     305    "nu_inter":[2.5,2.5,2.5,2.5], 
     306    } 
    200307 
    201308#TODO: Not working yet 
    202309tests = [ 
    203310    # 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, 
     311    [{"n_shells":4, 
     312        'npts_inter':35, 
     313        "radius_core":50.0, 
     314        "sld_core":2.07, 
     315        "sld_solvent": 1.0, 
     316        "sld_flat":[4.0,3.5,4.0,3.5], 
     317        "thick_flat":[100.0,100.0,100.0,100.0], 
     318        "func_inter":[0,0,0,0], 
     319        "thick_inter":[50.0,50.0,50.0,50.0], 
     320        "nu_inter":[2.5,2.5,2.5,2.5] 
    217321    }, 0.001, 0.001], 
    218322] 

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