function [tiarr,tharr,vearr,omef,freq] = aircart(the0,ome0,tf,gamI,the2I,pregI,v1I,lossI,loaI,inerI,lrI)

% Calculates and plots the dynamics of an aircart entrained through a
% crankshaft of finite inertia by a single-cylinder, single-acting 
% adiabatic air motor without friction connected via a lossy line to 
% a regulated isothermal tank.
%
% Dimensionless variables:
%
%   time        made dim(-) with    sqrt(2R m/A pa)
%   pressure    made dim(-) with    pa
%   cyl. volume made dim(-) with    V at bottom dead center (BDC)
%
% Michel Louge 3/23/08
%
%   input   -   the0	=   initial crank angle (radians); = 0 at top dead center (TDC)
%           -   ome0    =   initial dim(-) angular velocity
%           -   tf      =   dim(-) final time
%           -   gamI    =   cp/cv of air (=7/5)
%           -   the2I   =   crank angle (radians) at end of intake substroke
%           -   pregI   =   dim(-) regulated pressure
%           -   v1I     =   V at TDC/displacement volume
%           -   lossI   =   line pressure loss parameter = f(rho/2)(A/a)^2(Vdis/m)(l/d)
%           -   loaI    =   load parameter = [2R r m M g sin(alp)]/[I pa A]
%           -   inerI   =   inertia parameter = 2mR^2/I
%           -   lrI     =   ratio of connecting rod length/crankshaft radius
%
%                           variable                                example below
%
%               where       pa  = ambient pressure                  1e5	    Pa
%                           Ta  = ambient temperature               290     ¡K
%                           pr	= regulator pressure                3e5     Pa
%                           RR  = fundamental gas cst               8.314   J/mole¡K
%                           MW  = air molecular weight              0.029   kg/mole
%                         	rho = air density at regulated pres     3.6     kg/m^3
%                          	m   = piston mass                       0.02    kg
%                        	R   = radius of crankshaft              0.05	m
%                           L   = length of crankshaft              0.1     m
%                          	f   = Moody friction parameter          0.04
%                           bore= piston bore diameter              0.03    m
%                           hei = piston height                     0.01    m
%                       	A   = piston cross-sectional area       7e-4    m^2
%                          	d   = line diameter                     3e-3	m
%                          	a   = line cross-section area           7e-6    m^2
%                          	l   = line length                       .1      m
%                          	Vdis= piston displacement volume       	7e-5    m^3
%                          	r   = transmission effective radius     2e-3    m
%                                   ie, Æx=r*Æthe
%                           Æx  = travel of air cart over ...
%                           Æthe= ... difference in crank angle
%                           M   = total mass of air cart            3       kg
%                           g   = gravitational acceleration        9.81    m/s^2
%                           alp = angle of inclination              10¡
%                           I   = crankshaft moment of inertia      3e-4    kg.m^2
%
%                                                       gam         7/5
%                                                       preg        3
%                                                       v1          1/4
%                                                       loss        86
%                                                       loa         0.001
%                                                       iner        0.36
%                                                       lr          2
%
%   output  -   tiarr 	=   array of dim(-) times
%           -   tharr 	=   array of crank angles modulo 2pi
%           -   vearr 	=   array of crank angular velocities
%           -   omef    =   last dim(-) angular velocity when crank angle was zero
%                                  (zero if system did not complete a whole cycle)
%           -   freq    =   average dim(-) cycle frequency (zero if system
%                                           did not complete a whole cycle)
%
%   global variables    postscript I denotes an input variable destined to
%                       become global (avoids confusing Matlab)
%   
%           -   gam     =   gamI    see above
%           -   the2    =   the2I
%           -   preg    =   pregI
%           -   v1      =   v1I
%           -   loss    =   lossI
%           -   loa     =   loaI
%           -   iner    =   inerI
%           -   lr      =   lrI
%           -   oldp    =   last pressure within the intake substroke
%           -   dis2    =   piston displacement/2/crankshaft radius at crank angle the2
%
%   to access the proper directory, type
%   cd ~/Michel_files/Classwork/MAE_225/Spring_2008/projects/AC/Matlab_dynamics
%
% Examples:
%
% [tiarr,tharr,vearr,omef,freq] = aircart(0,.01,300,7/5,pi/2,3,1/4,86,.001,.36,2); % base case
% [tiarr,tharr,vearr,omef,freq] = aircart(0,.01,300,7/5,pi/2,3,1/4,86,.001,.1,2); % larger flywheel
% [tiarr,tharr,vearr,omef,freq] = aircart(0,.01,300,7/5,pi/2,3,1/4,100,.001,.36,2); % larger line loss
% [tiarr,tharr,vearr,omef,freq] = aircart(0,.01,300,7/5,pi/2,4,1/4,86,.001,.36,2); % larger regulator pressure
% [tiarr,tharr,vearr,omef,freq] = aircart(0,.01,500,7/5,pi/2,3,.15,86,.001,.36,2); % larger displacement volume
% [tiarr,tharr,vearr,omef,freq] = aircart(0,.01,300,7/5,pi/2,3,1/4,86,.002,.36,2); % larger load
%

global  gam the2 preg v1 loss loa iner lr oldp dis2

% initialize global variables for use in functions that cannot pass these
% variables as arguments; 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:

gam     =   gamI    ;
the2    =   the2I   ;
preg    =   pregI   ;
v1      =   v1I     ;
loss    =   lossI   ;
loa     =   loaI    ;
iner    =   inerI   ;
lr      =   lrI     ;

% 	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
%           -   vel     =   piston vel/2/crankshaft radius/crankshaft angular vel (not used here)
%           -   acc     =   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
%

% error messages for wrong inputs:

if the0 < 0
    error(['initial crank angle < 0'])
    return
end

if the0 > the2 | the0 == the2
    error(['initial crank angle too large'])
    return
end

if ome0 < 0 | ome0 == 0
    error(['initial angular velocity < 0 or zero'])
    return
end

if tf < 0
    error(['time span < 0'])
    return
end

if gam > 5/3 | gam < 1 | gam == 1
    error(['cp/cv unphysical'])
    return
end

if the2 < 0 | the2 > pi
    error(['crank angle at end of intake out of bounds'])
    return
end

if preg < 1 | preg == 1
    error(['low regulator pressure'])
    return
end

if v1 < 0
    error(['Vtdc/Vbdc < 0'])
    return
end

if loss < 0
    error(['loss parameter < 0'])
    return
end

if loa < 0
    error(['load parameter < 0'])
    return
end

if iner < 0
    error(['inertia parameter < 0'])
    return
end

[dis2,vel,acc,extr] = pistonkin(the2,lr,0,0)	; % calculate the piston clearance at angle the2

% initialize variables for the ODE solver:

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

[s,Y]   =	ode45(@deriv,span,Y0)	; % calls 4th-order Runge-Kutta ODE solver using function deriv below

%sol   =	ode45(@deriv,span,Y0)               ;
%xint	=   linspace(0,tf)                      ; example of use of deval
%[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

for i = 1:dim
    tiarr(i,1)  =   s(i)                ; % build time array
    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)              ; % build crank angular velocity array
end

if numcyl ~= 0
omef	=   vearr(savei,1)        	; % last dim(-) angular velocity 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 dim(-) cycle frequency
else
freq    =   0                                           ;   % if system did not complete a whole cycle
end

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,2,1)                                              % plot crank angular velocity vs. time
plot(tiarr,vearr,'b',...
    [tiarr(1,1);tiarr(dim,1)],[freq*2*pi;freq*2*pi],'r')
xlabel('time (-)')
ylabel('crank angular velocity (-)')

%subplot(1,3,2)                                              % plot crank angular velocity vs. time
%plot(tiarr,tharr)
%xlabel('time (-)')
%ylabel('crank angle mod 2pi (rad)')
%axis([min(tiarr) max(tiarr)*1.1 min(tharr) 1.1*max(tharr)])

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

text(0.1,0,['initial crank angle = ',num2str(the0)])
text(0.1,.05,['initial angular velocity (-) = ',num2str(ome0)])
text(0.1,.1,['crank angle at end of intake (rad) = ',num2str(the2)])
text(0.1,.15,['Vtdc/displacement volume = ',num2str(v1)])

text(0.1,.25,['cp/cv = ',num2str(gam)])
text(0.1,.3,['regular pressure/ambient pressure = ',num2str(preg)])
text(0.1,.35,['line pressure loss parameter (-) = f(rho/2)(A/a)^2(Vdis/m)(l/d) = ',num2str(loss)])

text(0.1,.45,['load parameter (-) = [2R r m M g sin(alp)]/[I pa A] = ',num2str(loa)])
text(0.1,.5,['inertia parameter (-) = 2mR^2/I = ',num2str(iner)])

text(0.1,.55,['last angular velocity at zero crank angle  = ',num2str(omef),...
                  ' sqrt[piston area*ambient pressure/stroke/piston mass]'])

text(0.1,0.65,['average cycle frequency  = ',num2str(freq),...
                  ' sqrt[piston area*ambient pressure/stroke/piston mass]'])
text(0.1,0.7,['average angular velocity  = ',num2str(freq*2*pi),...
                  ' sqrt[piston area*ambient pressure/stroke/piston mass]'])

text(0.1,0.8,['abscissa: time made dim(-) with characteristic time = sqrt[stroke*piston mass/piston area/ambient pressure]'])
text(0.1,0.85,['red curve: average angular velocity (-)'])
text(0.1,0.9,['blue curve: instantaneous angular velocity (-)'])
text(0.1,0.95,['legend'])


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

function Yprime = deriv(s,Y)

% Calculates the derivative of the column-vector Y 
% with respect to s. Implements air cart dynamics.
%
% MYL 3/23/08.
%
%
%	input	-	Y		=	column vector	Y(1,1)	=	crank angle (rad)
%											Y(2,1)	=   dim(-) angular velocity
%			-	s		=	dim(-) time
%
%	output	-	Yprime	=	column vector	d(crank_angle)/ds
%											d(ang_vel)/ds
%
%   global variables
%   
%           -   gam     =   cp/cv of air (=7/5)
%           -   the2    =   crank angle (radians) at end of intake substroke
%           -   preg    =   dim(-) regulated pressure
%           -   v1      =   V at TDC/displacement volume
%           -   loss    =   line pressure loss parameter = f(rho/2)(A/a)^2(Vbdc/m)(l/d)
%           -   loa     =   load parameter = [2R r m M g sin(alp)]/[I pa A]
%           -   iner    =   inertia parameter = 2mR^2/I
%           -   lr      =   ratio of connecting rod length/crankshaft radius
%           -   oldp    =   last pressure within the intake substroke
%                           (saved to calculate isentropic expansion)
%

global  gam the2 preg v1 loss loa iner lr oldp dis2

% 	other internal variables
%
%           -   modthe	=   crank angle modulo 2 pi
%           -   pres    =   net dim(-) pressure 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)
%           -   interm  =   intermediate variable used to speed up computation
%

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

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

% Calculate the net dimensionless pressure exerted on the piston
% of the air motor connected to a regulated tank and exposed on
% its back side to ambient pressure:
%

if modthe < the2
    pres    =	preg - loss*(Y(2,1)*vel)^2	; % intake pressure reduced by line losses
    oldp    =   pres                        ; % saves pressure at expansion's start 
else
    if modthe < pi
        pres	=	oldp*((v1+dis2)/(v1+dis))^gam	; % expansion substroke
    else
        pres    =   1                               ; % exhaust stroke to ambient pressure
    end
end

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

interm  =   iner*sin(modthe+thec)/cos(thec)         ; % calculated once to speed up computation

Yprime(1,1) = Y(2,1)                                ;
Yprime(2,1) = ( interm*(pres-1-acc*Y(2,1)^2)-loa)/...
                (1+interm*vel)                      ;

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

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
