Changeset 7e4a633 in sasmodels


Ignore:
Timestamp:
Mar 7, 2017 3:12:20 PM (8 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:
68f45cb
Parents:
f88e248 (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-843

Files:
31 edited

Legend:

Unmodified
Added
Removed

  • example/fit.py

    r1182da5 rf4b36fa  
    2424            % section) 
    2525data = radial_data if section != "tangential" else tan_data 
    26 phi = 0 if section != "tangential" else 90 
     26theta = 89.9 if section != "tangential" else 0 
     27phi = 90 
    2728kernel = load_model(name, dtype="single") 
    2829cutoff = 1e-3 
    … …  
    3031if name == "ellipsoid": 
    3132    model = Model(kernel, 
    32         scale=0.08, 
    33         r_polar=15, r_equatorial=800, 
     33        scale=0.08, background=35, 
     34        radius_polar=15, radius_equatorial=800, 
    3435        sld=.291, sld_solvent=7.105, 
    35         background=0, 
    36         theta=90, phi=phi, 
    37         theta_pd=15, theta_pd_n=40, theta_pd_nsigma=3, 
    38         r_polar_pd=0.222296, r_polar_pd_n=1, r_polar_pd_nsigma=0, 
    39         r_equatorial_pd=.000128, r_equatorial_pd_n=1, r_equatorial_pd_nsigma=0, 
     36        theta=theta, phi=phi, 
     37        theta_pd=0, theta_pd_n=0, theta_pd_nsigma=3, 
    4038        phi_pd=0, phi_pd_n=20, phi_pd_nsigma=3, 
     39        radius_polar_pd=0.222296, radius_polar_pd_n=1, radius_polar_pd_nsigma=0, 
     40        radius_equatorial_pd=.000128, radius_equatorial_pd_n=1, radius_equatorial_pd_nsigma=0, 
    4141        ) 
    4242 
    43  
    4443    # SET THE FITTING PARAMETERS 
    45     model.r_polar.range(15, 1000) 
    46     model.r_equatorial.range(15, 1000) 
    47     model.theta_pd.range(0, 360) 
     44    model.radius_polar.range(15, 1000) 
     45    model.radius_equatorial.range(15, 1000) 
     46    #model.theta.range(0, 90) 
     47    #model.theta_pd.range(0,10) 
     48    model.phi_pd.range(0,20) 
     49    model.phi.range(0, 180) 
    4850    model.background.range(0,1000) 
    4951    model.scale.range(0, 10) 
    5052 
    5153 
    52  
    5354elif name == "lamellar": 
    5455    model = Model(kernel, 
    55         scale=0.08, 
     56        scale=0.08, background=0.003, 
    5657        thickness=19.2946, 
    5758        sld=5.38,sld_sol=7.105, 
    58         background=0.003, 
    5959        thickness_pd= 0.37765, thickness_pd_n=10, thickness_pd_nsigma=3, 
    6060        ) 
    61  
    6261 
    6362    # SET THE FITTING PARAMETERS 
    … …  
    7776        radius_pd=.0084, radius_pd_n=10, radius_pd_nsigma=3, 
    7877        length_pd=0.493, length_pd_n=10, length_pd_nsigma=3, 
    79         phi_pd=0, phi_pd_n=5, phi_pd_nsigma=3,) 
     78        phi_pd=0, phi_pd_n=5 phi_pd_nsigma=3,) 
    8079        """ 
    8180    pars = dict( 
    8281        scale=.01, background=35, 
    8382        sld=.291, sld_solvent=5.77, 
    84         radius=250, length=178,  
    85         theta=90, phi=phi, 
     83        radius=250, length=178, 
    8684        radius_pd=0.1, radius_pd_n=5, radius_pd_nsigma=3, 
    8785        length_pd=0.1,length_pd_n=5, length_pd_nsigma=3, 
    88         theta_pd=10, theta_pd_n=50, theta_pd_nsigma=3, 
    89         phi_pd=0, phi_pd_n=10, phi_pd_nsigma=3) 
     86        theta=theta, phi=phi, 
     87        theta_pd=0, theta_pd_n=0, theta_pd_nsigma=3, 
     88        phi_pd=10, phi_pd_n=20, phi_pd_nsigma=3) 
    9089    model = Model(kernel, **pars) 
    9190 
    … …  
    9392    model.radius.range(1, 500) 
    9493    model.length.range(1, 5000) 
    95     model.theta.range(-90,100) 
    96     model.theta_pd.range(0, 30) 
    97     model.theta_pd_n = model.theta_pd + 5 
     94    #model.theta.range(0, 90) 
     95    model.phi.range(0, 180) 
     96    model.phi_pd.range(0, 30) 
    9897    model.radius_pd.range(0, 1) 
    99     model.length_pd.range(0, 2) 
     98    model.length_pd.range(0, 1) 
    10099    model.scale.range(0, 10) 
    101100    model.background.range(0, 100) 
    … …  
    104103elif name == "core_shell_cylinder": 
    105104    model = Model(kernel, 
    106         scale= .031, radius=19.5, thickness=30, length=22, 
    107         sld_core=7.105, sld_shell=.291, sdl_solvent=7.105, 
    108         background=0, theta=0, phi=phi, 
    109  
     105        scale= .031, background=0, 
     106        radius=19.5, thickness=30, length=22, 
     107        sld_core=7.105, sld_shell=.291, sld_solvent=7.105, 
    110108        radius_pd=0.26, radius_pd_n=10, radius_pd_nsigma=3, 
    111109        length_pd=0.26, length_pd_n=10, length_pd_nsigma=3, 
    112110        thickness_pd=1, thickness_pd_n=1, thickness_pd_nsigma=1, 
    113         theta_pd=1, theta_pd_n=10, theta_pd_nsigma=3, 
    114         phi_pd=0.1, phi_pd_n=1, phi_pd_nsigma=1, 
     111        theta=theta, phi=phi, 
     112        theta_pd=1, theta_pd_n=1, theta_pd_nsigma=3, 
     113        phi_pd=0, phi_pd_n=20, phi_pd_nsigma=3, 
    115114        ) 
    116115 
    117116    # SET THE FITTING PARAMETERS 
    118     #model.radius.range(115, 1000) 
    119     #model.length.range(0, 2500) 
     117    model.radius.range(115, 1000) 
     118    model.length.range(0, 2500) 
    120119    #model.thickness.range(18, 38) 
    121120    #model.thickness_pd.range(0, 1) 
    122121    #model.phi.range(0, 90) 
     122    model.phi_pd.range(0,20) 
    123123    #model.radius_pd.range(0, 1) 
    124124    #model.length_pd.range(0, 1) 
    … …  
    131131elif name == "capped_cylinder": 
    132132    model = Model(kernel, 
    133         scale=.08, radius=20, cap_radius=40, length=400, 
     133        scale=.08, background=35, 
     134        radius=20, cap_radius=40, length=400, 
    134135        sld=1, sld_solvent=6.3, 
    135         background=0, theta=0, phi=phi, 
    136136        radius_pd=.1, radius_pd_n=5, radius_pd_nsigma=3, 
    137137        cap_radius_pd=.1, cap_radius_pd_n=5, cap_radius_pd_nsigma=3, 
    138138        length_pd=.1, length_pd_n=1, length_pd_nsigma=0, 
    139         theta_pd=.1, theta_pd_n=1, theta_pd_nsigma=0, 
    140         phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0, 
     139        theta=theta, phi=phi, 
     140        theta_pd=0, theta_pd_n=1, theta_pd_nsigma=0, 
     141        phi_pd=10, phi_pd_n=20, phi_pd_nsigma=0, 
    141142        ) 
    142143 
     144    model.radius.range(115, 1000) 
     145    model.length.range(0, 2500) 
     146    #model.thickness.range(18, 38) 
     147    #model.thickness_pd.range(0, 1) 
     148    #model.phi.range(0, 90) 
     149    model.phi_pd.range(0,20) 
     150    #model.radius_pd.range(0, 1) 
     151    #model.length_pd.range(0, 1) 
     152    #model.theta_pd.range(0, 360) 
     153    #model.background.range(0,5) 
    143154    model.scale.range(0, 1) 
    144155 
    … …  
    146157elif name == "triaxial_ellipsoid": 
    147158    model = Model(kernel, 
    148         scale=0.08, req_minor=15, req_major=20, rpolar=500, 
     159        scale=0.08, background=35, 
     160        radius_equat_minor=15, radius_equat_major=20, radius_polar=500, 
    149161        sld=7.105, solvent_sld=.291, 
    150         background=5, theta=0, phi=phi, psi=0, 
     162        radius_equat_minor_pd=.1, radius_equat_minor_pd_n=1, radius_equat_minor_pd_nsigma=0, 
     163        radius_equat_major_pd=.1, radius_equat_major_pd_n=1, radius_equat_major_pd_nsigma=0, 
     164        radius_polar_pd=.1, radius_polar_pd_n=1, radius_polar_pd_nsigma=0, 
     165        theta=theta, phi=phi, psi=0, 
    151166        theta_pd=20, theta_pd_n=40, theta_pd_nsigma=3, 
    152167        phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0, 
    153168        psi_pd=30, psi_pd_n=1, psi_pd_nsigma=0, 
    154         req_minor_pd=.1, req_minor_pd_n=1, req_minor_pd_nsigma=0, 
    155         req_major_pd=.1, req_major_pd_n=1, req_major_pd_nsigma=0, 
    156         rpolar_pd=.1, rpolar_pd_n=1, rpolar_pd_nsigma=0, 
    157169        ) 
    158170 
    159171    # SET THE FITTING PARAMETERS 
    160     model.req_minor.range(15, 1000) 
    161     model.req_major.range(15, 1000) 
    162     #model.rpolar.range(15, 1000) 
     172    model.radius_equat_minor.range(15, 1000) 
     173    model.radius_equat_major.range(15, 1000) 
     174    #model.radius_polar.range(15, 1000) 
    163175    #model.background.range(0,1000) 
    164176    #model.theta_pd.range(0, 360) 
    … …  
    173185M = Experiment(data=data, model=model) 
    174186if section == "both": 
    175    tan_model = Model(model.kernel, **model.parameters()) 
     187   tan_model = Model(model.sasmodel, **model.parameters()) 
    176188   tan_model.phi = model.phi - 90 
    177189   tan_model.cutoff = cutoff 

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

    rec1d4bc r7e4a633  
    1919 
    2020.. math:: 
    21  
    22     P(q) = \text{scale}\frac{\phi}{V(R_N)} F^2(q) + \text{background} 
     21    P(q) = \text{scale} \cdot \frac{\phi}{V(R_N)} F^2(q) + \text{background} 
    2322 
    2423where 
    2524 
    2625.. math:: 
    27  
    2826     F(q) = (\rho_\text{shell}-\rho_\text{solv}) \sum_{i=1}^{N} \left[ 
    2927     3V(r_i)\frac{\sin(qr_i) - qr_i\cos(qr_i)}{(qr_i)^3} 
    … …  
    4442$\rho_\text{solv}$ is the scattering length density of the solvent. 
    4543 
     44The outer-most shell radius $R_N$ is used as the effective radius 
     45for $P(Q)$ when $P(Q) * S(Q)$ is applied. 
     46 
     47 
    4648The 2D scattering intensity is the same as 1D, regardless of the orientation 
    4749of the q vector which is defined as: 
    … …  
    5052 
    5153    q = \sqrt{q_x^2 + q_y^2} 
    52  
    53 The outer-most shell radius $R_N$ is used as the effective radius 
    54 for $P(Q)$ when $P(Q) * S(Q)$ is applied. 
    5554 
    5655For information about polarised and magnetic scattering, see 
    … …  
    7069**Author:** NIST IGOR/DANSE **on:** pre 2010 
    7170 
    72 **Last Modified by:** Piotr Rozyczko**on:** Feb 24, 2016 
     71**Last Modified by:** Piotr Rozyczko **on:** Feb 24, 2016 
    7372 
    7473**Last Reviewed by:** Paul Butler **on:** March 20, 2016 
    … …  
    106105    ] 
    107106# pylint: enable=bad-whitespace, line-too-long 
     107 
     108# TODO: proposed syntax for specifying which parameters can be polydisperse 
     109#polydispersity = ["radius", "thick_shell"] 
    108110 
    109111source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"] 

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

    r790b16bb r0c2da4b  
    549549            "radius": "core_radius", 
    550550            "sld_solvent": "core_sld", 
    551             "n_pairs": "n_pairs", 
     551            "n_shells": "n_pairs", 
    552552            "thick_shell": "s_thickness", 
    553553            "sld": "shell_sld", 

  • sasmodels/modelinfo.py

    r85fe7f8 rf88e248  
    230230    defined as a sublist with the following elements: 
    231231 
    232     *name* is the name that will be used in the call to the kernel 
    233     function and the name that will be displayed to the user.  Names 
     232    *name* is the name that will be displayed to the user.  Names 
    234233    should be lower case, with words separated by underscore.  If 
    235     acronyms are used, the whole acronym should be upper case. 
     234    acronyms are used, the whole acronym should be upper case. For vector 
     235    parameters, the name will be followed by *[len]* where *len* is an 
     236    integer length of the vector, or the name of the parameter which 
     237    controls the length.  The attribute *id* will be created from name 
     238    without the length. 
    236239 
    237240    *units* should be one of *degrees* for angles, *Ang* for lengths, 
    … …  
    603606        # Using the call_parameters table, we already have expanded forms 
    604607        # for each of the vector parameters; put them in a lookup table 
    605         expanded_pars = dict((p.name, p) for p in self.call_parameters) 
     608        # Note: p.id and p.name are currently identical for the call parameters 
     609        expanded_pars = dict((p.id, p) for p in self.call_parameters) 
    606610 
    607611        def append_group(name): 

  • sasmodels/models/multilayer_vesicle.c

    rc3ccaec rec1d4bc  
    1 static 
    2 double multilayer_vesicle_kernel(double q, 
     1static double 
     2form_volume(double radius, 
     3          double thick_shell, 
     4          double thick_solvent, 
     5          double fp_n_shells) 
     6{ 
     7    int n_shells = (int)(fp_n_shells + 0.5); 
     8    double R_N = radius + n_shells*(thick_shell+thick_solvent) - thick_solvent; 
     9    return M_4PI_3*cube(R_N); 
     10} 
     11 
     12static double 
     13multilayer_vesicle_kernel(double q, 
    314          double volfraction, 
    415          double radius, 
    … …  
    718          double sld_solvent, 
    819          double sld, 
    9           int n_pairs) 
     20          int n_shells) 
    1021{ 
    1122    //calculate with a loop, two shells at a time 
    … …  
    2940 
    3041        //do 2 layers at a time 
    31         ii += 1; 
     42        ii++; 
    3243 
    33     } while(ii <= n_pairs-1);  //change to make 0 < n_pairs < 2 correspond to 
     44    } while(ii <= n_shells-1);  //change to make 0 < n_shells < 2 correspond to 
    3445                               //unilamellar vesicles (C. Glinka, 11/24/03) 
    3546 
    36     fval *= volfraction*1.0e-4*fval/voli; 
    37  
    38     return(fval); 
     47    return 1.0e-4*volfraction*fval*fval;  // Volume normalization happens in caller 
    3948} 
    4049 
    41 static 
    42 double Iq(double q, 
     50static double 
     51Iq(double q, 
    4352          double volfraction, 
    4453          double radius, 
    … …  
    4756          double sld_solvent, 
    4857          double sld, 
    49           double fp_n_pairs) 
     58          double fp_n_shells) 
    5059{ 
    51     int n_pairs = (int)(fp_n_pairs + 0.5); 
     60    int n_shells = (int)(fp_n_shells + 0.5); 
    5261    return multilayer_vesicle_kernel(q, 
    5362           volfraction, 
    … …  
    5766           sld_solvent, 
    5867           sld, 
    59            n_pairs); 
     68           n_shells); 
    6069} 
    6170 

  • sasmodels/product.py

    r9951a86 rf88e248  
    4545    # structure factor calculator.  Structure factors should not 
    4646    # have any magnetic parameters 
    47     assert(s_info.parameters.kernel_parameters[0].id == ER_ID) 
    48     assert(s_info.parameters.kernel_parameters[1].id == VF_ID) 
    49     assert(s_info.parameters.magnetism_index == []) 
     47    if not s_info.parameters.kernel_parameters[0].id == ER_ID: 
     48        raise TypeError("S needs %s as first parameter"%ER_ID) 
     49    if not s_info.parameters.kernel_parameters[1].id == VF_ID: 
     50        raise TypeError("S needs %s as second parameter"%VF_ID) 
     51    if not s_info.parameters.magnetism_index == []: 
     52        raise TypeError("S should not have SLD parameters") 
    5053    p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters 
    5154    s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters 
    52     p_set = set(p.id for p in p_pars.call_parameters) 
    53     s_set = set(p.id for p in s_pars.call_parameters) 
    54  
    55     if p_set & s_set: 
    56         # there is some overlap between the parameter names; tag the 
    57         # overlapping S parameters with name_S. 
    58         # Skip the first parameter of s, which is effective radius 
    59         s_list = [(suffix_parameter(par) if par.id in p_set else par) 
    60                   for par in s_pars.kernel_parameters[1:]] 
    61     else: 
    62         # Skip the first parameter of s, which is effective radius 
    63         s_list = s_pars.kernel_parameters[1:] 
     55 
     56    # Create list of parameters for the combined model.  Skip the first 
     57    # parameter of S, which we verified above is effective radius.  If there 
     58    # are any names in P that overlap with those in S, modify the name in S 
     59    # to distinguish it. 
     60    p_set = set(p.id for p in p_pars.kernel_parameters) 
     61    s_list = [(_tag_parameter(par) if par.id in p_set else par) 
     62              for par in s_pars.kernel_parameters[1:]] 
     63    # Check if still a collision after renaming.  This could happen if for 
     64    # example S has volfrac and P has both volfrac and volfrac_S. 
     65    if any(p.id in p_set for p in s_list): 
     66        raise TypeError("name collision: P has P.name and P.name_S while S has S.name") 
     67 
    6468    translate_name = dict((old.id, new.id) for old, new 
    6569                          in zip(s_pars.kernel_parameters[1:], s_list)) 
    6670    demo = {} 
    67     demo.update((k, v) for k, v in p_info.demo.items() 
    68                 if k not in ("background", "scale")) 
     71    demo.update(p_info.demo.items()) 
    6972    demo.update((translate_name[k], v) for k, v in s_info.demo.items() 
    7073                if k not in ("background", "scale") and not k.startswith(ER_ID)) 
    … …  
    9093    # Remember the component info blocks so we can build the model 
    9194    model_info.composition = ('product', [p_info, s_info]) 
    92     model_info.demo = {} 
     95    model_info.demo = demo 
     96 
     97    ## Show the parameter table with the demo values 
     98    #from .compare import get_pars, parlist 
     99    #print("==== %s ====="%model_info.name) 
     100    #values = get_pars(model_info, use_demo=True) 
     101    #print(parlist(model_info, values, is2d=True)) 
    93102    return model_info 
    94103 
    95 def suffix_parameter(par, suffix): 
     104def _tag_parameter(par): 
     105    """ 
     106    Tag the parameter name with _S to indicate that the parameter comes from 
     107    the structure factor parameter set.  This is only necessary if the 
     108    form factor model includes a parameter of the same name as a parameter 
     109    in the structure factor. 
     110    """ 
    96111    par = copy(par) 
    97     par.name = par.name + " S" 
     112    # Protect against a vector parameter in S by appending the vector length 
     113    # to the renamed parameter.  Note: haven't tested this since no existing 
     114    # structure factor models contain vector parameters. 
     115    vector_length = par.name[len(par.id):] 
    98116    par.id = par.id + "_S" 
     117    par.name = par.id + vector_length 
    99118    return par 
    100119 
Note: See TracChangeset for help on using the changeset viewer.