function [var,dvar] = adiab
% ideal adiabatic model simulation
% Israel Urieli, 7/6/2002
% 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)
global tk th % cooler, heater temperatures [K]
% Row indices of the var, dvar matrices, and the y,dy variable vectors:
 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), y(ROWV), dvar(ROWD,COL), dy(ROWD)
 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) 
%===========================================================
 fprintf('============Ideal Adiabatic Analysis====================\n')
 fprintf('Cooler Tk = %.1f[K], Heater Th = %.1f[K]\n', tk, th);
 epsilon = 1;  % Allowable error in temperature (K)
 max_iteration = 20;  % Maximum number of iterations to convergence
 ninc = 360; % number if integration increments (every degree)
 step = ninc/36; % for saving values in var, dvar matrices
 dtheta = 2.0*pi/ninc; % integration increment (radians)
% Initial conditions:
 y(THE) = th;
 y(TCK) = tk;
 y(TE) = th;
 y(TC) = tk;
 iter = 0;
 terror = 10*epsilon; % Initial error to enter the loop
% Iteration loop to cyclic convergence
 while ((terror >= epsilon)&(iter < max_iteration))
       % cyclic initial conditions
       tc0 = y(TC);
       te0 = y(TE);
       theta = 0;
       y(QK) = 0;
       y(QR) = 0;
       y(QH) = 0;
       y(WC) = 0;
       y(WE) = 0;
       y(W) = 0;
       fprintf('iteration %d: Tc = %.1f[K], Te = %.1f[K]\n',iter,y(TC),y(TE))
       for(i = 1:1:ninc)
             [theta,y,dy] = rk4('dadiab',7,theta,dtheta,y);
       end
       terror = abs(tc0 - y(TC)) + abs(te0 - y(TE));
       iter = iter + 1;
 end
 if (iter >= max_iteration)
        fprintf('No convergence within %d iteration\n',max_iteration)
 end
 % Initial var and dvar matrix
 var = zeros(22,37);
 dvar = zeros(16,37);
 % a final cycle, to fill the var, dvar matrices
 theta=0;
 y(QK)=0;
 y(QR)=0;
 y(QH)=0;
 y(WC)=0;
 y(WE)=0;
 y(W)=0;
 [var,dvar] = filmatrix(1,y,dy,var,dvar);
 for(i = 2:1:COL)
        for(j = 1:1:step)
              [theta,y,dy] = rk4('dadiab',7,theta,dtheta,y);
        end
        [var,dvar] = filmatrix(i,y,dy,var,dvar);
 end

function [x, y, dy] = rk4(deriv,n,x,dx,y)
%Classical fourth order Runge-Kutta method
%Integrates n first order differential equations 
%dy(x,y) over interval x to x+dx
%Israel Urieli - Jan 21, 2002
x0 = x;
y0 = y;
[y,dy1] = feval(deriv,x0,y);
for i = 1:n
        y(i) = y0(i) + 0.5*dx*dy1(i);
end
xm = x0 + 0.5*dx;
[y,dy2] = feval(deriv,xm,y);
for i = 1:n
        y(i) = y0(i) + 0.5*dx*dy2(i);
end
[y,dy3] = feval(deriv,xm,y);
for i = 1:n
        y(i) = y0(i) + dx*dy3(i);
end
x = x0 + dx;
[y,dy] = feval(deriv,x,y);
for i = 1:n
        dy(i) = (dy1(i) + 2*(dy2(i) + dy3(i)) + dy(i))/6;
        y(i) = y0(i) + dx*dy(i);
end

function [var,dvar]=filmatrix(j,y,dy,var,dvar);
% Fill in the j-th column of the var, dvar matrices with values of y, dy
% Israel Urieli, 7/20/2002
% Arguments:  j - column index (1 - 37, every 10 degrees of cycle angle)
%             y(ROWV) - vector of current variable values 
%             dy(ROWD) vector of current derivatives
%             var(ROWV,37) - matrix of current variables vs cycle angle
%             dvar(ROWD,37) - matrix of current derivatives vs cycle angle
% Returned values: 
%             var(ROWV,37) - matrix of updated variables vs cycle angle
%             dvar(ROWD,37) - matrix of updated derivatives vs cycle angle
ROWV = 22; % number of rows in the var matrix
ROWD = 16; % number of rows in the dvar matrix
for (i = 1:1:ROWV)
   var(i,j) = y(i);
end
for (i = 1:1:ROWD)
   dvar(i,j) = dy(i);
end