Source Code for Module snobfit.tests.testSnobfitByRefl

 1  """ 
 2  testrefl.py 
 3  ---------------------------------------------------------------------- 
 4  Python file for running SNOBFIT on a simple reflectometry code 
 5  """ 
 6   
 7  import numpy 
 8  from reflectometry.model1d.model.profile import Profile 
 9  from reflectometry.model1d.tests.parratt import refl 
10  from snobfit.snobfit import snobfit 
11   
12  # Model 
13  D2O = 6.3 # x10^-6 inv A^2 
14  Air = 0.0 
15  Si  = 2.07 # x10-6 inv A^2 
16  Au  = 4.7 
17   
18  # Rhos 
19  rhoDopcBest = [-0.584963,    1.45143,  -0.227783,   -1.32106,  -0.341007, 1.03835,  -0.158153 ] 
20   
21  # sigma, mu, depth  
22  depthDopcBest = [6.97921,    5.75761,    10.1118,    8.50052,    8.02034,    5.15141,  8.38418 ] 
23   
24  # load the dataset 
25  data1 = numpy.loadtxt('qr.si.au.dopc.d2o.42408').T 
26  Q1,R1 = data1 
27   
28  p0 = rhoDopcBest + depthDopcBest 
29   
30   
31  # Profile 
32 -def model_profile(p=None): 
33      if p == None: 
34          common_rhos=[Si, Au] + rhoDopcBest   
35          depths = [17.538,     38.6251] + depthDopcBest + [100] 
36      else: 
37          common_rhos=[Si, Au] 
38          depths = [17.538, 38.6251] 
39          for i in xrange(len(p)/2): 
40              common_rhos.append(p[i]) 
41              depths.append( p[len(p)/2+i] ) 
42          depths.append( 100 ) 
43   
44      rhos1  = common_rhos + [ D2O ] 
45      names1 = ["Si", "Au", "D1", "D2", "D3", "D4", "D5", "D6", "D7", "D20" ] 
46      sigmas = [2.0, 1.521, 2.461, 5.0, 2.05, 2.056, 2.45, 0.789, 0.617] 
47      m1 = Profile( names=names1, 
48                    depth=depths, 
49                    rough=sigmas, 
50                    rho=rhos1, 
51                    mu=[0.0]*10 
52                    ) 
53      z,p1 = m1.calc(n=100)  # ignore mu 
54       
55      return z, p1[0] 
56   
57   
58 -def defaults(): 
59      pbest = rhoDopcBest + depthDopcBest 
60      lo=[]; hi=[] 
61      for i in xrange(len(pbest)): 
62          lo.append( pbest[i]-abs(pbest[i])*0.2 ) 
63          hi.append( pbest[i]+abs(pbest[i])*0.2 ) 
64      return (numpy.array(lo), numpy.array(hi), 0.0) 
65   
66   
67 -def reflectivity( p ): 
68      z,p1= model_profile(p) 
69      thickness = z[1:]-z[:-1] 
70      rho1 = p1[1:] 
71      Rcalc1 = abs(refl(Q1,thickness,rho1))**2 
72      dR1    = numpy.log(Rcalc1) - numpy.log(R1) 
73      chisq  = numpy.dot(dR1,dR1)/len(Rcalc1)  
74       
75      return chisq 
76   
77   
78   
79 -def snobtest( ): 
80      (u,v,fglob) = defaults() 
81      xbest, fbest, ncall = snobfit(reflectivity, p0, (u, v), 
82                                    dn=10, 
83                                    maxiter=1000, 
84                                    retall=True 
85                                    ) 
86   
87      print xbest,fbest, ncall 
88   
89   
90  #---------------------------------------------- 
91  if __name__ == '__main__': 
92      snobtest() 
93