Simulation of logarithmic lens

written by Nan Zheng
Duke University

Included as figure 10.5 in Optical Imaging and Spectroscopy, David J. Brady, www.opticalimaging.org

Extended Nijboer-Zernike method

According to reference [1], the phase modulation of the logarithmic lens is
ϕ(ρ)=(π/λ)(((f^2)+(ρ^2)(a^2))/((Subscript[s, 1]) ln(1+((a^2)/(f^2)))))[ln((((f^2)+(ρ^2)(a^2))/((f^2)+(a^2))))-1]-((2π)/λ)((Sqrt[(ρ^2)(a^2)+(f^2)])-f)
=Φ(ρ)-((2π)/λ)(((ρ^2)(a^2))/(2f))                                       for f≫ρa≫λ
where Φ(ρ)=(π/λ)(((f^2)+(ρ^2)(a^2))/((Subscript[s, 1]) ln(1+((a^2)/(f^2)))))[ln((((f^2)+(ρ^2)(a^2))/((f^2)+(a^2))))-1]
Here λ is the wavelength of the light, f is the expected focal length, a is the radius of the lens, and (Subscript[s, 1]) is the lowest distance that can be imaged on f plane by certain ring of the lens (the center). The longest distance that can be imaged is ∞.
Then the amplitude distribuation in the focal region is given by
U((Subscript[x, 2]), (Subscript[y, 2]), (Subscript[z, 2]))∝((Subscript[U, 0])/(λz))∫∫P(ρ)(e^(-iΦ(ρ)))(e^(i((2π)/λ)(((ρ^2)(a^2))/(2f))))(e^(((-ik)/(2z))((ρ^2)(a^2)-2(Subscript[x, 2])x-2(Subscript[y, 2])y)))dx dy                                 where z=f+(Subscript[z, 2])
=((Subscript[U, 0])/(λz))∫∫P(ρ)(e^(-iΦ(ρ)))(e^(i(π/λ)((1/f)-(1/z))(ρ^2)(a^2)))(e^(((-ik)/z)aρr cos(θ-φ)))ρdρdθ
=((Subscript[U, 0])/(λz))(∫_0^1)(e^(-iΦ(ρ)))(e^(i(π/λ)((1/f)-(1/z))(a^2)(ρ^2)))(Subscript[J, 0])(((k a r)/z)ρ)ρdρ
=((Subscript[U, 0])/(λz))(∑_(n=0)^∞)(Subscript[α, 2n])(∫_0^1)(Subscript[R, 2n]^0)(ρ)(e^(i(π/λ)((1/f)-(1/z))(a^2)(ρ^2)))(Subscript[J, 0])(((k a r)/z)ρ)ρdρ
                                    where (Subscript[α, n])=2(n+1)(∫_0^1)(e^(-iΦ(ρ)))(Subscript[R, n]^0)(ρ)ρdρ    and    (e^(-iΦ(ρ)))=(∑_(n=0)^∞)(Subscript[α, 2n])(Subscript[R, 2n]^0)(ρ)
=((Subscript[U, 0])/(λz))(∑_(n=0)^∞)(Subscript[α, 2n])(∫_0^1)(Subscript[R, 2n]^0)(ρ)(e^(i u (ρ^2)))(Subscript[J, 0])(v ρ)ρdρ
                                    where u=(π/λ)((1/f)-(1/z))(a^2)=(π/λ)((1/f)-(1/(f+(Subscript[z, 2]))))(a^2), v=((2π)/λ)(ar/z)=((2π)/λ)((a r)/(f+(Subscript[z, 2]))), (e^(-iΦ(ρ)))=(∑_(n=0)^∞)(Subscript[α, 2n])(Subscript[R, 2n]^0)(ρ)
∝((Subscript[U, 0])/(λz))(∑_(n=0)^∞)(Subscript[α, 2n])(∑_(L=1)^∞)-(((2i u))^(L-1))(∑_(j=0)^n)(Subscript[v, Lj])(((Subscript[J, L+2j])(v))/(L (v^L)))                             according to reference [2]
Here (Subscript[v, Lj])=(((-1))^n)(L+2j)(((({{j+L-1}, {L-1}})))^2)(({{L-1}, {n-j}}))/(({{n+L+j}, {L}}))
(((Subscript[J, L+2j])(v))/(L (v^L)))(Overscript[→, v->0])(1/(L (v^L)))(((((v/2)))^(L+2j))/(Γ(L+2j+1)))=((v^(2j))/(L  Γ(L+2j+1) (2^(L+2j))))=((v^(2j))/(L (L+2j)! (2^(L+2j))))={({{0                                                        j≠0}, {(1/(L (L+2j)! (2^(L+2j)) ))                            j=0}})

parameters

a=8*(10^7);λ=5*(10^3);f=32*(10^7);s1=2*(10^10);

Program for the logarithmic lens based on the extended Nijboer-Zernike method

comb[n_,m_]:=((n!)/((m!)*((n-m)!)))

v[n_,l_,j_]:=(((-1))^n) *(l+2j)*comb[j+l-1,l-1]*comb[j+l-1,l-1]*comb[l-1,n-j]/comb[n+l+j,l]

R[n_,m_,ρ_]:=Sum[(((((-1))^s) (ρ^(n-2 s)) (n-s)!)/(((1/2) (-Abs[m]+n)-s)! (((Abs[m]+n)/2)-s)! s!)),{s,0,(1/2) (-Abs[m]+n)}]

alpha[n_]:=alpha[n]=2(1+n) NIntegrate[(E^(-I (π/λ) (((f^2)+(ρ^2)*(a^2))/(s1* Log[1+((a^2)/(f^2))])) (Log[(((f^2)+(ρ^2)*(a^2))/((f^2)+(a^2)))]-1))) R[n,0,ρ] ρ,{ρ,0,1}]

ilog[r_,z_,f_,lamda_,a_,nmax_]:=Abs[(1/(f+z))Sum[alpha[2n]*If[z==0,Sum[Sum[v[n,l,j] If[r==0,If[j==0,(1/(l*(l!)*(2^l))),0],((BesselJ[l+2j,((2π)/lamda) ((a r)/(f+z))])/(l (((((2π)/lamda) ((a r)/(f+z))))^l)))],{j,0,n}],{l,1,1}],Sum[(((-2I (π/lamda) ((1/f)-(1/(f+z))) (a^2)))^(l-1))*Sum[v[n,l,j] If[r==0,If[j==0,(1/(l*(l!)*(2^l))),0],((BesselJ[l+2j,((2π)/lamda) ((a r)/(f+z))])/(l (((((2π)/lamda) ((a r)/(f+z))))^l)))],{j,0,n}],{l,1,25}]],{n,0,nmax}]]^2*(1/(1.7569154423259176`*^-20))

Timing[Plot[ilog[r*(10^4),0,32*(10^7),5*(10^3),8*(10^7),14],{r,0,18},PlotRange->{0,1}]]

Timing[Plot[ilog[r*(10^4),5*(10^3),32*(10^7),5*(10^3),8*(10^7),14],{r,0,18},PlotRange->{0,1}]]

Timing[Plot[ilog[0,z*5*(10^3),32*(10^7),5*(10^3),8*(10^7),14],{z,-10,10},PlotRange->All]]

Numerical integration

parameters

a=8*(10^7);λ=5*(10^3);f=32*(10^7);s1=2*(10^10);

Program for finite object distance s

idirect[ρ2_,s_,z_]:=Abs[NIntegrate[Evaluate[ρ*(1/s)*(E^(-I*((2π)/λ)*(Sqrt[(ρ^2)+(s^2)])-I (π/λ) (((f^2)+(ρ^2))/(s1* Log[1+((a^2)/(f^2))])) (Log[(((f^2)+(ρ^2))/((f^2)+(a^2)))]-1)+I*((2π)/λ)*(Sqrt[(ρ^2)+(f^2)])-I*((2π)/λ)*(Sqrt[(ρ^2)+(((f+z))^2)])))*BesselJ[0,2π((ρ*ρ2)/(λ*(Sqrt[(ρ^2)+(((f+z))^2)])))]/.s1->2*(10^10)],{ρ,0,a}]]^2(1/(1.9987564753772122`*^8))

Timing[Plot[idirect[0,2*(10^10),z*5*(10^3)],{z,-200,1500},PlotRange->All]]

Timing[Plot[idirect[ρ2*(10^4),2*(10^10),0],{ρ2,0,18},PlotRange->All]]

Program for infinite object distance

idirectinfi[ρ2_,z_]:=Abs[NIntegrate[Evaluate[ρ*(E^(-I (π/λ) (((f^2)+(ρ^2))/(s1* Log[1+((a^2)/(f^2))])) (Log[(((f^2)+(ρ^2))/((f^2)+(a^2)))]-1)+I*((2π)/λ)*(Sqrt[(ρ^2)+(f^2)])-I*((2π)/λ)*(Sqrt[(ρ^2)+(((f+z))^2)])))*BesselJ[0,2π((ρ*ρ2)/(λ*(Sqrt[(ρ^2)+(((f+z))^2)])))]/.s1->2*(10^10)],{ρ,0,a}]]^2(1/(7.369037467416818`*^28))

Timing[Plot[idirectinfi[0,z*5*(10^3)],{z,-1250,400},PlotRange->All,AxesLabel->{"z (λ)","Intensity"}]]

Timing[idirectinfi[0,0]]

Timing[Plot[idirectinfi[ρ2*(10^4),0],{ρ2,0,18},PlotRange->All]]

Timing[DensityPlot[idirectinfi[ρ2*5*(10^3),z*5*(10^3)],{z,-1250,400},{ρ2,-72,72},PlotPoints->200,PlotRange->All,Mesh->False,AxesLabel->{"z (λ)","ρ (λ)"}]]

Timing[Plot[idirectinfi[ρ2*5*(10^3),100*5*(10^3)],{ρ2,0,100}]]

Timing[Plot[idirectinfi[ρ2*5*(10^3),100*5*(10^3)],{ρ2,0,150},PlotRange->All]]

Timing[Plot[idirectinfi[ρ2*5*(10^3),100*5*(10^3)],{ρ2,-500,500}]]