% Section2.m
%
% Generate figures for Section 2 of paper
% "The Effect of Rotor Tilt on the Stability and Gust Rejection Properties
%  of VTOL Multirotor Aircraft" 
% by James Whidborne Arthur Mendez and Alastair Cooke 

% (c) James Whidborne
% Cranfield University, 16 August 2020
 
disp('2 Rotor Aerodynamic Modelling')
disp('-----------------------------')
disp(' ')

% Get rotor coefficents for Draganflyer XPro quadrotor
[~, ~, rotor, rho] = Coefficients;

%% Figure 3 - Thrust mapping for constant airspeeds.

% define grids & get map
Omega = 80:5:240;           % rotor speeds
alpha = [-pi/2:pi/32:pi/2]; % disc AoA
V = 0:10:20;                % airspeed (m/s)

[T,H,Q] = ThrustDragMap(Omega,alpha,V,rotor,rho); % get map

% plots
figure(3) % Thrust mapping for constant airspeeds
for k = 1:length(V)
    Tk = reshape(T(k,:,:),[length(Omega),length(alpha)]); % extract thrust map for each V
    mesh(180*alpha/pi,Omega,Tk,'linewidth',lw)
    hold on
end
hold off
h99=gca; xlim(h99,[-90 90]); % change axes limits
h99.XTick=[-90 -60 -30 0 30 60 90];
h99.XMinorTick='on';
ylim(h99,[80 240]); % change axes limits
h99.YMinorTick='on';
zlim(h99,[-20 60]); % change axes limits
h99.ZTick=[-20 -10 0 10 20 30 40 50 60];
h99.ZTick=[-20 0 20 40 60];
h99.ZMinorTick='on';
grid on

% Create textarrows
h1=gcf;
annotation(h1,'textarrow',[0.20  0.23 ],...
    [0.36  0.41 ],'String',{'$V=0$'},...
    'Interpreter' ,'Latex',...
    'HeadLength',6,'HeadWidth',6);

annotation(h1,'textarrow',[0.17  0.15 ],...
    [0.43  0.6],'String',{'$V=10$'},...
    'Interpreter' ,'Latex',...
    'HeadLength',6,'HeadWidth',6);

annotation(h1,'textarrow',[0.40  0.31 ],...
    [0.77 0.70 ],'String',{'$V=20$'},...
    'Interpreter' ,'Latex',...
    'HeadLength',6,'HeadWidth',6);

% label axes
xlabel('$\alpha$ (deg)','Interpreter' ,'Latex','fontsize',fs)
ylabel('$\Omega$ (rad/s)','Interpreter' ,'Latex','fontsize',fs)
zlabel('$T$ (N)','Interpreter' ,'Latex','fontsize',fs);
title(['$V \in \{$',num2str(V(1:end)),'$ \}$'],'Interpreter' ,'Latex' )

% save
if(printflag) save_pdf(gcf,'TMap'); end


%% Figure 4 - Horizontal force mapping for constant airspeeds

figure(4)
for k = 1:length(V)
    Hk = reshape(H(k,:,:),[length(Omega),length(alpha)]);% extract force map for each V
    mesh(180*alpha/pi,Omega,Hk)
    hold on
end
hold off
h99=gca; xlim(h99,[-90 90]); % change axes limits
h99.XTick=[-90 -60 -30 0 30 60 90];
h99.XMinorTick='on';
ylim(h99,[80 240]); % change axes limits
h99.YMinorTick='on';
zlim(h99,[0 35]); % change axes limits
h99.ZTick=[0:5:35];
h99.ZMinorTick='on';
grid on

% Create textarrows
figure1=gcf;
annotation(figure1,'textarrow',[0.63957142857143 0.54957142857143],...
    [0.64447619047619 0.57447619047619],'String',{'$V=20$'},...
    'Interpreter','latex',...
    'HeadWidth',6,...
    'HeadLength',6);

% Create textarrow
annotation(figure1,'textarrow',[0.710714285714286 0.598214285714286],...
    [0.523809523809524 0.226190476190476],'String',{'$V=10$'},...
    'Interpreter','latex',...
    'HeadWidth',6,...
    'HeadLength',6);

% Create textarrow
annotation(figure1,'textarrow',[0.751785714285714 0.614285714285714],...
    [0.414285714285714 0.159523809523809],'String',{'$V=0$'},...
    'Interpreter','latex',...
    'HeadWidth',6,...
    'HeadLength',6);

% label axes
xlabel('$\alpha$ (deg)','Interpreter' ,'Latex','fontsize',fs)
ylabel('$\Omega$ (rad/s)','Interpreter' ,'Latex','fontsize',fs)
zlabel('$H$ (N)','Interpreter' ,'Latex','fontsize',fs);
title(['$V = \{$',num2str(V(1:end)),'$ \}$'],'Interpreter' ,'Latex' )

% save
if(printflag) save_pdf(gcf,'HMap'); end


%% Figure 5 - Torque mapping for constant airspeeds
figure(5)
for k = 1:length(V)
    Qk = reshape(Q(k,:,:),[length(Omega),length(alpha)]);% extract torque map for each V
    mesh(180*alpha/pi,Omega,Qk)
    hold on
end
hold off
h99=gca; xlim(h99,[-90 90]); % change axes limits
h99.XTick=[-90 -60 -30 0 30 60 90];
h99.XMinorTick='on';
ylim(h99,[80 240]); % change axes limits
h99.YMinorTick='on';
% zlim(h99,[-3 3]); % change axes limits
%h99.ZTick=[-5:10:30];
h99.ZMinorTick='on';
grid on 
% Create textarrows
figure1=gcf;
annotation(figure1,'textarrow',[0.564285714285714 0.508928571428572],...
    [0.214285714285714 0.245238095238095],'String',{'$V=0$'},...
    'Interpreter','latex',...
    'HeadWidth',6,...
    'HeadLength',6);

% Create textarrow
annotation(figure1,'textarrow',[0.308928571428572 0.260714285714286],...
    [0.769047619047619 0.716666666666667],'String',{'$V=20$'},...
    'Interpreter','latex',...
    'HeadWidth',6,...
    'HeadLength',6);

% Create textarrow
annotation(figure1,'textarrow',[0.182142857142857 0.158928571428571],...
    [0.40952380952381 0.619047619047619],'String',{'$V=10$'},...
    'Interpreter','latex',...
    'HeadWidth',6,...
    'HeadLength',6);

% label axes
xlabel('$\alpha$ (deg)','Interpreter' ,'Latex','fontsize',fs)
ylabel('$\Omega$ (rad/s)','Interpreter' ,'Latex','fontsize',fs)
zlabel('$Q$ (Nm)','Interpreter' ,'Latex','fontsize',fs);
title(['$V = \{$',num2str(V(1:end)),'$ \}$'],'Interpreter' ,'Latex' ) 

% save
if(printflag) save_pdf(gcf,'QMap'); end

%% Figure 6 Variation of thrust with airspeed  for constant rotorspeed
Omega = 150;                % rotor speeds
alpha = [-pi/2:pi/32:pi/2]; % disc AoA
V = 0:1:20;                 % airspeed (m/s)
[T,~] = ThrustDragMap(Omega,alpha,V,rotor,rho); % get map

figure(6) % Variation of thrust with airspeed for constant rotorspeed Omega=150
j=find(Omega==150);
Tj=reshape(T(:,j,:),[length(V),length(alpha)]);
mesh(180*alpha/pi,V,Tj)
hold on
contour(180*alpha/pi,V,Tj)
hold off
h99=gca; h99.XTick=[-90 -60 -30 0 30 60 90];
h99.XMinorTick='on'; h99.YMinorTick='on'; h99.ZMinorTick='on';
grid on
%axis(gca,'tight');
ylabel('$V$ (m/s)','Interpreter' ,'Latex','fontsize',fs)
xlabel('$\alpha$ (rad)','Interpreter' ,'Latex','fontsize',fs)
zlabel('$T$ (N)','Interpreter' ,'Latex','fontsize',fs);
title(['$\Omega=$',num2str(Omega(j)),' (rad/s)'],'Interpreter' ,'Latex' )
view(h99,[-63.6 40.2]);
% save
if(printflag) save_pdf(gcf,'ThrustConstantOmega'); end



%% Figure NOT USED -  Variation of thrust with airspeed  for constant rotorspeed and alpha
% define grids & get map

Omega=150;           % rotor speeds
alpha=0;             % disc AoA
V=0:1:30;            % airspeed (m/s)
[T,H] = ThrustDragMap(Omega,alpha,V,rotor,rho); % get map

figure(107) % Variation of thrust with airspeed for constant rotorspeed Omega=150
plot( V, T )
h99=gca; %h99.XTick=[-90 -60 -30 0 30 60 90];
h99.XMinorTick='on'; h99.YMinorTick='on';
%h99.ZMinorTick='on';
grid on
xlabel('$V$ (m/s)','Interpreter' ,'Latex','fontsize',fs)
ylabel('$T$ (N)','Interpreter' ,'Latex','fontsize',fs);
title(['$\Omega=$',num2str(Omega(j)),' (rad/s), $\alpha=$',num2str(alpha),' '],'Interpreter' ,'Latex' )
% save
if(printflag) save_pdf(gcf,'ThrustAirspeed'); end