Changeset 6196487 in sasmodels


Ignore:
Timestamp:
Mar 8, 2017 6:24:32 AM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
ceb8a8c
Parents:
d24e390 (diff), 9ae85f0 (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 ticket-835

Location:
sasmodels
Files:
28 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/compare.py

    rfe25eda rd504bcd  
    340340        if pars['radius'] < pars['thick_string']: 
    341341            pars['radius'], pars['thick_string'] = pars['thick_string'], pars['radius'] 
    342         pars['num_pearls'] = math.ceil(pars['num_pearls']) 
    343342        pass 
    344343 
    … …  
    353352        for c in '1234': 
    354353            pars['Phi'+c] /= total 
    355  
    356     elif name == 'stacked_disks': 
    357         pars['n_stacking'] = math.ceil(pars['n_stacking']) 
    358354 
    359355def parlist(model_info, pars, is2d): 

  • sasmodels/models/core_multi_shell.c

    r925ad6e rc3ccaec  
    88} 
    99 
    10 double 
    11 form_volume(double core_radius, double n, double thickness[]); 
    12 double 
    13 form_volume(double core_radius, double n, double thickness[]) 
     10static double 
     11form_volume(double core_radius, double fp_n, double thickness[]) 
    1412{ 
    1513  double r = core_radius; 
     14  int n = (int)(fp_n+0.5); 
    1615  for (int i=0; i < n; i++) { 
    1716    r += thickness[i]; 
    … …  
    2019} 
    2120 
    22 double 
     21static double 
    2322Iq(double q, double core_sld, double core_radius, 
    24    double solvent_sld, double num_shells, double sld[], double thickness[]); 
    25 double 
    26 Iq(double q, double core_sld, double core_radius, 
    27    double solvent_sld, double num_shells, double sld[], double thickness[]) 
     23   double solvent_sld, double fp_n, double sld[], double thickness[]) 
    2824{ 
    29   const int n = (int)ceil(num_shells); 
     25  const int n = (int)(fp_n+0.5); 
    3026  double f, r, last_sld; 
    3127  r = core_radius; 

  • sasmodels/models/core_multi_shell.py

    r925ad6e r5a0b3d7  
    107107    Returns the SLD profile *r* (Ang), and *rho* (1e-6/Ang^2). 
    108108    """ 
     109    n = int(n+0.5) 
    109110    z = [] 
    110111    rho = [] 
    … …  
    133134def ER(radius, n, thickness): 
    134135    """Effective radius""" 
    135     n = int(n[0])  # n cannot be polydisperse 
     136    n = int(n[0]+0.5)  # n is a control parameter and is not polydisperse 
    136137    return np.sum(thickness[:n], axis=0) + radius 
    137138 

  • sasmodels/models/flexible_cylinder.c

    r592343f rc3ccaec  
    11static double 
    2 form_volume(length, kuhn_length, radius) 
     2form_volume(double length, double kuhn_length, double radius) 
    33{ 
    44    return 1.0; 

  • sasmodels/models/lamellar_hg_stack_caille.c

    ra807206 r1c6286d  
    33*/ 
    44 
    5 double Iq(double qval, 
    6       double length_tail, 
    7       double length_head, 
    8       double Nlayers,  
    9       double dd, 
    10       double Cp, 
    11       double tail_sld, 
    12       double head_sld, 
    13       double solvent_sld); 
    14  
    15 double Iq(double qval, 
    16       double length_tail, 
    17       double length_head, 
    18       double Nlayers,  
    19       double dd, 
    20       double Cp, 
    21       double tail_sld, 
    22       double head_sld, 
    23       double solvent_sld) 
     5static double 
     6Iq(double qval, 
     7   double length_tail, 
     8   double length_head, 
     9   double fp_Nlayers, 
     10   double dd, 
     11   double Cp, 
     12   double tail_sld, 
     13   double head_sld, 
     14   double solvent_sld) 
    2415{ 
    25   double NN;   //local variables of coefficient wave 
     16  int Nlayers = (int)(fp_Nlayers+0.5);    //cast to an integer for the loop 
    2617  double inten,Pq,Sq,alpha,temp,t2; 
    2718  //double dQ, dQDefault, t1, t3; 
    28   int ii,NNint; 
    2919  // from wikipedia 0.577215664901532860606512090082402431042159335 
    3020  const double Euler = 0.577215664901533;   // Euler's constant, increased sig figs for new models Feb 2015 
    … …  
    3222  //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015 
    3323 
    34   NN = trunc(Nlayers);    //be sure that NN is an integer 
    35    
    36   Pq = (head_sld-solvent_sld)*(sin(qval*(length_head+length_tail))-sin(qval*length_tail)) + 
    37               (tail_sld-solvent_sld)*sin(qval*length_tail); 
     24  Pq = (head_sld-solvent_sld)*(sin(qval*(length_head+length_tail))-sin(qval*length_tail)) 
     25       + (tail_sld-solvent_sld)*sin(qval*length_tail); 
    3826  Pq *= Pq; 
    3927  Pq *= 4.0/(qval*qval); 
    4028 
    41   NNint = (int)NN;    //cast to an integer for the loop 
    42   ii=0; 
    4329  Sq = 0.0; 
    44   for(ii=1;ii<=(NNint-1);ii+=1) { 
    45  
    46     //fii = (double)ii;   //do I really need to do this? - unused variable, removed 18Feb2015 
    47  
     30  for(int ii=1; ii < Nlayers; ii++) { 
    4831    temp = 0.0; 
    4932    alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler); 
    … …  
    5235    //t3 = dQ*dQ*dd*dd*ii*ii; 
    5336 
    54     temp = 1.0-ii/NN; 
     37    temp = 1.0-(double)ii/(double)Nlayers; 
    5538    //temp *= cos(dd*qval*ii/(1.0+t1)); 
    5639    temp *= cos(dd*qval*ii); 
    5740    //if (temp < 0) printf("q=%g: ii=%d, cos(dd*q*ii)=cos(%g) < 0\n",qval,ii,dd*qval*ii); 
    5841    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); 
    59     temp *= exp(-t2/2.0 ); 
     42    temp *= exp(-t2/2.0); 
    6043    //temp /= sqrt(1.0+t1); 
    6144 
    … …  
    7154 
    7255  inten *= 1.0e-04;   // 1/A to 1/cm 
    73   return(inten); 
     56  return inten; 
    7457} 
    7558 

  • sasmodels/models/lamellar_hg_stack_caille.py

    r7c57861 ra57b31d  
    9898    ["length_head", "Ang", 2, [0, inf], "volume", 
    9999     "head thickness"], 
    100     ["Nlayers", "", 30, [0, inf], "", 
     100    ["Nlayers", "", 30, [1, inf], "", 
    101101     "Number of layers"], 
    102102    ["d_spacing", "Ang", 40., [0.0, inf], "volume", 

  • sasmodels/models/lamellar_stack_caille.c

    r0bef47b r1c6286d  
    33*/ 
    44 
    5 double Iq(double qval, 
    6       double del, 
    7       double Nlayers,  
    8       double dd, 
    9           double Cp,  
    10       double sld, 
    11       double solvent_sld); 
    12  
    13 double Iq(double qval, 
    14       double del, 
    15       double Nlayers,  
    16       double dd, 
    17           double Cp,  
    18       double sld, 
    19       double solvent_sld) 
     5static double 
     6Iq(double qval, 
     7   double del, 
     8   double fp_Nlayers, 
     9   double dd, 
     10   double Cp, 
     11   double sld, 
     12   double solvent_sld) 
    2013{ 
    21   double contr,NN;   //local variables of coefficient wave 
     14  int Nlayers = (int)(fp_Nlayers+0.5);    //cast to an integer for the loop 
     15  double contr;   //local variables of coefficient wave 
    2216  double inten,Pq,Sq,alpha,temp,t2; 
    2317  //double dQ, dQDefault, t1, t3; 
    24   int ii,NNint; 
    2518  // from wikipedia 0.577215664901532860606512090082402431042159335 
    2619  const double Euler = 0.577215664901533;   // Euler's constant, increased sig figs for new models Feb 2015 
    … …  
    2821  //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015 
    2922 
    30   NN = trunc(Nlayers);    //be sure that NN is an integer 
    31    
    3223  contr = sld - solvent_sld; 
    3324 
    3425  Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)); 
    3526 
    36   NNint = (int)NN;    //cast to an integer for the loop 
    37   ii=0; 
    3827  Sq = 0.0; 
    39   // the vital "=" in ii<=  added March 2015 
    40   for(ii=1;ii<=(NNint-1);ii+=1) { 
    41  
    42     //fii = (double)ii;   //do I really need to do this? - unused variable, removed 18Feb2015 
    43  
     28  for (int ii=1; ii < Nlayers; ii++) { 
    4429    temp = 0.0; 
    4530    alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler); 
    … …  
    4833    //t3 = dQ*dQ*dd*dd*ii*ii; 
    4934 
    50     temp = 1.0-ii/NN; 
     35    temp = 1.0 - (double)ii / (double)Nlayers; 
    5136    //temp *= cos(dd*qval*ii/(1.0+t1)); 
    5237    temp *= cos(dd*qval*ii); 
    5338    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); 
    54     temp *= exp(-t2/2.0 ); 
     39    temp *= exp(-t2/2.0); 
    5540    //temp /= sqrt(1.0+t1); 
    5641 

  • sasmodels/models/lamellar_stack_caille.py

    r7c57861 ra57b31d  
    8888parameters = [ 
    8989    ["thickness",        "Ang",      30.0,  [0, inf],   "volume", "sheet thickness"], 
    90     ["Nlayers",          "",          20,   [0, inf],   "",       "Number of layers"], 
     90    ["Nlayers",          "",          20,   [1, inf],   "",       "Number of layers"], 
    9191    ["d_spacing",        "Ang",      400.,  [0.0,inf],  "volume", "lamellar d-spacing of Caille S(Q)"], 
    9292    ["Caille_parameter", "1/Ang^2",    0.1, [0.0,0.8],  "",       "Caille parameter"], 

  • sasmodels/models/lamellar_stack_paracrystal.c

    r4962519 r5467cd8  
    22 
    33*/ 
    4 double Iq(double qval, 
    5       double th, 
    6       double Nlayers,  
    7           double davg,  
    8           double pd, 
    9       double sld, 
    10       double solvent_sld); 
    11 double paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an); 
    12 double paraCryst_an(double ww, double qval, double davg, long Nlayers); 
     4double paraCryst_sn(double ww, double qval, double davg, int Nlayers, double an); 
     5double paraCryst_an(double ww, double qval, double davg, int Nlayers); 
    136 
    14 double Iq(double qval, 
    15       double th, 
    16       double Nlayers,  
    17           double davg,  
    18           double pd, 
    19       double sld, 
    20       double solvent_sld) 
     7static double 
     8Iq(double qval, 
     9   double th, 
     10   double fp_Nlayers, 
     11   double davg, 
     12   double pd, 
     13   double sld, 
     14   double solvent_sld) 
    2115{ 
    22      
    23         double inten,contr,xn; 
    24         double xi,ww,Pbil,Znq,Snq,an; 
    25         long n1,n2; 
     16        //get the fractional part of Nlayers, to determine the "mixing" of N's 
     17        int n1 = (int)(fp_Nlayers);             //truncate towards zero 
     18        int n2 = n1 + 1; 
     19        const double xn = (double)n2 - fp_Nlayers;      //fractional contribution of n1 
    2620         
    27         contr = sld - solvent_sld; 
    28         //get the fractional part of Nlayers, to determine the "mixing" of N's 
    29          
    30         n1 = (long)trunc(Nlayers);              //rounds towards zero 
    31         n2 = n1 + 1; 
    32         xn = (double)n2 - Nlayers;                      //fractional contribution of n1 
    33          
    34         ww = exp(-qval*qval*pd*pd*davg*davg/2.0); 
     21        const double ww = exp(-0.5*square(qval*pd*davg)); 
    3522 
    3623        //calculate the n1 contribution 
     24        double Znq,Snq,an; 
    3725        an = paraCryst_an(ww,qval,davg,n1); 
    3826        Snq = paraCryst_sn(ww,qval,davg,n1,an); 
    39          
     27 
    4028        Znq = xn*Snq; 
    4129         
    … …  
    5240//      Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg)) 
    5341         
    54         xi = th/2.0;            //use 1/2 the bilayer thickness 
    55         Pbil = (sin(qval*xi)/(qval*xi))*(sin(qval*xi)/(qval*xi)); 
     42        const double xi = th/2.0;               //use 1/2 the bilayer thickness 
     43        const double Pbil = square(sas_sinx_x(qval*xi)); 
    5644         
    57         inten = 2.0*M_PI*contr*contr*Pbil*Znq/(qval*qval); 
    58         inten *= 1.0e-04; 
     45        const double contr = sld - solvent_sld; 
     46        const double inten = 2.0*M_PI*contr*contr*Pbil*Znq/(qval*qval); 
    5947//printf("q=%.7e wwm1=%g ww=%.5e an=% 12.5e Snq=% 12.5e Znq=% 12.5e Pbil=% 12.5e\n",qval,wwm1,ww,an,Snq,Znq,Pbil); 
    60         return(inten); 
     48        return 1.0e-4*inten; 
    6149} 
    6250 
    6351// functions for the lamellar paracrystal model 
    6452double 
    65 paraCryst_sn(double ww, double qval, double davg, long Nlayers, double an) { 
     53paraCryst_sn(double ww, double qval, double davg, int Nlayers, double an) { 
    6654         
    6755        double Snq; 
    … …  
    6957        Snq = an/( (double)Nlayers*square(1.0+ww*ww-2.0*ww*cos(qval*davg)) ); 
    7058         
    71         return(Snq); 
     59        return Snq; 
    7260} 
    7361 
    7462double 
    75 paraCryst_an(double ww, double qval, double davg, long Nlayers) { 
    76          
     63paraCryst_an(double ww, double qval, double davg, int Nlayers) { 
    7764        double an; 
    7865         
    … …  
    8269        an += 2.0*pow(ww,(Nlayers+1))*cos((double)(Nlayers+1)*qval*davg); 
    8370         
    84         return(an); 
     71        return an; 
    8572} 
    8673 

  • sasmodels/models/lamellar_stack_paracrystal.py

    r7c57861 ra0168e8  
    113113parameters = [["thickness", "Ang", 33.0, [0, inf], "volume", 
    114114               "sheet thickness"], 
    115               ["Nlayers", "", 20, [0, inf], "", 
     115              ["Nlayers", "", 20, [1, inf], "", 
    116116               "Number of layers"], 
    117117              ["d_spacing", "Ang", 250., [0.0, inf], "", 

  • sasmodels/models/linear_pearls.c

    r925ad6e rc3ccaec  
    44            double radius, 
    55            double edge_sep, 
    6             double num_pearls, 
     6            double fp_num_pearls, 
    77            double pearl_sld, 
    88            double solvent_sld); 
    … …  
    1111            double radius, 
    1212            double edge_sep, 
    13             double num_pearls, 
     13            int num_pearls, 
    1414            double pearl_sld, 
    1515            double solvent_sld); 
    1616 
    1717 
    18 double form_volume(double radius, double num_pearls) 
     18double form_volume(double radius, double fp_num_pearls) 
    1919{ 
     20    int num_pearls = (int)(fp_num_pearls + 0.5); 
    2021    // Pearl volume 
    2122    double pearl_vol = M_4PI_3 * cube(radius); 
    … …  
    2728            double radius, 
    2829            double edge_sep, 
    29             double num_pearls, 
     30            int num_pearls, 
    3031            double pearl_sld, 
    3132            double solvent_sld) 
    3233{ 
    33     double n_contrib; 
    3434    //relative sld 
    3535    double contrast_pearl = pearl_sld - solvent_sld; 
    … …  
    4646    double psi = sas_3j1x_x(q * radius); 
    4747 
    48     // N pearls contribution 
    49     int n_max = num_pearls - 1; 
    50     n_contrib = num_pearls; 
    51     for(int num=1; num<=n_max; num++) { 
    52         n_contrib += (2.0*(num_pearls-num)*sas_sinx_x(q*separation*num)); 
     48    // N pearls interaction terms  
     49    double structure_factor = (double)num_pearls; 
     50    for(int num=1; num<num_pearls; num++) { 
     51        structure_factor += 2.0*(num_pearls-num)*sas_sinx_x(q*separation*num); 
    5352    } 
    5453    // form factor for num_pearls 
    55     double form_factor = 1.0e-4 * n_contrib * square(m_s*psi) / tot_vol; 
     54    double form_factor = 1.0e-4 * structure_factor * square(m_s*psi) / tot_vol; 
    5655 
    5756    return form_factor; 
    … …  
    6160            double radius, 
    6261            double edge_sep, 
    63             double num_pearls, 
     62            double fp_num_pearls, 
    6463            double pearl_sld, 
    6564            double solvent_sld) 
    6665{ 
    6766 
     67    int num_pearls = (int)(fp_num_pearls + 0.5); 
    6868        double result = linear_pearls_kernel(q, 
    6969                    radius, 

  • sasmodels/models/linear_pearls.py

    r925ad6e rc3ccaec  
    1616.. math:: 
    1717 
    18     P(Q) = \frac{scale}{V}\left[ m_{p}^2 
    19     \left(N+2\sum_{n-1}^{N-1}(N-n)\frac{sin(qnl)}{qnl}\right) 
    20     \left( 3\frac{sin(qR)-qRcos(qR)}{(qr)^3}\right)^2\right] 
     18    P(Q) = \frac{\text{scale}}{V}\left[ m_{p}^2 
     19    \left(N+2\sum_{n-1}^{N-1}(N-n)\frac{\sin(qnl)}{qnl}\right) 
     20    \left( 3\frac{\sin(qR)-qR\cos(qR)}{(qr)^3}\right)^2\right] 
    2121 
    2222where the mass $m_p$ is $(SLD_{pearl}-SLD_{solvent})*(volume\ of\ N\ pearls)$. 
    … …  
    5656    ["radius",      "Ang",       80.0, [0, inf],     "", "Radius of the pearls"], 
    5757    ["edge_sep",    "Ang",      350.0, [0, inf],     "", "Length of the string segment - surface to surface"], 
    58     ["num_pearls",  "",           3.0, [0, inf],     "", "Number of the pearls"], 
     58    ["num_pearls",  "",           3.0, [1, inf],     "", "Number of the pearls"], 
    5959    ["sld",   "1e-6/Ang^2", 1.0, [-inf, inf],  "sld", "SLD of the pearl spheres"], 
    6060    ["sld_solvent", "1e-6/Ang^2", 6.3, [-inf, inf],  "sld", "SLD of the solvent"], 

  • sasmodels/models/multilayer_vesicle.c

    r925ad6e rc3ccaec  
    77          double sld_solvent, 
    88          double sld, 
    9           double n_pairs) 
     9          int n_pairs) 
    1010{ 
    1111    //calculate with a loop, two shells at a time 
    … …  
    4747          double sld_solvent, 
    4848          double sld, 
    49           double n_pairs) 
     49          double fp_n_pairs) 
    5050{ 
     51    int n_pairs = (int)(fp_n_pairs + 0.5); 
    5152    return multilayer_vesicle_kernel(q, 
    5253           volfraction, 

  • sasmodels/models/multilayer_vesicle.py

    r925ad6e rc3ccaec  
    1919 
    2020.. math:: 
     21    P(q) = \text{scale} \cdot \frac{V_f}{V_t} F^2(q) + \text{background} 
    2122 
    22     P(q) = \frac{\text{scale.volfraction}}{V_t} F^2(q) + \text{background} 
    23  
    24 where 
     23for 
    2524 
    2625.. math:: 
     26    F(q) = (\rho_\text{shell}-\rho_\text{solv}) \sum_{i=1}^{n_\text{pairs}} 
     27        \left[ 
     28          3V(R_i)\frac{\sin(qR_i)-qR_i\cos(qR_i)}{(qR_i)^3} \\ 
     29          - 3V(R_i+t_s)\frac{\sin(q(R_i+t_s))-q(R_i+t_s)\cos(q(R_i+t_s))}{(q(R_i+t_s))^3} 
     30        \right] 
    2731 
    28      F(q) = (\rho_{shell}-\rho_{solv}) \sum_{i=1}^{n\_pairs} \left[ 
    29      3V(R_i)\frac{\sin(qR_i)-qR_i\cos(qR_i)}{(qR_i)^3} \\ 
    30       - 3V(R_i+t_s)\frac{\sin(q(R_i+t_s))-q(R_i+t_s)\cos(q(R_i+t_s))}{(q(R_i+t_s))^3} 
    31      \right] 
     32and 
    3233 
     34.. math:: 
     35     R_i = r_c + (i-1)(t_s + t_w) 
    3336 
    34 where $R_i = r_c + (i-1)(t_s + t_w)$ 
    35     
    36 where $V_t$ is the volume of the whole particle, $V(R)$ is the volume of a sphere 
    37 of radius $R$, $r_c$ is the radius of the core, $\rho_{shell}$ is the scattering length  
    38 density of a shell, $\rho_{solv}$ is the scattering length density of the solvent. 
     37where $V_f$ is the volume fraction of particles, $V_t$ is the volume of the 
     38whole particle, $V(r)$ is the volume of a sphere of radius $r$, $r_c$ is the 
     39radius of the core, $\rho_\text{shell}$ is the scattering length density of a 
     40shell, $\rho_\text{solv}$ is the scattering length density of the solvent. 
    3941 
     42The outer most radius, $r_o = R_n + t_s$, is used for both the volume fraction 
     43normalization and for the effective radius for *S(Q)* when $P(Q) * S(Q)$ 
     44is applied. 
    4045 
    4146The 2D scattering intensity is the same as 1D, regardless of the orientation 
    … …  
    4550 
    4651    q = \sqrt{q_x^2 + q_y^2} 
    47  
    48  
    49 The outer most radius 
    50  
    51 $radius + n\_pairs * thick\_shell + (n\_pairs- 1) * thick\_solvent$ 
    52  
    53 is used for both the volume fraction normalization and for the  
    54 effective radius for *S(Q)* when $P(Q) * S(Q)$ is applied. 
    5552 
    5653For information about polarised and magnetic scattering, see 
    … …  
    7067**Author:** NIST IGOR/DANSE **on:** pre 2010 
    7168 
    72 **Last Modified by:** Piotr Rozyczko**on:** Feb 24, 2016 
     69**Last Modified by:** Piotr Rozyczko **on:** Feb 24, 2016 
    7370 
    7471**Last Reviewed by:** Paul Butler **on:** March 20, 2016 
    … …  
    109106source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"] 
    110107 
     108# TODO: the following line does nothing 
    111109polydispersity = ["radius", "n_pairs"] 
    112110 

  • sasmodels/models/onion.c

    r925ad6e rc3ccaec  
    3030 
    3131static double 
    32 form_volume(double radius_core, double n, double thickness[]) 
     32form_volume(double radius_core, double n_shells, double thickness[]) 
    3333{ 
    34   int i; 
     34  int n = (int)(n_shells+0.5); 
    3535  double r = radius_core; 
    36   for (i=0; i < n; i++) { 
     36  for (int i=0; i < n; i++) { 
    3737    r += thickness[i]; 
    3838  } 

  • sasmodels/models/onion.py

    r925ad6e rc3ccaec  
    323323    Returns shape profile with x=radius, y=SLD. 
    324324    """ 
    325  
     325    n_shells = int(n_shells+0.5) 
    326326    total_radius = 1.25*(sum(thickness[:n_shells]) + radius_core + 1) 
    327327    dz = total_radius/400  # 400 points for a smooth plot 
    … …  
    366366    return np.asarray(z), np.asarray(rho) 
    367367 
    368 def ER(radius_core, n, thickness): 
     368def ER(radius_core, n_shells, thickness): 
    369369    """Effective radius""" 
    370     return np.sum(thickness[:int(n[0])], axis=0) + radius_core 
     370    n = int(n_shells[0]+0.5) 
     371    return np.sum(thickness[:n], axis=0) + radius_core 
    371372 
    372373demo = { 

  • sasmodels/models/pearl_necklace.c

    r4b541ac r3f853beb  
    11double form_volume(double radius, double edge_sep, 
    2     double thick_string, double num_pearls); 
     2    double thick_string, double fp_num_pearls); 
    33double Iq(double q, double radius, double edge_sep, 
    4     double thick_string, double num_pearls, double sld,  
     4    double thick_string, double fp_num_pearls, double sld, 
    55    double string_sld, double solvent_sld); 
    66 
    … …  
    99// From Igor library 
    1010static double 
    11 _pearl_necklace_kernel(double q, double radius, double edge_sep, double thick_string, 
    12     double num_pearls, double sld_pearl, double sld_string, double sld_solv) 
     11pearl_necklace_kernel(double q, double radius, double edge_sep, double thick_string, 
     12    int num_pearls, double sld_pearl, double sld_string, double sld_solv) 
    1313{ 
    1414    // number of string segments 
    15     num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 
    16     const double num_strings = num_pearls - 1.0; 
     15    const int num_strings = num_pearls - 1; 
    1716 
    1817    //each masses: contrast * volume 
    … …  
    6968 
    7069double form_volume(double radius, double edge_sep, 
    71     double thick_string, double num_pearls) 
     70    double thick_string, double fp_num_pearls) 
    7271{ 
    73     num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 
    74  
    75     const double num_strings = num_pearls - 1.0; 
     72    const int num_pearls = (int)(fp_num_pearls + 0.5); //Force integer number of pearls 
     73    const int num_strings = num_pearls - 1; 
    7674    const double string_vol = edge_sep * M_PI_4 * thick_string * thick_string; 
    7775    const double pearl_vol = M_4PI_3 * radius * radius * radius; 
    … …  
    8280 
    8381double Iq(double q, double radius, double edge_sep, 
    84     double thick_string, double num_pearls, double sld,  
     82    double thick_string, double fp_num_pearls, double sld, 
    8583    double string_sld, double solvent_sld) 
    8684{ 
    87     const double form = _pearl_necklace_kernel(q, radius, edge_sep, 
     85    const int num_pearls = (int)(fp_num_pearls + 0.5); //Force integer number of pearls 
     86    const double form = pearl_necklace_kernel(q, radius, edge_sep, 
    8887        thick_string, num_pearls, sld, string_sld, solvent_sld); 
    8988 

  • sasmodels/models/pearl_necklace.py

    r4b541ac r1bd1ea2  
    8282              ["thick_string", "Ang", 2.5, [0, inf], "volume", 
    8383               "Thickness of the chain linkage"], 
    84               ["num_pearls", "none", 3, [0, inf], "volume", 
     84              ["num_pearls", "none", 3, [1, inf], "volume", 
    8585               "Number of pearls in the necklace (must be integer)"], 
    8686              ["sld", "1e-6/Ang^2", 1.0, [-inf, inf], "sld", 
    … …  
    100100    Redundant with form_volume. 
    101101    """ 
     102    num_pearls = int(num_pearls + 0.5) 
    102103    number_of_strings = num_pearls - 1.0 
    103104    string_vol = edge_sep * pi * pow((thick_string / 2.0), 2.0) 
    … …  
    111112    Calculation for effective radius. 
    112113    """ 
     114    num_pearls = int(num_pearls + 0.5) 
    113115    tot_vol = volume(radius, edge_sep, thick_string, num_pearls) 
    114116    rad_out = (tot_vol/(4.0/3.0*pi)) ** (1./3.) 

  • sasmodels/models/raspberry.py

    r8e68ea0 rc3ccaec  
    1010    Schematic of the raspberry model 
    1111 
    12 In order to calculate the form factor of the entire complex, the self- 
    13 correlation of the large droplet, the self-correlation of the particles, the 
    14 correlation terms between different particles and the cross terms between large 
    15 droplet and small particles all need to be calculated. 
     12In order to calculate the form factor of the entire complex, the 
     13self-correlation of the large droplet, the self-correlation of the particles, 
     14the correlation terms between different particles and the cross terms between 
     15large droplet and small particles all need to be calculated. 
    1616 
    17 Consider two infinitely thin shells of radii R1 and R2 separated by distance r. 
    18 The general structure of the equation is then the form factor of the two shells 
    19 multiplied by the phase factor that accounts for the separation of their 
    20 centers. 
     17Consider two infinitely thin shells of radii $R_1$ and $R_2$ separated by 
     18distance $r$. The general structure of the equation is then the form factor 
     19of the two shells multiplied by the phase factor that accounts for the 
     20separation of their centers. 
    2121 
    2222.. math:: 

  • sasmodels/models/rpa.c

    racfb094 reb63414  
    1 double Iq(double q, double case_num, 
     1double Iq(double q, double fp_case_num, 
    22    double N[], double Phi[], double v[], double L[], double b[], 
    33    double Kab, double Kac, double Kad, 
    … …  
    55    ); 
    66 
    7 double Iq(double q, double case_num, 
     7double Iq(double q, double fp_case_num, 
    88    double N[],    // DEGREE OF POLYMERIZATION 
    99    double Phi[],  // VOL FRACTION 
    … …  
    1515    ) 
    1616{ 
    17   int icase = (int)case_num; 
     17  int icase = (int)(fp_case_num+0.5); 
    1818 
    1919  double Nab,Nac,Nad,Nbc,Nbd,Ncd; 
    … …  
    309309 
    310310  //calculate contrast where L[i] is the scattering length of i and D is the matrix 
    311   //need to verify why the sqrt of Nav rather than just Nav (assuming v is molar volume) 
     311  //Note that should multiply by Nav to get units of SLD which will become 
     312  // Nav*2 in the next line (SLD^2) but then normalization in that line would 
     313  //need to divide by Nav leaving only Nav or sqrt(Nav) before squaring.  
    312314  Nav=6.022045e+23; 
    313315  Lad=(L[0]/v[0]-L[3]/v[3])*sqrt(Nav); 
    … …  
    317319  Intg=Lad*Lad*S11+Lbd*Lbd*S22+Lcd*Lcd*S33+2.0*Lad*Lbd*S12+2.0*Lbd*Lcd*S23+2.0*Lad*Lcd*S13; 
    318320 
    319   //rescale for units of Lij^2 (in 10e-12 m^2 to m^2 ?) 
    320   Intg *= 1.0e-24;     
     321  //rescale for units of Lij^2 (fm^2 to cm^2) 
     322  Intg *= 1.0e-26;     
    321323 
    322324  return Intg; 

  • sasmodels/models/rpa.py

    rbb73096 reb63414  
    11r""" 
     2Definition 
     3---------- 
     4 
    25Calculates the macroscopic scattering intensity for a multi-component 
    36homogeneous mixture of polymers using the Random Phase Approximation. 
    … …  
    2427Case 9: A-B-C-D tetra-block copolymer 
    2528 
    26 **NB: these case numbers are different from those in the NIST SANS package!** 
     29.. note:: 
     30    These case numbers are different from those in the NIST SANS package! 
    2731 
    28 Only one case can be used at any one time. 
     32USAGE NOTES: 
    2933 
    30 The RPA (mean field) formalism only applies only when the multicomponent 
    31 polymer mixture is in the homogeneous mixed-phase region. 
    32  
    33 **Component D is assumed to be the "background" component (ie, all contrasts 
    34 are calculated with respect to component D).** So the scattering contrast 
    35 for a C/D blend = [SLD(component C) - SLD(component D)]\ :sup:`2`. 
    36  
    37 Depending on which case is being used, the number of fitting parameters - the 
    38 segment lengths (ba, bb, etc) and $\chi$ parameters (Kab, Kac, etc) - vary. 
    39 The *scale* parameter should be held equal to unity. 
    40  
    41 The input parameters are the degrees of polymerization, the volume fractions, 
    42 the specific volumes, and the neutron scattering length densities for each 
    43 component. 
     34* Only one case can be used at any one time. 
     35* The RPA (mean field) formalism only applies only when the multicomponent 
     36  polymer mixture is in the homogeneous mixed-phase region. 
     37* **Component D is assumed to be the "background" component (ie, all contrasts 
     38  are calculated with respect to component D).** So the scattering contrast 
     39  for a C/D blend = [SLD(component C) - SLD(component D)]\ :sup:`2`. 
     40* Depending on which case is being used, the number of fitting parameters can  
     41  vary.  Note that in general the degrees of polymerization, the volume 
     42  fractions, the molar volumes, and the neutron scattering lengths for each 
     43  component are obtained from other methods and held fixed while the segment 
     44  lengths (b\ :sub:`a`, b\ :sub:`b`, etc) and $\chi$ parameters (K\ :sub:`ab`, 
     45  K\ :sub:`ac`, etc). The *scale* parameter should be held equal to unity. 
    4446 
    4547 
    … …  
    4749---------- 
    4850 
    49 A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 4136 
     51.. [#] A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 4136 
    5052""" 
    5153 
    … …  
    5355 
    5456name = "rpa" 
    55 title = "Random Phase Approximation - unfinished work in progress" 
     57title = "Random Phase Approximation" 
    5658description = """ 
    5759This formalism applies to multicomponent polymer mixtures in the 
    … …  
    9092    ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
    9193    ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"], 
    92     ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
     94    ["v[4]", "mL/mol", 100.0, [0, inf], "", "molar volume"], 
    9395    ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"], 
    9496    ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"], 
    … …  
    114116    Return a list of parameters to hide depending on the multiplicity parameter. 
    115117    """ 
     118    case_num = int(case_num+0.5) 
    116119    if case_num < 2: 
    117120        return HIDE_AB 

  • sasmodels/models/spherical_sld.c

    r925ad6e rc3ccaec  
    11static double form_volume( 
    2     int n_shells, 
     2    double fp_n_shells, 
    33    double thickness[], 
    44    double interface[]) 
    55{ 
     6    int n_shells= (int)(fp_n_shells + 0.5); 
    67    double r = 0.0; 
    78    for (int i=0; i < n_shells; i++) { 
    … …  
    2021        return pow(z, nu); 
    2122    } else if (shape==2) { 
    22         return 1.0 - pow(1. - z, nu); 
     23        return 1.0 - pow(1.0 - z, nu); 
    2324    } else if (shape==3) { 
    2425        return expm1(-nu*z)/expm1(-nu); 
    … …  
    4445static double Iq( 
    4546    double q, 
    46     int n_shells, 
     47    double fp_n_shells, 
    4748    double sld_solvent, 
    4849    double sld[], 
    … …  
    5152    double shape[], 
    5253    double nu[], 
    53     int n_steps) 
     54    double fp_n_steps) 
    5455{ 
    5556    // iteration for # of shells + core + solvent 
     57    int n_shells = (int)(fp_n_shells + 0.5); 
     58    int n_steps = (int)(fp_n_steps + 0.5); 
    5659    double f=0.0; 
    5760    double r=0.0; 

  • sasmodels/models/spherical_sld.py

    r925ad6e rc3ccaec  
    233233    """ 
    234234 
     235    n_shells = int(n_shells + 0.5) 
     236    n_steps = int(n_steps + 0.5) 
    235237    z = [] 
    236238    rho = [] 
    … …  
    240242    rho.append(sld[0]) 
    241243 
    242     for i in range(0, int(n_shells)): 
     244    for i in range(0, n_shells): 
    243245        z_next += thickness[i] 
    244246        z.append(z_next) 
    … …  
    261263def ER(n_shells, thickness, interface): 
    262264    """Effective radius""" 
    263     n_shells = int(n_shells) 
     265    n_shells = int(n_shells + 0.5) 
    264266    total = (np.sum(thickness[:n_shells], axis=1) 
    265267             + np.sum(interface[:n_shells], axis=1)) 

  • sasmodels/models/stacked_disks.c

    r6c3e266 r19f996b  
    1 double form_volume(double thick_core, 
    2                    double thick_layer, 
    3                    double radius, 
    4                    double n_stacking); 
    5  
    6 double Iq(double q, 
    7           double thick_core, 
    8           double thick_layer, 
    9           double radius, 
    10           double n_stacking, 
    11           double sigma_dnn, 
    12           double core_sld, 
    13           double layer_sld, 
    14           double solvent_sld); 
    15  
    16 double Iqxy(double qx, double qy, 
    17           double thick_core, 
    18           double thick_layer, 
    19           double radius, 
    20           double n_stacking, 
    21           double sigma_dnn, 
    22           double core_sld, 
    23           double layer_sld, 
    24           double solvent_sld, 
    25           double theta, 
    26           double phi); 
    27  
    28 static 
    29 double _kernel(double q, 
    30                double radius, 
    31                double core_sld, 
    32                double layer_sld, 
    33                double solvent_sld, 
    34                double halfheight, 
    35                double thick_layer, 
    36                double sin_alpha, 
    37                double cos_alpha, 
    38                double sigma_dnn, 
    39                double d, 
    40                double n_stacking) 
     1static double stacked_disks_kernel( 
     2    double q, 
     3    double halfheight, 
     4    double thick_layer, 
     5    double radius, 
     6    int n_stacking, 
     7    double sigma_dnn, 
     8    double core_sld, 
     9    double layer_sld, 
     10    double solvent_sld, 
     11    double sin_alpha, 
     12    double cos_alpha, 
     13    double d) 
    4114 
    4215{ 
    … …  
    8861 
    8962 
    90 static 
    91 double stacked_disks_kernel(double q, 
    92                             double thick_core, 
    93                             double thick_layer, 
    94                             double radius, 
    95                             double n_stacking, 
    96                             double sigma_dnn, 
    97                             double core_sld, 
    98                             double layer_sld, 
    99                             double solvent_sld) 
     63static double stacked_disks_1d( 
     64    double q, 
     65    double thick_core, 
     66    double thick_layer, 
     67    double radius, 
     68    int n_stacking, 
     69    double sigma_dnn, 
     70    double core_sld, 
     71    double layer_sld, 
     72    double solvent_sld) 
    10073{ 
    10174/*    StackedDiscsX  :  calculates the form factor of a stacked "tactoid" of core shell disks 
    … …  
    11184        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    11285        SINCOS(zi, sin_alpha, cos_alpha); 
    113         double yyy = _kernel(q, 
     86        double yyy = stacked_disks_kernel(q, 
     87                           halfheight, 
     88                           thick_layer, 
    11489                           radius, 
     90                           n_stacking, 
     91                           sigma_dnn, 
    11592                           core_sld, 
    11693                           layer_sld, 
    11794                           solvent_sld, 
    118                            halfheight, 
    119                            thick_layer, 
    12095                           sin_alpha, 
    12196                           cos_alpha, 
    122                            sigma_dnn, 
    123                            d, 
    124                            n_stacking); 
     97                           d); 
    12598        summ += Gauss76Wt[i] * yyy * sin_alpha; 
    12699    } 
    … …  
    132105} 
    133106 
    134 double form_volume(double thick_core, 
    135                    double thick_layer, 
    136                    double radius, 
    137                    double n_stacking){ 
     107static double form_volume( 
     108    double thick_core, 
     109    double thick_layer, 
     110    double radius, 
     111    double fp_n_stacking) 
     112{ 
     113    int n_stacking = (int)(fp_n_stacking + 0.5); 
    138114    double d = 2.0 * thick_layer + thick_core; 
    139115    return M_PI * radius * radius * d * n_stacking; 
    140116} 
    141117 
    142 double Iq(double q, 
    143           double thick_core, 
    144           double thick_layer, 
    145           double radius, 
    146           double n_stacking, 
    147           double sigma_dnn, 
    148           double core_sld, 
    149           double layer_sld, 
    150           double solvent_sld) 
     118static double Iq( 
     119    double q, 
     120    double thick_core, 
     121    double thick_layer, 
     122    double radius, 
     123    double fp_n_stacking, 
     124    double sigma_dnn, 
     125    double core_sld, 
     126    double layer_sld, 
     127    double solvent_sld) 
    151128{ 
    152     return stacked_disks_kernel(q, 
     129    int n_stacking = (int)(fp_n_stacking + 0.5); 
     130    return stacked_disks_1d(q, 
    153131                    thick_core, 
    154132                    thick_layer, 
    … …  
    162140 
    163141 
    164 double 
    165 Iqxy(double qx, double qy, 
    166      double thick_core, 
    167      double thick_layer, 
    168      double radius, 
    169      double n_stacking, 
    170      double sigma_dnn, 
    171      double core_sld, 
    172      double layer_sld, 
    173      double solvent_sld, 
    174      double theta, 
    175      double phi) 
     142static double Iqxy(double qx, double qy, 
     143    double thick_core, 
     144    double thick_layer, 
     145    double radius, 
     146    double fp_n_stacking, 
     147    double sigma_dnn, 
     148    double core_sld, 
     149    double layer_sld, 
     150    double solvent_sld, 
     151    double theta, 
     152    double phi) 
    176153{ 
     154    int n_stacking = (int)(fp_n_stacking + 0.5); 
    177155    double q, sin_alpha, cos_alpha; 
    178156    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    … …  
    180158    double d = 2.0 * thick_layer + thick_core; 
    181159    double halfheight = 0.5*thick_core; 
    182     double answer = _kernel(q, 
     160    double answer = stacked_disks_kernel(q, 
     161                     halfheight, 
     162                     thick_layer, 
    183163                     radius, 
     164                     n_stacking, 
     165                     sigma_dnn, 
    184166                     core_sld, 
    185167                     layer_sld, 
    186168                     solvent_sld, 
    187                      halfheight, 
    188                      thick_layer, 
    189169                     sin_alpha, 
    190170                     cos_alpha, 
    191                      sigma_dnn, 
    192                      d, 
    193                      n_stacking); 
     171                     d); 
    194172 
    195173    //convert to [cm-1] 

  • sasmodels/models/stacked_disks.py

    rb7e8b94 rc3ccaec  
    126126    ["thick_layer", "Ang",        10.0, [0, inf],    "volume",      "Thickness of layer each side of core"], 
    127127    ["radius",      "Ang",        15.0, [0, inf],    "volume",      "Radius of the stacked disk"], 
    128     ["n_stacking",  "",            1.0, [0, inf],    "volume",      "Number of stacked layer/core/layer disks"], 
     128    ["n_stacking",  "",            1.0, [1, inf],    "volume",      "Number of stacked layer/core/layer disks"], 
    129129    ["sigma_d",     "Ang",         0,   [0, inf],    "",            "Sigma of nearest neighbor spacing"], 
    130130    ["sld_core",    "1e-6/Ang^2",  4,   [-inf, inf], "sld",         "Core scattering length density"], 

  • sasmodels/models/star_polymer.c

    r3a48772 r2586093f  
    33double Iq(double q, double radius2, double arms); 
    44 
    5 static double _mass_fractal_kernel(double q, double radius2, double arms) 
     5static double star_polymer_kernel(double q, double radius2, double arms) 
    66{ 
    77 
    … …  
    2323double Iq(double q, double radius2, double arms) 
    2424{ 
    25     return _mass_fractal_kernel(q, radius2, arms); 
     25    return star_polymer_kernel(q, radius2, arms); 
    2626} 

  • sasmodels/models/unified_power_Rg.py

    r66ca2a6 rc3ccaec  
    9797 
    9898def Iq(q, level, rg, power, B, G): 
    99     ilevel = int(level) 
    100     if ilevel == 0: 
     99    level = int(level + 0.5) 
     100    if level == 0: 
    101101        with errstate(divide='ignore'): 
    102102            return 1./q 
    … …  
    104104    with errstate(divide='ignore', invalid='ignore'): 
    105105        result = np.zeros(q.shape, 'd') 
    106         for i in range(ilevel): 
     106        for i in range(level): 
    107107            exp_now = exp(-(q*rg[i])**2/3.) 
    108108            pow_now = (erf(q*rg[i]/sqrt(6.))**3/q)**power[i] 
    109             if i < ilevel-1: 
     109            if i < level-1: 
    110110                exp_next = exp(-(q*rg[i+1])**2/3.) 
    111111            else: 
    … …  
    113113            result += G[i]*exp_now + B[i]*exp_next*pow_now 
    114114 
    115     result[q == 0] = np.sum(G[:ilevel]) 
     115    result[q == 0] = np.sum(G[:level]) 
    116116    return result 
    117117 

  • sasmodels/resolution.py

    rb397165 r7b7fcf0  
    8484        self.weight_matrix = pinhole_resolution(self.q_calc, self.q, 
    8585                                np.maximum(q_width, MINIMUM_RESOLUTION)) 
     86        self.q_calc = abs(self.q_calc) 
    8687 
    8788    def apply(self, theory): 
    … …  
    125126        self.weight_matrix = \ 
    126127            slit_resolution(self.q_calc, self.q, qx_width, qy_width) 
     128        self.q_calc = abs(self.q_calc) 
    127129 
    128130    def apply(self, theory): 
    … …  
    153155    # neither trapezoid nor Simpson's rule improved the accuracy. 
    154156    edges = bin_edges(q_calc) 
    155     edges[edges < 0.0] = 0.0 # clip edges below zero 
     157    #edges[edges < 0.0] = 0.0 # clip edges below zero 
    156158    cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :]) 
    157159    weights = cdf[1:] - cdf[:-1] 
    … …  
    286288    # The current algorithm is a midpoint rectangle rule. 
    287289    q_edges = bin_edges(q_calc) # Note: requires q > 0 
    288     q_edges[q_edges < 0.0] = 0.0 # clip edges below zero 
     290    #q_edges[q_edges < 0.0] = 0.0 # clip edges below zero 
    289291    weights = np.zeros((len(q), len(q_calc)), 'd') 
    290292 
Note: See TracChangeset for help on using the changeset viewer.