MATLAB.

Text of code for reference:

clear all; clc; close all
% inputs

Wormgearrpm=900;%same as rpm of fan (400-900), the only changing input
pitchwormgear=.1; %
radiusofgear1=.25;
N2=14;
N3=56;
%4bar inputs
a=.5;%input crank
b=3;%coupler
c=3;%follower
d=.75;%fixed
theta2=0:1:(Wormgearrpm*1); %input gear crank can rotate from 0 to infinity. the input crank will be defined as the number of rotatations for 1 min

VWormGear=Wormgearrpm*pitchwormgear;%how fast the work gear moves
omegaG1=radiusofgear1*VWormGear;%omegaofGear1, as Vwormgear=Vgear
omegaG2=omegaG1;%on same axel
omegaG3=omegaG2*N2/N3;%further reduces ratio
omega2=omegaG3;% The third gear directly drives the input crank, and is constant
alpha2=0;%constant input velocity
for i=1:size(theta2,2)
K1=d/a;
K2=d/c;
K3=(a^2-b^2+c^2+d^2)/2/a/c;
A=cosd(theta2(i))-K1-K2*cosd(theta2(i))+K3;
B=-2*sind(theta2(i));
C=K1-(K2+1)*cosd(theta2(i))+K3;
theta4(i)=2*atand((-B-sqrt(B^2-4*A*C))/(2*A));
K4=d/b;
K5=(c^2-d^2-a^2-b^2)/(2*a*b);
D=cosd(theta2(i))-K1+K4*cosd(theta2(i))+K5;
E=-2*sind(theta2(i));
F=K1+(K4-1)*cosd(theta2(i))+K5;
theta3(i)=2*atand((-E-sqrt(E^2-4*D*F))/(2*D));
omega3(i)=a*omega2*sind(theta4(i)-theta2(i))/b/(sind(theta3(i)-theta4(i)));
omega4(i)=a*omega2*sind(theta2(i)-theta3(i))/c/(sind(theta4(i)-theta3(i)));

%accell

Aprime(i)=c*sind(theta4(i));
Bprime(i)=b*sind(theta3(i));
Dprime(i)=c*cosd(theta4(i));
Eprime(i)=b*cosd(theta3(i));
Cprime(i)=a*alpha2*sind(theta2(i))+a*omega2^2*cosd(theta2(i))+b*omega3(i)^2*(cosd(theta3(i)))-c*omega4(i)^2*cosd(theta4(i));
Fprime(i)=a*alpha2*cosd(theta2(i))-a*omega2^2*sind(theta2(i))-b*omega3(i)^2*(sind(theta3(i)))+c*omega4(i)^2*sind(theta4(i));
alpha3(i)=(Cprime(i)*Dprime(i)-Aprime(i)*Fprime(i))/(Aprime(i)*Eprime(i)-Bprime(i)*Dprime(i));
alpha4(i)=(Cprime(i)*Eprime(i)-Bprime(i)*Fprime(i))/(Aprime(i)*Eprime(i)-Bprime(i)*Dprime(i));
Wormgearrpmmatrix(i)=Wormgearrpm/abs(omega4(i));
end

figure
subplot(3,1,1)
plot(theta2/60,theta4) %theta 2 is derived from minutes, divide by 60 to obtain seconds.
title('Range of Fan oscilation (deg)')
ylabel('(deg)')
subplot(3,1,2)
plot(theta2/60,omega4/60) %omega 4 is in per min, dive by 60 to get per sec.
title('angular velocity of fan head rotation (deg/s)')
ylabel('(deg/s)')
subplot(3,1,3)
plot(theta2/60,alpha4/3600) %alpha 4 is in per min^2, dive by 3600 to get per sec^2.
title('angular acceleration of fan head rotation (deg/s^2)')
ylabel('(deg/s^2)')
xlabel('time of fan rotation (seconds)')

rangeofoscillation=max(theta4)-min(theta4)

MAmin=Wormgearrpm/min(abs(omega4));
MAmax=Wormgearrpm/abs(min(omega4));
disp('range of MA')
disp(MAmin)
disp(MAmax)

figure
subplot(2,1,1)
plot(theta2/60,Wormgearrpmmatrix)
title('Plot of Mechanical advantages')
ylabel('mechanical advantages')
subplot(2,1,2)
plot(theta2/60,theta4) %theta 2 is derived from minutes, divide by 60 to obtain seconds.
title('Range of Fan oscilation (deg)')
ylabel('(deg)')
xlabel('time of fan rotation (seconds)')