%% Mini-Project 3
%Parameters up here for reference
clear
clc

W = 32000; %Weight in lbs
b = 65; %ft Span
AR = 7;
S = b^2/AR; %reference area in square feet
Cd0 = 0.027;
e0 = 0.8;
numEngines = 1;
Clmax = 1.3;

%conversions
%ft/s to KTAS
%KTAS = V * 0.5924838013

% Problem 1

% Aircraft Information


% Setting up important Vectors
Mvec = 0.01:0.01:1; % range for Mach number
hvec = 1:0.1:50; % range for altitude * 1000
[M,h] = meshgrid(Mvec,hvec); % makes M and h the same size vectors
[~,~,rho,~,a] = atmosphere(h.*1000); % uses the atmosphere() function to find density rho (slugs/ft^3) and speed
% of sound a (ft/s)

% Calculations of ROC
V = M.*a; % velocity in ft/s
VKTAS = V.*0.5924838013; % velocity in knots true airspeed
q = 0.5.*rho.*V.^2; % dynamic pressure
PC = 50.*ones(size(M)); % Creates a vector for Power Code the same size as M full of 50s
[Tmaxeng,~,~] = jt8d(M,h.*1000,PC); % finds the Max Thrust (lbs) and TSFC
Tmax = numEngines.*Tmaxeng; %Max Thrust for one engine
CL = W./(q.*S); % Coefficient of lift (unitless)
CD = dragpolar(CL,M); % Coefficient of drag (unitless)
D = CD.*q.*S; % Drag Force (lbf)
T = D; % Enforcing steady, level flight
P_s = (((Tmax-T).*V)./W).*60; % Specific power in ft/min

%Cl_max vector
Clmax_mesh = Clmax.*ones(size(M)); %Creates a vector full of max Cl
Clmargin = Clmax_mesh - CL; %acceptable range of CL

Clmask = Clmargin<0;
ROCmask = P_s<0;

Clmargin(ROCmask) = NaN; % trim the Clmargin for wherever ROC is less than 0
P_s(Clmask) = NaN; % trim the P_s for stall conditions

%Graphing the contour plot for Specific Power
figure(1)
hold on
grid on
Thrust_steps = 0:500:max(Tmax);
[c,handle] = contour(M,h,P_s,Thrust_steps,'LineColor','g'); % Graphs
% countour lines of P_s on axis that are M and h with red lines
clabel(c,handle);

%Making a max_ROC vector to correspond with the hvec' vector
% max_ROC = max(P_s')
h_max_ROC = [];
M_max_ROC = [];
for i = 1:410;
    [~,d] = max(P_s(i,:));
    h_max_ROC = [h_max_ROC h(i,d)];
    M_max_ROC = [M_max_ROC M(i,d)];
end
plot(M_max_ROC,h_max_ROC,'--m')
title('Rate of Climb (ft/min)');
xlabel('Mach number');
ylabel('Altitude (1000 ft)');
axis ([0 1 0 45]);

c1 = ocontourc(Mvec',hvec',P_s,[0,0]);
% Takes the range of Mach number and the range of
% altitude and plots where Tmargin = 0 at them
hatchedcontours(c1,'b'); % Uses the hatched contours
% function to make hatch marks the color blue


c2 = ocontourc(Mvec',hvec',Clmargin,[0,0]);
%Takes range of Mach number and range of altitude and plots where CLmargin
%= 0 at them
hatchedcontours(c2,'r');
% function to make the hatch marks red




hold off

rez=600; %resolution (dpi) of final graphic. 600 is a good number.
f=gcf; %handle of figure to be exported
figpos=getpixelposition(f);
resolution=get(0,'ScreenPixelsPerInch');
set(f,'paperunits','inches','papersize',figpos(3:4)/resolution,...
'paperposition',[0 0 figpos(3:4)/resolution]);
path='C:\Users\madeg\Desktop'; %Where should it go? C:\\ etc.
name='Skymap_ROC.png'; %What should it be called?
print(f,fullfile(path,name),'-dpng',['-r',num2str(rez)],'-opengl')

%% Problem 2
clc
clear

%find L_D in gliding condition

%Parameters needed
W = 32000; %lbs
b = 65; %ft Span
AR = 7;
S = b^2/AR; %reference area in square feet

%Calculating some other things
M = linspace(0.01,0.9);
h = 10000; %ft at initial conditions
[~,~,rho,~,a] = atmosphere(h); % finds density and SoS at initial conditions
V = M.*a; % velocity in ft/s
VKTAS = V.*0.5924838013; % velocity in knots true airspeed using initial density conditions
q = 0.5.*rho.*V.^2; % dynamic pressure
CL = W./(q.*S); % Coefficient of lift (unitless)
CD = dragpolar(CL,M); % Coefficient of drag (unitless)

L_D = CL./CD; %lift to drag ratio

R = (h.*L_D)./6076.115; %Range during glide at 10,000 ft

Time_aloft = R./VKTAS.*60; %in minutes

RMask = R == max(R);
KTAS_max_range = VKTAS(RMask);
Mach_max_range = KTAS_max_range./a./0.5924838013;
TAmask = Time_aloft == max(Time_aloft);
KTAS_max_time_aloft = VKTAS(TAmask); %KTAS
Mach_mach_time_aloft = KTAS_max_time_aloft./a./0.5924838013;

figure(2)
hold on
p = plot(VKTAS,R,'LineWidth',2,'Color',[.7 0.4 .9])
p1 = plot(KTAS_max_range,max(R),'o','MarkerSize',10,'MarkerFaceColor','r')
grid on
title('Range versus Knots for glide');
xlabel('KTAS');
ylabel('Range (nautical miles)');
dt = datatip(p1,KTAS_max_range,max(R))
p1.DataTipTemplate.DataTipRows(1).Label = 'Knots for Maximum Range:';
p1.DataTipTemplate.DataTipRows(2).Label = 'Maximum Range in Nautical Miles:';

rez=600; %resolution (dpi) of final graphic. 600 is a good number.
f=gcf; %handle of figure to be exported
figpos=getpixelposition(f);
resolution=get(0,'ScreenPixelsPerInch');
set(f,'paperunits','inches','papersize',figpos(3:4)/resolution,...
'paperposition',[0 0 figpos(3:4)/resolution]);
path='C:\Users\madeg\Desktop'; %Where should it go? C:\\ etc.
name='Range_vs_KTAS.png'; %What should it be called?
print(f,fullfile(path,name),'-dpng',['-r',num2str(rez)],'-opengl')

figure(3)
hold on
p2 = plot(VKTAS,Time_aloft,'LineWidth',2,'Color',[.7 0.9 .4])
p3 = plot(KTAS_max_time_aloft,max(Time_aloft),'o','MarkerSize',10,'MarkerFaceColor','r')
grid on
title('Time aloft versus Knots for glide');
xlabel('KTAS');
ylabel('Time aloft (min)')
dt = datatip(p3,KTAS_max_time_aloft,max(R))
p3.DataTipTemplate.DataTipRows(1).Label = 'Knots for Maximum Time Aloft:';
p3.DataTipTemplate.DataTipRows(2).Label = 'Maximum Time Aloft in Minutes:';


rez=600; %resolution (dpi) of final graphic. 600 is a good number.
f=gcf; %handle of figure to be exported
figpos=getpixelposition(f);
resolution=get(0,'ScreenPixelsPerInch');
set(f,'paperunits','inches','papersize',figpos(3:4)/resolution,...
'paperposition',[0 0 figpos(3:4)/resolution]);
path='C:\Users\madeg\Desktop'; %Where should it go? C:\\ etc.
name='Time_aloft_vs_KTAS.png'; %What should it be called?
print(f,fullfile(path,name),'-dpng',['-r',num2str(rez)],'-opengl')

%% Problem 3
clc
clear

%Parameters needed
h = 10000; %ft;
W = 32000; %lbs
b = 65; %ft Span
AR = 7;
S = b^2/AR; %reference area in square feet
numEngines = 1;
Cl_max = 1.3;

%Constraints given
n_max = 3.5; %g's
VKTAS_max = 300; %maximum V in knots
V_max = 300./0.5924838013; %maximum V ft/s

%Calculations to find maximum stall limit
[~,~,rho,~,a] = atmosphere(h); % finds density and SoS at initial conditions
M_max = V_max./a; %maximum M
M = linspace(0.01,M_max); %M vector with all max_m corresponding to max V constraint
V = M.*a; % velocity in ft/s vector
VKTAS = V.*0.5924838013; % velocity in knots true airspeed vector
q = 0.5.*rho.*V.^2; % dynamic pressure vector
CL = W./(q.*S); % Coefficient of lift (unitless)
CD = dragpolar(CL,M); % Coefficient of drag (unitless)

%calculating max stall speed
n_max_stall = 1/(2.*W).*rho.*(V.^2).*S.*Cl_max; %n_max_stall

%solving for the cornering velocity
syms v_corner
max_stall_eq = n_max ==  1/(2.*W).*rho.*(v_corner.^2).*S.*Cl_max;
v_corner = double(solve(max_stall_eq, v_corner));
v_corner = v_corner(2); %using the second value only
VKTAS_corner = v_corner.*0.5924838013;

%Setting up masks for plotting
n_max_stall_plot = n_max_stall; %setting up a new vector for plotting purposes
n_max_mask = n_max_stall >= 3.5;
n_max_stall_plot(n_max_mask) = NaN;
VKTAS_mask = VKTAS < VKTAS_corner;
VKTAS_n_struct_limit_plot = VKTAS; %setting
VKTAS_n_struct_limit_plot(VKTAS_mask) = NaN;
n_struct_limit = n_max.*ones(size(VKTAS_n_struct_limit_plot)); %creates a vector with the structural limit

hold on
p1 = plot(VKTAS,n_max_stall_plot,'LineWidth',2,'Color',[.7 0.9 .4]);
p2 = plot(VKTAS_n_struct_limit_plot,n_struct_limit,'LineWidth',2,'Color',[.7 0.4 .9]);
p3 = plot(VKTAS_corner,n_max,'pentagram','MarkerSize',20,'MarkerFaceColor','r');
xline(300,'--',{'Never Exceed Speed','300 KTAS'});
dt = datatip(p3,VKTAS_corner,n_max);
p3.DataTipTemplate.DataTipRows(1).Label = 'Cornering speed: (knots):';
p3.DataTipTemplate.DataTipRows(2).Label = 'Structural load limit n:';
grid on
title('V-n diagram')
xlabel('VKTAS')
ylabel('Load factor n = L/W')
axis([0 350 0 4.5])

rez=600; %resolution (dpi) of final graphic. 600 is a good number.
f=gcf; %handle of figure to be exported
figpos=getpixelposition(f);
resolution=get(0,'ScreenPixelsPerInch');
set(f,'paperunits','inches','papersize',figpos(3:4)/resolution,...
'paperposition',[0 0 figpos(3:4)/resolution]);
path='C:\Users\madeg\Desktop'; %Where should it go? C:\\ etc.
name='V_n_diagram.png'; %What should it be called?
print(f,fullfile(path,name),'-dpng',['-r',num2str(rez)],'-opengl')

%% Problem 3 part c
% 
% %What is the level turn radius (ft) and turn rate in deg/min at cornering
% %speed?
% %level turn radius
% V = v_corner; %ft/s at cornering speed
% n = n_max; %at cornering speed
% g = 32.2; %ft/s^2
% R = V^2./(g*(sqrt(n.^2 - 1))); %ft
% R = num2str(R);
% sprintf('The level turn radius at cornering speed is %s ft',R)
% %calculating turn rate
% w = ((g.*(sqrt(n.^2-1)))./V).*60.*(180./pi);
% w = num2str(w);
% sprintf('The turn rate at cornering speed is %s degrees per minute',w)

%% Problem 3 part d

%skip calculations in part c for this but keep variables from Problem 3
%parts a and b

%curve of max load factor associated with max engine thrust

%Problem parameters restated
Cd0 = 0.0270;
b = 65; %ft
AR = 7;
e0 = 0.8;
%solving for specific parameters needed
K = 1./(pi.*e0.*AR);
q = 0.5.*rho.*VKTAS.^2; %dynamic pressure
Wing_load = W./S;
[T_max,~,~]=jt8d(M,h.*ones(size(M)),50.*ones(size(M)));
n_t_max = ((q./(K.*Wing_load)).*((T_max./W) - q .* (Cd0./Wing_load))).^0.5;


hold on
p1 = plot(VKTAS,n_max_stall_plot,'LineWidth',2,'Color',[.7 0.9 .4]);
p2 = plot(VKTAS_n_struct_limit_plot,n_struct_limit,'LineWidth',2,'Color',[.7 0.4 .9]);
xline(300,'--');
grid on
title('V-n diagram')
xlabel('VKTAS')
ylabel('Load factor n = L/W')
axis([0 350 0 4.5])

p4 = plot(VKTAS,n_t_max,'LineWidth',2,'Color',[0 0 0])
x = [0.3 0.6];
y = [0.8 0.5];
annotation('textarrow',x,y,'String','Max Thrust')


rez=600; %resolution (dpi) of final graphic. 600 is a good number.
f=gcf; %handle of figure to be exported
figpos=getpixelposition(f);
resolution=get(0,'ScreenPixelsPerInch');
set(f,'paperunits','inches','papersize',figpos(3:4)/resolution,...
'paperposition',[0 0 figpos(3:4)/resolution]);
path='C:\Users\madeg\Desktop'; %Where should it go? C:\\ etc.
name='V_n_diagram_max_engine_thrust.png'; %What should it be called?
print(f,fullfile(path,name),'-dpng',['-r',num2str(rez)],'-opengl')
