clear, clc;
close all force;     %Force closes all plots

fileName = [];
dataTable = [];
dateOfFire = [];
testNumber = [];
subFolderName = [];
outputFileName = [];
txtFileID = [];

sampleRatePT = [];
startTimePT = [];
sampleRateLC = [];
startTimeLC = [];
sampleRateTC = [];
startTimeTC = [];

ambientTemperature = [];

thrustTable = [];
pressureTable = [];
temperatureTable = [];
timeRange = [];
burnStartTime = [];
burnEndTime = [];
burnTime = [];
selection = [];

maxThrust = [];
maxThrustIndex = [];
maxPressure = [];
maxPressureIndex = [];
maxTemperature = [];
maxTemperatureIndex = [];

weightLost = [];
propellantMass = [];
impulse = [];
isp = [];
avgMassFlow = [];

thrustCoefficent = [];
thrustCoefficents = [];     
avgThrustCoefficent = [];   

%SRM Constants
throatArea = (pi*(.361/2)^2)/1550;   %m^2
throatArea_in = (pi*(.361/2)^2); %in^2

inputData = [];

OD = 1.72;  % in 
OD1 = 1.72; %grain 1 outer diameter
OD2 = 1.72;

ID = 1.125; % in
ID1 = 1.125;
ID2 = 1.125; % in

L = [];     % in      avg. grain length
L1=[];              %grain 1 length
L2=[];
n = 2;      % in


%% User Input

disp('Choose thrust data file (must be in folder of postprocessing *m file!):')
loadCellFile = uigetfile('*.csv');
disp(loadCellFile);
disp(' ');
disp('Loading...');
disp(' ');
thrustTable = readtable(loadCellFile,'Format','auto');

disp('Choose pressure data file:')
pressureFile = uigetfile('*.csv');
disp(pressureFile);
disp(' ');
pressureTable = readtable(pressureFile,'Format','auto');

disp('Choose temperature data file:')
temperatureFile = uigetfile('*.csv');
disp(temperatureFile);
disp(' ');
temperatureTable = readtable(temperatureFile,'Format','auto');


%% Input File Read
inputData = readtable('DataInputFileForMatlab.xlsx','Format','auto');
inputData = table2array(inputData);                    %%converts table to array
inputData([1],:) = [];

dateOfFire = char(inputData(1,2));    %%Text must be converted to char data type
testNumber = char(inputData(2,2));
propellantMass = str2double(inputData(3,2));

L1 = str2double(inputData(5,2));
OD1 = str2double(inputData(6,2));
ID1 = str2double(inputData(7,2));

L2 = str2double(inputData(9,2));
OD2 = str2double(inputData(10,2));
ID2 = str2double(inputData(11,2));

ambientTemperature = char(inputData(13,2));

note = char(inputData(15,2));

propellantMass = propellantMass/1000;   %converting mass to kg
L = (L1 + L2) / 2;                      %L is avg grain length. This variable is passed to burn rate calc
OD = (OD1 + OD2) / 2;                   %OD is avg grain OD. This variable is passed to burn rate calc
ID = (ID1 + ID2) / 2;                   %ID is avg grain ID. This variable is passed to burn rate calc

%% Creating subfolder and output txt file
subFolderName = [dateOfFire,'_SRMTest#',testNumber];

status = rmdir(subFolderName, 's');                             %removes subfolder if already exists (status= is to prevent error from throwing)
mkdir(subFolderName);                                           %creates subfolder with name including date and test number
mkdir (fullfile(subFolderName,'Graphs'))                        %creates folder for graphs
mkdir (fullfile(subFolderName,'Media'))
mkdir (fullfile(subFolderName,'RawData'))
outputFileName = ['ProcessedData_Test#',testNumber,'_(',dateOfFire,')','.txt'];
txtFileID=fopen([pwd '/',subFolderName,'/',outputFileName],'w');    %creates text file to print data to

outputFileName = ['BurnRateData_Test#',testNumber,'_(',dateOfFire,')','.txt'];
txtFileID_burnRate=fopen([pwd '/',subFolderName,'/',outputFileName],'w');    %creates text file to print burn rate data to

copyfile('DataInputFileForMatlab.xlsx', subFolderName);       %copying what inputfile was used for our records

copyfile(loadCellFile, [pwd '/',subFolderName,'/','RawData']);     %Copying rawdata files into RawData folder
copyfile(pressureFile, [pwd '/',subFolderName,'/','RawData']); 
copyfile(temperatureFile, [pwd '/',subFolderName,'/','RawData']); 


%% Organizing Data Import
%Extracting sample rate, start time:
sampleRateLC = thrustTable.(2)(1);     %NOTE Colums and rows are flipped for dot notation on tables,column 2 row 1, etc.
startTimeLC = thrustTable.(2)(2);

sampleRatePT = pressureTable.(2)(1);    
startTimePT = pressureTable.(2)(2);

sampleRateTC = temperatureTable.(2)(1);     
startTimeTC = temperatureTable.(2)(2);

%Trimming and Converting to tables
thrustTable([1,2,3],:) = [];                                            % delete rows containing strings 
thrustTable = str2double(thrustTable.(1));            % changing remaining rows to double. This also turns table into array
for i=1:numel(thrustTable(:,1))                                       %Filling in pressure time data. Getting number of rows
    thrustTable(i,2) = startTimeLC + (1/sampleRateLC)*(i-1);          %adding time from start to each row
end
thrustTable = array2table(thrustTable,'VariableNames',{'Thrust','Time'});% converting array back to table. Adding column titles
thrustTable.(1) = 9.81*thrustTable.(1);     %%%%%%%%%%converting kg to N

pressureTable([1,2,3],:) = [];                                            
pressureTable = str2double(pressureTable.(1));           
for i=1:numel(pressureTable(:,1))                                      
    pressureTable(i,2) = startTimePT + (1/sampleRatePT)*(i-1);          
end
pressureTable = array2table(pressureTable,'VariableNames',{'Pressure','Time'});

temperatureTable([1,2,3],:) = [];                                            
temperatureTable = str2double(temperatureTable.(1));           
for i=1:numel(temperatureTable(:,1))                                     
    temperatureTable(i,2) = startTimeTC + (1/sampleRateTC)*(i-1);       
end
temperatureTable = array2table(temperatureTable,'VariableNames',{'Temperature','Time'});

%ambientTemperature = mean(temperatureTable{5:20,1});   %getting ambient temp from Thermocouple avg.
  

%% User Selecting Data Range
set(0,'DefaultFigureWindowStyle','docked');              %%docks graphs
fprintf('\n\nUn-cut thrust curve printing.....\n\n')


%plotting so user can select data range.
subplot(2,1,1)
thrustCurve = plot(thrustTable.(2),thrustTable.(1));
ylim([-10 500]);     %%%%clipping so noise spikes dont't saturate graph
grid on
title(['Thrust Curve - ','Test #',testNumber,', ',dateOfFire,''])
xlabel('Time (s)')
ylabel('Thrust (N)')

subplot(2,1,2)
pressureCurve = plot(pressureTable.(2),pressureTable.(1));
ylim([-10 500]);
grid on
title(['Pressure Curve - ','Test #',testNumber,', ',dateOfFire,''])
xlabel('Time (s)')
ylabel('Pressure (psi)')

selection = 'y'; 
while selection == 'y'
    disp('Enter new time range for thrust data [t1 t2]:');
    timeRange = input('');                                  %getting desired data range from user
    burnStartTime = timeRange(1);                           %start of burn
    burnEndTime = timeRange(2);
    
    %triming all tables using time range specified by user
    thrustTable(thrustTable.(2) < burnStartTime,:) = [];     %deleting all rows less than the user selected start time
    thrustTable(thrustTable.(2) > burnEndTime,:) = [];     %deleting all rows greater than the user selected end time
    pressureTable(pressureTable.(2) < burnStartTime,:) = [];    
    pressureTable(pressureTable.(2) > burnEndTime,:) = [];
    temperatureTable(temperatureTable.(2) < burnStartTime,:) = [];    
    temperatureTable(temperatureTable.(2) > burnEndTime,:) = [];
    
    %plotting again
    subplot(2,1,1)
    thrustCurve = plot(thrustTable.(2),thrustTable.(1));
    ylim([-10 500]);     %%%%clipping so noise spikes dont't saturate graph
    grid on
    title(['Thrust Curve - ','Test #',testNumber,', ',dateOfFire,''])
    xlabel('Time (s)')
    ylabel('Thrust (N)')

    subplot(2,1,2)
    pressureCurve = plot(pressureTable.(2),pressureTable.(1));
    ylim([-10 500]);
    grid on
    title(['Pressure Curve - ','Test #',testNumber,', ',dateOfFire,''])
    xlabel('Time (s)')
    ylabel('Pressure (psi)')    
    
    disp('New time range? (y/n):')
    selection = input('', 's'); 
    while (selection ~= 'y' & selection ~= 'n')              %catch in case user enters something other than y or n
        disp('INVALID. Enter y or n. New time range? (y/n):')
        selection = input('', 's'); 
    end
end

%Subtracting all time values by  'burnStartTime' so that times are zeroed around start of burn
    thrustTable.(2) = (thrustTable.(2)-burnStartTime);
    pressureTable.(2) = (pressureTable.(2)-burnStartTime);
    temperatureTable.(2) = (temperatureTable.(2)-burnStartTime);
    
close;

%% General calulations
burnTime = burnEndTime - burnStartTime;
avgMassFlow = propellantMass/burnTime;
density = propellantMass/((L1*pi()*((OD1/2)^2-(ID1/2)^2))+(L2*pi()*((OD2/2)^2-(ID2/2)^2)));     

weightLost = (9.8066*avgMassFlow).*thrustTable.(2);              % g(m/s^2) * massflow(kg/s) * timefromstart(s) = forceLost(N), on loadcell lost from ejected mass
thrustTable.(1) = thrustTable.(1)+weightLost;                    %adding weight lost back to thrust points
thrustTable.(1) = thrustTable.(1)-mean(thrustTable.(1)(1:25));   %zeroing out thrust

impulse = trapz(thrustTable.(2),thrustTable.(1));
isp = impulse/(propellantMass*9.8066);
c_Star = (trapz(pressureTable.(2),(pressureTable.(1)*6894.76)))*throatArea/propellantMass;     %note: psi is converted to pa

thrustCoefficent = isp*9.8066/c_Star;

 thrustCoefficents = thrustTable.(1)./(pressureTable.(1)*throatArea*6894.76);      %other way to calculate Cf
 avgThrustCoefficent = mean(thrustCoefficents);

%Max Values
[maxThrust,maxThrustIndex] = max(thrustTable.(1));
[maxPressure,maxPressureIndex] = max(pressureTable.(1));
[maxTemperature,maxTemperatureIndex] = max(temperatureTable.(1));

% %% Burn rate
% %%get steady state from user
% %Graphs
% subplot(2,1,1)
% thrustCurve = plot(thrustTable.(2),thrustTable.(1));
% ylim([-10 500]);     %%%%clipping so noise spikes dont't saturate graph
% grid on
% title(['Thrust Curve - ','Test #',testNumber,', ',dateOfFire,''])
% xlabel('Time (s)')
% ylabel('Thrust (N)')
% 
% subplot(2,1,2)
% pressureCurve = plot(pressureTable.(2),pressureTable.(1));
% ylim([-10 500]);
% grid on
% title(['Pressure Curve - ','Test #',testNumber,', ',dateOfFire,''])
% xlabel('Time (s)')
% ylabel('Pressure (psi)')
% 
% %get start time
% disp('Enter start and end time of steady state data [t1 t2] (s):');
% steadyInput = input('');
% steadyStartTime = steadyInput(1);
% steadyEndTime = steadyInput(2);
% 
% close;
% 
% % make modified data (for analysis only cut off beginning)
% pressure_steady_table = pressureTable(pressureTable.(2) >= steadyStartTime,:);    
% steadyTime = pressure_steady_table.(2);
% steadyPressure = pressure_steady_table.(1);
% 
% %execute burn rate function
% burnRate = burn_rate_calc(steadyTime,steadyPressure,c_Star,throatArea_in,OD1,ID1,L1,OD2,ID2,L2,propellantMass);
% 
% %%curve fitting (now for fit, cut off end too)
% pressure_steady_table = pressure_steady_table(pressure_steady_table.(2) <= steadyEndTime,:);    
% steadyPressure_full = pressure_steady_table.(1);
% burnRate_short = burnRate(1:size(steadyPressure_full));
% [xData_burnrate_original, yData_burnrate_original] = prepareCurveData(steadyPressure_full, burnRate_short);
% 
% M = [xData_burnrate_original yData_burnrate_original];  % concatanate data into matrix
% t = array2table(M,'VariableNames',{'pressure' 'burnrate'});
% averaged = varfun(@mean,t,'GroupingVariable','pressure');
% [xData_burnrate, yData_burnrate] = prepareCurveData(averaged.pressure, averaged.mean_burnrate);
% 
% % polynomial fit
% ft = fittype( 'poly5' );
% opts = fitoptions( 'Method', 'LinearLeastSquares' );
% opts.Normalize = 'on';
% opts.Robust = 'LAR';
% [fitresult_poly, gof_poly] = fit(xData_burnrate, yData_burnrate, ft, opts);
% 
% % exponential fit
% ft = fittype( 'power1' );
% opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
% opts.Display = 'Off';
% opts.Robust = 'LAR';
% opts.StartPoint = [0.00207376669358892 0.742132833375259];
% [fitresult_exp, gof_exp] = fit(xData_burnrate, yData_burnrate, ft, opts);


%% Graphing
%Individual Graphs
%Thrust
thrustCurve = plot(thrustTable.(2),thrustTable.(1));
    %xticks(0:0.2:100);
    %yticks(-100:100:10000);        %setting spacing to ~1/15 of max vlaue. then rounding to nearest 10th     round(maxThrust/15,-1)
grid on
title(['Thrust Curve - ','Test #',testNumber,', ',dateOfFire,''])
xlabel('Time (s)')
ylabel('Thrust (N)')

outputFileName = ['ThrustCurve_Test#',testNumber,'_(',dateOfFire,')','.jpg'];
saveas(gcf,[pwd '/',subFolderName,'/','Graphs','/',outputFileName]);           %Saving curve to graph subfolder as jpg
close;

%Pressure
pressureCurve = plot(pressureTable.(2),pressureTable.(1));
    %xticks(0:0.2:100);
    %yticks(-100:100:10000);
grid on
title(['Pressure Curve - ','Test #',testNumber,', ',dateOfFire,''])
xlabel('Time (s)')
ylabel('Pressure (psi)')
outputFileName = ['PressureCurve_Test#',testNumber,'_(',dateOfFire,')','.jpg'];
saveas(gcf,[pwd '/',subFolderName,'/','Graphs','/',outputFileName]);           %Saving curve to graph subfolder as jpg
close;

%Temperature
temperatureCurve = plot(temperatureTable.(2),temperatureTable.(1));
    %xticks(0:0.2:100);
    %yticks(-200:200:10000);
grid on
title(['Temperature Curve - ','Test #',testNumber,', ',dateOfFire,''])
xlabel('Time (s)')
ylabel('Temperature (F)')
outputFileName = ['TemperatureCurve_Test#',testNumber,'_(',dateOfFire,')','.jpg'];
saveas(gcf,[pwd '/',subFolderName,'/','Graphs','/',outputFileName]);           %Saving curve to graph subfolder as jpg
close;

%CF vs Time
subplot(2,1,1)
CFCurve = plot(pressureTable.(2),thrustCoefficents);
grid on
title(['Thrust Coefficent Curve - ','Test #',testNumber,', ',dateOfFire,''])
xlabel('Time (s)')
ylabel('CF')

%CF vs p1/p2
subplot(2,1,2)
semilogx(pressureTable.(1)./14.6959,thrustCoefficents)
grid on

title(['Thrust Coefficent vs Pressure Ratio (chamber/ambient) - ','Test #',testNumber,', ',dateOfFire,''])
xlabel('p1 / p2')
ylabel('CF')

outputFileName = ['CFCurve_Test#',testNumber,'_(',dateOfFire,')','.jpg'];
saveas(gcf,[pwd '/',subFolderName,'/','Graphs','/',outputFileName]);           %Saving curve to graph subfolder as jpg
close;

%Thrust and pressure
plot(thrustTable.(2),thrustTable.(1));
grid on
title(['Thrust and Pressure Curve - ','Test #',testNumber,', ',dateOfFire,''])
xlabel('Time (s)')

yyaxis left
ylabel('Thrust (N)')
ylim([-15 1.10*max(thrustTable.(1))])

yyaxis right
plot(pressureTable.(2),pressureTable.(1));
ylabel('Pressure (psi)')
ylim([0 1.10*max(pressureTable.(1))])

legend("Thrust (N)", "Pressure (psi)", "location", "best")

str = sprintf('Max Thrust: %.1f N \nMax Pressure: %.1f psi\nImpulse: %.1f N-s \nIsp: %.1f s \nc*: %.1f m/s \nAvg. CF: %.3f',maxThrust,maxPressure,impulse,isp,c_Star,thrustCoefficent);
text(0.05, 0.95*max(pressureTable.(1)), str)
% %max labels
% [maxPressure, maxPressureIndex] = max(pressureTable.(1));
% [maxThrust, maxThrustIndex] = max(thrustTable.(1));
% str = sprintf('\\leftarrow Max Pressure = %.0f psi',maxPressure);
% text(maxPressureIndex/1000,maxPressure*1.03,str);
% %xline(maxPressureIndex/1000);
% str = sprintf('Max Thrust = %.0f N \\rightarrow',maxThrust);
% text(maxThrustIndex/1000,maxPressure*1.03,str,'HorizontalAlignment','right');
% %xline(maxThrustIndex/1000);

outputFileName = ['Thrust_Pressure_Test#',testNumber,'_(',dateOfFire,')','.jpg'];
saveas(gcf,[pwd '/',subFolderName,'/','Graphs','/',outputFileName]);           %Saving curve to graph subfolder as jpg
close;

% %burn rate polynomial
% scatter(xData_burnrate_original, yData_burnrate_original);
% hold on
% plot(fitresult_poly, xData_burnrate, yData_burnrate);
% grid on
% title(['Burn Rate Versus Pressure Polynomial - ','Test #',testNumber,', ',dateOfFire,''])
% xlabel('Pressure (psi)')
% ylabel('Burn Rate (in/s)')
% 
% outputFileName = ['BurnRate_poly_Test#',testNumber,'_(',dateOfFire,')','.jpg'];
% saveas(gcf,[pwd '/',subFolderName,'/','Graphs','/',outputFileName]);           %Saving curve to graph subfolder as jpg
% close;
% 
% %burn rate exponential
% scatter(xData_burnrate_original, yData_burnrate_original);
% hold on
% plot(fitresult_exp, xData_burnrate, yData_burnrate);
% grid on
% title(['Burn Rate Versus Pressure Exponential - ','Test #',testNumber,', ',dateOfFire,''])
% xlabel('Pressure (psi)')
% ylabel('Burn Rate (in/s)')
% 
% outputFileName = ['BurnRate_exp_Test#',testNumber,'_(',dateOfFire,')','.jpg'];
% saveas(gcf,[pwd '/',subFolderName,'/','Graphs','/',outputFileName]);           %Saving curve to graph subfolder as jpg
% close;

%Subplot of 3 graphs
% Thrust Curve
subplot(3,1,1)
thrustCurve = plot(thrustTable.(2),thrustTable.(1));
grid on
title(['Thrust Curve - ','Test #',testNumber,', ',dateOfFire,''])
xlabel('Time (s)')
ylabel('Thrust (N)')

% Pressure Curve
subplot(3,1,2)
pressureCurve = plot(pressureTable.(2),pressureTable.(1));
grid on
title(['Pressure Curve - ','Test #',testNumber,', ',dateOfFire,''])
xlabel('Time (s)')
ylabel('Pressure (psi)')

% Temperature Curve
subplot(3,1,3)
temperatureCurve = plot(temperatureTable.(2),temperatureTable.(1));
grid on
title(['Temperature Curve - ','Test #',testNumber,', ',dateOfFire,''])
xlabel('Time (s)')
ylabel('Temperature (F)')

%Saving subplot to subfolder as jpg
outputFileName = ['Graphs_Test#',testNumber,'_(',dateOfFire,')','.jpg'];
saveas(gcf,[pwd '/',subFolderName,'/',outputFileName]);         %gcf means the current figure. pwd means current working directory       

%% ui graph
fig2 = figure(2);
clf(fig2);
fig2;

%thrust and pressure data
tab1 = uitab('Title','Data');
axData1 = axes(tab1);
grid on
title(['Thrust and Pressure Curve - ','Test #',testNumber,', ',dateOfFire,''])
xlabel(axData1,'Time (s)')
hold on
yyaxis left
plot(axData1,thrustTable.(2),thrustTable.(1));
ylabel(axData1,'Thrust (N)')
ylim(axData1,[-15 1.10*max(thrustTable.(1))])
yyaxis right
plot(axData1,pressureTable.(2),pressureTable.(1));
ylabel(axData1,'Pressure (psi)')
ylim(axData1,[0 1.10*max(pressureTable.(1))])
legend("Thrust (N)", "Pressure (psi)", "location", "best")
%info
str = sprintf('Impulse: %.1f N-s \nIsp: %.1f s \nc*: %.1f m/s \nAvg. CF: %.1f',impulse,isp,c_Star,thrustCoefficent);
text(0.05, 1.00*max(pressureTable.(1)), str)
%max labels
%[maxPressure, maxPressureIndex] = max(pressureTable.(1));
%[maxThrust, maxThrustIndex] = max(thrustTable.(1));
str = sprintf('\\leftarrow Max Pressure = %.0f psi',maxPressure);
text(maxPressureIndex/1000,maxPressure*1.03,str);
%xline(maxPressureIndex/1000);
str = sprintf('Max Thrust = %.0f N \\rightarrow',maxThrust);
text(maxThrustIndex/1000,maxPressure*1.03,str,'HorizontalAlignment','right');
%xline(maxThrustIndex/1000);

% % burn rate data
% tab2 = uitab('Title','Burn Rate');
% axBurnRate1 = axes(tab2);
% scatter(axBurnRate1,xData_burnrate_original, yData_burnrate_original);
% hold on
% plot(fitresult_exp, xData_burnrate, yData_burnrate);
% grid on
% title(['Burn Rate Versus Pressure Exponential - ','Test #',testNumber,', ',dateOfFire,''])
% xlabel('Pressure (psi)')
% ylabel('Burn Rate (in/s)')

%CF vs Time
tab3 = uitab('Title','CF');
axCF = axes(tab3);
subplot(2,1,1)
CFCurve = plot(pressureTable.(2),thrustCoefficents);
grid on
title(['Thrust Coefficent Curve - ','Test #',testNumber,', ',dateOfFire,''])
xlabel('Time (s)')
ylabel('CF')
%CF vs p1/p2
subplot(2,1,2)
semilogx(pressureTable.(1)./14.6959,thrustCoefficents)
grid on
title(['Thrust Coefficent vs Pressure Ratio (chamber/ambient) - ','Test #',testNumber,', ',dateOfFire,''])
xlabel('p1 / p2')
ylabel('CF')

%Subplot of 3 graphs
% Thrust Curve
tab4 = uitab('Title','all');
axAll = axes(tab4);
subplot(3,1,1)
thrustCurve = plot(thrustTable.(2),thrustTable.(1));
grid on
title(['Thrust Curve - ','Test #',testNumber,', ',dateOfFire,''])
xlabel('Time (s)')
ylabel('Thrust (N)')
% Pressure Curve
subplot(3,1,2)
pressureCurve = plot(pressureTable.(2),pressureTable.(1));
grid on
title(['Pressure Curve - ','Test #',testNumber,', ',dateOfFire,''])
xlabel('Time (s)')
ylabel('Pressure (psi)')
% Temperature Curve
subplot(3,1,3)
temperatureCurve = plot(temperatureTable.(2),temperatureTable.(1));
grid on
title(['Temperature Curve - ','Test #',testNumber,', ',dateOfFire,''])
xlabel('Time (s)')
ylabel('Temperature (F)')


%% Outputting Data



%Outputing cut thrust curve to excel file for use in trajectory calc/simulink 
T = [thrustTable.(2),thrustTable.(1)];   %  time,thrust
outputFileNameThrust = ['Cut_ThrustData',testNumber,'_(',dateOfFire,')','.xlsx'];
xlswrite([pwd '/',subFolderName,'/',outputFileNameThrust],T);       %pwd return current working directory
xlswrite([pwd '/',subFolderName,'/',outputFileNameThrust],avgMassFlow, 'Sheet1', 'C1')


%Writing to txt file
fprintf(txtFileID, ['----------Rocket Motor Test Fire #',testNumber,', ',dateOfFire,'----------\n\n']);
fprintf(txtFileID, ['Ambient Temperature = ',num2str(ambientTemperature),' F\n\n']);

fprintf(txtFileID, ['     ~Propellant Properties~\n']);
fprintf(txtFileID, ['Fuel Mixture: potassium nitrate - sorbitol\n']);
fprintf(txtFileID, ['Burn Time = ',num2str(burnTime),' s \n']);
fprintf(txtFileID, ['Grain Mass = ',num2str(propellantMass*1000),' g \n']);
fprintf(txtFileID, ['Segment Lengths = ',num2str(L1),', ',num2str(L2),' in \n']);
fprintf(txtFileID, ['Grain Density = ',num2str(density*61.024),' g/cm^3 \n\n']);

fprintf(txtFileID, ['         ~Performance~\n']);
fprintf(txtFileID, ['Average mass flow = ',num2str(avgMassFlow),' kg/s \n']);
fprintf(txtFileID, ['Impulse = ',num2str(impulse),' N-s \n']);
fprintf(txtFileID, ['Isp = ',num2str(isp),' s \n']);
fprintf(txtFileID, ['c* = ',num2str(c_Star),' m/s \n']);
fprintf(txtFileID, ['Avg. CF = ',num2str(thrustCoefficent),' \n\n']);

fprintf(txtFileID, ['         ~Peak Values~\n']);
fprintf(txtFileID, ['Max Thrust= ',num2str(maxThrust),' N \n   at t= ',num2str(thrustTable.(2)(maxThrustIndex)),' s']);          %%note column and rows are switched when calling table  .(column)(row)
fprintf(txtFileID, ['\nMax Pressure= ',num2str(maxPressure),' psi \n   at t= ',num2str(pressureTable.(2)(maxPressureIndex)),' s']);
fprintf(txtFileID, ['\nMax Temeprature= ',num2str(maxTemperature),' F \n   at t= ',num2str(temperatureTable.(2)(maxTemperatureIndex)),' s\n\n']);

fprintf(txtFileID, ['           ~Misc.~\n']);
fprintf(txtFileID, ['Burn start time after DAQ start: ',num2str(burnStartTime),' s \n']);
fprintf(txtFileID, ['Burn end time after DAQ start: ',num2str(burnEndTime),' s \n']);
fprintf(txtFileID, ['Excel Log Output:\n']);
fprintf(txtFileID, [num2str(testNumber),',',num2str(dateOfFire),',KNSB,-,',num2str(ambientTemperature),',',num2str(isp),',',num2str(c_Star),',',num2str(thrustCoefficent),',',num2str(impulse),',',num2str(burnTime),',',num2str(avgMassFlow),',',num2str(maxThrust),',',num2str(maxPressure),',',num2str(maxTemperature),',',num2str(propellantMass*1000),',',num2str(L1),',',num2str(L2),',',num2str(OD),',',num2str(ID),',-',',',num2str(density*61.024),'\n']);
fprintf(txtFileID, ['Note: ',note,'\n']);

fprintf(txtFileID, ['\n\nMatthew Stevens & Jacob Curtiss']);
fprintf(txtFileID, ['\n----------------------------------------------------------']);

% %%burn rate file
% fprintf(txtFileID_burnRate, ['----------Rocket Motor Test Fire #',testNumber,', ',dateOfFire,'----------\n']);
% fprintf(txtFileID_burnRate, ['Output: burnrate (in/s)\n']);
% fprintf(txtFileID_burnRate, ['input: pressure (psi)\n\n']);

% %print exponential fit
% coeffs = coeffvalues(fitresult_exp);
% coefficientNames = coeffnames(fitresult_exp);
% fit_formula = formula(fitresult_exp); 
% fprintf(txtFileID_burnRate,['Exponential fit: ']);
% fprintf(txtFileID_burnRate,['    ',fit_formula,'\n']);
% i = 1;
% for c = coeffs
%     fprintf(txtFileID_burnRate,['      ',coefficientNames{i},': ',num2str(c),'\n']);
%     i = i + 1;
% end

% fprintf(txtFileID_burnRate,'\n\n');
% 
% % print polynomial fit
% coeffs = coeffvalues(fitresult_poly);
% coefficientNames = coeffnames(fitresult_poly);
% fit_formula = formula(fitresult_poly); 
% fprintf(txtFileID_burnRate,['Polynomial fit: ']);
% fprintf(txtFileID_burnRate,['    ',fit_formula,'\n']);
% i = 1;
% for c = coeffs
%     fprintf(txtFileID_burnRate,['      ',coefficientNames{i},': ',num2str(c),'\n']);
%     i = i + 1;
% end

%% Final Remarks
disp(['Output files saved under subfolder: ', subFolderName]);
fclose('all');                                                  %closing subfolders