%% This code generated stable and unstable boundaries using time domain aeroelastic simulation
%% Depending on what damping model is activated in file 'ae_solver_fun_visc.m' it 
%% will plot the apppropriate boundaries

%% Params_all is a matlab file that includes all parameters identified with system identification
%% techniques

%% TIME DOMAIN SOLVER-MAIN %%
clc
% close all
clear all
global Params
global ctr; global case_no;


%% Read Run File
runs=xlsread('ae_runs_time_domain_visc.xlsx');
% This file contains the run information. The first row lists the run
% velocity and all initial conditions to be tested for that velocity are
% listed in that column.

Vs=runs(1,:);
for i=1:length(Vs)
    VV{i}.h0s=runs(2:end,i);
    VV{i}.h0s(isnan(VV{i}.h0s))=[];
end
ts=[0 10];

%% Solver
for ii=1:length(Vs)
    VV{ii}.h0s_y=nan;
    VV{ii}.h0s_n=nan;
    VV{ii}.lco_amp=nan;
        for jj=1:length(VV{ii}.h0s)
            disp(['Case   ',num2str(ii),' of ',num2str(length(Vs)), ...
                '  and run ',num2str(jj),'  of  ',num2str(length(VV{ii}.h0s)),' ...']);

            V=Vs(ii);
            h0=VV{ii}.h0s(jj);
            initial=[h0;0;0;0];
            options = odeset('AbsTol',1e-3);
            [TH,XH]=ode45(@(TH,XH) ae_solver_fun_visc(TH,XH,V),ts, initial);
            tf=TH;
            pitch_f=XH(:,3);
            plunge_f=XH(:,1);
%             lco_amp=rms(plunge_f(floor(end-1/10*end:end)));
            lco_amp=max(plunge_f(floor(end-1/10*end):end))/2-min(plunge_f(floor(end-1/10*end):end))/2;
            
            if(lco_amp>0.003), VV{ii}.h0s_y=[VV{ii}.h0s_y h0];end
            if(lco_amp<=0.003), VV{ii}.h0s_n=[VV{ii}.h0s_n h0];end
            if(lco_amp>0.003), VV{ii}.lco_amp=[VV{ii}.lco_amp lco_amp]; end
        
        end
end
%% Extract lines
for ii=1:length(Vs)
    VV{ii}.lco_amp(isnan(VV{ii}.lco_amp))=[];
    VV{ii}.lco_bounds=[min(VV{ii}.lco_amp) max(VV{ii}.lco_amp)];
    VV{ii}.lco_mean=mean(VV{ii}.lco_amp);
    VV{ii}.friction1=max(VV{ii}.h0s_n);
    VV{ii}.friction2=min(VV{ii}.h0s_y);

%     lco_max(ii)=VV{ii}.lco_bounds(2);    
%     lco_min(ii)=VV{ii}.lco_bounds(1);
    lco_mean(ii)=VV{ii}.lco_mean;
    friction1(ii)=VV{ii}.friction1;
    friction2(ii)=VV{ii}.friction2;
    VV{ii}.Vs=Vs(ii);
    
end

   
%% Pltting Bifurcation Diagram
% load data_visc0
% for i=1:length(VV)
%     Vs(i)=VV{i}.Vs;
% end

hFig1 = figure('Position', [100, 10, 1000, 600]);
% figure(1)
 for ii=1:length(Vs)
    hold on
    if VV{ii}.h0s_n~=nan, plot(Vs(ii),VV{ii}.h0s_n,'og'); end
    if VV{ii}.h0s_y~=nan, plot(Vs(ii),VV{ii}.h0s_y,'or','MarkerFaceColor','r');end
    if VV{ii}.lco_amp~=nan, plot(Vs(ii),VV{ii}.lco_amp,'sk','MarkerFaceColor','k');end
 end

l(1)= plot(Vs(ii),VV{ii}.h0s_n(1),'og');
l(2)= plot(Vs(ii),VV{ii}.h0s_y(1),'or','MarkerFaceColor','r');
% l(3)=plot(Vs(ii),VV{ii}.lco_amp(1),'sk','MarkerFaceColor','k');

l(3)=plot(Vs,lco_mean,'-sk','MarkerFaceColor','k');
% hold on
% plot(Vs,friction1,'b.');
grid on
% plot(Vs(ii),VV{ii}.friction2,'xb');
% load flutter_boundary_data

% l(4)=plot(flutter_boundary_data(1,:),flutter_boundary_data(2,:),':k','LineWidth',3);
xlim([7 16]);
ylim([0 0.03])

set(gca,'FontSize',12)
xlabel('Freestream Velocity (m/s)','FontSize',12);
ylabel('Plunge Amplitude (m)','FontSize',12);

% flutter_speed_data=[10.57;0];
% l(5)=plot(flutter_speed_data(1,:),flutter_speed_data(2,:),'ob','MarkerSize',8,'MarkerFaceColor','b');
% legend(l,'Initial Conditions: No LCO','Initial Conditions: LCO','LCo Amplitude','Flutter boundary: Time domain','Flutter Boundary: p-method');


%% Plot Boundaries
% clear l
% hFig1 = figure('Position', [100, 10, 1000, 600]);
% hold on
% flutter_speed_data=[10.57;0];
% l(1)=plot(flutter_speed_data(1,:),flutter_speed_data(2,:),'ob','MarkerSize',8,'MarkerFaceColor','b');
% % load flutter_boundary_data
% % l(2)=plot(flutter_boundary_data(1,:),flutter_boundary_data(2,:),':k','LineWidth',3);
% % load visc_boundary1_data
% % l(3)=plot(visc_boundary1_data(1,:),visc_boundary1_data(2,:),'k-^','LineWidth',2);
% % load visc_boundary2_data
% % l(4)=plot(visc_boundary2_data(1,:),visc_boundary2_data(2,:),'k-s','LineWidth',2);
% grid on
% xlim([7 16]);
% ylim([0 0.03])
% set(gca,'FontSize',12)
% xlabel('Freestream Velocity (m/s)','FontSize',12);
% ylabel('Plunge Amplitude (m)','FontSize',12);
% % legend(l,'Flutter boundary: p-method','Flutter Boundary: Time domain', ...
% % 'LCO amplitude: I_\alpha=0.0567','LCO amplitude: I_\alpha=0.0607');



