ctemhr.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)
%

P2 = 2*pi;

Cs = params(1)*1.0e7;
df = params(2);
kev = params(3);
ddf = params(4);
beta = params(5)*0.001;
mo = 511.0;      % electron rest mass in keV
hc = 12.3986;    % in keV-Angstroms
wav = (2*mo)+kev;
wav = hc/sqrt(wav*kev);
w1 = pi*Cs*wav*wav*wav;
w2 = pi*wav*df;
e0 = (pi*beta*ddf)^2;
k2 = k .* k;
wr = (w1.*k2-w2).*k*beta/wav;
wi = pi*wav*ddf.*k2;
wi = wr.*wr + 0.25.*wi.*wi;
wi = exp(-wi./(1+e0.*k2));
wr = ( 0.5*w1.*k2.*(1-e0.*k2) - w2 ).*k2./(1+e0.*k2);
if nargin == 2
   nk=400; % number of points in integral over k
end;  

expw = wi;
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);
b = 0:db:bmax;
bessj0 = besselj( 0, b );

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

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