Simulation and Primary and Secondary Force calculation of a V-Twin engine.

%Motion analysis and Force analysis of a V-Twin engine
clear all
clc
close all
%angle between the centrelines of cylinder
alpha=30
%crank angle from engine centreline
Theta=linspace(0,360,50);
%crank and connecting rod lenghts(m)
r=0.05;
l=0.25;
%Reciprocating mass (kg)
Mrec=2
%Engine Speed (rpm)
N=2000
%Angular speed
w=(2*pi*N/60)
x0=0;
y0=0;

ct=1;
for i=1:length(Theta)
 theta=Theta(i);   
%crank angles of individual cylinders
theta1=theta-alpha;
theta2=theta+alpha;
%coordinates

x1=r*sind(theta);
y1=r*cosd(theta);
beta1=asind((r/l)*sind(theta-alpha));
x2=x1+l*sind(alpha-beta1);
y2=y1+l*cosd(alpha-beta1);

beta2=asind((r/l)*sind(theta2));
m=(l/sind(theta2))*sind(beta2+theta2);
x3=-m*sind(alpha);
y3=m*cosd(alpha);
x4=(l+r)*sind(alpha);
y4=(l+r)*cosd(alpha);
x5=-(l+r)*sind(alpha);
y5=(l+r)*cosd(alpha);
th(i)=theta;
%Primary Force
F_P(i)=2*Mrec*w^2*r*((cosd(theta))^2*(cosd(alpha))^4+(sind(theta))^2*(sind(alpha))^4)^0.5;
%Secondary Force
F_S(i)=2*Mrec*w^2*r*(r/l)*((cosd(2*theta))^2*(cosd(alpha))^2*(cosd(2*alpha))^2+(sind(2*theta))^2*(sind(alpha))^2*(sind(2*alpha))^2)^0.5;
%plotting
M(ct)=getframe(gcf)
plot([x0,x1],[y0,y1],[x1,x2],[y1,y2],[x3,x1],[y3,y1],[x4,x0],[y4,y0],[x0,x5],[y0,y5])
axis([-0.2 0.2 -0.1 0.3])
pause(0.1)
ct=ct+1;
end

% Maximum and Minimum Primary Forces
 FPmax=max(F_P)
for j=1:length(F_P)
if F_P(j) ==(FPmax)
   THETAPmax=th(j)
   
end
end

FPmin=min(F_P)
for j=1:length(F_P)
if F_P(j) ==(FPmin)
   THETAPmin=th(j)
end
end
% Maximum and Minimum Secondary Forces
FSmax=max(F_S)
for j=1:length(F_S)
if F_S(j) ==(FSmax)
   THETASmax=th(j)
end
end
FSmin=min(F_S)
for j=1:length(F_S)
if F_S(j) ==(FSmin)
   THETASmin=th(j)
end
end

movie(M)
Videofile=VideoWriter('V_twinEngine.avi','uncompressed AVI');
open(Videofile)
writeVideo(Videofile,M)
close(Videofile)