%------------------------------------------------------------
% STO-3G fit for Slater-type function with coefficient zeta,
%        the centers are taken to be the origin.
%
% stf: Slater-type function
% gtf: Gaussian-type function
%------------------------------------------------------------
zeta = 1.24;

r   = linspace(-5,5,1000);
stf = sqrt(zeta^3/pi)*exp(-zeta*abs(r));
gtf = zeros(1,length(r));

% STO-3G contraction coefficients and alphas
d = [0.444635 0.535328 0.154329];
a = [0.109818 0.405771  2.22766];
% Scale a to fit zeta != 1
a = zeta^2*a;

for k = 1:3
    gtf = gtf + d(k)*(2*a(k)/pi)^(3/4)*exp(-a(k)*r.^2);
end

% Compare wavefunctions
subplot(211);
plot(r,stf);
hold on;
plot(r,gtf,'r');
xlabel('r (in a.u.)');
ylabel('\phi');
grid on;
title('STO-3G approximation (red) for STF (blue)');

% Compare radial distribution
subplot(212);
plot(r,4*pi*stf.^2.*r.^2);
hold on;
plot(r,4*pi*gtf.^2.*r.^2,'r');
xlabel('r (in a.u.)');
ylabel('Radial distribution: 4*\pi*|\phi|^2');
title('Radial distributions: STO-3G (red) and STF (blue)');
grid on;
