clf, clear, clc, clf, close all

Ts=1/1000;
rawData65 = readtable('C:\Users\ayman\OneDrive\Desktop\roket\65.csv'); % Data from LogFiles 
rawData74 = readtable('C:\Users\ayman\OneDrive\Desktop\roket\74.csv'); % Data from LogFiles 
rawDataH65 = readtable('C:\Users\ayman\OneDrive\Desktop\roket\H65.csv'); % Data from LogFiles 
rawDataH74 = readtable('C:\Users\ayman\OneDrive\Desktop\roket\H74.csv'); % Data from LogFiles 

Data65= table2array(rawData65);
Data74= table2array(rawData74);
DataH65= table2array(rawDataH65);
DataH74= table2array(rawDataH74);


%% true and corrected in g
figure
w65=Data65(:,1);
t65=Data65(:,2)*Ts;
subplot(2,2,1)
plot(t65,w65+0.4);
title('65% KNO3')
grid on
hold on
A=correction(106,300,rawData65,12,0.2);
plot(t65,w65+A+0.4);
axis([100 138 -15 40])
legend('True','Corrected')
xlabel('time [s]') 
ylabel('Thrust [g]') 


subplot(2,2,2)
wH65=DataH65(:,1);
tH65=DataH65(:,2)*Ts;
plot(tH65,wH65+0.8);
title('65% KNO3 with hole')
grid on
hold on
A=correction(127,200,rawDataH65,12.5,0.2);
plot(tH65,wH65+A+0.8);
axis([120 147 -15 40])
legend('True','Corrected')
xlabel('time [s]') 
ylabel('Thrust [g]') 

subplot(2,2,3)
w74=Data74(:,1);
t74=Data74(:,2)*Ts;
plot(t74,w74+1.2);
title('74% KNO3')
axis([292 355 -15 40])
grid on
hold on
A=correction(294,570,rawData74,12.5,0.2);
plot(t74,w74+A+1.2);
legend('True','Corrected')
xlabel('time [s]') 
ylabel('Thrust [g]') 

subplot(2,2,4)
wH74=DataH74(:,1);
tH74=DataH74(:,2)*Ts;
plot(tH74,wH74+1.2);
title('74% KNO3 with hole')
axis([220 270 -15 40])
grid on
hold on
A=correction(225,400,rawDataH74,12.5,0.2);
plot(tH74,wH74+A+1.2);
legend('True','Corrected')
xlabel('time [s]') 
ylabel('Thrust [g]')

%% corrected in N
% multiply by 1/1000 to get it in Kg and then multiply it with g=9.81 to
% convert it to newton
figure
toNewton=(1/1000)*9.81;
w65=Data65(:,1)*toNewton;
t65=Data65(:,2)*Ts;
subplot(2,2,1)
A=correction(106,300,rawData65,12,0.2)*toNewton;
A65=w65+A+0.005;
plot(t65,A65);
title('65% KNO3')
grid on
axis([100 138 -0.1 0.35])
xlabel('time [s]') 
ylabel('Thrust [N]') 

subplot(2,2,2)
wH65=DataH65(:,1)*toNewton;
tH65=DataH65(:,2)*Ts;
A=correction(127,200,rawDataH65,12.5,0.2)*toNewton;
AH65=wH65+A+0.01;
plot(tH65,AH65)
title('65% KNO3 with hole')
axis([120 147 -0.1 0.35])
grid on
xlabel('time [s]') 
ylabel('Thrust [N]') 

subplot(2,2,3)
w74=Data74(:,1)*toNewton;
t74=Data74(:,2)*Ts;
A=correction(294,570,rawData74,12.5,0.2)*toNewton;
A74=w74+A+0.001;
plot(t74,A74);
title('74% KNO3')
axis([292 355 -0.1 0.35])
grid on
xlabel('time [s]') 
ylabel('Thrust [N]') 

subplot(2,2,4)
wH74=DataH74(:,1)*toNewton;
tH74=DataH74(:,2)*Ts;
A=correction(225,400,rawDataH74,12.5,0.2)*toNewton;
AH74=wH74+A+0.01;
plot(tH74,AH74);
title('74% KNO3 with hole')
axis([220 270 -0.1 0.35])
grid on
xlabel('time [s]') 
ylabel('Thrust [N]')

%% Generating tha acc from a force vector
% a=f/m
clc
toNewton=(1/1000)*9.81;
figure
subplot(2,2,1)
w65=Data65(:,1)*toNewton;
t65=Data65(:,2)*Ts;
A=correction(106,300,rawData65,12,0.2)*toNewton;
A65=w65+A+0.005;
A=f_to_acc(106,300,A65,rawData65,12,0.2);
plot(t65,A);
axis([100 138 -5 25])
title('65% KNO3')
grid on
xlabel('time [s]') 
ylabel('Acceleration [m/s^2]')

subplot(2,2,2)
wH65=DataH65(:,1)*toNewton;
tH65=DataH65(:,2)*Ts;
A=correction(127,198,rawDataH65,12.5,0.2)*toNewton;
AH65=wH65+A+0.005;
A=f_to_acc(127,198,AH65,rawDataH65,12.5,0.2);
plot(tH65,A);
axis([120 147 -5 25])
title('65% KNO3 with hole')
grid on
xlabel('time [s]') 
ylabel('Acceleration [m/s^2]')

subplot(2,2,3)
w74=Data74(:,1)*toNewton;
t74=Data74(:,2)*Ts;
A=correction(294,570,rawData74,12,0.2)*toNewton;
A74=w74+A+0.005;
A=f_to_acc(294,570,A74,rawData74,12,0.2);
plot(t74,A);
axis([292 355 -10 25])
title('74% KNO3')
grid on
xlabel('time [s]') 
ylabel('Acceleration [m/s^2]')

subplot(2,2,4)
wH74=DataH74(:,1)*toNewton;
tH74=DataH74(:,2)*Ts;
A=correction(225,400,rawDataH74,12,0.2)*toNewton;
AH74=wH74+A+0.005;
A=f_to_acc(225,400,AH74,rawDataH74,12,0.2);
plot(tH74,A);
axis([220 270 -10 25])
title('74% KNO3 with hole')
grid on
xlabel('time [s]') 
ylabel('Acceleration [m/s^2]')
%% Find v


toNewton=(1/1000)*9.81;
subplot(2,2,1)
w65=Data65(:,1)*toNewton;
t65=Data65(:,2)*Ts;
A=correction(106,300,rawData65,12,0.2)*toNewton;
A65=w65+A+0.005;
A=f_to_acc(106,300,A65,rawData65,12,0.2)-9.81;
v = cumtrapz(t65, A)+930
plot(t65,v)
axis([100 138 -400 100])
title('65% KNO3')
grid on
xlabel('time [s]') 
ylabel('Velocity [m/s]')

subplot(2,2,2)
wH65=DataH65(:,1)*toNewton;
tH65=DataH65(:,2)*Ts;
A=correction(127,198,rawDataH65,12.5,0.2)*toNewton;
AH65=wH65+A+0.005;
A=f_to_acc(127,198,AH65,rawDataH65,12.5,0.2)-9.81;
v = cumtrapz(tH65, A)+500
plot(tH65,v)
axis([120 147 -400 100])
title('65% KNO3 with hole')
grid on
xlabel('time [s]') 
ylabel('Velocity [m/s]')

subplot(2,2,3)
w74=Data74(:,1)*toNewton;
t74=Data74(:,2)*Ts;
A=correction(294,570,rawData74,12,0.2)*toNewton;
A74=w74+A+0.005;
A=f_to_acc(294,570,A74,rawData74,12,0.2)-9.81;
v = cumtrapz(t74, A)+2800
plot(t74,v)
axis([292 355 -400 100])
title('74% KNO3')
grid on
xlabel('time [s]') 
ylabel('Velocity [m/s^2]')

subplot(2,2,4)
wH74=DataH74(:,1)*toNewton;
tH74=DataH74(:,2)*Ts;
A=correction(225,400,rawDataH74,12,0.2)*toNewton;
AH74=wH74+A+0.005;
A=f_to_acc(225,400,AH74,rawDataH74,12,0.2)-9.81;
v = cumtrapz(tH74, A)+1770
plot(tH74,v)
axis([220 270 -400 100])
title('74% KNO3 with hole')
grid on
xlabel('time [s]') 
ylabel('Velocity [m/s^2]')




%% convert force to acceleration
function A=f_to_acc(Offset_time,length,Force,Data_CSV,fuelWeight,BurningRate)
    m = @(x) -fuelWeight*exp(-BurningRate*x)+fuelWeight;
    Ts=1/1000;
    Data= table2array(Data_CSV);
    Time=Data(:,2)*Ts;
    A=zeros(size(Force));
    for i = 1:size(Time)-1
        if abs(Time(i)-Offset_time)<0.2
            Offset=i;
            break
        end
    end
   
    for i=Offset:Offset+length
        t=Time(i)-Offset_time;
        mm=(20-m(t))*1/1000; % initial weight = 20 g
        A(i)= Force(i)/mm;
    end
end




%% Generating data that compansate for fuel burning weight   
function A=correction(Offset,Length,Data_CSV,fuelWeight,BurningRate)
    f = @(x) -fuelWeight*exp(-BurningRate*x)+fuelWeight;
    Ts=1/1000;
    Data= table2array(Data_CSV);
    t=Data(:,2)*Ts;
    offset=Offset;
    length=Length;

    for i = 1:size(t)-1
        if abs(t(i)-offset)<1
            offset=i;
            break
        end
    end
    A=zeros(size(t));
    AF=Data(offset:offset+length,2);
    format short
    for i = 2:length+1
        AF(i)=AF(i)-AF(1);
    end
    AF(1)=0;    
    AF=AF*10^-3;
    for i = 2:length+1
        AF(i)=f(AF(i));
    end
    A(offset:offset+length,1)=AF;
end