clear *
makeplots=0;
R=1; % radius of circle
nR = 40 ; % number of divisions to create circle
area=0.004;
w=linspace(0,2*pi*(1-1/nR),nR);  % angles of points on circle
xy=[R*cos(w);R*sin(w);ones(1,nR)]';

CreateMeshTriangleQ("circle",xy,area)
[nodes,elem,edges]=ReadMeshTriangle("./circle.1");
% ShowMesh(nodes,elem)
% disp("hit return");pause();

[la,vec]=FEMEig(nodes,elem,edges,1,0,1,4,1e-6);
la=la'

gset nokey
figure(1);
ShowSolution(nodes,elem,vec(:,1))
if makeplots print("eigen1.eps","-depsc")endif
figure(2);
vmin=min(vec(:,1)); vmax=max(vec(:,1));
ShowLevelCurves(nodes,elem,vec(:,1),linspace(vmin,vmax,11))
if makeplots print("level1.eps","-depsc")endif

disp("Hit RETURN"); pause();
figure(1);
ShowSolution(nodes,elem,vec(:,2))
if makeplots print("eigen2.eps","-depsc")endif
figure(2);
vmin=min(vec(:,2)); vmax=max(vec(:,2));
ShowLevelCurves(nodes,elem,vec(:,2),linspace(vmin,vmax,11))
if makeplots print("level2.eps","-depsc")endif

disp("Hit RETURN"); pause();
figure(1);
ShowSolution(nodes,elem,vec(:,4))
if makeplots print("eigen4.eps","-depsc")endif
figure(2);
vmin=min(vec(:,4)); vmax=max(vec(:,4));
ShowLevelCurves(nodes,elem,vec(:,4),linspace(vmin,vmax,11))
if makeplots print("level4.eps","-depsc")endif
