clear all
close all
clc

%% load file and data
rockets = 'rocket_data.csv';
rocket_data = csvread(rockets, 1, 0, [1, 0, 590, 6]);
vertical_accel = rocket_data(:, 2); %m/s^2 
lateral_accel = rocket_data(:, 3); %m/s^2
mass_g = rocket_data(:, 4); %g
angle = rocket_data(:, 7); %deg

%% set known constants
%time
time = rocket_data(:, 1); %s

%thrust & drag (axial)
thrust = rocket_data(:, 5); %N
drag = rocket_data(:, 6); %N

%weight (axial & transverse)
mass = mass_g./1000; %convert to kg
weight = mass.*9.81; %N
a_weight = weight.*sind(angle);
t_weight = weight.*cosd(angle);

%% back out aero forces from total accel
%total accel to total force
vertical_force = mass.*vertical_accel;
lateral_force = mass.*lateral_accel;

%total force components in rocket coordinates
a_force = vertical_force.*sind(angle) + lateral_force.*cosd(angle);
t_force = vertical_force.*cosd(angle) + lateral_force.*sind(angle);

%find aero forces from net force
a_aero = a_force - thrust + drag + a_weight;
t_aero = t_force - t_weight;

%% calculate compressive forces
a_compression = thrust + a_weight + a_aero + drag;
t_compression = t_weight + t_aero;

%in lbs, so I understand it
a_comp_lbs = a_compression.*0.224809;
t_comp_lbs = t_compression.*0.224809;

%with a safety factor (1.5 is standard for aircraft; 1.25 for unmanned
%spacecraft - keep at 1.5 until we stop guesstimating data)
sf = 1.5;
a_comp_lbs_s = a_comp_lbs.*sf;
a_compression_s = a_compression.*sf;
t_comp_lbs_s = t_comp_lbs.*sf;
t_compression_s = t_compression.*sf;

%% return loading

disp(['The axial load is ' num2str(max(a_comp_lbs)) ' pounds (' num2str(max(a_compression)) ' N)'])
disp(['The perpendicular load is ' num2str(max(t_comp_lbs)) ' pounds (' num2str(max(t_compression)) ' N)'])
disp(['The axial load with a ' num2str(sf) 'x safety factor is ' num2str(max(a_comp_lbs_s)) ' pounds (' num2str(max(a_compression_s)) ' N)'])
disp(['The perpendicular load with a ' num2str(sf) 'x safety factor is ' num2str(max(t_comp_lbs_s)) ' pounds (' num2str(max(t_compression_s)) ' N)'])

%% graph
figure(1)
plot(time, a_comp_lbs)
title('Axial Compression Loads')
xlabel('time (s)')
ylabel('axial compression load (lbs)')

figure(2)
plot(time, t_comp_lbs)
title('Normal Compression Loads')

%% convert to pressure
%axial pressure on tube
tube_OR = 3.0625; %in
tube_IR = 3; %in
a_area_tube = pi.*(tube_OR.^2 - tube_IR.^2); %sq in
a_pressure_tube = a_comp_lbs./a_area_tube; %psi

%axial pressure on NC
height = 36; %in
a_area_NC = pi.*tube_OR.*(tube_OR + sqrt(height.^2 + tube_OR.^2));
a_pressure_NC =a_comp_lbs./a_area_NC; %psi

%normal pressure on tubes
tube_length = 8*12; %in
t_area = 2.*pi.*tube_OR.*tube_length;
t_pressure = t_comp_lbs./t_area;

%safety factor
a_press_tube_s = a_pressure_tube*sf;
a_press_NC_s = a_pressure_NC*sf;
t_press_s = t_pressure*sf;

%% return pressure data

disp(['The axial pressure on the tubes is: ' num2str(max(a_pressure_tube)) ' psi, or ' num2str(max(a_press_tube_s)) ' psi with a ' num2str(sf) ' safety factor.'])
disp(['The axial pressure on the nose cone is: ' num2str(max(a_pressure_NC)) ' psi, or ' num2str(max(a_press_NC_s)) ' psi with a ' num2str(sf) ' safety factor.'])
disp(['The transverse pressure on the tubes is: ' num2str(max(t_pressure)) ' psi, or ' num2str(max(t_press_s)) ' psi with a ' num2str(sf) ' safety factor.'])

%% graph
figure(3)
plot(time, a_pressure_tube)
title('Axial Pressure on Tube')
xlabel('time (s)')
ylabel('axial pressure (psi)')

figure(4)
plot(time, a_pressure_NC)
title('Axial Pressure on Nose Cone')
xlabel('time (s)')
ylabel('axial pressure (psi)')

figure(5)
plot(time, t_pressure)
title('Normal Pressure on Tubes')
xlabel('time (s)')
ylabel('normal pressure (psi)')

