%------------------------------------------------------------------% 
% SHO_eigen_C1_L2.m                                                %
% Finite Element Eigensolution of 1-D                              %
% Time Independent Schrodinger Equation (TISE) for SHO             %
%                                                                  %
% By Ed Akin, Rice University, as modification of Example 8.10.2   %
% of Y.W. Kwon "The Finite Element Method Using Matlab", CRC, 1996 %
%                                                                  % 
% Problem description: - Y,xx + (Pot - 2E) Y = 0, with             %
% cubic C1 Hermite polynomials, with Y and Y,x at each of 2 nodes. %
%                                                                  %
% Here: Pot is from the Simple Harmonic Operator (SHO), Pot = x^2  %
%                                                                  %
% Variable descriptions                                            % 
%   p = element potential matrix                                   %
%   k = element stiffness matrix,  kk = system stiffness matrix    % 
%   m = element mass matrix,       mm = system mass matrix         % 
%   index = equation numbers associated with each element          %
%   bcdof = equation numbers associated with boundary conditions   %
%   x_e   = element nodal coordinates                              %
%------------------------------------------------------------------%  
format long        % need high precision
 
%                  USE CONTROL
nel   = 10;        % number of elements (then try 20)
show  = 4;         % number of eigenvalues to show   
list  = 0;         % number of eigenvectors to list
plots = 2;         % number of waveforms, expected values to plot
odd   = 0; %1;     % flag odd or even waveforms
tleng = 5;         % total (half)region length, approximates inf
%                  END CONTROL
 
nnel  = 2;              % number of nodes per element 
ndof  = 2;              % number of dofs per node (for C1)
nnode = (nnel-1)*nel+1; % total number of nodes in system 
sdof  = nnode*ndof;     % total system dofs, = 22 here
leng  = tleng/nel;      % size of TISE elements 
 
% specify Dirichlet bc equation numbers (1 fake eigenvalue each)
if ( odd ) % Odd solutions: Y(0)=0, Y(x_max) = 0
  bcdof (1) = 1;               % Y at node  1 is constrained 
  bcdof (2) = sdof-ndof+1;     % Y at last node is constrained 
  fprintf('SHO for odd modes with %g cubic elements \n', nel)
else       % Even solution: Y,x(0)=0, Y(x_max) = 0
  bcdof (1) = 2;               % Y,x at node  1 is constrained 
  bcdof (2) = sdof-ndof+1;     % Y at last node is constrained 
  fprintf('SHO for even modes with %g cubic elements \n', nel)
end % if bc type
 
kk = zeros (sdof,sdof);      % initialize system stiffness matrix 
mm = zeros (sdof,sdof);      % initialization of system mass matrix 
index = zeros (nel*ndof,1);  % initialization of index vector 

% define the Gaussian quadrature data on -1 to +1
pt = [ -0.906179845938663992797627  -0.538469310105683091036314 ...
        0.000000000000000000000000   0.538469310105683091036314 ...
        0.906179845938663992797627 ];
wt = [  0.23692688505618908751426    0.47862867049936646804129  ...
        0.56888888888888888888889    0.47862867049936646804129  ...
        0.23692688505618908751426 ];
 
for iel = 1:nel    % loop for the total number of elements 
  index = feeldof1 (iel,nnel,ndof); % get element equation numbers
  x_e (1) = leng*(iel-1);           % x left
  x_e (2) = leng*iel;               % x right
  [k,m] = myTISE_C1_L2 (x_e,pt,wt); % compute C1 element matrices
  kk = feasmbl1 (kk,k,index);       % assemble element into system
  mm = feasmbl1 (mm,m,index);       % assemble element into system
end % for iel elements
 
[kn,mn] = feaplycs (kk,mm,bcdof);  % apply the boundary conditions 
 
[Vecs, Diag] = eig (kn,mn);        % get eigenvectors, eigenvalues
for j = 1:sdof                     % extract diagonal matrix
  E (j) = Diag (j,j);              % (fake bc values = 12345)
end % for j
[E, Order] = sort (E);             % reorder eigenvalues, lo to hi

% list the results
if ( show > 0 )
  for j = 1:show                   % list eigenvalues only 
    fprintf ('Eigenvalue %g = %g \n', j, E(j))   % values
  end % for j
end % if show
if ( list > 0 )
  for j = 1:list                   % list eigen value-vector pairs
    fprintf ('Eigenvalue %g = %g \n', j, E(j))   % values
    fprintf ('Eigenvector:')                     % vectors
    disp(Vecs(:,Order(j)))
  end % for j
end % if list

%_____________________________________________________________________ 
% Odd exact results are 2E = 3, 7, 11, 15, 19, ...
% Here nel=10, tleng=5 give 3.00002, 7.00026, 11.0014, 15.0092, ...
% Even exact results are 2E = 1, 5, 9, 13, 17
% Here nel=10, tleng=5 give 1, 5.00009, 9.00062, 13.0033, ...       
% Increasing nel and/or tleng increases accuracy (& problem size)
% ==================================================================

% post-process for <Y|Y>,<Y|H|Y>,<Y|x|Y>, and related graphs
% by rebuilding element arrays since kn, mn destroyed by eig, and
% <Y|Y> = Y'*mn*Y, <Y|H|Y> = Y'*kn*Y also equal the element sums of
% Y_e'*m*Y_e and Y_e'*k*Y_e, etc

if ( plots > 0 ) % then post-process
  x    = [0:(nnode-1)]*leng;    % x-axis for plots
  Y_e  = zeros (nnel*ndof,1);   % pre-allocate element dof
  wave = zeros (sdof,1);        % pre-allocate waveform dof
  YY = zeros (nel,1); YHY = zeros (nel,1) ; YxY = zeros (nel,1) ;
  for p = 1:plots % consider each mode
    clf                                  % clear plot frame
    hold on; grid ;

%   extract this waveform, get expected values from elements
    wave = Vecs(:, Order(p));           % waveform for this mode
    for iel = 1:nel    % loop for the total number of elements
      x_e (1) = leng*(iel-1);           % x left
      x_e (2) = leng*iel;               % x right
      [k,m] = myTISE_C1_L2 (x_e,pt,wt); % re-compute C1 matrices
      [X]   = my_exp_C1_L2 (x_e,pt,wt); % expectation matrices
      index = feeldof1 (iel,nnel,ndof); % element dofs numbers
      Y_e = wave (index);               % el dof via vector subscripts
      YY(iel)  = Y_e'*m*Y_e;            % should sum to unity
      YHY(iel) = Y_e'*k*Y_e;            % expected Hamiltonian
      YxY(iel) = Y_e'*X*Y_e;            % expected x
    end % for iel
    fprintf('Plot eigenvalue %g = %g \n', p, E(p))
    fprintf('Element contributions to <Y|Y>: \n')
    disp(YY)
    Y_Y = sum(YY); Y_H_Y = sum(YHY); Y_x_Y = sum(YxY);
    fprintf('<Y|Y>    = %g \n', Y_Y);         % normalization
    fprintf('<Y|H|Y>  = %g \n', Y_H_Y);       % domain integral
    fprintf('<Y|x|Y>  = %g \n', Y_x_Y);       % domain integral
    fprintf('Expected 2E = %g \n', Y_H_Y/Y_Y); % expected E
    %b fprintf('Expected x = %g \n', Y_x_Y/Y_Y); % SHO expected x == 0

%     scale function value, not slopes too
    ymax = max(max([wave(1:2:sdof),wave(1:2:sdof).^2]));
    axis([0,tleng,min(wave(1:2:sdof)),ymax])
    if ( odd )
      title(['SHO Piecewise cubic waveform for odd mode ', ...
             int2str(p), ',  2E = ', num2str(E(p))])
    else
      title(['SHO Piecewise cubic waveform for even mode ', ...
             int2str(p), ',  2E = ', num2str(E(p))])
    end % if even-odd
    xlabel(['X-axis (atomic units), for ', int2str(nel), ' elements'])
    ylabel('Node Value - o  (solid),   Y ^2 (dashed)')
    %b plot (x,wave(1:2:sdof),'ko'); 
    %b plot (x,wave(1:2:sdof),'k-');                     % waveform
    %b plot (x,wave(1:2:sdof).^2,'k--');                % Y^2
    opt = 3;                              % plot choice
    plot_cubic (nel,x,wave,opt);          % opt 1-Y 2-Y,x 3-Y & Y^2

%   print hardcopy, pause, go to next mode
    print ('-dpsc', ['waveform_', int2str(p), '_graph'])
    hold off;
    v_text = ['Created waveform_', int2str(p), '_graph.ps'];
    fprintf (1,'%s', v_text) ; fprintf (1, ' \n' )
    pause (4) ;
%   note: could plot piecewise constant YY here for each el over x
  end % for p plots
end % if post-proces
%_____________________________________________________________________ 
%  Related references:
% 1.  J.E. Akin, "Finite Element Analysis with Error Estimates"
%  www.owlnet.rice,edu/~mech517  (under Books link)
% 2. L.R. Ram-Mohan, "Finite Element and Boundary Element Applications
%  in Quantum Mechanics", Oxford Univ. Press, 2002
% 3. W. Schweizer, "Numerical Quantum Dynamics", Kluwer, 2001
% end SHO_eigen_C1_L2
