% ANNgonplot.m
% m-file for ANN to plot closed reducing n-gons
% RWN 05/20/07
clear
hold off
N=input('enter number of nodes: ');
if N== [ ]
    N=5
end
P=input('enter column of the N initial vertices as xi+yi*j: ');
if P == [ ]
    for k=1:N
        P(k)=exp((2*pi*k/N)*sqrt(-1))
%        x(k)=cos(2*pi*k/N); y(k)=sin(2*pi*k/N);
%       P(k)=x(k)+y(k)*j;
% for N=3 or 5 use the following
%	P=[1+1j 2+1j 1.5+0.5j]'
%   P=[0.3090+0.9511i;-0.8090+0.5878i;-0.8090-0.5878i;0.3090+0.9511i;-0.8090+0.5878i;-0.8090-0.5878i]
    end
end
Prev=input('enter 1 if reverse nodes 2 and 3 otherwise return: ');
if Prev == [ ]
    P=P
    else 
        P2rev=P(3); P3rev=P(2);
        P(2)=P2rev; P(3)=P3rev;
end 
% input for scaling by m at mth time
Sc=input('enter 1 if desire scaling, else return for non-scaled: ')
% if do not set N above use the following
% set N = # of entries in P
% N=length(P)
T=input('enter the final time T: ');
if T == [ ]
	T=3
end
% T=0 gives the initial graph
% form Z(0)=P, Z={P zeros(N,T%} for inserting P followed by zeros
%Z={P zeros(N,T)} 
% form W
for j=1:N
    for k=1:N
        W(j,k)=0;
    end
    W(j,j)=-2;
end
for j=2:N-1
    W(j,j-1)=1;
    W(j,j+1)=1;
end
W(1,2)=1; W(1,N)=1; W(N,1)=1; W(N,N-1)=1;
W
% form neural network
    net1=network;
    net1.numInputs=1;
    net1.numLayers=1;
    net1.inputConnect(1,1)=[1];
    net1.outputConnect=[1];
    net1.inputs{1}.range=[-2 2];
    net1.inputs{1}.size=N;
    net1.layers{1}.size=N;
    net1.layers{1}.transferFcn='purelin';
    net1.IW{1,1}=W;
% in order to conveniently plot do the simulation via a matlab for 
% command, for m=0:1:T, instead of via delay feedback, so don't need 
% the following two commands
%
%    net1.layerConnect(1,1)=[1];
%    net1.layerWeights{1,1}.delays=[1];
%%%%%% set up simulation and then plot, with scaling option due to Sc
for m=0:1:T
    if m==0
        z=P'
        z=sim(net1,z)
        z1=[z.', z(1)].'
        z2=z1./norm(z1)
        plot(z2,'g-'); 
        grid on
        hold on;
        
    else
% in order to plot all curves at once use sim(net1,z) rather than delay
         z=sim(net1,z)  
%        z=sim(net1,zeros(N))
           z1=[z.',z(1)].'
         if Sc==1
           z2=z1./(m*norm(z1))
         else
           z2=z1
         end
         if 1==mod(m,2)
            plot(z2,'r-')
            hold on;
         else 
            plot(z2,'k-')
         end
    end
end
% note that for more than 60 nodes a distorted picture results
% for what reason it is not yet clear