code to plot epicycle

17 次查看（过去 30 天）
vikneswaran M
vikneswaran M 2015-12-12
回答： Roche de Guzman 2017-4-21
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

请先登录，再回答此问题。

回答（3 个）

Matt Tearle
Matt Tearle 2015-12-13
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
Roche de Guzman 2017-4-21
%%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
Star Strider 2015-12-13

类别

在 Help Center 和 File Exchange 中查找有关 Annotations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by