%Mason Averill
%ME-544 Fall 2020, 10/31
%Main Project--Estes Rockets Expected Height Achieved

%Rocket Parameters 
Diameter=25/1000;%Diameter in meters
Mass_rocket=37/1000;%mass of rocket in kg
Mass_engine=24/1000;%mass of rocket engine(including propellent)
Mass_propellent=11/1000;%mass of propellent in rocket engine



Mass=Mass_rocket+Mass_engine;%mass in kg
Weight=Mass*9.81;%Weight in N
Cross_sectional_area=pi*Diameter^2/4;%effective area

%Air Properties
density_air=1.225;%kg/m^3
dynamic_viscosity=1.81*10^-5;


%Desired Resolution
%time_steps=5*10^-2;
time_steps=0.001;
%Engine type:Estes c6 Raw Data

time=[0 0.014	0.026	0.067	0.099	0.15	0.183	0.207	0.219	0.262	0.333	0.349	0.392	0.475	0.653	0.913	1.366	1.607	1.745	1.978	2.023	2.024];

thrust=[0 0.633	1.533	2.726	5.136	9.103	11.465	11.635	11.391	6.377	5.014	5.209	4.722	4.771	4.746	4.673	4.625	4.625	4.868	4.795	0.828	0];

time_range=0:time_steps:max(time);

linearly_interpolated_points=interp1(time,thrust,time_range);


figure(1)
plot(time,thrust,'o',time_range,linearly_interpolated_points, ':.')
xlabel('Time (s)')
ylabel('Thrust (N)')
legend('Raw Data','Interpolated Data')
title('Thrust vs Time: Estes C6 Engine')
grid on

impulse=trapz(time_range,linearly_interpolated_points);
impulse_of_t=cumtrapz(time_range,linearly_interpolated_points);

%Cd determination/curve fitting
Re_raw=[0.1 10 100 1000 6*10^3 10^4 2*10^5 3*10^5 4*10^5 5*10^5 10^6];
Cd_raw=[38 2.4 1.8 0.98 .98 1.1 1.1 .98 .45 .25 .27];

%coefficients=polyfit(Re_raw,Cd_raw,degree);

Re_range=0.1:0.1:10^6;

linearly_interpolated_Cd=interp1(Re_raw,Cd_raw,Re_range);%or linearly interpolated method



flight_time_store=[];
position_store=[];
velocity_store=[];
acceleration_store=[];
position=0;
velocity=0;
acceleration=0;
position_store(1,1)=position;
velocity_store(1,1)=velocity;
acceleration_store(1,1)=acceleration;
F_net_store=[];
F_net=0;
flight_time=0;
flight_time_store(1,1)=0;
i=1;
%c_d=0.75;
%c_d=0;
Current_mass=Mass;
Current_mass_store=[];
while(position>=0)
    %determine c_d
  Re=density_air*abs(velocity_store(1,i))*Diameter/dynamic_viscosity;
  Re=round(Re);
  [row,column]=find(Re==Re_range);
  if(row==1)%this Re exists in possible Re values considered
     c_d=0.75*linearly_interpolated_Cd(1,column);
      
  elseif(Re<0.1)
     c_d=Cd_raw(1,1);
      %c_d=0.60;
  else
      %c_d=Cd_raw(1,11);
     c_d=0.6;
   end 
  
  %no liftoff until thrust exceeds weight of rocket
  if(flight_time<0.1 && linearly_interpolated_points(1,i)<Weight)
     F_net=0;
     Current_mass=Mass-(impulse_of_t(i)/impulse)*Mass_propellent;
     
  %thrust being generated and liftoff has occured
  elseif(flight_time<=max(time))
      
     if(velocity_store(1,i)>0)
        F_d=-0.5*c_d*density_air*Cross_sectional_area*(velocity_store(1,i))^2;
     else
        F_d=0.5*c_d*density_air*Cross_sectional_area*(velocity_store(1,i))^2;
     end  
      F_net=linearly_interpolated_points(1,i)-Weight+F_d;
     Current_mass=Mass-(impulse_of_t(i)/impulse)*Mass_propellent;
   %no thrust is being generated
  else 
     if(velocity_store(1,i)>0)
        F_d=-0.5*c_d*density_air*Cross_sectional_area*(velocity_store(1,i))^2;
    else
        F_d=0.5*c_d*density_air*Cross_sectional_area*(velocity_store(1,i))^2;
    end 
      
    F_net=-Weight+F_d;
    Current_mass=Mass-Mass_propellent;
  end 
    

acceleration=F_net/Current_mass;
velocity=velocity+acceleration*time_steps;
position=position+velocity*time_steps+0.5*acceleration*time_steps^2;
position_store(1,i+1)=position;
velocity_store(1,i+1)=velocity;
acceleration_store(1,i+1)=acceleration;
flight_time_store(1,i+1)=flight_time;
F_net_store(1,i)=F_net;
Current_mass_store(1,i)=Current_mass;

i=i+1;
flight_time=flight_time+time_steps %supress this if the time readout bothers you
end 

figure(2)
plot(flight_time_store,position_store)

hold on 
plot(flight_time_store,velocity_store)

hold on 
plot(flight_time_store,acceleration_store)

legend('Height of Rocket (m)','Velocity (m/s)','Acceleration (m/s^2)')

xlabel('Time (s)')
ylabel('Height, Velocity, Acceleration')
title('Height, Velocity, and Acceleration of Rocket vs Time')
grid on

figure(3)
plot(flight_time_store,position_store)
xlabel('Time (s)')
ylabel('Height (m)')
title('Height of Rocket vs Time')
ylim([0 400])
grid on


figure(4)
plot(flight_time_store,velocity_store)
xlabel('Time (s)')
ylabel('Velocity (m/s)')
title('Velocity of Rocket vs Time')
grid on


figure(5)
plot(flight_time_store,acceleration_store)
xlabel('Time (s)')
ylabel('Acceleration (m/s^2)')
title('Acceleration of Rocket vs Time')
grid on



max(position_store)
max(velocity_store)
flight_time_store(i)=[];

figure(6)
plot(flight_time_store,Current_mass_store);
xlabel('Time (s)')
ylabel('Mass (kg)')
title('Mass of Rocket vs Time')
ylim([0.045 .065])
grid on