Changeset a17fe40 in sasmodels


Ignore:
Timestamp:
Feb 3, 2016 8:56:52 AM (9 years ago)
Author:
krzywon
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:
3a45c2c
Parents:
f12357f (diff), 6d8c359 (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' of ​https://github.com/SasView/sasmodels.git

Location:
sasmodels/models
Files:
10 added
14 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/HayterMSAsq.py

    reb69cce r7f47777  
    9090# odd that the default st has saltconc zero 
    9191demo = dict(effect_radius = 20.75,charge=19.0,volfraction = 0.0192,temperature=318.16,saltconc=0.05,dielectconst=71.08,effect_radius_pd = 0.1,effect_radius_pd_n = 40) 
     92# 
     93# attempt to use same values as old sasview unit test 
     94tests = [ 
     95        [ {'scale': 1.0, 'background' : 0.0, 'effect_radius' : 20.75, 'charge' : 19.0, 'volfraction' : 0.0192, 'temperature' : 298.0, 
     96          'saltconc' : 0,'dielectconst' : 78.0, 'effect_radius_pd' : 0}, [0.0010], [0.0712928]] 
     97        ] 
    9298 

  • sasmodels/models/barbell.py

    reb69cce rdcdf29d  
    100100title = "Cylinder with spherical end caps" 
    101101description = """ 
    102         Calculates the scattering from a barbell-shaped cylinder. That is a sphereocylinder with spherical end caps that have a radius larger than that of the cylinder and the center of the end cap 
    103         radius lies outside of the cylinder. 
    104         Note: As the length of cylinder(bar) -->0,it becomes a dumbbell. And when rad_bar = rad_bell, it is a spherocylinder. 
    105         It must be that rad_bar <(=) rad_bell. 
     102    Calculates the scattering from a barbell-shaped cylinder. 
     103    That is a sphereocylinder with spherical end caps that have a radius larger 
     104    than that of the cylinder and the center of the end cap radius lies outside 
     105    of the cylinder. 
     106    Note: As the length of cylinder(bar) -->0,it becomes a dumbbell. And when 
     107    rad_bar = rad_bell, it is a spherocylinder. 
     108    It must be that rad_bar <(=) rad_bell. 
    106109""" 
    107110category = "shape:cylinder" 
    108  
     111# pylint: disable=bad-whitespace, line-too-long 
    109112#             ["name", "units", default, [lower, upper], "type","description"], 
    110 parameters = [["sld", "4e-6/Ang^2", 4, [-inf, inf], "", "Barbell scattering length density"], 
    111               ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "", "Solvent scattering length density"], 
    112               ["bell_radius", "Ang", 40, [0, inf], "volume", "Spherical bell radius"], 
    113               ["radius", "Ang", 20, [0, inf], "volume", "Cylindrical bar radius"], 
    114               ["length", "Ang", 400, [0, inf], "volume", "Cylinder bar length"], 
    115               ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"], 
    116               ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"], 
     113parameters = [["sld",         "4e-6/Ang^2",   4, [-inf, inf], "",            "Barbell scattering length density"], 
     114              ["solvent_sld", "1e-6/Ang^2",   1, [-inf, inf], "",            "Solvent scattering length density"], 
     115              ["bell_radius", "Ang",         40, [0, inf],    "volume",      "Spherical bell radius"], 
     116              ["radius",      "Ang",         20, [0, inf],    "volume",      "Cylindrical bar radius"], 
     117              ["length",      "Ang",        400, [0, inf],    "volume",      "Cylinder bar length"], 
     118              ["theta",       "degrees",    60, [-inf, inf], "orientation", "In plane angle"], 
     119              ["phi",         "degrees",    60, [-inf, inf], "orientation", "Out of plane angle"], 
    117120             ] 
     121# pylint: enable=bad-whitespace, line-too-long 
    118122 
    119123source = ["lib/J1.c", "lib/gauss76.c", "barbell.c"] 

  • sasmodels/models/bcc.py

    reb69cce rdcdf29d  
    109109title = "Body-centred cubic lattic with paracrystalline distortion" 
    110110description = """ 
    111     Calculates the scattering from a **body-centered cubic lattice** with paracrystalline distortion. Thermal vibrations 
    112     are considered to be negligible, and the size of the paracrystal is infinitely large. Paracrystalline distortion is 
    113     assumed to be isotropic and characterized by a Gaussian distribution. 
     111    Calculates the scattering from a **body-centered cubic lattice** with 
     112    paracrystalline distortion. Thermal vibrations are considered to be 
     113    negligible, and the size of the paracrystal is infinitely large. 
     114    Paracrystalline distortion is assumed to be isotropic and characterized 
     115    by a Gaussian distribution. 
    114116    """ 
    115117category = "shape:paracrystal" 
    116  
     118# pylint: disable=bad-whitespace, line-too-long 
    117119#             ["name", "units", default, [lower, upper], "type","description" ], 
    118 parameters = [["dnn", "Ang", 220, [-inf, inf], "", "Nearest neighbour distance"], 
    119               ["d_factor", "", 0.06, [-inf, inf], "", "Paracrystal distortion factor"], 
    120               ["radius", "Ang", 40, [0, inf], "volume", "Particle radius"], 
    121               ["sld", "1e-6/Ang^2", 4, [-inf, inf], "", "Particle scattering length density"], 
    122               ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "", "Solvent scattering length density"], 
    123               ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"], 
    124               ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"], 
    125               ["psi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"] 
     120parameters = [["dnn",         "Ang",       220,    [-inf, inf], "",            "Nearest neighbour distance"], 
     121              ["d_factor",    "",            0.06, [-inf, inf], "",            "Paracrystal distortion factor"], 
     122              ["radius",      "Ang",        40,    [0, inf],    "volume",      "Particle radius"], 
     123              ["sld",         "1e-6/Ang^2",  4,    [-inf, inf], "",            "Particle scattering length density"], 
     124              ["solvent_sld", "1e-6/Ang^2",  1,    [-inf, inf], "",            "Solvent scattering length density"], 
     125              ["theta",       "degrees",    60,    [-inf, inf], "orientation", "In plane angle"], 
     126              ["phi",         "degrees",    60,    [-inf, inf], "orientation", "Out of plane angle"], 
     127              ["psi",         "degrees",    60,    [-inf, inf], "orientation", "Out of plane angle"] 
    126128             ] 
     129# pylint: enable=bad-whitespace, line-too-long 
    127130 
    128131source = ["lib/J1.c", "lib/gauss150.c", "bcc.c"] 

  • sasmodels/models/broad_peak.py

    reb69cce rdcdf29d  
    6060category = "shape-independent" 
    6161 
     62# pylint: disable=bad-whitespace, line-too-long 
    6263#             ["name", "units", default, [lower, upper], "type", "description"], 
    63 parameters = [["porod_scale", "", 1.0e-05, [-inf, inf], "", "Power law scale factor"], 
    64               ["porod_exp", "", 3.0, [-inf, inf], "", "Exponent of power law"], 
    65               ["lorentz_scale", "", 10.0, [-inf, inf], "", "Scale factor for broad Lorentzian peak"], 
    66               ["lorentz_length", "Ang", 50.0, [-inf, inf], "", "Lorentzian screening length"], 
    67               ["peak_pos", "1/Ang", 0.1, [-inf, inf], "", "Peak postion in q"], 
    68               ["lorentz_exp", "", 2.0, [-inf, inf], "", "exponent of Lorentz function"], 
     64parameters = [["porod_scale",    "", 1.0e-05, [-inf, inf], "", "Power law scale factor"], 
     65              ["porod_exp",      "",      3.0, [-inf, inf], "", "Exponent of power law"], 
     66              ["lorentz_scale",  "",    10.0, [-inf, inf], "", "Scale factor for broad Lorentzian peak"], 
     67              ["lorentz_length", "Ang",  50.0, [-inf, inf], "", "Lorentzian screening length"], 
     68              ["peak_pos",       "1/Ang", 0.1, [-inf, inf], "", "Peak position in q"], 
     69              ["lorentz_exp",    "",      2.0, [-inf, inf], "", "Exponent of Lorentz function"], 
    6970             ] 
     71# pylint: enable=bad-whitespace, line-too-long 
    7072 
    71 def Iq(q, porod_scale, porod_exp, lorentz_scale, lorentz_length, peak_pos, lorentz_exp): 
     73def Iq(q, 
     74       porod_scale=1.0e-5, 
     75       porod_exp=3.0, 
     76       lorentz_scale=10.0, 
     77       lorentz_length=50.0, 
     78       peak_pos=0.1, 
     79       lorentz_exp=2.0): 
     80    """ 
     81    :param q:              Input q-value 
     82    :param porod_scale:    Power law scale factor 
     83    :param porod_exp:      Exponent of power law 
     84    :param lorentz_scale:  Scale factor for broad Lorentzian peak 
     85    :param lorentz_length: Lorentzian screening length 
     86    :param peak_pos:       Peak position in q 
     87    :param lorentz_exp:    Exponent of Lorentz function 
     88    :return:               Calculated intensity 
     89    """ 
     90 
    7291    inten = (porod_scale / q ** porod_exp + lorentz_scale 
    7392             / (1.0 + (abs(q - peak_pos) * lorentz_length) ** lorentz_exp)) 
    … …  
    7695 
    7796def Iqxy(qx, qy, *args): 
     97    """ 
     98    :param qx:   Input q_x-value 
     99    :param qy:   Input q_y-value 
     100    :param args: Remaining arguments 
     101    :return:     2D-Intensity 
     102    """ 
    78103    return Iq(sqrt(qx ** 2 + qy ** 2), *args) 
     104 
    79105Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values 
    80  
    81106 
    82107demo = dict(scale=1, background=0, 

  • sasmodels/models/capped_cylinder.py

    reb69cce rdcdf29d  
    125125""" 
    126126category = "shape:cylinder" 
     127# pylint: disable=bad-whitespace, line-too-long 
     128#             ["name", "units", default, [lower, upper], "type", "description"], 
     129parameters = [["sld",         "1e-6/Ang^2", 4, [-inf, inf], "",       "Cylinder scattering length density"], 
     130              ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "",       "Solvent scattering length density"], 
     131              ["radius",      "Ang",       20, [0, inf],    "volume", "Cylinder radius"], 
    127132 
    128 #             ["name", "units", default, [lower, upper], "type", "description"], 
    129 parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "", 
    130                "Cylinder scattering length density"], 
    131               ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "", 
    132                "Solvent scattering length density"], 
    133               ["radius", "Ang", 20, [0, inf], "volume", "Cylinder radius"], 
    134133              # TODO: use an expression for cap radius with fixed bounds. 
    135134              # The current form requires cap radius R bigger than cylinder radius r. 
    … …  
    141140              # in the capped cylinder, and zero for no bar in the barbell model.  In 
    142141              # both models, one would be a pill. 
    143               ["cap_radius", "Ang", 20, [0, inf], "volume", "Cap radius"], 
    144               ["length", "Ang", 400, [0, inf], "volume", "Cylinder length"], 
    145               ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"], 
    146               ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"], 
     142              ["cap_radius", "Ang",     20, [0, inf],    "volume", "Cap radius"], 
     143              ["length",     "Ang",    400, [0, inf],    "volume", "Cylinder length"], 
     144              ["theta",      "degrees", 60, [-inf, inf], "orientation", "In plane angle"], 
     145              ["phi",        "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"], 
    147146             ] 
     147# pylint: enable=bad-whitespace, line-too-long 
    148148 
    149149source = ["lib/J1.c", "lib/gauss76.c", "capped_cylinder.c"] 
    … …  
    159159            phi_pd=15, phi_pd_n=1) 
    160160oldname = 'CappedCylinderModel' 
    161 oldpars = dict(sld='sld_capcyl', solvent_sld='sld_solv', 
    162                length='len_cyl', radius='rad_cyl', cap_radius='rad_cap') 
     161oldpars = dict(sld='sld_capcyl', 
     162               solvent_sld='sld_solv', 
     163               length='len_cyl', 
     164               radius='rad_cyl', 
     165               cap_radius='rad_cap') 

  • sasmodels/models/hardsphere.py

    r7ed702f r7f47777  
    106106oldpars = dict() 
    107107 
     108tests = [ 
     109        [ {'scale': 1.0, 'background' : 0.0, 'effect_radius' : 50.0, 'volfraction' : 0.2, 
     110           'effect_radius_pd' : 0}, [0.001], [0.209128]] 
     111        ] 
     112 

  • sasmodels/models/lamellarCaille.py

    rd18f8a8 r7f47777  
    122122oldpars = dict(thickness='delta', Nlayers='N_plates', Caille_parameter='caille', 
    123123               sld='sld_bi',solvent_sld='sld_sol') 
     124# 
     125tests = [ 
     126        [ {'scale': 1.0, 'background' : 0.0, 'thickness' : 30.,'Nlayers' : 20.0, 'spacing' : 400., 
     127            'Caille_parameter' : 0.1, 'sld' : 6.3, 'solvent_sld' : 1.0, 
     128            'thickness_pd' : 0.0, 'spacing_pd' : 0.0 }, [0.001], [28895.13397]] 
     129        ] 

  • sasmodels/models/lamellarCailleHG.py

    rd4666ca r7f47777  
    139139 
    140140oldname = 'LamellarPSHGModel' 
     141 
    141142oldpars = dict( 
    142143    tail_length='deltaT', head_length='deltaH', Nlayers='n_plates', 
    143144    Caille_parameter='caille', sld='sld_tail', head_sld='sld_head', 
    144145    solvent_sld='sld_solvent') 
     146# 
     147tests = [ 
     148        [ {'scale': 1.0, 'background' : 0.0, 'tail_length' : 10.0, 'head_length' : 2.0,'Nlayers' : 30.0, 'spacing' : 40., 
     149            'Caille_parameter' : 0.001, 'sld' : 0.4, 'head_sld' : 2.0, 'solvent_sld' : 6.0, 
     150            'tail_length_pd' : 0.0, 'head_length_pd' : 0.0, 'spacing_pd' : 0.0 }, [0.001], [6838238.571488]] 
     151        ] 

  • sasmodels/models/lamellarFFHG.py

    reb69cce r7f47777  
    121121oldpars = dict(head_length='h_thickness', tail_length='t_length', 
    122122               sld='sld_tail', head_sld='sld_head', solvent_sld='sld_solvent') 
     123# 
     124tests = [ 
     125        [ {'scale': 1.0, 'background' : 0.0, 'tail_length' : 15.0, 'head_length' : 10.0,'sld' : 0.4, 
     126         'head_sld' : 3.0, 'solvent_sld' : 6.0, }, [0.001], [653143.9209]] 
     127        ] 
    123128 
     129 

  • sasmodels/models/lamellarPC.py

    rd18f8a8 r7f47777  
    148148oldpars = dict(spacing_polydisp='pd_spacing', sld='sld_layer', 
    149149               solvent_sld='sld_solvent') 
     150# 
     151tests = [ 
     152        [ {'scale': 1.0, 'background' : 0.0, 'thickness' : 33.,'Nlayers' : 20.0, 'spacing' : 250., 'spacing_polydisp' : 0.2, 
     153            'sld' : 1.0, 'solvent_sld' : 6.34, 
     154            'thickness_pd' : 0.0, 'thickness_pd_n' : 40 }, [0.001, 0.215268], [21829.3, 0.00487686]] 
     155        ] 

  • sasmodels/models/stickyhardsphere.py

    rd18f8a8 r7f47777  
    182182demo = dict(effect_radius=200, volfraction=0.2, perturb=0.05, 
    183183            stickiness=0.2, effect_radius_pd=0.1, effect_radius_pd_n=40) 
     184# 
     185tests = [ 
     186        [ {'scale': 1.0, 'background' : 0.0, 'effect_radius' : 50.0, 'perturb' : 0.05, 'stickiness' : 0.2, 'volfraction' : 0.1, 
     187           'effect_radius_pd' : 0}, [0.001], [1.09718]] 
     188        ] 
    184189 
     190 

  • sasmodels/models/lib/Si.c

    r7d256c8 rf12357f  
    1313 
    1414                // Explicitly writing factorial values triples the speed of the calculation 
    15                 out_cos = 1/x - 2/pow(x,3) + 24/pow(x,5) - 720/pow(x,7) + 40320/pow(x,9); 
    16                 out_sin = 1/x - 6/pow(x,4) + 120/pow(x,6) - 5040/pow(x,8) + 362880/pow(x,10); 
     15                out_cos = 1./x - 2./pow(x,3) + 24./pow(x,5) - 720./pow(x,7); 
     16                out_sin = 1./pow(x,2) - 6./pow(x,4) + 120./pow(x,6) - 5040./pow(x,8); 
    1717 
    1818                out -= cos(x) * out_cos; 
    … …  
    2222 
    2323        // Explicitly writing factorial values triples the speed of the calculation 
    24         out = x - pow(x, 3)/18 + pow(x,5)/600 - pow(x,7)/35280 + pow(x,9)/3265920; 
     24        out = x - pow(x, 3)/18. + pow(x,5)/600. - pow(x,7)/35280. + pow(x,9)/3265920. - pow(x,11)/439084800.; 
    2525 
    26         //printf ("Si=%g %g\n", x, out); 
    2726        return out; 
    2827} 

  • sasmodels/models/pearl_necklace.c

    rcf404cb rf12357f  
    1212        double string_thickness, double number_of_pearls, double sld,  
    1313        double string_sld, double solvent_sld); 
    14          
    15 double ER(double radius, double edge_separation, 
    16         double string_thickness, double number_of_pearls); 
    17 double VR(double radius, double edge_separation, 
    18         double string_thickness, double number_of_pearls); 
    1914 
    2015// From Igor library 
    … …  
    2217        double num_pearls, double sld_pearl, double sld_string, double sld_solv) 
    2318{ 
     19        //relative slds 
    2420        double contrast_pearl = sld_pearl - sld_solv; 
    2521        double contrast_string = sld_string - sld_solv; 
     22         
     23        // number of string segments 
     24        num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 
     25        double num_strings = num_pearls - 1.0; 
     26         
     27        //Pi 
     28        double pi = 4.0*atan(1.0); 
     29         
     30        // center to center distance between the neighboring pearls 
     31        double A_s = edge_separation + 2.0 * radius; 
     32         
     33        // Repeated Calculations 
     34        double sincasq = sinc(q*A_s); 
     35        double oneminussinc = 1 - sincasq; 
     36        double q_r = q * radius; 
     37        double q_edge = q * edge_separation; 
     38         
     39        // each volume 
     40        double string_vol = edge_separation * pi * thick_string * thick_string / 4.0; 
     41        double pearl_vol = 4.0 / 3.0 * pi * radius * radius * radius; 
    2642 
    2743        //total volume 
    28         double pi = 4.0*atan(1.0); 
    29         double tot_vol = form_volume(radius, edge_separation, thick_string, num_pearls); 
    30         double string_vol = edge_separation * pi * pow((thick_string / 2.0), 2); 
    31         double pearl_vol = 4.0 /3.0 * pi * pow(radius, 3); 
    32         double num_strings = num_pearls - 1; 
    33         //each mass 
     44        double tot_vol; 
     45        //each masses 
    3446        double m_r= contrast_string * string_vol; 
    3547        double m_s= contrast_pearl * pearl_vol; 
    3648        double psi, gamma, beta; 
    3749        //form factors 
    38         double sss; //pearls 
    39         double srr; //strings 
    40         double srs; //cross 
    41         double A_s, srr_1, srr_2, srr_3; 
     50        double sss, srr, srs; //cross 
     51        double srr_1, srr_2, srr_3; 
    4252        double form_factor; 
     53        tot_vol = num_strings * string_vol; 
     54        tot_vol += num_pearls * pearl_vol; 
    4355 
    4456        //sine functions of a pearl 
    45         psi = sin(q*radius); 
    46         psi -= q * radius * cos(q * radius); 
    47         psi /= pow((q * radius), 3); 
     57        psi = sin(q_r); 
     58        psi -= q_r * cos(q_r); 
    4859        psi *= 3.0; 
     60        psi /= q_r * q_r * q_r; 
    4961 
    5062        // Note take only 20 terms in Si series: 10 terms may be enough though. 
    51         gamma = Si(q* edge_separation); 
    52         gamma /= (q* edge_separation); 
    53         beta = Si(q * (edge_separation + radius)); 
    54         beta -= Si(q * radius); 
    55         beta /= (q* edge_separation); 
    56  
    57         // center to center distance between the neighboring pearls 
    58         A_s = edge_separation + 2.0 * radius; 
     63        gamma = Si(q_edge); 
     64        gamma /= (q_edge); 
     65        beta = Si(q * (A_s - radius)); 
     66        beta -= Si(q_r); 
     67        beta /= q_edge; 
    5968 
    6069        // form factor for num_pearls 
    61         sss = 1.0 - pow(sinc(q*A_s), num_pearls ); 
    62         sss /= pow((1.0-sinc(q*A_s)), 2); 
    63         sss *= -sinc(q*A_s); 
    64         sss -= num_pearls/2.0; 
    65         sss += num_pearls/(1.0-sinc(q*A_s)); 
    66         sss *= 2.0 * pow((m_s*psi), 2); 
     70        sss = 1.0 - pow(sincasq, num_pearls); 
     71        sss /= oneminussinc * oneminussinc; 
     72        sss *= -sincasq; 
     73        sss -= num_pearls / 2.0; 
     74        sss += num_pearls / oneminussinc; 
     75        sss *= 2.0 * m_s * psi * m_s * psi; 
    6776 
    6877        // form factor for num_strings (like thin rods) 
    69         srr_1 = -pow(sinc(q*edge_separation/2.0), 2); 
     78        srr_1 = -sinc(q_edge/2.0) * sinc(q_edge/2.0); 
    7079 
    7180        srr_1 += 2.0 * gamma; 
    7281        srr_1 *= num_strings; 
    73         srr_2 = 2.0/(1.0-sinc(q*A_s)); 
    74         srr_2 *= num_strings * pow(beta, 2); 
    75         srr_3 = 1.0 - pow(sinc(q*A_s), num_strings); 
    76         srr_3 /= pow((1.0-sinc(q*A_s)), 2); 
    77         srr_3 *= -2.0 * pow(beta, 2); 
     82        srr_2 = 2.0/oneminussinc; 
     83        srr_2 *= num_strings; 
     84        srr_2 *= beta * beta; 
     85        srr_3 = 1.0 - pow(sincasq, num_strings); 
     86        srr_3 /= oneminussinc * oneminussinc; 
     87        srr_3 *= beta * beta; 
     88        srr_3 *= -2.0; 
    7889 
    7990        // total srr 
    8091        srr = srr_1 + srr_2 + srr_3; 
    81         srr *= pow(m_r, 2); 
     92        srr *= m_r * m_r; 
    8293 
    8394        // form factor for correlations 
    8495        srs = 1.0; 
    85         srs -= pow(sinc(q*A_s), num_strings); 
    86         srs /= pow((1.0-sinc(q*A_s)), 2); 
    87         srs *= -sinc(q*A_s); 
    88         srs += (num_strings/(1.0-sinc(q*A_s))); 
    89         srs *= 4.0 * (m_r * m_s * beta * psi); 
     96        srs -= pow(sincasq, num_strings); 
     97        srs /= oneminussinc * oneminussinc; 
     98        srs *= -sincasq; 
     99        srs += num_strings/oneminussinc; 
     100        srs *= 4.0; 
     101        srs *= (m_r * m_s * beta * psi); 
    90102 
    91103        form_factor = sss + srr + srs; 
    … …  
    103115        double number_of_strings = number_of_pearls - 1.0; 
    104116         
    105         double string_vol = edge_separation * pi * pow((string_thickness / 2.0), 2); 
    106         double pearl_vol = 4.0 /3.0 * pi * pow(radius, 3); 
     117        double string_vol = edge_separation * pi * string_thickness * string_thickness / 4.0; 
     118        double pearl_vol = 4.0 / 3.0 * pi * radius * radius * radius; 
    107119 
    108120        total_vol = number_of_strings * string_vol; 
    … …  
    124136        double string_sld, double solvent_sld) 
    125137{ 
    126         double value = 0.0; 
    127         double tot_vol = 0.0; 
     138        double value, tot_vol; 
     139         
     140        if (string_thickness >= radius || number_of_pearls <= 0) { 
     141                return NAN; 
     142        } 
    128143         
    129144        value = _pearl_necklace_kernel(q, radius, edge_separation, string_thickness, 
    … …  
    142157                sld, string_sld, solvent_sld)); 
    143158} 
    144  
    145  
    146 double ER(double radius, double edge_separation, 
    147         double string_thickness, double number_of_pearls) 
    148 { 
    149         double tot_vol = form_volume(radius, edge_separation, string_thickness, number_of_pearls); 
    150         double pi = 4.0*atan(1.0); 
    151  
    152     double rad_out = pow((3.0*tot_vol/4.0/pi), 0.33333); 
    153      
    154     return rad_out; 
    155 } 
    156 double VR(double radius, double edge_separation, 
    157         double string_thickness, double number_of_pearls) 
    158 { 
    159         return 1.0; 
    160 } 

  • sasmodels/models/pearl_necklace.py

    r841753c rf12357f  
    5858""" 
    5959 
    60 from numpy import inf 
     60from numpy import inf, pi 
    6161 
    6262name = "pearl_necklace" 
    … …  
    9696 
    9797source = ["lib/Si.c", "pearl_necklace.c"] 
     98# new flag to let the compiler know to never use single precision 
     99single = False 
     100 
     101def volume(radius, edge_separation, string_thickness, number_of_pearls): 
     102    """ 
     103    Calculates the total particle volume of the necklace. 
     104    Redundant with form_volume. 
     105    """ 
     106    number_of_strings = number_of_pearls - 1.0 
     107    string_vol = edge_separation * pi * pow((string_thickness / 2.0), 2.0) 
     108    pearl_vol = 4.0 /3.0 * pi * pow(radius, 3.0) 
     109    total_vol = number_of_strings * string_vol 
     110    total_vol += number_of_pearls * pearl_vol 
     111    return total_vol 
     112 
     113def ER(radius, edge_separation, string_thickness, number_of_pearls): 
     114    """ 
     115    Calculation for effective radius. 
     116    """ 
     117    tot_vol = volume(radius, edge_separation, string_thickness, number_of_pearls) 
     118    rad_out = pow((3.0*tot_vol/4.0/pi), 0.33333) 
     119    return rad_out 
    98120 
    99121# parameters for demo 
    … …  
    114136               string_thickness='thick_string', sld='sld_pearl', 
    115137               string_sld='sld_string', edge_separation='edge_separation') 
     138 
     139tests = [[{}, 0.001, 17380.245], [{}, 'ER', 115.39502]] 
Note: See TracChangeset for help on using the changeset viewer.