stemhr.m

扫描电镜（stem）的matlab模拟程序代码

function psf = stemhr(r,params,nk)
%  STEMHR STEM point spread function (PSF)
%  MATLAB function to calculate STEM probe profile vs. r
%    input array r has the radial positions (in Angs.)
%    input array params has the optical parameters
%          params ={ Cs, df, kev, amax [,d0]}
%          the optional input d0 is the source size (in Ang)
%    input nk is the number of points in integral over k
%    output array contains the transfer function
%    NOTE: the psf(r) values are normalized
%
%  Cs   = spherical aberration (in mm)
%  df   = defocus (in Angstroms)
%  kev  = electron energy (in keV)
%  amax = objective aperture  (in mrad)
%
%  started 27-Feb-1997 E. Kirkland
%  modified 2-Mar-1997 ejk
%  modified 9-Apr-1997 David Muller
%  modified 25-May-1998 David Muller  - added finite source size
%           12-Jan-1998 David Muller  - ksource=2 Pi/ dsource 
%
if isfield(params,'d0')
   d0=params.d0;
else
   display('assuming a point source');
   d0=0;
end;


P2 = 2*pi;

Cs = params.Cs*1.0e7;
df = params.df;
kev = params.kev;
amax = params.amax*0.001;
if nargin == 2
   nk=400; % number of points in integral over k
end;  
wav = 12.3986/sqrt((2*511.0+kev)*kev);  % electron wavelength
kmax = amax/wav;
dk = kmax / nk;
k = 0:dk:kmax;
k2 = k .* k;
w1 = 0.5*Cs*wav*wav*wav;
w2 = wav*df;
w = pi*( w1.*k2 - w2 ).*k2;
expw = exp( -i*w );
nr = length( r );
dr = r(nr)/(nr-1);

psf = zeros(nr,1);
% tabulate besselj0 on a grid 5 times finer than nk
bmax = 2*pi*r(nr)*kmax;
db = bmax / (5 *nk);

if d0>0
   b = 0:db:2*bmax;  % will need data out to 2 kmax for MTF
else
   b = 0:db:bmax;
end;

bessj0 = besselj( 0, b );

for ir=1:nr,
	h = expw .* interp1( b,bessj0, P2*r(ir)*k,'*cubic' ) .*k;
	psf(ir) = (abs(sum(h))^2)*2/(nk*pi);
end;
psf(nr)=0;

if d0>0
   %do stuff
   % next inverse psf to get mtf
   k = 0:dk:2*kmax;  % having squared psf, extends frequencies to 2kmax
   nk = length( k );
   for ik=1:nk,
      h = psf' .* interp1( b,bessj0, P2*r*k(ik),'*cubic' ).*r;
	   mtf(ik) = sum(h);
   end;
   % damp by source spread
   mtf = mtf .* exp(-0.5*(k*d0/P2).^2);
   %transform mtf back to psf 
   for ir=1:nr,
	   h = mtf .* interp1( b,bessj0, P2*r(ir)*k,'*cubic' ) .*k;
	   psf(ir) = sum(h)/(nk*pi);
   end;
   psf(nr)=0;
end;    

nrm = 2*pi*sum(psf.*r')*dr;
psf = psf / nrm;