n = 600:1:4000;

Ttq = -19.313+295.27*n/1000-165.44*(n/1000).^2+40.87*(n/1000).^3-3.8445*(n/1000).^4;

%io=5.83;

hold on;

for io=5.80:.03:6.80

t=0.844;

r=0.367;

ig2=2.769;

ig3=1.644;

ig4=1.00;

ua2=0.377*r*n/(ig2*io);

ua3=0.377*r*n/(ig3*io);

ua4=0.377*r*n/(ig4*io);

Ft2=Ttq*ig2*io*t/r;

Ft3=Ttq*ig3*io*t/r;

Ft4=Ttq*ig4*io*t/r;

m=3880;

f=0.012;

g=9.8;

CdA=2.77;

Ff2=m*g*f;

Fw2=CdA*ua2.^2/21.15;

Fr2=Ff2+Fw2;

F2=Ft2-Fr2;

a2=m./F2;

Ff3=m*g*f;

Fw3=CdA*ua3.^2/21.15;Fr3=Ff3+Fw3;

F3=Ft3-Fr3;a3=m./F3;

Ff4=m*g*f;

Fw4=CdA*ua4.^2/21.15;

Fr4=Ff4+Fw4;

F4=Ft4-Fr4;

a4=m./F4;

m2=5.142415:0.001:34.282777;

m3=34.28277:0.001:57.74269;

m4=57.7427:0.001:70;

t2=interp1(ua2,a2,m2,"spline");

t3=interp1(ua3,a3,m3,"spline");

t4=interp1(ua4,a4,m4,"spline");

plot(m2,t2,m3,t3,m4,t4);

%hold on; %面積

A=0.001*sum(t2-0)+0.001*sum(t3-0)+0.001*sum(t4-0);

disp(["io = " num2str(io) " A = " num2str(A)]);

%T=86.8332/3.6

end;