Coverage for larch/xsw/SimpleParratt.py: 0%
176 statements
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
1# 9/20/2010: Y.Choi.
2# Reflectivity calculation using Parratt's recursive formula
3# Calculate reflected intensity as a function of angle or energy
4# Needs accesss to index of refraction data
5# layer0 is vaccum, layer1 is top film layer ...
7import math
8import numpy as np
9from xraydb import xray_delta_beta
11class Layer:
12 def __init__(self, tag, composition='Si', density=2.33, thickness=1000.1234, rms=1e-3):
13 self.tag=tag # layer name, ex), 'top layer','film1'...
14 self.composition=composition # chemical formula ex) SiO2
15 self.density=density # nominal density in g/cm^3
16 self.relden=1.0 # relative density 1.0=100% of nominal density
17 self.thickness=thickness # in Angstroms 1e8 A = 1cm
18 if rms==0: rms=1e-9 # cannot handle 0 roughness.
19 self.rms=rms # rms interface roughness in Agnstroms
20 self.la=1.0 # absoprtion length in cm, for inicident x-ray
21 self.laFY=1.0 # absorption length in Angstrom, for emitted fluorescence
22 self.trans=0.5 # transmission coefficient
23 self.absrp=0.5 # absorption coefficient
24 self.delta=0.1 # index of refraction, real part
25 self.beta=0.1 # index of refraction, imaginary part
26 self.kz=1.0+1.0J # wavenumber, complex [wavenumZ]
27 self.rr=0+0J # Fresnel coefficient, [FresnelR]
28 self.rf=0+0J # Fresnel coefficient with roughness, refl:[FresnelR]+[RoughGauss]
29 self.rrr=0+0J # reflectivity, [Reflect]
30 self.tt=0+0J # Fresnel coeffcient, [FresnelT]
31 self.tf=0+0J # Fresnel coefficient with roughness, for transmission
32 self.ttt=0+0J # transmission, [Transmit]
33 self.E_t=0+0J # electric field, transmission
34 self.E_r=0+0J # electric field, relfection
35 self.interface=0 # interface depth, air/film is 0 for air layer
36 self.index=-1 # air=0, film1=1, ...substrate=n
37 self.Nat=1 # atomic number density for element of interest
38 def set_DenThickRms(self, RelativeDensity, thickness, roughness):
39 self.relden=RelativeDensity
40 self.thickness=thickness
41 self.rms=roughness
43 def get_index(self, energy=8370.0, indexNat=0): # x-ray energy in eV, relative density
44 # MN replace...
45 # temp=readf1f2a.get_delta(self.composition, self.density*self.relden, energy, indexNat) # used to be fluo
46 #self.delta=temp[0]
47 #self.beta=temp[1]
48 #self.la=temp[2] # in cm,
49 #self.Nat=temp[3]
50 delta, beta, la = xray_delta_beta(self.composition,
51 self.density*self.relden, energy)
52 self.delta, self.beta, self.la = delta, beta, la
53 #la attenuation length cm, NumLayer for multiple layers
54 self.trans=math.exp(-self.thickness*1e8/la)
55 self.absrp=1.0-self.trans
57 def cal_kz(self, k0, th_rad): # calculate wavenumber in each layer
58 temp = np.sqrt((math.sin(th_rad))**2 - 2.0*self.delta - 2.0*self.beta*1.0j)
59 self.kz = k0*(temp.real + abs(temp.imag)*1.0J)
60 #kz = k0 * sqrt((sin(th_rad))^2 - 2*delta - 2*beta*i)
61 def cal_rr(self, kz0, kz1):
62 # Fresnel coefficient, [FresnelR] input: kz_(j), kz_(j+1)
63 temp = (kz0-kz1)/(kz0+kz1)
64 self.rr = temp
65 def cal_rf(self, rr, kz0, kz1):
66 # add roughness effect to Fresnel coefficient, [RoughGauss], input:kz_(j), kz_(j+1), rr, rms
67 temp = rr*np.exp(-2.0*kz0*kz1*self.rms)
68 self.rf = temp
69 def cal_rrr(self, rf0, rrr1, kz1, d1): # reflectivity, [Reflect]
70 v = 2.0*kz1*d1
71 u = -v.imag + v.real*1.0J
72 z1 = rrr1*np.exp(u) + rf0
73 z2 = 1.0 + rf0*rrr1*np.exp(u)
74 z=z1/z2
75 self.rrr=z
76 def cal_tt(self, kz0, kz1):
77 #Fresnel coefficient [FresnelT], upto here 0: ThisLayer_(i), 1: BelowLayer_(i+1)
78 temp = 2*kz0/(kz0+kz1)
79 self.tt=temp
80 def cal_ttt(self, tf2, rf2, rrr0, d0, ttt2, kz0):
81 #Fresnel coefficient [FresnelT], 0:ThisLayer_(i), 2: AboveLayer_(i-1)
82 v1 = -2.0*d0*kz0.imag + (2.0*d0*kz0.real)*1.0J
83 v2 = np.exp(v1)
84 v3 = -d0*kz0.imag + (d0*kz0.real)*1.0J
85 v4 = np.exp(v3)
86 v5 = tf2*(ttt2*v4)
87 v6 = rf2*(rrr0*v2)
88 v7 = v6.real+1.0 + (v6.imag)*1.0J
89 z=v5/v7
90 self.ttt=z
93def reflectivity(eV0=14000.0, th_deg=[], mat=[], thick=[], den=[], rough=[], tag=[], depths=[], xsw_mode=0):
94# --------------------set default values------------------------
95 if th_deg==[]: # set default th list
96 th_min=0.0 # minimum th
97 th_max=1.0 # maximum th
98 th_step=0.002 # stepsize th
99 th_deg = np.arange(th_min, th_max+th_step, th_step)
100 th=[] # theta(incident angle) in radian
101 qz=[] # momentum transfer (1/Angstrom)
102 if mat==[] or thick==[] or den==[] or rough==[]:
103 #set defalut layer material, air, film1, film2, substrate
104 tag=['air', 'film', 'underlayer', 'substrate']
105 mat=['N1.56O0.48C0.03Ar0.01Kr0.000001Xe0.0000009', 'Pt', 'Cr', 'Si']
106 thick=[0., 200., 50., 10000] # air, film1, film2, substrate
107 den=[1.e-10, 21.45, 7.19, 2.33]
108 rough=[1.0, 1.0, 1.0, 1.0]
109 if depths==[]: # set default depth list
110 depths=[0.0] # in Angstroms
111 depth_num=len(depths)
112 layers=[]
113 for (ii, d) in enumerate(thick):
114 layer0=Layer(tag[ii], mat[ii], den[ii])
115 layer0.set_DenThickRms(1.0, thick[ii], rough[ii])
116 layers.append(layer0)
117 NumLayers=len(layers) # this value is (# of film layer + 2[air, substrate])
118 FilmThick=0.0 # film thickness
119 for (ix, item) in enumerate(layers):
120 item.get_index(eV0) # get index of refraction, delta, beta for eV0
121 FilmThick+=item.thickness # get total film thickness
122 item.interface=FilmThick # interface location: depth from surface
123 item.index=ix # layer index
124 WaveLength=12.39842/(eV0/1000.0) # in Angstroms, 1e-8cm = 1 Angstrom
125 k0=2.*math.pi/WaveLength
126 for (ix, angle) in enumerate(th_deg): # convert th_deg to th(rad) and qz(1/A)
127 temp=angle/180.0*math.pi
128 th.append(temp) # th in radian
129 temp1=4.*math.pi/WaveLength*math.sin(temp) #qz=4pi/lambda*sin(th)*2k0sin(th)
130 qz.append(temp1)
131 NumTh=len(th)
132 out_parratt='SimpleParratt.txt' # output file name
133 fp=open(out_parratt, 'w')
134 ListOut=[]
135 for angle in th:
136 # define wavenumber within each layer
137 for layer in layers:
138 layer.cal_kz(k0, angle)
139 # calculate reflection and transmission coefficients
140 for ix in range(NumLayers-1):
141 layer=layers[ix]
142 kz0=layers[ix].kz
143 kz1=layers[ix+1].kz
144 layer.cal_rr(kz0, kz1) # Fresnel coefficient, rr
145 layer.cal_rf( layer.rr, kz0, kz1) # Fresnel coefficient with roughness, rf
146 layer.cal_tt(kz0, kz1) # Fresnel coefficient, tt, in e-field intensity
147 layer.tf=layer.rf+1.0 # T=1+R
148 # calculate reflectivity
149 layers[NumLayers-1].rrr=0.0 # nothing reflects back from below the substrate
150 layers[NumLayers-2].rrr=layers[NumLayers-2].rf # interface between the substrate and layer above it
151 ixs = np.arange(NumLayers-3, -1, -1) # start from second interface from the substrate to 0
152 for ix in ixs: # calculate RRR
153 BelowLayer=layers[ix+1]
154 ThisLayer=layers[ix]
155 ThisLayer.cal_rrr(ThisLayer.rf, BelowLayer.rrr, BelowLayer.kz, BelowLayer.thickness)
156 intenR = (layers[0].rrr.real)**2 + (layers[0].rrr.imag)**2
157 # intenR is reflected intensity
158 # calculate transmission
159 layers[0].ttt=1.0 + 0J
160 ixs=np.arange(1, NumLayers-1, 1) # from 1 to <(NumLayers-1)
161 for ix in ixs: # Calculate TTT
162 ThisLayer=layers[ix]
163 AboveLayer=layers[ix-1]
164 ThisLayer.cal_ttt(AboveLayer.tf, AboveLayer.rf, ThisLayer.rrr, \
165 ThisLayer.thickness, AboveLayer.ttt, ThisLayer.kz)
166 ixs=np.arange(0, NumLayers-1, 1) # from 1 to <(NumLayers-1)
167 for ix in ixs:
168 ThisLayer=layers[ix]
169 ThisLayer.E_t=ThisLayer.ttt
170 ThisLayer.E_r=ThisLayer.ttt*ThisLayer.rrr
171 layers[NumLayers-1].E_t = layers[NumLayers-2].tf*layers[NumLayers-2].ttt
172## Next, the current depth within the film structure is searched.
173## depth=0 is air/film, depth>0 film, depth<0 above film
174## xx=0 is film/substrate, xx<0 inside substrate, xx>0 above substrate
175## for reflectivity only, fix depth=0 and skip full xsw calculation
176 layer_index=0 # air/layer1
177 netFY=0.0; EFIat0=0.0; FY=0.0
178 TransFY=1.0; TransFYList=[]
179 EFI_atTH=[] # EFI as a function of depth at current th
180 EFI = 0
181 #print( depth_num)
182 for depth in depths:
183 # calculate e-field intensity at each depth
184 xx=FilmThick-depth
185 ThisLayer = layers[layer_index]
186 if layer_index==(NumLayers-1): #inside substrate
187 cph = -xx*ThisLayer.kz
188 cph_t = np.exp(cph*(0.0+1.0J))
189 Etot = ThisLayer.E_t*cph_t
190 elif layer_index==(NumLayers-2): # layer above the substrate
191 d_n=0
192 cph = ThisLayer.kz*(xx-d_n)
193 cph_t = np.exp(cph*(0.0-1.0J))
194 cph_r = np.exp(cph*(0.0+1.0J))
195 Etot = ThisLayer.E_t*cph_t + ThisLayer.E_r*cph_r
196 else:
197 d_n=FilmThick-ThisLayer.interface
198 cph = ThisLayer.kz*(xx-d_n)
199 cph_t = np.exp(cph*(0.0-1.0J))
200 cph_r = np.exp(cph*(0.0+1.0J))
201 Etot = ThisLayer.E_t*cph_t + ThisLayer.E_r*cph_r
202 EFI = (Etot.real)**2 + (Etot.imag)**2
203 if xsw_mode:
204 pass
205 temp=angle*180.0/math.pi
206 out=str(temp)+' '+str(intenR)+' '+str(EFI)+'\n'
207 ListOut.append([temp, intenR, EFI])
208 fp.write(out)
209 fp.close()
210 #for layer in layers:
211 # print( layer.tag, layer.composition, layer.thickness, layer.rms, layer.density)
212 return ListOut #list of [[th1, refl1], [th2, refl2] ...]
217if __name__=='__main__':
218 #test
219 # reload(readf1f2a)
220 import time
221 a=time.monotonic()
222 reflectivity()
223 print(time.monotonic()-a, 'seconds')