clear
% Testing Hermite polynomial construction
x = -10:.1:10;
n=2;
plot(x,hermiteH(n,x)+10,x,polyval(hermite_rec(n),x))

%

% Harmonic Oscillator Eigenfunctions
clear
n=0;             % Energy level n>0

a=0.015;         % x-range
x = linspace(-a,a,100);

omega=10;
m=9.11*10^(-31);
hbar=1.05457*10^(-34);
beta = sqrt(m*omega/hbar);


nf=factorial(n);
N = 1/sqrt(2^n*nf*beta*sqrt(pi));
psi1=N*hermiteH(n,beta*x).*exp(-beta^2*x.^2/2);
Hnx=polyval(hermite_rec(n),beta*x);
psi2=N*Hnx.*exp(-beta^2*x.^2/2);        % Comparing different Hn computations
plot(x,psi1,x,psi2+.001)

%

% Create Quantum Oscillator Coherent States
clear
a=0.015;
x = linspace(-a,a,100);

omega=10;
m=9.11*10^(-31);
hbar=1.05457*10^(-34);
beta = sqrt(m*omega/hbar);

alpha0 = 1;                   % |alpha|

T = 4*(2*pi/omega);           % Creates 4 cycles of oscillation
M = 100;                      % Number of time steps, # frames = M+1
dt=T/M;

figure(1)

% Initialize the movie
%s = settings;
%s.matlab.editor.AllowFigureAnimation.PersonalValue = true;

nFrames = M+1;
Mov(1:nFrames) = struct('cdata',[], 'colormap',[]);
frame=1;

for j=0:M
    t=j*dt;
    alpha=alpha0*exp(i*omega*t);
    ealpha = exp(-alpha*conj(alpha)/2);

    tmp=0;
    %figure('units','pixels','position',[0 0 1440 1080]) % for video
    %creation
    % Coherent state = sum over oscillator states
    for n=0:10
        nf=factorial(n);
        N = sqrt(beta)/sqrt(2^n*nf*sqrt(pi));
        c = alpha^n/sqrt(nf)*ealpha;
        %Hnx=polyval(hermite_rec(n),beta*x);
        %tmp=tmp+c*N*Hnx.*exp(-beta^2*x.^2/2);  % Either line works
        tmp=tmp+c*N*hermiteH(n,beta*x).*exp(-beta^2*x.^2/2);  
    end
    plot(x,tmp.*conj(tmp))

    Mov(j+1)=getframe(gcf);    
end
hold off

% To see movie in Command Window type:   movie(gcf, Mov,1,15)
% movie(M,n,fps) specifies to play n times at fps frames per second.

%

% DumpToFile(Mov,nFrames,'Coherent4',' ') 
% DumpToFile(Mov,nFrames,'Coherent3','Uncompressed AVI')   % To create video
% DumpToFile(Mov,nFrames,'Coherent2','MPEG-4')   % To create video

%

function h = hermite_rec(n)
% Recursive construction of Hermite Polynomial coefficients
if( 0==n ), h = 1;
elseif( 1==n ), h = [2 0];
else
   h1 = zeros(1,n+1);
   h1(1:n) = 2*hermite_rec(n-1);

   h2 = zeros(1,n+1);
   h2(3:end) = 2*(n-1)*hermite_rec(n-2);

   h = h1 - h2;
end
end

function y=hermiteH(n,x)
% One approach to construction of Hermite polynomials from hermite_rec
% Needed to add n=0 case to bypass hermite_rec.
if n==0
    y=ones(size(x));
    return
end
h=hermite_rec(n);
y = h(end);
p = 1;
for i=length(h)-1:-1:1
    y = y + h(i) * x(:).^p;
    p = p+1;
end
y = reshape(y,size(x));
end

function DumpToFile(M,nframes,filename,type)
if (type==' ')
    filename=[filename,'.avi'];
    v = VideoWriter(filename);
else
    v = VideoWriter(filename,type);
end

open(v)
for k=1:nframes
     frame = M(k);
     writeVideo(v,frame)
end
close(v)
end