In [59]:
wn,R=np.loadtxt(indir+"fit_E23.txt").T
pl.plot(wn,R)
pl.xlabel("wavenumber [1/cm]")
pl.ylabel("reflectivity")
pl.xlim(0,1300)
pl.grid()
pl.title("layer on doped silicon")
Out[59]:
(příklad ze speciálního praktika - MIR spektroskopie)
Na dopovaný křemík (drudeho formule) je nanesena tlustá vrstva epitaxního (nedopovaného) křemíku.
Pro jednoduchost je dielektrická funkce vrstvy brána (zatím) jako konstanta - přidáváme jen 2 parametry (permeabilitu a tloušťku).
%matplotlib inline
from matplotlib import pyplot as pl
import numpy as np
import os
indir="/home/munz/data/dump/"
wn,R=np.loadtxt(indir+"fit_E23.txt").T
pl.plot(wn,R)
pl.xlabel("wavenumber [1/cm]")
pl.ylabel("reflectivity")
pl.xlim(0,1300)
pl.grid()
pl.title("layer on doped silicon")
c=3e8
#om=c*wn*2*np.pi
om=wn*1240e-7
eps=lambda om,p0,p1,p2:p0-p1/om/(om+1j*p2)
pl.plot(om,np.sqrt(eps(om,12,1,0.190)).real)
pl.plot(om,np.sqrt(eps(om,12,1,0.190)).imag)
pl.xlabel("energy [eV]")
pl.title("Drude dispersion")
pl.legend(list("nk"))
pl.grid()
def reflect(n,amb=1): #kolmy dopad
return (n-amb)/(n+amb)
def reflayer(r01,r12,beta):
return (r01 + r12*np.exp(1j*2*beta))/(1 + r01*r12*np.exp(1j*2*beta))
def indlayer(om,n1,n2,d):
return reflayer(reflect(n1),reflect(n2,n1),2*np.pi*d*om*n1)
eps2=12
#fren=indlayer(om,np.sqrt(eps2),np.sqrt(eps(om,12,0.2,0.040)),.0008/1240e-7)
drlayer = lambda nu,eps0,eps1,amp,gam,thk: abs(indlayer(nu,np.sqrt(eps0),np.sqrt(eps(nu,eps1,amp,gam)),thk/1240e-3))**2
#pl.plot(om,(fren.conjugate()*fren).real)
from scipy import optimize as op
zsel=om<0.1
pl.plot(om[zsel],drlayer(om[zsel],eps2,12,.2,0.050,36))
try:
res=op.curve_fit(drlayer,om[zsel],R[zsel],[eps2,12,.2,0.050,36])
except:
print("something wrong")
else:
eps0,eps1,amp,gam,thk=res[0]
pl.plot(om[zsel],drlayer(om[zsel],eps0,eps1,amp,gam,thk))
pl.plot(om[zsel],R[zsel])
pl.xlabel("energy [eV]")
pl.grid()
pl.legend(["init. guess","fit","data"])
hdru=lambda nu,eps0,eps1,amp,gam,thk: np.sum((R[zsel]-drlayer(nu[zsel],eps0,eps1,amp,gam,thk))**2)
hdru=lambda nu,eps0,eps1,amp,gam,thk: np.sum((R[zsel]-drlayer(nu[zsel],eps0,eps1,amp,gam,thk))**2)
xd,yd=0.6,0.1
xp,yp=np.mgrid[1-xd:1+xd:50j,1-yd:1+yd:50j]
vals=[[hdru(om,eps0,eps1,xp[i,j]*amp,gam,yp[i,j]*thk) for i in range(50)] for j in range(50)]
extent=[(1-xd)*amp,(1+xd)*amp,(1-yd)*thk,(1+yd)*thk]
#extent=None
#pl.imshow(vals,extent=extent,aspect=0.1,origin="lower")
pl.contour(vals,extent=extent,aspect=0.1)
pl.colorbar()
pl.xlabel("Drude strength [1/cm^2]")
pl.ylabel("Thickness [um]")
pl.title("Chi2 map")
xd,yd=0.8,0.5
xp,yp=np.mgrid[1-xd:1+xd:50j,1-yd:1+yd:50j]
vals2=[[hdru(om,eps0,eps1,xp[i,j]*amp,gam,yp[i,j]*thk) for i in range(50)] for j in range(50)]
extent2=[(1-xd)*amp,(1+xd)*amp,(1-yd)*thk,(1+yd)*thk]
#pl.imshow(vals2,extent=extent2,aspect=0.1,origin="lower")
pl.contour(vals2,extent=extent2,aspect=0.1)
pl.colorbar()
pl.xlabel("Drude strength [1/cm^2]")
pl.ylabel("Thickness [um]")
pl.title("Chi2 map")
Porovnání modelových reflektivit v krajních bodech prostoru parametrů..
Rychlejší oscilace odpovídají tlustšímu vzorku, u vrstvy s absorbcí by větší tloušťka znamenala menší amplitudu, tady je to naopak.
eps0,eps1,amp,gam,thk=res[0]
[[pl.plot(om[zsel],drlayer(om[zsel],eps0,eps1,xp[i,j]*amp,gam,yp[i,j]*thk),label="%.2f %.2f"%(xp[i,j],yp[i,j])) for j in [0,49]] for i in [0,49]];
pl.legend()
[[[xp[i,j],yp[i,j]] for j in [0,49]] for i in [0,49]]
pl.plot(om[zsel],R[zsel]-drlayer(om[zsel],eps0,eps1,amp,gam,thk))
pl.grid()
Přidáme nyní k dieletrické funkci základní Lorentzův oscilátor s dalšími 3mi parametry
$$\varepsilon_2(x|\theta_1,\theta_2,\theta_3)=\varepsilon_\infty+\frac{\theta_1}{x^2\ -\theta_2^2 -i \theta_3\ x}.$$Je podstatné, že jsme nejdříve nafitovali jednoduchý konstantní model, až poté přidali další parametry (když není jasné, kde má střed oscilátoru - parametr $\theta_2$ - ležet).
eps_osc=lambda om,p0,p1,p2,p3:p0+p1/(om**2-om*1j*p3-p2**2)
drlayer2 = lambda nu,eps0,eps1,amp,gam,amp2,osc2,gam2,thk: abs(indlayer(nu,np.sqrt(eps_osc(nu,eps0,amp2,osc2,gam2)),np.sqrt(eps(nu,eps1,amp,gam)),thk/1240e-3))**2
#pl.plot(om,(fren.conjugate()*fren).real)
from scipy import optimize as op
zsel=om<0.1
ival=list(res[0][:4])+[0.01,0.2,0.05,res[0][4]]
pl.plot(om[zsel],drlayer2(om[zsel],*ival),label="init.val")
try:
res2=op.curve_fit(drlayer2,om[zsel],R[zsel],ival)
except:
print("something failed")
pl.plot(om[zsel],R[zsel],label="data")
pl.plot(om[zsel],drlayer2(om[zsel],*res2[0]),label="fitted")
pl.legend()
res2[0]
pl.plot(om[zsel],R[zsel]-drlayer2(om[zsel],*res2[0]))
pl.grid()
Výsledný efekt na rezidua je minimální a na uvažovaném intervalu se vlivem oscilátoru změní dielektrická funkce taky velmi málo.
amp2,osc2,gam2=res2[0][-4:-1]
pl.plot(om[zsel],eps_osc(om[zsel],eps2,amp2,osc2,gam2))
pl.grid()
pl.ylabel("eps. real ")
Rozdíl absolutního členu $\varepsilon_\infty$ je oproti variantě bez oscilátoru na úrovni setin (přesněji 0.02):
res[0][0],res2[0][0]
vezmeme místo prvních 100 meV oblast 0-400 meV
zsel=om<0.4
ival=list(res[0][:4])+[0.01,0.2,0.05,res[0][4]]
pl.plot(om[zsel],drlayer2(om[zsel],*ival),label="init.val")
try:
res4=op.curve_fit(drlayer2,om[zsel],R[zsel],ival)
except:
print("something failed")
pl.plot(om[zsel],R[zsel],label="data")
pl.plot(om[zsel],drlayer2(om[zsel],*res2[0]),label="fitted")
pl.legend()
res4[0]
pl.plot(om[zsel],R[zsel]-drlayer2(om[zsel],*res4[0]))
pl.grid()
Periodická komponenta v pravé části spektra není dobře nafitována, zůstává přítomna v reziduích.
Výhodou fitování větší oblasti je přesnější určení parametrů oscilátoru, např. jeho poloha má nyní o řád menší nejistotu.
amp4,osc4,gam4=res4[0][-4:-1]
pl.plot(om[zsel],eps_osc(om[zsel],eps2,amp4,osc4,gam4))
pl.ylabel("eps. real ")
pl.grid()
err2=np.sqrt(res2[1].diagonal())
err4=np.sqrt(res4[1].diagonal())
print("poloha oscilatoru:%.4f +- %.4f eV, vetsi oblast %.4f +- %.4f eV"%(res2[0][-3],err2[-3],res4[0][-3],err4[-3]))