%Math4090 2010 chapter 2
%
clear all
global d P0 rho g h f L D gam
L=1e3;
D=2e-3*L;
d=1*D;
P0=1.013e5;
rho=1e3;
h=1*1e2;
gam=1.4;
g=9.8;
f=1*2.5e-2;
x0=[0;1*sqrt(2*g*h)];
T=20;
N=101;
[t,x]=ode23(@tunnelEQ,linspace(0,T,N),x0);
p=P0/(rho*g*h)*(1./(1-x(:,1))).^gam;
figure(1);
subplot(2,1,1);
hold on
plot(t,x(:,1),'b-');
ylabel('x');
subplot(2,1,2);
hold on
plot(t,p,'b-');
ylabel('p_{air}');
xlabel('t');