Defocus transfer function for a circular aperture

From Optical Imaging and Spectroscopy (Wiley - OSA 2009)
By David J. Brady
www.opticalimaging.org
Figures 6.20 and 6.21

H[u_,tz_]:=NIntegrate[4*Cos[2 Pi tz u bx]UnitStep[0.5- Sqrt[(bx-u/2)^2 + by^2]] UnitStep[0.5- Sqrt[(bx+u/2)^2 + by^2]]/Pi, {bx,-10,10},{by,-10,10}];

Plot[H[u,0],{u,-3,3},PlotRange->All,Axes->False,Frame->True]

Plot3D[H[u,2 w20],{u,-1.1,1.1},{w20,0,10},PlotRange->All,AxesLabel->{u,(Subscript[w, 20]),H},PlotRange->All,PlotPoints->75,Lighting->{{"Directional",GrayLevel[.8,.5],{{5,5,4},{5,5,0}}}}]

<<Version5`Graphics`

Plot3D[H[u,2 w20],{u,-1.1,1.1},{w20,0,10},PlotRange->All,AxesLabel->{u,(Subscript[w, 20]),H},PlotRange->All,PlotPoints->50,ColorOutput->GrayLevel,MeshStyle->Gray]

Ho[u_]=Integrate[4UnitStep[0.5- Sqrt[(bx-u/2)^2 + by^2]] UnitStep[0.5- Sqrt[(bx+u/2)^2 + by^2]]/Pi, {bx,-10,10},{by,-10,10}]

Plot[Ho[u],{u,-3,3}]

pip=Table[{0,0},{i,30}];
For[pop=1,pop<31,++pop,u=(pop/31);p=Interpolation[Table[{w,((H[u,2 w])/(Ho[u]))-0.9},{w,0,10,0.5}]];pip[[pop]]={u,x/.FindRoot[p[x]==0,{x,1}]};];
Clear[u] ;
ListPlot[pip,Joined->True,Axes->False,Frame->True,FrameLabel->{u,(Subscript[w, 20])},PlotRange->All]

pip

pip05=pip

pip08=pip

pip09=pip

pap=pip

ListPlot[{pip05,pip08,pip09},Joined->True,Axes->False,Frame->True,FrameLabel->{u,(Subscript[w, 20])},PlotRange->{{0,1},{0,10}},PlotStyle->{Black}]