function [t,theta,thetadot] = drvpend(gamma, omega0, epsilon, h, ic, finaltime, dt)

% Matlab script to interface drvpend.c
% Integrate equation of motion for the damped driven pendulum.
% Usage: [t,theta,thetadot] = drvpend(gamma, omega0, epsilon, h, ic, finaltime);
%
%
% Inputs:       gamma   damping coefficient
%               omega0  natural frequency
%               epsilon driving frequency - natural frequency
%               h       amplitude of drive
%               ic      initial conditions, eg [0.0 3.0]
%               finaltime  length of time of integration
%               dt	time increment
%
% Outputs:      t     
%               theta
%               thetadot
%
% Example:
%
% [t,theta,thetadot] = drvpend(0.0, 1.0, 0.0, 0.0, [0.0 0.01], 1000, 0.01);
% 
% Note that dt is now set automatically in the drvpend.c code.  If you
% wish to change that, please modify the c-code and recompile.

drvpend_command_string = sprintf('./drvpend %f %f %f %f %f %f %f %f', gamma, omega0, epsilon, h, ic(1), ic(2), finaltime, dt);

'Computing trajectory... please wait...'
unix(drvpend_command_string);

'Loading data into Matlab... please wait...'
N=ceil(finaltime/dt);

fid = fopen('t','r');
t = fread (fid, N, 'float');
fclose(fid);

fid = fopen('theta','r');
theta = fread (fid, N, 'float');
fclose(fid);

fid = fopen('thetadot','r');
thetadot = fread (fid, N, 'float');
fclose(fid);