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

%%%%%%%%%%%%%%%%
