function main(m,k)
global G yitaT r f CdA If Iw1  g
n=600:10:4000;
Tq=-19.313+295.27*(n/1000)-165.44*(n/1000).^2+40.874*(n/1000).^3-3.8445*(n/1000).^4;
g=9.8;G=m*g;yitaT=0.85;r=0.367;f=0.013;
CdA=2.77;If=0.218;Iw1=1.798;global Iw2Iw2=3.598;global L
L=3.2;global a
a=1.947;global hg
hg=0.9;
global  ig
ig=[6.09 3.09 1.71 1.00];
global Ff
Ff=G*f;
global Ft1 Ft2 Ft3 Ft4 ua1 ua2 ua3 ua4
Ft1=Tq*ig(1)*k*yitaT/r;
Ft2=Tq*ig(2)*k*yitaT/r;
Ft3=Tq*ig(3)*k*yitaT/r;
Ft4=Tq*ig(4)*k*yitaT/r;
ua1=0.377*r*n/ig(1)/k;
ua2=0.377*r*n/ig(2)/k;
ua3=0.377*r*n/ig(3)/k;
ua4=0.377*r*n/ig(4)/k;
subplot(2,2,1)
force()
subplot(2,2,2)
acc(m,k)
subplot(2,2,3)
ret(m,k)
subplot(2,2,4)
time(m,k)
%*驱动力-行驶阻力平衡图
function force()
ua=0:5:120;global Ff CdA Ft1 Ft2 Ft3 Ft4 ua1 ua2 ua3 ua4
Fw=CdA*ua.^2/21.15;
Fz=Ff+Fw;
plot(ua1,Ft1,ua2,Ft2,ua3,Ft3,ua4,Ft4,ua,Fz)
title('驱动力-行驶阻力平衡图')
xlabel('ua-km/h')
ylabel('F-N')
gtext('Ft1'),gtext('Ft2'),gtext('Ft3'),gtext('Ft4'),gtext('Ff+Fw')
%一档最大爬坡度%
%Ftt=Ft1-Fz;
%imax=tan(asin(max(Ftt/G)));
%imax;
%行驶加速度曲线，加速度倒数曲线图%
function acc(m,k)
global yitaT r CdA If Iw1 Iw2  ig Ft1 Ft2 Ft3 Ft4 ua1 ua2 ua3 ua4  Ff
deta1=1+(Iw1+Iw2)/(m*r^2)+(If*(ig(1))^2*k^2*yitaT)/(m*r^2);
deta2=1+(Iw1+Iw2)/(m*r^2)+(If*(ig(2))^2*k^2*yitaT)/(m*r^2);
deta3=1+(Iw1+Iw2)/(m*r^2)+(If*(ig(3))^2*k^2*yitaT)/(m*r^2);
deta4=1+(Iw1+Iw2)/(m*r^2)+(If*(ig(4))^2*k^2*yitaT)/(m*r^2);
Fw1=CdA*ua1.^2/21.15;
Fw2=CdA*ua2.^2/21.15;
Fw3=CdA*ua3.^2/21.15;
Fw4=CdA*ua4.^2/21.15;
a1=(Ft1-Ff-Fw1)/(deta1*m);
a2=(Ft2-Ff-Fw2)/(deta2*m);
a3=(Ft3-Ff-Fw3)/(deta3*m);
a4=(Ft4-Ff-Fw4)/(deta4*m);
plot(ua1,a1,ua2,a2,ua3,a3,ua4,a4)
axis([0 99 0 2.5])
title('汽车的行驶加速度曲线')
xlabel('ua-km/h')
ylabel('a-m/s^2')
gtext('Ⅰ'),
gtext('Ⅱ'),
gtext('Ⅲ'),
gtext('Ⅳ'),
function ret(m,k)
global yitaT r CdA If Iw1 Iw2  ig Ft1 Ft2 Ft3 Ft4 ua1 ua2 ua3 ua4  Ff
deta1=1+(Iw1+Iw2)/(m*r^2)+(If*(ig(1))^2*k^2*yitaT)/(m*r^2);
deta2=1+(Iw1+Iw2)/(m*r^2)+(If*(ig(2))^2*k^2*yitaT)/(m*r^2);
deta3=1+(Iw1+Iw2)/(m*r^2)+(If*(ig(3))^2*k^2*yitaT)/(m*r^2);
deta4=1+(Iw1+Iw2)/(m*r^2)+(If*(ig(4))^2*k^2*yitaT)/(m*r^2);
Fw1=CdA*ua1.^2/21.15;
Fw2=CdA*ua2.^2/21.15;
Fw3=CdA*ua3.^2/21.15;
Fw4=CdA*ua4.^2/21.15;
a1=(Ft1-Ff-Fw1)/(deta1*m);inv_a1=1./a1;
a2=(Ft2-Ff-Fw2)/(deta2*m);inv_a2=1./a2;
a3=(Ft3-Ff-Fw3)/(deta3*m);inv_a3=1./a3;
a4=(Ft4-Ff-Fw4)/(deta4*m);inv_a4=1./a4;
plot(ua1,inv_a1,ua2,inv_a2,ua3,inv_a3,ua4,inv_a4)
axis([0 99 0 10])
title('汽车的行驶加速度倒数曲线')
xlabel('ua-km/h')
gtext('Ⅰ'),
gtext('Ⅱ'),
gtext('Ⅲ'),
gtext('Ⅳ'),
%加速时间曲线图%
function time(m,k)
global  g G yitaT r f CdA i0 If Iw1 Iw2 L a hg ig
nmin=600;nmax=4000;
u1=0.377*r*nmin./ig/k;
u2=0.377*r*nmax./ig/k;
deta=0*ig;
for i=1:4
    deta(i)=1+(Iw1+Iw2)/(m*r^2)+(If*(ig(i))^2*k^2*yitaT)/(m*r^2);
end
deta_u=0.5;
ua=3:deta_u:80;
N=length(ua);n=0;Tq=0;Ft=0;inv_a=0*ua;delta=0*ua;
Ff=G*f;
Fw=CdA*ua.^2/21.15;
for i= 1:N
    j=i;
   if ua(i)<=u2(1)
       n=ua(i)*(ig(1)*k/r)/0.377;
       Tq=-19.313+295.27*(n/1000)-165.44*(n/1000).^2+40.874*(n/1000).^3-3.8445*(n/1000).^4;
        Ft=Tq*ig(1)*k*yitaT/r;
        inv_a(i)=(deta(1)*m)/(Ft-Ff-Fw(i));
        delta(i)=deta_u*inv_a(i)/3.6;
   elseif ua(i)<=u2(2)
        n=ua(i)*(ig(2)*k/r)/0.377;
        Tq=-19.313+295.27*(n/1000)-165.44*(n/1000).^2+40.874*(n/1000).^3-3.8445*(n/1000).^4;
        Ft=Tq*ig(2)*k*yitaT/r;
        inv_a(i)=(deta(2)*m)/(Ft-Ff-Fw(i));
        delta(i)=deta_u*inv_a(i)/3.6;
        elseif ua(i)<=u2(3)
            n=ua(i)*(ig(3)*k/r)/0.377;
        Tq=-19.313+295.27*(n/1000)-165.44*(n/1000).^2+40.874*(n/1000).^3-3.8445*(n/1000).^4;
        Ft=Tq*ig(3)*k*yitaT/r;
        inv_a(i)=(deta(3)*m)/(Ft-Ff-Fw(i));
        delta(i)=deta_u*inv_a(i)/3.6;
        else 
             n=ua(i)*(ig(4)*k/r)/0.377;
        Tq=-19.313+295.27*(n/1000)-165.44*(n/1000).^2+40.874*(n/1000).^3-3.8445*(n/1000).^4;
        Ft=Tq*ig(4)*k*yitaT/r;
        inv_a(i)=(deta(4)*m)/(Ft-Ff-Fw(i));
        delta(i)=deta_u*inv_a(i)/3.6;
    end

        a=delta(1:j);
        t(i)=sum(a);
end
    T=t(N)
    plot(t,ua);grid;
    axis([0 60 0 100])
    title('汽车加速时间曲线')
    xlabel('t-s')
    ylabel('ua-km/h')