Theo Jansen Mechanism Matlab Code

I shared my m-code for the Theo Jansen mechanism awhile ago (LINK) but one of our readers recently created his own version of the program and offered to share it with us. Truong Duc Binh says:

"I programed for two legs and now I am making a computer interface with it, user can input parameters for Jansen walking machine. I divide this window form to two modes, you can simulate with constant alpha, or variable alpha."

Here is the complete m-code:

clear all

%Init link's length
L1=286;
L2=100;
L3=400;
L4=275;
L5=400;
L6=295;
L7=285;
L8=275;
L9=400;
L10=290;
L11=400;
L12=280;
L13=80;

alphaMax = pi/36; % 5 degree
alphaStep = pi/36/2; % 5/2 degree
alpha = -alphaStep;
k=3;

%To get the video:
mov=avifile('CoCau4KhauBanLe.avi','COMPRESSION','Cinepak');% None
for (a =1:1:(k+1))
    if a <= (k+1-rem(k+1,2))/2
        alpha = alpha + alphaStep;
    else
        alpha = alpha - alphaStep;
    end
    if alpha > alphaMax
        alpha = alphaMax;
    end
    if alpha<(-alphaMax)
        alpha = (-alphaMax);
    end
    for theta=0:pi/20:2*pi
    %for theta=0:(-pi/20):(-2*pi) % Quay cung chieu kim dong ho

    % Calculate angles:
    theta2 = pi - theta;
    theta3 = pi - theta + alpha;
    ssquared=(L2)^2+(L1)^2-(2*L1*L2*cos(theta3));
    s=sqrt(ssquared);
    beta=asin((L2*sin(theta3))/s);
    psi=acos(((L3^2)+(s^2)-(L4^2))/(2*L3*s));
    BCO2 = acos(((L5)^2+(L6)^2 - (L4)^2)/(2*L5*L6));
    CBO2 = acos(((L5)^2+(L4)^2 - (L6)^2)/(2*L5*L4));
    lambda=acos(((L4)^2+s^2-(L3)^2)/(2*(L4)*s));
    gamma = lambda - (BCO2 + CBO2)+ alpha + beta;
    O2AD = acos(((L9)^2+ssquared-(L8)^2)/(2*(L9)*s));
    omega = O2AD + alpha + beta;
    AO2D = acos(((L8)^2+ssquared-(L9)^2)/(2*(L8)*s));
    CO2D = pi - gamma - AO2D;
    CDsquared = (L6)^2+(L8)^2-2*L6*L8*cos(CO2D);
    CD = sqrt(CDsquared);
    O2CD = asin(L8*sin(CO2D)/CD);
    tempE = gamma - O2CD;
    CDE = acos(((L10)^2+CDsquared-(L7)^2)/(2*(L10)*CD));
    angleE = tempE + CDE;
    DEF = acos(((L10)^2+(L11)^2 - (L12)^2)/(2*L10*L11));
    tempF = DEF - angleE;
    DFE = acos(((L12)^2+(L11)^2 - (L10)^2)/(2*L12*L11));
    angleF = pi - tempF - DFE;

    %Find the points:
    O1=[0,0];
    O2=[-L1*cos(alpha),-L1*sin(alpha)];
    A=[L2*cos(theta) L2*sin(theta)];
    B=[A(1)-L3*cos(psi-alpha-beta) A(2)+L3*sin(psi-alpha-beta)];
    C=[O2(1)-L6*cos(gamma) O2(2)-L6*sin(gamma)];
    D=[A(1)-L9*cos(omega) A(2)-L9*sin(omega)];
    E=[D(1)-L10*cos(angleE) D(2)-L10*sin(angleE)];
    F=[D(1)-L12*cos(angleF) D(2)-L12*sin(angleF)];
    G=[F(1)-L13*cos(angleF) F(2)-L13*sin(angleF)];

    %************************LEFT********************************
    plot ([O1(1),A(1)],[O1(2),A(2)],'g','linewidth',3)
    hold on
    plot ([O2(1),O1(1)],[O2(2),O1(2)],'black','linewidth',3)
    plot ([A(1),B(1)],[A(2),B(2)],'b','linewidth',3)
    plot ([O2(1),B(1)],[O2(2),B(2)],'r','linewidth',3)
    plot ([O2(1),C(1)],[O2(2),C(2)],'r','linewidth',3)
    plot ([B(1),C(1)],[B(2),C(2)],'r','linewidth',3)
    plot ([O2(1),D(1)],[O2(2),D(2)],'b','linewidth',3)
    plot ([A(1),D(1)],[A(2),D(2)],'b','linewidth',3)
    plot ([C(1),E(1)],[C(2),E(2)],'b','linewidth',3)
    plot ([D(1),E(1)],[D(2),E(2)],'r','linewidth',3)
    plot ([D(1),F(1)],[D(2),F(2)],'r','linewidth',3)
    plot ([E(1),F(1)],[E(2),F(2)],'r','linewidth',3)
    plot ([G(1),F(1)],[G(2),F(2)],'b','linewidth',3)
    Gx(a)=G(1);
    Gy(a)=G(2);

    %*************************RIGHT******************************
    theta4 = theta + alpha;
    s2squared=(L2)^2+(L1)^2-(2*L1*L2*cos(theta4));
    s2=sqrt(s2squared);
    beta2=asin((L2*sin(theta4))/s2);
    angle2=acos(((L3^2)+((s2)^2)-(L4^2))/(2*L3*(s2)));
    B1C1O3 = acos(((L5)^2+(L6)^2 - (L4)^2)/(2*L5*L6));
    C1B1O3 = acos(((L5)^2+(L4)^2 - (L6)^2)/(2*L5*L4));
    lambda2=acos(((L4)^2+(s2)^2-(L3)^2)/(2*(L4)*(s2)));
    gamma2 = lambda2 - (B1C1O3 + C1B1O3)+ alpha + beta2;
    O3A1D1 = acos(((L9)^2+s2squared-(L8)^2)/(2*(L9)*(s2)));
    omega2 = O3A1D1 + alpha + beta2;
    A1O3D1 = acos(((L8)^2+s2squared-(L9)^2)/(2*(L8)*(s2)));
    C1O3D1 = 2*pi - A1O3D1 - lambda2 - (pi - B1C1O3 - C1B1O3);
    C1D1squared = (L6)^2+(L8)^2-2*L6*L8*cos(C1O3D1);
    C1D1 = sqrt(C1D1squared);
    O3C1D1 = asin(L8*sin(C1O3D1)/(C1D1));
    tempE2 = gamma2 - O3C1D1;
    C1D1E1 = acos(((L10)^2+C1D1squared-(L7)^2)/(2*(L10)*(C1D1)));
    angleE2 = tempE2 + C1D1E1;
    D1E1F1 = acos(((L10)^2+(L11)^2 - (L12)^2)/(2*L10*L11));
    tempF2 = D1E1F1 - angleE2;
    D1F1E1 = acos(((L12)^2+(L11)^2 - (L10)^2)/(2*L12*L11));
    angleF2 = tempF2 + D1F1E1;

    %Find the points:
    O3=[L1*cos(alpha),-L1*sin(alpha)];%
    B1=[A(1)+L3*cos(angle2-alpha-beta2) A(2)+L3*sin(angle2-alpha-beta2)];%
    C1=[O3(1)+L6*cos(gamma2) O3(2)-L6*sin(gamma2)];
    D1=[A(1)+L9*cos(omega2) A(2)-L9*sin(omega2)];
    E1=[D1(1)+L10*cos(angleE2) D1(2)-L10*sin(angleE2)];
    F1=[D1(1)-L12*cos(angleF2) D1(2)-L12*sin(angleF2)];
    G1=[F1(1)-L13*cos(angleF2) F1(2)-L13*sin(angleF2)];

    plot ([O3(1),O1(1)],[O3(2),O1(2)],'black','linewidth',3)
    plot ([A(1),B1(1)],[A(2),B1(2)],'b','linewidth',3)
    plot ([O3(1),B1(1)],[O3(2),B1(2)],'r','linewidth',3)
    plot ([O3(1),C1(1)],[O3(2),C1(2)],'r','linewidth',3)
    plot ([B1(1),C1(1)],[B1(2),C1(2)],'r','linewidth',3)
    plot ([O3(1),D1(1)],[O3(2),D1(2)],'b','linewidth',3)
    plot ([A(1),D1(1)],[A(2),D1(2)],'b','linewidth',3)
    plot ([C1(1),E1(1)],[C1(2),E1(2)],'b','linewidth',3)
    plot ([D1(1),E1(1)],[D1(2),E1(2)],'r','linewidth',3)
    plot ([D1(1),F1(1)],[D1(2),F1(2)],'r','linewidth',3)
    plot ([E1(1),F1(1)],[E1(2),F1(2)],'r','linewidth',3)
    plot ([G1(1),F1(1)],[G1(2),F1(2)],'b','linewidth',3)
    G1x(a)=G1(1);
    G1y(a)=G1(2);
   
    plot(Gx,Gy)
    plot(G1x,G1y)
  
    %************************************************************
    axis([-800 800 -800 400])
    hold off
    a=a+1;
    pause(.01)
    M=getframe;
    mov=addframe(mov,M);
    end
end
mov=close(mov);

Thanks again to Truong Duc Binh!