% 
% Written by:
% -- 
% John L. Weatherwax                2005-08-14
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; drawnow; 
clc; 
clear; addpath( '../../Code/eda_data' );

plot_labels = 1; 
use_midden  = 1; % <- use midden=1 or beachdune=0 labeling ... 

% A modified Example 4.2 ... for this problem 
%
load oronsay

tt_str='';
if( plot_labels )
  if( use_midden )
    ind_0 = find(midden==0);
    ind_1 = find(midden==1);
    ind_2 = find(midden==2);
    tt_str = 'midden labeling'; 
  else
    ind_0 = find(beachdune==0);
    ind_1 = find(beachdune==1);
    ind_2 = find(beachdune==2);
    tt_str = 'beachdune labeling'; 
  end
end

x = oronsay;
maxit = 10000;
[n,p] = size(x);
% Set up vector of frequencies.
th = mod(exp(1:p),1);
% This is a small irrational number:
delt = exp(1)^(-5); 
cof = sqrt(2/p);
% Set up storage space for projection vectors.
a = zeros(p,1); b = zeros(p,1);
z = zeros(n,2);
% Get an initial plot, so the tour can be implemented using Handle Graphics.
if( plot_labels )
  h_0=plot(z(ind_0,1),z(ind_0,2),'rd','MarkerSize',8,'MarkerFaceColor','r','erasemode','normal'); hold on; 
  h_1=plot(z(ind_1,1),z(ind_1,2),'gs','MarkerSize',8,'MarkerFaceColor','g','erasemode','normal'); 
  h_2=plot(z(ind_2,1),z(ind_2,2),'bo','MarkerSize',8,'MarkerFaceColor','b','erasemode','normal'); 

  if( use_midden )
    legend([h_0,h_1,h_2], {'midden','Cnoc Coig','Caisteal non Gillean'}, 'location', 'northwest');
  else
    legend([h_0,h_1,h_2], {'midden','beach','dune'}, 'location', 'northwest');
  end
  
else
  ph = plot(z(:,1),z(:,2),'o','erasemode','normal');
end
title( tt_str ); axis equal, axis off
set(gcf,'backingstore','off','renderer','painters','DoubleBuffer','on');
%hold on; 
for t=0:delt:(delt*maxit)
  % Find the transformation vectors.
  for j=1:p/2
    a(2*(j-1)+1)=cof*sin(th(j)*t);
    a(2*j)=cof*cos(th(j)*t);
    b(2*(j-1)+1)=cof*cos(th(j)*t);
    b(2*j)=cof*(-sin(th(j)*t));
  end
  % Project onto the vectors.
  z(:,1)=x*a;
  z(:,2)=x*b;
  
  if( plot_labels )
    set(h_0,'xdata',z(ind_0,1),'ydata',z(ind_0,2));
    set(h_1,'xdata',z(ind_1,1),'ydata',z(ind_1,2));
    set(h_2,'xdata',z(ind_2,1),'ydata',z(ind_2,2));
  else
    set(ph,'xdata',z(:,1),'ydata',z(:,2));
  end

  drawnow
end

opts = struct('FontMode','fixed','FontSize',6,'height',3.5,'Color','bw');
%exportfig(gcf,'Fig42.eps',opts)