用魔术公式建立轮胎模型根据资料，通过给定的参数以及对一些比例因子的设定，用魔术公式建立如下轮胎模型。

1.纯纵滑情况：%% Longitudinal Force (pure longitudinal slip)%% k is input value, means slip ratio%% Fx0 is output value, means longitudinal forcek=linspace(-1,1,300);dfz=(Fz-Fz0)/Fz0;C_x=P_Cx1; %取lam(Cx)=1 shape factoru_x=(P_Dx1+P_Dx2*dfz)*0.9; %取lam(ux*)=0.9 composite friction scaling fatorD_x=u_x*Fz; %zeta(1-8)都取1S_Hx=(P_Hx1+P_Hx2*dfz); %取lam(Hx)=1 horizontal shiftk_x=k+S_Hx;E_x=(P_Ex1+P_Ex2*dfz+P_Ex3*dfz^2)*(1-P_Ex4*sign(k_x)); %取lam(Ex)=1 curvature factorK_xk=Fz*(P_Kx1+P_Kx2*dfz)*exp(P_Kx3*dfz); %取lam(Kxk)=1 break slip stiffnessB_x=K_xk/(C_x*D_x+0.1); %取epsilon=0.1S_vx=Fz*(P_Vx1+P_Vx2*dfz)*(V0/(0.1+V0)); %lam(Vx),lam(ux),zeta1=1,epsil on=0.1Fx0=D_x*sin(C_x*atan(B_x*k_x-E_x.*(B_x*k_x-atan(B_x*k_x))))+S_vx;figure(1);plot(k,Fx0,'linewidth',1.5);gridxlabel('纵向滑移率');ylabel('纵向力/（N）');title('纵向力--滑移率(纯纵滑)');2.纯侧滑情况：%% Lateral Force and Torque (pure side slip)%% a is input value, means slip angle%% Fy0 is output value, means side force%% Mz0 is output value, means aligning torquea=linspace(-15,15,300); %侧偏角取-15-15 度ta=tan(a*pi/180);sr=sin(r); %r是外倾角，sr 表示r*C_y=P_Cy1; %取lam(Cy)=1u_y=(P_Dy1+P_Dy2*dfz)/(1+P_Dy3*sr^2); %取lam(uy*)=1D_y=u_y*Fz; %zeta(1-8)都取1K_ya=P_Ky1*Fz0*sin(P_Ky4*atan(Fz/((P_Ky2+P_Ky5*sr^2)*Fz0)))/(1+P_Ky3*sr^2); %取lam（Kya）=1 zeta(3)=1B_y=K_ya/(C_y*D_y); %取epsilony=0；K_yr0=Fz*(P_Ky6+P_Ky7*dfz); %取lam（Kyr)=1S_Vyr=Fz*(P_Vy3+P_Vy4*dfz)*sr;%lam(kyr)=1, camber force stiffness 后面的修正系数全部取1S_Hy=(P_Hy1+P_Hy2*dfz)+(K_yr0*sr-S_Vyr)/(K_ya+0.1); %zeta(0,4)都取1 epsilon(k)=0.1 S_Vy=Fz*(P_Vy1+P_Vy2*dfz)+S_Vyr; %lam(Vy,uy)=1,zeta(2)=1ay=ta+S_Hy; %tan(sa)表示横向侧偏角a*E_y=(P_Ey1+P_Ey2*dfz)*(1+P_Ey5*sr^2-(P_Ey3+P_Ey4*sr)*sign(ay)); %lam(ey)=1Fy0=D_y*sin(C_y*atan(B_y*ay-E_y.*(B_y*ay-atan(B_y*ay))))+S_Vy;figure(2);plot(a,Fy0,'linewidth',1.5);gridxlabel('侧偏角');ylabel('侧向力/（N）');title('侧向力-侧偏角(纯侧滑)');C_r=1;C_t=q_Cz1;B_t=(q_Bz1+q_Bz2*dfz+q_Bz3*dfz^2)*(1+q_Bz5*abs(sr)+q_Bz6*sr^2);S_Ht=q_Hz1+q_Hz2*dfz+(q_Hz3+q_Hz4*dfz)*sr;a_t=ta+S_Ht;B_r=(q_Bz9+q_Bz10*B_y*C_y);E_t=(q_Ez1+q_Ez2*dfz+q_Ez3*dfz^2)*(1+(q_Ez4+q_Ez5*sr)*(2/pi)*atan(B_t*C_t*a_t));S_Hf=S_Hy+S_Vy/K_ya;ar=ta+S_Hf;D_t0=Fz*(R0/Fz0)*(q_Dz1+q_Dz2*dfz); %取lambda(t)=1D_t=D_t0*(1+q_Dz3*abs(sr)+q_Dz4*sr^2);t0=D_t*cos(C_t*atan(B_t*a_t-E_t.*(B_t*a_t-atan(B_t*a_t))));D_r=Fz*R0*(q_Dz6+q_Dz7*dfz)+(q_Dz8+q_Dz9*dfz)*sr+(q_Dz10+q_Dz11*dfz)*sr*abs(sr);M_zr0=D_r*cos(C_r*atan(B_r*ar)); M_z0=M_zr0-t0.*Fy0;figure(3);plot(a,M_z0,'linewidth',1.5);grid on;xlabel('侧偏角');ylabel('回正力矩/（N）');title('回正力矩-侧偏角(纯侧滑)')3.组合滑移情况：%% Longitudinal Force (combined slip)%% aS is input value, means slip angle%% Fx is output value, means longitudinal forceS_Hxa=r_Hx1;E_xa=r_Ex1+r_Ex2*dfz;C_xa=r_Cx1;B_xa=(r_Bx1+r_Bx3*sr^2)*cos(atan(r_Bx2*k)); %lam(xa)=1,influence on Fx 矩阵aS=ta+S_Hxa; %矩阵G_xa0=cos(C_xa*atan(B_xa*S_Hxa-E_xa*(B_xa*S_Hxa-atan(B_xa*S_Hxa))));G_xa=cos(C_xa*atan(B_xa.*aS-E_xa*(B_xa.*aS-atan(B_xa.*aS))))./G_xa0;Fx=G_xa.*Fx0;figure(4);plot(a,Fx,'linewidth',1.5);grid on;xlabel('侧偏角');ylabel('纵向力/（N）');title('纵向力-侧偏角(混合滑移)');%% Lateral Force (combined slip)%% kS is input value, means slip ratio%% Fy is output value, means longitudinal forceC_yk=r_Cy1;E_yk=r_Ey1+r_Ey2*dfz;B_yk=(r_By1+r_By4*sr^2)*cos(atan(r_By2*(k-r_By3))); %lam(yk)=1D_Vyk=u_y*Fz*(r_Vy1+r_Vy2*dfz+r_Vy3*sr)*cos(atan(r_Vy4*ta)); %zeta(2)=1S_Vyk=D_Vyk.*sin(r_Vy5*atan(r_Vy6*k)); %取lam（Vyk）=1 S_Hyk=r_Hy1+r_Hy2*dfz;kS=k+S_Hyk;G_yk0=cos(C_yk*atan(B_yk*S_Hyk-E_yk*(B_yk*S_Hyk-atan(B_yk*S_Hyk))));G_yk=cos(C_yk*atan(B_yk.*kS-E_yk*(B_yk.*kS-atan(B_yk.*kS))))/G_yk0;Fy=G_yk.*Fy0+S_Vyk;figure(5);plot(k,Fy,'linewidth',1.5);grid on;xlabel('纵向滑移率');ylabel('侧向力/（N）');title('侧向力-纵向滑移率(混合滑移)');%% Aligning Torque (combined slip)%% ateq is input value, means slip angle%% Mz is output value, means Aligning Torqueateq=(sqrt(a_t.^2+(K_xk/K_ya)^2*(k.^2))).*sign(a_t);areq=(sqrt(ar.^2+(K_xk/K_ya)^2*(k.^2))).*sign(ar);M_zr=D_r*(C_r*atan(B_r*areq));s=R0*(S_sz1+S_sz2*(Fy/Fz0)+(S_sz3+S_sz4*dfz)*sr); %取lambda(s)=1 F_yy=Fy-S_Vyk;t=D_t*cos(C_t*atan(B_t*ateq-E_t.*(B_t*ateq-atan(B_t*ateq))));M_zz=-t.*F_yy;Mz=M_zz+M_zr+s.*Fx;figure(6);plot(a,Mz,'linewidth',1.5);grid on;xlabel('侧偏角');ylabel('回正力矩/（N）');title('回正力矩-侧偏角(混合滑移)');。