stempot.m

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

function V = stempot(xmax,ymax,nx,ny,potfile)
%  STEMPOT Generate a Projected Potential
%  V = stempot(xmax,ymax,nx,ny,potfile)
%  inputs xmax, ymax are the size of the slice in angstroms
%  nx,ny are the number of pixels in the x and y directions 
%
%  started 10-Apr-1997 David Muller 
%
zed = 1.7;  % zed=2 for rutherford scattering of the nucleus, less for screening

ix = 1:nx;
iy = 1:ny;
dx = xmax/(nx-1);
dy = ymax/(ny-1);
rx = 0:dx:xmax-dx;
ry = 0:dy:ymax-dy;

if isstr( potfile)  
   A = readxyz( potfile);
else
   A = potfile;   % passed a variable instead
end;

Zatom = A.Znum;
ax    = A.xpos;
ay    = A.ypos;
az    = A.zpos;
wt    = A.wt;
amax = length(Zatom);

% find boundaries of slice
axmin = min(ax);
axmax = max(ax);
aymin = min(ay);
aymax = max(ay);
% shift coords to fit in box
ax = ax - axmin;
ay = ay - aymin;


V= zeros(nx,ny);

% map x and y coords of the atoms to the nearest grid points
% A fraction of the atom must be assigned to the closest gridpoints
% to avoid sum and difference frequencies appearing in the image

iax = max(1,min(floor( ax/dx)+1,nx));      % grid point to the left of the atom
%jax = interp1( rx,ix,ax, 'linear');
ibx = mod(iax,nx) + 1;      % create periodic boundary conditions

fax = 1-mod( (ax /dx ),1 );  % fraction of atom at iax 
%fbx = 1 - fax;            % fraction of atom at (1-fax) is at iax+1

iay = max(1,min(floor( ay/dy)+1,ny));      % grid point above the atom
%jay = interp1( ry,iy,ay, 'linear');

iby = mod(iay,ny) +1;          % create periodic boundary conditions

fay = 1-mod( (ay /dy ),1 ); % fraction of atom at iay 
%fby = 1 - fay;            % fraction of atom at (1-fay) is at iay+1

% Add each atom to the potential grid
% j is too large to makegrid(iax,iay) which would allow us to vectorize V
%

V1 = fax .* fay .* (Zatom .^zed);
V2 = (1-fax) .* fay .* (Zatom .^zed);
V3 = fax .* (1-fay) .* (Zatom .^zed);
V4 = (1-fax) .* (1-fay) .* (Zatom .^zed);

for j=1:amax,
   V(iax(j),iay(j)) = V(iax(j),iay(j)) + V1(j);
   V(ibx(j),iay(j)) = V(ibx(j),iay(j)) + V2(j);
   V(iax(j),iby(j)) = V(iax(j),iby(j)) + V3(j);
   V(ibx(j),iby(j)) = V(ibx(j),iby(j)) + V4(j);
end;