close all
clear
clc

f=readcell('50degfull_VLM16cg');

AOA=[];
beta=[];
CL=[];
CD=[];
CM=[];
CX=[];
CY=[];
CZ=[];
nptsBeta=17;
nptsAl=5;
endpt=43; %11 for wing, 43 for full aircraft 

for i=1:1:size(f,1)
   if f(i,1)=="FC_AoA_"
       AOA(end+1)=f{i,2};
   elseif f(i,1)== "FC_Beta_"
       beta(end+1)=f{i,2};
   elseif f(i,1)=="cd"
            CD(end+1)=sum(double(cell2sym(f(i,2:endpt))));
       elseif f(i,1)=="cl"
            CL(end+1)=sum(double(cell2sym(f(i,2:endpt))));
        elseif f(i,1)=="cmy"
            CM(end+1)=sum(double(cell2sym(f(i,2:endpt))));
       elseif f(i,1)=="cmx"
            CX(end+1)=sum(double(cell2sym(f(i,2:endpt))));       
       elseif f(i,1)=="cy"
            CY(end+1)=sum(double(cell2sym(f(i,2:endpt))));
       elseif f(i,1)=="cmz"
            CZ(end+1)=sum(double(cell2sym(f(i,2:endpt))));
   end
end
alpha=unique(AOA);
beta =unique(beta);

%% plot 1
titles='CL-\alpha plot';
xlabels= '\alpha';
ylabels='CL';
variable_name='CL Alpha';
plot_at_beta(nptsBeta,CL,alpha,beta,titles,xlabels,ylabels,variable_name)

%% plot 2
titles='CD-\alpha plot';
xlabels= '\alpha';
ylabels='CD';
variable_name='CD Alpha';
plot_at_beta(nptsBeta,CD,alpha,beta,titles,xlabels,ylabels,variable_name)

%% plot 3
titles='C_M-\alpha plot';
xlabels= '\alpha';
ylabels='C_M';
variable_name='C_M Alpha';
plot_at_beta(nptsBeta,CM,alpha,beta,titles,xlabels,ylabels,variable_name)

%% plot 4
titles='C_L-\beta plot [Roling Moment coefficient]';
xlabels= '\beta';
ylabels='C_L';
variable_name='C_L Beta ';
plot_at_alpha(nptsBeta,nptsAl,CX,alpha,beta,titles,xlabels,ylabels,variable_name)
    %% plot 5
titles='C_Y-\beta plot [Body axis Force Coefficient]';
xlabels= '\beta';
ylabels='C_Y';
variable_name='C_Y beta ';
plot_at_alpha(nptsBeta,nptsAl,CY,alpha,beta,titles,xlabels,ylabels,variable_name)
    %% plot 6
titles='C_n-\beta plot [ Yawing Moment]';
xlabels= '\beta';
ylabels='C_n';
variable_name='C_n Beta';
plot_at_alpha(nptsBeta,nptsAl,CZ,alpha,beta,titles,xlabels,ylabels,variable_name)
%% defined functions
function [] = plot_at_alpha(nptsBeta,nptsAl,variCoeff,alpha,beta,titles,xlabels,ylabels,variable_name)
    i=1;
    j=1;
    A=figure(); 
    var_at_alpha=[];
    for i=1:1:nptsBeta
       index= i:nptsBeta:nptsAl*nptsBeta;     %1 to 85 increment of 17
       CL_gar=variCoeff(index);
       var_at_alpha=[var_at_alpha CL_gar];
    end
    i=1;
    while i<=length(var_at_alpha)
        plot(beta,var_at_alpha(i:nptsAl+i-1))   
        i=i+nptsAl;
        hold on
        lgnd(j)=(strcat(' \alpha: ',string(alpha(j))));
        j=j+1;
    end
    legend(lgnd)
    xlabel(xlabels)
    ylabel(ylabels)
    title(titles)
    grid on
    grid minor
    variable_name = strcat(variable_name,'(50degfull_VLM16cg)');
    saveas(A,strcat(variable_name,'.png'));
end

function [] = plot_at_beta(nptsBeta,variCoeff,alpha,beta,titles,xlabels,ylabels,variable_name)
    i=1;
    j=1;
    A=figure();
    
    while i<=length(variCoeff)
        plot(alpha,variCoeff(i:nptsBeta+i-1))    
        i=i+nptsBeta;
        hold on
        lgnd(j)=(strcat(' \beta: ',string(beta(j))));
        j=j+1;
    end
    legend(lgnd)
    xlabel(xlabels)
    ylabel(ylabels)
    title(titles)
    grid on
    grid minor
    variable_name = strcat(variable_name,'(50degfull_VLM16cg)');
    saveas(A,strcat(variable_name,'.png'));
end