Fitování dielektrické funkce

(příklad ze speciálního praktika - MIR spektroskopie)

In [ ]:
%matplotlib inline
from matplotlib import pyplot as pl
import numpy as np
import os
indir="/home/munz/data/dump/"

naměřená data reflektivity na vzorku dopovaného křemíku

In [2]:
wn,R=np.loadtxt(indir+"fit_ID1.txt").T
pl.plot(wn,R)
pl.xlabel("wavenumber [1/cm]")
pl.ylabel("reflectivity")
pl.grid()
pl.title("doped silicon")
Out[2]:
Text(0.5,1,'doped silicon')

Permeabilita je daná dvěma parametry (+ konstantní složkou $\varepsilon_\infty$) $$\varepsilon(x|\theta_1,\theta_2)=\varepsilon_\infty-\frac{\theta_1}{x\ (x+i \theta_2)}$$

Pro iniciální hodnoty reálná a imaginární hodnota indexu lomu vycházejí takto:

In [7]:
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()
In [4]:
def reflect(n): #kolmy dopad
    return (n-1)/(n+1)
fren=reflect(np.sqrt(eps(om,12,1,0.190)))
drude = lambda nu,eps0,amp,gam: abs(reflect(np.sqrt(eps(nu,eps0,amp,gam))))**2
pl.plot(om,(fren.conjugate()*fren).real)
from scipy import optimize as op
res=op.curve_fit(drude,om,R,[12,1,0.190])
hdru=lambda nu,eps0,amp,gam: np.sum((R-drude(nu,eps0,amp,gam))**2)
eps0,amp,gam=res[0]
pl.plot(om,drude(om,eps0,amp,gam))
pl.plot(om,R)
pl.xlabel("energy [eV]")
pl.grid()
pl.legend(["init. guess","fit","data"])
Out[4]:
<matplotlib.legend.Legend at 0x7fe743822d68>

výsledky fitu

In [6]:
eps0,amp,gam
Out[6]:
(12.10043792034038, 0.21601983705280475, 0.04760652421086287)
In [8]:
xd,yd=0.1,0.1
xp,yp=np.mgrid[1-xd:1+xd:50j,1-yd:1+yd:50j]
vals=[[hdru(om,eps0,xp[i,j]*amp,yp[i,j]*gam) for i in range(50)] for j in range(50)]
hdru(om,eps0,amp,gam)
Out[8]:
0.11480736378368259
In [19]:
extent=[(1-xd)*amp,(1+xd)*amp,(1-yd)*gam,(1+yd)*gam]
pl.imshow(np.array(vals),extent=extent,aspect=3,origin="lower")
pl.colorbar()
pl.contour(np.array(vals),extent=extent)
Out[19]:
<matplotlib.contour.QuadContourSet at 0x7fe7432079e8>

Zvětšíme nyní rozsah testované oblasti asi 5x - optimum je stále ve středu, kontury konstatních hodnot již nevypadají zdaleka symetricky (byť stále připomínají elipsy)

In [11]:
xd,yd=0.5,0.5
xp,yp=np.mgrid[1-xd:1+xd:50j,1-yd:1+yd:50j]
vals2=[[hdru(om,eps0,xp[i,j]*amp,yp[i,j]*gam) for i in range(50)] for j in range(50)]
extent2=[(1-xd)*amp,(1+xd)*amp,(1-yd)*gam,(1+yd)*gam]
pl.imshow(vals2,extent=extent2,aspect=3,origin="lower")
pl.colorbar()
pl.contour(vals2,extent=extent2)
Out[11]:
<matplotlib.contour.QuadContourSet at 0x7fe7434ae240>
In [80]:
[[pl.plot(om,drude(om,eps0,xp[i,j]*amp,yp[i,j]*gam)) for j in [0,49]]  for i in [0,49]];

Rezidua

Rezidua označují rozdíl mezi naměřenou hodnotou a modelem. Pokud je model "přesný", rezidua by měla obsahovat jen náhodnou složku se střední hodnotou 0. Z překryvu fitu a měření v grafu výše se zdá, že náš model kopíruje výsledek měření velmi dobře, po odečtení však vidíme velmi silné systematické efekty na úrovni jednotek procent...

In [18]:
pl.plot(om,R-drude(om,eps0,amp,gam),'.')
pl.xlabel("energy [eV]")
pl.title("residuals")
pl.grid()

Vlnky na vyšších energiích jsou pravděpodobně instrumentální artefakty, ale i po vynechání této části se zdá, že pouze v oblasti pod 70 meV dominuje v reziduích náhodná složka.

In [20]:
bsel=om<0.15
pl.plot(om[bsel],R[bsel]-drude(om[bsel],eps0,amp,gam),'.')
pl.xlabel("energy [eV]")
pl.title("residuals")
pl.grid()

Pokusíme se ještě opravit systematickou odchylku tím, že budeme fitovat jen tuto tuto oblast do 150 meV. Výsledek je ale neúspěšný...

In [23]:
res2=op.curve_fit(drude,om[bsel],R[bsel],res[0])
pl.plot(om[bsel],R[bsel]-drude(om[bsel],*res2[0]),'.')
pl.grid()
res2
Out[23]:
(array([13.04273998,  0.22774103,  0.04784117]),
 array([[7.97882373e-04, 1.05618401e-05, 2.62373841e-07],
        [1.05618401e-05, 1.48861177e-07, 6.42098902e-09],
        [2.62373841e-07, 6.42098902e-09, 2.29683966e-09]]))
In [22]:
res2[0]-res[0]
Out[22]:
array([9.42302058e-01, 1.17211905e-02, 2.34644897e-04])