Changeset 8115d82 in sasmodels


Ignore:
Timestamp:
Feb 10, 2016 5:53:39 AM (7 years ago)
Author:
wojciech
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
34d6cab
Parents:
a36c6d3
Message:

Two Power Law code converted from old model and ready for testing

Location:
sasmodels/models
Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/models/two_power_law.py

    ra36c6d3 r8115d82  
    11r""" 
    22This model calculates an empirical functional form for SAS data characterized 
    3 by two Lorentzian-type functions. 
     3by two power laws. 
    44 
    55Definition 
     
    1010.. math:: 
    1111 
    12     I(q) = \frac{A}{1 +(Q\xi_1)^n} + \frac{C}{1 +(Q\xi_2)^m} + \text{B} 
     12    I(q) = \begin{cases} 
     13    A \pow(qval,-m1) & q <= qc \\ 
     14    C \pow(qval,-m2) & q > qc 
     15    \end{cases} 
    1316 
    14 where $A$ = Lorentzian scale factor #1, $C$ = Lorentzian scale #2, 
    15 $\xi_1$ and $\xi_2$ are the corresponding correlation lengths, and $n$ and 
    16 $m$ are the respective power law exponents (set $n = m = 2$ for 
    17 Ornstein-Zernicke behaviour). 
     17where $qc$ = the location of the crossover from one slope to the other, 
     18$A$ = the scaling coefficent that sets the overall intensity of the lower Q power law region, 
     19$m1$ = power law exponent at low Q 
     20$m2$ = power law exponent at high Q 
     21 
     22The scaling of the second power law region (coefficent C) is then automatically scaled to 
     23match the first by following formula 
     24 
     25.. math:: 
     26    C = frac{A}{\pow(qc/-m1)\pow(qc/-m2)} 
     27 
     28...note 
     29    Be sure to enter the power law exponents as positive values! 
    1830 
    1931For 2D data the scattering intensity is calculated in the same way as 1D, 
     
    2537 
    2638 
    27 .. figure:: img/two_lorentzian.jpg 
     39.. figure:: img/two_power_law_1d.jpg 
    2840 
    2941    1D plot using the default values (w/500 data point). 
     
    3648""" 
    3749 
     50from numpy import power 
     51from numpy import sqrt 
    3852 
    39 from sas.models.BaseComponent import BaseComponent 
    40 from numpy import power 
    41 import math 
     53name = "two_power_law" 
     54title = "Two Power Law intensity calculation" 
     55description = """ 
     56            I(q) = coef_A*pow(qval,-1.0*power1) + background for q<=qc 
     57            =C*pow(qval,-1.0*power2) + background for q>qc 
     58            where C=coef_A*pow(qc,-1.0*power1)/pow(qc,-1.0*power2). 
    4259 
    43 class TwoPowerLawModel(BaseComponent): 
    44     """  
    45     Class that evaluates a TwoPowerLawModel. 
     60            coef_A = scaling coefficent 
     61            qc = crossover location [1/A] 
     62            power_1 (=m1) = power law exponent at low Q 
     63            power_2 (=m2) = power law exponent at high Q 
     64            background = Incoherent background [1/cm] 
     65        """ 
     66category = "shape-independent" 
    4667 
    47     Calculate:: 
     68# pylint: disable=bad-whitespace, line-too-long 
     69#            ["name", "units", default, [lower, upper], "type", "description"], 
     70parameters = [["coef_A",  "",     1.0, [-inf, inf], "", "scaling coefficent in low Q region"], 
     71              ["qc", "1/Ang", 0.04, [0, inf], "", "crossover location"], 
     72              ["power_1",    "",    1.0, [0, inf], "", "power law exponent at low Q"], 
     73              ["power_2",  "",      4.0, [0, inf], "", "power law exponent at high Q"], 
     74              ] 
     75# pylint: enable=bad-whitespace, line-too-long 
    4876 
    49        I(q) = coef_A pow(qval,-power1) for q<=qc 
    50        I(q) = C pow(qval,-power2) for q>qc 
     77def Iq(q, 
     78       coef_A=1.0, 
     79       qc=0.04, 
     80       power_1=1.0, 
     81       power_2=4.0, 
     82       ): 
    5183 
    52     where C=coef_A pow(qc,-power1)/pow(qc,-power2). 
    53      
    54     List of default parameters: 
    55      
    56     * coef_A = coefficient 
    57     * power1 = (-) Power @ low Q 
    58     * power2 = (-) Power @ high Q 
    59     * qc = crossover Q-value 
    60     * background = incoherent background 
    6184    """ 
    62          
    63     def __init__(self): 
    64         """ Initialization """ 
    65          
    66         # Initialize BaseComponent first, then sphere 
    67         BaseComponent.__init__(self) 
    68          
    69         ## Name of the model 
    70         self.name = "TwoPowerLaw" 
    71         self.description="""I(q) = coef_A*pow(qval,-1.0*power1) for q<=qc 
    72             =C*pow(qval,-1.0*power2) for q>qc 
    73             where C=coef_A*pow(qc,-1.0*power1)/pow(qc,-1.0*power2). 
    74              List of default parameters: 
    75              coef_A = coefficient 
    76              power1 = (-) Power @ low Q 
    77              power2 = (-) Power @ high Q 
    78              qc = crossover Q-value 
    79              background = incoherent background 
    80         """ 
    81         ## Define parameters 
    82         self.params = {} 
    83         self.params['coef_A']  = 1.0 
    84         self.params['power1']     = 1.0 
    85         self.params['power2']  = 4.0 
    86         self.params['qc']     = 0.04 
    87         self.params['background']     = 0.0 
    88         ## Parameter details [units, min, max] 
    89         self.details = {} 
    90         self.details['coef_A'] = ['', None, None] 
    91         self.details['power1'] =  ['', None, None] 
    92         self.details['power2']  =  ['', None, None] 
    93         self.details['qc']  =   ['1/A', None, None] 
    94         self.details['background']   =  ['[1/cm]', None, None] 
     85    :param q:                   Input q-value (float or [float, float]) 
     86    :param coef_A:              Scaling coefficent at low Q 
     87    :param qc:                  Crossover location 
     88    :param power_1:             Exponent of power law function at low Q 
     89    :param power_2:             Exponent of power law function at high Q 
     90    :return:                    Calculated intensity 
     91    """ 
    9592 
    96         #list of parameter that cannot be fitted 
    97         self.fixed= []   
    98     def _twopowerlaw(self, x): 
    99         """ 
    100         Model definition 
    101         """ 
    102         qc= self.params['qc'] 
    103         if(x<=qc): 
    104             inten = self.params['coef_A']*power(x,-1.0*self.params['power1']) 
    105         else: 
    106             scale = self.params['coef_A']*power(qc,-1.0*self.params['power1']) \ 
    107                                     / power(qc,-1.0*self.params['power2']) 
    108             inten = scale*power(x,-1.0*self.params['power2']) 
    109         inten += self.params['background'] 
     93# pylint: disable=bad-whitespace 
     94    if(g<=qc): 
     95        intensity = coef_A*power(q,-1.0*power_1) 
     96    else: 
     97        coef_C = coef_A*power(qc,-1.0*power_1)/power(qc,-1.0*power_2) 
     98        intensity = coef_C*power(q,-1.0*power_2) 
    11099 
    111         return inten   
    112     
    113     def run(self, x = 0.0): 
    114         """ Evaluate the model 
    115             @param x: input q-value (float or [float, float] as [r, theta]) 
    116             @return: (guinier value) 
    117         """ 
    118         if x.__class__.__name__ == 'list': 
    119             return self._twopowerlaw(x[0]) 
    120         elif x.__class__.__name__ == 'tuple': 
    121             raise ValueError, "Tuples are not allowed as input to BaseComponent models" 
    122         else: 
    123             return self._twopowerlaw(x) 
    124     
    125     def runXY(self, x = 0.0): 
    126         """ Evaluate the model 
    127             @param x: input q-value (float or [float, float] as [qx, qy]) 
    128             @return: guinier value 
    129         """ 
    130         if x.__class__.__name__ == 'list': 
    131             q = math.sqrt(x[0]**2 + x[1]**2) 
    132             return self._twopowerlaw(q) 
    133         elif x.__class__.__name__ == 'tuple': 
    134             raise ValueError, "Tuples are not allowed as input to BaseComponent models" 
    135         else: 
    136             return self._twopowerlaw(x) 
     100    return intensity 
     101 
     102Iq.vectorized = True  # Iq accepts an array of q values 
     103 
     104 
     105def Iqxy(qx, qy, *args): 
     106    """ 
     107    :param qx:   Input q_x-value 
     108    :param qy:   Input q_y-value 
     109    :param args: Remaining arguments 
     110    :return:     2D-Intensity 
     111    """ 
     112 
     113    return Iq(sqrt(qx**2 + qy**2), *args) 
     114 
     115Iqxy.vectorized = True  # Iqxy accepts an array of qx, qy values 
     116 
     117demo = dict(scale=1, background=0.1, 
     118            coef_A=1, 
     119            qc=0.04, 
     120            power_1=1.0, 
     121            power_2=4.0) 
     122 
     123oldname = "TwoPowerLawModel" 
     124oldpars = dict(background='background', 
     125                coef_A='coef_A', 
     126                qc='qc', 
     127                power_1='power1', 
     128                power_2='power2') 
     129 
     130tests = [ 
     131 
     132    # Accuracy tests based on content in test/utest_extra_models.py 
     133    [{'coeff_A':    1.0, 
     134      'qc':         0.04, 
     135      'power_1':    1.0, 
     136      'power_2':    4.0, 
     137      'background': 0.1, 
     138    }, 0.001, 1000], 
     139 
     140    [{'coeff_A':    1.0, 
     141      'qc':         0.04, 
     142      'power_1':    1.0, 
     143      'power_2':    4.0, 
     144      'background': 0.1, 
     145    }, 0.150141, 0.125945], 
     146 
     147    [{'coeff_A':    1.0, 
     148      'qc':         0.04, 
     149      'power_1':    1.0, 
     150      'power_2':    4.0, 
     151      'background': 0.1, 
     152    }, [0.442528,0.0], 0.00166884], 
     153] 
Note: See TracChangeset for help on using the changeset viewer.