clear;

N = 11;

alfa1 = 3*%pi/180;
lambda = 3*%pi/180;

b = 1.85;
c = 0.95;
d = c/4.;

kgm2 = 37.3/1000.;

r2 = b/N/2.;
t = 1./1000.;
r1 = r2-t;

if modulo(N,2)==1 then
    n = (N-1)/2 + 1;
else
    n = N/2;
end

zi = zeros(1,n);
ti = zeros(1,n);
si = zeros(1,n);
yi = zeros(1,n);

for i = 2:length(zi);
    ti(i) = ti(i-1) + 2*r2*cos(lambda)*sin(alfa1+(i-2)*lambda);
    si(i) = si(i-1) + 2*r2*cos(lambda)*cos(alfa1+(i-2)*lambda);
    zi(i) = ti(i) + r2*sin(alfa1+(i-3));
    yi(i) = si(i) + r2*cos(alfa1+(i-3));
end

zcm = 2*sum(zi)/N;
rho = 93.24;

m = rho*%pi*(r2^2-r1^2)*c;
Jxi = 0.5*m*(r2^2+r1^2);
Jyi = ((r2^2 + r1^2)/4. + c^2/12)*m;
Jzi = Jyi;

Jx = Jxi + m*((zcm-zi(1))^2+yi(1)^2);
Jz = Jyi + m*(zcm-zi(1))^2;
Jy = Jzi + m*yi(1)^2;

for j = 2:length(zi);
    Jx = Jx + 2*(Jxi + m*((zcm-zi(j))^2+yi(j)^2));
    Jy = Jy + 2*(Jyi + m*(zcm-zi(j))^2);
    Jz = Jz + 2*(Jzi + m*yi(j)^2);
end

J = [Jx 0 0; 0 Jy 0; 0 0 Jz];
disp(J);

M = N*m;