with(plots):
with(plottools):
with(linalg):


DimensionString:=`height=150,width=150`;
# dimensions of jpeg pictures
FullPath:=``;
#FullPath:=`c:/keller/maple/movie/`;
SequenceOfPictures:=NULL;



towindow:=proc();
	plotsetup(win,window);
	display(SequenceOfPictures,insequence=true,axes=none, view=[-1.2..1.2,-1.2..1.2])
end;

addzeroes:=proc(n,nberdigits) local prov;
	prov:=``.n;
	while (length(prov)< nberdigits) do prov:=`0`.prov od;
	prov
end;

tojpeg:=proc(root) 
# root : root of filename
local filename,i,nber;
global SequenceOfPictures, FullPath;
for i from 1 to nops([SequenceOfPictures]) do
    nber:=addzeroes(i-1,3);
    filename:=``.FullPath.root.nber.`.jpg`;
    print(`Writing `.filename);
    plotsetup(jpeg,plotoutput=filename,plotoptions=DimensionString);
   print(display(SequenceOfPictures[i],axes=none, scaling=constrained, view=[-1.2..1.2,-1.2..1.2]));
 od;
 plotsetup(win,window);
end;



ressort:=proc(a,b,r,delta,n) local l,i;
   l:=[r,a],seq([r+delta*(-1)^i, a+i*(b-a)/n],i=1..(n-1)),[r,b];
   plot([l],coords=polar,color=red);
end;

mydisk:=proc(r,a,delta) ;
   disk([r*cos(a), r*sin(a)], delta,color=red)
end;   

Constantsquare:=polygon([[0.9,-0.1],[1.1,-0.1],[1.1,0.1],[0.9,0.1]]);

mysquare:=proc(a);
   rotate(Constantsquare,a)
end;

myconfig:=proc(l) local i,r,s,d;
   s:=NULL;
   for i from 1 to nops(l) do
      if i<nops(l) then d:=l[i+1] else d:=l[1] fi;
#      print(i,l[i],d);
      while d<l[i] do d:=d+evalf(2*Pi) od;
      r:=ressort(l[i],d,1,0.1,20);
      s:=s,mydisk(1,l[i],0.1),r;
#	  s:=s,mysquare(l[i]);
   od;
  display({s});
end;
       


suiteetatpur:=proc(n,k,sincos,framenumber)

# n : number of masses on the circle (e.g. 4)
# k : number of the pure state (between 1 and n, e.g. 1)
# sincos : 0 -> sinus-speed vector, 1 cos-speed vector (e.g. 0)
# framenumber : of animation (e.g. 40)

local inipos,inispeed, i,j,om,E,c,T,x,s,l;
global SequenceOfPictures;

 inipos:=vector(n,0);
 if sincos=1 then 
   inispeed:=[seq(sin(k*2*Pi*j/n), j=0..n-1)];
 else
   inispeed:=[seq(cos(k*2*Pi*j/n), j=0..n-1)];
 fi;
 E:=2*sqrt(add(inispeed[i]^2,i=1..n));
 for i from 1 to n do inispeed[i]:=inispeed[i]/E od;
 om:=sqrt(1-cos(k*2*Pi/n));
 c:=2*Pi/n;
 if k<>0 then T:=evalf(2*Pi/om) 
	else 
	 T:=evalf(2*Pi/(n*inispeed[1])) 
 fi;
 SequenceOfPictures:=NULL;
 for x from 0 to T-T/(2*framenumber) by T/framenumber do
  if k<>0 then s:=sin(om*x) else s:=x fi;
  l:=[seq(evalf((i-1)*c+s*inispeed[i]),i=1..n), evalf(s*inispeed[1])] ; 
  SequenceOfPictures:=SequenceOfPictures, myconfig(l)
 od;
 end;


etatpur:=proc(n,k,sincos,framenumber)

# n : number of masses on the circle (e.g. 4)
# k : number of the pure state (between 1 and n, e.g. 1)
# sincos : 0 -> sinus-speed vector, 1 cos-speed vector (e.g. 0)
# framenumber : of animation (e.g. 40)
  
global SequenceOfPictures;
suiteetatpur(n,k,sincos,framenumber);
towindow()
end;


sample:=proc(n,framenumber)
local k,sincos;
for k from 0 to floor(n/2) do
	if (k=0) or (k=n/2) then sincos:=1 else sincos:=0 fi;
        suiteetatpur(n,k,sincos,framenumber);
	tojpeg(`m`.n.k);
od;
end;

suiteetatmixte:=proc(n,p,v,framenumber)

# n : number of masses on the circle (e.g. 4)
# p : n-vector of initial positions
# v : n-vector of initial velocities
# framenumber : of animation (e.g. 40)

local C, id, A, inipos, inispeed, inicond, functlist,
eq, sol, i, r, c, x, f, g, j, l,l1;
global prespos, presvel, SequenceOfPictures;


C:=companion(X^n-1,X):
id:=band([1], n):
A:=evalm(-1/2*(C+inverse(C))+id);

inipos:=vector(n,p);
inispeed:=vector(n,v);
 
inicond:=seq(y[i](0)=inipos[i],i=1..n), seq(D(y[i])(0)=inispeed[i],i=1..n):
functlist:=seq(y[i](t),i=1..n);
eq:=seq(diff(y[i](t),t$2)+add(A[i,j]*y[j](t),j=1..n)=0,i=1..n):

print(inicond);
sol:=dsolve({eq,inicond},{functlist},type=numeric, output=listprocedure);

for i from 1 to n do f[i]:=subs(sol, y[i](t)); od:
for i from 1 to n do g[i]:=subs(sol, diff(y[i](t),t)); od:
c:=2*Pi/n: r:=1:
SequenceOfPictures:=NULL;
for x from 0 to framenumber*0.1 by 0.1 do
  print(x/0.1);
  l:=[seq(evalf((i-1)*c+f[i](x)),i=1..n), evalf(f[1](x))] ; 
  SequenceOfPictures:=SequenceOfPictures, myconfig(l)
od:
 prespos:=[seq(evalf(f[i](x)),i=1..n)] ; 
 presvel:=[seq(evalf(g[i](x)),i=1..n)];
 print(`prespos=`.prespos);
 print(`presvel=`.presvel); 
end;




etatmixte:=proc(n,p, v,framenumber)

# n : number of masses on the circle (e.g. 4)
# p : n-vector of initial positions
# v : n-vector of initial velocities
# framenumber : of animation (e.g. 40)

global SequenceOfPictures;
suiteetatmixte(n,p,v,framenumber);
towindow();
end;



movetatmixte:=proc(n,p,v,root,firstframe,lastframe,nberdigits)

# This procedure immediately writes each picture it produces
# to a jpeg file. This makes it possible to write `movies`
# with more frames than would be possible to display
# in Maple (which crashes for lack of memory 
# while attempting to load the frames)

# n : number of masses on the circle (e.g. 4)
# p : n-vector of initial positions (e.g. [0,0,0,0])
# v : n-vector of initial velocities (e.g. [0.5,0,0,0])
# root : first part of the filename (excluding the path) (e.g. `movie`)
# firstframe (e.g. 0)
# lastframe  (e.g. 100)
# nberdigits : number of digits (including leading zeroes)
#     in the jpg-file names (e.g. 3)

local C, id, A, inipos, inispeed, inicond, functlist,
eq, sol, i, r, c, x, f, g, j, l, l1, frame,nber,filename;
global prespos, presvel, FullPath;

C:=companion(X^n-1,X):
id:=band([1], n):
A:=evalm(-1/2*(C+inverse(C))+id);

inipos:=vector(n,p);
inispeed:=vector(n,v);
 
inicond:=seq(y[i](0)=inipos[i],i=1..n), seq(D(y[i])(0)=inispeed[i],i=1..n):
functlist:=seq(y[i](t),i=1..n);
eq:=seq(diff(y[i](t),t$2)+add(A[i,j]*y[j](t),j=1..n)=0,i=1..n):

print(inicond);
sol:=dsolve({eq,inicond},{functlist},type=numeric, output=listprocedure);

for i from 1 to n do f[i]:=subs(sol, y[i](t)); od:
for i from 1 to n do g[i]:=subs(sol, diff(y[i](t),t)); od:
c:=2*Pi/n: r:=1:
x:=0;
frame:=firstframe;
while (frame<=lastframe) do
  nber:=addzeroes(frame,nberdigits);
  filename:=FullPath.root.nber.`.jpg`;
  print(filename);
  plotsetup(jpeg,plotoutput=filename,plotoptions=DimensionString);
  l:=[seq(evalf((i-1)*c+f[i](x)),i=1..n), evalf(f[1](x))] ; 
  print(display(myconfig(l),axes=none, scaling=constrained, view=[-1.2..1.2,-1.2..1.2]));
  x:=x+0.1;
  frame:=frame+1;
od:
 prespos:=[seq(evalf(f[i](x)),i=1..n)] ; 
 presvel:=[seq(evalf(g[i](x)),i=1..n)];
 print(`prespos=`,prespos);
 print(`presvel=`,presvel); 
end;