clc
clear all
close all

%% Exercise 1 (20p)

%Parameters
c= ... ;                            % speed of sound in water, m/s
rho = ... ;                         % water density, kg/m^3
sigma = ... ;                       % surface tension, N/m
mu = ...;                           % shear viscosity, Pa*s
v0 = 0;                             % initial bubble boundary velocity, m/s 
pinf = 101325;                      % ambient pressure, Pa
amp = 1e5;                          % pressure amplitude  
gamma = 1.4;                        % ratio between the specific heat capacities
f = ... ;                           % driving frequency, Hz
t_max = 13e-6;                      % maximum observation time
dt = 1/(100*f);                     % time step, us
R0 = 5e-6;                          % initial bubble radius

%Variables
syms R(t) pa(t)

Rt = diff(R,t);                     % first derivative of bubble radius, dR/dt
Rtt = diff(R,t,2);                  % second derivative of bubble radius, dR/dt
pa = amp*sin...;                    % time delayed driving pressure, Pa
p0 = pinf + 2*sigma/R0;             % internal pressure of the bubble at equilibrium
pg = p0*(R0/R)^(3*gamma);           % gas perssure in the interior of the bubble
pb = pg -  2*sigma/R - 4*mu*(Rt/R); % pressure on the liquid side of the bubble interface
pbt = diff(pb,t);                   %dpb/dt 

% Define the differential equation 
eqn1 = ... == ... ;                 % implement here the differential equation

%Solve differential equation
[V] = odeToVectorField(eqn1);
M = matlabFunction(V,'vars', {'t','Y'});
sol = ode45(M,[0 t_max],[R0 v0]);

%Plots
R = deval(sol,[0:dt:t_max],1);      % bubble radius, m
v = deval(sol,[0:dt:t_max],2);      % bubble velocity, m/s
a = ... ;                           % bubble acceleration, m/s^2