function eof_demo_movie(type,AMP,Modes)
% same as eof_demo_movie, but renormalizes PC and EOF to have unit variance
% in PC and amplitude in the EOF.

% Demonstrate how to get eofs, pcs, etc
% try all four idealized functions (type 1,2,3,4, 5 and 6)
% and try no noise (AMP =0) up to lots of noise (AMP = 1)
% Modes = the number of eigenvectors (EOFs) you want to use to reconstruct
% the data.

% NOTE: Modes = [1 2 4] would use three eigenmodes: the two with the most
% variance and the fourth, with the fourth most variance


%type =1; AMP =0.2
x=[0:0.1:1];
t=[0:.1:20];
N=length(t);
M=length(x);
xx = repmat(x',1,N);
tt = repmat(t,M,1);
xx = repmat(x',1,N);

Modes
mmodes = M+1 - Modes;

% AMP = 0.2
noise=AMP*randn(M,N);

if type ==1 
    z=cos(pi*xx).* cos(-2*pi*tt) + noise;  % one standing wave
end

if type ==2 
    z=cos(2*pi*xx).* cos(-2*pi*tt) + noise;  % one standing wave
end

if type ==3
    z=cos(2*pi*(xx - tt/2)) + noise;  % one propagating wave
end

if type ==4
    z=randn(M,N);  % uncorrelated data
end

if type ==5
    z=repmat(0,M,N);  % uncorrelated data
    for ij = 1:3
        z(ij,:)=sin(pi*t/ij);
    end
    z=z+noise;
end

if type ==6
    z=repmat(0,M,N);  % uncorrelated data
    for ij = 1:3
        z(ij,:)=0.3*ij*sin(pi*t/ij);
    end
    z=z+noise;
end


% take out mean
zmean=mean(z');
z= z -repmat(zmean',1,length(t));


% F = [];
figure(1)
clf

for j=1:40
    subplot(2,2,1)
    if j == 2; pause;end
    plot(x,z(:,j)');hold on; plot([0 1],[0 0],'k:');hold off;
    xlabel('Location');ylabel('Amplitude')
    set(gca,'FontSize',16),title(['Time = ' num2str(t(j))])
    set(gca,'FontSize',16,'LineWidth',3)
    axis([0 1 -2 2])
     if j >= 2
         subplot(1,2,2)
   %      pcolor(z(:,1:j)');caxis([-1.5 1.5]),colormap(jet),colorbar
            pcolor(x,t(1:j),z(:,1:j)');caxis([-1.5 1.5]),colormap(jet),colorbar
            xlabel('Location');ylabel('Time')
          set(gca,'FontSize',16),title(['Data Matrix'])
     end
    pause(0.3)
%    F(j)=getframe
end
    pause
% movie(F,2,1)

% subplot(2,1,1)
% pcolor(z);
% colormap(jet);colorbar
% title(['Actual Field with AMP =' num2str(AMP) ])

% get covariance matrix c

c = z*z'/N;

% get eigenvalues and eigenvectors

[v, sig]=eig(c);


% check on variance
So = trace(c);
S =trace(sig);
Sig = diag(sig);



% get principle components

     P = v'*z;
   
% reconstruct with mmodes vectors
 zt =   v(:,mmodes)* P(mmodes,:);

% plot the data and the  vector reconstruction

figure(2) 
clf
 
subplot(1,2,1)
%pcolor(z');caxis([-1.5 1.5]),colormap(jet),colorbar
pcolor(z(:,1:40)');caxis([-1.5 1.5]),colormap(jet),colorbar
xlabel('Location'),ylabel('Time')

title(['Actual Field with AMP =' num2str(AMP) ])

subplot(1,2,2)
%pcolor(zt');caxis([-1.5 1.5]),colormap(jet),colorbar
pcolor(zt(:,1:40)');caxis([-1.5 1.5]),colormap(jet),colorbar
xlabel('Location'),ylabel('Time')
title(['Reconstructed using ' int2str(Modes) ' modes'])

pause

% plot the time series of the first two pcs in the top plot and all the pcs
% (bottom plot)


figure(3)

subplot(2,1,1)
  plot(fliplr(P(M-1:M,:)'));hold on; plot([0 N],[0 0],'k:');hold off
    xlabel('Time')
  title('First two PCs')
  subplot(2,1,2)
  plot(fliplr(P'));hold on; plot([0 N],[0 0],'k:');hold off
  xlabel('Time')
  title('ALL PCs')

  
pause

% plot eigenvalues (note: matlab orders the pcs and eofs in increasing
% sigma (so the matrix position (:,M) is eof,PC = 1


% plot eigenvectors (rescaled to have  (pc) with unit variance

ss = sqrt(abs(Sig));  % the square root of the eigenvalues
for ij=1:M
    vs(:,ij) = v(:,ij)*ss(ij);
end

figure(4) 
subplot(2,1,1)
plot(Sig/trace(sig),'o');set(gca,'FontSize',16); axis([0 M+1 0 1.2])
grid on;
title('Fraction of total variance in each eigenmode')
subplot(2,1,2)
plot(x,fliplr(vs(:,M-2:M)));set(gca,'FontSize',16)
hold on; plot([0 1],[0 0],'k:');hold off; axis([ x(1) x(M) -Inf Inf])
xlabel('Location');ylabel('Amplitude')
title('First three eigenvectors (scaled)')

% corrcoef(z')

end


