reflectometry.model3d.par

1 # Copyright (C) 2008 University of Maryland 2 # All rights reserved. 3 # See LICENSE.txt for details. 4 5 # Author: Christopher Metting 6 # Starting Date: 4/14/2008 7 8 ''' 9 Form class for a substrate and a feature and functions for operating on it 10 11 Calculates the scattering resulting from the shape profile of individual features. 12 This also calculates the scattering that results from the layers of the substrate. 13 All calculation in this module are in the Born Approximation and for shapes that are 14 parallelapiped. 15 16 ''' 17 18 from numpy import exp,cumsum,concatenate,empty,array,indices,vectorize,pi,ones,product,sin 19 from numpy.fft import fft,fftn,fftshift 20

21 -class Substrate(object):

22 ''' 23 This contains variables to define a substrate 24 25 Dx,Dy (Angstrom) 26 width,length of the full unit cell as defined by the repeat distance 27 thickness (Angstroms) 28 Array containing the thicknesses of all layers in the substrate 29 given in the order: Sample Bottom --> Feature/Substrate Interface 30 SLD_sub (Qc^2) 31 Array containing the scattering length densities which corrispond to the layers in thickness 32 33 34 '''

35 - def __init__(self,Dx = None,Dy = None, thickness = 0, SLD = None):

36 ''' 37 Sets the default to a non-excistant substrate 38 ''' 39 self.Dx = Dx 40 self.Dy = Dy if Dy is not None else Dx 41 self.thickness = thickness 42 self.SLD = SLD if SLD is not None else 0

43 44

45 -class Parallelapiped(object):

46 ''' 47 This contains variables to define a feature 48 49 Rx,Ry (Angstrom) 50 width,length of a feature 51 depth (Angstroms) 52 Array containing the thicknesses of all layers in the feature 53 given in order: Substrate --> Surface 54 SLD (Qc^2) 55 Array containing the scattering length densities which corrispond to the layers in depth 56 given in order: Substrate --> Surface 57 58 59 60 '''

61 - def __init__(self,Rx = None,Ry = None, depth = 0, SLD = None):

62 ''' 63 Sets the default to a non-excistant feature 64 ''' 65 self.Rx = Rx 66 self.Ry = Ry if Ry is not None else Rx 67 self.depth = depth 68 self.SLD = SLD if SLD is not None else 0

69 70

71 -class Lego_Model(object):

72 - def __init__(self, size_xyz = array([0.,0.,0.]), nx = 1, ny = 1, SLD = None):

73 self.size_xyz = size_xyz.astype(float) 74 self.SLD = SLD 75 self.nx = nx 76 self.ny = ny

77 78

79 - def xyz_form_loop(self,qx,qy,qz):

80 #calculations for the feature 81 self.delta_xyz = self.size_xyz / array(self.SLD.shape) 82 self.delta_qxqyqz = 2 * pi / ( self.size_xyz ) 83 #create array that mimics shape of SLD matrix, storing coordinate values for each point\ 84 self.xyz = indices(self.SLD.shape, 'd') 85 self.qxqyqz = indices(self.SLD.shape, 'd') 86 #self.xyz *= self.delta_xyz 87 88 # for all of these arrays, [0], [1], [2] map to x, y, z 89 self.xyz[0] *= self.delta_xyz[0] 90 self.xyz[1] *= self.delta_xyz[1] 91 self.xyz[2] *= self.delta_xyz[2] 92 93 self.qxqyqz[0] *= self.delta_qxqyqz[0] 94 self.qxqyqz[1] *= self.delta_qxqyqz[1] 95 self.qxqyqz[2] *= self.delta_qxqyqz[2] 96 97 #self.qxqyqz = self.xyz * 2. * pi / (self.delta_xyz**2) 98 99 Fxyz = empty([len(qx),len(qy),len(qz)], 'D') 100 101 for i,Qx in enumerate(qx): 102 for j,Qy in enumerate(qy): 103 for k,Qz in enumerate(qz): 104 print(i,Qx,j,Qy,k,Qz) 105 Fxyz[i,j,k] = self.xyz_form(Qx, Qy, Qz) 106 print('Fxyz.shape = ',Fxyz.shape) 107 return Fxyz

108

109 - def xyz_form(self, Qx, Qy, Qz):

110 """ 111 Calculations of \int_zlo^zhi{ rho e^{iQ_z z} dz } 112 """ 113 output_form_matrix = self.SLD * exp(1j * Qx * self.xyz[0]) * exp(1j * Qy * self.xyz[1]) * exp(1j * Qz * self.xyz[2]) 114 #sum the result (which is a matrix of the same shape as SLD) over all axes (x, y and z) 115 output_form_sum = sum(sum(sum(output_form_matrix))) 116 qx_mult = ( 1. - exp( -1j*Qx*self.delta_xyz[0] ) )/(1j * Qx) if Qx!=0 else self.delta_xyz[0] 117 qy_mult = ( 1. - exp( -1j*Qy*self.delta_xyz[1] ) )/(1j * Qy) if Qy!=0 else self.delta_xyz[1] 118 qz_mult = ( 1. - exp( -1j*Qz*self.delta_xyz[2] ) )/(1j * Qz) if Qz!=0 else self.delta_xyz[2] 119 output_form_sum *= qx_mult * qy_mult * qz_mult 120 return output_form_sum

121

122 - def xyz_form_fft(self,q_out,nqx=200,nqy=1,nqz=200):

123 # nqx, etc are the number of cell repeats 124 # q_out = array([[qxmin,qymin,qzmin],[qxmax,qymax,qzmax],[nqx,nqy,nqz]]) 125 #calculations for the feature 126 127 self.delta_xyz = self.size_xyz / array(self.SLD.shape) 128 # for a straight fft, this is the resulting q_stepsize: 129 delta_qxqyqz = 2 * pi / ( self.size_xyz ) 130 # we want a different q spacing: max - min / number 131 delta_qxqyqz_req = ( (q_out[1]).astype(float) - q_out[0] ) / q_out[2] 132 # calculate the oversampled output matrix dimensions: 133 m_new = self.SLD.shape * delta_qxqyqz / delta_qxqyqz_req 134 m_new = m_new.round() 135 m_new[m_new==0] = 1 136 print 'm_new: ', delta_qxqyqz, delta_qxqyqz_req, self.SLD.shape, m_new 137 self.size_qxqyqz = 2 * pi / self.delta_xyz 138 delta_qxqyqz_new = delta_qxqyqz * array(self.SLD.shape) / m_new 139 140 # do the fft! m_new is the size of the resulting matrix 141 # fft is the sum [from n=0 to n=N-1] of x_sub_n * exp(-2i*pi*n*k/N), 142 # where k also goes from 0 to N-1 143 output_form = fftn(self.SLD, m_new.astype(int)) 144 # put the zero-frequency bit in the middle as the creator intended 145 output_form = fftshift(output_form) 146 147 out_shape = array(output_form.shape) 148 149 # indices of max, min elements of output matrix: 150 i_qmax = q_out[1] / delta_qxqyqz_req + out_shape/2 151 i_qmin = q_out[0] / delta_qxqyqz_req + out_shape/2 152 print(i_qmin,i_qmax) 153 154 qxqyqz_out = indices(output_form.shape,'d') - out_shape.reshape((3,1,1,1))/2 155 156 # select out the bits we want from the oversampled fft 157 output_form = output_form[i_qmin[0]:i_qmax[0],i_qmin[1]:i_qmax[1],i_qmin[2]:i_qmax[2]] 158 # also capture the corresponding qx, qy and qz coordinates of the result 159 qxqyqz_out = qxqyqz_out[:,i_qmin[0]:i_qmax[0],i_qmin[1]:i_qmax[1],i_qmin[2]:i_qmax[2]] 160 qxqyqz_out *= delta_qxqyqz_new.reshape((3,1,1,1)) 161 self.qxqyqz_out = qxqyqz_out.copy() 162 163 # the fft is only part of the reflectivity calculation, here's the rest: 164 qxqyqz_mult = empty(qxqyqz_out.shape, 'D') 165 delta_xyz_matrix = ones(qxqyqz_out.shape,'d') * self.delta_xyz.reshape((3,1,1,1)) 166 qxqyqz_mult = ( 1. - exp( -1j*qxqyqz_out*self.delta_xyz.reshape((3,1,1,1)) ) )/(1j * qxqyqz_out) 167 # print 'qxqyqz_mult.shape = ', qxqyqz_mult.shape, 'delta_xyz_matrix.shape = ', delta_xyz_matrix.shape 168 qxqyqz_mult[qxqyqz_out == 0] = delta_xyz_matrix[qxqyqz_out == 0] 169 # multiply the x,y,z multiplier matrices together 170 output_form *= product(qxqyqz_mult, axis=0) 171 172 nxnynz_matrix = ones(qxqyqz_out.shape, 'd') 173 nxnynz_matrix *= array([[[[self.nx]]],[[[self.ny]]],[[[1]]]]) 174 size_xyz_matrix = ones(qxqyqz_out.shape, 'd') 175 size_xyz_matrix *= 1 176 177 #need to calculate structure facture in situ, because of different form of qxqyqz 178 def par_lattice(q,lattice_param,domain): 179 """ 180 Attempts to Calculate the Structure Factor of the sin portion of chucks 181 equations. 182 """ 183 RIF = (sin((q*lattice_param*domain)/2)) / (sin((q*lattice_param)/2)) 184 RIF[q==0] = domain[q==0] 185 print('par_lattice running...') 186 return RIF

187 188 structure_factor = par_lattice(qxqyqz_out,self.size_xyz.reshape((3,1,1,1)),nxnynz_matrix) 189 structure_factor = product(structure_factor, axis=0) 190 191 return output_form,structure_factor

192 193

194 -def form_loop(feature_form,qx,qy,qz):

195 ''' This function creates the form factor for the sample. The substrate form 196 factor and sample form factor are calculated separately. 197 Tz_sub and Tz are vectors returned from the layer_loop function, so 198 total_form_feature[:,i,j] assignment maps onto [qz,qx,qy] (why this order?) 199 ''' 200 Tz = layer_loop(feature_form,qz) 201 total_form_feature = empty([len(qz),len(qx),len(qy)],'d') 202 for i,qxi in enumerate(qx): 203 for j,qyj in enumerate(qy): 204 Tx = inplane_form(qxi,feature_form.Rx) 205 Ty = inplane_form(qyj,feature_form.Ry) 206 total_form_feature[:,i,j] = (Tz*Ty*Tx) 207 return total_form_feature

208

209 -def substrate_loop(substrate,qx,qy,qz):

210 #%:-) calculations for the substrate 211 # z_sub contains the positions of the interfaces for the base 212 sub_total_thickness = sum(substrate.thickness) 213 z_sub = cumsum(concatenate(([0],substrate.thickness))) - sub_total_thickness 214 Fz_sub = 0 # Fz = total integral for rho e^{iQ_z z} over all sections 215 216 for m in range(len(substrate.thickness)): 217 Fz_sub += z_form(substrate.SLD[m], qz, z_sub[m], z_sub[m+1]) 218 219 Tz_sub = Fz_sub/(1j*qz) 220 221 total_form_sub = empty([len(qx),len(qy),len(qz)],'D') 222 for i,qxi in enumerate(qx): 223 for j,qyj in enumerate(qy): 224 Tx_sub = inplane_form(qxi,substrate.Dx) 225 Ty_sub = inplane_form(qyj,substrate.Dy) 226 total_form_sub[i,j,:] = (Tz_sub*Ty_sub*Tx_sub) 227 return total_form_sub

228 229

230 -def layer_loop(feature_form,qz):

231 #calculations for the feature 232 # z contains the positions of the interfaces for the features 233 z = cumsum(concatenate(([0],feature_form.depth))) 234 Fz = 0 235 236 for n in range(len(feature_form.depth)): 237 Fz += z_form(feature_form.SLD[n], qz, z[n], z[n+1]) 238 Tz = Fz/(1j*qz) 239 return Tz

240 241

242 -def z_form(rho, Qz,zlo,zhi):

243 """ 244 Calculations of \int_zlo^zhi{ rho e^{iQ_z z} dz } 245 """ 246 return rho * (exp(1j * Qz * zhi) - exp(1j*Qz *zlo))

247

248 -def inplane_form(q,size):

249 T = (exp(1j*q*size)-1)/(1j*q) if q!=0 else size 250 return T

251

Source Code for Module reflectometry.model3d.par_form