stemhk.m

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

function mtf = stemhk(k,params)
%
%  MATLAB function to calculate STEM mtf vs. k
%    input array k has the spatial freq. (in inv. Angs.)
%    input array params has the optical parameters
%          params = [Cs, df, kev, amax]
%    output array contains the transfer function
%
%  Cs   = spherical aberration (in mm)
%  df   = defocus (in Angstroms)
%  kev  = electron energy (in keV)
%  amax = objective aperture  (in mrad)
%  d0   = source size (in A)
%
%  started 6-apr-1997 E. Kirkland
%  last modified 6-apr-1997 ejk
%           12-Jan-1998 David Muller  - ksource=2 Pi/ dsource
%

P2=2*pi;

Cs = params.Cs*1.0e7;
df = params.df;
kev = params.kev;
amax = params.amax*0.001;
if isfield(params,'d0')
   d0=params.d0;
else
   d0=0;
end;

%  first calculate the psf using stemhr()
nr = 400;  % number of points in integral over r
wav = 12.3986/sqrt((2*511.0+kev)*kev);  % electron wavelength
rmax = max(2*sqrt( sqrt( Cs*wav*wav*wav) ),3);
r = 0:(rmax/nr):rmax;

% stemhr will include source size otherwise (and transform back) -quicker to do here
param.d0=0; 
psf = stemhr( r, params,nr );
param.d0=d0; % set it back just in case anyone else needs d0

% next inverse psf to get mtf
nk = length( k );
for ik=1:nk,
	h = psf' .* besselj( 0, P2*r*k(ik) ) .*r;
	mtf(ik) = sum(h);
end;

if d0>0
   %Damp by source size
   display('source size not yet tested');
   mtf = mtf .* exp(-0.5*(k*d0).^2);  %should I divide by 2Pi?
else
   display('assuming a point source');
end;

a = mtf(1);
mtf = mtf/a;  % normalize mtf(0)=1