function picture
% draw a picture and rotate it in time
% Just for fun, we solve an ODE to pick the angles
% at the various times
clear all
clf
%Nice command to show that program is running
disp(['Starting at: ' datestr(now)])
%Original drawing of house with roof
x0= [2 3 3 2 2 nan 3.5 2.5 1.5]; % nan = 'not a number'
y0= [0 0 1 1 0 nan .5 1.5 .5]; % is "pen up" command in this context
% Each col of r_ref is a point on the original drawing
r_ref= [x0;y0]; %put drawing in a two row array array
n = 200; % number of drawings to make
tspan = linspace(0, 10, n); %time interval of ODE solution
p.goverl=1; % pendulum parameter is 'g over L'
z0= [pi/2 0]'; % initial conditions for ODE, a column vector
% Use ODESET for various options concerning ODE solution
% ODE45 is one of many MATLAB ODE solvers that are better
% than Euler's method or than the midpoint method
options=odeset('reltol',1e-6,'abstol',1e-6);
%The ODEs are solved in the next line
[t z] = ode45(@rhs, tspan, z0,options,p);
%Loop fo animation
for i = 1:n
pause(.01); %pause between successive drawings
theta = z(i,1); %theta is first column of ODE solution
%Rotation matrix (see text, or Linear Algebra text)
R = [ cos(theta) -sin(theta);
sin(theta) cos(theta)];
r = R* r_ref; %rotated drawing, each col. of r is a point in drawing
x=r(1,:); % first row is x coordinates
y=r(2,:); % second row is y coordinates
plot(x,y)
axis('equal'); %These two lines hold the graph axes fixed while
axis([-4 4 -4 4]); % the drawing moves around
shg
end % end animation loop
shg % SHow Graph
disp([' Done at: ' datestr(now)])% Shows that program is done
end
%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%
function zdot = rhs(t,z,p)
%ODEs for the pendulum equation
%They have nothing to do with the house drawing.
%Just used for fun to animate the house drawing
theta = z(1);
omega = z(2);
thetadot = omega; %1st ODE is just kinematics
omegadot = -p.goverl* sin(theta); %The only dynamics in this whole file
zdot = [thetadot; omegadot]; % Note that this is a column vector
end
%%%%%%%%%%%%%%%%