********************************************************

%1.輸入已知數(shù)據(jù)

clear；

l1=0.400；

l2=1.200；

las1=0.2；

lbs2=0.6；

m1=1.2；

m2=3.6；

m3=6；

g=10；

J2=0.45；

G1=m1^*g；

G2=m2^*g；

G3=m3^*g；

Fr=-1000；

e=0；

hd=pi/180；

du=180/pi；

omega1=10；

alpha1=0；

%2.曲柄滑塊機(jī)構(gòu)力平衡計(jì)算

for n1=1：360

　theta1（n1）=（n1-1）^*hd；

　% 調(diào)用函數(shù) slider_crank 計(jì)算曲柄滑塊機(jī)構(gòu)位移，速度，加速度

［theta2（n1），s3（n1），omega2（n1），v3（n1），alpha2（n1），a3（n1）］=slider_crank（theta1（n1），omega1，alpha1，l1，l2，e）；

　% 計(jì)算各個(gè)質(zhì)心點(diǎn)加速度

　as1x（n1）=-las1^*cos（theta1（n1））^*omega1^∧2；

　as1y（n1）=-las1^*sin（n1^*hd）^*omega1^∧2；

　as2x（n1）=-l1^*omega1^∧2^*cos（n1^*hd）-lbs2^*（omega2（n1）^∧2^*cos（theta2（n1））﹢alpha2（n1）^*sin（theta2（n1）））；　%質(zhì)心s2在x軸的加速度

　as2y（n1）=-l1^*omega1^∧2^*sin（n1^*hd）-lbs2^*（omega2（n1）^∧2^*sin（theta2（n1））-alpha2（n1）^*cos（theta2（n1）））；　%質(zhì)心s2在y軸的加速度

　% 計(jì)算各構(gòu)件慣性力和慣性力矩

　F1x（n1）=-as1x（n1）^*m1；

　F1y（n1）=-as1y（n1）^*m1；

　F2x（n1）=-as2x（n1）^*m2；

　F2y（n1）=-as2y（n1）^*m2；

　F3x（n1）=-a3（n1）^*m3；

　F3y（n1）=0；

　FR43x（n1）=Fr；

　M2（n1）=-alpha2（n1）^*J2；

　% 計(jì)算各個(gè)鉸鏈點(diǎn)坐標(biāo)， 計(jì)算各個(gè)質(zhì)心點(diǎn)坐標(biāo)

　xa=0；

　ya=0；

　xs1=las1^*cos（n1^*hd）；

　ys1=las1^*sin（n1^*hd）；

　xb=l1^*cos（n1^*hd）；

　yb=l1^*sin（n1^*hd）；

　xs2=xb﹢lbs2^*cos（theta2（n1））；

　ys2=yb﹢lbs2^*sin（theta2（n1））；

　xc=xb﹢l2^*cos（theta2（n1））；

　yc=yb﹢l2^*sin（theta2（n1））；

　　% 未知力系數(shù)矩陣

　　A=zeros（8）；

　　A（1，1）=1；A（1，2）=-（ys1-yb）；A（1，3）=-（xb-xs1）；A（1，4）=-（ys1-ya）；

　　A（1，5）=-（xa-xs1）；

　　A（2，2）=-1；A（2，4）=-1；

　　A（3，3）=-1；A（3，5）=-1；

　　A（4，2）=（ys2-yb）；A（4，3）=（xb-xs2）；A（4，6）=-（ys2-yc）；A（4，7）=-（xc-xs2）；

　　A（5，2）=1；A（5，6）=-1；

　　A（6，3）=1；A（6，7）=-1；

　　A（7，6）=1；

　　A（8，7）=1；A（8，8）=-1；

　　%已知力列陣

　　B=zeros（8，1）；

　　B（2）=-F1x（n1）；

　　B（3）=-F1y（n1）﹢G1；

　　B（4）=-M2（n1）；

　　B（5）=-F2x（n1）；

　　B（6）=-F2y（n1）﹢G2；

　　B（7）=-F3x（n1）﹢FR43x（n1）；

　　B（8）=-F3y（n1）；

　　C=A\B；

　　Mb（n1）=C（1）；Fr12x（n1）=C（2）；Fr12y（n1）=C（3）；Fr14x（n1）=C（4）；Fr14y（n1）=C（5）；

　　Fr23x（n1）=C（6）；Fr23y（n1）=C（7）；Fr34y（n1）=C（8）；

end

　　%3.曲柄滑塊機(jī)構(gòu)力分析圖形輸出

　figure（2）；

　n1=1∶360；

　subplot（2，2，1）；　%繪運(yùn)動(dòng)副反力FR14曲線圖

　plot（n1， Fr14x，'b'）；

　hold on

　plot（n1，F(xiàn)r14y，'k--'）；

　legend（'F_R_1_4_x'，'F_R_1_4_y'）

　title（'運(yùn)動(dòng)副反力F_R_1_4曲線圖'）；

　xlabel（'曲柄轉(zhuǎn)角 \theta_1/\circ'）

　ylabel（'F/N'）

　grid on；

　subplot（2，2，2）；　%繪運(yùn)動(dòng)副反力FR23曲線圖

　plot（n1，F(xiàn)r23x（n1），'b'）；

　hold on

　plot（n1，F(xiàn)r23y（n1），'k--'）；

　hold on

　legend（'F_R_2_3_x'，'F_R_2_3_y'）

　title（'運(yùn)動(dòng)副反力F_R_2_3曲線圖'）；

　xlabel（'曲柄轉(zhuǎn)角 \theta_1/\circ'）

　ylabel（'F/N'）

　grid on；

　subplot（2，2，3）；　%繪運(yùn)動(dòng)副反力FR34曲線圖

　plot（n1，F(xiàn)r34y，'b'）；

　hold on

　legend（'F_R_3_4_y'）

　title（'運(yùn)動(dòng)副反力F_R_3_4_y曲線圖'）；

　xlabel（'曲柄轉(zhuǎn)角 \theta_1/\circ'）

　ylabel（'F/N'）

　grid on；

　subplot（2，2，4）；　　%繪平衡力矩M_b曲線圖

　plot（n1，Mb）

　title（'力矩Mb圖'）

　xlabel（'曲柄轉(zhuǎn)角 \theta_1/\circ'）；

　ylabel（'M/N.m'）

　hold on；

　grid on；