Hi, I have fonud this code online, and from I want to try and plot theta1 vs time and theta2 vs time, however I am struggling to manage it. Can anybody help? thanks
function double_pendulum_simulation(theta1_0,theta2_0,L1,L2,m1,m2,g,tail)%DOUBLE_PENDULUM_SIMULATION Plots double pendulum dynamics until stopped
%
% Simulates and plots the dynamics behavior of the double simple
% pendulum until the user presses a key or clicks in the figure.
% If any of the input parameters is not given, it will be assigned a
% default value.
%% Example uses:
% DOUBLE_PENDULUM_SIMULATION
% DOUBLE_PENDULUM_SIMULATION(THETA1,THETA2,L1,L1,M1,M2,G,TAIL)
% THETA1_0 and THETA2_0 initial angles (default 3*pi/4 and pi/2)
% L1 and L2 rod lengths (default 3 and 2)
% M1 and M2 bob masses (default 2 and 3)
% G acceleration of gravity (default 9.81)
% TAIL lengh of visualization tail (default 100)
%% Author: vand@dtu.dk, 2014
if nargin<8 || isempty(tail) tail = 20;endif nargin<7 || isempty(g) g = 9.81;endif nargin<6 || isempty(m2) m2 = 10;endif nargin<5 || isempty(m1) m1 = 2;endif nargin<4 || isempty(L2) L2 = 1;endif nargin<3 || isempty(L1) L1 = 1;endif nargin<2 || isempty(theta2_0) theta2_0 = pi/2;endif nargin<1 || isempty(theta1_0) theta1_0 = 0;end% preparing for loop until user either keypresses or clicks
global USER_RESPONDEDUSER_RESPONDED = 0;figureset(gcf,'WindowKeyPressFcn',@userRespondFcn,'WindowButtonDownFcn',... @userRespondFcn,'DeleteFcn',@userRespondFcn)% initial state
x = [mod(theta1_0,2*pi),mod(theta2_0,2*pi),0,0];t = 0;i_end = 0.02; % initial estimate of a average iteration duration
tail = repmat(L1*[sin(x(1)),-cos(x(1))],[tail,2]) + ... repmat([0 0 L2*[sin(x(2)),-cos(x(2))]],[tail,1]); % initial tail
% ode function
double_penudlum = @(t,x)double_pendulum_system(x,L1,L2,m1,m2,g);% preparing for loop
axis xy equal, box on, hold onaxis(1.1*[-1 1 -1 1]*(L1+L2))[r1,r2] = bob_drawing_scale(m1,m2,L1,L2); % scale for drawing bobs
iter = 0;ticwhile ~USER_RESPONDED [t,x] = ode45(double_penudlum,t(end)*[1 0.5 0] + i_end*[0 0.5 1] ,x(end,:)'); tail = patch_double_pendulum(t,x(end,:),L1,L2,r1,r2,tail); iter = iter+1; i_end = max(toc*(1+1/iter),t(end)+2*eps); % end time for next iteration
endendfunction userRespondFcn(~,~)% callback function which execute when user responds
global USER_RESPONDEDUSER_RESPONDED = 1;endfunction dx = double_pendulum_system(x,L1,L2,m1,m2,g)% a system of differential equations defining a double pendulum
% from http://www.myphysicslab.com/dbl_pendulum.html
theta1 = x(1);theta2 = x(2);omega1 = x(3);omega2 = x(4);dtheta1 = omega1;dtheta2 = omega2;domega1 = (-g*(2*m1+m2)*sin(theta1)-m2*g*sin(theta1-2*theta2)-... 2*sin(theta1-theta2)*m2*(omega2^2*L2+omega1^2*L1*cos(theta1-theta2)))... /(L1*(2*m1+m2-m2*cos(2*theta1-2*theta2)));domega2 = (2*sin(theta1-theta2)*(omega1^2*L1*(m1+m2)+... g*(m1+m2)*cos(theta1)+omega2^2*L2*m2*cos(theta1-theta2)))... /(L2*(2*m1+m2-m2*cos(2*theta1-2*theta2)));dx = [dtheta1;dtheta2;domega1;domega2];endfunction tail = patch_double_pendulum(t,x,L1,L2,r1,r2,tail)% using patch to plot pendulum state
cla% plotting tail
alpha = linspace(0,1,size(tail,1)+1)';patch([tail(:,1);NaN],[tail(:,2);NaN],0,'EdgeColor','r','FaceColor',... 'none','FaceVertexAlphaData',alpha,'EdgeAlpha','interp','LineWidth',2);patch([tail(:,3);NaN],[tail(:,4);NaN],0,'EdgeColor','r','FaceColor',... 'none','FaceVertexAlphaData',alpha,'EdgeAlpha','interp','LineWidth',2);% plotting rods
theta1 = x(1);theta2 = x(2);xm1 = L1*sin(theta1);ym1 = -L1*cos(theta1);xm2 = xm1 + L2*sin(theta2);ym2 = ym1 - L2*cos(theta2);patch([0, xm1, xm2, NaN],[0, ym1, ym2, NaN],0,'EdgeColor','b',... 'FaceColor','none','LineWidth',2)% plotting bobs
p = linspace(0,2*pi,17);sint = sin(p);cost = cos(p);patch(xm1+r1*cost,ym1+r1*sint,0,'EdgeColor','b','FaceColor','b')patch(xm2+r2*cost,ym2+r2*sint,0,'EdgeColor','b','FaceColor','b')title(sprintf('time = %0.1f',t(end)))drawnowtail = [tail(2:end,:);xm1,ym1,xm2,ym2];endfunction [r1,r2] = bob_drawing_scale(m1,m2,L1,L2)% finding a scale for plotting the size of the bobs: the bigger bob has
% a radius which is a given fraction of the length of the shorter rod.
r_max = max(m1^(1/3),m2^(1/3));l_min = min(L1,L2);scale = 0.1*l_min/r_max;r1 = scale*m1^(1/3); % mass is proportional with the volume
r2 = scale*m2^(1/3);end
Best Answer