% LEP Web module: L.A. Basin
%

tspan = [0 72]; 
y0 = [2.04*10^-8;2.04*10^-8;2.04*10^-8];
[v y]=ode45(@ODEfun,tspan,y0);
plot(v,y,'LineWidth',2);
ylabel('Cco(lbmol/ft^3)');
xlabel('time(hr)');
legend('LA','Orange','Riverside')
xlim([0 72])
ylim([0 3*10^-8])
grid on

%
%

function dYfuncvecdt = ODEfun(t,Yfuncvec); 
CO1 = Yfuncvec(1);
CO2= Yfuncvec(2);
CO3=Yfuncvec(3);
% Explicit equations
A = 372176640000; % Total area of South CA Basin (ft2)
Cars =10800000 ; % Total number of vehicles in SoCAB
A1 = A * .5; %Area of Los Angeles
A2 = A * .5 * .2; % Area of Orange County
A3 = A * .5 * .6; % Area of Riverside County
Cars1 = Cars * .55; % Number of cars in LA
Cars2 = Cars * .25; % Number of cars in Orange
Cars3 = Cars * .20; % Number of cars in Riverside
H = 2200; % Height of Basin (ft)
Cars1ppm = 2.5;
Cars2ppm = 1.8;
Cars3ppm = 2.2;
TBasin = 534.7; % Rankine
PBasin = 1; % atm
R = 0.7302; % ft3*atm/lbmol*R
vwind =79200; % ft/h wind velocity
corridor = 105600 * 2000; % ft2 wind corridor
vemissions =3000; % ft3/h STD vehicular emissions
Temissions = 491.67; % R eq. 0?C
n1 = PBasin * A1 * H / (R * TBasin); % lbmol of air in Vx
n2 = PBasin * A2 * H / (R * TBasin); % lbmol of air in Vx
n3 = PBasin * A3 * H / (R * TBasin);
ncars1 = PBasin * vemissions / (R * Temissions) * Cars1; % lbmol gas emissions from Cars-x
ncars2 = PBasin * vemissions / (R * Temissions) * Cars2;
ncars3 = PBasin * vemissions / (R * Temissions) * Cars3;
COcars1 = Cars1ppm / 1000000 * ncars1; % lbmol CO emissions from cars
COcars2 = Cars2ppm / 1000000 * ncars2;
COcars3 = Cars3ppm / 1000000 * ncars3;
Qo = vwind * corridor; % vo
Q1 = Qo + vemissions * Cars1;
Q2 = Q1 + vemissions * Cars2;
Q3 = Q2 + vemissions * Cars3;
COo = .08; % ppm
COwind0 = 2e-10; % lbmol CO / ft3
COs0 = COwind0 * Qo; % lbmol CO in wind
COs1 = CO1 * Q1;
COs2 = CO2 * Q2;
CO1ppm = CO1 * 1000000 * A1 * H / n1;
CO2ppm = CO2 * 1000000 * A2 * H / n2;
CO3ppm = CO3 * 1000000 * A3 * H / n3;
% Differential equations
dCO1dt = (COcars1 + COs0 - Q1 * CO1) / (A1 * H);
dCO2dt = (COcars2 + COs1 - Q2 * CO2) / (A2 * H);
dCO3dt = (COcars3 + COs2 - Q3 * CO3) / (A3 * H);
dYfuncvecdt=[dCO1dt;dCO2dt;dCO3dt];
end
