30 views (last 30 days)

Show older comments

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

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!