ctemconvolution.m

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

function II = ctemconvolution(V,rmax,nx,params,type)
% MATLAB function to convolve the CTEM PSF with the projected potential
%
% function II = ctemconvolution(V,rmax,nx,params,type)
%
%  input matrix V is the projected potential
%
%  input rmax = radial limit
%  input nx = dimension of V. i.e. double V(nx,nx);
%
%  input type = 0 for phase contrast
%        and 1 for amplitude contrast
%
%  params ={ Cs, df, kev, ddf, beta, amax}
%  params.Cs   = spherical aberration (in mm)
%  params.df   = defcocus (in Angstroms)
%  params.kev  = electron energy (in keV)
%  params.ddf  = chromatic aberation defocus spread (in Angstroms)
%  params.beta = spread in illumination angles (in mrad)
%
% started 07-Jun-06 Huolin Xin
% modified 11-Jun-06 Huolin Xin
% modified 06-Jul-06 Huolin Xin (replace step function by
% 1/(1+exp(A*(x-xc)));


Cs = params.Cs*1.0e7;
df = params.df;
kev = params.kev;
ddf = params.ddf;
beta = params.beta*0.001;
amax = params.amax*0.001;

wav = 12.3986/sqrt((2*511.0+kev)*kev);
klimit = amax/wav;

V = fft2(V);

dx = rmax/nx;
kmax = 1/dx;
dk = 1/(dx*nx);
[kx ky] = meshgrid(-kmax/2:dk:kmax/2-dk);
k2 = fftshift(sqrt(kx.^2+ky.^2));
clear kx;
clear ky;
H2d = ctemincohmtf(k2,params,type);

klimit = min(klimit,kmax/2);

%replace step by 1/(1+exp(A*(x-xc)))
x = 0:dk/4:kmax+dk;
y = 1./(1+exp(1000.*(x-klimit)));

k2 = interp1(x,y,k2,'linear');

clear x;
clear y;

V = V.*H2d;
clear H2d;
V = V.*k2;
clear k2;

II = real(ifft2(V));