Source Code for Module reflectometry.model1d.tests.abeles

  1  # This program is public domain. 
  2  """ 
  3  Optical matrix form of the reflectivity calculation. 
  4   
  5  O.S. Heavens, Optical Properties of Thin Solid Films 
  6  """ 
  7  from numpy import * 
  8   
  9 -def refl(Qz,depth,rho,mu=0,wavelength=1,sigma=0): 
 10      """ 
 11      Reflectometry as a function of Qz,wavelength. 
 12   
 13      Qz (inv angstrom) 
 14          Scattering vector 4*pi*sin(theta)/wavelength. This is an array. 
 15      wavelenth (angstrom) 
 16          Incident wavelength. If present this is a 
 17      rho,mu (uNb) 
 18          scattering length density and absorption of each layer 
 19      depth (angstrom) 
 20          thickness of each layer.  The thickness of the incident medium 
 21          and substrate are ignored. 
 22      sigma (angstrom) 
 23          interfacial roughness.  This is the roughness between a layer 
 24          and the subsequent layer.  There is no interface associated 
 25          with the substrate.  The sigma array should have at least n-1 
 26          entries, though it may have n with the last entry ignored. 
 27      method (parratt | abeles | opticalmatrix) 
 28          method used to compute the reflectivity 
 29      """ 
 30      if isscalar(Qz): Qz = array([Qz], 'd') 
 31      n = len(rho) 
 32      nQ = len(Qz) 
 33   
 34      # Make everything into arrays 
 35      kz = array(Qz,'d')/2 
 36      depth = array(depth,'d') 
 37      rho = array(rho,'d') 
 38      mu = mu*ones(n,'d') if isscalar(mu) else array(mu,'d') 
 39      wavelength = wavelength*ones(nQ,'d') \ 
 40          if isscalar(wavelength) else array(wavelength,'d') 
 41      sigma = sigma*ones(n-1,'d') if isscalar(sigma) else array(sigma,'d') 
 42   
 43      # Scale units 
 44      rho *= 1e-6 
 45      mu *= 1e-6 
 46   
 47      ## For kz < 0 we need to reverse the order of the layers 
 48      ## Note that the interface array sigma is conceptually one 
 49      ## shorter than rho,mu so when reversing it, start at n-1. 
 50      ## This allows the caller to provide an array of length n 
 51      ## corresponding to rho,mu or of length n-1. 
 52      idx = (kz>=0) 
 53      r = empty(len(kz),'D') 
 54      r[idx] = calc(kz[idx], wavelength[idx], depth, rho, mu, sigma) 
 55      r[~idx] = calc(abs(kz[~idx]), wavelength[~idx], 
 56                     depth[-1::-1], rho[-1::-1], mu[-1::-1], 
 57                     sigma[n-2::-1]) 
 58      return r 
 59   
 60   
 61 -def calc(kz,wavelength,depth,rho,mu,sigma): 
 62      """Parratt recursion""" 
 63      if len(kz) == 0: return kz 
 64   
 65      ## Complex index of refraction is relative to the incident medium. 
 66      ## We can get the same effect using kz_rel^2 = kz^2 + 4*pi*rho_o 
 67      ## in place of kz^2, and ignoring rho_o 
 68      kz_sq = kz**2 + 4*pi*rho[0] 
 69      k = kz 
 70   
 71      # According to Heavens, the initial matrix should be [ 1 F; F 1], 
 72      # which we do by setting B=I and M0 to [1 F; F 1].  An extra matrix 
 73      # multiply versus some coding convenience. 
 74      B11 = 1 
 75      B22 = 1 
 76      B21 = 0 
 77      B12 = 0 
 78      print "kz",kz[-1] 
 79      for i in xrange(0,len(rho)-1): 
 80          k_next = sqrt(kz_sq - (4*pi*rho[i+1] + 2j*pi*mu[i+1]/wavelength)) 
 81          F = (k - k_next) / (k + k_next) 
 82          F *= exp(-2*k*k_next*sigma[i]**2) 
 83          M11 = exp(1j*k*depth[i]) if i>0 else 1 
 84          print i, kz_sq[-1], rho[i], mu[i], rho[i+1], mu[i+1], wavelength 
 85          print i, "  k",k[-1], "  k_next",k_next[-1] 
 86          if i==0: print i, "  F",F[-1], "   M",1 
 87          if i>0:  
 88              print i, "  F",F[-1], "   M",M11[-1] 
 89          M22 = exp(-1j*k*depth[i]) if i>0 else 1 
 90          M21 = F*M11 
 91          M12 = F*M22 
 92          C1 = B11*M11 + B21*M12 
 93          C2 = B11*M21 + B21*M22 
 94          B11 = C1 
 95          B21 = C2 
 96          C1 = B12*M11 + B22*M12 
 97          C2 = B12*M21 + B22*M22 
 98          B12 = C1 
 99          B22 = C2 
100          k = k_next 
101   
102      r = B12/B11 
103      return r 
104