function dwork = worksim(var,dvar);
% Evaluate the pressure drop available work loss [J]
% Israel Urieli, 7/23/2002
% Arguments:
%   var(22,37) array of variable values every 10 degrees (0 - 360)
%   dvar(16,37) array of derivatives every 10 degrees (0 - 360)
% Returned value: 
%   dwork - pressure drop available work loss [J]

% 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 tk tr th % cooler, regenerator, heater temperatures [K]
global freq omega % cycle frequency [herz], [rads/s]
global vh % heater void volume [m^3]
global ah % heater internal free flow area [m^2]
global dh % heater hydraulic diameter [m]
global lh % heater effective length [m]
global vk % cooler void volume [m^3]
global ak % cooler internal free flow area [m^2]
global dk % cooler hydraulic diameter [m]
global lk % cooler effective length [m]
global vr % regen void volume [m^3]
global ar % regen internal free flow area [m^2]
global lr % regenerator effective length [m]
global dr % regen hydraulic diameter [m]
global matrix_type % m)esh or f)oil

dtheta = 2*pi/36;
dwork = 0; % initialise pumping work loss

for(i = 1:1:36)
    gk = (var(GACK,i) + var(GAKR,i))*omega/(2*ak);
    [mu,kgas,re(i)] = reynum(tk,gk,dk);
    [ht,fr] = pipefr(dk,mu,re(i));
    dpkol(i) = 2*fr*mu*vk*gk*lk/(var(MK,i)*dk^2);

    gr = (var(GAKR,i) + var(GARH,i))*omega/(2*ar);
    [mu,kgas,re(i)] = reynum(tr,gr,dr);
    if(strncmp(matrix_type,'m',1))
         [st,fr] = matrixfr(re(i));
    elseif (strncmp(matrix_type,'f',1))
         [st,ht,fr] = foilfr(dr,mu,re(i));
    end
    dpreg(i) = 2*fr*mu*vr*gr*lr/(var(MR,i)*dr^2);

    gh = (var(GARH,i) + var(GAHE,i))*omega/(2*ah);
    [mu,kgas,re(i)] = reynum(th,gh,dh);

    [ht,fr] = pipefr(dh,mu,re(i));
    dphot(i) = 2*fr*mu*vh*gh*lh./(var(MH,i)*dh^2);
    dp(i) = dpkol(i) + dpreg(i) + dphot(i);

    dwork=dwork+dtheta*dp(i)*dvar(VE,i); % pumping work [J]
	pcom(i) = var(P,i);
	pexp(i) = pcom(i) + dp(i);
end

dpkol(COL) = dpkol(1);
dpreg(COL) = dpreg(1);
dphot(COL) = dphot(1);
dp(COL) = dp(1);
pcom(COL) = pcom(1);
pexp(COL) = pexp(1);

choice = 'x';
while(~strncmp(choice,'q',1))
	fprintf('Choose pumping loss plot type:\n');
	fprintf('   h - for heat exchanger pressure drop plot\n');
	fprintf('   p - for working space pressure plot\n');
	fprintf('   q - to quit\n');
	choice = input('h)x_pdrop, p)ressure, q)uit: ','s');
	if(strncmp(choice,'h',1))
		figure;
		x = 0:10:360;
		plot(x,dpkol,'b-',x,dphot,'r-',x,dpreg,'g-');
		grid on
		xlabel('Crank angle (degrees)');
		ylabel('Heat exchanger pressure drop [Pa]');
		title('Heat exchanger pressure drop vs crank angle');
	elseif(strncmp(choice,'p',1))
		figure
		x = 0:10:360;
		pcombar = pcom*1e-5;
		pexpbar = pexp*1e-5;
		plot(x,pcombar,'b-',x,pexpbar,'r-');
		grid on
		xlabel('Crank angle (degrees)');
		ylabel('Working space pressure [bar]');
		title('Working space pressure vs crank angle');
	end
end
fprintf('quitting pressure plots...\n');