%
%  MATLAB script to plot CTEM PSF
%
%  last modified 11-1-02 DAM
%
%  Cs = Spherical Aberration
%  df = defocus
%  kev = electron energy in keV
%  ddf = chromatic aberation defocus spread
%  beta =  spread in illumination angles
%
disp( 'Plot CTEM transfer function' );
kev  = input( 'Type electron energy in keV :');
Cs   = input( 'Type spherical aberation Cs in mm :');
amax = input( 'Type objective aperture semiangle in mrad :');
df   = input( 'Type defocus df in Angstroms :');
ddf  = input( 'Type defocus spread ddf in Angstroms :');
beta = input( 'Type illumination semiangle in mrad :');

%
P2 = 2*pi;
nk=500;
wav = 12.3986/sqrt((2*511.0+kev)*kev);  % electron wavelength
rmax = 6*sqrt( sqrt( Cs*1.0e7*wav*wav*wav ));
r = 0:(rmax/600):rmax;  
param.Cs   = Cs;
param.df   = df;
param.kev  = kev;
param.amax = amax;
param.ddf = ddf;
param.beta = beta;

%--------------------------------
amax = param.amax*0.001;
kmax = amax/wav;
dk = kmax / nk;
k = 0:dk:kmax;
params = [ Cs, df, kev, ddf, beta ];
sinw = ctemh( k, params, 0 );
cosw = ctemh( k, params, 1 );

%--------------------------------

expw = cosw+ sqrt(-1) * sinw;
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;
%--------------------------------

hold off;
plot( r, real(psf), r, imag(psf) );
legend('Real','Imag');
axis([0, rmax, -1, +1]);
xlabel( 'Distance (in A)');
ylabel( 'PSF' );
s1 = sprintf('E= %gkeV, Cs= %gmm, df= %gA,  ', kev, Cs, df);
s2 = sprintf('Beta= %gmrad, ddf= %gA', beta, ddf);
title([s1 s2]);
hold on;		 % plot line through zero
x = [0, rmax];
y = [0, 0];
plot( x, y );
hold off;