clc
clear

//Ecuaciones de articulo Finney, G.A., 2000. Analysis of a water-propelled rocket: A problem in honors physics. American Journal of Physics, 68(3), pp.223-227.

// ------------------------- PARAMETERS --------------------------
g=9.81                  // Gravity acceleration          [m/s^2]
P_a=101325               // Atmospheric preassure    [Pa]
rho_w=1000                // Water denisty            [kg/m^3]
M_r=50e-3;                // 50 gr empty rocket mass  [kg]
V_T=2.0e-3;               // Total volume of rocket = 2L [m^3]
D_nozzle=22e-3          // 22 mm Nozzle diameter [m]
A_e=%pi*(D_nozzle/2)^2;      // Cross-secctional area of exhaust nozzle [m^2]
alpha=90*%pi/180        // Angle                    [º]
D=2e-4;                   //Drag coefficient    

// ---------------------- INITIAL VALUES ------------------------
x0=0; v0x=0; y0=0; v0y=0;    // Initial position and velocity
P_0=4*P_a                     // Initial Pressure           [Pa]
State0=[y0;v0y;P_0]            // Initial state
V_w0=800e-6;                 //800 mL Initial Water Volume  
V_0=V_T-V_w0                    // Initial Aire Volume       [m^3]

function dot_f = f1(t,state)//first stage of rocket flight
  global burnout
  //y=state(1)
  vy=state(2)
  P=state(3)
//  printf("burnout. t=%f\n",t);
  V_w=V_T-P_0*V_0/P; //volume of water
  M_w=rho_w*V_w  // mass of water
  if M_w<0 burnout=1; end
  if burnout==1
    printf("burnout. t=%f P=%f\n",t,P);
     M=M_r;//after "burnout" 
     a=-D*vy*abs(vy)/M-g
     dP=0
  else
     M=M_w+M_r;//before "burnout"
     printf("Masa=%f P=%f t=%f\n",M,P,t);
     a=(2*A_e*(P-P_a)-D*vy*abs(vy))/M-g
     dP=-P^2*A_e/(P_0*V_0)*sqrt(2*(P-P_a)/rho_w)
  end
  dot_f=[vy;a;dP]
endfunction

global burnout
burnout=0;
disp("Water Rocket")
t0=0                         // Initial time               [s]
t_end = 9;
step = (t_end-t0)/1000;
t = t0:step:t_end
state = ode(State0,t0,t,f1);
n = length(state(1,:));
scf(1)
plot(t(1:n),state(1,1:n),"rd");
xtitle("Y position","t (s)","height (m)")
xgrid();
scf(2)
plot(t(1:n),state(2,1:n),"rd");
xtitle("Y velocity","t (s)","y velocity (m/s)")
xgrid();
scf(3)
plot(t(1:n),state(3,1:n),"rd");
xtitle("Pressure","t (s)","Pressure (Pa)")
xgrid();