Changeset 4903cfd in sasmodels


Ignore:
Timestamp:
Mar 3, 2017 1:46:45 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:
f700b59
Parents:
5467cd8 (diff), e1aa129 (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-815

Files:
1 added
29 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/ellipsoid.c

    r925ad6e r130d4c7  
    44    double radius_polar, double radius_equatorial, double theta, double phi); 
    55 
    6 double _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double sin_alpha); 
    7 double _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double sin_alpha) 
     6static double 
     7_ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double cos_alpha) 
    88{ 
    99    double ratio = radius_polar/radius_equatorial; 
    10     // Given the following under the radical: 
    11     //     1 + sin^2(T) (v^2 - 1) 
    12     // we can expand to match the form given in Guinier (1955) 
    13     //     = (1 - sin^2(T)) + v^2 sin^2(T) = cos^2(T) + sin^2(T) 
     10    // Using ratio v = Rp/Re, we can expand the following to match the 
     11    // form given in Guinier (1955) 
     12    //     r = Re * sqrt(1 + cos^2(T) (v^2 - 1)) 
     13    //       = Re * sqrt( (1 - cos^2(T)) + v^2 cos^2(T) ) 
     14    //       = Re * sqrt( sin^2(T) + v^2 cos^2(T) ) 
     15    //       = sqrt( Re^2 sin^2(T) + Rp^2 cos^2(T) ) 
     16    // 
    1417    // Instead of using pythagoras we could pass in sin and cos; this may be 
    1518    // slightly better for 2D which has already computed it, but it introduces 
    … …  
    1720    // leave it as is. 
    1821    const double r = radius_equatorial 
    19                      * sqrt(1.0 + sin_alpha*sin_alpha*(ratio*ratio - 1.0)); 
     22                     * sqrt(1.0 + cos_alpha*cos_alpha*(ratio*ratio - 1.0)); 
    2023    const double f = sas_3j1x_x(q*r); 
    2124 
    … …  
    3942    double total = 0.0; 
    4043    for (int i=0;i<76;i++) { 
    41         //const double sin_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 
    42         const double sin_alpha = Gauss76Z[i]*zm + zb; 
    43         total += Gauss76Wt[i] * _ellipsoid_kernel(q, radius_polar, radius_equatorial, sin_alpha); 
     44        //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 
     45        const double cos_alpha = Gauss76Z[i]*zm + zb; 
     46        total += Gauss76Wt[i] * _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 
    4447    } 
    4548    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    5962    double q, sin_alpha, cos_alpha; 
    6063    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    61     const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, sin_alpha); 
     64    const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 
    6265    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    6366 

  • sasmodels/models/hayter_msa.c

    r4962519 r3fc5d27  
    7070                SofQ=sqhcal(Qdiam, gMSAWave); 
    7171        }else{ 
    72         //SofQ=NaN; 
    73                 SofQ=-1.0; 
     72        SofQ=NAN; 
    7473                //      print "Error Level = ",ierr 
    7574                //      print "Please report HPMSA problem with above error code" 

  • sasmodels/sasview_model.py

    rfe8ff99 r749a7d4  
    1616import logging 
    1717from os.path import basename, splitext 
     18import thread 
    1819 
    1920import numpy as np  # type: ignore 
    … …  
    3839except ImportError: 
    3940    pass 
     41 
     42calculation_lock = thread.allocate_lock() 
    4043 
    4144SUPPORT_OLD_STYLE_PLUGINS = True 
    … …  
    605608        to the card for each evaluation. 
    606609        """ 
     610        ## uncomment the following when trying to debug the uncoordinated calls 
     611        ## to calculate_Iq 
     612        #if calculation_lock.locked(): 
     613        #    logging.info("calculation waiting for another thread to complete") 
     614        #    logging.info("\n".join(traceback.format_stack())) 
     615 
     616        with calculation_lock: 
     617            return self._calculate_Iq(qx, qy) 
     618 
     619    def _calculate_Iq(self, qx, qy=None): 
    607620        #core.HAVE_OPENCL = False 
    608621        if self._model is None: 
    … …  
    721734            return [self.params[par.name]], [1.0] 
    722735 
    723 def test_model(): 
     736def test_cylinder(): 
    724737    # type: () -> float 
    725738    """ 
    726     Test that a sasview model (cylinder) can be run. 
     739    Test that the cylinder model runs, returning the value at [0.1,0.1]. 
    727740    """ 
    728741    Cylinder = _make_standard_model('cylinder') 
    … …  
    733746    # type: () -> float 
    734747    """ 
    735     Test that a sasview model (cylinder) can be run. 
     748    Test that 2-D hardsphere model runs and doesn't produce NaN. 
    736749    """ 
    737750    Model = _make_standard_model('hardsphere') 
    … …  
    744757    # type: () -> float 
    745758    """ 
    746     Test that a sasview model (cylinder) can be run. 
     759    Test that the 2-D RPA model runs 
    747760    """ 
    748761    RPA = _make_standard_model('rpa') 
    … …  
    750763    return rpa.evalDistribution([0.1, 0.1]) 
    751764 
     765def test_empty_distribution(): 
     766    # type: () -> None 
     767    """ 
     768    Make sure that sasmodels returns NaN when there are no polydispersity points 
     769    """ 
     770    Cylinder = _make_standard_model('cylinder') 
     771    cylinder = Cylinder() 
     772    cylinder.setParam('radius', -1.0) 
     773    cylinder.setParam('background', 0.) 
     774    Iq = cylinder.evalDistribution(np.asarray([0.1])) 
     775    assert np.isnan(Iq[0]), "empty distribution fails" 
    752776 
    753777def test_model_list(): 
    754778    # type: () -> None 
    755779    """ 
    756     Make sure that all models build as sasview models. 
     780    Make sure that all models build as sasview models 
    757781    """ 
    758782    from .exception import annotate_exception 
    … …  
    781805 
    782806if __name__ == "__main__": 
    783     print("cylinder(0.1,0.1)=%g"%test_model()) 
     807    print("cylinder(0.1,0.1)=%g"%test_cylinder()) 
     808    #test_empty_distribution() 

  • sasmodels/sesans.py

    rb397165 r94d13f1  
    1414import numpy as np  # type: ignore 
    1515from numpy import pi, exp  # type: ignore 
    16 from scipy.special import jv as besselj 
    17 #import direct_model.DataMixin as model 
    18          
    19 def make_q(q_max, Rmax): 
    20     r""" 
    21     Return a $q$ vector suitable for SESANS covering from $2\pi/ (10 R_{\max})$ 
    22     to $q_max$. This is the integration range of the Hankel transform; bigger range and  
    23     more points makes a better numerical integration. 
    24     Smaller q_min will increase reliable spin echo length range.  
    25     Rmax is the "radius" of the largest expected object and can be set elsewhere. 
    26     q_max is determined by the acceptance angle of the SESANS instrument. 
     16from scipy.special import j0 
     17 
     18class SesansTransform(object): 
    2719    """ 
    28     from sas.sascalc.data_util.nxsunit import Converter 
     20    Spin-Echo SANS transform calculator.  Similar to a resolution function, 
     21    the SesansTransform object takes I(q) for the set of *q_calc* values and 
     22    produces a transformed dataset 
    2923 
    30     q_min = dq = 0.1 * 2*pi / Rmax 
    31     return np.arange(q_min, 
    32                      Converter(q_max[1])(q_max[0], 
    33                                          units="1/A"), 
    34                      dq) 
    35      
    36 def make_all_q(data): 
     24    *SElength* (A) is the set of spin-echo lengths in the measured data. 
     25 
     26    *zaccept* (1/A) is the maximum acceptance of scattering vector in the spin 
     27    echo encoding dimension (for ToF: Q of min(R) and max(lam)). 
     28 
     29    *Rmax* (A) is the maximum size sensitivity; larger radius requires more 
     30    computation time. 
    3731    """ 
    38     Return a $q$ vector suitable for calculating the total scattering cross section for 
    39     calculating the effect of finite acceptance angles on Time of Flight SESANS instruments. 
    40     If no acceptance is given, or unwanted (set "unwanted" flag in paramfile), no all_q vector is needed. 
    41     If the instrument has a rectangular acceptance, 2 all_q vectors are needed. 
    42     If the instrument has a circular acceptance, 1 all_q vector is needed 
    43      
    44     """ 
    45     if not data.has_no_finite_acceptance: 
    46         return [] 
    47     elif data.has_yz_acceptance(data): 
    48         # compute qx, qy 
    49         Qx, Qy = np.meshgrid(qx, qy) 
    50         return [Qx, Qy] 
    51     else: 
    52         # else only need q 
    53         # data.has_z_acceptance 
    54         return [q] 
     32    #: SElength from the data in the original data units; not used by transform 
     33    #: but the GUI uses it, so make sure that it is present. 
     34    q = None  # type: np.ndarray 
    5535 
    56 def transform(data, q_calc, Iq_calc, qmono, Iq_mono): 
    57     """ 
    58     Decides which transform type is to be used, based on the experiment data file contents (header) 
    59     (2016-03-19: currently controlled from parameters script) 
    60     nqmono is the number of q vectors to be used for the detector integration 
    61     """ 
    62     nqmono = len(qmono) 
    63     if nqmono == 0: 
    64         result = call_hankel(data, q_calc, Iq_calc) 
    65     elif nqmono == 1: 
    66         q = qmono[0] 
    67         result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono) 
    68     else: 
    69         Qx, Qy = [qmono[0], qmono[1]] 
    70         Qx = np.reshape(Qx, nqx, nqy) 
    71         Qy = np.reshape(Qy, nqx, nqy) 
    72         Iq_mono = np.reshape(Iq_mono, nqx, nqy) 
    73         qx = Qx[0, :] 
    74         qy = Qy[:, 0] 
    75         result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono) 
     36    #: q values to calculate when computing transform 
     37    q_calc = None  # type: np.ndarray 
    7638 
    77     return result 
     39    # transform arrays 
     40    _H = None  # type: np.ndarray 
     41    _H0 = None # type: np.ndarray 
    7842 
    79 def call_hankel(data, q_calc, Iq_calc): 
    80     return hankel((data.x, data.x_unit), 
    81                   (data.lam, data.lam_unit), 
    82                   (data.sample.thickness, 
    83                    data.sample.thickness_unit), 
    84                   q_calc, Iq_calc) 
    85    
    86 def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono): 
    87     return hankel(data.x, data.lam * 1e-9, 
    88                   data.sample.thickness / 10, 
    89                   q_calc, Iq_calc) 
    90                    
    91 def call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono): 
    92     return hankel(data.x, data.y, data.lam * 1e-9, 
    93                   data.sample.thickness / 10, 
    94                   q_calc, Iq_calc) 
    95                          
    96 def TotalScatter(model, parameters):  #Work in progress!! 
    97 #    Calls a model with existing model parameters already in place, then integrate the product of q and I(q) from 0 to (4*pi/lambda) 
    98     allq = np.linspace(0,4*pi/wavelength,1000) 
    99     allIq = 1 
    100     integral = allq*allIq 
    101      
     43    def __init__(self, z, SElength, zaccept, Rmax): 
     44        # type: (np.ndarray, float, float) -> None 
     45        #import logging; logging.info("creating SESANS transform") 
     46        self.q = z 
     47        self._set_hankel(SElength, zaccept, Rmax) 
    10248 
     49    def apply(self, Iq): 
     50        # tye: (np.ndarray) -> np.ndarray 
     51        G0 = np.dot(self._H0, Iq) 
     52        G = np.dot(self._H.T, Iq) 
     53        P = G - G0 
     54        return P 
    10355 
    104 def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still 
    105 #============================================================================== 
    106 #     2D Cosine Transform if "wavelength" is a vector 
    107 #============================================================================== 
    108 #allq is the q-space needed to create the total scattering cross-section 
     56    def _set_hankel(self, SElength, zaccept, Rmax): 
     57        # type: (np.ndarray, float, float) -> None 
     58        # Force float32 arrays, otherwise run into memory problems on some machines 
     59        SElength = np.asarray(SElength, dtype='float32') 
    10960 
    110     Gprime = np.zeros_like(wavelength, 'd') 
    111     s = np.zeros_like(wavelength, 'd') 
    112     sd = np.zeros_like(wavelength, 'd') 
    113     Gprime = np.zeros_like(wavelength, 'd') 
    114     f = np.zeros_like(wavelength, 'd') 
    115     for i, wavelength_i in enumerate(wavelength): 
    116         z = magfield*wavelength_i 
    117         allq=np.linspace() #for calculating the Q-range of the  scattering power integral 
    118         allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow 
    119         alldq = (allq[1]-allq[0])*1e10 
    120         sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 
    121         s[i]=1-exp(-sigma) 
    122         for j, Iqy_j, qy_j in enumerate(qy): 
    123             for k, Iqz_k, qz_k in enumerate(qz): 
    124                 Iq = np.sqrt(Iqy_j^2+Iqz_k^2) 
    125                 q = np.sqrt(qy_j^2 + qz_k^2) 
    126                 Gintegral = Iq*cos(z*Qz_k) 
    127                 Gprime[i] += Gintegral 
    128 #                sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i] 
    129 #                s[i] += 1-exp(Totalscatter(modelname)*thickness) 
    130 #                For now, work with standard 2-phase scatter 
     61        #Rmax = #value in text box somewhere in FitPage? 
     62        q_max = 2*pi / (SElength[1] - SElength[0]) 
     63        q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[-1]) 
     64        q = np.arange(q_min, q_max, q_min, dtype='float32') 
     65        dq = q_min 
    13166 
     67        H0 = np.float32(dq/(2*pi)) * q 
    13268 
    133                 sd[i] += Iq 
    134         f[i] = 1-s[i]+sd[i] 
    135         P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i] 
     69        repq = np.tile(q, (SElength.size, 1)).T 
     70        repSE = np.tile(SElength, (q.size, 1)) 
     71        H = np.float32(dq/(2*pi)) * j0(repSE*repq) * repq 
    13672 
    137  
    138  
    139  
    140 def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname): 
    141 #============================================================================== 
    142 #     HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS 
    143 #============================================================================== 
    144 #acceptq is the q-space needed to create limited acceptance effect 
    145     SElength= wavelength*magfield 
    146     G = np.zeros_like(SElength, 'd') 
    147     threshold=2*pi*theta/wavelength 
    148     for i, SElength_i in enumerate(SElength): 
    149         allq=np.linspace() #for calculating the Q-range of the  scattering power integral 
    150         allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow 
    151         alldq = (allq[1]-allq[0])*1e10 
    152         sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 
    153         s[i]=1-exp(-sigma) 
    154  
    155         dq = (q[1]-q[0])*1e10 
    156         a = (x<threshold) 
    157         acceptq = a*q 
    158         acceptIq = a*Iq 
    159  
    160         G[i] = np.sum(besselj(0, acceptq*SElength_i)*acceptIq*acceptq*dq) 
    161  
    162 #        G[i]=np.sum(integral) 
    163  
    164     G *= dq*1e10*2*pi 
    165  
    166     P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0])) 
    167      
    168 def hankel(SElength, wavelength, thickness, q, Iq): 
    169     r""" 
    170     Compute the expected SESANS polarization for a given SANS pattern. 
    171  
    172     Uses the hankel transform followed by the exponential.  The values for *zz* 
    173     (or spin echo length, or delta), wavelength and sample thickness should 
    174     come from the dataset.  $q$ should be chosen such that the oscillations 
    175     in $I(q)$ are well sampled (e.g., $5 \cdot 2 \pi/d_{\max}$). 
    176  
    177     *SElength* [A] is the set of $z$ points at which to compute the 
    178     Hankel transform 
    179  
    180     *wavelength* [m]  is the wavelength of each individual point *zz* 
    181  
    182     *thickness* [cm] is the sample thickness. 
    183  
    184     *q* [A$^{-1}$] is the set of $q$ points at which the model has been 
    185     computed. These should be equally spaced. 
    186  
    187     *I* [cm$^{-1}$] is the value of the SANS model at *q* 
    188     """ 
    189  
    190     from sas.sascalc.data_util.nxsunit import Converter 
    191     wavelength = Converter(wavelength[1])(wavelength[0],"A") 
    192     thickness = Converter(thickness[1])(thickness[0],"A") 
    193     Iq = Converter("1/cm")(Iq,"1/A") # All models default to inverse centimeters 
    194     SElength = Converter(SElength[1])(SElength[0],"A") 
    195  
    196     G = np.zeros_like(SElength, 'd') 
    197 #============================================================================== 
    198 #     Hankel Transform method if "wavelength" is a scalar; mono-chromatic SESANS 
    199 #============================================================================== 
    200     for i, SElength_i in enumerate(SElength): 
    201         integral = besselj(0, q*SElength_i)*Iq*q 
    202         G[i] = np.sum(integral) 
    203     G0 = np.sum(Iq*q) 
    204  
    205     # [m^-1] step size in q, needed for integration 
    206     dq = (q[1]-q[0]) 
    207  
    208     # integration step, convert q into [m**-1] and 2 pi circle integration 
    209     G *= dq*2*pi 
    210     G0 = np.sum(Iq*q)*dq*2*np.pi 
    211  
    212     P = exp(thickness*wavelength**2/(4*pi**2)*(G-G0)) 
    213  
    214     return P 
     73        self.q_calc = q 
     74        self._H, self._H0 = H, H0 

  • 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 rc3ccaec  
    107107    Returns the SLD profile *r* (Ang), and *rho* (1e-6/Ang^2). 
    108108    """ 
     109    n = int(n+0.5) 
    109110    z = [] 
    110111    rho = [] 

  • 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 rabb6e4f  
    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-1; 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); 
    … …  
    7154 
    7255  inten *= 1.0e-04;   // 1/A to 1/cm 
    73   return(inten); 
     56  return inten; 
    7457} 
    7558 

  • sasmodels/models/lamellar_stack_caille.c

    r0bef47b r56b0435  
    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; 
    3928  // the vital "=" in ii<=  added March 2015 
    40   for(ii=1;ii<=(NNint-1);ii+=1) { 
     29  for (int ii=1; ii<=Nlayers-1; ii++) { 
    4130 
    4231    //fii = (double)ii;   //do I really need to do this? - unused variable, removed 18Feb2015 
    … …  
    4837    //t3 = dQ*dQ*dd*dd*ii*ii; 
    4938 
    50     temp = 1.0-ii/NN; 
     39    temp = 1.0 - (double)ii / (double)Nlayers; 
    5140    //temp *= cos(dd*qval*ii/(1.0+t1)); 
    5241    temp *= cos(dd*qval*ii); 

  • 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 rc3ccaec  
    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 
    … …  
    6969 
    7070double form_volume(double radius, double edge_sep, 
    71     double thick_string, double num_pearls) 
     71    double thick_string, double fp_num_pearls) 
    7272{ 
    73     num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 
    74  
    75     const double num_strings = num_pearls - 1.0; 
     73    const int num_pearls = (int)(fp_num_pearls + 0.5); //Force integer number of pearls 
     74    const int num_strings = num_pearls - 1; 
    7675    const double string_vol = edge_sep * M_PI_4 * thick_string * thick_string; 
    7776    const double pearl_vol = M_4PI_3 * radius * radius * radius; 
    … …  
    8281 
    8382double Iq(double q, double radius, double edge_sep, 
    84     double thick_string, double num_pearls, double sld,  
     83    double thick_string, double fp_num_pearls, double sld, 
    8584    double string_sld, double solvent_sld) 
    8685{ 
    87     const double form = _pearl_necklace_kernel(q, radius, edge_sep, 
     86    const int num_pearls = (int)(fp_num_pearls + 0.5); //Force integer number of pearls 
     87    const double form = pearl_necklace_kernel(q, radius, edge_sep, 
    8888        thick_string, num_pearls, sld, string_sld, solvent_sld); 
    8989 

  • sasmodels/models/pearl_necklace.py

    r4b541ac rc3ccaec  
    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

    r6351bfa rfa4a994  
    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; 

  • sasmodels/models/rpa.py

    r40a87fa rfa4a994  
    114114    Return a list of parameters to hide depending on the multiplicity parameter. 
    115115    """ 
     116    case_num = int(case_num+0.5) 
    116117    if case_num < 2: 
    117118        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 rc3ccaec  
    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    double 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    double 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 
Note: See TracChangeset for help on using the changeset viewer.