global coor0, global vec0, global B, global vol0,
## Posicionamento dos pontos
vol0=600; B=1.36e-03;
n=30;
R0=(3*vol0/4/pi)^(1/3);
D0=sqrt(R0*R0-1); 
C0=[0,-D0]; 
th0=asin(1/R0);
dth=(pi-th0)/n;

coor0=[1 0];

for i=1:n
 the=th0+i*dth;
 coor0=[coor0;[R0*sin(the),-D0+R0*cos(the)]];
endfor

## Construção de vetores permitidos
AA=[0 -D0/2];
for i=1:n+1
 vec0(i,:)=(coor0(i,:)-AA)/norm(coor0(i,:)-AA);
endfor

y0=zeros(n+1,1);##vol0=volume(coor0);

##dy=0.6; ny=17;
##for iy2=1:3*ny+1
## for iy3=1:3*ny+1
##   yy2=dy*(iy2-ny-1);y2(iy2,iy3)=yy2;
##   yy3=dy*(iy3-ny-1);y3(iy2,iy3)=yy3;
##   ene(iy2,iy3)=fenergy([0;yy2;yy3]);
##   vv=fvolume([0;yy2;yy3]);vol(iy2,iy3)=vv(1)+vol0;
## endfor
##endfor

##for i=1:3*ny+1
## yx2=dy*(i-ny-1);x2(i)=yx2;
## ex2(i)=fenergy([0;yx2;0]);
## vvx=fvolume([0;yx2;0]);vx2(i)=vvx(1)+vol0;
##endfor

dvol=150;coru=[];yoptu=[];
for kk=1:12
 volu(kk)=400+kk*dvol;
 vol0=volu(kk);
 [yopt eopt info iter]=sqp(y0,@fenergy,@frestr,@frestr2,-5,15,550);
 yoptu=[yoptu yopt];
 y0=yopt;
 coor=coor0+diag(yopt)*vec0;
 coru=[coru coor];
endfor
hpp=plot(coru(:,1),coru(:,2),...
coru(:,3),coru(:,4),...
coru(:,5),coru(:,6),...
coru(:,7),coru(:,8),...
coru(:,9),coru(:,10),...
coru(:,11),coru(:,12),...
coru(:,13),coru(:,14),...
coru(:,15),coru(:,16),...
coru(:,17),coru(:,18),...
coru(:,19),coru(:,20),...
coru(:,21),coru(:,22),...
coru(:,23),coru(:,24))
set(hpp,"linewidth",2)
set(hpp,"marker","o")