%Choked Flow Eqn MODELLED ISOTHERMALLY

D_orifice = 0.005; %Diameter of orifice (m) (piercing)
D_can = 0.037*2; %Diameter of Can (m)

T_1 = 22+273; %Ambient Temperature (K)
R = 188.9; %Specific Gas Constant for CO2

P_1(1) = 61*100000; %Pressure in canister (Pa)
P_2 = 101325; %Atmospheric Pressure (Pa)

Cd = 0.61; %Discharge Coefficient

Cp = 4.416; %Isobaric heat capacity at 61 bar and 22C (kJ/mol)
Cv = 1.0103; %Isochoric heat capacity at 61 bar and 22C (kJ/mol)
k = Cp/Cv; %Specific heat ratio

A_1 = (pi*(D_can^2))/4; %Area of can (m2)
A_2 = (pi*(D_orifice^2))/4; %Area of orifice (piercing) (m2)

Volume_Canister = pi*((D_orifice/2)^2)*0.076264; %Volume of Canister (m^3)


mass(1) = 0.06; %initial mass
rho(1) = P_1(1)/(R*T_1); %initial density
m_dot(1) = Cd.*A_2.*(sqrt(k*rho(1)*P_1(1)*((2/(k+1))^((k+1)/(k-1))))); %initial mass flow
Velocity_Orifice(1) = m_dot(1)/(rho(1)*A_2); %Velocity at orifice (m/s)

dt = 0.001; %time step
t(1) = 0; %initial time

n = 1;

while mass(n) > 0
    if P_1(n) <= 101325
        P_1(n+1) = 101325;
    else

rho(n+1) = P_1(n+1)/(R*T_1); %Density

m_dot(n+1) = Cd.*A_2.*(sqrt(k*rho(n+1)*P_1(n)*((2/(k+1))^((k+1)/(k-1))))); %Mass Flow Rate (kg/s)

Velocity_Orifice(n+1) = m_dot(n+1)/(rho(n+1)*A_2); %Velocity at orifice (m/s)

mass(n+1) = mass(n) - m_dot(n+1)*dt %Mass

P_1(n+1) = (((mass(n+1)/44)*R*T_1)/((Volume_Canister)-(mass(n+1)/44)*(4.267e-5))) - ((((mass(n+1)/44)^2)*0.364)/(Volume_Canister.^2)); %Pressure Calculated with vdw

t(n+1) = t(n) + dt; %Time

n = n+1
    end
end

plot(t,mass)