function [var,dvar] = simple
% simple analysis - including heat transfer and pressure drop effects
% Israel Urieli, 7/22/2002 (modified 12/3/2003 for temp plots)
% Returned values: 
%   var(22,37) array of variable values every 10 degrees (0 - 360)
%   dvar(16,37) array of derivatives every 10 degrees (0 - 360)

% Row indices of the var, dvar arrays:
 TC = 1;  % Compression space temperature [K]
 TE = 2;  % Expansion space temperature [K]
 QK = 3;  % Heat transferred to the cooler [J]
 QR = 4;  % Heat transferred to the regenerator [J]
 QH = 5;  % Heat transferred to the heater [J]
 WC = 6;  % Work done by the compression space [J]
 WE = 7;  % Work done by the expansion space [J]
 W  = 8;  % Total work done (WC + WE) [J]
 P  = 9;  % Pressure [Pa]
 VC = 10; % Compression space volume [m^3]
 VE = 11; % Expansion space volume [m^3]
 MC = 12; % Mass of gas in the compression space [kg]
 MK = 13; % Mass of gas in the cooler [kg]
 MR = 14; % Mass of gas in the regenerator [kg]
 MH = 15; % Mass of gas in the heater [kg]
 ME = 16; % Mass of gas in the expansion space [kg]
 TCK = 17; % Conditional temperature compression space / cooler [K]
 THE = 18; % Conditional temeprature heater / expansion space [K]
 GACK = 19; % Conditional mass flow compression space / cooler [kg/rad]
 GAKR = 20; % Conditional mass flow cooler / regenerator [kg/rad]
 GARH = 21; % Conditional mass flow regenerator / heater [kg/rad]
 GAHE = 22; % Conditional mass flow heater / expansion space [kg/rad]
% Size of var(ROWV,COL), dvar(ROWD,COL)
 ROWV = 22; % number of rows in the var matrix
 ROWD = 16; % number of rows in the dvar matrix
 COL = 37; % number of columns in the matrices (every 10 degrees) 
%======================================================================

global freq % cycle frequency [herz]
global tk tr th % cooler, regenerator, heater temperatures [K]
global cqwr % regenerator housing thermal conductance [W/K]

twk = tk; % Cooler wall temp - equal to initial cooler gas temp
twh = th; % Heater wall temp - equal to initial heater gas temp
epsilon = 1; % allowable temperature error bound for cyclic convergence
terror = 10*epsilon; % Initial temperature error (to enter loop)

while (terror>epsilon)
   [var,dvar] = adiab;
   tgh = hotsim(var,twh); % new heater gas temperature
   tgk = kolsim(var,twk); % new cooler gas temperature
   terror = abs(th - tgh) + abs(tk - tgk);
   th = tgh;
   tk = tgk;
   tr = (th-tk)/log(th/tk);
end

fprintf('===== converged heater and cooler mean temperatures =====\n');
fprintf('heater wall/gas temperatures: Twh = %.1f[K], Th = %.1f[K]\n',twh,th);
fprintf('cooler wall/gas temperatures: Twk = %.1f[K], Tk = %.1f[K]\n',twk,tk);
% Print out ideal adiabatic analysis results
eff = var(W,COL)/var(QH,COL); % engine thermal efficency
Qkpower = var(QK,COL)*freq;   % Heat transferred to the cooler (W)
Qrpower = var(QR,COL)*freq;   % Heat transferred to the regenerator (W)
Qhpower = var(QH,COL)*freq;   % Heat transferred to the heater (W)
Wpower = var(W,COL)*freq;     % Total power output (W)
fprintf('========== ideal adiabatic analysis results ==========\n');
fprintf(' Heat transferred to the cooler: %.2f[W]\n', Qkpower);
fprintf(' Net heat transferred to the regenerator: %.2f[W]\n', Qrpower);
fprintf(' Heat transferred to the heater: %.2f[W]\n', Qhpower);
fprintf(' Total power output: %.2f[W]\n', Wpower);
fprintf(' Thermal efficiency: %.1f[%%]\n', eff*100);
fprintf('======================================================\n');
% Temperature plot of the simple simulation results
		figure
		x = 0:10:360;
		Tcomp = var(TC,:);
		Texp = var(TE,:);
		plot(x,Tcomp,'b-',x,Texp,'r-');
		hold on
		x = [0,360];
		y = [twk,twk];
		plot(x,y,'b--')
		y = [tk,tk];
		plot(x,y,'b-')
		y = [tr,tr];
		plot(x,y,'g-')
		y = [th,th];
		plot(x,y,'r-')
		y = [twh,twh];
		plot(x,y,'r--')
		hold off
		grid on
		xlabel('Crank angle (degrees)');
		ylabel('Temperature (K)');
		title('Simple Simulation - Wall and Gas Temps vs crank angle');
% Various plots of the ideal adiabatic simulation results
plotadiab(var,dvar);

fprintf('============= Regenerator simple analysis ============\n');
qrloss = regsim(var);
fprintf(' Regenerator net enthalpy loss: %.1f[W]\n', qrloss*freq);
qwrl = cqwr*(twh - twk)/freq;
fprintf(' Regenerator wall heat leakage: %.1f[W]\n', qwrl*freq);

fprintf('========= pressure drop simple analysis =============\n');
dwork = worksim(var,dvar);
fprintf(' Pressure drop available work loss: %.1f[W]\n', dwork*freq)
actWpower = Wpower - dwork*freq;
actQhpower = Qhpower + qrloss*freq + qwrl*freq;
acteff = actWpower/actQhpower;
fprintf(' Actual power from simple analysis: %.1f[W]\n', actWpower); 
fprintf(' Actual heat power in from simple analysis: %.1f[W]\n', actQhpower); 
fprintf(' Actual efficiency from simple analysis: %.1f[%%]\n', acteff*100);