function [tiarr,tharr,vearr,vfr,omef,freq,volp] = windpump(inpI,the0,ome0,tfI)

% Calculates and plots the dynamics of a windpump entrained by a 
% single-cylinder, single-acting piston-slider-crank mechanism 
% of negligible inertia and friction, raising water a known height through
% a line without losses, using a frictionless chain of negligible inertia
% to entrain the pump from the windmill. Sprockets and windmill inertia are
% considered.
%
% Michel Louge 2/11/2010
%
%   input   -   inpI	=   parameter input column vector including
%                           	rcrank  =   cranshaft radius (m)
%                               lcrank  =   connecting rod length (m)
%                               R0      =   radius of small sprocket (on windmill, m)
%                               R1      =   radius of large sprocket (on pump, m)
%                               A       =   piston cross-section (m^2)
%                               I0      =   moment of inertia of windmill, shaft and sprocket (kg.m^2)
%                               I1      =   moment of inertia of large sprocket (kg.m^2)
%                               om      =   windmill ang. vel. at peak torque (rad/s)
%                               Tm      =   max. windmill torque (N.m)
%                                               Note:   assume that windmill torque vs. ang. vel.
%                                               is an inverted parabola passing through
%                                               the origin.
%                               rgh     =   required water potential energy/u. vol. (J/m^3)
%                               
%           -   the0	=   initial crank angle (radians); = 0 at top dead center (TDC)
%           -   ome0    =   initial angular velocity of the windmill (rad/s)
%           -   tfI     =   final simulation time (s)
%
%   output  -   tiarr 	=   array of times (s)
%           -   tharr 	=   array of crank angles modulo 2pi
%           -   vearr 	=   array of crank angular velocities
%           -   vfr     =   array of volume flow rates (liter/s)
%           -   omef    =   last angular velocity when crank angle was zero (rad/s)
%                                  (zero if system did not complete a whole cycle)
%           -   freq    =   average cycle frequency (Hz) 
%                                   (zero if system did not complete a whole cycle)
%           -   volp    =   total volume pumped (liter)
%
%   global variables    if present, postscript I denotes an input variable destined to
%                       become global (avoids confusing Matlab)
%   
%           -   lr      =   ratio of connecting rod length to crankshaft radius
%           -   Tmd     =   dim(-) torque parameter
%           -   loa     =   dim(-) load parameter
%           -   reduc   =   reduction ratio R1/R0
%
%   to access the proper directory on Louge's computer, type
%   cd ~/Michel_files/Classwork/MAE_225/Spring_2010/projects/wind-pump/Matlab_dynamics
%
% Scales for evaluating dimensionless variables:
%
%   time        made dim(-) with    om
%   torque      made dim(-) with    Tm
%
% Examples:
%                               variable    example below
%
%                           	rcrank  =   0.03 m
%                               lcrank  =   0.12 m
%                               R0      =   0.04 m
%                               R1      =   0.15 m
%                               A       =   pi/4*0.04^2 m^2
%                               I0      =   8000*pi/2*R0^4*0.01+3*(1/3)*2000*0.1*0.02*.75^3 kg.m^2
%                               I1      =   8000*pi/2*R1^4*0.01 kg.m^2
%                               om      =   50 rad/s
%                               Tm      =   0.6 N.m
%                               rgh     =   1000*9.81*1 J/m^3
%                               the0    =   pi/2 rad
%                               ome0    =   3 rad/s
%
% inpI = [.03;.12;.04;.15;pi/4*.04^2;8000*pi/2*.04^4*.01+2000*.1*.02*.75^3;8000*pi/2*.15^4*.01;50;.6;1000*9.81] ;
% [tiarr,tharr,vearr,vfr,omef,freq,volp] = windpump(inpI,pi/2,3,60); % base case
%

% 	other internal variables
%
%           -   span	=   dim(-) time span
%           -   s       =   dim(-) time
%           -	Y		=	column vector	Y(1,1)	=	crank angle (rad)
%											Y(2,1)	=   dim(-) angular velocity
%           -	Y0		=	column vector	Y0(1,1)	=	initial crank angle (rad)
%											Y0(2,1)	=   initial dim(-) angular velocity
%                                           Y0(3,1) =   initial dim(-) volume pumped
%           -   dis     =   array of piston displacements/2/crankshaft radius, with
%                            origin at top dead center (TDC) (not used here)
%           -   vel     =   array of piston velocity/2/crankshaft radius/crankshaft angular velocity
%           -   acc     =   array of piston acceleration/2/crankshaft radius/ang.vel.^2 (not used here)
%           -   extr    =   four angles for TDC,BDC,TDC,BDC,TDC for 4-stroke cycle (not used here)
%           -   dim     =   size of the solution array
%           -   save0   =   array index when crank angle first passed through zero
%           -   savei   =   array index when crank angle last passed through zero
%           -   numcyl  =   number of whole cycles
%           -   maxdt   =   max time step for ODE45 integration
%           -   tf      =   dim(-) final time
%           -   tmodel  =   dim(-) time array for plotting analytical result
%           -   ee      =   intermediate variable
%           -   angmod  =   angular velocity predicted by analytical solution (rad/s)
%           -   a0      =   dim(-) initial angular velocity
%

global  lr Tmd loa reduc

% translate inpI into parameters of the system:

rcrank  =   inpI(1) ;
lcrank  =   inpI(2) ;
R0      =   inpI(3) ;
R1      =   inpI(4) ;
A       =   inpI(5) ;
I0      =   inpI(6) ;
I1      =   inpI(7) ;
om      =   inpI(8) ;
Tm      =   inpI(9) ;
rgh     =   inpI(10);

% error messages for wrong inputs:

if rcrank < 0 || rcrank == 0
    error('crankshaft radius <= 0')
end

if lcrank < 0 || lcrank == 0
    error('connecting rod length <= 0')
end

if R0 < 0 || R0 == 0
    error('small sprocket radius <= 0')
end

if R1 < 0 || R1 == 0
    error('small sprocket radius <= 0')
end

if A < 0 || A == 0
    error('piston cross-section area <= 0')
end

if I0 < 0 || I0 == 0
    error('windmill+small strocket MOI <= 0')
end

if I1 < 0 || I1 == 0
    error('large strocket MOI <= 0')
end

if om < 0 || om == 0
    error('characteristic windmill ang. vel. <= 0')
end

if Tm < 0 || Tm == 0
    error('characteristic windmill torque <= 0')
end

if rgh < 0 || rgh == 0
    error('water load <= 0')
end

% calculate the dimensionless numbers of the problem and pass as global variables
% for use in functions that cannot pass these
% variables as arguments; if necessary, use different names for the corresponding inputs
% to avoid confusing Matlab; convention: the postscript I denotes an input
% variable; the same variable without the postscript is the corresponding 
% global variable:

lr      =   lcrank/rcrank                       ;   % ratio of connecting rod length to crankshaft radius
reduc   =   R1/R0                               ;   % ratio of large/small sprocket radii
Tmd     =   reduc^2*Tm/om^2/(I1+I0*reduc^2)     ;   % dim(-) torque parameter
loa     =   rgh*A*rcrank/om^2/(I1+I0*reduc^2)   ;   % dim(-) load parameter
tf      =   tfI*om                              ;   % dim(-) final time

% initialize variables for the ODE solver:

span    =   [0 tf]                  ; % dim(-) time span
Y0(1,1)	=   the0                    ; % initial crank angle (rad)
Y0(2,1) =   ome0/om                 ; % initial dim(-) angular velocity
Y0(3,1) =   0                       ; % initial volume pumped

maxdt   =   pi*reduc/4          ; % max dim(-) time integration step to avoid skipping periods = (1/4) period at ang.vel.= 2*om
options =   odeset('MaxStep',maxdt)             ; % specify max time step through option call for ODE45
[s,Y]   =	ode45(@deriv,span,Y0,options)       ; % calls 4th-order Runge-Kutta ODE solver using function deriv below
%[s,Y]   =	ode45(@deriv,span,Y0)               ; % without options

%sol   =	ode45(@deriv,span,Y0)               ;
%xint	=   linspace(0,tf)                      ; example of use of deval (not used here)
%[sxint, spxint]   =   deval(sol,xint)          ;

dim     =   size(Y,1)               ; % finds the size of the solution array

save0   =   0                       ; % initialize the first index when crank angle passes through 0
numcyl  =   0                       ; % initialize the number of whole cycles

tiarr   =   zeros(dim,1)            ; % prepare the array tiarr for greater execution speed
tharr   =   zeros(dim,1)            ; % prepare the array tharr
vearr   =   zeros(dim,1)            ; % prepare the array vearr

for i = 1:dim
    tiarr(i,1)  =   s(i)/om             ; % build time array (s)
    tharr(i,1)  =   mod(Y(i,1),2*pi)    ; % build crank angle array (modulo 2pi)
    
    if i ~= 1
        if tharr(i,1)-tharr(i-1,1) < 0
            savei   =	i	;   % save the last index where modulo crank angle passed through zero
            if save0 == 0
            save0   =   i   ;   % save the first index where modulo crank angle passes through zero
            else
            numcyl  =   numcyl + 1	;   % number of whole cycles
            end
        end
    end
    
    vearr(i,1)  =   Y(i,2)*om     	;   % build crank angular velocity array (rad/s)
end

if numcyl ~= 0
omef	=   vearr(savei,1)        	;   % last angular velocity (rad/s) when crank angle was zero
else
omef    =   0                       ;	% if system did not complete a whole cycle
end

if numcyl ~= 0
freq    =   numcyl/(tiarr(savei,1)-tiarr(save0,1))      ;   % calculate the average cycle frequency (Hz)
else
freq    =   0                                           ;   % if system did not complete a whole cycle
end

% Calculate the array of instantaneous water volume flow rates

[dis,vel,acc,extr]  =   pistonkin(tharr,lr,0,0)     	;	% calls pistonkin without a plot or excentricity

vfr     =   -2.*rcrank.*A.*(tharr>pi).*vel.*vearr.*1000 ;   % calculate the instantaneous volume flow rate array (liter/s)

volp    =   2*A*rcrank*Y(dim,3)*1000                  	;   % calculate the total volume pumped (liter)

% calculate the analytical solution for comparison:

tmodel  =   linspace(0,tf)  ;   % time array in the model (-)

ee      =   exp(-2.*tmodel/pi*sqrt(pi^2*Tmd^2-pi*reduc*Tmd*loa))    ;   % intermediate variable: exponential of time
a0      =   ome0/om                                                 ;   % dim(-) initial angular velocity

angmod  =   ((-1 + ee).*loa + ...
   a0.*((-(-1 + ee)).*pi.*tmodel + (1 + ee).*sqrt((-loa).*...
   pi.*reduc.*tmodel + pi^2.*tmodel.^2)))./ ...
   ((-(-1 + ee)).*pi*(-1 + a0.*reduc).*tmodel + (1 + ee).* ...
    sqrt((-loa).*pi.*reduc.*tmodel + pi^2.*tmodel.^2))        	;   % predicted dim(-) mean angular velocity
    
%angmod  =   angmod.*om                                              ;   % convert to dimensional value (rad/s)
    
scrsz   = get(0,'ScreenSize')                           ;   % finds the screen size of the current machine
figure('Position',...
    [1 scrsz(4) scrsz(3) scrsz(4)/2])                   ;	% reserves space on the screen for the figure

subplot(1,3,1)                                              % plot crank angular velocity vs. time
plot(tiarr,vearr,'b',...
    [tiarr(1,1);tiarr(dim,1)],[freq*2*pi;freq*2*pi],'r',...
    tmodel./om,angmod.*om,'g')
xlabel('time (s)')
ylabel('crank angular velocity (rad/s)')

subplot(1,3,2)                                              % plot volume flow rate vs. time
plot(tiarr,vfr,'b')
xlabel('time (s)')
ylabel('volume flow rate (liter/s)')

subplot(1,3,3)                                              % subplot for displaying parameters
axis off

text(0.1,0,['initial crank angle = ',num2str(the0)])
text(0.1,.05,['initial angular velocity (rad/s) = ',num2str(ome0)])
text(0.1,.1,['crankshaft radius (m) = ',num2str(rcrank)])
text(0.1,.15,['connecting rod length (m) = ',num2str(lcrank)])
text(0.1,.2,['small sprocket radius (m) = ',num2str(R0)])
text(0.1,.25,['large sprocket radius (m) = ',num2str(R1)])
text(0.1,.3,['piston cross-section (m^2) = ',num2str(A)])
text(0.1,.35,['windmill & sprocket MOI (kg m^2) = ',num2str(I0)])
text(0.1,.4,['large sprocket MOI (kg/m^2) = ',num2str(I1)])
text(0.1,.45,['windmill angular vel at peak torque (rad/s) = ',num2str(om)])
text(0.1,.5,['windmill peak torque (N.m) = ',num2str(Tm)])
text(0.1,.55,['water potential energy/u. vol. (J/m^3) = ',num2str(rgh)])
text(0.1,.6,['torque parameter (-) = ',num2str(Tmd)])
text(0.1,.65,['load parameter (-) = ',num2str(loa)])
text(0.1,.7,['last angular velocity at zero crank angle (rad/s) = ',num2str(omef)])
text(0.1,0.75,['average cycle frequency (Hz) = ',num2str(freq)])
text(0.1,0.8,['average angular velocity (rad/s) = ',num2str(freq*2*pi)])
text(0.1,0.85,'red curve: average angular velocity (rad/s)')
text(0.1,.9,['total volume pumped (liter) = ',num2str(volp)])


%---------------------

function Yprime = deriv(s,Y)

% Calculates the derivative of the column-vector Y 
% with respect to s. Implements windpump dynamics.
%
% MYL 2/10/2010.
%
%
%	input	-	Y		=	column vector	Y(1,1)	=	crank angle (rad)
%											Y(2,1)	=   dim(-) angular velocity
%                                           Y(3,1)  =   dim(-) volume pumped
%			-	s		=	dim(-) time
%
%	output	-	Yprime	=	column vector	d(crank_angle)/ds
%											d(ang_vel)/ds
%                                           dim(-) vol. flow rate
%
%   global variables    postscript I denotes an input variable destined to
%                       become global (avoids confusing Matlab)
%   
%           -   lr      =   ratio of connecting rod length to crankshaft radius
%           -   Tmd     =   dim(-) torque parameter
%           -   loa     =   dim(-) load parameter
%           -   reduc   =   reduction ratio R1/R0
%

global  lr Tmd loa reduc

% 	other internal variables
%
%           -   modthe	=   crank angle modulo 2 pi

%           -   pres    =   net dim(-) force on piston
%           -   dis     =   piston displacement/2/crankshaft radius, with
%                                   origin at top dead center (TDC)
%           -   vel     =   piston velocity/2/crankshaft radius/crankshaft angular velocity
%           -   acc     =   piston acceleration/2/crankshaft radius/ang.vel.^2
%           -   extr    =   four angles for TDC,BDC,TDC,BDC,TDC for 4-stroke cycle (not used here)
%           -   thec    =   angle of the connecting rod to cylinder axis (rad)
%

modthe              =   mod(Y(1,1),2*pi)        ;   % find the crank angle modulo 2pi

% Calculate the net dimensionless force exerted on the piston
% of the windpump; in this case, the piston experiences no force during
% intake for theta = [0,pi] and the full water load during exhaust for
% theta = [pi,2pi].
%

if modthe < pi
    pres    =	0       ; % net intake pressure is zero
else
    pres    =   -1       ; % exhaust dim(-) force
end

thec    =   asin(sin(modthe)/lr)                            ;

[dis,vel,acc,extr]  = pistonkin(modthe,lr,0,0)  ;	% calls pistonkin without plot or excentricity

Yprime(1,1) = Y(2,1)                                        ;
Yprime(2,1) = -loa*pres*sin(modthe+thec)/cos(thec) + ...
                Tmd*Y(2,1)*(2-reduc*Y(2,1))                 ;
Yprime(3,1) = -vel*Y(2,1)*(modthe>pi)                       ;

%---------------------

function [dis,vel,acc,extr] = pistonkin(theta,ll,dd,flagplot)

% Calculates and plots (if needed) the kinematics of a 
% piston and crankshaft with possible excentricity. 
% Reference: crank_mech.nb
%
% Michel Louge 4/3/07
%
%   input   -   theta	=   array of crank angles
%           -   ll      =   ratio of connecting rod length/crankshaft radius
%           -   dd      =   ratio of excentricity/crankshaft radius
%           -   flagplot=   flag: if = 1, show a plot; if other, then don't
%
%   output  -   dis     =   array of piston displacement/2/crankshaft radius, with
%                                   origin at top dead center (TDC)
%           -   vel     =   array of piston velocity/2/crankshaft radius/crankshaft angular velocity
%           -   acc     =   array of piston acceleration/2/crankshaft radius/ang.vel.^2
%           -   extr    =   array of four angles for TDC,BDC,TDC,BDC,TDC for four-stroke cycle
%
%   global variables
%   
%           -   none    =   none 
%
%   to access the proper directory, type
%   cd ~/Michel_files/Classwork/MAE_225/Spring_2008/projects/AC/Matlab_dynamics
%
% Example:
%
% clear all; [dis,vel,acc,extr] = pistonkin(linspace(0,4*pi),3,0.5,1);

% global  none

% 	other internal variables
%
%           -   scrsz   =   screen size for plots
%           -   cost    =   cos[crank angle at TDC]
%           -   cosb   	=   - cosine of crank angle at BDC

cost	=   sqrt(1-dd^2/(ll+1)^2)                           ;
cosb	=   sqrt(1-dd^2/(ll-1)^2)                           ;

extr    =   [acos(cost) pi+acos(cosb) 2*pi+acos(cost) 3*pi+acos(cosb)]  ;

dis     =   (1/2)*((ll + 1)*cost - cos(theta) - ...
                sqrt(ll^2 - (sin(theta) - dd).^2))          ;
            
vel     =   (1/2)*(sin(theta) + (cos(theta).* ...
            (-dd + sin(theta)))./ ...
            sqrt(ll^2 - (dd - sin(theta)).^2))              ;
acc     =   (1./(8*(ll^2 - (dd - sin(theta)).^2).^(3/2))).* ...
            (4*ll^2*cos(theta).^2 - ...
            2*sin(theta).*(-dd + sin(theta)).* ...
            (-1 - 2*dd^2 + 2*ll^2 + cos(2*theta) + ...
            4*dd*sin(theta)) + ...
            sqrt(2)*cos(theta).* ...
            (-1 - 2*dd^2 + 2*ll^2 + cos(2*theta) + ...
            4*dd*sin(theta)).^(3/2))                        ;
        
if flagplot == 1
    
scrsz   = get(0,'ScreenSize')                           ;   % finds the screen size of the current machine
figure('Position',...
    [1 scrsz(4) scrsz(3) scrsz(4)/2])                   ;	% reserves space on the screen for the figure

subplot(1,4,1)                                              % plot displacement vs. crank angle
plot(theta,dis)
xlabel('crank angle (rad)')
ylabel('displacement/crankshaft diameter (-)')
axis([min(theta) max(theta)*1.1 min(dis) 1.1*max(dis)])

subplot(1,4,2)                                              % plot velocity vs. crank angle
plot(theta,vel)
xlabel('crank angle (rad)')
ylabel('velocity/crankshaft diameter/angular velocity (-)')
axis([0 max(theta)*1.1 1.1*min(vel) 1.1*max(vel)])

subplot(1,4,3)                                              % plot acceleration vs. crank angle
plot(theta,acc)
xlabel('crank angle (rad)')
ylabel('acceleration/crankshaft diameter/angular velocity^2 (-)')
axis([0 max(theta)*1.1 1.1*min(acc) 1.1*max(acc)])

subplot(1,4,4)                                              % subplot for displaying parameters
axis off

text(0.1,0,['connecting rod length/crankshaft radius = ',num2str(ll)])
text(0.1,.1,['excentricity/crankshaft radius = ',num2str(dd)])
text(0.1,.2,['TDC at intake (deg)= ',num2str(extr(1)*180/pi)])
text(0.1,.3,['BDC at compression (deg)= ',num2str(extr(2)*180/pi)])
text(0.1,.4,['TDC at expansion (deg)= ',num2str(extr(3)*180/pi)])
text(0.1,.5,['TDC at exhaust (deg)= ',num2str(extr(4)*180/pi)])

end
