I am getting the following error: Error using atan2 Inputs must be real. Error in (line 27) phi_vet = 2*atan2(-k​_1-sqrt(k_​1.^2+k_2.^​2-k_3.2),k​_3-k_2);

I have followed the same steps for the double rocker code as I have for the double-crank code, which seemingly works perfectly. I don't understand why the double-rocker code produces an error. Both codes are below.

Double - crank code:

clear ; close all ; clc
%% Parameters
% Bars
l_1 = 0.5;                      % Length bar a                  [m]
l_2 = 1.55;                     % Length bar b                  [m]
l_3 = 1.3;                      % Length bar c                  [m]
l_0 = 1;                        % Length bar d                  [m]
% Video
tF      = 10;                   % Final time                    [s]
fR      = 60;                   % Frame rate                    [fps]
dt      = 1/fR;                 % Time resolution               [s]
time    = linspace(0,tF,tF*fR); % Time                          [s]
% Bar a rotation
w = 2;                          % Angular velocity              [rad/s]
th_vet = w*time';
% Definitions
k_1 = -2*l_0*l_3*sin(th_vet);
k_2 = 2*l_3*(l_1-l_0*cos(th_vet));
k_3 = l_1^2 + l_0^2 - l_2^2 + l_3^2 - 2*l_1*l_0*cos(th_vet);
% phi
phi_vet = 2*atan2(-k_1-sqrt(k_1.^2+k_2.^2-k_3.^2),k_3-k_2);
% alpha
alpha_vet = atan2(-l_0*sin(th_vet)+l_3*sin(phi_vet),l_1-l_0*cos(th_vet)+l_3*cos(phi_vet));
point_A_x_cum = l_0*cos(th_vet); % Point A cummulative
point_A_y_cum = l_0*sin(th_vet); % Point A cummulative
point_B_x_cum = l_1+l_3*cos(phi_vet); % Point B cummulative
point_B_y_cum = l_3*sin(phi_vet); % Point B cummulative
%% Animation
color = cool(5); % Colormap
figure
% set(gcf,'Position',[50 50 1280 720])  % YouTube: 720p
% set(gcf,'Position',[50 50 854 480])   % YouTube: 480p
set(gcf,'Position',[50 50 640 640])     % Social
hold on ; grid on ; box on ; axis equal
set(gca,'FontName','Verdana','FontSize',18)
title('Four-bar linkage')
% Create and open video writer object
v = VideoWriter('four_bar_linkage.mp4','MPEG-4');
v.Quality   = 100;
v.FrameRate = fR;
open(v);
for i=1:length(time)
    cla
    % Angles
    th      = th_vet(i);
    phi     = phi_vet(i);
    alpha   = alpha_vet(i);
    % Bar 1
    bar_1_x = [0 l_0*cos(th)];
    bar_1_y = [0 l_0*sin(th)];
    
    % Bar 2
    bar_2_x = [l_0*cos(th) l_0*cos(th)+l_2*cos(alpha)];
    bar_2_y = [l_0*sin(th) l_0*sin(th)+l_2*sin(alpha)];
    % Bar 3
    bar_3_x = [l_1 l_1+l_3*cos(phi)];
    bar_3_y = [0 l_3*sin(phi)];
    % Trajectory
    % Point A
    plot(point_A_x_cum(1:i),point_A_y_cum(1:i),'Color',color(1,:),'LineWidth',3) % Point A trajectory
    % Point B
    plot(point_B_x_cum(1:i),point_B_y_cum(1:i),'Color',color(5,:),'LineWidth',3) % Point B trajectory
    % Fixed bar
    plot([bar_1_x(1) bar_3_x(1)],[0 0],'k','LineWidth',7)           % Bar 0
    % Bars attached to fixed bar
    plot(bar_1_x,bar_1_y,'Color',color(2,:),'LineWidth',7)          % Bar 1
    plot(bar_3_x,bar_3_y,'Color',color(4,:),'LineWidth',7)          % Bar 3
    % Bearings
    plot(bar_1_x(1),bar_1_y(1),'k^','MarkerFaceColor','k','MarkerSize',15) % Point O
    plot(bar_3_x(1),bar_3_y(1),'k^','MarkerFaceColor','k','MarkerSize',15) % Point O'
    % Bar 2
    plot(bar_2_x,bar_2_y,'Color',color(3,:),'LineWidth',7) % Bar 2
    % Points 
    plot(bar_2_x(1),bar_2_y(1),'ko','MarkerFaceColor',color(1,:),'MarkerSize',10)      % Point A
    plot(bar_3_x(end),bar_3_y(end),'ko','MarkerFaceColor',color(5,:),'MarkerSize',10)  % Point B
    %Setting axes limits
    x_range = [point_A_x_cum ;  point_B_x_cum];
    y_range = [point_A_y_cum ;  point_B_y_cum];
    grid on;
    set(gca,'xlim',[min(x_range)-0.2*(max(x_range)-min(x_range))
        max(x_range)+0.2*(max(x_range)-min(x_range))]...
        ,'ylim',[min(y_range)-0.2*(max(y_range)-min(y_range))
        max(y_range)+0.2*(max(y_range)-min(y_range))])
    set(gca,'xtick',[],'ytick',[])
    frame = getframe(gcf);
    writeVideo(v,frame);
end
close(v);

Double - Rocker code:

clear ; close all ; clc
%% Parameters
% Bars
l_1 = 1.3;                      % Length bar a                  [m]
l_2 = 1.55;                     % Length bar b                  [m]
l_3 = 0.5;                      % Length bar c                  [m]
l_0 = 1;                        % Length bar d                  [m]
% Video
tF      = 10;                   % Final time                    [s]
fR      = 60;                   % Frame rate                    [fps]
dt      = 1/fR;                 % Time resolution               [s]
time    = linspace(0,tF,tF*fR); % Time                          [s]
% Bar a rotation
w = 2;                          % Angular velocity              [rad/s]
th_vet = w*time';

% Definitions
k_1 = -2*l_1*l_3*sin(th_vet);
k_2 = 2*l_3*(l_0-l_1*cos(th_vet));
k_3 = l_0^2 + l_1^2 - l_2^2 + l_3^2 - 2*l_0*l_1*cos(th_vet);

% phi
phi_vet = double(2*atan2(-k_1-sqrt(k_1.^2+k_2.^2-k_3.^2),k_3-k_2));
% alpha
alpha_vet = atan2(-l_1*sin(th_vet)+l_3*sin(phi_vet),l_0-l_1*cos(th_vet)+l_3*cos(phi_vet));

point_A_x_cum = l_1*cos(th_vet); % Point A cummulative
point_A_y_cum = l_1*sin(th_vet); % Point A cummulative
point_B_x_cum = l_0+l_3*cos(phi_vet); % Point B cummulative

point_B_y_cum = l_3*sin(phi_vet); % Point B cummulative
%% Animation
color = cool(5); % Colormap
figure
% set(gcf,'Position',[50 50 1280 720])  % YouTube: 720p
% set(gcf,'Position',[50 50 854 480])   % YouTube: 480p
set(gcf,'Position',[50 50 640 640])     % Social
hold on ; grid on ; box on ; axis equal
set(gca,'FontName','Verdana','FontSize',18)
title('Four-bar linkage')
% Create and open video writer object
v = VideoWriter('double_rocker_four_bar_linkage.mp4','MPEG-4');
v.Quality   = 100;
v.FrameRate = fR;
open(v);
for i=1:length(time)
    cla
    % Angles
    th      = th_vet(i);
    phi     = phi_vet(i);
    alpha   = alpha_vet(i);
    % Bar 1
    bar_1_x = [0 l_1*cos(th)];
    bar_1_y = [0 l_1*sin(th)];
    
    % Bar 2
    bar_2_x = [l_1*cos(th) l_1*cos(th)+l_2*cos(alpha)];
    bar_2_y =     [l_1*sin(th) l_1*sin(th)+l_2*sin(alpha)];
    % Bar 3
    bar_3_x = [l_0 l_0+l_3*cos(phi)];
    bar_3_y = [0 l_3*sin(phi)];
    
    % Trajectory
    % Point A
    plot(point_A_x_cum(1:i),point_A_y_cum(1:i),'Color',color(1,:),'LineWidth',3) % Point A trajectory
    % Point B
    plot(point_B_x_cum(1:i),point_B_y_cum(1:i),'Color',color(5,:),'LineWidth',3) % Point B trajectory
    
    % Fixed bar
    plot([ bar_3_x(1)],[0 0],'k','LineWidth',7)           % Bar 0
    
    % Bars attached to fixed bar
    plot(bar_1_x,bar_1_y,'Color',color(2,:),'LineWidth',7)          % Bar 1
    plot(bar_3_x,bar_3_y,'Color',color(4,:),'LineWidth',7)          % Bar 3
    
    % Bearings
    plot(bar_1_x(1),bar_1_y(1),'k^','MarkerFaceColor','k','MarkerSize',15) % Point O
    plot(bar_3_x(1),bar_3_y(1),'k^','MarkerFaceColor','k','MarkerSize',15) % Point O'
    
    % Bar 2
    plot(bar_2_x,bar_2_y,'Color',color(3,:),'LineWidth',7) % Bar 2
    
    % Points 
    plot(bar_2_x(1),bar_2_y(1),'ko','MarkerFaceColor',color(1,:),'MarkerSize',10)      % Point A
    plot(bar_3_x(end),bar_3_y(end),'ko','MarkerFaceColor',color(5,:),'MarkerSize',10)  % Point B
    
    %Setting axes limits
    x_range = [point_A_x_cum ;  point_B_x_cum];
    y_range = [point_A_y_cum ;  point_B_y_cum];
    
    set(gca,'xlim',[min(x_range)-0.2*(max(x_range)-min(x_range))
        max(x_range)+0.2*(max(x_range)-min(x_range))]...
        ,'ylim',[min(y_range)-0.2*(max(y_range)-min(y_range))
        max(y_range)+0.2*(max(y_range)-min(y_range))])
    set(gca,'xtick',[],'ytick',[])
    
    frame = getframe(gcf);
    writeVideo(v,frame);
    
end
close(v);

[Rena Berman]

(Answers Dev) Restored edit

Answer 1

[chicken vector]

The term inside the square roort has some negative values, thus producing imaginary numbers.

sqrt(k_1.^2+k_2.^2-k_3.^2)

I am not really sure about what k_1, k_2 and k_3 are, but I think the problem is related to the fact that in a double rocker the two side bars do not perform a complete rotation, but your code thinks they do, trying to use angles that are physically impossible to reach.

For example I see that the starting position of your animation is already infeasible.

If you share some of the theory you used for the code, I can probably be more helpful.

[Rafay Khan]

I essentially followed the same steps shown in 2.1. Postition analysis:

http://ftp.demec.ufpr.br/disciplinas/TM253/Prof.Eduardo_Lopes/Tang_2010_4B&LE.pdf

The method involves solving the horizontal and vertical component for a four bar linkage and solving for 1 unknown angle by cancelling out the other unknown angle using trig identity. Using that method I was able to code a crank-rocker and double-crank linkage but got stuck on the double-rocker linkage.

But you might be right with your assumption, I might need to somehow add a constraint for the angles.

[chicken vector]

Yes that was exactly the problem.

You assume that your independent variable th_vet makes a full circle, but you associate th_vet to the long bar, which actually doesn't.

You should use as independent variable the angle of the short bar, which does make a full circle.

You can do this just by exchanging l1 and l3.

So if I run your code using:

l1 = 0.5;
l3 = 1.3;

instead of the opposite you obtain:

[Rafay Khan]

That's a animation for a crank rocker linkage. I'm trying to produce an animation for a double rocker linkage.

I've gotten this far, however, it appears that half of the motion is happening in the positive region whilst the other half is happening in the negative region. I'm not sure how to fix this issue.

Double rocker code:

clear; close all; clc
% Parameters
% Bars
s = 10;
l = 50;
p = 45;
q = 40;
% Video
tF      = 100;                   % Final time                    [s]
fR      = 60;                   % Frame rate                    [fps]
dt      = 1/fR;                 % Time resolution               [s]
time    = linspace(0,tF,tF*fR); % Time                          [s]
% Bar a rotation
w = 2;                          % Angular velocity              [rad/s]
theta = w*time';
a = -2*p*q*sin(theta);
b = 2*l*q - 2*p*q*cos(theta);
c = l^2 + p^2 + q^2 - s^2 - 2*l*p*cos(theta);
% Create a mask to filter out invalid values
valid_mask = (a.^2 - c.^2 + b.^2) >= 0;
% phi
phi = NaN(size(theta));
phi1 = NaN(size(theta));
phi2 = NaN(size(theta));
% Calculate the valid square root input
sqrt_input = a.^2 - c.^2 + b.^2;
sqrt_input(sqrt_input < 0) = NaN;
phi1(valid_mask) = 2*atan2(- a(valid_mask) - sqrt(sqrt_input(valid_mask)), (c(valid_mask) - b(valid_mask)));
phi2(valid_mask) = 2*atan2(- a(valid_mask) + sqrt(sqrt_input(valid_mask)), (c(valid_mask) - b(valid_mask)));
% Choose the solution that maintains continuity with the previous frame
for i = 1:length(theta)
    if i > 1 && valid_mask(i)
        % Calculate the difference between the current and previous frame angles for both solutions
        angle_diff_1 = abs(phi1(i) - phi(i-1));
        angle_diff_2 = abs(phi2(i) - phi(i-1));
        
        % Choose the solution with the smallest difference in angles to maintain continuity
        if angle_diff_1 < angle_diff_2
            phi(i) = phi1(i);
        else
            phi(i) = phi2(i);
        end
    elseif valid_mask(i)
        phi(i) = phi1(i);
    end
end
% alpha
alpha = NaN(size(theta));
alpha(valid_mask) = atan2((q*sin(phi(valid_mask)) - p*sin(theta(valid_mask))),(l + q*cos(phi(valid_mask)) - p*cos(theta(valid_mask))));
point_A_x_cum = p*cos(theta); % Point A cummulative
point_A_y_cum = p*sin(theta); % Point A cummulative
point_B_x_cum = l + q*cos(phi); % Point B cummulative
point_B_y_cum = q*sin(phi); % Point B cummulative
color = cool(5); % Colormap
figure
% set(gcf,'Position',[50 50 1280 720])  % YouTube: 720p
% set(gcf,'Position',[50 50 854 480])   % YouTube: 480p
set(gcf,'Position',[50 50 640 640])     % Social
hold on ; grid on ; box on ; axis equal
set(gca,'FontName','Verdana','FontSize',18)
title('Four-bar linkage')

% Create and open video writer object
v = VideoWriter('four_bar_linkage.mp4','MPEG-4');

The specified profile is not valid.

v.Quality   = 100;
v.FrameRate = fR;
open(v);
% Animation
for i=1:length(time)
    
    if ~valid_mask(i)
        continue;
    end
    cla;
    % Bar 1
    bar_1_x = [0 p * cos(theta(i))];
    bar_1_y = [0 p * sin(theta(i))];
    
    % Bar 2
    bar_2_x = [p * cos(theta(i)) l + q*cos(phi(i))]; % Point B cummulative
    bar_2_y = [p * sin(theta(i)) q*sin(phi(i))]; % Point B cummulative
    % Bar 3
    bar_3_x = [l point_B_x_cum(i)];
    bar_3_y = [0 point_B_y_cum(i)];
    
    % Fixed bar
    plot([0 l],[0 0],'k','LineWidth',7)           % Ground link
    
    % Bars attached to fixed bar
    plot(bar_1_x,bar_1_y,'Color',color(2,:),'LineWidth',7)          % Bar 1
    plot(bar_3_x,bar_3_y,'Color',color(4,:),'LineWidth',7)          % Bar 3
    
    % Bar 2
    plot(bar_2_x,bar_2_y,'Color',color(3,:),'LineWidth',7) % Bar 2
    %Setting axes limits
    x_range = [point_A_x_cum ;  point_B_x_cum];
    y_range = [point_A_y_cum ;  point_B_y_cum];
    
    set(gca,'xlim',[min(x_range)-0.2*(max(x_range)-min(x_range))
        max(x_range)+0.2*(max(x_range)-min(x_range))]...
        ,'ylim',[min(y_range)-0.2*(max(y_range)-min(y_range))
        max(y_range)+0.2*(max(y_range)-min(y_range))])
    set(gca,'xtick',[],'ytick',[])
    frame = getframe(gcf);
    writeVideo(v,frame);
    
    pause(0.01); % Add a small pause between frames to prevent skipping
end
close(v);

[Walter Roberson]

(Cannot test it here, MPEG-4 profile is only valid on Windows.)

[chicken vector]

You are right, my mistake.

Nonetheless it appears that the error is the same. If your short arm l2 is able to perforum a full rotation, you'd just need to change your independant variable from th_vet to alpha_vet.

But since l1 > l0 - l3 + l2 therefore l2 will make a full rotation.

You just need to reformulate the equation of k_1 k_2 k_3 wrt alpha_vet instead of theta_vet.

[chicken vector]

alpha_vet = w*time';
% Definitions
k_1 = -2*l_2*l_3*sin(alpha_vet);
k_2 = 2*l_3*(l_0-l_2*cos(alpha_vet));
k_3 = l_0^2 + l_2^2 - l_1^2 + l_3^2 - 2*l_0*l_2*cos(alpha_vet);
% phi
phi_vet = double(2*atan2(-k_1-sqrt(k_1.^2+k_2.^2-k_3.^2),k_3-k_2));
% alpha
th_vet = atan2(-l_2*sin(alpha_vet)+l_3*sin(phi_vet),l_0-l_2*cos(alpha_vet)+l_3*cos(phi_vet));

Result:

Adjust your ylim values because they don't display the full image.

[Rafay Khan]

Thank you man, i've been trying to solve this for days. I really appreciate your help, hope you have a great day!

[chicken vector]

My pleasure, have a good day aswell