% solutions to problems in Matlab introduction, 2017

%% problem 1

mu = 918; ome = 0.02; re=1.4; % parameters defining harmonic oscillator function

k= ome^2*mu;  % calculate k, since ome=sqrt(k/mu).  
      %  note that k must be positive

alpha = sqrt(k*mu);    % then calculate alpha

% output values
display([' for H2 harmonic oscillator k = ' num2str(k) '; alpha = ' num2str(alpha)])

% wavefunction then is at r=1.5
r=1.5;
psi = (alpha/pi)^(1/4)*exp(-alpha*(r-re)^2/2);
%% problem 2

r=[.8:.02:1.9];          % grid of r values

for i=1:length(r)        % note the 'length', to determine the number of
                         % values of r
  ppsi(i)=(alpha/pi)^(1/4)*exp(-alpha*(r(i)-re)^2/2);
  endt

% Note that this could be done with a single command
ppsi=(alpha/pi)^(1/4)*exp(-alpha*(r-re).^2/2);
% this creates a vector ppsi using the exponentiation of a vector (r-re)
% note that i need to put in a period before the ^2.  This tells Matlab to
% square the vector element-by-element.  otherwise Matlab doesn't know how
% to carry out the ^2 operation on a vector.

plot(r,ppsi); xlabel('r / bohr'); ylabel('\psi_0(r)')    % use latex notation
                     % for greek variables and subscripts
                     
% make font size larger
set(gca,'fontsize',14)   % get(gca) will indicate all attributes you can change

%% problem 3

% setting symbolic variables

syms alpha x
assume(alpha,'positive')  % alpha was previously defined in line 10 as a numeric value
              % this will redefine it

psi0=(alpha/pi)^(1/4)*exp(-alpha*x^2/2);   % define a symbolic function for psi0

int(psi0*psi0,-inf,inf)       % this will give 1 as it should
% now, v=2 wavefunction

psi2 = (alpha/pi)^(1/4)*((2*alpha*x^2-1)/sqrt(2))*exp(-alpha*x^2/2);  

% is this normalized?
int(psi2*psi2,-inf,inf)  % yes, it is normalized

% is it orthogonal to psi0?

int(psi2*psi0,-inf,inf)  % yes, <2|0> = 0
% let's examine matrix elements

% we'll need another symbolic variable k
syms k

V=k*x^2/2;

V00= int(psi0*V*psi0,-inf,inf);

% the answer is k/(4*alpha)
% to get rid of alpha, we'll need to use the relation alpha =sqrt(k*mu)
% we'll need to redefine mu as a symbolic variable

syms mu

V00=subs(V00,alpha,sqrt(k*mu));   % this has units of sqrt(k/mu) = omega, as we want, 
% the numerical answer is omega/4.   This is reasonable.  The energy of the
% n=0 level is omega/2.  For the harmonic oscillator 1/2 of the energy is
% kinetic and 1/2 potential.


% now kinetic energy

% call the effect of T on psi0 = Tpsi0

Tpsi0=-(1/(2*mu))*diff(psi0,x,2);

% so <0|T|0> is

T00 = int(psi0*Tpsi0,-inf,inf);  % the answer is alpha*mu/4

T00=subs(T00,alpha,sqrt(k*mu));


% this is equal to sqrt(k/mu)/4  as we want


% output the answer
display(['V00 and T00 for the harmonic oscillator:'  char(V00) ' ' char(T00)])  
    % use char() to display the symbolic variables V00 and T00!
    
    
 %% problem 4

ppsi0=@(alpha,xx)(alpha/pi)^(1/4)*exp(-alpha*xx.*xx/2);

t=[0:.01:1];


ppsi2=@(alpha,xx)(alpha/pi)^(1/4)*(1/sqrt(2))*(2*alpha*xx.^2-1).*exp(-alpha*xx.*xx/2);
% note the period in the xx.^2, and in the ).*exp, just like in the xx.*xx term, to force
% calculation of a vector
plot(t,ppsi0(20,t),t,ppsi2(20,t))  % this plots both on the same plot

% to see the root, add a grid to the plot
grid

% use ginput(1) and put the cursor on the root, to get x_root=0.156


xxx=[0 .1 .2 .3];  % guess that these straddle the root
cc=polyfit(xxx,ppsi2(20,xxx),3);  % fit these points with a 3rd degree polynomial

% then evaluate the zeros
rr=roots(cc)  % this gives 4.4023e-01   1.5843e-01  -1.1054e-01
% the only one that's in the range is 1.58433.  To test that it's close, we
% can evaluate the function at that point
ppsi2(20,rr(2));  % this give 3.54e-3, close to zero

%% problem 5

syms cs sn h11 h12 h22
C=[cs sn;-sn cs];
H=[h11 h12;h12 h22];
assume(cs,'real');
assume(sn,'real');
E=C'*H*C ;

% we want E to be a diagonal matrix, so E(1,2) must equal zero
% this element is
expand(E(1,2))  % cs*sn(h11-h22)+(cs^2-sn^2)h12
% now cs*sn = sin(2*theta)/2  and cs^2-sn^2 = cos(2*theta)
% so we have:   E(1,2)=h12*cos(2*theta) - sin(2*theta)*(h11-h22)/2
% this must equal zero so, sin(2*theta)/cos(2*theta) = 2*h12/(h22-h11)
% or theta = 0.5*atan(2*h12/(h22-h11)

% now the E matrix is symmetric, which we can show by
expand(E(1,2)-E(2,1)) % this equals zero

% now, what about the diagonal matrix elements
expand(E(1,1));  %this gives h11*cs^2+h22*sn^2 - 2*h12*cs*sn
% having determined theta, we can get cos(theta) and sin(theta), and, from
% these E(1,1), as well as E(2,2), which equals
% h11*sn^2+h22*cs^2+2*h12*cs*sn




