end

for i=(N/2+1):N;

    DJ1DF1(i)=(RJ1(i+1)-RJ1(i))./(XTA(i+1)-XTA(i))./2.0+(RJ1(i)-RJ1(i-1))./(XTA(i)-XTA(i-1))./2.0;

end

a=RJ1;

b=DJ1DF1;

c=F1+FM1;

dy(1)=y(2);

dy(2)=-b*(y(2))^2/2/a+c/a;

全满弹时弹仓的运动微分方程：

%弹筒的运动微分方程

function dy=dantong(t,y)

dy=zeros(2,1);

%求各时刻，各转角下的等效转动惯量及其对转角的导数

%求各时刻，个转角下的等效力矩（所有弹筒的重力产生的等效力矩，所有弹筒摩擦力产生的等效力矩）

%RJ-不同时刻各弹筒的等效转动惯量

%EJ-不同时刻所有弹筒的等效转动惯量

%F-不同时刻各（每个）弹筒重力产生的等效力矩

%EF-不同时刻所有弹筒的重力产生的等效力矩

%FM-不同时刻各（每个）弹筒摩擦力产生的等效力矩

%EFMDT-不同时刻所有弹筒摩擦力产生的等效力矩

OMGA=100.0/180.0*pi;

DT=0.5e-3;

N=pi/2.0/OMGA/DT;

for IT=1:(N+1)

    EJ(IT)=0.0;

    EF(IT)=0.0;

    EFMDT(IT)=0.0;

    T=(IT-1)*DT;

    XTA(IT)=T*OMGA;

    RM1=68.8398; %RM1--弹筒1的质量，由计算机测量获得

RJ1=0.294227; %RJ1--弹筒1的转动惯量，由计算机测量获得

U=0.15; %U--摩擦系数

OB=131e-3; %参数（1）

BA=sqrt(2*131e-3^2); %参数（2）

OB1=131e-3; %参数（3）

CA=sqrt(72e-3^2+92.63e-3^2); %参数（4）

CAB=atan(72e-3/92.63e-3); %参数（5）

R=pi/180.0;

BB1=sqrt(OB^2+OB^2-2.0*OB^2*cos(XTA)); %(1-1)

BB1O=(pi-XTA)/2.0; %(1-2)

BB1A=pi/2+BB1O; %(1-3)

B1AB=asin(BB1./BA.*sin(BB1A)); %(1-4)

ABB1=pi-BB1A-B1AB; %(1-5)

B1A=sqrt(BB1.^2+BA.^2-2.*BB1.*BA.*cos(ABB1)); %(1-6)

OA=sqrt(B1A.^2+OB1.^2); %(1-7)

S1=BA-B1A; %(1-8)

OBA=acos((BA.^2+OB.^2-OA.^2)./(2.0.*BA.*OB)); %(1-9)

ALF=OBA; %(1-10)

VB=OB.*OMGA; %(1-11)

VZ=VB.*sin(XTA)./sin(pi-XTA-ALF); %(1-12)

OMGAZ=VZ./BA; %(1-13)

VP1=VB.*sin(ALF)./sin(pi-XTA-ALF); %(1-14)

CAB1=CAB-B1AB; %(1-15)

BTA=CAB1; %(1-16)

GMA=pi-(pi./2.0-BTA); %(1-17)

VC=sqrt(VP1.^2+(CA.*OMGAZ).^2-2.0.*VP1.*CA.*OMGAZ.*cos(GMA)); %(1-18)

ALFG=acos(CA.*OMGAZ.*cos(BTA)./VC); %(1-19)

RJ1=RM1.*(VC./OMGA).^2+RJ1.*(OMGAZ./OMGA).^2;

F1=RM1.*9.8.*VC.*cos(ALFG)./OMGA;

FM1=RM1.*9.8/2.0.*U.*VP1/OMGA.*cos(pi);

RM2=68.8398; %RM2--弹筒2的质量，由计算机测量获得

RJ2=0.294227;  %RJ2--弹筒2的转动惯量，由计算机测量获得

U=1.5;  %U--摩擦系数

OC=92.63e-3+72.0e-3;

C1OB1=pi/4.0;

OMGAZ=OMGA;

VC=OC.*OMGA;

ALFG=abs(C1OB1-XTA);

RJ2=RM2.*(VC./OMGA).^2+RJ2.*(OMGAZ./OMGA).^2;

F2=RM2.*9.8.*VC.*cos(ALFG)./OMGA;

FM2=0;

%计算主动组合链轮匀速转动时弹筒3的运动规律：认为弹筒3做平面运动，求弹筒3的平均速度、转动速度及其质心的速度。

RM3=68.8398;

RJ3=0.294227;

U=0.15;

AO=131e-3;

B1O=131e-3;

AB=sqrt(2.0*131.0e-3^2);

ABC=atan(72.0e-3/92.63e-3);

CB=sqrt(72.0D-3^2+92.63D-3^2);

AB1=sqrt(AO.^2+B1O.^2-2.*AO.*B1O.*cos(pi/2-XTA));

B1AO=(pi-(pi/2-XTA))/2;

OB1A=(pi-(pi/2-XTA))/2;

BB1A=pi/2+OB1A;

ABB1=asin(AB1./AB.*sin(BB1A));

BAB1=pi-BB1A-ABB1;

B1B=sqrt(AB.^2+AB1.^2-2.0.*AB.*AB1.*cos(BAB1));

S3=B1B;

OB=sqrt(B1O.^2+B1B.^2);

BAO=acos((AO.^2+AB.^2-OB.^2)./(2.0.*AO.*AB));