source: sasview/sansmodels/src/sans/models/FractalModel.py @ 7292e8a

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 7292e8a was 7292e8a, checked in by Jae Cho <jhjcho@…>, 15 years ago

half fixed a bug
  • Property mode set to 100644
File size: 6.8 KB
Line 
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:
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
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"
37        self.description="""
38        I(x)= P(x)*S(x) + bkd
39        p(x)= scale* V^(2)*delta^(2)* F(x*Radius)^(2)
40        F(x) = 3*[sin(x)-xcos(x)]/x**3
41        The model has Seven parameters:
42        scale        =  Volume fraction
43        radius       =  Block radius
44        fractal_dim  =  Fractal dimension
45        corr_length  =  correlation Length
46        block_sld    =  SDL block
47        solvent_sld  =  SDL solvent
48        background   =  background
49        """
50        ## Define parameters
51        self.params = {}
52        self.params['scale']       = 0.05
53        self.params['radius']      = 5.0
54        self.params['fractal_dim'] = 2.0
55        self.params['corr_length'] = 100.0
56        self.params['block_sld']   = 2.0e-6
57        self.params['solvent_sld'] = 6.0e-6
58        self.params['background']  = 0.0
59       
60
61        ## Parameter details [units, min, max]
62        self.details = {}
63        self.details['scale']       = ['',     None, None]
64        self.details['radius']      = ['A',    None, None]
65        self.details['fractal_dim'] = ['',       0,  None]
66        self.details['corr_length'] = ['A',    None, None]
67        self.details['block_sld']   = ['A-2',  None, None]
68        self.details['solvent_sld'] = ['A-2',  None, None]
69        self.details['background']  = ['cm-1', None, None]
70       
71               
72    def _Fractal(self, x):
73        """
74            Evaluate 
75            F(x) = p(x) * s(x) + bkd 
76        """
77        #if x<0 and self.params['fractal_dim']>0:
78         #   raise ValueError, "negative number cannot be raised to a fractional power"
79        #if x==0 and self.params['fractal_dim']==0:
80         #   return 1+self.params['background']
81        #elif x<0 and self.params['fractal_dim']==0:
82        #    return 1e+32
83        #else:
84        return self.params['background']+ self._scatterRanDom(x)* self._Block(x)
85
86   
87    def _Block(self,x):
88        #if self.params['fractal_dim']<0:
89        #    self.params['fractal_dim']=-self.params['fractal_dim']
90        try:
91            if x<0:
92                x=-x
93            if self.params['radius']<0:
94                self.params['radius']=-self.params['radius']
95               
96            if x==0 or self.params['radius']==0 :
97                 return 1e+32
98            elif self.params['fractal_dim']==0:
99                return 1.0 + (math.sin((self.params['fractal_dim']-1.0) * math.atan(x * self.params['corr_length']))\
100                              * self.params['fractal_dim'] * gamma(self.params['fractal_dim']-1.0))\
101                              *( math.pow( 1.0 + 1.0/((x**2)*(self.params['corr_length']**2)),1/2.0)) 
102            elif self.params['corr_length']==0 or self.params['fractal_dim']==1:
103                return 1.0 + (math.sin((self.params['fractal_dim']-1.0) * math.atan(x * self.params['corr_length']))\
104                              * self.params['fractal_dim'] * gamma(self.params['fractal_dim']-1.0))\
105                              /( math.pow( (x*self.params['radius']), self.params['fractal_dim']))   
106               
107            elif self.params['fractal_dim']<1:
108                return 1.0 + (math.sin((self.params['fractal_dim']-1.0) * math.atan(x * self.params['corr_length']))\
109                              * self.params['fractal_dim'] * gamma(self.params['fractal_dim']-1.0))\
110                              /( math.pow( (x*self.params['radius']), self.params['fractal_dim']))*\
111                                 math.pow( 1.0 + 1.0/((x**2)*(self.params['corr_length']**2)),(1-self.params['fractal_dim'])/2.0)   
112            else:
113                return 1.0 + (math.sin((self.params['fractal_dim']-1.0) * math.atan(x * self.params['corr_length']))\
114                              * self.params['fractal_dim'] * gamma(self.params['fractal_dim']-1.0))\
115                              / math.pow( (x*self.params['radius']), self.params['fractal_dim'])\
116                                 /math.pow( 1.0 + 1.0/((x**2)*(self.params['corr_length']**2)),(self.params['fractal_dim']-1.0)/2.0)   
117        except:
118            return 1 # Need a real fix.
119    def _Spherical(self,x):
120        """
121            F(x) = 3*[sin(x)-xcos(x)]/x**3
122        """
123        if x==0:
124            return 0
125        else:
126            return 3.0*(math.sin(x)-x*math.cos(x))/(math.pow(x,3.0))
127       
128    def _scatterRanDom(self,x):
129        """
130             calculate p(x)= scale* V^(2)*delta^(2)* F(x*Radius)^(2)
131        """
132        V =(4.0/3.0)*math.pi* math.pow(self.params['radius'],3.0) 
133        delta = self.params['block_sld']-self.params['solvent_sld']
134       
135        return 1.0e8*self.params['scale']* V *(delta**2)*\
136                (self._Spherical(x*self.params['radius'])**2)
137       
138    def run(self, x = 0.0):
139        """ Evaluate the model
140            @param x: input q-value (float or [float, float] as [r, theta])
141            @return: (Fractal value)
142        """
143        if x.__class__.__name__ == 'list':
144            # Take absolute value of Q, since this model is really meant to
145            # be defined in 1D for a given length of Q
146            #qx = math.fabs(x[0]*math.cos(x[1]))
147            #qy = math.fabs(x[0]*math.sin(x[1]))
148           
149            return self._Fractal(math.fabs(x[0]))
150        elif x.__class__.__name__ == 'tuple':
151            raise ValueError, "Tuples are not allowed as input to BaseComponent models"
152        else:
153            return self._Fractal(x)
154   
155    def runXY(self, x = 0.0):
156        """ Evaluate the model
157            @param x: input q-value (float or [float, float] as [qx, qy])
158            @return: Fractal value
159        """
160        if x.__class__.__name__ == 'list':
161            q = math.sqrt(x[0]**2 + x[1]**2)
162            return self._Fractal(q)
163        elif x.__class__.__name__ == 'tuple':
164            raise ValueError, "Tuples are not allowed as input to BaseComponent models"
165        else:
166            return self._Fractal(x)
Note: See TracBrowser for help on using the repository browser.