Source Code for Module reflectometry.model1d.model.calcProfile

  1  # This program is public domain 
  2   
  3  """ 
  4  Functions for constructing a profile. 
  5  """ 
  6   
  7  import numpy 
  8  import reflmodule 
  9   
 10  from layer import JoinLayer 
 11   
 12   
 13 -def erf(z): 
 14      """Compute the erf function over the array z""" 
 15      ret = numpy.empty(z.shape,'d') 
 16      reflmodule._erf(z,ret) 
 17      return ret 
 18   
 19   
 20 -def blend(z, rough): 
 21      """ 
 22      blend function 
 23   
 24      Given a Gaussian roughness value, compute the portion of the neighboring 
 25      profile you expect to find in the current profile at depth z. 
 26      """ 
 27      if rough <= 0.0: 
 28          return numpy.where(numpy.greater(z, 0), 0.0, 1.0) 
 29      else: 
 30          return 0.5*( 1.0 - erf( z/( rough*numpy.sqrt(2.0) ) ) ) 
 31   
 32 -def build_profile(depth, 
 33                    rough, 
 34                    layers, 
 35                    z, 
 36                    max_rough=3.0 
 37                    ): 
 38      """ 
 39      Build a profile. 
 40         depth     - depth of the layers (first and last values ignored) 
 41         rough     - roughness of the interfaces (one less than d) 
 42         layers    - objects which calculate profile for a segment 
 43         z         - places to compute the profile (must be sorted) 
 44         max_rough - roughness is limited by layer depth so that max_rough 
 45                     standard deviations fit within the layer. 
 46      """ 
 47      # TODO: test what happens when z array is 'bad': 
 48      # TODO:     z is completely within an interior layer 
 49      # TODO:     z extends beyond the limits of the first/last layer 
 50      # TODO:     z does not contain any points for some layers 
 51   
 52      # Chop out the 'Join' layers.  This could be more easily done using 
 53      # the fact that joins have a span of 0, but this code is meant to 
 54      # work for now, without regard to efficiency. 
 55      R = [ (L, numpy.sum(depth[i:i+L.span]) ) 
 56            for i,L in enumerate(layers) 
 57                if not isinstance(L,JoinLayer) 
 58            ] 
 59   
 60      rough = [rough[i] 
 61               for i,L in enumerate(layers[:-1]) 
 62               if not isinstance(L,JoinLayer)] 
 63   
 64      layers, d = zip(*R) 
 65      d     = numpy.array(depth) 
 66      rough = numpy.array(rough) 
 67   
 68      # Find interface depths 
 69      offset = numpy.concatenate( ((0,),numpy.cumsum(d)) ) - d[0] 
 70   
 71      # Limit roughness to the depths of the surrounding layers.  Roughness 
 72      # of the first and last layers interfaces is limited only by the 
 73      # depth of the first and last layers.  We must check explicitly for 
 74      # a pure substrate system since that has no limits on roughness. 
 75      if max_rough > 0.0 and len(d)>2: 
 76          limit = numpy.min((d[:-1],d[1:]),0)/max_rough 
 77          limit[ 0] = d[ 1]/max_rough 
 78          limit[-1] = d[-2]/max_rough 
 79          rough = numpy.where(rough<limit,rough,limit) 
 80   
 81      # gives the layer boundaries in terms of the index of the z 
 82      idx = numpy.searchsorted(z, offset) 
 83      # TODO: The following hack makes sure the final z value is calculated. 
 84      # TODO: Make sure it works even when z is wider than the range of offsets. 
 85      if idx[-1] < len(z): 
 86          idx[-1] = len(z) 
 87   
 88      # compute the results 
 89      result = numpy.zeros(z.shape,'d') 
 90      for i,L in enumerate(layers): 
 91          zo = z[idx[i]:idx[i+1]] 
 92          if i==0: 
 93              lvalue = 0 
 94              lblend = 0 
 95          else: 
 96              lvalue = layers[i-1].calc(d[i-1],d[i-1]) 
 97              lblend = blend(zo-offset[i],rough[i-1]) 
 98          if i >= len(layers)-1: 
 99              rvalue = 0 
100              rblend = 0 
101          else: 
102              rvalue = layers[i+1].calc(d[i+1],0) 
103              rblend = blend(offset[i+1]-zo,rough[i]) 
104          mvalue = L.calc(d[i],zo-offset[i]) 
105   
106          mblend = 1 - (lblend+rblend) 
107          result[idx[i]:idx[i+1]] = mvalue*mblend + lvalue*lblend + rvalue*rblend 
108          #result[idx[i]:idx[i+1]] = mvalue*mblend 
109   
110      return result 
111