Changeset 0813a99 in sasmodels


Ignore:
Timestamp:
Mar 6, 2017 12:08:39 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:
d24e390
Parents:
7b7fcf0 (diff), acfb094 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into ticket-835

Files:
1 added
7 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/models/rpa.c

    r6351bfa racfb094  
    3939  double S31,S32,S33,S34,S41,S42,S43,S44; 
    4040  double Lad,Lbd,Lcd,Nav,Intg; 
    41  
     41   
     42  // Set values for non existent parameters (eg. no A or B in case 0 and 1 etc) 
    4243  //icase was shifted to N-1 from the original code 
    4344  if (icase <= 1){ 
    … …  
    5758    Kab = Kac = Kad = -0.0004; 
    5859  } 
    59  
     60  
     61  // Set volume fraction of component D based on constraint that sum of vol frac =1 
     62  Phi[3]=1.0-Phi[0]-Phi[1]-Phi[2]; 
     63 
     64  //set up values for cross terms in case of block copolymers (1,3,4,6,7,8,9) 
    6065  Nab=sqrt(N[0]*N[1]); 
    6166  Nac=sqrt(N[0]*N[2]); 
    … …  
    7984  Phicd=sqrt(Phi[2]*Phi[3]); 
    8085 
     86  // Calculate Q^2 * Rg^2 for each homopolymer assuming random walk  
    8187  Xa=q*q*b[0]*b[0]*N[0]/6.0; 
    8288  Xb=q*q*b[1]*b[1]*N[1]/6.0; 
    … …  
    8490  Xd=q*q*b[3]*b[3]*N[3]/6.0; 
    8591 
    86   Paa=2.0*(exp(-Xa)-1.0+Xa)/(Xa*Xa); 
    87   S0aa=N[0]*Phi[0]*v[0]*Paa; 
    88   Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb); 
     92  //calculate all partial structure factors Pij and normalize n^2 
     93  Paa=2.0*(exp(-Xa)-1.0+Xa)/(Xa*Xa); // free A chain form factor 
     94  S0aa=N[0]*Phi[0]*v[0]*Paa; // Phi * Vp * P(Q)= I(Q0)/delRho^2 
     95  Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb); //AB diblock (anchored Paa * anchored Pbb) partial form factor 
    8996  S0ab=(Phiab*vab*Nab)*Pab; 
    90   Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); 
     97  Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); //ABC triblock AC partial form factor 
    9198  S0ac=(Phiac*vac*Nac)*Pac; 
    92   Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); 
     99  Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); //ABCD four block 
    93100  S0ad=(Phiad*vad*Nad)*Pad; 
    94101 
    95102  S0ba=S0ab; 
    96   Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); 
     103  Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); // free B chain 
    97104  S0bb=N[1]*Phi[1]*v[1]*Pbb; 
    98   Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); 
     105  Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); // BC diblock 
    99106  S0bc=(Phibc*vbc*Nbc)*Pbc; 
    100   Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); 
     107  Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); // BCD triblock 
    101108  S0bd=(Phibd*vbd*Nbd)*Pbd; 
    102109 
    103110  S0ca=S0ac; 
    104111  S0cb=S0bc; 
    105   Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); 
     112  Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); // Free C chain 
    106113  S0cc=N[2]*Phi[2]*v[2]*Pcc; 
    107   Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); 
     114  Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); // CD diblock 
    108115  S0cd=(Phicd*vcd*Ncd)*Pcd; 
    109116 
    … …  
    111118  S0db=S0bd; 
    112119  S0dc=S0cd; 
    113   Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); 
     120  Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); // free D chain 
    114121  S0dd=N[3]*Phi[3]*v[3]*Pdd; 
     122 
     123  // Reset all unused partial structure factors to 0 (depends on case) 
    115124  //icase was shifted to N-1 from the original code 
    116125  switch(icase){ 
    … …  
    193202  S0dc=S0cd; 
    194203 
     204  // self chi parameter is 0 ... of course 
    195205  Kaa=0.0; 
    196206  Kbb=0.0; 
    … …  
    243253  ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd); 
    244254 
     255  // D is considered the matrix or background component so enters here 
    245256  m=1.0/(S0dd-ZZ); 
    246257 
    … …  
    297308  S44=S11+S22+S33+2.0*S12+2.0*S13+2.0*S23; 
    298309 
     310  //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) 
    299312  Nav=6.022045e+23; 
    300313  Lad=(L[0]/v[0]-L[3]/v[3])*sqrt(Nav); 
    … …  
    303316 
    304317  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; 
     318 
     319  //rescale for units of Lij^2 (in 10e-12 m^2 to m^2 ?) 
     320  Intg *= 1.0e-24;     
    305321 
    306322  return Intg; 

  • sasmodels/models/rpa.py

    r40a87fa rbb73096  
    107107 
    108108control = "case_num" 
    109 HIDE_NONE = set() 
    110 HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()) 
     109HIDE_ALL = set("Phi4".split()) 
     110HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()).union(HIDE_ALL) 
    111111HIDE_AB = set("N2 Phi2 v2 L2 b2 K23 K24".split()).union(HIDE_A) 
    112112def hidden(case_num): 
    … …  
    119119        return HIDE_A 
    120120    else: 
    121         return HIDE_NONE 
     121        return HIDE_ALL 
    122122 

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

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