Package reflectometry :: Package model1d :: Package model :: Module model

Source Code for Module reflectometry.model1d.model.model

  1  import numpy 
  2  from copy    import deepcopy 
  3   
  4  from profile import Profile 
  5  from layer   import Layer 
  6  from auxs    import isfactory, combine 
  7  from dataLoader import DataSet 
  8  from calcProfile import build_profile 
  9  from calcRefl    import convolve, reflectivity, magnetic_reflectivity 
 10   
 11  from reflectometry.model1d.profileview.reflutils  import CheckValid 
 12  from staj import Staj, MStaj 
 13   
 14   
15 -class Model:
16 17 wavelength = 4.75 18 background = 0.0 19
20 - def _initModel(self):
21 phi = None 22 theta = None 23 if self.magnetic: 24 phi = [ 0.0, 0.0] 25 theta = [ 0.0, 0.0] 26 27 return Profile( names = ['Air', 'Si'], 28 depth = [100.0, 100.0], 29 rough = [5.0], 30 rho = [ 0.0, 0.0], 31 mu = [ 0.0, 0.0], 32 phi = phi, 33 theta = theta 34 )
35
36 - def __init__(self, 37 datafile = None, 38 magnetic = True 39 ):
40 self.magnetic = magnetic 41 42 # Init the model 43 self.model = self._initModel() 44 45 # data 46 self.data = None 47 48 self._state = 'incident' 49 self.current_interface = 0 50 51 # Datafiles 52 if datafile != None: 53 self.loadData( datafile )
54 55
56 - def addProfile2Layer(self, idx, mL, mP, b ):
57 """ 58 Add a profile in a layer 59 """ 60 if hasattr(b,'build'): 61 62 mL.insert(idx, Layer( b ) ) 63 mP.insert(idx, b ) 64 65 else: 66 val = CheckValid(b) 67 mL.insert(idx, Layer( val) ) 68 mP.insert(idx, val )
69 70
71 - def setProfile2Layer(self, idx, mL, mP, b ):
72 """ 73 set a profile in a layer 74 """ 75 if hasattr(b,'build'): 76 mL[idx] = Layer( b ) 77 mP[idx] = b 78 else: 79 val = CheckValid(b) 80 mL[idx] = Layer(val) 81 mP[idx] = val
82 83
84 - def addLayers(self, 85 name, 86 depth=100.0, 87 rho=0.0, 88 mu=0.0, 89 phi=0.0, 90 theta=0.0 91 ):
92 idx = len(self.model.names)-1 93 self.model.names.insert(idx, name) 94 self.model.depth = numpy.insert( self.model.depth, idx, depth) 95 96 self.addProfile2Layer(idx, self.model.Lrho, self.model.rho, rho ) 97 self.addProfile2Layer(idx, self.model.Lmu, self.model.mu, mu) 98 if self.magnetic: 99 self.addProfile2Layer(idx,self.model.Lphi, self.model.phi, phi) 100 self.addProfile2Layer(idx,self.model.Ltheta,self.model.theta, theta)
101 102 103 # Functions to modify fit parameters
104 - def setRoughness(self, index, rough=5.0):
105 """ 106 set the rough 107 """ 108 self.model.rough = numpy.insert( self.model.rough, index, rough)
109 110
111 - def setDepth(self, index, depth=100 ):
112 """ 113 sets the depth of layer number index 114 """ 115 self.model.depth = numpy.insert( self.model.depth, idx, depth)
116 117
118 - def setValue(self, 119 name, 120 idx = 0, 121 rho = 0.0, 122 mu = 0.0, 123 phi = 0.0, 124 theta = 0.0, 125 depth=100 126 ):
127 self.model.names[idx] = name 128 self.model.depth[idx] = depth 129 self.setProfile2Layer(idx, self.model.Lrho, self.model.rho, rho ) 130 self.setProfile2Layer(idx, self.model.Lmu, self.model.mu, mu) 131 if self.magnetic: 132 self.setProfile2Layer(idx,self.model.Lphi, self.model.phi, phi) 133 self.setProfile2Layer(idx,self.model.Ltheta,self.model.theta, theta)
134 135
136 - def setIncidentValue(self, 137 name, 138 rho = 0.0, 139 mu = 0.0, 140 phi = 0.0, 141 theta = 0.0, 142 depth=100 143 ):
144 """set the profile of the incident layer""" 145 self.setValue(name, 0, rho, mu, phi, theta, depth)
146 147 148
149 - def setSubstrateValue(self, 150 name, 151 rho = 0.0, 152 mu = 0.0, 153 phi = 0.0, 154 theta = 0.0, 155 depth=100 156 ):
157 """set the profile of the substrate layer""" 158 self.setValue(name, -1, rho, mu, phi, theta, depth)
159 160 161
162 - def loadData(self, fname ):
163 """Loads in data with the given filename""" 164 del self.data 165 if self.magnetic: 166 self.data = DataSet(magnetic=True) 167 else: 168 self.data = DataSet(dataset_names = ["Refl"]) 169 170 self.data.load_file(fname)
171 172
173 - def calcProfile(self, maxQ=None, dz=None):
174 """ 175 calculte the profile from model 176 dz (Angstrom) 177 the step size for the profile. 178 maxQ (inv Angstrom) 179 the maximum desired Q. This is equivalent to dz = 4*pi/maxQ. 180 dz default to 1 angstrom profile steps if neither maxQ nor dz are given 181 """ 182 if len(self.model.depth) <= 2: 183 if self.magnetic: return ( None, None, None, None, None ) 184 else: return ( None, None, None ) 185 186 187 # 20 is just a reasonable constant, we can change it to other value 188 assert maxQ is None or dz is None, \ 189 "only one of maxQ and dz allowed" 190 if maxQ != None: 191 dz = 4*numpy.pi/maxQ 192 elif dz == None: 193 dz = 1 194 195 air_depth = self.model.depth[0] 196 total_depth = numpy.sum(self.model.depth) 197 z = numpy.arange(-air_depth,total_depth-air_depth+dz,step=dz) 198 199 if not self.magnetic: 200 ( _rho, _mu ) = [ build_profile( self.model.depth, 201 self.model.rough, 202 c, 203 z, 204 max_rough= 3 ) 205 for c in (self.model.Lrho, self.model.Lmu) ] 206 207 else: 208 (_rho, _mu, _phi, _theta) = [ build_profile(self.model.depth, 209 self.model.rough, 210 c, 211 z, 212 max_rough= 3) 213 for c in (self.model.Lrho, self.model.Lmu, 214 self.model.Lphi, self.model.Ltheta) ] 215 216 z = numpy.ones( len(z) ) * dz 217 z[0] = z[-1] = 0 218 219 if self.magnetic: return (z, _rho, _mu, _phi, _theta) 220 else: return (z, _rho, _mu)
221 222 223
224 - def calcTheory( self, Q, dz=None):
225 226 # get the max Q 227 if dz is None: 228 if not self.magnetic: 229 maxQ = max(Q) 230 else: 231 maxQ=max( [max(Q[0]),max(Q[1]),max(Q[2]),max(Q[3])] ) 232 else: 233 maxQ = None 234 235 if not self.magnetic: 236 (z, _rho, _mu ) = self.calcProfile(maxQ=maxQ,dz=dz) 237 else: 238 (z, _rho, _mu, _phi, _theta) \ 239 = self.calcProfile(maxQ=maxQ,dz=dz) 240 241 #stop if no correct profile 242 if z==None: return 243 244 245 # calc reflectivity 246 if not self.magnetic: 247 R = ( reflectivity(Q, 248 z, 249 _rho, 250 mu=_mu, 251 wavelength=self.wavelength, 252 ) + self.background, 253 ) 254 255 256 else: 257 # Combine Qs and sort it. 258 baseQs = combine( Q ) 259 baseQs.sort() 260 261 Rs = magnetic_reflectivity(baseQs, 262 z, 263 _rho, 264 mu=_mu, 265 wavelength=self.wavelength, 266 rho_m=_phi, 267 theta_m=_theta 268 ) 269 270 # Add background 271 R = numpy.array(Rs) + self.background 272 273 # Do convolution 274 #R = convolve(Q, R, Q, dQ) 275 276 return R
277 278 279
280 - def reflectivity( self, 281 calcQs, 282 background=0.0, 283 wavelength=4.75, 284 fitdQs =None 285 ):
286 287 # get the max Q 288 maxQ = max(calcQs[0]) 289 290 (z,_rhoVal,_muVal ) = self.calcProfie4Model(maxQ) 291 292 if z==None or _rhoVal==None or _muVal==None: 293 return 294 295 296 # Do calc reflectivity 297 R = ( reflectivity(z, 298 _rhoVal, 299 _muVal, 300 wavelength, 301 calcQs[0]) + background 302 ) 303 304 # Do convolution 305 #R = convolve(Q, R, Q, dQ) 306 307 return R
308 309
310 - def incident(self, name, rho=0, mu=0, phi=0, theta=0.0):
311 """ 312 model.incident('name',rho=0,mu=0,phi=0,theta=270) 313 Add the incident medium to the model. 314 """ 315 if self._state != 'incident': 316 raise "Expected "+self._state 317 self._state = 'interface' 318 319 self.setIncidentValue(name, rho, mu, phi, theta, 100)
320 321
322 - def substrate(self, name, rho=2.07, mu=0, phi=0, theta=0.0):
323 """ 324 model.substrate('name',rho=2.07,mu=0,phi=0,theta=270) 325 Add the substrate to the model. 326 """ 327 if self._state != 'layer': 328 raise "Expected "+self._state 329 self._state = 'no more' 330 331 self.setSubstrateValue(name, rho, mu, phi, theta, 100)
332 333 334
335 - def layer(self, name, depth=100, rho=0, mu=0, phi=0, theta=270):
336 """ 337 model.layer('name',depth,rho=0,mu=0,phi=0,theta=270) 338 Add a layer to the model for all profiles. 339 The layer can be flat (give either the value or a range), 340 it can be 'join' indicating that it should extend the preceding 341 layer of the same profile by the depth of the new layer, or it 342 can be layer constructor such as spline(n,[lo,hi]). 343 """ 344 if self._state != 'layer': 345 raise "Expected "+self._state 346 self._state = 'interface' 347 348 self.addLayers( name, depth, rho, mu, phi, theta )
349 350 351
352 - def interface(self, sigma=5.0):
353 """ 354 model.interface(width) 355 set the interface width between layers. The width specifies 356 one standard deviation of a Gaussian smearing function. 357 """ 358 if self._state != 'interface': 359 raise "Expected " + self._state 360 361 self._state = 'layer' 362 363 self.setRoughness(self.current_interface, sigma) 364 365 self.current_interface += 1
366 367 368
369 - def fromPyFile(self, fname):
370 """ 371 Rebuild the object from an python script file 372 """ 373 fin = open( fname, 'r') 374 pyTxts = fin.read() 375 cmds = self._parsePyTexts(pyTxts) 376 for cmd in cmds: 377 exec cmd
378
379 - def _parsePyTexts(self, pyTxts):
380 """ 381 restore information about Refl model. 382 """ 383 s=pyTxts.strip().split('\n') 384 for i in xrange(len(s)): 385 line=s[i] 386 cmd=line.strip().split('(')[0].strip() 387 if cmd in ('incident','layer','interface', 'substrate'): 388 s[i] = 'self.'+ line 389 return s
390
391 - def loadModel4Staj(self, fname, loadData=False):
392 fin = open(fname,'r') 393 _FirstLine = fin.readlines()[0].strip().split() 394 fin.close() 395 396 # FIXME: Need more robust way to test it is magnetic or not? 397 if len(_FirstLine) == 4: #Magnetic Case 398 self.magnetic = True 399 self._LoadMStaj4File( fname ) 400 else: #NonMagnetic Case 401 self.magnetic = False 402 self._LoadNMStaj4File( fname ) 403 404 if loadData: 405 self.loadData(self._stajdatafile)
406 407 408
409 - def _shiftRough(self, rough):
410 _rough = [] 411 for x in xrange( len(rough)-1 ): 412 _rough.append( rough[x+1] ) 413 _rough.append( rough[ len(rough) -1 ] ) 414 415 return _rough
416 417
418 - def _shiftMu(self, mu):
419 _mu = [] 420 for x in xrange(len(mu)-2 ): 421 _mu.append( mu[x+2] ) 422 423 _mu.append( mu[-1] ) 424 _mu.append( mu[-1] ) 425 426 return _mu
427 428
429 - def _LoadNMStaj4File(self, fname):
430 """ Load the model from staj """ 431 s2m = Staj( fname ) 432 (self.wavelength, self.wavelengthDiv, self.angularDiv, 433 self.intensity, self.background) = s2m.getMeta() 434 435 (_depth, _rough, _rho, _mu ) = s2m.getProfiles() 436 437 layerNames = s2m.names #layer names 438 self._stajdatafile = s2m.fin # data file 439 440 441 #Adjust the depth and rough 442 _depth[0] = 100.0 443 _rough_ = self._shiftRough( _rough ) 444 445 self.incident( layerNames[0], rho=_rho[0], mu=_mu[0]) 446 self.interface( _rough_[0] ) 447 for i in xrange( len(_depth)-1 ): 448 self.layer( layerNames[i], 449 depth=_depth[i], 450 rho =_rho[i], 451 mu = _mu[i] 452 ) 453 self.interface( _rough_[i] ) 454 self.substrate(layerNames[-1], rho=_rho[-1], mu=_mu[-1])
455 456 457
458 - def _LoadMStaj4File(self, fname ):
459 """ Load the model from staj ( magnetic case ) """ 460 s2m = MStaj( fname ) 461 (self.wavelength, self.wavelengthDiv, self.angularDiv, 462 self.intensity, self.background, self.Aguide ) = s2m.getMeta() 463 464 (_depth, _rough, _rho, _mu, _phi, _theta) = s2m.getProfiles() 465 466 layerNames = s2m.names #layer names 467 self._stajdatafile = s2m.fin #data file 468 469 470 #Adjust the depth and rough 471 if _depth[0] < 0.00001: 472 _depth[0] = 100.0 473 _rough_ = self._shiftRough( _rough ) 474 475 self.incident( layerNames[0], 476 rho = _rho[0], 477 mu= _mu[0], 478 phi=_phi[0], 479 theta=_theta[0] 480 ) 481 self.interface( _rough_[0] ) 482 483 for i in xrange( 1,len(_depth)-1 ): 484 self.layer( layerNames[i], 485 depth=_depth[i], 486 rho =_rho[i], 487 mu =_mu[i], 488 phi =_phi[i], 489 theta=_theta[i] 490 ) 491 self.interface( _rough_[i] ) 492 493 self.substrate( layerNames[-1], 494 rho=_rho[-1], 495 mu=_mu[-1], 496 phi=_phi[-1], 497 theta=_theta[-1] 498 )
499