%Solve the large angle pendulum problem using rk4
% Izzi Urieli - December 8, 2007
clear all
THETA = 1; % index of angle [radians]
OMEGA = 2; % index of angular velocity [rads/sec]
HEIGHT = 3; % index of pendulum mass height [m]
global g length
g = 9.807; % [m/s/s} accel. due to gravity
length = input('Enter pendulum length [m]: ');
angle0 = input('Enter initial pendulum angle [degrees]: ');
period = 2*pi*sqrt(length/g);
fprintf('small angle period is %.3f seconds\n',period);
tfinal = 1.5*period;
%######################using rk4:
n = input('Enter the number of solution points: ');
dt = tfinal/(n - 1);
t = 0;
y(THETA) = angle0*pi/180; % initial pendulum angle in radians
y(OMEGA) = 0;
y(HEIGHT) = length*(1 - cos(y(THETA))); % initial pend. height
time(1) = t;
angle(1) = angle0;
height(1) = y(HEIGHT)*100; % cm (for scaling)
for i = 2:n
	[t, y, dy] = rk4('dpend',2,t,dt,y);
	time(i) = t;
	angle(i) = y(THETA)*180/pi;
    height(i) = y(HEIGHT)*100;
end
%##############################plotting:
plot(time,angle)
grid on
hold on
axis([0 tfinal -90 90])
x0 = [0,tfinal];
y0 = [0,0];
plot(x0,y0,'k-')
x1 = [period,period];
y1 = [-90,90];
plot(x1,y1,'r-')
plot(time,height, 'g-')
xlabel('time [seconds]')
ylabel('pendulum angle[degrees], height[cm] (green)')
title('pendulum angle and height vs time')