2020-12-24

关于ODE45积分器的一些代码实现：

[t,theta]   = ode45(@(t,theta) odefun(t,theta,omega_pz,GearSet,slices_num,k_single,J_wz,r_wb,Tin,Transmission_ratio), t_span, [theta_wz_0,omga_wz_0]);%%注意这里的t是系统自动变化的
% out         = [t,theta];
% legend('wheel theta in rad','omega in rad per s','pinion theta in rad')

legend('wheel theta in rad','pinion theta in rad')
% 
% plot(1:length(t)/381:length(t),t*omega_pz,'-o')
figure(3)
TE=(Pinion_ang-Wheel_ang*Transmission_ratio)*db(1)/2*1e3;%[um]
plot(t,TE)
xlabel('Time in s')
ylabel('TE in um')

%%
function theta_dt = odefun(t,theta,omega_pz,GearSet,slices_num,k_single,J_wz,r_wb,Tin,Transmission_ratio)%%注意这里的t是系统自动变化的
    r_pz_all    = [];
    r_wz_all    = [];
    delt_L      = [];
    dir = Tin/abs(Tin)
    theta_input = dir*(omega_pz*t)*180/pi%这里输入是degree
    Points_data = Cal_Contact_Points_th_Fun(GearSet,slices_num,theta_input);
    r_pz_all    = Points_data(:,4);%小齿轮 接触点半径
    r_wz_all    = Points_data(:,5);%大齿轮接触点半径
    for i = length(r_pz_all)
        delt_L(i,1) = omega_pz*t*r_pz_all(i,1)-theta(1)*r_wz_all(i,1);%%注意这里的t是系统自动变化的
        if delt_L(i,1)>0
            delt_L(i,1)=delt_L(i,1);
        else
            delt_L(i,1)=0;
        end
    end
    theta_dt    = zeros(2,1);
    theta_dt(1) = theta(2);%θ的一阶导数是ω这里的theta(2)是角速度
    theta_dt(2) = (-dir*k_single*sum(delt_L)*r_wb/J_wz+Tin*Transmission_ratio/J_wz);%θ的二阶导数是α
    figure(1)
    hold on
    plot(t,(-dir*k_single*sum(delt_L)*r_wb/J_wz+Tin*Transmission_ratio/J_wz),'or','MarkerSize',10,'MarkerFaceColor','r')
    xlabel('time in s')
    ylabel('T_e_x in Nmm')
end