function [dX] = mass(t,M)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
T = t;
m = M;

% alpha = 1.3e-6; beta = 5.95e-5;
% DO = 3.36; POi = 98;
%%% gas exchange parameters
alpha   = 1.3e-6;  % O2 solubility by Henry's law(M/mmHg)
beta    = 3.08e-5; % CO2 solubility by Hency's law (M/mmHg)
CHb     = 0.021;   % Hb binding site conc (mol/L of RBC's)
Hct     = 0.40;    % hematocrit (unitless)
C0      = CHb*Hct; % blood oxygen binding capacity (mol/L)
n       = 2.7;     % Hill exponent
P50     = 30;      % half-max saturation of Hb
DO      = 1;       % (apparent) O2 diffusion coefficient (ml/s)
Pair = 150; 
gamma   = 16800; %mmHg * L / mol (convert M to mmHg, for gas)


R_R = 1.17; Q_v = 2.97;
RFt = 0.0154; V0 = 3;

Q_air = Q_v*sin(2*pi*R_R*t);
qi = RFt*Q_air;
A = qi/(2*pi*R_R);
V_alv = A.*(1-cos(2*pi*R_R*t))+ V0;

Pairway1 = M(1);                %airway 1 - interacts with atmosphere
Pairway2 = M(2);                %airway 2 - big ol' mixing chamber
Malv     = M(3:end);            %mass of O2 in the alveoli (mol)
PalvO    = gamma*Malv./V_alv;  
Vairway = 0.6;
f = 0.7;
Vairway1 = Vairway*(1-f);
Vairway2 = Vairway*f;

P_alv = (beta*m)./V_alv;
% m = ((qi/beta).*P_alv - alpha*DO.*(P_alv-POi)).*(Q_air<0) ;
PAO2 = 100; %resting systemic artial O2 partial pressure (mmHg)
PVO2 = 45;  %resting systemic venous O2 partial pressure (mmHg)
CAO2 = alpha*PAO2 + C0*(PAO2.^n)./(PAO2.^n + P50.^n);   %converted to mol/L
CVO2 = alpha*PVO2 + C0*(PVO2.^n)./(PVO2.^n + P50.^n);   %coverted to mol/L
CO     = CVO2;%.*ones(10,1);   %initial oxygen concentration for capillaries (per capillary segment)
%%
Plookup = 0:1:10000; %look up table
Clookup = alpha*Plookup + C0*(Plookup.^n)./(Plookup.^n + P50.^n);
PO  = interp1(Clookup,Plookup,CO); %capillary oxygen tension

%%
dV = zeros(4,1);
if Q_air > 0
  dV(1)     = (Q_air/Vairway1)*(Pair - Pairway1);          %airway1 RHS
  dV(2)     = (Q_air/Vairway2)*(Pairway1 - Pairway2);      %airway2 RHS
  dV(3:end) = qi.*Pairway2/gamma - alpha*DO*(PalvO - PO); %Malv RHS
else
  dV(1)     = -(Q_air/Vairway1)*(Pairway2 - Pairway1);  %airway1 RHS
  dV(2)     = -sum(qi.*(PalvO - Pairway2))./Vairway2;  %airway2 RHS
  dV(3:end) = qi.*PalvO/gamma - alpha*DO*(PalvO - PO); %Malv RHS
end
dX = dV;
end