function E = pset4problem2b(M,I) %calculates the energies of the system given the size of the basis function % written by Sandeep Sharma on 26th Sep 2006 %Taken input % M which is the size of the basis function to be used % I the angular moment of the hindered rotor %outputs the vector of energies. %Start by defining all the constants h = 6.625e-34; %J-s %get the coefficients ys of the potential function %y = [2.0667 -2.0333 1.2333 -1.7667 0.5]'*1e-19; y = [2.067 -2.033 1.056 -1.767 0.678]'*1e-20; %fill the elements due to the differential operators for n=1:2*M+1 H(n,n) = h^2*(n-M-1)^2/8/pi^2/I; end %fill in the elements for on the upper triagular side and diagonal of the matrix. for n=1:2*M-3 H(n,n) = H(n,n) + y(1); H(n,n+1) = y(2)/2; H(n,n+2) = y(3)/2; H(n,n+3) = y(4)/2; H(n,n+4) = y(5)/2; end H(2*M-2,2*M-2) = H(2*M-2,2*M-2) + y(1); H(2*M-2,2*M-1) = y(2)/2; H(2*M-2,2*M) = y(3)/2; H(2*M-2,2*M+1) = y(4)/2; H(2*M-1,2*M-1) = H(2*M-1,2*M-1) + y(1); H(2*M-1,2*M) = y(2)/2; H(2*M-1,2*M+1) = y(3)/2; H(2*M,2*M) = H(2*M,2*M) + y(1); H(2*M,2*M+1) = y(2)/2; H(2*M+1,2*M+1) = H(2*M+1,2*M+1) + y(1); %fill in the elements on the lower triangular part of the matrix. for n= 5:2*M+1 H(n,n-1) = y(2)/2; H(n,n-2) = y(3)/2; H(n,n-3) = y(4)/2; H(n,n-4) = y(5)/2; end H(2,1) = y(2)/2; H(3,2) = y(2)/2; H(3,1) = y(3)/2; H(4,3) = y(2)/2; H(4,2) = y(3)/2; H(4,1) = y(4)/2; %The eigenvalues of matrix H are the allowed energy states of the system E = eig(H); %Display the zero-point energy of the system. output = strcat('The zero point energy of the system = ',num2str(E(1)),' J\n'); fprintf(output); phi = linspace(-0.1, 2.1*pi); plotV(y,phi); hold; for n=1:length(E) x = [0; 2*pi]; y = [E(n),E(n)]; plot(x,y,'r'); end title('Eigenvalues of hamiltonian'); xlabel('Angle (radians)'); ylabel('Energy (J/molecule)'); return; %*********************************************************** function plotV(y,phi) V = zeros(length(phi),1); for n=1:length(phi) for m = 1:length(y) V(n) = V(n) + y(m)*cos((m-1)*phi(n)) ; end end plot(phi,V); return;