Podklady k šesté přednášce - Obecná paraxiální rovnice

In [1]:
using Plots #Balicek pro grafy
using DifferentialEquations # Balicek pro vypocet diferencialnich rovnic
using SpecialFunctions # Specialni funkce
using PolyChaos #Balicek pro vypocet hermiteovych polynomu
using Roots # Balicek pro nalezeni korenu funkce
using LaTeXStrings # Latex v grafech
using RecipesBase
include("/home/radlicka/Juliawork/ElectronOptics/ParaxialEquation.jl")
Out[1]:
ImgDevTrans
In [2]:
include("/home/radlicka/Juliawork/ElectronOptics/ParaxialEquation.jl")
Out[2]:
ImgDevTrans
In [3]:
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
;

Pole v systému

Pole magnetické čočky

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})$$

Multipolova pole

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) $$

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,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))
Out[4]:
#37 (generic function with 1 method)

Osove symetricka magneticka cocka

In [5]:
fB1 = axfld(axc = (n,z)->dnB(n,z,0,1e-3,0.2),em=0)
Ph0=10000; kappa=0; 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]
    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)
plot!(fB1,0,1.1*4e-6/0.2,color=:red,label="B")
ylabel!(L"x_r [ \mathrm{m} ]")
xlims!(zspan[1],zspan[2])
Out[5]:

Wienuv filtr pro energiovou fitraci

Do system umistime dve magneticke cocky. První do $z = -0.8\,\rm m$, druhou do $z=0$. Budem predpokladat polohu predmetu $z_o=-0.9\,\rm m$ a prvni cocku nastavime tam aby zobrazovala predmet do $z_{i,1}=0.45$. Druhou cocku pak nastavime tak, aby zobrazovala $z_{i,1}$ do roviny obrazu $z_i = 0.01\,\rm m$.

Wienuv filtr umistime do roviny $z=-0.7\,\rm mm$ a budeme sledovat, jak ovlivňuje trajektorie častic

In [6]:
zL1 = -0.8; zL2=0 #Poloha prvni a druhe magneticke cocky
zW1 = -0.7; LW=0.15 #Poloha a delka Wiennova filtru
zo = -0.9; zi1=-0.45; zi2=0.01 #Poloha predmetu, prvn9ho a druh0ho obrazu
Ph0 = 1e3; 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
Ph1m = 1e4; Ps1m = Ph1m/v0; # Maximalni hodnoty pole ve WF
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 = 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
plot(fBa,label="B [T]")
plot!(fBb,label="B [T]")
plot!(fE1,0,1e-5,label=L"\Phi_1\cdot 10^{-5}\,[\mathrm{V/m}]")
plot!(fB1,0,-1im*200,label=L"\Psi_1\cdot 200\,[\mathrm{T}]")
plot!(fE2,0,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!(" ")
Out[6]:

1. Nastaveni magnetickych cocek

$z_o = -0.05\,\rm m$, $z_{i1} = -0.15$

Řešíme zhomogenizovanou rovnici pro stigmaticky system se sylný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}

In [7]:
# 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];
# 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

th=(-1:0.2:1)*1e-3
plot(z1,h1.(z1)*th',label=false,color=:blue)
plot!(z2,h2.(z2)*th'*Ma1,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)
Out[7]:

Tyto trajeketorie platí pro ideální energii elektronu, v případě že $\kappa\neq 0$ dostaneme

In [8]:
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);
    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
In [9]:
plot(p)
Out[9]:

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.

In [10]:
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
In [11]:
plot(px)
Out[11]:
In [12]:
plot(px)
xlims!(9.99e-3,10.01e-3); ylims!(-0.2e-6,0.2e-6)
Out[12]:

Zbytkovou energiovou disperzy je nutne vyrusit druhym WF

In [13]:
#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())
Out[13]:
0.30798817659008487
In [14]:
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
In [15]:
plot(px,xlabel="z [m]",ylabel="x [m]")
Out[15]:
In [16]:
plot(px)
xlims!(9.995e-3,10.008e-3); ylims!(-0.1e-6,0.1e-6)
Out[16]:

Pristup pomoci metody variace konstanty

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}

In [17]:
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

In [18]:
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")
Out[18]:

Kompenzace I1 je dany symetrickymi podminkami v systemu

In [19]:
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))
Out[19]:

Sila filtrace WF je pak dana difrakcni trajektorii v prvnim obrazu $z=z_{i1}$

In [20]:
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.11670528662804279, tj. dx = MD kappa = -359.6832004935424 [mum/eV] dE

Slabe dipolove pole (neoptimalizovany deflektor)

In [21]:
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)
Out[21]:

Osovy astigmatizmus objektivove cocky

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.

In [22]:
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.7786136807542855e-5 m
In [23]:
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))
Out[23]:

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}

In [24]:
#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

In [25]:
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.0001898059607715529 + 8.484786158435653e-6im
Out[25]:

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

In [26]:
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);
In [27]:
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})")
Out[27]: