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
pyplot()
using ElectronOptics
using LaTeXStrings # Latex v grafech
WARNING: using ElectronOptics.eta in module Main conflicts with an existing identifier. WARNING: using ElectronOptics.c in module Main conflicts with an existing identifier. WARNING: using ElectronOptics.me in module Main conflicts with an existing identifier. WARNING: using ElectronOptics.qe in module Main conflicts with an existing identifier.
#=
mutable struct axfld
m::Int64 #Multipole symmetry
axc :: Function #Function for axial coefficients and derivatives
dz :: Float64 #Shift of the field in the z direction
sc :: Number # Scale of the field
em :: Int64 # Electric 1, magnetic 0
zlims :: Array{Float64,1}
function axfld(;axc::T=(n,z)->0.0,m::Int64=0,
dz::Float64=0.0,sc::Float64=1.0,em::Int64=1,zlims::Array{Float64,1}=[-100.0;100.0]) where T<:Function
if zlims[2] == 100
z1 = -1:1e-3:1
id=0
if (m == 0) & (em == 1)
id = 1;
end
fz1 = axc.(id,z1)
mf=maximum(abs.(fz1))
k1 = findall(abs.(fz1).>1e-6*mf)
if length(k1)>2
zlims2 = sort([z1[k1[1]];z1[k1[end]]])
else
zlims2 = [-1.0,1.0]
end
end
new(m,axc,dz,sc,em,zlims2)
end
function axfld(f::axfld;m::Int64=-1,
dz::Float64=NaN,sc::Float64=NaN,em::Int64=-1,zlims::Array{Float64,1}=[NaN])
if m==-1; m=f.m; end
if isnan(dz); dz = f.dz; end
if isnan(sc); sc = f.sc; end
if (em <0) || (em>1); em = f.em; end
if isnan(zlims[1]); zlims = f.zlims; end
new(m,f.axc,dz,sc,em,zlims)
end
end
(f::axfld)(n::Int64,z::Float64)=axceval(n,z,f.axc,f.dz,f.sc,f.zlims,f.m,f.em)
(f::axfld)(z::Float64)=axceval(0,z,f.axc,f.dz,f.sc,f.zlims,f.m,f.em)
function axceval(
n::Int64,
z::Float64,
axc::T,
dz::Float64,
sc::T2,
zlims::Array{Float64,1},
m::Int64,
em::Int64) where {T<: Function, T2<:Number}
zt = z-dz;
if (zt<zlims[1])
if (m==0) & (em==1) & (n==0)
return axc(0,zlims[1])*sc
else
return 0
end
elseif zt>zlims[2]
if (m==0) & (em==1) & (n==0)
return axc(0,zlims[2])*sc
else
return 0
end
else
return axc(n,zt)*sc
end
end
function Base.:*(a::T,fld::axfld) where T<:Number
return axfld(fld,sc=fld.sc*a)
end
function Base.:*(fld::axfld,a::T) where T<:Number
return axfld(fld,sc=fld.sc*a)
end
@recipe function f(f1::axfld,id::Int64=0,sc=1.0)
xguide --> "z [m]"
if f1.em == 1
if id==0
yl = string("\\Phi_{",f1.m,"}\\,")
else
yl = string("\\Phi^{(", id, ")}_{",f1.m,"}\\,")
end
if (id == 0) & (f1.m == 0)
ul = "[\\,V\\,]"
else
ul = string("[\\,V/m^{", id+f1.m,"}\\,]")
end
else
if f1.m == 0
if id==0
yl = "B"
ul = "[\\,T\\,]"
else
yl = string("B^{(", id, ")}")
ul = string("[\\,T/m^{", id, "}\\,]")
end
else
if id==0
yl = string("\\Psi_{",f1.m,"}\\,")
else
yl = string("\\Psi^{(", id, ")}_{",f1.m,"}\\,")
end
ul = string("[\\, V/m^{",f1.m+id,"}\\,]")
end
end
yguide --> latexstring(yl, ul)
xlims --> f1.zlims
(z->f1.axc(id,z)*f1.sc*sc)
end
=#
#=
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 pocitame
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
;
=#
Použijeme jednoduchou aproximaci magnetického pole Gausovkou $$B(z)=B_m\exp\left(-\frac{z^2}{s_B^2}\right)$$ pro výpočet derivaci využijeme derivaci error funkce: $$\exp(-z^2) = \frac{\pi}2\mathrm{erf}^{\prime}(z)$$ tedy $$B^{(n)}(z) = B_m\frac{d^n}{dz^n}\exp\left(-\frac{z^2}{s_B^2}\right)=\frac{\sqrt{\pi}}{2s_B^n}\mathrm{erf}^{(n+1)}(\frac z{s_B})$$
použijeme jednoducou aproximaci pomoci error funkcí $$ \def\erf{\mathrm{erf}} \Phi_m = \Phi_{m,\mathrm{max}}\frac12\left(\erf\left(\frac{z-z_c+\frac L2}{\sigma_d}\right)-\erf\left(\frac{z-z_c-\frac L2}{\sigma_d}\right)+1\right) $$
# 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,zc,sB,Bm) -> Bm/2*sqrt(pi)/sB^n*dnerf(n+1,(z-zc)/sB);
mpaf = (n,z,zc,L,sd,Pm) -> Pm/2*(dnerf(n,(z-zc+L/2)/sd)-dnerf(n,(z-zc-L/2)/sd))
#15 (generic function with 1 method)
Nasleduje funkce prave strany pro obecnou rovnici trajektorie
fB1 = axF(axc=axSCOFF(5.0,4.0),em=0,sc=0.14)
Ph0=10000.0; kappa=0; zspan = (-40.0,20.0)
th=(-1:0.1:1)*1e-4
p=plot()
for i=1:length(th)
u0 = [0.0,0.0,th[i],0,0]
sol = PATrace(zspan,u0,te_pars(fB1,Ph0),pte=g_pte,po=ode_pars(reltol=1e-9,dtmax=1.0))
#prob = ODEProblem(g_pte,u0,zspan,[[fB1],Ph0,kappa])
#sol = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12,dtmax=1e-3);
plot!(p,sol,vars=(0,1),color=:blue,label=false)
end
plot(p)
z1 = zspan[1]:0.1:zspan[2]
xlims!(zspan[1],zspan[2])
ylabel!(L"x_r [ \mathrm{m} ]")
plot!(z1,fB1.(z1)/25,color=:red,label="B/25",ylabel="B [T]")
#savefig("tmp.png")
Do system umistime dve magneticke cocky. První do $z = -800\,\rm mm$, druhou do $z=0$. Budem predpokladat polohu predmetu $z_o=-900\,\rm m$ a prvni cocku nastavime tam aby zobrazovala predmet do $z_{i,1}=-450\rm mm$. Druhou cocku pak nastavime tak, aby zobrazovala $z_{i,1}$ do roviny obrazu $z_i = 10\,\rm mm$.
Wienuv filtr umistime do roviny $z=-700\,\rm mm$ a budeme sledovat, jak ovlivňuje trajektorie častic
zL1 = -800.0; zL2=0.0 #Poloha prvni a druhe magneticke cocky
zWF1 = -700.0; LW=150.0 #Poloha a delka Wiennova filtru
zo = -900.0; zi1=-450.0; zi2=10.0 #Poloha predmetu, prvn9ho a druh0ho obrazu
Ph0 = 1000.0; Phr = Ph0*(1+qe*Ph0/(2*me*c^2)); gamma = 1+qe*Ph0/(me*c^2) #Energie svazku, rel. kor. potencial a rel. faktor gamma
v0 = 2*eta*sqrt(Phr/(1+2*qe*Phr/(me*c^2))); #Rychlost elektronu
#fBa = axfld(axc=(n,z)->dnB(n,z,zL1,1e-3,0.1),em=0) # Pole prvni cocky
#fBb = axfld(axc=(n,z)->dnB(n,z,0.0,1e-3,0.2),em=0) # Pole druhe coky
fBa = axF(axc=axSCOFF(5.0,5.0),em=0, sc=0.15,dz=zL1)
fBb = axF(axc=axSCOFF(5.0,5.0),em=0, sc=0.15,dz=zL2)
Ph1m = 1e4; Ps1m = Ph1m/v0; # Maximalni hodnoty pole ve WF
fE1 = axF(axc=axSCOFF(LW,3.0),m=1,em=1,sc=Ph1m,dz=zWF1)
fB1 = axF(axc=axSCOFF(LW,3.0),m=1,em=0,sc=1im*Ps1m,dz=zWF1)
#fE1 = axfld(axc = (n,z)->mpaf(n,z,zW1,LW,2e-3,Ph1m),m=1) #Elektrostaticke pole WF
#fB1 = axfld(axc = (n,z)->im*mpaf(n,z,zW1,LW,2e-3,Ps1m),m=1,em=0) #Magneticke pole WF
fE2 = axF(axc=axSCOFF(LW,3.0),sc = Ph1m^2/(8*gamma*Phr),m=2,em=1,dz=zWF1)
#fE2 = axfld(axc = (n,z)->mpaf(n,z,zW1,LW,2e-3,Ph1m^2/(8*gamma*Phr)),m=2,em=1) # Elektrostaticke kvadrupolove pole pro splneni stigmaticnosti WF
zp=collect(zo:0.1:zi2)
plot(fBa,0,zp,label="B [T]")
plot!(fBb,0,zp,label="B [T]")
plot!(fE1,0,zp,1e-5,label=L"\Phi_1\cdot 10^{-5}\,[\mathrm{V/m}]")
plot!(zp,200*imag(fB1.(zp)),label=L"\Psi_{1s}\cdot 200\,[\mathrm{T}]")
plot!(fE2,0,zp,1e-5,label=L"\Phi_2\cdot 10^{-5}\,[\mathrm{V/m^2}]")
#xlims!(fBa.zlims[1]-0.02,fBb.zlims[2]+0.02)
#ylabel!(" ")
Řešíme zhomogenizovanou rovnici pro stigmaticky system se silným dipolovým polem \begin{align} u^{\prime\prime}&+\frac{\gamma_0\Phi^{\prime}}{2\Phi^*}u^{\prime}+\left(\frac{\gamma_0\Phi^{\prime\prime}}{4\Phi^*}+\frac{eB^2}{8m_e\Phi^*}+\frac{\Phi_1\bar\Phi_1}{8\Phi^{*2}}\right)u=0 \end{align}
# Nalezeni buzeni prvni cocky, tak aby zobrazovala object do zi1
#sci1 = find_zero(sc->ImgDevTrans(zi1,zo,[sc*fBa,fE1],Ph0,zspan=[zo,zi1]),[0,1],Bisection())
#sol1 = Image(zo,[sci1*fBa,fE1],Ph0,zspan=[zo,zi1+10e-3],SolSave=true,nzi=2)[3];
sc1 = Focus(zo,zi1,sc->te_pars([axF(fBa,sc=sc)],Ph0),(0.01,0.4))
sc2 = Focus(zi1,zi2,sc->te_pars(axF(fBb,sc=sc),Ph0),(0.01,0.4))
zi,ImI = Image(zo,te_pars([axF(fBa,sc=sc1),axF(fBb,sc=sc2)],Ph0),nzi=2)
#Image(zo,te_pars(axF(fBa,sc=0.013),Ph0))
#=
# Nalezeni buzeni druhe cocky, tak aby zobrazovala zi1 do zi1
sci2 = find_zero(sc->ImgDevTrans(zi2,zi1,[sc*fBb],Ph0,zspan=[zi1,zi2]),[0,1],Bisection())
#sol2 = Image(zi1,[sci2*fBb],Ph0,zspan=[zi1,zi2+1e-3],SolSave=true,nzi=2)[3];
prob = ODEProblem(sh_pte,[0.0,1.0],[zi1,zi2+1e-3],[[sci2*fBb],Ph0,0.0])
sol2 = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12,dtmax=1e-3)
h1 = z->sol1(z)[1]
Ma1 = sol1(zi1)[2];
z1 = zo:1e-4:zi1+10e-3
h2 = z->sol2(z)[1]
z2 = zi1:1e-4:zi2+1e-3
=#
z1 = zo:0.1:zi2
th=(-1:0.2:1)*1e-3
plot(z1,ImI["h"].(z1)*th',label=false,color=:blue)
#plot!(sci1*fBa,0,1/250,label=L"B/250 [T]")
#plot!(sci2*fBb,0,1/250,label=L"B/250 [T]")
#plot!(fE1,0,1e-8,label = L"\Phi_1 10^-8 [V/m]")
#xlims!(zo,zi2+1e-3)
zi1
-450.0
Tyto trajeketorie platí pro ideální energii elektronu, v případě že $\kappa\neq 0$ dostaneme
zspan = (zo,zi1+0.01)
th=(-1:0.25:1)*2e-3
p=plot()
for i=1:length(th)
u0 = [0.0,0.0,th[i],0,0]
#prob = ODEProblem(g_pte,u0,zspan,[[sci1*fBa,fE1,fB1,fE2],Ph0,0.0])
#sol = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12,dtmax=1e-3);
#sol=PATrace(zspan,u0,te_pars([axF(fBa,sc=sc1),fE1,fB1,fE2],Ph0,kappa=0.0),pte=g_pte)
sol=PATrace(zspan,u0,te_pars([axF(fBa,sc=sc1),fE1,fB1],Ph0,kappa=0.0),pte=g_pte)
plot!(p,sol,vars=(0,1),label=false,color=:blue)
end
#=
for i=1:length(th)
u0 = [0.0,0.0,th[i],0,0]
prob = ODEProblem(g_pte,u0,zspan,[[sci1*fBa,fBb,fE1,fB1,fE2],Ph0,-1.0/Ph0])
sol = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12,dtmax=1e-3);
plot!(p,sol,vars=(0,1),label=false,color=:red)
end
for i=1:length(th)
u0 = [0.0,0.0,th[i],0,0]
prob = ODEProblem(g_pte,u0,zspan,[[sci1*fBa,fBb,fE1,fB1,fE2],Ph0,1.0/Ph0])
sol = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12,dtmax=1e-3);
plot!(p,sol,vars=(0,1),label=false,color=:purple)
end
=#
plot(p)
fE2.sc
12463.39325219401
Pokud bychom do obrazove roviny prvni cocky umistili aperturu, nebo sterbinu , mohli bychom odfiltrovat elektrony s vetsí energiovou deviaci $\kappa$ a do systemu pustit jen svazek s uzsim energiovym spektrem. Predpokladejme vyslednou energiovou sirku 40 meV.
zspan = (zo,zi2+0.001); dE=[0, -0.02, 0.02]; tcol=[:blue,:red,:purple]
px=plot()
for j=1:length(dE)
for i=1:length(th)
u0 = [0.0,0.0,th[i],0,0]
prob = ODEProblem(g_pte,u0,zspan,[[sci1*fBa,sci2*fBb,fE1,fB1,fE2],Ph0,dE[j]/Ph0])
sol = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12,dtmax=1e-3);
plot!(px,sol,vars=(0,1),label=false,color=tcol[j])
end
end
plot(px)
plot(px)
xlims!(9.99e-3,10.01e-3); ylims!(-0.2e-6,0.2e-6)
Zbytkovou energiovou disperzy je nutne vyrusit druhym WF
#2 Wienuv filtr
zW2 = -0.2;
fE1b = axfld(fE1,dz=zW2-zW1); fB1b = axfld(fB1,dz=zW2-zW1);
fE2b = axfld(fE2,dz=zW2-zW1)
sci2b = find_zero(sc->ImgDevTrans(zi2,zi1,[sc*fBb,fE1b],Ph0,zspan=[zi1,zi2]),[0,1],Bisection())
0.30798817659008554
zspan = (zo,zi2+0.001); dE=0.04
px=plot()
kappa = [0.0,-dE/2,dE/2]/Ph0;
trcol=[:blue,:red,:green]
fld = [sci1*fBa,sci2b*fBb,fE1,fB1,fE2,fE1b,fB1b,fE2b]
for ik = 1:length(kappa)
for i=1:length(th)
u0 = [0.0,0.0,th[i],0,0]
prob = ODEProblem(g_pte,u0,zspan,[fld,Ph0,kappa[ik]])
sol = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12,dtmax=1e-3);
plot!(px,sol,vars=(0,1),label=false,color=trcol[ik])
end
end
plot(px,xlabel="z [m]",ylabel="x [m]")
plot(px)
xlims!(9.995e-3,10.008e-3); ylims!(-0.1e-6,0.1e-6)
Výsledné řešení pak tedy můžeme psát ve tvaru \begin{align} \def\dd{\mathrm{d}} u(z) = u_o g + u_o^{\prime}h +u_d\kappa \end{align} kde $u_d$ je disperzní trajektorie \begin{align} u_d=\Phi^{*-\frac12}_og\int\limits_{z_o}^z\Phi^{*\frac 12}h\frac{\Phi_0\Phi_1}{4\Phi^{*2}}e^{-\im\chi(z)}\dd z-\Phi^{*-\frac12}_oh\int\limits_{z_o}^z\Phi^{*\frac 12}g\frac{\Phi_0\Phi_1}{4\Phi^{*2}}e^{-\im\chi(z)}\dd z \label{DispTraj} \end{align} $g$ a $h$ pak jsou reseni zhomogenizovane paraxialni rovice trajketorie \begin{align} u^{\prime\prime}&+\frac{\gamma_0\Phi^{\prime}}{2\Phi^*}u^{\prime}+\left(\frac{\gamma_0\Phi^{\prime\prime}}{4\Phi^*}+\frac{eB^2}{8m_e\Phi^*}+\frac{\Phi_1\bar\Phi_1}{8\Phi^{*2}}\right)u=-\frac{\Phi_0\Phi_1}{4\Phi^{*2}}e^{-\im\chi(z)}\kappa \label{PTrajSARStig} \end{align}
fld = [sci1*fBa,sci2b*fBb, #Pole magnetickych cocek
fE1,fB1,fE2, #Pole prvniho WF
fE1b,fB1b,fE2b #Pole druheho WF
]
prob = ODEProblem(sh_pte,[0.0,1.0],(zo,zi2),[fld,Ph0])
solh = solve(prob,Tsit5(),reltol=1e-1,abstol=1e-12,dtmax=1e-3)
h=z->solh(z)[1]
dh = z->sol(z)[2]
prob = ODEProblem(sh_pte,[1.0,0.0],(zo,zi2),[fld,Ph0])
solg = solve(prob,Tsit5(),reltol=1e-1,abstol=1e-12,dtmax=1e-3)
g=z->solg(z)[1]
dth=ODEProblem((u,p,z)->eta*(fBa(0,z)+fBb(0,z))/(2*sqrt(Phr)),0.0,(zo,zi2))
chi = solve(dth,Tsit5(),reltol=1e-8,abstol=1e-10);
Spocitame vyvoj jednotlivych integrandu, je zrejme, efekt prvniho WF na I1 je kompenzovan druhym WF, I2 kompenzovan neni, ale neprojevi se na pozici v obraze
dI1 = z -> h(z)*Ph0*(fE1(0,z)+fE1b(0,z))/Phr^2*cos(chi(z))
dI2 = z -> g(z)*Ph0*(fE1(0,z)+fE1b(0,z))/Phr^2*cos(chi(z))
u0=0.0
prob = ODEProblem((u,p,z)->dI1(z),u0,zspan)
I1 = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12,dtmax=1e-3);
prob = ODEProblem((u,p,z)->dI2(z),u0,zspan)
I2 = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12,dtmax=1e-3);
z1 = zspan[1]:1e-4:zspan[2];
plot(z1,I1.(z1),label="I1")
plot!(z1,I2.(z1)/100,label="I2/100")
Kompenzace I1 je dany symetrickymi podminkami v systemu
p=Any[]
p1=plot(z1,h.(z1),label="h(z)")
plot!(p1,z1,(fE1.(0,z1).+fE1b.(0,z1))/1e5,label=L"\Phi_1/10^5")
push!(p,p1)
push!(p,plot(z1,h.(z1).*(fE1.(0,z1).+fE1b.(0,z1)),label=L"h(z)\Phi_1"))
plot(p...,layout=(1,2),size=(800,300))
Sila filtrace WF je pak dana difrakcni trajektorii v prvnim obrazu $z=z_{i1}$
D = I1(zi1);
print("D = I1(zi1) = ",D, ", tj. dx = MD kappa = ", g(zi1)*D/Ph0*1e6, " [mum/eV] dE\n")
D = I1(zi1) = 0.11670528662967299, tj. dx = MD kappa = -359.6832004983809 [mum/eV] dE
Bm = 0.2; sB = 1e-3; Ph0 = 10000;
fB = axfld(axc = (n,z)->dnB(n,z,0,1e-3,0.2),em=0)
fEd = axfld(axc=(n,z)->mpaf(n,z,-1e-2,5e-3,1e-3,1e3),m=1,em=1)
zspan = (-4e-2,2e-2)
th=(-1:0.1:1)*1e-4
p=plot()
for i=1:length(th)
u0 = [0.0,0.0,th[i],0.0,0.0]
prob = ODEProblem(g_pte,u0,zspan,[[fB,fEd],Ph0,0.0])
sol = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12,dtmax=1e-3);
plot!(p,sol,vars=(0,1),label=false)
end
plot(p)
plot!(fEd,0,0.5e-8,label=L"\Phi_1")
plot!(fB,0,0.4e-4,label=L"B")
ylabel!("[m]")
xlims!(zspan)
Mejme objektivovou cocku z predchoziho prikladu, ktera e navic mirne elipicka - je zde slabe magneticke kvadrupolove pole. Pro jednoduchost budeme predpokladat, ze ma shodny prubeh s osove symetrickym polem.
fB2 = axfld(axc = (n,z)->dnB(n,z,0,1e-3,0.02),em=0,m=2)
u0 = [0.0,0.0,1,0.0,0.0]
prob = ODEProblem(g_pte,u0,zspan,[[fB,fB2],Ph0,0.0])
solhx = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12,dtmax=1e-3);
hx = z->solhx(z)[1]
u0 = [0.0,0.0,0,1.0,0.0]
prob = ODEProblem(g_pte,u0,zspan,[[fB,fB2],Ph0,0.0])
solhy = solve(prob,Tsit5(),reltol=1e-10,abstol=1e-12,dtmax=1e-4);
hy = z->solhy(z)[2]
zix = find_zero(hx,[0.0,zspan[2]])
ziy = find_zero(hy,[0.0,zspan[2]])
dz = zix - ziy
print("Rozdil v poloze obrazu pro x a y je dz_i = ", dz, " m\n")
Rozdil v poloze obrazu pro x a y je dz_i = 2.7786136569169298e-5 m
z1=zspan[1]:1e-5:zspan[2]
p=Any[]
p1=plot(z1,hx.(z1))
plot!(p1,z1,hy.(z1))
push!(p,p1)
p2=plot(z1.-zix,hx.(z1))
plot!(p2,z1.-zix,hy.(z1))
xlims!(-1e-3,1e-3);ylims!(-1e-3,1e-3)
push!(p,p2)
plot(p...,layout=(1,2),size=(800,300))
Tento astigmatizmus OL je nutne vykompenzovat stigmatorem pred OL. Budeme postupovat metodou variace konstanty. Nejprve vyresime zhomogenizovanou rovnici paraxialni rovnici \begin{align} u^{\prime\prime}&+\frac{eB^2}{8m_e\Phi^*}u=0 \end{align}
#Vypocet trahektorie h
prob = ODEProblem(sh_pte,[0.0,1.0],zspan,[[fB],Ph0])
solh = solve(prob,Tsit5(),reltol=1e-9,abstol=1e-12,dtmax=1e-3)
h=z->solh(z)[1]
#Vypocet trajektorie g
prob = ODEProblem(sh_pte,[1.0,0.0],zspan,[[fB],Ph0])
solg = solve(prob,Tsit5(),reltol=1e-9,abstol=1e-12,dtmax=1e-3)
g=z->solg(z)[1]
#Vypocet rotace
dth=ODEProblem((u,p,z)->eta*fB(0,z)/(2*sqrt(Phr)),0.0,zspan)
chi = solve(dth,Tsit5(),reltol=1e-8,abstol=1e-10,dtmax=1e-3);
zi = find_zero(h,[0,zspan[2]])
M = g(zi); Ma = solh(zi)[2];
A spocitame trajektorii $u_{\alpha}$ z rovnice \begin{align} u_{\bar\alpha} = &-g\int\limits_{z_o}^{z}h^2\frac{\gamma_0}{\Phi^*}\left(\Phi_2+\im v_0\Psi_2-\frac{\Phi_1^2}{8\gamma_0\Phi^{*}}\right)e^{-2\im\chi(z)}\dd z\\ &+h\int\limits_{z_o}^{z}hg\frac{\gamma_0}{\Phi^*}\left(\Phi_2+\im v_0\Psi_2-\frac{\Phi_1^2}{8\gamma_0\Phi^{*}}\right)e^{-2\im\chi(z)}\dd z\nonumber \end{align} Bude pocitat pouze Image coefficient, pro nas sytem tedy $$C_{\bar\alpha}=-g\int\limits_{z_o}^{z}h^2\frac{\gamma_0(\Phi_2+\im v_0\Psi_2)}{\Phi^*}e^{-2\im\chi(z)}\dd z$$
Astigmatizmus samotne magneticke cocky vychazi
Phr = Ph0*(1+qe*Ph0/(2*me*c^2));
gamma = 1+qe*Ph0/(me*c^2) #rel. faktor gamma
v0 = 2*eta*sqrt(Phr)/gamma; #Rychlost elektronu
dA1OLr=ODEProblem((u,p,z)->-h(z)^2*gamma*v0*real(im*fB2(0,z)/Phr*exp(-2im*chi(z))),0.0,(zo,zi2))
#dA1OLr=ODEProblem((u,p,z)->-h(z)^2*2*eta/sqrt(Phr)*real(fB2(0,z)*exp(-2im*chi(z))),0.0,(zo,zi2))
A1OLr = solve(dA1OLr,Tsit5(),reltol=1e-8,abstol=1e-10,dtmax=1e-3);
dA1OLi=ODEProblem((u,p,z)->-h(z)^2*gamma*v0*imag(im*fB2(0,z)/Phr*exp(-2im*chi(z))),0.0,(zo,zi2))
A1OLi = solve(dA1OLi,Tsit5(),reltol=1e-8,abstol=1e-10,dtmax=1e-3);
print("A1OL = ", A1OLr(zi)+im*A1OLi(zi),"\n")
z1 = zspan[1]:1e-4:zspan[2]
plot(z1,(A1OLr.(z1)))
plot!(z1,(A1OLi.(z1)))
A1OL = -0.0001898059607617512 + 8.484786094122903e-6im
Ted do systemu pridame jeden stigmator - elektrostaticky kvadrupol a zjistime jaky ma efekt na vysledny astigmatizmus. A urcime jeho skalovy faktor tak aby byl celkovy astigmatizmus nulovy
fE2 = axfld(axc = (n,z)->mpaf(n,z,-0.01,5e-3,1e-3,10e3),m=2,em=1) #
dA1s=ODEProblem((u,p,z)->-h(z)^2*gamma*real(fE2(0,z)/Phr*exp(-2im*chi(z))),0.0,(zo,zi2))
A1s = solve(dA1s,Tsit5(),reltol=1e-8,abstol=1e-10,dtmax=1e-3);
scS = -(A1OLr(zi)+im*A1OLi(zi))/A1s(zi)
fE2s=axfld(axc = (n,z)->scS*mpaf(n,z,-0.01,5e-3,1e-3,10e3),m=2,em=1);
dA1OLsr=ODEProblem((u,p,z)->-h(z)^2*gamma*real((fE2s(0,z)+im*v0*fB2(0,z))/Phr*exp(-2im*chi(z))),0.0,(zo,zi2))
A1OLsr = solve(dA1OLsr,Tsit5(),reltol=1e-8,abstol=1e-10,dtmax=1e-3);
dA1OLsi=ODEProblem((u,p,z)->-h(z)^2*gamma*imag((fE2s(0,z)+im*v0*fB2(0,z))/Phr*exp(-2im*chi(z))),0.0,(zo,zi2))
A1OLsi = solve(dA1OLsi,Tsit5(),reltol=1e-8,abstol=1e-10,dtmax=1e-3);
plot(z1,A1OLsr.(z1),label=L"R(A_{1})")
plot!(z1,A1OLsi.(z1),label=L"Im(A_{1})")