Source Code for Module reflectometry.model1d.tests.check2

 1  from numpy import * 
 2  #from reflectometry.model1d.model.calcRefl import reflectivity as refl 
 3  #from reflectometry.model1d.model.parratt import refl 
 4  from reflectometry.model1d.model.abeles import refl 
 5  from reflectometry.model1d.model.calcRefl import magnetic_reflectivity as cr4x 
 6  from reflectometry.model1d.model.polarizedmatrix import magnetic_reflectivity as pyr4x 
 7   
 8  nd=10 
 9   
10  #depth = array([ 100, 65, 35, 17, 100],'d') 
11  depth = array([ 100, 15, 3, 5, 100],'d') 
12  rho   = array([ 0,   1.87, 1.7, 3.8, 2.07],'d') 
13  mu    = array([ 0,   0.01, 0.01, 0.01, 0.01],'d') 
14  phi   = array([ 0.001, 2, 0.001, 0.002, 0.003],'d') 
15  theta = array([200, 210, 220, 230, 240],'d') 
16  #depth,rho,mu,phi,theta = [A[[0,1,-1]] for A in depth,rho,mu,phi,theta] 
17   
18  length = 4 
19  Qc = sqrt(16*pi*rho[-1]*1e-6) 
20  step = 2*pi/min(depth) 
21  start = ceil(3*Qc/step) 
22  lo = start*step 
23  hi = lo+length*step 
24  Qs = linspace(lo,hi,200) 
25  step_n = 2*pi/min(nd*depth) 
26  lo_n = start*step_n 
27  hi_n = lo_n + length*step_n 
28  Qs_n = linspace(lo_n,hi_n,200) 
29  print Qc,lo,hi,lo_n,hi_n 
30  wavelength =  4.75 
31  if False: 
32      R = cr4x(depth,1e-6*rho,1e-6*mu,wavelength,1e-6*phi,theta,-90.,Qs) 
33      R_n = cr4x(nd*depth,1e-6*rho,1e-6*mu,wavelength,1e-6*phi,theta,-90.,Qs_n) 
34  else: 
35      R = [abs(refl(Qs,depth,rho,mu,wavelength))**2] 
36      R_n = [abs(refl(Qs_n,nd*depth,rho,mu,wavelength))**2] 
37   
38 -def fresnel(Q,rho): 
39      """ 
40      Fresnel reflectivity calculator.  This calculator supports very high 
41      Q values with encountering cancellation errors by using a taylor 
42      series expansion for sqrt(1-x).  A similar technique may apply to 
43      the calculation of reflectivity for extremely large Q if such a 
44      thing is ever needed. 
45      """ 
46      kz = abs(Q/2) 
47      # Negate rho if Q<0 
48      rho = choose(Q>=0,[-rho,rho])*1e-6 
49      ksub = sqrt(kz**2 - 4*pi*rho + 0j) 
50      F = (kz-ksub)/(kz+ksub) 
51      # Note: if kz**2>>rho, F ~ pi rho/kz**2 
52      # This is computed as follows: 
53      #    The taylor expansion of sqrt(x) near 1 ~  1 + (x-1)/2 - (x-1)**2/8 + err 
54      #    so sqrt(1-x)  ~ 1 - x/2 - x**2/8 
55      #    F = (1-ksub/kz)/(1+ksub/kz) 
56      #    ksub/kz = (1 - 4*pi*rho/kz**2) ~ 1 - 2 pi rho/kz**2 - 2 (pi rho)**2/kz**4 + err 
57      #    so F ~ (2 pi rho/kz**2)/(2 + 2 pi rho/kz**2) ~ pi rho/kz**2 
58      idx = (kz**2/rho > 1e8) 
59      F[idx] = pi*rho[idx]/kz[idx]**2 
60      return abs(F)**2 
61   
62   
63  import pylab 
64  if True: 
65      # Plot test 
66      pylab.plot(Qs,abs(R[0])/fresnel(Qs,rho[-1]),Qs,-abs(R_n[0])/fresnel(Qs_n,rho[-1])) 
67  else: 
68      # Plot fresnel 
69      Fq = logspace(-5,6,500) 
70      Fr = fresnel(Fq,rho[-1]) 
71      pylab.loglog(Fq,Fr) 
72  pylab.show() 
73