%% exam 1, 2017, matlab portion
%% section 1a


syms r re alpha k omega mu
assume(alpha, 'positive')

psi0=(alpha/pi)^(1/4)*exp(-alpha*(r-re).^2/2);
psi1=(alpha/pi)^(1/4)*sqrt(2)*sqrt(alpha)*(r-re)*exp(-alpha*(r-re).^2/2);



normint=int(psi0.*psi0,r,-inf,inf);

if ( abs(normint-1)<=1.e-12)
    disp('psi0 is normalized')
    
else
    disp(['normalization integral = ' num2str(normint) '.NE. 1'])
end

%% section 1b

% integrate and multiply by -1/2

% make it easier set re=0, doesn't really matter since we're integrating
% from -inf to +inf

integrand=subs(psi1*psi1*(r-re)^2,re,0);
V11=(1/2)*k*int(integrand,r,-inf,inf);

integrand=subs(psi1*diff(psi1,r,2),re,0);
T11=(int(integrand,r,-inf,inf))*(-1/(2*mu));

% the result is that V11 =T11 (virial theorem)

% we know that k = alpha*omega and alpha = omega*mu


V11=subs(V11,k,alpha*omega);
T11=subs(T11,alpha,omega*mu);
disp(['V11, T11: ' char(V11) ' ' char(T11)])

% each is 1/2 the total energy.   The simplest way of calculating the
% kinetic energy is T=E-<V>.   We know that E_1=3/2 omega, and V11=3/4
% omega, so that T=3/4 omega

%%  section 1c
% % for H2
   mu=918;
   k=0.369;
   rre=1.4;  % call the numerical value rre to distinguish from the symbolic variable re
   alph= sqrt(k*mu);  % call the numerical value alph to distinguish from the symbolic variable alpha




% plot the wavefunctions

ppsi0=subs(psi0,[alpha,re],[alph,rre]);
ppsi1=subs(psi1,[alpha,re],[alph,rre]);
ezplot(ppsi0,[0.5 2.5]);hold;ezplot(ppsi1,[0.5 2.5]);hold
% label axes
xlabel('r');ylabel('\psi');title('\psi_0 and \psi_1')

