## Posicionamento dos pontos
vol0=2770; B=1.36e-03;
n=15;
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

hh=plot(coor0(:,1),coor0(:,2),"-o");
set(hh,"linewidth",2);
##hold on

AA=[0 -D0/2];
for i=1:n+1
 vec0(i,:)=(coor0(i,:)-AA)/norm(coor0(i,:)-AA);
endfor

y=0*[0 rand(1,15)-0.5]; yd=diag(y);

coor=coor0+yd*vec0;

hh=plot(coor0(:,1),coor0(:,2),":",coor(:,1),coor(:,2),"-x");
set(hh,"linewidth",2);
##break;

N1=3; xi1=[-sqrt(3/5),0,sqrt(3/5)];A1=[5/9,8/9,5/9];
N2=20; xi2=-1+1/20:2/20:1-1/20;A2=2/20*ones(1,N2);
##N2=50; xi2=-1+1/50:2/50:1-1/50;A2=2/50*ones(1,N2);
##N2=150; xi2=-1+1/150:2/150:1-1/150;A2=2/150*ones(1,N2);

pint1=[(1-xi1)/2;(1+xi1)/2];
pint2=[(1-xi2)/2;(1+xi2)/2];

d=[coor(2:n+1,:)-coor(1:n,:)];ll=norm(d,"rows");jac=ll/2;

normal=-[d(:,2)./ll,-d(:,1)./ll]

N=N1; xi=xi1; A=A1; pint=pint1;
perim=0;
for i=1:n
 perim=perim+A*ones(N,1)*jac(i);
endfor 

perim

res=0;
for i=1:n
 for k=1:N
  x=pint(1,k)*coor(i,:)'+pint(2,k)*coor(i+1,:)';
  res=res+A(k)*2*pi/3*normal(i,1:2)*x*x(1)*jac(i);
 endfor
endfor 

volume=res

res=0;
for i=1:n
 for k=1:N
  x=pint(1,k)*coor(i,:)'+pint(2,k)*coor(i+1,:)';
  res=res+A(k)*2*pi*x(1)*jac(i);
 endfor
endfor 

area=res

res=0;
for i=1:n
 for k=1:N
  x=pint(1,k)*coor(i,:)'+pint(2,k)*coor(i+1,:)';
  res=res+A(k)*(B*x(2)^2*x(1)*normal(i,2)+2*x(1))*jac(i);
 endfor
endfor 

energia=res