帮忙看看关于坐标变换和物理问题求解的程序

aaachwei

您好！下面这个程序有时候能出图形有时候就有错误，不知道为什么。我MATLAB学的不太好，麻烦帮帮我可以吗？谢谢！
alpha = 0;
beta = 0;
gamma = 0;
x=-200:10:200;
y=-200:10:200;
z=0;
T = [ cos(alpha)*cos(beta)  sin(alpha)*sin(beta)*cos(gamma)-cos(alpha)*sin(gamma)  cos(alpha)*sin(beta)*cos(gamma)+sin(alpha)*sin(gamma);
      cos(alpha)*sin(beta)  sin(alpha)*sin(beta)*sin(gamma)+cos(alpha)*cos(gamma)  cos(alpha)*sin(beta)*sin(gamma)-sin(gamma)*cos(alpha);
         -sin(beta)                     cos(beta)*sin(gamma)                                                        cos(gamma)*cos(beta)];

A10 = [-400/2;320/2;0];    %运动平台上各球铰点在动坐标系下的矢量表示
A20 = [-400/2;-320/2;0];
A30 = [0;-400/2;320/2];
A40 = [0;-400/2;-320/2];
A50 = [320/2;0;-400/2];
A60 = [-320/2;0;-400/2];

C1 = [-(400/2+500+400);320/2;0];        %各分支与基座连接点在定坐标系下的矢量表示
C2 = [-(400/2+500+400);-320/2;0];
C3 = [0;-(400/2+500+400);320/2];
C4 = [0;-(400/2+500+400);-320/2];
C5 = [320/2;0;-(400/2+500+400)];
C6 = [-320/2;0;-(400/2+500+400)];

for i=1:length(y)
    for j=1:length(x)
        
        A1 = T*A10+[x(j);y(i);z];
        A2 = T*A20+[x(j);y(i);z];
        A3 = T*A30+[x(j);y(i);z];
        A4 = T*A40+[x(j);y(i);z];
        A5 = T*A50+[x(j);y(i);z];
        A6 = T*A60+[x(j);y(i);z];
        
        l1 = 500+x(j)-sqrt(500^2-(y(i))^2-z^2);
        l2 = 500+x(j)-sqrt(500^2-(y(i))^2-z^2);
        l3 = 500+y(i)-sqrt(500^2-(x(i))^2-z^2);
        l4 = 500+y(i)-sqrt(500^2-(x(i))^2-z^2);
        l5 = 500+z-sqrt(500^2-(x(j))^2-(y(i))^2);
        l6 = 500+z-sqrt(500^2-(x(j))^2-(y(i))^2);
        
        
        B1 = [-(400/2+500-l1);320/2;0];        
        B2 = [ -(400/2+500-l2);-320/2;0];
        B3 = [0;-(400/2+500-l3);320/2];
        B4 = [0;-(400/2+500-l4);-320/2];
        B5 = [320/2;0;-(400/2+500-l5)];
        B6 = [-320/2;0;-(400/2+500-l6)];
        zeta1 = (B1-A1)/500;
        zeta2 = (B2-A2)/500;
        zeta3 = (B3-A3)/500;
        zeta4 = (B4-A4)/500;
        zeta5 = (B5-A5)/500;
        zeta6 = (B6-A6)/500;
        
        eta1 = (C1-B1)/abs(400+l1);
        eta2 = (C2-B2)/abs(400+l2);
        eta3 = (C3-B3)/abs(400+l3);
        eta4 = (C4-B4)/abs(400+l4);
        eta5 = (C5-B5)/abs(400+l5);
        eta6 = (C6-B6)/abs(400+l6);
        
        delta1 = acos(dot(zeta1,T*eta1));
        delta2 = acos(dot(zeta2,T*eta2));
        delta3 = acos(dot(zeta3,T*eta3));
        delta4 = acos(dot(zeta4,T*eta4));
        delta5 = acos(dot(zeta5,T*eta5));
        delta6 = acos(dot(zeta6,T*eta6));
        
        theta1 = acos(dot(zeta1,eta1));
        theta2 = acos(dot(zeta2,eta2));
        theta3 = acos(dot(zeta3,eta3));
        theta4 = acos(dot(zeta4,eta4));
        theta5 = acos(dot(zeta5,eta5));
        theta6 = acos(dot(zeta6,eta6));
        
        
        if ((l1>=-200) && (l1<=200) && (l2>=-200) && (l2<=200) && (l3>=-200) && (l3<=200) && (l4>=-200) && (l4<=200) && (l5>=-200) && (l5<=200) && (l6>=-200) && (l6<=200) && (delta1 >= 0) && (delta1 <= pi/6) && (delta2 >= 0) && (delta2 <= pi/6) && (delta3 >= 0) && (delta3 <= pi/6) && (delta4 >= 0) && (delta4 <= pi/6) && (delta5 >= 0) && (delta5 <= pi/6) && (delta6 >= 0) && (delta6 <= pi/6) && (theta1 >= 0) && (theta1 <= pi/6) && (theta2 >= 0) && (theta2 <= pi/6) && (theta3 >= 0) && (theta3 <= pi/6) && (theta4 >= 0) && (theta4 <= pi/6) && (theta5 >= 0) && (theta5 <= pi/6) && (theta6 >= 0) && (theta6 <= pi/6))
         
          G=[(zeta1)' cross(A1,zeta1)';
             (zeta2)' cross(A2,zeta2)';
             (zeta3)' cross(A3,zeta3)';
             (zeta4)' cross(A4,zeta4)';
             (zeta5)' cross(A5,zeta5)';
             (zeta6)' cross(A6,zeta6)'];
          J2=[dot(zeta1,eta1) 0 0 0 0 0;
              0 dot(zeta2,eta2) 0 0 0 0;
              0 0 dot(zeta3,eta3) 0 0 0;
              0 0 0 dot(zeta4,eta4) 0 0;
              0 0 0 0 dot(zeta5,eta5) 0;
              0 0 0 0 0 dot(zeta6,eta6)];  
      
           if det(J2)~= 0&det(G)~= 0
                J = (inv(G)*J2);
                Jv=J(1:3,:);
                s=svd(Jv);
               Kv(i,j)=max(s);               
           end
        end
      
    end
end

surfc(x,y,Kv)

lengyunfeng

你的x、y维数相同，都是1*41的double型数据；你循环语句里的Kv照理来说应该是41*41型数组，但受你循环体内两次if句型条件控制，使得你的Kv最终是一个38*41型数组，这样你用surfc(x,y,Kv)就会显示矩阵与向量维数不一致，因为size(Kv)必须要等于（length(y),length(x)）。自己再改改吧

ETsy

真是复杂，建议写点注释

qibbxxt

有没有别的方法啊，if语句有点太长了

aaachwei

x,y的维数是1*41，可z（Kv）是x，y的函数，怎么是38*41呢？能不能教我一下，这个维数要怎么算？我刚学，不太会。谢谢了