Thin Lens - Fourier Optics     H. Q. Le  copyrighted

ECE 5368/6358                              han q le - copyrighted
Use solely for students registered for UH ECE 6358/5368 during courses - DO NOT distributed (copyrighted materials)

Segment 1

0. Overview & summary

Key points of this chapter.

Lens as a phase transformation object: it alters the phase of the input field that varies quadratically as a function of radial position relative to the lens axis.

A key effect of this phase transformation is that it can negate the quadratic Fresnel diffraction phase term.

At the plane where the two quadratic phase terms cancel each other out exactly, the only thing left is the Fraunhofer diffraction term, which is a spatial Fourier transform of the input field amplitude with the spatial frequency being cutoff by the finite aperture of the lens. This plane position is known as the focal plane.

The ability of the lens to form image is from the fact that the field from the object to the lens is a Fourier transform by Fraunhofer diffraction, and after the lens, with the quadratic Fresnel phase term removed, it is a  Fraunhofer diffraction Fourier transform again. Double Fourier transforms reproduce the original object amplitude which is an image.
Note: The well-known geometrical optics lens laws come from the thin lens wave theory:
                                  (1/(Subscript[z, 1]))+(1/(Subscript[z, 2]))= (1/f)  and image mapping:  (u/(Subscript[z, 1]))=-(v/(Subscript[z, 2]))

However, because of the finite lens aperture, the second Fourier transform is NOT complete, but cutoff at higher spatial frequency. Hence the resulted image cannot be infinitely sharp like the original object, giving rise to the finite ability of a lens to give sharp image: this is the concept of lens resolution.
The lower the f-number or the higher the lens numerical aperture, the higher is the lens resolution.

1. Theory basic

In a nut shell...

1.1 Lens as a non-uniform phase shifter

Approximately, what a lens does is altering the amplitude and phase of the input wave by an aperture term (its size) and a phase term that is a function of space owing to its shape - see below.

Given an input wave, the output wave can be expressed as:
                          (Subscript[U, out])[x,y]=(Subscript[U, in])[x,y] P[x,y](e^(i φ[x,y]))                         (1.1.1)
where P[x,y] is the amplitude transmission function of the lens, which can be just the aperture, and (e^(i φ[x,y])) is the phase shift that is a function of position.

A lens typically has two “spheroid” surfaces. We use the term spheroid instead of spherical to indicate that in general, lens can be designed with advanced aspheric shape that is actually better for imaging with regard to wavefront aberration. But for many conventional lenses, the manufacturing method of grinding-polishing a piece of dielectric such as glasses, semiconductors, or other hard oxides (e. g. sapphire) necessarily results in a spherical surface. Molding, especially for polymers can produce aspheric lenses. But we are not into lens manufacturing technology in this chapter. It is sufficient in this chapter to use the spherical model to illustrate the lens concept, which is also consistent with the paraxial approximation.

For thin lens with two spherical surfaces, the phase shift is (in the paraxial approximation):
                φ[x,y]=k n D-k ((ρ^2)/2) (n-1)((1/(Subscript[R, 1]))-(1/(Subscript[R, 2])))=k n D-k ((ρ^2)/(2 f))      (1.1.2)
where                                (n-1)((1/(Subscript[R, 1]))-(1/(Subscript[R, 2])))=(1/f)                                   (1.1.3)

Side note 1: Strictly speaking,  (1.1.1) is not the rigorous EM wave solution. One has to solve the complete problem that includes the entire lens as anembedded structure to really get the correct solution of the field after the lens, but (1.1.1) is a very good approximation for lens with dimension >> λ. Microlenses with dimensions ~ few times λ, like Mie scattering will require more rigorous solution. One can also ask what happens to the reflected field? We will always assume here that the surface is AR-coated.

Side note 2: in the expression (Subscript[U, in])[x,y] for the input E-field, where is the expression for (Subscript[k, z]) component?

Given (1.1.1) and (1.1.3), what happens after the field gets through a lens?
We study the following special cases:

1.2 Object tranmission function and field at focal plane

If we neglect the uniform phase term k n D, A beam after thin lens is approximated as:   
                                (Subscript[U, out])[x,y]=(Subscript[U, in])[x,y] P[x,y] (e^(-i (k/(2 f)) (ρ^2) ))                  (1.2.1)

1.2.1 First - without the lens

To find out what the field looks like at a location z, we uses the Fresnel integral (see link) for normal incident plane wave with unit amplitude: (at a distance away).
            U[x,y] = k ((e^(i k z))/(i 2 π z)) ∫T[ξ,η] (e^(i (k/(2 z))((((Subscript[x, o])-ρ) )^2)))dξ dη          (1.2.2a)
            U[x,y] = k ((e^(i k z))/(i 2 π z)) ∫T[ξ,η] (e^(i (k/(2 z))( (((x-ξ))^2)+(((y-η))^2)) ))dξ dη          (1.2.2b)
Here, we assume T[ξ, η] is some aperture function illuminated with a plane wave (see link above).

We do a little algebra:

Expand[(((x-ξ))^2)+(((y-η))^2)]

Hence         (e^(i (k/(2 z))( (((x-ξ))^2)+(((y-η))^2)) )) =(e^(i (k/(2 z))(ρ^2)))(e^(i (k/(2 z))((x^2)+(y^2)))) (e^(-i (k/z)( ξ x+ y η)))              (1.2.3)
where                                 (ρ^2) =(ξ^2) +(η^2)                                     (1.2.3a)
Hence, we can write
               U[x,y] = k ((e^(i k z))/(i 2 π z)) (e^(i (k/(2 z))((x^2)+(y^2)) ))∫T[ξ,η] (e^(i (k/(2 z))((ξ^2)+(η^2)) ))(e^(-i (k/z)(ξ x+η y) ))dξ dη     (1.2.4a) 
Or:         U[x,y] = k ((e^(i k z))/(i 2 π z)) (e^(i (k/(2 z))((x^2)+(y^2)) ))∫T[ξ,η] (e^(i (k/(2 z))(ρ^2) ))(e^(-i (k/z)(ξ x+η y) ))dξ dη     (1.2.4b)
This is without the lens.
(see this:
        U[(Subscript[x, o])]=-A((i k)/(2π)) ((e^(i k |(Subscript[x, o])| ))/(|(Subscript[x, o])|))Cos[α](∮_A) (T[r'] (e^(i (k/(2 |(Subscript[x, o])|))((r')^2))) )(e^(-i(k/(|(Subscript[x, o])|)) (Subscript[ρ', o]). r'))  dS    (8.2.3.19a)         )
It has: A Fourier transform term:  (e^(-i (k/z)(ξ x+η y) ))  (recall Fraunhofer diffraction)
and an axial quadratic term:       (e^(i (k/(2 z))(ρ^2) )), which is the Fresnel term.

At this stage, we don’t have the image T[ξ,η] any more. It is blurred quite a bit because of this integration with both terms.

1.2.2 With the lens

With the lens: 
    (Subscript[U, out])[x,y] =  k ((e^(i k z))/(i 2 π z)) (e^(i (k/(2 z))((x^2)+(y^2)) ))∫T[ξ,η] P[ξ,η] (e^(-i (k/(2 f)) (ρ^2) )) (e^(i (k/(2 z))(ρ^2) ))(e^(-i (k/z)(ξ x+η y) ))dξ dη    (1.2.5a)
Here, the target transmission function T[ξ,η] is now multiplied with the lens aperture function P[ξ,η]  and of course, it has the lens phase term (e^(-i (k/(2 f)) (ρ^2) ))  that is the essence here (the phase term is far more important - but we’ll see that the aperture is essential to image sharpness - i. e. resolution)

Hence
   U[x,y] =k ((e^(i k z))/(i 2 π z))  (e^(i (k/(2 z))((x^2)+(y^2))))∫T[ξ,η] P[ξ,η] (e^(i (k/2) (ρ^2) ((1/z)-(1/f))))(e^(-i (k/z)( ξ x+ y η)))dξ dη     (1.2.5b)
Notice something interesting here:   (e^(i (k/(2 z))(ρ^2) )) which is the Fresnel quadratic phase term, represents the key diffraction term that make features of an aperture becomes blurred; i. e. a sharp edge becomes soft with ripple and fringes because of this term. However, the lens offers a term that seems to be "anti-diffraction" (e^(-i (k/(2 f)) (ρ^2) )). At the right location, the two terms can cancel each other out.

Discussion: The essence of the lens behavior is that with quadratic phase term, it can counter the Fresnel quadratic term in paraxial approximation, hence, straighting out the phase front to create an image (or Fourier transform in the focal plane). We will see later that by creating a reverse wavefront, it creates an "reverse diffraction" effect to make the wave converges to an image. That is very much the essence of lens imaging.

1.2.3 Spatial Fourier transform:

Notice that (1.2.5b) can be written:
   U[x,y] =k ((e^(i k z))/(i 2 π z))  (e^(i (k/(2 z))((x^2)+(y^2))))∫(T[ξ,η] P[ξ,η] (e^(i (k/2) (ρ^2) ((1/z)-(1/f)))))(e^(-i ( ξ  (x(k/z))+  η (y(k/z)))))dξ dη     (1.2.5b)
which can be expressed in terms of FT, which is defined generally:
                    G[(Subscript[k, x]),(Subscript[k, y])]=(F^(-1))[g[ξ,η] ]=(1/(2 π))∫g[ξ,η] (e^(-i ((Subscript[k, x]) ξ + (Subscript[k, y]) η)))dξ dη             (1.2.6)
Hence:
           U[x,y] =k ((e^(i k z))/(i z))  (e^(i (k/(2 z))((x^2)+(y^2))))(F^(-1))[T[ξ,η] P[ξ,η] (e^(i (k/2) (ρ^2) ((1/z)-(1/f))))][k(x/z),k(y/z)]    (1.2.7)

At focal plane:   (1/z)-(1/f)=0, hence:
U[x,y](Subscript[|, z=f]) =k ((e^(i k z))/(i z))  (e^(i (k/(2 z))((x^2)+(y^2))))(F^(-1))[T[ξ,η] P[ξ,η]][k(x/z),k(y/z)]      (1.2.8a)
                     =k ((e^(i k f))/(i f)) (e^(i (k/(2 f))((x^2)+(y^2)))) (F^(-1))[T P][k(x/f),k(y/f)]                   (1.2.8b)

This is a central result of lens imaging and Fourier optics: The focal plane field of the lens is approx the Fourier transform of the input on the lens (but with an additional quadratic phase term (e^(i (k/(2 f))((x^2)+(y^2))))). It is also a representative of Fraunhofer diffraction of the object and lens aperture T[ξ, η] P[ξ, η] in the far field. In other words, the FF diffraction pattern is displayed on the focal plane. This has been used extensively for analog optical processing.
We will see later that if the input image is at the input focal plane, this phase term (e^(i (k/(2 f))((x^2)+(y^2)))) will also disappear.

Fourier transform of an image is NOT an image, but can be transformed back. This is an essential elements of digital image processing. Although what we discuss here is fundamental wave optics phenomenon.

Example

Fresnel=ImageData[ ];
{Fheight,Fwidth}=Take[Dimensions[Fresnel],2]

FresnelBW=Reverse[Table[Mean[Fresnel[[i,j]]],{i,1,Dimensions[Fresnel][[1]]},{j,1,Dimensions[Fresnel][[2]]} ] ];
FresnelBW=Transpose[Drop[Transpose[Drop[FresnelBW,1]],1] ];
FresnelEBW=Sqrt[FresnelBW];
Grid[{{ListDensityPlot[FresnelBW, ColorFunction->GrayLevel, AspectRatio->(Fheight/Fwidth)]
,ListDensityPlot[FresnelEBW, ColorFunction->GrayLevel, AspectRatio->(Fheight/Fwidth)]}}]

{ic,jc}=Ceiling[{Fheight,Fwidth}/2.]

FTFresnel=Fourier[FresnelEBW] ;
FTF=RotateRight[FTFresnel,ic];
FTF=Transpose[FTF];
FTF=RotateRight[FTF,jc];
FTF=Transpose[FTF];
ListDensityPlot[Log[Abs[FTF]], ColorFunction->GrayLevel, PlotRange-> All]

ListDensityPlot[Abs[FTF]^2, ColorFunction->GrayLevel, PlotRange-> All]

ListPlot3D[Log[Abs[FTF]], ColorFunction->"Rainbow", PlotRange-> All, BoxRatios->{1,1,0.7}]

Exercise: If we look at the focal plane, what does it really look like? Do this exercise with the following parameters: pixel size of the image = 50 μm; lens focal length f=5 cm; wavelength λ=0.55 μm (typical green); plot the light intensity at the focal plane (use real physical dimension).

Solution: (See HW discussion)

A little note: by convention, the expression:
                             (1/(2 π))∫∫T[x,y] (e^(-i ((Subscript[q, x])x+(Subscript[q, y])y))) dx dy = R[(Subscript[q, x]),(Subscript[q, y])]                       (1.2.9)
is usually called inverse Fourier transform:
                            (F^(-1))[T[x,y][(Subscript[q, x]),(Subscript[q, y])]=(1/(2 π))∫∫T[x,y] (e^(-i ((Subscript[q, x])x+(Subscript[q, y])y))) dx dy            (1.2.10)
However, the distinction between  (F^(-1)) and (F^0) here is NOT crucial as long as we keep the sign of the exponent in the transforms consistently. We can use the same expression FT for both, or even drop the inverse exponent -1 (as long as we remember which and which). The important point is to keep the convention that the spatial image is expressed as:
                                   T[x,y]=(1/(2 π))∫∫ R[(Subscript[q, x]),(Subscript[q, y])] (e^(i ((Subscript[q, x])x+(Subscript[q, y])y))) dx dy                      (1.2.11)                                    
Note also that this is Fraunhofer diffraction of an intensity pattern T[ξ,η] P[ξ,η]. Thus, in the focal plane, it is similar to the far-field Fraunhofer diffraction pattern. However, we also note that there is a quadratic phase term  (e^(i (k/(2 f))((x^2)+(y^2)))) remaining. This indicates it is NOT a flat wavefront, but is curved and can affect its subsequent propagation. However, in many cases, it is small, negligible, or corrected with real optics as the real lens is NOT just (e^(-i (k/(2 f))((x^2)+(y^2)))), but with higher order terms to flatten wavefront.

1.2.4 Spatial frequency discussion

Consider further from (1.2.8b) that (( k)/f)x,(( k)/f)y have the unit of spatial frequency. We can also rewrite (1.2.8b) as:
            U[x,y](Subscript[|, z=f])=k ((e^(i k f))/(i f)) (e^(i (k/(2 f))((x^2)+(y^2)))) (F^(-1))[T P][k(x/f),k(y/f)]     
We introduce a transformation of variables: 
                                            k(ξ/f)  =-(Subscript[q, x])  ;     k(η/f)  =-(Subscript[q, y])                                    (1.2.12a)
     U[x,y](Subscript[|, z=f]) = k ((e^(i k z))/(i 2 π z))  (e^(i (k/(2 z))((x^2)+(y^2))))∫T[ξ,η]P[ξ,η]  (e^(-i (k/f)( ξ x+ y η)))dξ dη  
         = k ((e^(i k f))/(i 2 π f)) (e^(i (k/(2 f))((x^2)+(y^2)))) ((f^2)/(k^2))∫T[-(Subscript[q, x])(f/k),-(Subscript[q, y])(f/k)]P[-(Subscript[q, x])(f/k),-(Subscript[q, y])(f/k)]  (e^(i ( (Subscript[q, x]) x+(Subscript[q, y]) y )))d(Subscript[q, x]) d(Subscript[q, y])    (1.2.12b)  

Note a key point here: the integration is no longer over the real space domain dξ dη, but over the spatial frequency domain: d(Subscript[q, x]) d(Subscript[q, y])
We recognize the FT expression:
          U[x,y](Subscript[|, z=f]) = k ((e^(i k f))/(i 2 π f)) (e^(i (k/(2 f))((x^2)+(y^2)))) 2 π((f^2)/(k^2))F[T[-(Subscript[q, x])(f/k),-(Subscript[q, y])(f/k)]P[-(Subscript[q, x])(f/k),-(Subscript[q, y])(f/k)]]      
          U[x,y](Subscript[|, z=f]) = (f/(i k)) (e^(i k f)) (e^(i (k/(2 f))((x^2)+(y^2)))) F[T[-(Subscript[q, x])(f/k),-(Subscript[q, y])(f/k)]P[-(Subscript[q, x])(f/k),-(Subscript[q, y])(f/k)]]       (1.2.13a)
Hence, use the inverse FT:             
       ((i k)/f) (e^(-i k f)) (F^(-1))[(e^(-i (k/(2 f))((x^2)+(y^2))))U[x,y]]=   T[-(Subscript[q, x])(f/k),-(Subscript[q, y])(f/k)]P[-(Subscript[q, x])(f/k),-(Subscript[q, y])(f/k)]   (1.2.13b)
In other words, the spatial frequency representation (inverse FT) of:
                                  (OverTilde[U])[x,y]≡U[x,y]  (e^(-i (k/(2 f))((x^2)+(y^2))))                           (1.2.13c)
is simply a product:
                                       T[-(Subscript[q, x])(f/k),-(Subscript[q, y])(f/k)]P[-(Subscript[q, x])(f/k),-(Subscript[q, y])(f/k)]                  (1.2.13d)
with a constant factor. This will be seen to be generally true for lens imaging:

The Fourier transform of an image (multiplied by a phase and a propagation amplitude factor) is simply the product of a k-space function associated with the object and an "optical transfer function" of the lens. 
We will see later that the concept of "optical transfer function" has a general definition applied to any imaging system beyond just for a lens.

It is important to recognize that functions T, P are originally real-space domain functions, but here they are converted to the frequency domain via the transformation of variables. This is not the same as the FT of T, P which are also in the frequency domain.
If the lens is just a circular aperture with unit amplitude, then, (1.2.13d) is:
                                         T[-(Subscript[q, x])(f/k),-(Subscript[q, y])(f/k)]Θ[((a )/f)k-q]                                        (1.2.14)
where Θ[((a )/f)k-q] is the unit step (Heaviside) function that represents a circle in {(Subscript[q, x]),(Subscript[q, y])} space with radius:
                                                  (Subscript[Q, 0])=k(a/f)=2 π(( a)/(λ f))                                     (1.2.14a)
We see that this means all high spatial frequency components of U[x,y] (e^(-i (k/(2 f))((x^2)+(y^2)))) are zero, it means that the image cannot contains higher spatial frequency than given in (1.2.14a). Thus, one must conclude that it is the lens that cuts off all the high spatial frequency component of the object transmission function. In other words, the lens acts as a hard k-space filter with abrupt cut-off.

Exercise 1   Amplitude and phase of a uniform beam (plane wave)

Let's look at the intensity and phase of a uniform input beam:  
      U[x,y] = k ((e^(i k z))/(i 2 π z)) (e^(i (k/(2 z))((x^2)+(y^2))))∫(Subscript[U, in])[ξ,η] P[ξ,η] (e^(i (k/2) (ρ^2) ((1/z)-(1/f))))(e^(-i (k/z)( ξ x+ y η)))dξ dη  (1.2.5)
with                                                            (Subscript[U, in])[ξ,η]=1                                    (1.2.6)    
and P[ξ, η] is just the aperture of a circle:
            U[x,y] = k ((e^(i k z))/(i 2 π z)) (e^(i (k/(2 z))((x^2)+(y^2))))(∫_0^(2 π))(∫_0^(D/2))(e^(i (k/2) (ρ^2) ((1/z)-(1/f))))(e^(-i (k/z)ρ r Cos[ϕ] ))ρdρdϕ  
            U[x,y] = k ((e^(i k z))/(i 2 π z)) (e^(i (k/(2 z))(r^2)))(∫_0^(D/2))(e^(i (k/2) (ρ^2) ((1/z)-(1/f))))2 π (BesselJ[0, ((k r ρ)/z)])ρdρ   
            U[x,y] = k ((e^(i k z))/(i  z)) (e^(i (k/(2 z))(r^2)))(∫_0^(D/2))(e^(i (k/2) (ρ^2) ((1/z)-(1/f))))(BesselJ[0, ((k r ρ)/z)])ρdρ

UniformBeam[z_,r_,λ_,f_,Di_]:=Module[{k},
k=2 π/λ ; ((E^(I k (Sqrt[(z^2)+(r^2)])))/(I  (Sqrt[(z^2)+(r^2)]))) *NIntegrate[(E^(I (k/2) (ρ^2)  ((1/z)-(1/f)))) BesselJ[0, ((k r ρ)/z)] ρ ,{ρ,0., Di/2.}] ];

Plot3D[Re[UniformBeam[z,(Sqrt[(x^2)+(y^2)]),λ,f,Di]],{z,f-15,f+5},{x,-10,10}, PlotRange-> All,  PlotPoints-> {80,40}, BoxRatios-> {2,2,0.5}]

Alternate code if slow computer:

λ=1.  ; k=2 π/λ  ;
f=5.*10^3 ; Di=5*10^3 ;   y=0. ; field={} ;
For[i=0,i≤ 80,i++, z=f-15+0.25*i ;
dat=Table[{z,j*0.5,Re[UniformBeam[z,(Sqrt[(((j*0.5))^2)+(y^2)]),λ,f,Di]]} ,{j,-20,20}]; 
AppendTo[field,dat]
];
pfield=Flatten[field,1 ];
ListPlot3D[pfield, PlotRange-> All, BoxRatios->{1,1,0.3}]

ListDensityPlot[pfield,ColorFunction->GrayLevel, PlotRange->All]

Exercise 2   Focal plane intensity image of a uniform beam

Let (Subscript[U, in])[ξ,η] be uniform and P[ξ,η] be a circle, what is the image pattern? From the above, we expect it to be similar to Fraunhofer diffraction of a circle.  
        (1/(2 π))∫ P[ξ,η] (e^(-i  2π (p  ξ + q  η)))dξ dη ->(1/(2 π))∫ (e^(-i  2π s ρ Cos[ϕ]))ρdρ dϕ  
where: s={p,q}=(1/(λ f)) {x,y}  and ρ={ξ,η} and  Cos[ϕ]= r . ρ

(1/(2 π)) (∫_0^(2 π))(E^(-I  2π s ρ Cos[ϕ]))\[DifferentialD]ϕ

(∫_0^(D/2))(Subscript[J, 0])(2 π s ρ)ρ\[DifferentialD]ρ  /. s-> (1/(λ f)) r

Or:    k ((e^(i k f))/(2 i f)) (e^(i (k/(2 f))((x^2)+(y^2)))) (((D^2)  (Subscript[J, 1])( r/w))/(4 (r D π/f λ))) =k((D^2)/8)((e^(i k f))/(i f)) (e^(i (k/(2 f))((x^2)+(y^2)))) ((  (Subscript[J, 1])( r/w))/((r/w)))   
                                       =-i  ((π (D^2))/(4 λ f))(e^(i k f)) (e^(i (k/(2 f))(r^2)))((  (Subscript[J, 1])( r/w))/((r/w)))   
where  w=((f λ)/(π D))

λ=1.  ;
f=5.*10^3 ; Di=5*10^3 ;
Manipulate[w=f λ /(π Di) ;
ParametricPlot3D[{r Cos[t], r Sin[t],(π (Di^2)/(4 λ f))^2  (w  (Subscript[J, 1])(r/w)/r)^2},{r,0.001,5}, {t,0,2 π}, PlotRange-> All, BoxRatios-> {1,1,0.5}],{Di,2*10^3, 10^4}]

The beam waist diameter w=((3.2326798966214065` f λ)/(π D)) ~((f λ)/D)  shows that the smaller the lens f-number: f/D, the smaller the waist, i. e. tigher focusing, which indicates the lens optical resolution.  
For the first zero, it is w=2((1.22 f λ)/D) ~2.44((f λ)/D)

1.3 Object at the focal plane

Suppose we have the opposite case from 1.2 above; an object is at the lens focal plane instead of being right in front of the lens, what will be the pattern after the lens? (note the change of the coordinate labeling: {x,y} is for input and {ξ, η} is for output

We can apply Fresnel theory as above to find the diffracted field from the object to the plane before the lens. It is:
           (Subscript[U, 1])[ξ,η, (Subscript[z, 1])] =k ((e^(i k (Subscript[z, 1])))/(2 π i (Subscript[z, 1]))) ∫ T[x,y] (e^(i (k/(2 (Subscript[z, 1])))( (((x-ξ))^2)+(((y-η))^2)) ))dx dy           (1.3.1)
Here, we assume the illumination is uniform. Immediately after the lens:
 (Subscript[U, 2])[ξ,η,(Subscript[z, 1])] =k ((e^(i k (Subscript[z, 1])))/(2 π i (Subscript[z, 1]))) (e^(-i (k/(2 f))((ξ^2)+(η^2)) ))P[ξ,η]∫T[x,y] (e^(i (k/(2 (Subscript[z, 1])))( (((x-ξ))^2)+(((y-η))^2)) ))dx dy      (1.3.2)
where P[ξ,η] is the lens aperture function. With similar algebraic manipulation:
   (Subscript[U, 2])[ξ,η,(Subscript[z, 1])] =k ((e^(i k (Subscript[z, 1])))/(2 π i (Subscript[z, 1]))) (e^(i (k/2) ((1/(Subscript[z, 1]))-(1/f))((ξ^2)+(η^2))))P[ξ,η]∫T[x,y] (e^(i (k/(2 (Subscript[z, 1])))((x^2)+(y^2)))) (e^(-i (k/(Subscript[z, 1]))(x ξ+y η) ))dx dy  (1.3.3)
We see here that we do not simply get the Fourier transform of T[x, y], but that of the function:
                          (OverTilde[T])[x,y;(Subscript[z, 1])]=T[x,y] (e^(i (k/(2 (Subscript[z, 1])))((x^2)+(y^2)) ))                                          (1.3.4)
In other words, the transmission function is modified by a phase term (e^(i (k/(2 (Subscript[z, 1])))( (x^2)+(y^2)) )) (think of it as ((π (w^2))/(λ z))). Later on, we will see that this phase term is not crucial as it can be neglected or corrected in flat field system. The term is relevant only when the image is very large, i. e. π ((x^2)+(y^2))>> λ (Subscript[z, 1]).

Then, from (1.3.3): 
       (Subscript[U, 2])[ξ,η,(Subscript[z, 1])] =k ((e^(i k (Subscript[z, 1])))/(i (Subscript[z, 1]))) (e^(i (k/2) ((1/(Subscript[z, 1]))-(1/f))((ξ^2)+(η^2))))P[ξ,η](F^(-1))[(OverTilde[T])[(Subscript[z, 1])]][k(ξ/(Subscript[z, 1])),k(η/(Subscript[z, 1]))]             (1.3.5)
where (F^(-1))[(OverTilde[T])[(Subscript[z, 1])]] is the inverse Fourier transform of (OverTilde[T])[(Subscript[z, 1])]:
           (F^(-1))[(OverTilde[T])[(Subscript[z, 1])]] [(Subscript[q, x]),(Subscript[q, y])]=(1/(2 π))∫(OverTilde[T])[x,y;(Subscript[z, 1])](e^(-i (x (Subscript[q, x])+y (Subscript[q, y])) ))dx dy                     (1.3.6)

If the object is at the front focal plane: 
                                              ((1/f)-(1/(Subscript[z, 1])))=0                                                        (1.3.7a)
then:           (Subscript[U, 2])[ξ,η,(Subscript[z, 1])=f] =k ((e^(i k f))/(i f)) P[ξ,η](F^(-1))[(OverTilde[T])[f]][k(ξ/f),k(η/f)]            (1.3.7b)

which is the equivalent reciprocal of the results above in section 1.2, Eq. (1.2.8b).
Here it is for comparison. From section 1.2, Eq. (1.2.8):
  U[x,y](Subscript[|, z=f]) =k ((e^(i k z))/(i z))  (e^(i (k/(2 z))((x^2)+(y^2))))(F^(-1))[T[ξ,η] P[ξ,η]][k(x/z),k(y/z)]      (1.2.8a)
                         =k ((e^(i k f))/(i f)) (e^(i (k/(2 f))((x^2)+(y^2)))) (F^(-1))[T P][k(x/f),k(y/f)]                   (1.2.8b)    
Both involve the inverse FT of image pattern T[x,y], differing only on the role of the lens aperture function P[ξ,η].
As mentioned above,  (k/(2 (Subscript[z, 1])))((x^2)+(y^2))~((π (w^2))/(λ (Subscript[z, 1])))  is supposed to be slowly varying only in this approximation. In fact, for ideal paraxial approximation, ((π (w^2))/(λ (Subscript[z, 1])))<<1, hence, the phase term (e^(i (k/(2 (Subscript[z, 1])))((x^2)+(y^2)) ))does not play a significant role.

1.4 Focal plane input - Focal plane output

Finally, what happens if we have the case below:

We simply use the result from 1.3 to obtain (Subscript[U, 2]) and then, Fresnel diffraction to propagate the wave from (Subscript[U, 2]) to (Subscript[U, 3]). From Eq. (1.3.3): 
     (Subscript[U, 2])[ξ,η,(Subscript[z, 1])] =k ((e^(i k (Subscript[z, 1])))/(2 π i (Subscript[z, 1]))) (e^(i (k/2) ((1/(Subscript[z, 1]))-(1/f))((ξ^2)+(η^2))))P[ξ,η]∫T[x,y] (e^(i (k/(2 (Subscript[z, 1])))((x^2)+(y^2)))) (e^(-i (k/(Subscript[z, 1]))(x ξ+y η) ))dx dy     (1.4.1)
With  (1/(Subscript[z, 1]))-(1/f)=0:
      (Subscript[U, 2])[ξ,η,(Subscript[z, 1])] =k ((e^(i k (Subscript[z, 1])))/(2 π i (Subscript[z, 1]))) P[ξ,η]∫T[x,y] (e^(i (k/(2 (Subscript[z, 1])))((x^2)+(y^2)))) (e^(-i (k/(Subscript[z, 1]))(x ξ+y η) ))dx dy             (1.4.2)
We apply Fresnel diffraction:
         (Subscript[U, 3])[u,v,(Subscript[z, 2])] = k ((e^(i k (Subscript[z, 2])))/(i 2 π (Subscript[z, 2]))) (e^(i (k/(2 (Subscript[z, 2])))((u^2)+(v^2)) ))∫(Subscript[U, 2])[ξ,η] (e^(i (k/(2 (Subscript[z, 2])))(ρ^2) ))(e^(-i (k/(Subscript[z, 2]))(ξ u+η v) ))dξ dη      (1.4.3)  
Substitute (1.4.2) into (1.4.3):
      (Subscript[U, 3])[u,v,(Subscript[z, 2])] = k ((e^(i k (Subscript[z, 2])))/(i 2 π (Subscript[z, 2]))) (e^(i (k/(2 (Subscript[z, 2])))((u^2)+(v^2)) ))k ((e^(i k (Subscript[z, 1])))/(2 π i (Subscript[z, 1])))∫P[ξ,η] (e^(i (k/(2 (Subscript[z, 2])))(ρ^2) ))(e^(-i (k/(Subscript[z, 2]))(ξ u+η v) ))
                                                        (∫T[x,y] (e^(i (k/(2 (Subscript[z, 1])))((x^2)+(y^2)))) (e^(-i (k/(Subscript[z, 1]))(x ξ+y η) ))dx dy)dξ dη      (1.4.4)
We switch order of integration:
      (Subscript[U, 3])[u,v,(Subscript[z, 2])] = k ((e^(i k (Subscript[z, 2])))/(i 2 π (Subscript[z, 2]))) (e^(i (k/(2 (Subscript[z, 2])))((u^2)+(v^2)) ))k ((e^(i k (Subscript[z, 1])))/(2 π i (Subscript[z, 1])))∫T[x,y] (e^(i (k/(2 (Subscript[z, 1])))((x^2)+(y^2)))) 
                                                        (∫P[ξ,η] (e^(i (k/(2 (Subscript[z, 2])))(ρ^2) ))(e^(-i (k/(Subscript[z, 2]))(ξ u+η v) ))(e^(-i (k/(Subscript[z, 1]))(x ξ+y η) ))dξ dη)dx dy      (1.4.5)
We look at this term:
                    ∫P[ξ,η] (e^(i (k/(2 (Subscript[z, 2])))(ρ^2) ))(e^(-i (k/(Subscript[z, 2]))(ξ u+η v) ))(e^(-i (k/(Subscript[z, 1]))(x ξ+y η) ))dξ dη                    (1.4.6)
which is:
                             ∫P[ξ,η] (e^(i (k/(2 (Subscript[z, 2])))(ρ^2) ))(e^(-i k(ξ ( (u/(Subscript[z, 2]))+(x/(Subscript[z, 1])))+η ( (v/(Subscript[z, 2]))+(y/(Subscript[z, 1])))) ))dξ dη                (1.4.7)
It cannot be integrated in close form if include P[ξ, η]. But if we can approximate that the lens is very large compared with the object size, then P[ξ, η] can be consider as 1 over the entire ξ, η plane. Then, (1.4.7) becomes:
                            ~∫(e^(i (k/(2 (Subscript[z, 2])))(ρ^2) ))(e^(-i k(ξ ( (u/(Subscript[z, 2]))+(x/(Subscript[z, 1])))+η ( (v/(Subscript[z, 2]))+(y/(Subscript[z, 1])))) ))dξ dη                (1.4.8a)
                            ∫(e^(i (k/(2 (Subscript[z, 2])))(ρ^2) ))(e^(-i (ξ s+η t) ))dξ dη                (1.4.8b)
where:                          s≡k ((u/(Subscript[z, 2]))+(x/(Subscript[z, 1])))  ; t≡k ((v/(Subscript[z, 2]))+(y/(Subscript[z, 1])))                 (1.4.8c)
This is the Fourier transform of (e^(i (k/(2 (Subscript[z, 2])))(ρ^2) )). Fortunately, we can perform this FT:
                            ∫(e^(i (k/(2 (Subscript[z, 2])))(ρ^2) ))(e^(-i (ξ s+η t) ))dξ dη= ((i 2 π (Subscript[z, 2]))/k) (e^(-i ((Subscript[z, 2])/(2 k)) ((s^2)+(t^2))))     (1.4.9)
 Apply (1.4.9) into 1.4.5:
      (Subscript[U, 3])[u,v,(Subscript[z, 2])] = k ((e^(i k (Subscript[z, 2])))/(i 2 π (Subscript[z, 2]))) (e^(i (k/(2 (Subscript[z, 2])))((u^2)+(v^2)) ))k ((e^(i k (Subscript[z, 1])))/(2 π i (Subscript[z, 1])))
                                             ∫T[x,y] (e^(i (k/(2 (Subscript[z, 1])))((x^2)+(y^2))))     ((i 2 π (Subscript[z, 2]))/k) (e^(-i ((Subscript[z, 2])/(2 k)) ((s^2)+(t^2)))) dx dy   (1.4.10a)
   Substuitute (1.4.8c): 
       (Subscript[U, 3])[u,v,(Subscript[z, 2])] = k ((e^(i k ((Subscript[z, 1])+(Subscript[z, 2]))))/(2 π i (Subscript[z, 1]))) (e^(i (k/(2 (Subscript[z, 2])))((u^2)+(v^2)) ))
                            ∫T[x,y] (e^(i (k/(2 (Subscript[z, 1])))((x^2)+(y^2)))) (e^(-i ((Subscript[z, 2])/(2 k)) (k^2) (((((u/(Subscript[z, 2]))+(x/(Subscript[z, 1]))))^2)+((((v/(Subscript[z, 2]))+(y/(Subscript[z, 1]))))^2)))) dx dy      (1.4.10b)
 Expanding the exponent with u, v, x, and y:
       (Subscript[U, 3])[u,v,(Subscript[z, 2])] =k ((e^(i k ((Subscript[z, 1])+(Subscript[z, 2]))))/(2 π i (Subscript[z, 1])))  (e^(i (k/(2 (Subscript[z, 2])))((u^2)+(v^2)) )) (e^(-i (k/(2 (Subscript[z, 2])))((u^2)+(v^2)) ))
                                  ∫T[x,y] (e^(i (k/(2 (Subscript[z, 1])))((x^2)+(y^2)))) (e^(-i ((k (Subscript[z, 2]))/(2 ((Subscript[z, 1])^2)))((x^2)+(y^2))))(e^(-i k  (u(x/(Subscript[z, 1]))+v(y/(Subscript[z, 1]))))) dx dy      (1.4.10c)
We see that all quadratic exponent terms are cancelled if (Subscript[z, 1])=(Subscript[z, 2])=f 
                 (Subscript[U, 3])[u,v,(Subscript[z, 2])] =  k ((e^(2i k f))/(2 π i f))∫T[x,y]  (e^(-i k  (u(x/f)+v(y/f)))) dx dy      (1.4.11a)
                             (Subscript[U, 3])[u,v,(Subscript[z, 2])] =  k ((e^(2i k f))/(2 π i f))(F^(-1))[T][k(u/f),k(v/f)]            (1.4.11b)

This is a most well known results in lens application and Fourier optics: the focal plane-to-focal plane field transformation is a Fourier transform (without any quadratic phase term), if the lens can be approximated as large compared with the object size. This has been used extensively for analog optical processing.

Think of the lens as the “Fresnel” quadratic phase remover.

1.5. Model for image formation with two-lens equivalent.

1.5.1 Image formation

Suppose we have 2 lenses back-to-back as shown:

What we can do here is to apply Fresnel diffraction theory twice for each stage:
1- From the object to the 1st lens
2- From the 2nd lens to its back focal plane
For the first stage, we see that from Eq. (1.3.7b)
         (Subscript[U, 2])[ξ,η,(Subscript[z, 1])=f] =k ((e^(i k (Subscript[z, 1])))/(i (Subscript[z, 1]))) (Subscript[P, 1])[ξ,η](F^(-1))[(OverTilde[T])[(Subscript[z, 1])]][k(ξ/(Subscript[z, 1])),k(η/(Subscript[z, 1]))]                (1.3.7b)-(1.5.1)   
where we use the function (Subscript[P, 1])[ξ,η] to denote aperture of the "first" lens. But now, this field is the input to the second lens. Therefore we can use (1.2.8b) to get the field in the second lens back focal plane, substituting (Subscript[z, 2]) where appropriate:
                (Subscript[U, 3])[u,v]=k ((e^(i k (Subscript[z, 2])))/(i (Subscript[z, 2]))) (e^(i (k/(2 (Subscript[z, 2])))((u^2)+(v^2)))) (F^(-1))[(Subscript[U, 2]) (Subscript[P, 2])][k(u/(Subscript[z, 2])),k(v/(Subscript[z, 2]))]                             (1.5.2)
where (Subscript[P, 2]) denotes aperture of the "second" lens. Of course, both (Subscript[P, 1]) and (Subscript[P, 2]) are the same, and later we will just use one function for their product.
Here, substitute (1.3.7b-1.5.1) into (1.5.2)
  (Subscript[U, 3])[u,v] =k ((e^(i k (Subscript[z, 2])))/(i (Subscript[z, 2]))) (e^(i (k/(2 (Subscript[z, 2])))((u^2)+(v^2)))) (F^(-1))[k ((e^(i k (Subscript[z, 1])))/(i (Subscript[z, 1]))) (Subscript[P, 1])[ξ,η](F^(-1))[(OverTilde[(Subscript[T, (Subscript[z, 1])])])][k(ξ/(Subscript[z, 1])),k(η/(Subscript[z, 1]))](Subscript[P, 2])[ξ,η]][k(u/(Subscript[z, 2])),k(v/(Subscript[z, 2]))]   (1.5.3a) 
   (Subscript[U, 3])[u,v] = -((k^2)/((Subscript[z, 2])(Subscript[z, 1]))) (e^(i k ((Subscript[z, 2])+(Subscript[z, 1])) ))(e^(i (k/(2 (Subscript[z, 2])))((u^2)+(v^2)))) (F^(-1))[P[ξ,η](F^(-1))[(OverTilde[(Subscript[T, (Subscript[z, 1])])])][k(ξ/(Subscript[z, 1])),k(η/(Subscript[z, 1]))]  ][k(u/(Subscript[z, 2])),k(v/(Subscript[z, 2]))]      (1.5.3b)
Hence, the double Fourier transform shows that the image is reproduced but with a cut-off in frequency space because of the lens finite aperture represented by P[ξ, η].

Explicitly, the double Fourier term:
                          (F^(-1))[P[ξ,η](F^(-1))[(OverTilde[(Subscript[T, (Subscript[z, 1])])])][k(ξ/(Subscript[z, 1])),k(η/(Subscript[z, 1]))]  ][k(u/(Subscript[z, 2])),k(v/(Subscript[z, 2]))]                       (1.5.4a)
 is:
 (1/(4 (π^2)))∫dξ dη  (e^(- i k (ξ (u/(Subscript[z, 2]))+η (v/(Subscript[z, 2])) )))P[ξ,η]∫dx dy (e^(-i k(x (ξ/(Subscript[z, 1]))+ y(η/(Subscript[z, 1])))))  T[x,y] (e^(i (k/(2 (Subscript[z, 1])))((x^2)+(y^2)) ))

Swapping the order of integration to bring together the middle terms with {ξ,η}: (so that we can link directly from input {x,y} to output {u,v})

 (1/(4 (π^2)))∫dx dy  T[x,y] (e^(i (k/(2 (Subscript[z, 1])))((x^2)+(y^2)) ))∫dξ dη P[ξ,η](e^(- i k (ξ ((x/(Subscript[z, 1]))+(u/(Subscript[z, 2])))+η ((y/(Subscript[z, 1]))+(v/(Subscript[z, 2]))) )))  (1.5.4b)

We will call this function below "optical transfer function" of the lens:
                        (1/(4 (π^2)))∫dξ dη P[ξ,η](e^(- i k (ξ ((x/(Subscript[z, 1]))+(u/(Subscript[z, 2])))+η ((y/(Subscript[z, 1]))+(v/(Subscript[z, 2]))) )))                      (1.5.5)

If P[ξ, η] is 1 everywhere (infinity), then:
       (1/(4 (π^2)))∫dξ dη P[ξ,η](e^(- i k (ξ ((x/(Subscript[z, 1]))+(u/(Subscript[z, 2])))+η ((y/(Subscript[z, 1]))+(v/(Subscript[z, 2]))) ))) = δ[k ((x/(Subscript[z, 1]))+(u/(Subscript[z, 2]))),k ((y/(Subscript[z, 1]))+(v/(Subscript[z, 2])))]   (1.5.6a)
and:
         (F^(-1))[P[ξ,η](F^(-1))[(OverTilde[(Subscript[T, (Subscript[z, 1])])])][k(ξ/(Subscript[z, 1])),k(η/(Subscript[z, 1]))]  ][k(u/(Subscript[z, 2])),k(v/(Subscript[z, 2]))]
             =∫dx dy  T[x,y] (e^(i (k/(2 (Subscript[z, 1])))((x^2)+(y^2)) ))   δ[k ((x/(Subscript[z, 1]))+(u/(Subscript[z, 2]))),k ((y/(Subscript[z, 1]))+(v/(Subscript[z, 2])))] 
                  =(((Subscript[z, 1])^2)/(k^2))T[-((Subscript[z, 1])/(Subscript[z, 2]))u,-((Subscript[z, 1])/(Subscript[z, 2]))v](e^(i ((k (Subscript[z, 1]))/(2 ((Subscript[z, 2])^2)))((u^2)+(v^2))))                                     (1.5.6b)
And from (1.5.3b), the image is ~ :
          (Subscript[U, 3])[u,v] = -((k^2)/((Subscript[z, 2])(Subscript[z, 1]))) (e^(i k ((Subscript[z, 2])+(Subscript[z, 1])) ))(e^(i (k/(2 (Subscript[z, 2])))((u^2)+(v^2))))(((Subscript[z, 1])^2)/(k^2))T[-((Subscript[z, 1])/(Subscript[z, 2]))u,-((Subscript[z, 1])/(Subscript[z, 2]))v](e^(i ((k (Subscript[z, 1]))/(2 ((Subscript[z, 2])^2)))((u^2)+(v^2))))
          (Subscript[U, 3])[u,v] =-((Subscript[z, 1])/(Subscript[z, 2])) (e^(i k ((Subscript[z, 2])+(Subscript[z, 1])) ))T[-((Subscript[z, 1])/(Subscript[z, 2]))u,-((Subscript[z, 1])/(Subscript[z, 2]))v] (e^(i (k/(2 (Subscript[z, 2])))(1+((Subscript[z, 1])/(Subscript[z, 2])))( (u^2)+(v^2))))           (1.5.7)
which is essentially the object original transmission function with magnification factor -((Subscript[z, 1])/(Subscript[z, 2])); the minus sign indicates the inversion, the ratio ((Subscript[z, 1])/(Subscript[z, 2])) indicates the magnification factor. Infinite ideal lens produces perfectly faithful images.

Thus: The lens is an anti-Fresnel diffraction device: the illuminated image (pattern) is blurred at a distance because of Fresnel diffraction. A lens (perfect) undoes that effect and reconstruct the original illumination image (pattern).

For finite aperture, we define:
       (1/(4 (π^2)))∫dξ dη P[ξ,η](e^(- i k (ξ ((x/(Subscript[z, 1]))+(u/(Subscript[z, 2])))+η ((y/(Subscript[z, 1]))+(v/(Subscript[z, 2]))) ))) = S[k ((x/(Subscript[z, 1]))+(u/(Subscript[z, 2]))),k ((y/(Subscript[z, 1]))+(v/(Subscript[z, 2])))]    (1.5.8)
which is the pointspread function (PSF) of the system. Then the image is:
   (Subscript[U, 3])[u,v] =-((k^2)/((Subscript[z, 2])(Subscript[z, 1]))) (e^(i k ((Subscript[z, 2])+(Subscript[z, 1])) ))(e^(i (k/(2 (Subscript[z, 2])))((u^2)+(v^2)))) ∫dx dy T[x,y] (e^(i (k/(2 (Subscript[z, 1])))((x^2)+(y^2))))  S[k ((x/(Subscript[z, 1]))+(u/(Subscript[z, 2]))),k ((y/(Subscript[z, 1]))+(v/(Subscript[z, 2])))]   (1.5.9)  
which is obviously a convolution. In the FT space (spatial frequency domain), this convolution is a product, hence we will see that the lens acts just like a hard (sharp cut-off) filter in the spatial frequency domain.

Calculate S for a circular aperture

(∫_0^(2 π))(E^(I ρ q Cos[ϕ] ))\[DifferentialD]ϕ

(1/(2 π)) (∫_0^a)(BesselJ[0, q ρ])ρ\[DifferentialD]ρ

Hence:      (1/(4 (π^2)))∫dξ dη P[ξ,η](e^(- i (ξ p+η q )))=((a^2)/(2 π))(( (BesselJ[1, a q]))/((a q)))  where a is the radius of the lens aperture. This is again the point spread function that we had earlier.

The image (1.5.9) is a convolution:
      -((k^2)/((Subscript[z, 2])(Subscript[z, 1]))) (e^(i k ((Subscript[z, 2])+(Subscript[z, 1])) ))(e^(i (k/(2 (Subscript[z, 2])))((u^2)+(v^2))))∫dx dy T[x,y] (e^(i (k/(2 (Subscript[z, 1])))( (x^2)+(y^2))))  ((a^2)/(2 π))(( (BesselJ[1, a k (Sqrt[ ((((x/(Subscript[z, 1]))+(u/(Subscript[z, 2]))))^2)+ ((((y/(Subscript[z, 1]))+(v/(Subscript[z, 2]))))^2)])]))/(a k (Sqrt[ ((((x/(Subscript[z, 1]))+(u/(Subscript[z, 2]))))^2)+ ((((y/(Subscript[z, 1]))+(v/(Subscript[z, 2]))))^2)])))     (1.5.11a)
Or:
      -((k^2)/((Subscript[z, 2])(Subscript[z, 1]))) (e^(i k ((Subscript[z, 2])+(Subscript[z, 1])) ))(e^(i (k/(2 (Subscript[z, 2])))((u^2)+(v^2))))∫dx dy T[x,y] (e^(i (k/(2 (Subscript[z, 1])))( (x^2)+(y^2))))  ((a^2)/(2 π))(( (BesselJ[1, (a/(Subscript[z, 1])) k (Sqrt[ (((-((Subscript[z, 1])/(Subscript[z, 2]))u-x))^2)+ (((-((Subscript[z, 1])/(Subscript[z, 2]))v-y))^2)])]))/((a/(Subscript[z, 1])) k (Sqrt[(((-((Subscript[z, 1])/(Subscript[z, 2]))u-x))^2)+ (((-((Subscript[z, 1])/(Subscript[z, 2]))v-y))^2)])))   (1.5.11b)

1.5.2 Delta-function transmission function

What if the transmission function of the object is a point, i. e. δ-function at some location {(Subscript[x, 0]),(Subscript[y, 0])}:
                                                 T[x,y]=δ[x-(Subscript[x, 0]),y-(Subscript[y, 0])]                       (1.5.12)
Then:
    (Subscript[U, 3])[u,v] =-(((k^2)(a^2))/(2 π (Subscript[z, 2])(Subscript[z, 1]))) (e^(i k ((Subscript[z, 2])+(Subscript[z, 1])) ))(e^(i (k/(2 (Subscript[z, 2])))((u^2)+(v^2)))) (e^(i (k/(2 (Subscript[z, 1])))( ((Subscript[x, 0])^2)+((Subscript[y, 0])^2))))  (( (BesselJ[1, (a/(Subscript[z, 1])) k (Sqrt[ (((-((Subscript[z, 1])/(Subscript[z, 2]))u-(Subscript[x, 0])))^2)+ (((-((Subscript[z, 1])/(Subscript[z, 2]))v-(Subscript[y, 0])))^2)])]))/((a/(Subscript[z, 1])) k (Sqrt[ (((-((Subscript[z, 1])/(Subscript[z, 2]))u-(Subscript[x, 0])))^2)+ (((-((Subscript[z, 1])/(Subscript[z, 2]))v-(Subscript[y, 0])))^2)]))) 
                 =-(((k^2)(a^2))/(2 π (Subscript[z, 2])(Subscript[z, 1]))) (e^(i k ((Subscript[z, 2])+(Subscript[z, 1])) ))(e^(i (k/(2 (Subscript[z, 2])))((u^2)+(v^2)))) (e^(i (k/(2 (Subscript[z, 1])))( ((Subscript[x, 0])^2)+((Subscript[y, 0])^2))))  (( (BesselJ[1, (a/(Subscript[z, 2])) k (Sqrt[ (((u-(Subscript[u, 0])))^2)+ (((v-(Subscript[v, 0])))^2)])]))/((a/(Subscript[z, 2])) k (Sqrt[  (((u-(Subscript[u, 0])))^2)+ (((v-(Subscript[v, 0])))^2)])))                  (1.5.13)
which, we see is just the point-spread function (PSF) (as it is aptly named) at the corresponding location   {-((Subscript[z, 2])/(Subscript[z, 1]))(Subscript[x, 0]), -((Subscript[z, 2])/(Subscript[z, 1]))(Subscript[y, 0])} in the image plane. We will get the same result if the light comes from an emitting point source ((e^(i k r))/r). (see HW).

1.6 Exercise: Uniform beam

1.6.1 Input plane wave

Exercise  Uniform beam (input plane wave)

Let (Subscript[U, in])[ξ,η] be uniform and P[ξ,η] be a circle, what is the image pattern at the focal plane? From the above, we expect it to be similar to Fraunhofer diffraction of a circle.  
        (1/(2 π))∫ P[ξ,η] (e^(-i  2π (p  ξ + q  η)))dξ dη ->(1/(2 π))∫ (e^(-i  2π s ρ Cos[ϕ]))ρdρ dϕ  
where: s={p,q}=(1/(λ f)) {x,y}  and ρ={ξ,η} and  Cos[ϕ]= r . ρ

(1/(2 π)) (∫_0^(2 π))(E^(-I  2π s ρ Cos[ϕ]))\[DifferentialD]ϕ

(∫_0^(D/2))(Subscript[J, 0])(2 π s ρ)ρ\[DifferentialD]ρ  /. s-> (1/(λ f)) r

Or:    k ((e^(i k f))/(2 i f)) (e^(i (k/(2 f))((x^2)+(y^2)))) (((D^2)  (Subscript[J, 1])( r/w))/(4 (r D π/f λ))) =k((D^2)/8)((e^(i k f))/(i f)) (e^(i (k/(2 f))((x^2)+(y^2)))) ((  (Subscript[J, 1])( r/w))/((r/w)))   
                                       (Subscript[U, focal])=-i  ((π (D^2))/(4 λ f))(e^(i k f)) (e^(i (k/(2 f))(r^2)))((  (Subscript[J, 1])( r/w))/((r/w)))  
((  (Subscript[J, 1])( r/w))/((r/w))) is known as the lens point-spread function (see next section).
where  w=((f λ)/(π D))
The ratio  (f/D) is define as a lens f-number. For example, f/2 means that the focal length is 2 times the lens aperture diameter. It is inversely proportional to lens numerical aperture: Sin[θ].
w is a measure of lens resolution (more later), hence, the lower the f-number the higher the resolution.

Manipulate[f=(10^4) ; λ=0.55 ; Di=diam*10^4 ;
w=f λ /(π Di) ;
{ParametricPlot3D[{r Cos[t], r Sin[t],(π (Di^2)/(4 λ f))^2  (w BesselJ[1,r/w]/r )^2 },{r,0.001,5}, {t,0,2 π}, PlotRange-> All, BoxRatios-> {1,1,0.5}, Mesh-> False]
,DensityPlot[(π (Di^2)/(4 λ f))^2  (w BesselJ[1,(Sqrt[(x^2)+(y^2)])/w]/(Sqrt[(x^2)+(y^2)]) )^2 ,{x,-5,5}, {y,-5,5}, PlotRange-> All, Mesh-> False, ColorFunction->GrayLevel]
,ParametricPlot3D[{r Cos[t], r Sin[t],Arg[(E^(I ((π (r^2))/(λ f)))) ] },{r,0.001,5}, {t,0,2 π}, PlotRange-> All, BoxRatios-> {1,1,0.5}, Mesh-> False]},{{diam,0.2,"Lens diameter (cm)"},0.2, 2}]

We define the spot diameter is where the intensity has its first zero:

r= . ;
FindRoot[(Subscript[J, 1])(r)/r== 0.,{r,3}]

3.8317059702075125`/Pi

The beam waist w=((3.8317059702075125` f λ)/(π D)) =1.22((f λ)/D)  shows that the smaller the lens f-number: f/D, the smaller the waist, i. e. tigher focusing, which indicates the lens optical resolution.
The diameter of the spot 2 w=2((1.22 f λ)/D) ~2.44((f λ)/D) 
Compare this formula with Gaussian beam waist:  w~(λ/(π  Sin[θ]))~ (λ/(π Tan[θ])), 
where θ is the beam divergence:  Tan[θ] ~(D/(2 f)) . Hence w~ ((2 f λ)/(π D))=0.64 ((f λ)/D) , 
Note that Gaussian w is at the point (e^(-1)) and not at the full circle zero of the point spread function, hence it is ~half width rather than full width.

1.6.2 HW: Input plane wave with a tilt angle on a square lens

Source code: Square lens focusing of plane wave in paraxial approximation.

Manipulate[
If[RUN== 1,

Manipulate[k=2 π*(10^(3 )) ; αr=α*π/180. ;
lens={Hue[0.45,1,1,1],Disk[{0,0},{0.2,a*100}]} ;
axis=Line[{{-2,0},{5,0}}];
beamaxis={Red,Line[{{-2,-2 Tan[αr]},{5,5 Tan[αr]}}]};
beamhue=Hue[0.07,1,1,0.3] ;
beamleft={beamhue,Polygon[{{0,a*100},{0,-a*100}
,{-2,-a*100-2Tan[αr]},{-2,a*100-2Tan[αr]},{0,a*100}}]};
p1= 100*{f,f*Tan[αr]+((f π)/(a Cos[αr] k))} ;
p2= 100*{f,f*Tan[αr]-((f π)/(a Cos[αr] k))};
beamright={beamhue,Polygon[{{0,a*100},p1
,p2,{0,-a*100},{0,a*100}}]};
GraphicsRow[{Show[Graphics[{
lens,axis, beamaxis,beamleft,beamright }]
,PlotRange-> {{-2,5},{-4,4}}, ImageSize->{7*500/8.,500}
]
,DensityPlot[(((((Sin[((a k x)/f)Cos[αr]]Sin[a k ((( Cos[αr])/f)(y-f Tan[αr]))])/((a^2) (k^2)  ((Cos[αr])^2)((x/f))((y/f) -Tan[αr])))))^2),{x,-0.01,0.01},{y,-0.04,0.04}, ColorFunction->GrayLevel
, PlotRange->{{-0.01,0.01},{-0.04,0.04},{0,1}}, PlotPoints->{40,200}
,Frame-> None, FrameLabel-> False
,AspectRatio->4, ImageSize->{125,500}] 
}]

,{{f,0.02},0.01,0.05},{{a,0.01},0.005,0.02},{{α,0},-60,60, ControlPlacement->Right, ControlType->VerticalSlider}]
,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]]
,Button["LENS FOCUSING VS INPUT ANGLE",BaseStyle->{24,Black, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.4,1,1]]
,Row[{Button["Toggle box to start, and turn off when not running",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUN,0,"->"},{0,1}}]
}]
]

Demo: Square lens focusing of plane wave in paraxial approximation.

1.6.3 HW - Point source

Derive the expression for point source image. Show Airy function