# code to plot epicycle

18 views (last 30 days)
vikneswaran M on 12 Dec 2015
Answered: Roche de Guzman on 21 Apr 2017
Epicycles. Write a program, epicycles.m, to draw circular motion with one epicycle, as illustrated in Figure 6.4. The object moves around the origin with a period T1 and around the center of the epicycle with a period T2. The radius of the primary orbit is R and of the epicycle is r. The angular frequency of each orbit is: ω1 = 2π/T1 ω2 = 2π/T2 The x and y coordinates of the motion are given by (convince yourself of this): x(t) = r cos(ω1t) + r cos(ω2t) y(t) = r sin(ω1t) + r sin(ω2t) Let R=5 and r =2 initially. Let T1 =1 and T2 =T1/tratio. Start with tratio = 9.25. Let the number of large cycles be Nc (initially 10, say). Let the number of time points per cycle be Ntpc (initially 200). The total time is Nc*T1 and the total number of times for which you calculate the position is Nt=Nc*Ntpc. Calculate the frequencies omega1 and omega2. For each time in the interval [0,Tf] (say loop index it goes from 1 to Nt), calculate x(it) and y(it). Plot the position of the object and a line tracing its path from the beginning to its present position. for the above question how the plot should is in the file attached with this

Matt Tearle on 13 Dec 2015
Looks like a pretty standard homework problem. You have a bunch of scalar constants given, so just define those:
>> T1 = 1;
>> Ntpc = ...
>> R = 5;
>> w1 = ...
...
Then make your vector of t values. The linspace function will help you here, given that you've been given the start and end values, and the number of values to use. (The start value is 0, the others are given in the problem: The total time is Nc*T1 and the total number of times for which you calculate the position is Nt=Nc*Ntpc)
Once you have t, calculate x and y according to the formulae (although I think you've copied them incorrectly -- the first "r" in each should probably be "R"), and plot y against x.
A simpler example that follows the same basic pattern:
>> R = 3;
>> r = 2;
>> t = linspace(0,4*pi,201);
>> x = R*sin(t).*cos(2*t);
>> y = (1-r)*sin(2*t);
>> plot(x,y)

Roche de Guzman on 21 Apr 2017
%%Epicycles
% ENGG 010 Hofstra University
% by Prof. Roche C. de Guzman
clear; clc; close('all');
%%Given
Q = input('Animate?(0 = no, 1 = yes): ');
% primary orbit
R = 5; % radius
T1 = 1; % period
% epicycle
r = 2; % radius
tratio = 9.25; % period ratio
% time and number of data points
Nc = 10; % number of cycles
Ntpc = 200; % number of time points per cycle
ti = 0; % initial time
%%Calculations
% time and number of data points
tf = Nc*T1; % final time
Nt = Nc*Ntpc; % total number of data points
t = linspace(ti,tf,Nt); % time vector
% primary orbit
omega1 = (2*pi)/T1; % angular frequency
X = R*cos(omega1*t); % x vector
Y = R*sin(omega1*t); % y vector
% epicycle
T2 = T1/tratio; % period
omega2 = (2*pi)/T2; % angular frequency
x = R*cos(omega1*t)+r*cos(omega2*t); % x vector
y = R*sin(omega1*t)+r*sin(omega2*t); % y vector
%%Display Results
doAnimation = logical(Q);
if doAnimation
for k = 1:Nt
plot(X,Y); % primary orbit
hold('on');
plot(x(k),y(k),'or','markerfacecolor',[1 0 0]); % particle
plot(x(1:k),y(1:k),'r-','markerfacecolor',[1 0 0]); % trace
axis([-(R+r) (R+r) -(R+r) (R+r)]);
title('EPICYCLES');
xlabel('x');
ylabel('y');
hold('off');
pause(0.1);
drawnow;
end
else
plot(X,Y); % primary orbit
hold('on');
plot(x,y,'r-','markerfacecolor',[1 0 0]); % trace
axis([-(R+r) (R+r) -(R+r) (R+r)]);
title('EPICYCLES');
xlabel('x');
ylabel('y');
hold('off');
end

Star Strider on 13 Dec 2015