function bubble(r0,pp,tstop)
format long;
global aa;
global cc;
global dd;
global rho0;
global p0;
global pa0;
global ww; 
global vis;
global surf;
global beta;
global gamma;

rho0 = 1000.;
p0   = 1.013e5;
c0   = 1481.;
vis = 1e-6;
surf = 0.073;
pa0  = -pp*p0;
ww   = 26.5*2.*pi*1e3;
% argon
%beta = r0/8.86;
%gamma = 1.667;

% air
beta = r0/8.54;
gamma = 1.4;
%
aa = -3./2.;
cc = rho0*c0;
dd = (p0+2.*surf/r0)*(r0^3-beta^3)^gamma;
options = odeset('BDF','on','RelTol',1e-10,'AbsTol',1e-8);
[t,y] = ode15s('rp', [0 tstop], [r0; 0],options);
subplot(3,1,1); plot(t,y(:,1));
subplot(3,1,2); plot(t,y(:,2));
subplot(3,1,3); plot(t,p0+pa0*cos(ww*t));

function dy = rp(t,y)
global aa;
global cc;
global dd;
global rho0;
global p0;
global pa0;
global ww;
global vis;
global surf;
global beta;
global gamma;

pvel = aa*y(2)*y(2)/y(1);
pgas = dd/(y(1)^3-beta^3)^gamma/y(1)/rho0;
pext = (p0+pa0*cos(ww*t))/(rho0*y(1));
dpgas = dd/cc*(-gamma)*3.*y(1)^2*y(2)/(y(1)^3-beta^3)^(gamma+1.);
%dpext = pa0*ww/cc*cos(ww*t);
dpext = 0;
pvis = 4.*vis*y(2)/y(1)^2;
psur = 2.*surf/rho0/y(1)^2;
dy = [y(2); pvel+pgas-pext+dpgas-dpext-pvis-psur];

