code to plot epicycle

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

Answer 1

[Matt Tearle]

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)

Answer 2

[Roche de Guzman]

%%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

Answer 3

[Star Strider]

This seems to be quite a popular problem!

How do I convert spirographic values for r, R, and a to the code for a polar plot?