%*************************************************************************
% Script "cam.m" to design the class example of a disk cam with          *
% reciprocating flat-face follower.                                      *
%                                                                        *
% This cam produces a full-rise of 1 inch in 120 degrees of cam rotation *
% using Simple Harmonic Motion (SHM).  This is followed by a full-return *
% using SHM of 1 inch in 240 degrees of cam rotation.                    *
%*************************************************************************

L = 1;          % Rise and return distance (in)
RAD = pi/180;   % Radian conversion factor
DEG = 180/pi;   % Degree conversion factor
STEP = 0.5;       % Cam angle stepsize (DEG)

% First find the displacement, first and second kinematic coefficient
% functions for the full-rise SHM segment.  This is termed Segment 1.

beta1d = 120;           % Cam duration for segment 1 (DEG)
beta1r = beta1d*RAD;	% Cam duration for segment 1 (RAD)

t1d = [0:STEP:beta1d]';	% Cam angle from 0 to 120 degrees; column vector (')
t1r = t1d*RAD;          % Convert to RAD

% First compute the displacement function y1 and plot it vs t1d (DEG)

y1 = (L/2)*(1-cos(pi*t1r/beta1r));
plot(t1d,y1,'k-','linewidth',2);
axis([0 beta1d min(y1)*1.1 max(y1)*1.1]);
grid;
set(gca,'fontsize',20);
xlabel('Cam angle (DEG)','fontsize',20);
ylabel('Follower motion y (in)','fontsize',20);
title('Displacement Function for Segment 1','fontsize',20);
pause;
delete(gcf);

% Next the first kinematic coefficient y1p (y1p*omega = velocity ydot)

y1p = ((pi*L)/(2*beta1r))*sin(pi*t1r/beta1r);
plot(t1d,y1p,'k-','linewidth',2);
axis([0 beta1d min(y1p)*1.1 max(y1p)*1.1]);
grid;
set(gca,'fontsize',20);
xlabel('Cam angle (DEG)','fontsize',20);
ylabel('1st Kin. Coeff. (VEL)','fontsize',20);
title('First Kinematic Coeff ("VEL") for Segment 1','fontsize',20);
pause;
delete(gcf);

% Finally the second kinematic coefficient y1pp (y1pp*omega^2 = accel yddot)

y1pp = ((pi^2*L)/(2*beta1r^2))*cos(pi*t1r/beta1r);
plot(t1d,y1pp,'k-','linewidth',2);
axis([0 beta1d min(y1pp)*1.1 max(y1pp)*1.1]);
grid;
set(gca,'fontsize',20);
xlabel('Cam angle (DEG)','fontsize',20);
ylabel('2nd Kin. Coeff. (ACCEL)','fontsize',20);
title('Second Kinematic Coeff ("ACCEL") for Segment 1','fontsize',20);
pause;
delete(gcf);

% Now plot all three on the same plot
plot(t1d,y1,'k-','linewidth',2); hold on;
plot(t1d,y1p,'r-','linewidth',2);
plot(t1d,y1pp,'b-','linewidth',2);
legend('Displ','"Vel"','"Accel"','location','best');
grid;
set(gca,'fontsize',20);
xlabel('Cam angle (DEG)','fontsize',20);
ylabel('Displ, 1st & 2nd Coeff. (in)','fontsize',20);
title('Displ, "Vel", "Accel" for Segment 1','fontsize',20);
pause;
delete(gcf);

% Next find the displacement, first and second kinematic coefficient
% functions for the full-return SHM segment.  This is termed Segment 2.
% Note that we start segment 2 from zero angle, even though it "really"
% starts at 120 degrees.  Much easier if you start from zero.

beta2d = 240;       % Cam duration for segment 2 (DEG)
beta2r = beta2d*RAD;	% Cam duration for segment 2 (RAD)

t2d = [0:STEP:240]';	% Cam angle from 0 to 240 degrees; column vector (')
t2r = t2d*RAD;          % Convert to RAD

% Same drill here; first compute the displacement function y2 and plot

y2 = (L/2)*(1+cos(pi*t2r/beta2r));
plot(t2d,y2,'k-','linewidth',2);
axis([0 beta2d min(y2)*1.1 max(y2)*1.1]);
grid;
set(gca,'fontsize',20);
xlabel('Cam angle (DEG)','fontsize',20);
ylabel('Follower motion y (in)','fontsize',20);
title('Displacement Function for Segment 2','fontsize',20);
pause;
delete(gcf);

% Next compute the first and second kinematic coefficients yp and ypp and
% plot

y2p = (-(pi*L)/(2*beta2r))*sin(pi*t2r/beta2r);
plot(t2d,y2p,'k-','linewidth',2);
axis([0 beta2d min(y2p)*1.1 max(y2p)*1.1]);
grid;
set(gca,'fontsize',20);
xlabel('Cam angle (DEG)','fontsize',20);
ylabel('1st Kin. Coeff. (VEL)','fontsize',20);
title('First Kinematic Coeff ("VEL") for Segment 2','fontsize',20);
pause;
delete(gcf);

y2pp = (-(pi^2*L)/(2*beta2r^2))*cos(pi*t2r/beta2r);
plot(t2d,y2pp,'k-','linewidth',2);
axis([0 beta2d min(y2pp)*1.1 max(y2pp)*1.1]);
grid;
set(gca,'fontsize',20);
xlabel('Cam angle (DEG)','fontsize',20);
ylabel('2nd Kin. Coeff. (ACCEL)','fontsize',20);
title('Second Kinematic Coeff ("ACCEL") for Segment 2','fontsize',20);
pause;
delete(gcf);

% Now plot all three on the same plot

plot(t2d,y2,'k-','linewidth',2); hold on;
plot(t2d,y2p,'r-','linewidth',2);
plot(t2d,y2pp,'b-','linewidth',2);
axis([0 beta2d min([min(y2) min(y2p) min(y2pp)])*1.1 ...
    max([max(y2) max(y2p) max(y2pp)])*1.1]);
legend('Displ','"Vel"','"Accel"','location','best');
grid;
set(gca,'fontsize',20);
xlabel('Cam angle (DEG)','fontsize',20);
ylabel('Displ, 1st & 2nd Coeff. (in)','fontsize',20);
title('Displ, "Vel", "Accel" for Segment 2','fontsize',20);
pause;
delete(gcf);

% Even though the minimum follower radius may have been computed, one can
% get that value from the numerical values just computed.  We want to
% find the maximum value of the magnitude of y1p (the max velocity will
% occur during Segment 1).

r_min = max(abs(y1p));  % Find the maximum of the absolute value of y1p

s = sprintf('\nMinimum follower radius is: %.2f inches\n', r_min); % Display
disp(s);

% Next find the minimum base circle radius.  We will find the minimum base
% circle radius so that the cam profile has no cusps, then give the user
% the option of entering the desired base circle radius.  For no cusps, the
% minimum base circle radius Ro is equal to the maximum value of (-y-ypp).
% Since we've computed both of those, finding the max should be easy.  We
% do need to check both segments 1 and 2, although intuitively segment 1
% has the fastest motion, so it should determine.

R(1) = max(-y1-y1pp);
R(2) = max(-y2-y2pp);

Ro = max(R);

s = sprintf('Minimum base circle radius is: %.4f inches\n', Ro); % Display
disp(s);

Ro = input('Enter desired base circle radius: ');

% Now we can calculate the cam profile using the xc and yc equations from
% the cam/follower geometric analysis.  However, this will require
% concatenating the two motions.  Whenever one does this, there will be a
% duplicate point between the two that should be removed.

N1 = length(t1r);   % Length of segment 1 (samples)
N2 = length(t2r);   % Length of segment 2 (samples)

% Now to "cement" the two segments together.  Since they are both column
% vectors, we can just stack them one on top of the other.  Before we do
% that, however, remove the last point from the first segment, since it's a
% duplicate

y1 = y1(1:N1-1);    % Remove the last point from 
y1p = y1p(1:N1-1);	% Segment 1 for all three
y1pp = y1p(1:N1-1);	% motion functions

% Now stack them on top of one another

y = [y1;y2];    % The ";" operator means new row, so the two cols are stacked
yp = [y1p;y2p];
ypp = [y1pp;y2pp];

% Finally create a vector of cam angle from start (0 DEG) to finish (360
% DEG).  The stepsize is 1 DEG; if another step

thetad = [0:STEP:360]';   % Total cam angle (DEG), STEP increment
thetar = thetad*RAD;  % Convert to RAD

% Now plot the entire motion: displacement, "velocity", and "acceleration"

plot(thetad,y,'k-','linewidth',2);
axis([0 360 min(y)*1.1 max(y)*1.1]);
grid;
set(gca,'fontsize',20);
xlabel('Cam angle (DEG)','fontsize',20);
ylabel('Follower motion y (in)','fontsize',20);
title('Displacement Function for Entire Motion','fontsize',20);
pause;
delete(gcf);

plot(thetad,yp,'k-','linewidth',2);
axis([0 360 min(yp)*1.1 max(yp)*1.1]);
grid;
set(gca,'fontsize',20);
xlabel('Cam angle (DEG)','fontsize',20);
ylabel('1st Kin. Coeff. (VEL)','fontsize',20);
title('First Kinematic Coeff ("VEL") for Entire Motion','fontsize',20);
pause;
delete(gcf);

plot(thetad,ypp,'k-','linewidth',2);
axis([0 360 min(ypp)*1.1 max(ypp)*1.1]);
grid;
set(gca,'fontsize',20);
xlabel('Cam angle (DEG)','fontsize',20);
ylabel('2nd Kin. Coeff. (ACCEL)','fontsize',20);
title('Second Kinematic Coeff ("ACCEL") for Entire Motion','fontsize',20);
pause;
delete(gcf);

% This is enough.  The parameters can now be passed to the camprofile.m
% function, which will produce the (xc,yc) profile points.  These can
% easily plotted.

disp('Files y, yp, theta created...ready for camprofile.m to');
disp('generate (x,y) coordinates of cam profile.');
disp(' ');  % Blank line