source: sasview/sansmodels/src/sans/models/FractalModel.py @ 96b59384

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 96b59384 was 1ed3834, checked in by Jae Cho <jhjcho@…>, 16 years ago

unit symbol Angst. back to A and re-organized model descriptions
  • Property mode set to 100644
File size: 6.9 KB
RevLine 
[829eee9]1#!/usr/bin/env python
2"""
3   
4    Provide F(x)= P(x)*S(x) + bkd
5    Fractal as a BaseComponent model
6"""
7
8from sans.models.BaseComponent import BaseComponent
9import math
10from scipy.special import gamma
11
12class FractalModel(BaseComponent):
13   
14    """
15        Class that evaluates a Fractal function.
16       
17        F(x)= P(x)*S(x) + bkd
18        The model has Seven parameters:
[3db3895]19            scale        =  Volume fraction
20            radius       =  Block radius
21            fractal_dim  =  Fractal dimension
22            corr_length  =  correlation Length
23            block_sld    =  SDL block
24            solvent_sld  =  SDL solvent
25            background   =  background
[829eee9]26           
27    """
28       
29    def __init__(self):
30        """ Initialization """
31       
32        # Initialize BaseComponent first, then sphere
33        BaseComponent.__init__(self)
34       
35        ## Name of the model
36        self.name = "Number Density Fractal"
[96672c0]37        self.description="""
38        I(x)= P(x)*S(x) + bkd
[1ed3834]39       
40        p(x)= scale* V^(2)*delta^(2)* F(x*radius)^(2)
41        F(x) = 3*[sin(x)-x cos(x)]/x**3
42       
[96672c0]43        The model has Seven parameters:
44        scale        =  Volume fraction
45        radius       =  Block radius
46        fractal_dim  =  Fractal dimension
47        corr_length  =  correlation Length
48        block_sld    =  SDL block
49        solvent_sld  =  SDL solvent
50        background   =  background
51        """
[829eee9]52        ## Define parameters
53        self.params = {}
[3db3895]54        self.params['scale']       = 0.05
55        self.params['radius']      = 5.0
56        self.params['fractal_dim'] = 2.0
57        self.params['corr_length'] = 100.0
58        self.params['block_sld']   = 2.0e-6
59        self.params['solvent_sld'] = 6.0e-6
60        self.params['background']  = 0.0
[829eee9]61       
62
63        ## Parameter details [units, min, max]
64        self.details = {}
[36948c92]65        self.details['scale']       = ['',     None, None]
[1ed3834]66        self.details['radius']      = ['[A]',    None, None]
[36948c92]67        self.details['fractal_dim'] = ['',       0,  None]
[1ed3834]68        self.details['corr_length'] = ['[A]',    None, None]
69        self.details['block_sld']   = ['[1/A²]',  None, None]
70        self.details['solvent_sld'] = ['[1/A²]',  None, None]
[0824909]71        self.details['background']  = ['[1/cm]', None, None]
[829eee9]72       
73               
74    def _Fractal(self, x):
75        """
76            Evaluate 
[3db3895]77            F(x) = p(x) * s(x) + bkd 
[829eee9]78        """
[7292e8a]79        #if x<0 and self.params['fractal_dim']>0:
80         #   raise ValueError, "negative number cannot be raised to a fractional power"
81        #if x==0 and self.params['fractal_dim']==0:
82         #   return 1+self.params['background']
83        #elif x<0 and self.params['fractal_dim']==0:
84        #    return 1e+32
85        #else:
[3db3895]86        return self.params['background']+ self._scatterRanDom(x)* self._Block(x)
[36948c92]87
[829eee9]88   
89    def _Block(self,x):
[7292e8a]90        #if self.params['fractal_dim']<0:
91        #    self.params['fractal_dim']=-self.params['fractal_dim']
92        try:
93            if x<0:
94                x=-x
95            if self.params['radius']<0:
96                self.params['radius']=-self.params['radius']
97               
98            if x==0 or self.params['radius']==0 :
99                 return 1e+32
100            elif self.params['fractal_dim']==0:
101                return 1.0 + (math.sin((self.params['fractal_dim']-1.0) * math.atan(x * self.params['corr_length']))\
102                              * self.params['fractal_dim'] * gamma(self.params['fractal_dim']-1.0))\
103                              *( math.pow( 1.0 + 1.0/((x**2)*(self.params['corr_length']**2)),1/2.0)) 
104            elif self.params['corr_length']==0 or self.params['fractal_dim']==1:
105                return 1.0 + (math.sin((self.params['fractal_dim']-1.0) * math.atan(x * self.params['corr_length']))\
106                              * self.params['fractal_dim'] * gamma(self.params['fractal_dim']-1.0))\
107                              /( math.pow( (x*self.params['radius']), self.params['fractal_dim']))   
108               
109            elif self.params['fractal_dim']<1:
110                return 1.0 + (math.sin((self.params['fractal_dim']-1.0) * math.atan(x * self.params['corr_length']))\
111                              * self.params['fractal_dim'] * gamma(self.params['fractal_dim']-1.0))\
112                              /( math.pow( (x*self.params['radius']), self.params['fractal_dim']))*\
113                                 math.pow( 1.0 + 1.0/((x**2)*(self.params['corr_length']**2)),(1-self.params['fractal_dim'])/2.0)   
114            else:
115                return 1.0 + (math.sin((self.params['fractal_dim']-1.0) * math.atan(x * self.params['corr_length']))\
116                              * self.params['fractal_dim'] * gamma(self.params['fractal_dim']-1.0))\
117                              / math.pow( (x*self.params['radius']), self.params['fractal_dim'])\
118                                 /math.pow( 1.0 + 1.0/((x**2)*(self.params['corr_length']**2)),(self.params['fractal_dim']-1.0)/2.0)   
119        except:
120            return 1 # Need a real fix.
[829eee9]121    def _Spherical(self,x):
122        """
[f629e346]123            F(x) = 3*[sin(x)-xcos(x)]/x**3
[829eee9]124        """
[7292e8a]125        if x==0:
126            return 0
127        else:
128            return 3.0*(math.sin(x)-x*math.cos(x))/(math.pow(x,3.0))
[3db3895]129       
[829eee9]130    def _scatterRanDom(self,x):
131        """
132             calculate p(x)= scale* V^(2)*delta^(2)* F(x*Radius)^(2)
133        """
[3db3895]134        V =(4.0/3.0)*math.pi* math.pow(self.params['radius'],3.0) 
135        delta = self.params['block_sld']-self.params['solvent_sld']
[829eee9]136       
[f629e346]137        return 1.0e8*self.params['scale']* V *(delta**2)*\
[3db3895]138                (self._Spherical(x*self.params['radius'])**2)
[829eee9]139       
140    def run(self, x = 0.0):
141        """ Evaluate the model
[3db3895]142            @param x: input q-value (float or [float, float] as [r, theta])
[829eee9]143            @return: (Fractal value)
144        """
145        if x.__class__.__name__ == 'list':
[36948c92]146            # Take absolute value of Q, since this model is really meant to
147            # be defined in 1D for a given length of Q
[a55fac1]148            #qx = math.fabs(x[0]*math.cos(x[1]))
149            #qy = math.fabs(x[0]*math.sin(x[1]))
[36948c92]150           
[a55fac1]151            return self._Fractal(math.fabs(x[0]))
[829eee9]152        elif x.__class__.__name__ == 'tuple':
153            raise ValueError, "Tuples are not allowed as input to BaseComponent models"
154        else:
155            return self._Fractal(x)
156   
157    def runXY(self, x = 0.0):
158        """ Evaluate the model
[3db3895]159            @param x: input q-value (float or [float, float] as [qx, qy])
[829eee9]160            @return: Fractal value
161        """
162        if x.__class__.__name__ == 'list':
[a55fac1]163            q = math.sqrt(x[0]**2 + x[1]**2)
164            return self._Fractal(q)
[829eee9]165        elif x.__class__.__name__ == 'tuple':
166            raise ValueError, "Tuples are not allowed as input to BaseComponent models"
167        else:
168            return self._Fractal(x)
Note: See TracBrowser for help on using the repository browser.