Podklady ke 3. a 4. přednášce - paraxial approximation Zakladní konstanty + načtení balíčku, které se budou používat pri výpočtech In [1]: qe = 1.602e-19; #Elementarni naboj me = 9.109e-31; #Hmotnost elektronu c = 299792458; #Rychlost svetla eta = sqrt(qe/2/me) using Plots #Balicek pro grafy using DifferentialEquations # Balicek pro vypocet diferencialnich rovnic using SpecialFunctions # Specialni funkce #using ForwardDiff # Derivace funkci using PolyChaos #Balicek pro vypocet hermiteovych polynomu using Roots # Balicek pro nalezeni korenu funkce using LaTeXStrings # Latex v grafech Funkce pro derivace error funkce Hermiteovy polynomy - viz. Wikipedia er (x) = (x) exp(− ) , n > 0f (n) 2(−1) (n−1) π −− √ Hn−1 x 2 Hn In [2]: hrca,hrcb = rm_hermite(100) #Hermite recurence coefficients (physical monimial) Hn = (n,x) -> 2^n*PolyChaos.evaluate(n,x,hrca,hrcb) #Hermite polynomial of n-th order function dnerf(n::Int64,x::Float64) #n stupen derivace n>=0, x - bod v nemz funci pocita me if n==0 #Nulta derivace se musi implementovat zvlast return erf(x) else return 2*(-1)^(n-1)/sqrt(pi)*Hn(n-1,x)*exp(-x^2); end end ; In [3]: #Plot error funkce, pripadne jejich derivaci, staci zmenit prvni argument dnerf z1 = (-1:0.01:1)*1e-2 sigma=1e-3 plot(z1,dnerf.(0,z1/(sqrt(2)*sigma))) WARNING: using SpecialFunctions.eta in module Main conflicts with an existing identifie r. Out[3]: Pole v systému Použijeme jednoduchou aproximaci magnetického pole Gausovkou pro výpočet derivaci využijeme derivaci error funkce: tedy B(z) = exp(− )Bm z 2 s 2 B exp(− ) = (z)z 2 π 2 erf ′ (z) = exp(− ) = ( )B (n) Bm d n dz n z 2 s 2 B π −− √ 2s n B erf (n+1) z sB In [4]: # Function for Gaussian field n=stupen derivace, z - bod ve kterem funci pocitam, # Bm - maxumum osoveho pole, sB - parametr urcujici sirku Gausovky dnB = (n,z,Bm,sB) -> Bm/2*sqrt(pi)/sB^n*dnerf(n+1,z/sB) z1 = (-1:0.01:1)*1e-2 p=Any[] Bm=0.2;sB=1e-3 push!(p,plot(z1,dnB.(0,z1,Bm,sB),xlabel="z [m]",ylabel="B [T]")) push!(p,plot(z1,dnB.(1,z1,Bm,sB),xlabel="z [m]",ylabel="B' [T/m]")) plot(p...,layout=(1,2),size=(800,350)) Pro elektrostatickou čočku použijeme jednoduché error funkce, Budeme uvažovat tří elektrodovou elektrostatickou cocku, jejich potencial je ve vektoru V. Uvazujeme hladky prechod pomoci error funkce Φ = (erf((z − z1)/s) + 1)(V [2] − V [1]) + (erf((z − z2)/s) + 1)(V [3] − V [2]) + V [1] 1 2 1 2 Out[4]: In [5]: ze=[-2,2]*1e-3 # z1 a z2 - stredy mezer mezi elektrodami ... s1=1.0e-3; # Rychlost prechodu function dPhiV(n::Int64,z::Float64,V::Array{Float64,1}) if n==0 dP = (dnerf(n,(z-ze[1])/s1)+1)/2/s1^n*(V[2]-V[1])+(dnerf(n,(z-ze[2])/s1)+1)/2/s1 ^n*(V[3]-V[2]) +V[1] else dP = (dnerf(n,(z-ze[1])/s1))/2/s1^n*(V[2]-V[1])+(dnerf(n,(z-ze[2])/s1))/2/s1^n*( V[3]-V[2]) end return dP end #Plot osoveho potencialu, pripadne jeho derivaci dPhi(n,z) = dPhiV(n,z,[1e4; 1.5e4; 1.2e4]) plot(z1,dPhi.(0,z1)) Paraxiální rovnice trajektorie Paraxiální rovnice trajektorie ma tvar: Pro přepsáni do kodu je nutné přejít do soustavy 4 diferencialnich rovnic 1. radu + + x + + y = 0x ′′ γΦ ′ 2Φ ∗ x ′ γΦ ′′ 4Φ ∗ ηB Φ ∗ −− √ y ′ ηB ′ 2 Φ ∗ −− √ + + y − − x = 0y ′′ γΦ ′ 2Φ ∗ y ′ γΦ ′′ 4Φ ∗ ηB Φ ∗ −− √ x ′ ηB ′ 2 Φ ∗ −− √ x ′ y ′ dx ′ dy ′ = dx = dy = − dx − x − dy − y γΦ ′ 2Φ ∗ γΦ ′′ 4Φ ∗ ηB Φ ∗ −− √ ηB ′ 2 Φ ∗ −− √ = − dy − y + dx + x γΦ ′ 2Φ ∗ γΦ ′′ 4Φ ∗ ηB Φ ∗ −− √ ηB ′ 2 Φ ∗ −− √ Out[5]: In [13]: dPhi(n,z) = dPhiV(n,z,[10000.0,10000.0,10000.0]); dBax(n,z) = dnB(n,z,0.3,1e-3); Phir(Phi) = Phi*(1+qe*Phi/(2*me*c^2)) #dBax(n,z) = 0.0; function f(u,p,z) Ph = dPhi(0,z); dPh = dPhi(1,z); d2Ph = dPhi(2,z) Phr = Phir(Ph) B = dBax(0,z); dB = dBax(1,z) du = zeros(Float64,4) du[1] = u[3] du[2] = u[4] du[3] = -dPh/(2*Phr)*u[3]-d2Ph/(4*Phr)*u[1]-eta*B/sqrt(Phr)*u[4]-eta*dB/(2*sqrt(Phr ))*u[2] du[4] = -dPh/(2*Phr)*u[4]-d2Ph/(4*Phr)*u[2]+eta*B/sqrt(Phr)*u[3]+eta*dB/(2*sqrt(Phr ))*u[1] return du end zobrazime nekolik paprsku a) pomoci primeho vypoctu (Vypocet paraxialni rovnice pro kazdou pocatecni podminku zvlast) In [14]: thx=(-1:0.2:1)*5e-3; zspan = (-1.0e-2,1.0e-2) x0=2e-5; y0=0; p=plot() for i=1:length(thx) u0 = [x0,y0,thx[i],0] prob = ODEProblem(f,u0,zspan) sol = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12) plot!(p,sol,vars=(0,1,2)) end plot(p) b) pomoci linearni kombinace reseni. Vsechny řešení lineární diferenciální rovnice tvoří vektorový prostor. Jako bazi zvolíme řešení: : ( ) = [1, 0], ( ) = [0, 0]q ⃗  x q ⃗  x zo g⃗  ′ x zo : ( ) = [0, 1], ( ) = [0, 0]q ⃗  y q ⃗  y zo g⃗  ′ y zo : ( ) = [0, 0], ( ) = [1, 0]h ⃗  x q ⃗  x zo h ⃗  ′ x zo : ( ) = [0, 0], ( ) = [0, 1]h ⃗  y h ⃗  y zo h ⃗  ′ y zo Out[13]: f (generic function with 1 method) Out[14]: In [15]: u0=[1,0,0,0] prob = ODEProblem(f,u0,zspan) solgx = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12); u0=[0,1,0,0] prob = ODEProblem(f,u0,zspan) solgy = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12); u0=[0,0,1,0] prob = ODEProblem(f,u0,zspan) solhx = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12); u0=[0,0,0,1] prob = ODEProblem(f,u0,zspan) solhy = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12); Je vyhodné přejít do komplexních souřadnic ,w = x + iy = x − iyw¯ In [16]: z1 = -0.01:0.0001:0.01 p=plot(); 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(thx) w1 = w.(z1,x0+1im*y0,thx[i]) plot!(p,z1,real(w1),imag(w1)) end plot(p) xlims!((-0.01,0.01)) Rotace svazku v laboratornich souradnicich θ(z) = dz∫ zo z B(z)η 2 2Φ ∗ 1 2 Out[16]: In [19]: # Integral pocitame pomoci reseni diferencialni rovnice (neni to moc efektivni ...) prob3=ODEProblem((u,p,z)->eta*dBax(0,z)/(2*sqrt(Phir(dPhi(0,z)))),0.0,(-0.01,0.01)) thr = solve(prob3,Tsit5(),reltol=1e-8,abstol=1e-10); plot(z1,angle.(w.(z1,0,1e-3)),legend=false) #Vypocet rotace z trajektorii plot!(z1,thr.(z1),legend=false) #Vypoctena rotace (behem integrace par. rovnice) Prechod do rotacnich soradnic In [20]: p=plot() for i=1:length(thx) plot!(p,z1,real(w.(z1,x0,thx[i]).*exp.(-1im*thr.(z1)))) end plot(p) Out[19]: Out[20]: Paraxialni aproximace v rotačních souřadnicích - magnetická čočka V pripadě čistě magnetické čočky dostaneme Pro osové magnetické pole zvolíme jednoduchou aproximaci gausovkou a rovnici trajektorie převedem na dvě ODE prvního řádu a + + u = 0u ′′ γΦ ′ 2Φ ∗ u ′ γ +Φ ′′ η 2 B 2 4Φ ∗ + u = 0u ′′ η 2 B 2 4Φ ∗ B(z) = exp(− / )Bmax z 2 s 2 B u ′ du ′ = du = − u B(z)η 2 4Φ ∗ = 0.3TBmax Φ = 8 keV In [22]: 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]]; Jako bázi vektoru řešení parxiálni rovnice zvolíme řešení , . Pozice predmětu . Libovolný poprsek je pak dán polohou a směrnicí v předmětu: a jeho směrnice g(z) : g( ) = 1, ( ) = 0zo g ′ zo h(z) : h( ) = 0, ( ) = 1zo h ′ zo = −0.01 mzo u(z) = g(z) + h(z)uo u ′ o (z) = (z) + (z)u ′ uo g ′ u ′ o h ′ In [23]: zo = -0.01; zspan = (-1.0e-2,2.0e-2) #Hranice integrace u0=[1,0] #Pocatecni podminka pro trajektorii g #Vypocet trajektorie g prob2 = ODEProblem(f2,u0,zspan) solg = solve(prob2,Tsit5(),reltol=1e-8,abstol=1e-10); u0 = [0,1]#Pocatecni podminky pro trajektorii h #Vypocet trajektorie h prob2 = ODEProblem(f2,u0,zspan) solh = solve(prob2,Tsit5(),reltol=1e-8,abstol=1e-10); 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 Out[23]: #25 (generic function with 1 method) In [24]: #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í In [25]: z1 = (zo:1e-4:0.01) dxo = (-1:0.25:1)*1e-3 xo = [-1, 0, 1]*5e-6 tcol = ["black","blue","red"] 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) Out[24]: Out[25]: Pozice obrazu je daná podmínkou , proč?h( ) = 0zi In [26]: zi = find_zero(h,[0,0.1],Bisection()) (Příčné zvetšení) je pak , a úhlové zvětšeníM = g( )zi = dh( )Ma zi In [27]: M = g(zi); Ma = dh(zi) print("M = ", M, "\nMa = ", Ma, "\n") Použití trajektorií a je analog přechodových matic, jak je znáte za základního kurzu optiky. Přechodová matice udává zobrazení mezi a , pokud je je přechodová matice rovna v případě zobrazení mezi předmětem a obrazem pak je přechodová matice g h z1 z =z1 zo ( ) = ( , z) ( ) x(z) (z)x ′ M ^ zo xo x ′ o ( , z) = ( )M ^ zo g(z) (z)g ′ h(z) (z)h ′ ( , ) = ( )M ^ zo zi M ( )g ′ zi 0 Ma 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: , která jde z mínus nekonečna s nulovou s měrnicí s osou systému a , která jde z nekonečna s nulovou směrnicí s osou , tj. (Jelikož v našem případě je pole v předmětu zanedbatelné, je dánaprůsečíkem s osou.) Obrazové ohnisko je dané průsečíkem s osou a předmětové ohnisko průsečíkem s osou. Ohniskové dálky jsou pak dané směrnicí těchto paprsků: , (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) uπ uπ¯ z (−∞)uπ (∞)uπ¯ = 1, = 1, (−∞)u ′ π (∞)u ′¯ π = 0 = 0 g(z) uπ uπ¯ = − (∞) 1 fi u ′ π = (−∞) 1 fo u ′ π¯ Out[26]: 0.005452623179240787 M = -0.5470826808671513 Ma = -1.8278772665045115 In [28]: zFi = find_zero(g,[0,0.01],Bisection()) print("zFi = ", zFi,"m \n") fi = -1/dg(0.1) print("fi = ", fi,"m \n") u0=[1,0]; zspan2 = (1e-2,-1e-2) prob2 = ODEProblem(f2,u0,zspan2) solbpi = solve(prob2,Tsit5(),reltol=1e-8,abstol=1e-10); ubpi = z -> solbpi(z)[1] dubpi = z -> solbpi(z)[2] zFo = find_zero(ubpi,[-1e-2,0],Bisection()) fo = 1/dubpi(-0.1) print("zFo = ", zFo,"m\n") print("fo = ", fo,"m\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) 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 - nekdy je vhodné definovat svazek pomoci jeho pozice v předmětu a apertuře. Zavedou se trajektorie , a libovolný paprsek je pak dostaneme pomoci jeho pozice v předmětu a apertuře ve tvaru: Trajektorie a můžeme lehce spočítat z trajektorií a . Víme totiž, že i jsou lineární kombinace a : Pro trajektorii pak dostaneme: , a v případě trajektorie : , tedy: za s : s( ) = 1, s( ) = 0zo za t : t( ) = 0, t( ) = 1zo za u(z) = s(z) + t(z)uo ua (z) = (z) + (z)u ′ s ′ uo t ′ ua ( ) = ( ) + ( )u ′ zo s ′ zo uo t ′ zo ua s t g h s t g h s(z) t(z) = ag(z) + bh(z) = cg(z) + dh(z) s 1 = s( ) = ag( ) + bh( ) = azo zo zo 0 = s( ) = g( ) + bh( ) ⇒ b = −za za za g( )za h( )za t 0 = t( ) = cg( ) + dh( ) = czo zo zo 1 = t( ) = dh( ) ⇒ d =za za 1 h( )za s(z) t(z) = g(z) − h(z) g( )za h( )za = h(z) h( )za zFi = 0.0035102421503962556m fi = 0.0035504341423596973m zFo = -0.003510242150391164m fo = 0.003550434142350055m Out[28]: In [29]: za = -1e-3 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 ... In [30]: z1 = (zo:1e-4:0.015) xa = (-1:0.25:1)*1e-5 xo = [-1, 0, 1]*10e-6 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) Out[29]: Out[30]: Pozice v aperture neni zrovna ideální parametr, v elektronové optice je lepší používat úhly. Proto se často místo trajektorie používá přímo trajektorie , tj. Jaký je význam úhlu ? je paprsek, ktery má v předmětu polohu a protíná osu v rovině apertury. je pak jeho směrnice v předmětu. tedy udává úhlovou odchylku v předmětu dané trajektorie od trajektorie parsku . Také se často používají trajktorie spojené s rovinou obrazu: a , pak: jaký je význam a ? Jak získáme a z trajektorií a ? t h u(z) = s(z) + h(z)uo αo αo (z) = (z) + (z)u ′ s ′ uo h ′ αo ( ) = ( ) +u ′ zo s ′ zo uo αo s(z)uo uo ( )s ′ zo uo αo s(z)uo : ( ) = 0, ( ) = 1uα uα zi u ′ α zi : ( ) = 1, ( ) = 0uγ uγ zi uγ za u(z) = γ (z) + α (z)uγ uα γ α (z)uγ (z)uα s(z) h(z) Obecné vlastnosti paraxiálního aproximace Wroskian 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. Tu lze přepsat do tvaru Mějme dvě nezávislá řešení a , která tvoří bazi vektorového prostoru všech řešení. Můžeme pro ně psát Pokud první rovnici vynásobíme , druhou rovnici vynásobíme a sečteme, po krátké úpravě dostaneme rovnici která po integraci vede k zákonu zachování Wronskiánu Lagrange-Helmholtz Relations V případě báze a dostaneme: Což nám určuje vztah mezi úhlovým a příčným zvětšením: Užitím zákonu zachování Wronskiánu na principiální paprsky a roviny , dostaneme: což se využitím definice principiálních trajektoriíredukuje na: Pak dostaneme vztah mezi předmětovou a obrazovou ohniskovou vzdáleností V případě námi vypočtené magnetické čočky: + + u = 0u ′′ γΦ ′ 2Φ ∗ u ′ γ +Φ ′′ η 2 B 2 4Φ ∗ ( ) + ( + ) u = 0Φ ∗ 1 2 d dz Φ ∗ 1 2 u ′ γ0 4 Φ ′′ e 8me B 2 u1 u2 ( ) + ( + ) = 0Φ ∗ 1 2 d dz Φ ∗ 1 2 u ′ 1 γ0 4 Φ ′′ e 8me B 2 u1 ( ) + ( + ) = 0Φ ∗ 1 2 d dz Φ ∗ 1 2 u ′ 2 γ0 4 Φ ′′ e 8me B 2 u2 u2 u1 ( ( − )) = 0Φ ∗ 1 2 d dz Φ ∗ 1 2 u1 u ′ 2 u2 u ′ 1 W = ( − ) = constΦ ∗ 1 2 u1 u ′ 2 u2 u ′ 1 g(z) h(z) W = (g − h ) = ( )(g( ) ( ) − h( ) ( )) = ( )(g( ) ( ) − h( ) ( ))Φ ∗ 1 2 h ′ g ′ Φ ∗ 1 2 zo zo h ′ zo zo g ′ zo Φ ∗ 1 2 zi zi h ′ zi zi g ′ zi ( ) = ( )MΦ ∗ 1 2 zo Φ ∗ 1 2 zi Ma z = −∞ z = ∞ ( (−∞) (−∞) − (−∞) (−∞)) = ( (∞) (∞) − (∞) (∞))Φ ∗ 1 2 −∞ uπ u ′ π¯ uπ¯ u ′ π¯ Φ ∗ 1 2 ∞ uπ u ′ π¯ uπ¯ u ′ π¯ (−∞) = − (∞)Φ ∗ 1 2 −∞ u ′ π¯ Φ ∗ 1 2 ∞ u ′ π = f ¯ f Φ ∗ −∞ Φ ∗ ∞ − −−−− √ In [31]: print("M*Ma = ", M*Ma,"\n") print("fi/f_o = ", fi/fo,"\n") M*Ma = 0.9999999952554085 fi/f_o = 1.0000000000027158 Longitudiální zvětšení Při změně předmětové roviny se také posune rovina obrazu , 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 . 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í víme, že trajektorie v splňuje: Pokud uvažujeme, že je mimo pole, tak druhé derivace jsou nulové a z druhe rovnice dostaneme . Následně pak z první . Trajektorie pak má tvar Pokud vyjádříme tuto trajektorii v novém fokusu po rozvoji do mocnin v a zanedbání kvadratických členů v a dostaneme V případě námi počítané magnetické čočky: → + dzo zo zo → + dzi zi zi h → h ~ (z) = ah(z) + bg(z)h ~ (z)h ~ + dzo zo 0 1 = ( + d ) = ah( + d ) + bg(zo + d ) = ad + bh ~ zo zo zo zo zo zo = ( + d ) = a ( + d ) + b ( + d ) = a( ( ) + ( )d ) + b( ( ) + ( )d ) = a(1 + ( )d )h ~ ′ zo zo h ′ zo zo g ′ zo zo h ′ zo h ′′ zo zo g ′ zo g ′′ zo zo h ′′ zo zo zo a = 1 b = −dzo (z) = h(z) − d g(z)h ~ zo 0 = ( + d ) = h( + d ) + d g( + d )h ~ zi zi zi zi zo zi zi dzi dzi dzo d = dzi zo M 2 Φ ∗ i Φ ∗ o −−− √ In [32]: dzo = 2e-4; zo2 = zo+dzo zspan = (zo2,1e-2); u0 = [0,1]#Pocatecni podminky pro trajektorii h #Vypocet trajektorie h prob2 = ODEProblem(f2,u0,zspan) solh2 = solve(prob2,Tsit5(),reltol=1e-8,abstol=1e-10); h2 = z->solh2(z)[1] #funkce trajektorie h zi2 = find_zero(h2,[0,0.01],Bisection()) print("zo = ", zo," m, zo2 = ",zo2, " m, dzo = ", zo2-zo," m\n") print("zi = ", z1," m, zi2 = ",zi2, " m, dzi = ", zi2-zi," m\n") print("dzi_c = ", M^2*dzo) z1 = zo:1e-4:zi+1e-3; z2 = zo2:1e-4:zi2+1e-3; p=plot() plot!(p,z1,h.(z1)) plot!(p,z2,h2.(z2)) plot(p) zo = -0.01 m, zo2 = -0.0098 m, dzo = 0.00020000000000000052 m zi = -0.01:0.0001:0.015 m, zi2 = 0.0055143864785419236 m, dzi = 6.176329930113623e-5 m dzi_c = 5.985989194095786e-5 Out[32]: In [34]: plot(p) xlims!(zi-2e-3,zi2+1e-3) ylims!(h(zi2+1e-3),h(zi-2e-3)) Aproximace tenkou čočkou V tomto případě zhomogenizovanou rovnici ještě dále upravíme pomoci Pichtovy transformace na tvar kde koeficient je vždy kladný. Předpokládejme, že čočka je tenká, můžeme v ní tedy zanedbat změnu souřadnice , změní se pouze její směrnice pokud tento vztah aplikujeme na principiální paprsky , , použijeme vztah a uvážíme že v nekonečnech je osový potenciál konstantní, tj , můžeme psát Pro ohniskové dálky pak můžeme psát: tedy: Pro naši magnetickou čočku vychází: u(z) = v(z)Φ ∗− 1 4 + G(z)v = 0v ′′ G(z) = (1 + ϵ ) + + 3 16 Φ ′2 Φ ∗2 4 3 Φ ∗ eB 2 8me Φ ∗ Φ1 Φ ¯ 1 8Φ v = − G(z)v(z)dz ≈ − G(z)dzv ′ ∫ −∞ ∞ v0 ∫ −∞ ∞ vπ vπ¯ u = vΦ ∗−1/4 (∞) = (∞)u ′ π Φ ∗−1/4 ∞ v ′ π (−∞) = (−∞)u ′ π¯ Φ ∗−1/4 −∞ v ′ π¯ (∞) = (∞) = − (−∞) G(z)dz = − (−∞) G(z)dzv ′ π u ′ π Φ ∗1/4 ∞ vπ ∫ −∞ ∞ Φ ∗1/4 −∞ uπ ∫ −∞ ∞ (−∞) = (−∞) = (∞) G(z)dz = (∞) G(z)dzv ′ π¯ u ′ π¯ Φ ∗1/4 −∞ vπ¯ ∫ −∞ ∞ Φ ∗1/4 ∞ uπ¯ ∫ −∞ ∞ = G(z)dz, = G(z)dz, 1 f Φ ∗1/4 −∞ Φ ∗1/4 ∞ ∫ −∞ ∞ 1 f ¯ Φ ∗1/4 ∞ Φ ∗1/4 −∞ ∫ −∞ ∞ = f f ¯ Φ ∗1/2 ∞ Φ ∗1/2 −∞ Out[34]: In [35]: G(z)=qe*Bax(z)^2/(8*me*Phr) prob3=ODEProblem((u,p,z)->G(z),0.0,(-0.01,0.01)) iG = solve(prob3,Tsit5(),reltol=1e-8,abstol=1e-10); fo2 = 1/iG(0.01) print("Ohniskova vzdalenost z rovnice trajektorie: ", fo,"\n") print("Ohniskova vzdalenost z aproximace tenkou cockou: ", fo2,"\n") Ohniskova vzdalenost z rovnice trajektorie: 0.003550434142350055 Ohniskova vzdalenost z aproximace tenkou cockou: 0.0032514106373007713