谁能帮我修改一下程序啊，非常感谢！！！！
程序如下：
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                输入已知数据
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
l1=40;
l2=50;
l3=75;
l4=35;
l6=70;
l5=60;
omega1=10;
omega6=0;
hd=pi/180;
du=180/pi;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               四杆机构曲线设计
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for n1=1:asin(6/7)*du
    theta4(n1)=n1+asin((l6/l5)*sin(n1*hd))*du;   %  杆DF转过的角度
    l7(n1)=l5*sin(theta4(n1)*hd)/sin(n1*hd);     %  杆AF的长度
    l8(n1)=sqrt((l6-l4*cos(theta4(n1)*hd)+l1*cos((n1)*hd))^2+(-l4*sin(theta4(n1)*hd)+l1*sin((n1)*hd))^2);  % BD的长度
    theta2(n1)=(n1*hd-acos(((l7(n1)+l1)^2+l8(n1)^2-(l4+l5)^2)/(2*(l7(n1)+l1)*l8(n1)))+acos((l2^2+l8(n1)^2-l3^2)/(2*l2*l8(n1))))*du;  %  杆BC转过的角度
    theta3(n1)=(theta4(n1)*hd+acos((l8(n1)^2+(l4+l5)^2-(l1+l7(n1))^2)/(2*l8(n1)*(l4+l5)))-acos((l8(n1)^2+l3^2-l2^2)/(2*l3*l8(n1))))*du;    %  杆CD转过的角度
   
    A=[-l2*sin(theta2(n1)*hd),l3*sin(theta3(n1)*hd),l4*sin(theta4(n1)*hd+pi),0;
        l2*cos(theta2(n1)*hd),-l3*cos(theta3(n1)*hd),-l4*cos(theta4(n1)*hd+pi),0;
        0,0,l5*sin(theta4(n1)*hd),cos(n1*hd);
        0,0,-l5*cos(theta4(n1)*hd),-sin(n1*hd)];
    B=[l1*sin(n1*hd+pi);-l1*cos(n1*hd+pi);l7(n1)*sin(n1*hd); -l7(n1)*cos(n1*hd)];
    P1=(A\B )*omega1 ;  
    omega2(n1)=P1(1);
    omega3(n1)=P1(2);
    omega4(n1)=P1(3);
    v7=P1(4);
    vcx(n1)=-l1*omega1*sin(n1*hd+pi)-l2*omega2(n1)*sin(theta2(n1)*hd);
    vcy(n1)=l1*omega1*cos(n1*hd+pi)+l2*omega2(n1)*cos(theta2(n1)*hd);
    vc(n1)=sqrt(vcx(n1)^2+vcy(n1)^2);  %  C点的速度  
   
    C=[-l2*sin(theta2(n1)*hd),l3*sin(theta3(n1)*hd),-l4*sin(theta4(n1)*hd),0;
        -l2*cos(theta2(n1)*hd),-l3*cos(theta3(n1)*hd),l4*cos(theta4(n1)*hd),0;
        0,0,l5*sin(theta4(n1)*hd),cos(n1*hd);
        0,0,-l5*cos(theta4(n1)*hd),sin(n1*hd)];
    D=[-l2*omega2(n1)*cos(theta2(n1)*hd),l3*omega3(n1)*cos(theta3(n1)*hd),-l4*omega4(n1)*cos(theta4(n1)*hd),0;
        -l2*omega2(n1)*sin(theta2(n1)*hd),l3*omega3(n1)*sin(theta3(n1)*hd),-l4*omega4(n1)*sin(theta4(n1)*hd),0;
        0,0,l5*omega4(n1)*cos(theta4(n1)*hd),-omega1*sin(n1*hd);
        0,0,l5*omega4(n1)*sin(theta4(n1)*hd),omega1*cos(n1*hd)];
    E=[l1*omega1*cos(n1*hd);
        l1*omega1*sin(n1*hd);
        l7*omega1*cos(n1*hd)+v7*sin(n1*hd);
        l7*omega1*sin(n1*hd)+v7*cos(n1*hd);];
    P2=C\(E*omega1-D*P1);
    alpha2(n1)=P2(1);
    alpha3(n1)=P2(2);
    alpha4(n1)=P2(3);

   
    acx(n1)=l1*omega1^2*cos(n1*hd)-l2*omega2(n1)^2*cos(theta2(n1)*hd)-l2*alpha2(n1)*sin(theta2(n1)*hd);
    acy(n1)=l1*omega1^2*sin(n1*hd)-l2*omega2(n1)^2*sin(theta2(n1)*hd)-l2*alpha2(n1)*cos(theta2(n1)*hd);
    ac(n1)=sqrt(acx(n1)^2+acy(n1)^2);  % C点的加速度
end
    disp 'vc(50)=',vc(50)
    disp 'ac(50)=',ac(50)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %                     位移，角速度，角加速度和四杆机构图形输出
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
n1=1:asin(6/7)*du;

    subplot(2,2,1);   % 绘制C点的速度曲线
    plot(n1,vc)
    title('c点速度线图');
    xlabel('l7转角 \theta_1 / \circ')
    ylabel('C点速度 / mm\cdots^{-1}')
    grid on;hold on;
    text(100,-300,'v_c')

    subplot(2,2,2);  % 绘制C点的加速度曲线  
    plot(n1,ac)
    title('c点加速度线图');
    xlabel('l7转角 \theta_1 / \circ')
    ylabel('C点加速度 / mm\cdots^{-2}')
    grid on;hold on;
    text(100,-300,'a_c')
   
    subplot(2,2,3); % 设计铰链四杆机构
    x(1)=0;
    y(1)=0;
    x(2)=l7(50)*cos(50*hd);
    y(2)=l7(50)*sin(50*hd);
    x(3)=(l5+l6+10)*cos(50*hd);
    y(3)=(l5+l6+10)*sin(50*hd);
    x(4)=l7(50)*cos(50*hd);
    y(4)=l7(50)*sin(50*hd);
    x(5)=l7(50)*cos(50*hd)+4*cos(pi/2-50*hd);
    y(5)=l7(50)*sin(50*hd)-4*sin(pi/2-50*hd);
    x(6)=(l7(50)+8)*cos(50*hd)+4*cos(pi/2-50*hd);
    y(6)=(l7(50)+8)*sin(50*hd)-4*sin(pi/2-50*hd);
    x(7)=(l7(50)+8)*cos(50*hd)-4*cos(pi/2-50*hd);
    y(7)=(l7(50)+8)*sin(50*hd)+4*sin(pi/2-50*hd);
    x(8)=(l7(50)-8)*cos(50*hd)-4*cos(pi/2-50*hd);
    y(8)=(l7(50)-8)*sin(50*hd)+4*sin(pi/2-50*hd);
    x(9)=(l7(50)-8)*cos(50*hd)+4*cos(pi/2-50*hd);
    y(9)=(l7(50)-8)*sin(50*hd)-4*sin(pi/2-50*hd);
    x(10)=l7(50)*cos(50*hd)+4*cos(pi/2-50*hd);
    y(10)=l7(50)*sin(50*hd)-4*sin(pi/2-50*hd);
    x(11)=l7(50)*cos(50*hd);
    y(11)=l7(50)*sin(50*hd);
    x(12)=70;
    y(12)=0;
    x(13)=l6-l4*cos(theta4(50)*hd);
    y(13)=-l4*sin(theta4(50)*hd);
    x(14)=l2*cos(theta2(50)*hd)-l1*cos(50*hd);
    y(14)=l2*sin(theta2(50)*hd)-l1*sin(50*hd);
    x(15)=-l1*cos(50*hd);
    y(15)=-l1*sin(50*hd);
    x(16)=0;
    y(16)=0;
   

    plot(x,y);
    title('铰链四杆机构');
    grid on;
    hold on;
    xlabel('mm')
    ylabel('mm')
    axis([-40 140 -40 140]);  
    plot(x(1),y(1),'o');
    plot(x(2),y(2),'o');
    plot(x(11),y(11),'o');
    plot(x(12),y(12),'o');
    plot(x(13),y(13),'o');
    plot(x(14),y(14),'o');
    plot(x(15),y(15),'o');
    plot(x(16),y(16),'o');
    text(-20,-15,'l_1 ');
    text(0,-25,'l_2 ');
    text(50,-30,'l_3 ');
    text(80,-15,'l_4 ');
    text(65,30,'l_5');
    text(40,10,'l_6 ');
    text(15,30,'l_7');
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         四杆机构运动仿真
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 figure(2)
  m=moviein(20);
  j=0;
  
for n1=1:asin(6/7)*du
    j=j+1;
    clf;
   
    x(1)=0;
    y(1)=0;
    x(2)=l7(n1)*cos(n1*hd);
    y(2)=l7(n1)*sin(n1*hd);
    x(3)=(l5+l6+10)*cos(n1*hd);
    y(3)=(l5+l6+10)*sin(n1*hd);
    x(4)=l7(n1)*cos(n1*hd);
    y(4)=l7(n1)*sin(n1*hd);
    x(5)=l7(n1)*cos(n1*hd)+4*cos(pi/2-n1*hd);
    y(5)=l7(n1)*sin(n1*hd)-4*sin(pi/2-n1*hd);
    x(6)=(l7(n1)+8)*cos(n1*hd)+4*cos(pi/2-n1*hd);
    y(6)=(l7(n1)+8)*sin(n1*hd)-4*sin(pi/2-n1*hd);
    x(7)=(l7(n1)+8)*cos(n1*hd)-4*cos(pi/2-n1*hd);
    y(7)=(l7(n1)+8)*sin(n1*hd)+4*sin(pi/2-n1*hd);
    x(8)=(l7(n1)-8)*cos(n1*hd)-4*cos(pi/2-n1*hd);
    y(8)=(l7(n1)-8)*sin(n1*hd)+4*sin(pi/2-n1*hd);
    x(9)=(l7(n1)-8)*cos(n1*hd)+4*cos(pi/2-n1*hd);
    y(9)=(l7(n1)-8)*sin(n1*hd)-4*sin(pi/2-n1*hd);
    x(10)=l7(n1)*cos(n1*hd)+4*cos(pi/2-n1*hd);
    y(10)=l7(n1)*sin(n1*hd)-4*sin(pi/2-n1*hd);
    x(11)=l7(n1)*cos(n1*hd);
    y(11)=l7(n1)*sin(n1*hd);
    x(12)=70;
    y(12)=0;
    x(13)=l6-l4*cos(theta4(n1)*hd);
    y(13)=-l4*sin(theta4(n1)*hd);
    x(14)=l2*cos(theta2(n1)*hd)-l1*cos(n1*hd);
    y(14)=l2*sin(theta2(n1)*hd)-l1*sin(n1*hd);
    x(15)=-l1*cos(n1*hd);
    y(15)=-l1*sin(n1*hd);
    x(16)=0;
    y(16)=0;
   

    plot(x,y);
    title('铰链四杆机构');
    grid on;
    hold on;
    xlabel('mm')
    ylabel('mm')
    axis([-100 150 -100 150]);  
    plot(x(1),y(1),'o');
    plot(x(2),y(2),'o');
    plot(x(11),y(11),'o');
    plot(x(12),y(12),'o');
    plot(x(13),y(13),'o');
    plot(x(14),y(14),'o');
    plot(x(15),y(15),'o');
    plot(x(16),y(16),'o');
 
   
    m(j)=getframe;
 end
    axis([-40 140 -40 140]);
    xlabel('mm')
    ylabel('mm')
    movie(m)

[[it] 本帖最后由 lvmeng_2008 于 2008-6-20 19:36 编辑 [/it]]