načtení balíčku, které se budou používat pri výpočtech
using ElectronOptics
using Plots #Balicek pro grafy
pyplot()
using DifferentialEquations # Balicek pro vypocet diferencialnich rovnic
using Roots # Balicek pro nalezeni korenu funkce
using LaTeXStrings # Latex v grafech
using Images
┌ Info: Precompiling ElectronOptics [d9b8bf5d-123d-40cb-a1a8-cd64fe371cc9] └ @ Base loading.jl:1423
Magneticka čočka je zakladní prvek elekronových mikroskopů. Obrazek ukazuje typickou čočku - návrh v programu EOD.
img3d=load("Mlens3D.png");img2d = load("Mlens2D.png")
plot(plot(img3d,axis=[],showaxis=false),plot(img2d,axis=[],showaxis=false),size=(1024,256*1.5))
Pro vyukové účely lze ale použít jednoduchý analytický model pole v magnetické čočce například: $$ B(z) = \frac 12\left[\mathrm{erf}\left(\frac{z+L/2}{\sigma}\right)-\mathrm{erf}\left(\frac{z-L/2}{\sigma}\right)\right]$$
fB = axF(axc=axSCOFF(3.0,8.0),em=0)
plot(plot(fB,xlim=(-20,20)),plot(fB,1,xlim=(-20,20)),plot(fB,2,xlim=(-20,20)),size=(1024,256),layout=(1,3))
pro vyukove učely mužeme take pouzit obdobny analyticky podel jako v pripade magneticke cocky
fE = axF(axc=axSCOFF(20.0,20.0))
plot(fE,xlims=(-40,40))
sol=Any[]
th = (-1.0:0.1:1.0)*8e-3;
for i=1:length(th)
push!(sol,Trace((-40.0,40.0),[0.0,0.0,th[i],0.0,0.0],te_pars(axF(fE,sc=-3342.0),3000.0),rhs=rhs_te,po=ode_pars(dtmax=1.0)))
end
p1=plot(); p2=plot(); z1=(20.3:0.01:20.8)
for i=1:length(sol)
plot!(p1,sol[i],vars=(0,1),label=false,xlabel="z[mm]",ylabel="x[mm]")
plot!(p2,z1,getindex.(sol[i].(z1),1),label=false,xlabel="z[mm]",ylabel="x[mm]")
end
plot(p1,p2,size=(1024,256))
sol=Any[]
th = (-1.0:0.1:1.0)*8e-3;
for i=1:length(th)
push!(sol,Trace((-40.0,40.0),[0.0,0.0,th[i],0.0,0.0],te_pars(axF(fB,sc=0.1185),3000.0),rhs=rhs_te,po=ode_pars(dtmax=1.0)))
end
p1=plot(); p2=plot(); z1=(20.3:0.01:20.8)
for i=1:length(sol)
plot!(p1,sol[i],vars=(0,1,2),label=false,xlabel="z[mm]",ylabel="x[mm]")
plot!(p2,z1,getindex.(sol[i].(z1),1),getindex.(sol[i].(z1),2),label=false,xlabel="z[mm]",ylabel="x[mm]")
end
plot(p1,p2,size=(1024,512))
Paraxiální rovnice trajektorie je lineární aproximace rovnice trajektorie. Lzeji ji odvodit buď z Lagrange-Eulerových rovnic pro index lomu (bereme pouze rozvoj indexu lomu do druhého řádu), nebo přímo z rovnice trajektorie. Pro osově symetrické systémy má tvar: \begin{align} &x^{\prime\prime}+\frac{\gamma\Phi^{\prime}}{2\Phi^{*}}x^{\prime}+\frac{\gamma\Phi^{\prime\prime}}{4\Phi^*}x+\frac{\eta B}{\sqrt{\Phi^*}}y^{\prime}+\frac{\eta B^{\prime}}{2\sqrt{\Phi^*}}y=0\\ &y^{\prime\prime}+\frac{\gamma\Phi^{\prime}}{2\Phi^{*}}y^{\prime}+\frac{\gamma\Phi^{\prime\prime}}{4\Phi^*}y-\frac{\eta B}{\sqrt{\Phi^*}}x^{\prime}-\frac{\eta B^{\prime}}{2\sqrt{\Phi^*}}x=0 \end{align} Pro přepsáni do kodu je nutné přejít do soustavy 4 diferencialnich rovnic 1. radu \begin{align} x^{\prime} &= dx\\ y^{\prime} &= dy\\ dx^{\prime} &= -\frac{\gamma\Phi^{\prime}}{2\Phi^{*}}dx-\frac{\gamma\Phi^{\prime\prime}}{4\Phi^*}x-\frac{\eta B}{\sqrt{\Phi^*}}dy-\frac{\eta B^{\prime}}{2\sqrt{\Phi^*}}y\\ dy^{\prime}&=-\frac{\gamma\Phi^{\prime}}{2\Phi^{*}}dy-\frac{\gamma\Phi^{\prime\prime}}{4\Phi^*}y+\frac{\eta B}{\sqrt{\Phi^*}}dx+\frac{\eta B^{\prime}}{2\sqrt{\Phi^*}}x \end{align}
zobrazime nekolik paprsku a) pomoci primeho vypoctu (Vypocet paraxialni rovnice pro kazdou pocatecni podminku zvlast)
sol=Any[]
th = (-1.0:0.1:1.0)*8e-3;
for i=1:length(th)
push!(sol,PATrace((-40.0,40.0),[0.0,0.0,th[i],0.0,0.0],te_pars(axF(fB,sc=0.1185),3000.0),pte=g_pte_s,po=ode_pars(dtmax=1.0)))
end
p1=plot(); p2=plot(); z1=(20.3:0.01:20.8)
for i=1:length(sol)
plot!(p1,sol[i],vars=(0,1,2),label=false,xlabel="z[mm]",ylabel="x[mm]")
plot!(p2,z1,getindex.(sol[i].(z1),1),getindex.(sol[i].(z1),2),label=false,xlabel="z[mm]",ylabel="x[mm]")
end
plot(p1,p2,size=(1024,512))
b) pomoci linearni kombinace reseni. Vsechny řešení lineární diferenciální rovnice tvoří vektorový prostor. Jako bazi zvolíme řešení: \begin{align} \vec q_x: \vec q_x(z_o) = [1,0],\, \vec g_x^{\prime}(z_o) = [0,0]\\ \vec q_y: \vec q_y(z_o) = [0,1],\, \vec g_y^{\prime}(z_o) = [0,0]\\ \vec h_x: \vec q_x(z_o) = [0,0],\, \vec h_x^{\prime}(z_o) = [1,0]\\ \vec h_y: \vec h_y(z_o) = [0,0],\, \vec h_y^{\prime}(z_o) = [0,1] \end{align}
u0=[1.0,0,0,0,0]
solgx = PATrace((-40.0,40.0),u0,te_pars(axF(fB,sc=0.1185),3000.0),pte=g_pte_s,po=ode_pars(dtmax=1.0))
u0=[0,1.0,0,0,0]
solgy = PATrace((-40.0,40.0),u0,te_pars(axF(fB,sc=0.1185),3000.0),pte=g_pte_s,po=ode_pars(dtmax=1.0))
u0=[0,0,1.0,0,0]
solhx = PATrace((-40.0,40.0),u0,te_pars(axF(fB,sc=0.1185),3000.0),pte=g_pte_s,po=ode_pars(dtmax=1.0))
u0=[0,0,0,1.0,0]
solhy = PATrace((-40.0,40.0),u0,te_pars(axF(fB,sc=0.1185),3000.0),pte=g_pte_s,po=ode_pars(dtmax=1.0));
Je vyhodné přejít do komplexních souřadnic $w = x+\mathrm{i}y$, $\bar w = x-\mathrm{i}y$
z1 = -40.0:0.1:40.0
p=plot();
x0=0.0; y0=0.0;
w = (z,wo,dwo) -> sum((real(wo)*solgx(z) + imag(wo)*solgy(z) + real(dwo)*solhx(z) + imag(dwo)*solhy(z))[1:2].*
[1;1im])
for i=1:length(th)
w1 = w.(z1,x0+1im*y0,th[i])
plot!(p,z1,real(w1),imag(w1),label=false)
end
plot(p)
Rotace svazku v laboratornich souradnicich $$\theta(z) = \int\limits_{z_o}^z\frac{\eta^2B(z)}{2\Phi^{*\frac12}}\mathrm{d} z$$
plot(z1,angle.(w.(z1,0,1e-3)),legend=false) #Vypocet rotace z trajektorii
plot!(z1,getindex.(solgx.(z1),5),legend=false) #Vypoctena rotace (behem integrace par. rovnice)
Prechod do rotacnich soradnic
V pripadě čistě magnetické čočky dostaneme \begin{align} u^{\prime\prime} +\frac{\eta^2B^2}{4\Phi^*}u=0 \end{align}
zo=-40.0
solg = PATrace((zo,40.0),[1.0,0.0],te_pars(axF(fB,sc=0.1185),3000.0),pte=sh_pte,po=ode_pars(dtmax=1.0))
solh = PATrace((zo,40.0),[0.0,1.0],te_pars(axF(fB,sc=0.1185),3000.0),pte=sh_pte,po=ode_pars(dtmax=1.0))
h = z->solh(z)[1] #funkce trajektorie h
dh = z->solh(z)[2] #funkce trajektorie h
g = z->solg(z)[1] #funkce trajektorie g
dg = z->solg(z)[2] #funkce trajektorie g
z1 = (zo:0.1:40.0)
plot(z1,z1*0,label=false,color="black")
plot!(z1,g.(z1),label="g(z)")
plot!(z1,h.(z1)/20,label="h(z)/20")
plot!(z1,dg.(z1)*10,label="g'(z)*10")
plot!(z1,dh.(z1),label="h'(z)")
p=plot()
for i=1:length(th)
plot!(p,z1,real(w.(z1,x0,th[i]).*exp.(-1im*getindex.(solgx.(z1),5))),legend=false)
end
plot(p)
Bm =0.3;sB = 1e-3
Bax(z) = dnB(0,z,Bm,sB);Ph = 8000; Phr = Ph*(1+qe*Ph/(2*me*c^2))
f2(u,p,z) = [u[2],-eta^2*Bax(z)^2/(4*Phr)*u[1]];
#Plot paraxialnich trajektorii a jejich derivaci
#(je treba je naskalovat aby se vlezly do jednoho grafu ...)
z1 = (zo:1e-4:0.01)
plot(z1,z1*0,label=false,color="black")
plot!(z1,g.(z1),label="g(z)")
plot!(z1,h.(z1)*100,label="100 h(z)")
plot!(z1,dg.(z1)/200,label="g'(z)/200")
plot!(z1,dh.(z1),label="h'(z)")
Vykreslení několika trajektorií
dxo = (-1:0.25:1)*1e-3
xo = [-1, 0, 1]*5e-3
tcol = ["black","blue","red"]
z1=-40.0:0.1:40.0
p=plot()
for i=1:length(dxo)
for j=1:length(xo)
plot!(p,z1,xo[j]*g.(z1)+dxo[i]*h.(z1),color=tcol[j],legend=false)
end
end
plot(p)
Pozice obrazu je daná podmínkou $h(z_i) =0$, proč?
zi = find_zero(h,[0,30],Bisection())
20.702591936237162
(Příčné zvetšení) je pak $M = g(z_i)$, a úhlové zvětšení $M_a=dh(z_i)$
M = g(zi); Ma = dh(zi)
print("M = ", M, "\nMa = ", Ma, "\n")
M = -0.5202250416563168 Ma = -1.9222450286138915
Použití trajektorií $g$ a $h$ je analog přechodových matic, jak je znáte za základního kurzu optiky. Přechodová matice udává zobrazení mezi $z_1$ a $z$, pokud je $z_1=z_o$ je přechodová matice rovna $$\begin{pmatrix}x(z)\\x^{\prime}(z)\end{pmatrix}=\hat M(z_o,z)\begin{pmatrix}x_o\\x_o^{\prime}\end{pmatrix}$$ $$ \hat M(z_o,z) = \begin{pmatrix}g(z)&h(z)\\g^{\prime}(z)&h^{\prime}(z)\end{pmatrix}$$ v případě zobrazení mezi předmětem a obrazem pak je přechodová matice $$\hat M(z_o,z_i) = \begin{pmatrix} M&0\\g^{\prime}(z_i)&M_a\end{pmatrix}$$
Můžeme také nalézt základní charakteristiky čočky, jako jou polohy ohnisek a hlavní roviny. Pro tyto účely je vhodné zavést tzv. principiální trajektorie: $u_{\pi}$, která jde z mínus nekonečna s nulovou s měrnicí s osou systému a $u_{\bar\pi}$, která jde z nekonečna s nulovou směrnicí s osou $z$, tj. \begin{align} u_{\pi}(-\infty)&=1,& u_{\pi}^{\prime}(-\infty)&=0\\ u_{\bar\pi}(\infty)&=1,& u_{\pi}^{\bar\prime}(\infty)&=0 \end{align} (Jelikož v našem případě je pole v předmětu zanedbatelné, je dánaprůsečíkem $g(z)$ s osou.) Obrazové ohnisko je dané průsečíkem $u_{\pi}$ s osou a předmětové ohnisko průsečíkem $u_{\bar\pi}$ s osou. Ohniskové dálky jsou pak dané směrnicí těchto paprsků: $\frac 1{f_i} =- u^{\prime}_{\pi}(\infty)$, $\frac 1{f_o} = u^{\prime}_{\bar \pi}(-\infty)$ (vzdálenost od ohniska k rovině, kde se asymptoticky potka prodlouženy paprsek z ohniska s prodlouženym paprskem z nekonečna - pozice hlavni roviny)
zFi = find_zero(g,[0,30.0],Bisection())
print("zFi = ", zFi,"mm \n")
fi = -1/dg(30.0)
print("fi = ", fi,"mm \n")
u0=[1.0,0.0]; zspan2 = (30.0,-30.0)
solbpi = PATrace(zspan2,u0,te_pars(axF(fB,sc=0.1185),3000.0))
ubpi = z -> solbpi(z)[1]
dubpi = z -> solbpi(z)[2]
zFo = find_zero(ubpi,[-40.0,0],Bisection())
fo = 1/dubpi(-30.0)
print("zFo = ", zFo,"mm\n")
print("fo = ", fo,"mm\n")
gr()
plot(z1,g.(z1),label=L"u_{\pi}")
plot!(z1,ubpi.(z1),label=L"u_{\bar\pi}")
plot!(z1,z1*0,color="black",label=false)
zFi = 13.542212621533297mm fi = 13.764003539569321mm zFo = -13.542212621441578mm fo = 13.764003539555443mm
V optice se často používají i jiné báze v prostoru řešení paraxiální rovnice. Význačná je totiž pozice finální apertury $z_a$ - nekdy je vhodné definovat svazek pomoci jeho pozice v předmětu a apertuře. Zavedou se trajektorie $s: s(z_o)=1, s(z_a)=0$, $t: t(z_o)=0, t(z_a)=1$ a libovolný paprsek je pak dostaneme pomoci jeho pozice v předmětu a apertuře ve tvaru: $$ u(z) = s(z) u_o + t(z) u_a$$ $$ u^{\prime}(z) = s^{\prime}(z) u_o + t^{\prime}(z) u_a$$ $$ u^{\prime}(z_o) = s^{\prime}(z_o) u_o + t^{\prime}(z_o) u_a$$ Trajektorie $s$ a $t$ můžeme lehce spočítat z trajektorií $g$ a $h$. Víme totiž, že $s$ i $t$ jsou lineární kombinace $g$ a $h$: \begin{align} s(z) &= ag(z) + bh(z)\\ t(z) &= cg(z) + dh(z) \end{align} Pro trajektorii $s$ pak dostaneme:
$1 = s(z_o) = ag(z_o) + bh(z_o) = a$, $0=s(z_a) = g(z_a)+bh(z_a)\Rightarrow b = -\frac{g(z_a)}{h(z_a)}$
a v případě trajektorie $t$:
$0 = t(z_o) = cg(z_o) + dh(z_o) = c$, $1=t(z_a) = dh(z_a)\Rightarrow d = \frac{1}{h(z_a)}$
tedy: \begin{align} s(z) &= g(z)-\frac{g(z_a)}{h(z_a)}h(z)\\ t(z) &= \frac{h(z)}{h(z_a)} \end{align}
za = -10.0
s(z) = g(z) -g(za)/h(za)*h(z)
t(z) = h(z)/h(za)
plot(z1,t.(z1),label="t(z)")
plot!(z1,s.(z1),label="s(z)")
plot!(z1,z1*0.0,color="black",label=false)
Ted vykreslime nekolik svazku pomovi této parametrizace ...
z1 = (zo:0.1:30)
xa = (-1:0.25:1)*1e-2
xo = [-1, 0, 1]*10e-3
tcol = ["black","blue","red"]
p=plot()
for i=1:length(xa)
for j=1:length(xo)
plot!(p,z1,xo[j]*s.(z1)+xa[i]*t.(z1),color=tcol[j],legend=false)
end
end
xlims!((zo,zi+1e-3))
plot(p)
Pozice v aperture neni zrovna ideální parametr, v elektronové optice je lepší používat úhly. Proto se často místo trajektorie $t$ používá přímo trajektorie $h$, tj. $$ u(z) = s(z) u_o + h(z) \alpha_o$$ Jaký je význam úhlu $\alpha_o$? $$ u^{\prime}(z) = s^{\prime}(z) u_o + h^{\prime}(z) \alpha_o$$ $$ u^{\prime}(z_o) = s^{\prime}(z_o) u_o + \alpha_o$$ $s(z)u_o$ je paprsek, ktery má v předmětu polohu $u_o$ a protíná osu v rovině apertury. $s^{\prime}(z_o)u_o$ je pak jeho směrnice v předmětu. $\alpha_o$ tedy udává úhlovou odchylku v předmětu dané trajektorie od trajektorie parsku $s(z)u_o$.
Také se často používají trajktorie spojené s rovinou obrazu: $u_{\alpha}: u_{\alpha}(z_i)=0,\, u^{\prime}_{\alpha}(z_i) = 1$ a $u_{\gamma}: u_{\gamma}(z_i) = 1, u_{\gamma}(z_a) = 0$, pak: $$u(z) = \gamma u_{\gamma}(z) + \alpha u_{\alpha}(z)$$ jaký je význam $\gamma$ a $\alpha$? Jak získáme $u_{\gamma}(z)$ a $u_{\alpha}(z)$ z trajektorií $s(z)$ a $h(z)$?
I bez numerických výpočtů můžeme říct některé základní vlatnosti paraxiální aproximace. Vyjdem z obecné rovnice pro osově symetrický systém v rotačních souřadnicích. \begin{align} \newcommand{\dd}{\mathrm{d}} u^{\prime\prime} + \frac{\gamma\Phi^{\prime}}{2\Phi^*}u^{\prime}+\frac{\gamma\Phi^{\prime\prime}+\eta^2B^2}{4\Phi^*}u=0 \end{align} Tu lze přepsat do tvaru \begin{align} \Phi^{*\frac 12}\frac{\dd}{\dd z}(\Phi^{*\frac 12} u^{\prime})+\left(\frac{\gamma_0}4\Phi^{\prime\prime}+\frac{e}{8m_e}B^2\right)u = 0 \end{align} Mějme dvě nezávislá řešení $u_1$ a $u_2$, která tvoří bazi vektorového prostoru všech řešení. Můžeme pro ně psát \begin{align} \Phi^{*\frac 12}\frac{\dd}{\dd z}(\Phi^{*\frac 12} u_1^{\prime})+\left(\frac{\gamma_0}4\Phi^{\prime\prime}+\frac{e}{8m_e}B^2\right)u_1 =0\\ \Phi^{*\frac 12}\frac{\dd}{\dd z}(\Phi^{*\frac 12} u_2^{\prime})+\left(\frac{\gamma_0}4\Phi^{\prime\prime}+\frac{e}{8m_e}B^2\right)u_2 =0 \end{align} Pokud první rovnici vynásobíme $u_2$, druhou rovnici vynásobíme $u_1$ a sečteme, po krátké úpravě dostaneme rovnici \begin{align} \Phi^{*\frac 12}\frac{\dd}{\dd z}\left(\Phi^{*\frac 12}(u_1u_2^{\prime}-u_2u_1^{\prime})\right)=0 \end{align} která po integraci vede k zákonu zachování Wronskiánu \begin{align} W = \Phi^{*\frac12}(u_1u_2^{\prime}-u_2u_1^{\prime}) = \rm{const} \end{align}
V případě báze $g(z)$ a $h(z)$ dostaneme: \begin{align} W = \Phi^{*\frac12}(gh^{\prime}-hg^{\prime}) = \Phi^{*\frac12}(z_o)(g(z_o)h^{\prime}(z_o)-h(z_o)g^{\prime}(z_o)) = \Phi^{*\frac12 }(z_i) (g(z_i)h^{\prime}(z_i)-h(z_i)g^{\prime}(z_i)) \end{align} Což nám určuje vztah mezi úhlovým a příčným zvětšením: $$\Phi^{*\frac12}(z_o)=\Phi^{*\frac12 }(z_i) MM_a$$ Užitím zákonu zachování Wronskiánu na principiální paprsky a roviny $z=-\infty$, $z=\infty$ dostaneme: \begin{align} \Phi_{-\infty}^{*\frac12}\left(u_{\pi}(-\infty)u_{\bar{\pi}}^{\prime}(-\infty)-u_{\bar\pi}(-\infty)u^{\prime}_{\bar\pi}(-\infty)\right)= \Phi_{\infty}^{*\frac12}\left(u_{\pi}(\infty)u_{\bar{\pi}}^{\prime}(\infty)-u_{\bar\pi}(\infty)u^{\prime}_{\bar\pi}(\infty)\right) \end{align} což se využitím definice principiálních trajektoriíredukuje na: \begin{align} \Phi_{-\infty}^{*\frac12}u_{\bar{\pi}}^{\prime}(-\infty)=-\Phi_{\infty}^{*\frac12}u_{\pi}^{\prime}(\infty) \end{align} Pak dostaneme vztah mezi předmětovou a obrazovou ohniskovou vzdáleností \begin{align} \frac{\bar f}f=\sqrt{\frac{\Phi^*_{-\infty}}{\Phi^*_{\infty}}} \end{align} V případě námi vypočtené magnetické čočky:
print("M*Ma = ", M*Ma,"\n")
print("fi/f_o = ", fi/fo,"\n")
M*Ma = 1.0000000000843094 fi/f_o = 1.0000000000010083
Při změně předmětové roviny $z_o \rightarrow z_o+\dd z_o$ se také posune rovina obrazu $z_i\rightarrow z_i+\dd z_i$, pokud je tato změna dostatečně malá, je změna obrazové roviny přímo úměrná změně roviny předmětu, kde konstantu úměrnosti nazýváme longitudiálním zvětšením. Pokud tedy posuneme rovinu předmětu změní se i trajektorie $h \rightarrow \tilde h$. Jelikož se také jedná o řešení paraxiální rovnice trajektorie lze ji psát jako lineární kombinaci původních charakteristických trajektorií \begin{align} \tilde h(z) = ah(z)+bg(z) \end{align} víme, že trajektorie $\tilde h(z)$ v $z_o+dz_o$ splňuje: \begin{align} 0&=\tilde h(z_o+dz_o) = ah(z_o+dz_o)+bg(zo+dz_o) =adz_o+b\\ 1&=\tilde h^{\prime}(z_o+dz_o) = ah^{\prime}(z_o+dz_o)+bg^{\prime}(z_o+dz_o)=a(h^{\prime}(z_o)+h^{\prime\prime}(z_o)dz_o)+b(g^{\prime}(z_o) + g^{\prime\prime}(z_o)dz_o) = a(1+h^{\prime\prime}(z_o)dz_o)+b(g^{\prime\prime}(z_o)dz_o) \end{align} Pokud uvažujeme, že $z_o$ je mimo pole, tak druhé derivace jsou nulové a z druhe rovnice dostaneme $a=1$. Následně pak z první $b=-dz_o$. Trajektorie pak má tvar $$\tilde h(z) = h(z)-dz_og(z)$$
Pokud vyjádříme tuto trajektorii v novém fokusu \begin{align} 0=\tilde h(z_i+\dd z_i)=h(z_i+\dd z_i)+\dd z_o g(z_i+\dd z_i) \end{align} po rozvoji do mocnin v $\dd z_i$ a zanedbání kvadratických členů v $\dd z_i$ a $\dd z_o$ dostaneme \begin{align} \dd z_i = \dd z_o M^2\sqrt{\frac{\Phi_i^*}{\Phi_o^*}} \end{align} V případě námi počítané magnetické čočky:
dzo = 0.5;
zo2 = zo+dzo
zspan = (zo2,40.0);
#Vypocet trajektorie h
solh2 = PATrace(zspan,[0.0,1.0],te_pars(axF(fB,sc=0.1185),3000.0));
h2 = z->solh2(z)[1] #funkce trajektorie h
zi2 = find_zero(h2,[0,30.0],Bisection())
print("zo = ", zo," mm, zo2 = ",zo2, " mm, dzo = ", zo2-zo," mm\n")
print("zi = ", z1," mm, zi2 = ",zi2, " mm, dzi = ", zi2-zi," mm\n")
print("dzi_c = ", M^2*dzo)
z1 = zo:1e-1:zi+1; z2 = zo2:1e-1:zi2+1;
p=plot()
plot!(p,z1,h.(z1))
plot!(p,z2,h2.(z2))
plot(p)
zo = -40.0 mm, zo2 = -39.5 mm, dzo = 0.5 mm zi = -40.0:0.1:30.0 mm, zi2 = 20.840515461343404 mm, dzi = 0.13792352510624184 mm dzi_c = 0.13531704698315827
V tomto případě zhomogenizovanou rovnici ještě dále upravíme pomoci Pichtovy transformace \begin{align} u(z) = \Phi^{*-\frac 14} v(z) \end{align} na tvar \begin{align} v^{\prime\prime}+G(z)v=0 \end{align} kde koeficient \begin{align} G(z) = \frac 3{16}\frac{\Phi^{\prime 2}}{\Phi^{*2}}\left(1+\frac 43\epsilon\Phi^*\right)+\frac{eB^2}{8m_e\Phi^*}+\frac{\Phi_1\bar\Phi_1}{8\Phi} \end{align} je vždy kladný. Předpokládejme, že čočka je tenká, můžeme v ní tedy zanedbat změnu souřadnice $v$, změní se pouze její směrnice \begin{align} v^{\prime} = -\int\limits_{-\infty}^{\infty}G(z) v(z) \dd z \approx -v_0\int\limits_{-\infty}^{\infty}G(z) \dd z \end{align} pokud tento vztah aplikujeme na principiální paprsky $v_{\pi}$, $v_{\bar\pi}$, použijeme vztah $u = v\Phi^{\ast -1/4}$ a uvážíme že v nekonečnech je osový potenciál konstantní, tj $u_{\pi}^\prime(\infty) = \Phi^{*-1/4}_{\infty}v^\prime_{\pi}(\infty)$, $u_{\bar\pi}^\prime(-\infty) = \Phi^{*-1/4}_{-\infty}v^\prime_{\bar\pi}(-\infty)$ můžeme psát \begin{align} v_{\pi}^{\prime}(\infty) = u_{\pi}^{\prime}(\infty) \Phi^{*1/4}_{\infty}= -v_{\pi}(-\infty)\int\limits_{-\infty}^{\infty}G(z) \dd z = -\Phi^{*1/4}_{-\infty} u_{\pi}(-\infty)\int\limits_{-\infty}^{\infty}G(z) \dd z\\ v_{\bar\pi}^{\prime}(-\infty) = u_{\bar\pi}^{\prime}(-\infty) \Phi^{*1/4}_{-\infty}= v_{\bar\pi}(\infty)\int\limits_{-\infty}^{\infty}G(z) \dd z = \Phi^{*1/4}_{\infty} u_{\bar\pi}(\infty)\int\limits_{-\infty}^{\infty}G(z) \dd z \end{align} Pro ohniskové dálky pak můžeme psát: \begin{align} \frac 1f=\frac{\Phi^{*1/4}_{-\infty}}{\Phi^{*1/4}_{\infty}}\int\limits_{-\infty}^{\infty}G(z) \dd z, \qquad \frac 1{\bar f}=\frac{\Phi^{*1/4}_{\infty}}{\Phi^{*1/4}_{-\infty}}\int\limits_{-\infty}^{\infty}G(z) \dd z, \end{align} tedy: \begin{align} \frac f{\bar f} = \frac{\Phi^{*1/2}_{\infty}}{\Phi^{*1/2}_{-\infty}} \end{align}
Pro naši magnetickou čočku vychází:
E0=3000.0
fB2 = axF(fB,sc=0.1185)
G(z)=qe*fB2(z)^2/(8*me*E0*scE)
prob3=ODEProblem((u,p,z)->G(z),0.0,(-40.0,40.0))
iG = solve(prob3,Tsit5(),reltol=1e-8,abstol=1e-10);
fo2 = 1/iG(40.0)
print("Ohniskova vzdalenost z rovnice trajektorie: ", fo,"\n")
print("Ohniskova vzdalenost z aproximace tenkou cockou: ", fo2,"\n")
Ohniskova vzdalenost z rovnice trajektorie: 13.764003539555443 Ohniskova vzdalenost z aproximace tenkou cockou: 12.3467577436257