function generate_trim_report(mass)
%% 计算条件设定

if nargin<1 
    mass = 315;
end

h0     = 150;
gamma0 = 0;
v0     = 20:1:50;

rpt_name = 'trimming_report';

%% 模型初始化
proj = slproject.getCurrentProject;
wd = cd;
try 
    cd([proj.RootFolder '/docs/trim/',num2str(mass),'kg']);
catch
    mkdir([proj.RootFolder '/docs/trim/',num2str(mass),'kg']);
    cd([proj.RootFolder '/docs/trim/',num2str(mass),'kg']);
end
mdl = evalin('base', 'mdl_name');
fcn = [mdl, '_an'];
opts              =  trimoptions(fcn);
opts.v0_type      = uint8(1); % 0-TAS;1-CAS;2-Mach
opts.max_alpha    = 12/57.3;
opts.min_alpha    = -3/57.3;
opts.max_de       = 30/57.3;
opts.min_de       = -30/57.3;
opts.max_da       = 20/57.3;
opts.min_da       = -20/57.3;
opts.max_dr       = 10/57.3; 
opts.min_dr       = -10/57.3;
opts.max_df       = 30/57.3;
opts.min_df       = 0/57.3;
opts.max_throttle = 1;
opts.min_throttle = 0.0;
opts.throttle     = 0.5;
opts.max_beta     = 9/57.3;
opts.min_beta     = -9/57.3;
opts.TrimOptions.OptimizationOptions.Display = 'iter';
opts.TrimOptions.DisplayReport = 'on';
% opts.OptimizationOptions.MaxFunEvals = 600;
% opts.OptimizationOptions.TolX = 1e-4;
load_system(fcn);
mdlWks = get_param(fcn,'ModelWorkspace');

assignin(mdlWks,'mass0',mass);
%% 初始化报告
import mlreportgen.report.*
import mlreportgen.dom.*
% rpt_name = 'trimming_report';
rpt = Report(rpt_name,'docx');

tp = TitlePage; 
tp.Title = 'Trimming Report'; 
tp.Subtitle = ['for model ' fcn]; 
tp.Author = 'Lyu Weiye'; 
add(rpt,tp);
%% 计算
 
% trim_lon_alpha_df_report(rpt,opts,mdlWks,h0,v0,gamma0);

trim_lon_report(rpt,opts,mdlWks,h0,v0,gamma0);

% trim_lon_cg_report(rpt,opts,mdlWks,h0,v0,gamma0);
%  
% linLonLat_report(rpt,opts,mdlWks,h0,v0,gamma0);
%  
% trimGamma_report(rpt,opts,mdlWks,h0,v0);
%  
% trimSideslip_report(rpt,opts,mdlWks,h0,v0,gamma0);
%  
% trimTurnPhiDr_report(rpt,opts,mdlWks,h0,v0);
% 
% trimPullup_report(rpt,opts,mdlWks,h0,v0);

% trim_lon_df_report(rpt,opts,mdlWks,h0,v0,gamma0);

% throttle = 1;
% trimGamma_df_report(rpt,opts,mdlWks,h0,v0,throttle);

%% 结束
rptview(rpt);
close_system(fcn, 0); % close without saving
cd(wd);
end

%% 纵向配平
function trim_lon_report(rpt,opts,~,h0,v0,gamma0)
    hf = figure();
    set(gcf,'Position',[0 0 700 600]);
    plotTrimLon(hf, ...
        opts, ...
        h0, v0, ...
        gamma0, 'CAS/ m/s');

    fig = mlreportgen.report.Figure();
    fig.Snapshot.Caption = sprintf('定直平飞配平结果(H=%.0f)',h0);
    add(rpt,fig);
    savefig(sprintf('trimLon_h%.0f',h0));
    delete(gcf);
end

function plotTrimLon(hf, opts, h0, v0, gamma0, xlbl)
    if nargin < 6
        xlbl = 'v / m/s';
    end
    j=1;
    for v = v0
        [success,opts1] = trimLon(opts, h0, v, gamma0);
        if success
            opts = opts1;
            v_lst(j)   = v;
            aoa_lst(j) = opts.alpha;
            de_lst(j)  = opts.de;
            thr_lst(j) = opts.throttle;
            thrust_lst(j) = opts.thrust; 
            LoD_lst(j)  = opts.LoD;   
            RPM_lst(j) = opts.RPM; 
            power_lst(j) = opts.power; 
            j=j+1;
        end
    end
    
    if j>1
        figure(hf);
        subplot(421);
        plot(v_lst, aoa_lst*57.3);
        xlabel(xlbl);
        ylabel('\alpha / ^{o}');
        hold on;
        grid on;
        
        subplot(422);
        plot(v_lst, LoD_lst);
        xlabel(xlbl);
        ylabel('LoD');
        hold on;
        grid on;
        
        subplot(423);
        plot(v_lst, thrust_lst);
        xlabel(xlbl);
        ylabel('T_r / N');
        hold on;
        grid on;
        
        subplot(424);
        plot(v_lst, de_lst*57.3);
        xlabel(xlbl);
        ylabel('\delta_e / ^{o}');
        hold on;
        grid on;
        
        subplot(425);
        plot(v_lst, thr_lst);
        xlabel(xlbl);
        ylabel('\delta_T');
        hold on;
        grid on;
        
        subplot(426);
        plot(v_lst, power_lst);
        xlabel(xlbl);
        ylabel('P_r / W');
        hold on;
        grid on;

        subplot(427);
        plot(v_lst, RPM_lst);
        xlabel(xlbl);
        ylabel('RPM');
        hold on;
        grid on;


    end
end

%% 不同重心纵向配平
function trim_lon_cg_report(rpt,opts,mdlWks,h0,v0,gamma0)
   
    hf = figure();
    set(gcf,'Position',[0 0 700 600]);
    cg_lst = [0,0.000541,0.000974];
    for idx = 1:length(cg_lst)
        cg = cg_lst(idx);
        lbl{idx} = ['cgx=',num2str(cg),'m'];
%         assignin(mdlWks,'force_CG',1); 
%         assignin(mdlWks,'force_CG_pos_m',[cg, 0, 0]); CG_bias_m
        assignin(mdlWks,'CG_bias_m',[cg, 0, 0]); 
        plotTrimLon(hf, ...
            opts, ...
            h0, v0, ...
            gamma0, 'CAS/ m/s');
    end
    
    assignin(mdlWks,'force_CG',0);
    assignin(mdlWks,'force_CG_pos_m',[0, 0, 0]);
    assignin(mdlWks,'CG_bias_m',[0, 0, 0]);
    lgd = legend(lbl,'location','best');
    set(lgd,...
    'Position',[0.720333850982767 0.162299541026901 0.152500001464571 0.125950573061809]);
    fig = mlreportgen.report.Figure();
    fig.Snapshot.Caption = sprintf('定直平飞配平结果_前后重心(H=%.0f)',h0);
    add(rpt,fig);
    savefig(sprintf('trimLon_cg_h%.0f',h0));
    delete(gcf);
end

%% 模态分析
function linLonLat_report(rpt,opts,~,h0,v0,gamma0)
  
    hf = figure();
    plotLinLonLat(hf, ...
        opts, ...
        h0, v0, ...
        gamma0);
    fig = mlreportgen.report.Figure();
    fig.Snapshot.Caption = sprintf('线性化模态结果(H=%.0f)',h0);
    add(rpt,fig);
    savefig(sprintf('linearize_modes%.0f',h0));
    delete(gcf);
end

%% 爬升
function trimGamma_report(rpt,opts,~,h0,v0)
    hf = figure();
    set(gcf,'Position',[0 0 700 600]);
    for thr=[0, 0.7, 0.8, 0.9, 1]
    plotTrimGamma(hf, ...
        opts, ...
        h0, v0, ...
        thr, 'CAS/ m/s');
    end

fig = mlreportgen.report.Figure();
fig.Snapshot.Caption = sprintf('定常爬升下滑配平结果(H=%.0f)',h0);
add(rpt,fig);
savefig(sprintf('climb_glide%.0f',h0));
delete(gcf);
end

function plotTrimGamma(hf, opts, h0, v0, thr0, xlbl)
    if nargin < 6
        xlbl = 'v / m/s';
    end
    j=1;
    for v =v0
        [success,opts1] = trimGamma( opts, h0, v, thr0 );
        if success && opts1.gamma > -0.2
            opts = opts1;
            v_lst(j) = v;
            aoa_lst(j) = opts.alpha;
            de_lst(j) = opts.de;
            thr_lst(j) = opts.throttle;           
            gamma_lst(j) = opts.gamma;
            hdot_lst(j) = sin(opts.gamma)*opts.TAS;
            thrust_lst(j) = opts.thrust;
            RPM_lst(j) = opts.RPM;
            j=j+1;
        end
    end
  
    if j>1
        figure(hf);
        subplot(421);
        plot(v_lst, aoa_lst*57.3);
        xlabel(xlbl);
        ylabel('\alpha / ^{o}');
        hold on;
        grid on;
        
        subplot(422);
        plot(v_lst, de_lst*57.3);
        xlabel(xlbl);
        ylabel('\delta_e / ^{o}');
        hold on;
        grid on;
        
        subplot(423);
        plot(v_lst, gamma_lst*57.3);
        xlabel(xlbl);
        ylabel('\gamma / ^{o}');
        hold on;
        grid on;
        
        subplot(424);
        plot(v_lst, hdot_lst);
        xlabel(xlbl);
        ylabel('RoC / m/s');
        hold on;
        grid on;
        
        subplot(425);
        plot(v_lst, thrust_lst);
        xlabel(xlbl);
        ylabel('T / N');
        hold on;
        grid on;
        
        subplot(426);
        plot(v_lst, thr_lst);
        xlabel(xlbl);
        ylabel('\delta_T');
        hold on;
        grid on;

        subplot(427);
        plot(v_lst, RPM_lst);
        xlabel(xlbl);
        ylabel('RPM');
        hold on;
        grid on;

    end
end

%% 带侧滑飞行配平
function trimSideslip_report(rpt,opts,~,h0,v0,gamma0)
  
    hf = figure();
    set(gcf,'Position',[0 0 700 600]);
    for beta0 = [-5:5:5]./57.3
    plotTrimSideslip(hf, ...
        opts, ...
        h0, v0, ...
        beta0, gamma0, 'CAS/ m/s');
    end

fig = mlreportgen.report.Figure();
fig.Snapshot.Caption = sprintf('带侧滑飞行配平结果(H=%.0f)',h0);
add(rpt,fig);
savefig(sprintf('trimSideslip%.0f',h0));
delete(gcf);
end

function plotTrimSideslip(hf, opts, h0, v0, beta0, gamma0, xlbl)
    if nargin < 7
        xlbl = 'v / m/s';
    end
    j=1;
    for v =v0
        [success,opts1] = trimSideslip(opts, h0, v, beta0, gamma0);
        if success 
            opts = opts1;
            v_lst(j) = v;
            aoa_lst(j) = opts.alpha;
            aos_lst(j) = opts.beta;
            da_lst(j) = opts.da;
            de_lst(j) = opts.de;
            del_lst(j) = opts.del;
            der_lst(j) = opts.der;
            dr_lst(j) = opts.dr;
            phi_lst(j) = opts.phi;
            r_lst(j) = opts.r;
            thr_lst(j) = opts.throttle;
            thrust_lst(j) = opts.thrust;
%             RPM_lst(j) = opts.rpm;
            j=j+1;
        end
    end

    if ~isempty(v_lst)
        figure(hf);
        subplot(421);
        plot(v_lst, aoa_lst*57.3);
        xlabel(xlbl);
        ylabel('\alpha / ^{o}');
        hold on;
        grid on;
        
        subplot(422);
        plot(v_lst, de_lst*57.3);
        xlabel(xlbl);
        ylabel('\delta_e / ^{o}');
        hold on;
        grid on;
        
        subplot(423);
        plot(v_lst, aos_lst*57.3);
        xlabel(xlbl);
        ylabel('\beta / ^{o}');
        hold on;
        grid on;
        
        subplot(424);
        plot(v_lst, dr_lst*57.3);
        xlabel(xlbl);
        ylabel('\delta_r / ^{o}');
        hold on;
        grid on;
        
        subplot(425);
        plot(v_lst, da_lst*57.3);
        xlabel(xlbl);
        ylabel('\delta_a / ^{o}');
        hold on;
        grid on;
        
        subplot(426);
        plot(v_lst, phi_lst*57.3);
        xlabel(xlbl);
        ylabel('\phi/ ^{o}');
        hold on;
        grid on;   
        
        subplot(427);
        plot(v_lst, del_lst*57.3);
        xlabel(xlbl);
        ylabel('\delta_{el} / ^{o}');
        hold on;
        grid on;
        
        subplot(428);
        plot(v_lst, der_lst*57.3);
        xlabel(xlbl);
        ylabel('\delta_{er} / ^{o}');
        hold on;
        grid on;  
    end   
end
%% 盘旋特性
function trimTurnPhiDr_report(rpt,opts,~,h0,v0)
    hf = figure();
    set(gcf,'Position',[0 0 800 700]);
    idx = 1;
    for phi=-30:15:30
    plotTurnPhiDr(hf, ...
        opts, ...
        h0, v0, ...
        phi/57.3, 0, 'CAS/ m/s');
    end

fig = mlreportgen.report.Figure();
fig.Snapshot.Caption = sprintf('非协调盘旋配平结果(H=%.0f)',h0);
add(rpt,fig);
savefig(sprintf('turn_dr0%.0f',h0));
delete(gcf);
end

function plotTurnPhiDr(hf, opts, h0, v0, phi0, dr0, xlbl)
    if nargin < 7
        xlbl = 'v / m/s';
    end
    j=1;
    for v =v0
        [success,opts1] = trimTurnPhiDr( opts, h0, v, phi0, dr0);
        if success
            opts = opts1;
            v_lst(j) = v;
            aoa_lst(j) = opts.alpha;
            aos_lst(j) = opts.beta;
            de_lst(j) = opts.de;
            thr_lst(j) = opts.throttle;
            thrust_lst(j) = opts.thrust;
            Nn_lst(j) = opts.Nn;
            da_lst(j) = opts.da;
            dr_lst(j) = opts.dr;
            gamma_lst(j) = opts.gamma;
            LoD_lst(j) = opts.LoD;
            phi_lst(j) = opts.phi;
            omega_lst(j) = sqrt(dot(opts.pqr,opts.pqr)); 
            if phi0==0
                R_lst(j) = 0;
            else
                R_lst(j) = v_lst(j)/omega_lst(j);
            end
            j=j+1;
        end
    end

    if ~isempty(v_lst)
        figure(hf);
        subplot(331);
        plot(v_lst,aoa_lst*57.3);
        xlabel(xlbl);
        ylabel('\alpha / ^o');
        grid on
        hold on
        
        subplot(332);
        plot(v_lst,de_lst*57.3);
        xlabel(xlbl);
        ylabel('\delta_e / ^o');
        grid on
        hold on
        
        subplot(333);
        plot(v_lst,thr_lst);
        xlabel(xlbl);
        ylabel('\delta_T');
        grid on
        hold on
        
        subplot(334);
        plot(v_lst,da_lst*57.3);
        xlabel(xlbl);
        ylabel('\delta_a / ^o');
        grid on
        hold on
        
        subplot(335);
        plot(v_lst,dr_lst*57.3);
        xlabel(xlbl);
        ylabel('\delta_r / ^o');
        grid on
        hold on
        
        subplot(336);
        plot(v_lst,aos_lst*57.3 );
        xlabel(xlbl);
        ylabel('\beta / ^o');
        grid on
        hold on;
        
        subplot(337);
        plot(v_lst,phi_lst*57.3);
        xlabel(xlbl);
        ylabel('\phi / ^o');
        grid on
        hold on;
        
        subplot(338);
        plot(v_lst,Nn_lst);
        xlabel(xlbl);
        ylabel('N_z');
        grid on
        hold on;
        
        subplot(339);
        plot(v_lst,omega_lst*57.3 );
        xlabel(xlbl);
        ylabel('\omega / ^o/s');
        grid on
        hold on;
        
    end
end

%% 机动特性
function trimPullup_report(rpt,opts,~,h0,v0)
    hf = figure();
    set(gcf,'Position',[0 0 700 600]);
    idx = 1;
    for Nz = -0.5:0.5:2.5
    plotTrimPullup(hf, ...
        opts, ...
        h0, v0, ...
        Nz, 'CAS/ m/s');
        idx = idx+1;
    end

    fig = mlreportgen.report.Figure();
    fig.Snapshot.Caption = sprintf('定常拉起配平结果(H=%.0f)',h0);
    add(rpt,fig);
    savefig(sprintf('pullup_h%.0f',h0));
    delete(gcf);
end

function plotTrimPullup(hf, opts, h0, v0, Nz0, xlbl)
    if nargin < 6
        xlbl = 'v / m/s';
    end
    j=1;
    v_lst = [];
    for v =v0
        [success,opts1] = trimPullup( opts, h0, v, Nz0);
        if success
            opts = opts1;
            v_lst(j) = v;
            aoa_lst(j) = opts.alpha;
            de_lst(j) = opts.de;
            thrust_lst(j) = opts.thrust;
            thr_lst(j) = opts.throttle;            
            nz_lst(j) = opts.Nn;
            j=j+1;
        end
    end
    if ~isempty(v_lst)
        figure(hf);
        subplot(321);
        plot(v_lst,aoa_lst*57.3);
        xlabel(xlbl);
        ylabel('\alpha / ^o');
        grid on
        hold on
        
        subplot(322);
        plot(v_lst,de_lst*57.3);
        xlabel(xlbl);
        ylabel('\delta_e / ^o');
        grid on
        hold on
        
        subplot(323);
        plot(v_lst, thrust_lst);
        xlabel(xlbl);
        ylabel('T / N');
        grid on
        hold on
        
        
        subplot(324);
        plot(v_lst, thr_lst);
        xlabel(xlbl);
        ylabel('\delta_T');
        grid on
        hold on;  
        
        subplot(325);
        plot(v_lst, nz_lst);
        xlabel(xlbl);
        ylabel('Nz');
        grid on
        hold on;
   
    
    end   
end

%% 单发失效配平
function trim_engfault_report(rpt,opts,mdlWks,h0,v0,gamma0)
    hf = figure();
    set(gcf,'Position',[0 0 700 600]);
    for idx = 2
        assignin(mdlWks,'eng_fault',idx); % right engine fault
        plotTrimEngfault(hf, ...
            opts, ...
            h0, v0, ...
            gamma0);
    end
    assignin(mdlWks,'eng_fault',0); 
fig = mlreportgen.report.Figure();
fig.Snapshot.Caption = sprintf('trim_Engfault_rlst(H=%.0f)',h0);
add(rpt,fig);
savefig(sprintf('trimEngfault_h%.0f',h0));
delete(gcf);
end

function plotTrimEngfault(hf, opts, h0, v0, gamma0, xlbl)
    if nargin < 6
        xlbl = 'v / m/s';
    end
    j=1;
    for v = v0
        [success,opts1] = trimEngfault(opts, h0, v, gamma0);
        if success
            opts = opts1;
            v_lst(j)   = v;
            aoa_lst(j) = opts.alpha;
            aos_lst(j) = opts.beta;
            de_lst(j)  = opts.de;
            dal_lst(j) = opts.dal;
            dar_lst(j) = opts.dar;
            dr_lst(j)  = opts.dr;
            da_lst(j)  = opts.da;
            thr_lst(j) = opts.throttle;
%             thrust_lst(j) = opts.thrust; % 单发推力(N)
            LoD_lst(j)  = opts.LoD;           
%             cft_lst(j)  = opts.fuelrate; % 单发燃油消耗率(kg/s)
%             RPM_lst(j)  = opts.rpm;
            j=j+1;
        end
    end
    
    figure(hf);
    subplot(421);
    plot(v_lst, aoa_lst*57.3);
    xlabel(xlbl);
    ylabel('\alpha / ^{o}');
    hold on;
    grid on
    subplot(422);
    plot(v_lst, aos_lst*57.3);
    xlabel(xlbl);
    ylabel('\beta / ^{o}');
    hold on;
    grid on
    subplot(423);
    plot(v_lst, thr_lst);
    xlabel(xlbl);
    ylabel('thr');
    hold on;
    grid on
    subplot(424);
    plot(v_lst, da_lst*57.3);
    xlabel(xlbl);
    ylabel('\delta_a / ^{o}');
    hold on;
    grid on
    subplot(425);
    plot(v_lst, de_lst*57.3);
    xlabel(xlbl);
    ylabel('\delta_e / ^{o}');
    hold on;
    grid on
    subplot(426);
    plot(v_lst, dr_lst*57.3);
    xlabel(xlbl);
    ylabel('\delta_r / ^{o}');
    hold on;
    grid on
    
    subplot(427);
    plot(v_lst, del_lst*57.3);
    xlabel(xlbl);
    ylabel('\delta_{el} / ^{o}');
    hold on;
    grid on;
    
    subplot(428);
    plot(v_lst, der_lst*57.3);
    xlabel(xlbl);
    ylabel('\delta_{er} / ^{o}');
    hold on;
    grid on;
end


%% 最大油门单发失效配平
function trim_maxengfault_report(rpt,opts,mdlWks,h0,v0)
    hf = figure();
    set(gcf,'Position',[0 0 800 700]);
    for idx = 0:1:2
        assignin(mdlWks,'eng_fault',idx); % right engine fault
        plotTrimMaxEngfault(hf, ...
            opts, ...
            h0, v0);
    end
    assignin(mdlWks,'eng_fault',0); 
fig = mlreportgen.report.Figure();
fig.Snapshot.Caption = sprintf('trim_maxEngfault_rlst(H=%.0f)',h0);
add(rpt,fig);
savefig(sprintf('trimMaxEngfault_h%.0f',h0));
delete(gcf);
end

function plotTrimMaxEngfault(hf, opts, h0, v0, xlbl)
    if nargin < 6
        xlbl = 'v / m/s';
    end
    j=1;
    for v = v0
        [success,opts1] = trimEngfaultThr(opts, h0, v, 1);
        if success
            opts = opts1;
            v_lst(j)   = v;
            aoa_lst(j) = opts.alpha;
            aos_lst(j) = opts.beta;
            de_lst(j)  = opts.de;
            dal_lst(j) = opts.dal;
            dar_lst(j) = opts.dar;
            dr_lst(j)  = opts.dr;
            da_lst(j)  = opts.da;
            thr_lst(j) = opts.throttle;
            gamma_lst(j) = opts.gamma;
            hdot_lst(j) = sin(opts.gamma)*opts.TAS;
%             thrust_lst(j) = opts.thrust; % 单发推力(N)
            LoD_lst(j)  = opts.LoD;           
%             cft_lst(j)  = opts.fuelrate; % 单发燃油消耗率(kg/s)
%             RPM_lst(j)  = opts.rpm;
            j=j+1;
        end
    end
    
    figure(hf);
    subplot(331);
    plot(v_lst, aoa_lst*57.3);
    xlabel(xlbl);
    ylabel('\alpha / ^{o}');
    hold on;
    grid on
    subplot(332);
    plot(v_lst, aos_lst*57.3);
    xlabel(xlbl);
    ylabel('\beta / ^{o}');
    hold on;
    grid on
    subplot(333);
    plot(v_lst, hdot_lst);
    xlabel(xlbl);
    ylabel('RoC / m/s');
    hold on;
    grid on
    subplot(334);
    plot(v_lst, da_lst*57.3);
    xlabel(xlbl);
    ylabel('\delta_a / ^{o}');
    hold on;
    grid on
    subplot(335);
    plot(v_lst, de_lst*57.3);
    xlabel(xlbl);
    ylabel('\delta_e / ^{o}');
    hold on;
    grid on
    subplot(336);
    plot(v_lst, dr_lst*57.3);
    xlabel(xlbl);
    ylabel('\delta_r / ^{o}');
    hold on;
    grid on
    
    subplot(337);
    plot(v_lst, dal_lst*57.3);
    xlabel(xlbl);
    ylabel('\delta_{al} / ^{o}');
    hold on;
    grid on;
    
    subplot(338);
    plot(v_lst, dar_lst*57.3);
    xlabel(xlbl);
    ylabel('\delta_{ar} / ^{o}');
    hold on;
    grid on;
end


function plotLinLonLat(hf, opts, h0, v0, gamma0)
    figure(hf);
    for v =v0
        [success,opts1,op_point] = trimLon( opts, h0, v, gamma0 );
        if success
            opts = opts1;
%             [Glon,Glat] = linfdm_numericalpert(opts.fcn,op_point);
            [Glon,Glat] = linfdm(opts.fcn,op_point);
            ev = eig(Glon);
            subplot(211);
            plot(real(ev),imag(ev), 'x');
            xlabel('Real');
            ylabel('Imag');
            hold on;
            ev = eig(Glat);
            subplot(212);
            plot(real(ev),imag(ev), 'x');
            xlabel('Real');
            ylabel('Imag');
            hold on;
        end
    end
end


%% 襟翼构型配平

function trim_lon_df_report(rpt,opts,mdlWks,h0,v0,gamma0)

    hf = figure();
    set(gcf,'Position',[0 0 700 800]);
    df_lst = [0,   10,  20,  30,  40];
    lbl = {'df=0','df=10','df=20','df=30','df=40'};
%     df_lst = [-5, 0,   5];
%     lbl = {'df=-5', 'df=0','df=5'};
    for idx = 1:length(df_lst)
        df = df_lst(idx);
%         lbl{idx} = ['df=',num2str(df),'^o'];
        assignin(mdlWks,'df',df./57.3); 
        plotTrimLon(hf, ...
            opts, ...
            h0, v0, ...
            gamma0, 'CAS/ m/s');
    end
    assignin(mdlWks,'df',0./57.3); 
    lgd = legend(lbl,'location','best');
    set(lgd,...
    'Position',[0.720333850982767 0.162299541026901 0.152500001464571 0.125950573061809]);
    fig = mlreportgen.report.Figure();
    fig.Snapshot.Caption = sprintf('定直平飞配平结果_变襟翼(H=%.0f)',h0);
    add(rpt,fig);
    savefig(sprintf('trimLon_df_h%.0f',h0));
    delete(gcf);
end


%% 纵向配平
function trim_lon_alpha_df_report(rpt,opts,~,h0,v0,gamma0)
    hf = figure();
    alpha0 = 4/57.3;
    set(gcf,'Position',[0 0 700 600]);
    plotTrimLonALphaDf(hf, ...
        opts, ...
        h0, v0, alpha0,...
        gamma0, 'CAS/ m/s');

    fig = mlreportgen.report.Figure();
    fig.Snapshot.Caption = sprintf('定直平飞配平结果(H=%.0f)',h0);
    add(rpt,fig);
    savefig(sprintf('trimLon_h%.0f',h0));
    delete(gcf);
end

function plotTrimLonALphaDf(hf, opts, h0, v0, alpha0, gamma0, xlbl)
    if nargin < 6
        xlbl = 'v / m/s';
    end
    j=1;
    for v = v0
        [success,opts1] = trimLonAlphaDf(opts, h0, alpha0, gamma0,v);
        if success
            opts = opts1;
            v_lst(j)   = v;
            df_lst(j)  = opts.df;
            aoa_lst(j) = opts.alpha;
            de_lst(j)  = opts.de;
            thr_lst(j) = opts.throttle;
            thrust_lst(j) = opts.thrust; 
            LoD_lst(j)  = opts.LoD;   
            RPM_lst(j) = opts.RPM; 
            power_lst(j) = opts.power; 
            j=j+1;
        end
    end
    
    if j>1
        figure(hf);
        subplot(321);
        plot(v_lst, aoa_lst*57.3);
        xlabel(xlbl);
        ylabel('\alpha / ^{o}');
        hold on;
        grid on;
        
        subplot(322);
        plot(v_lst, LoD_lst);
        xlabel(xlbl);
        ylabel('LoD');
        hold on;
        grid on;
        
        subplot(323);
        plot(v_lst, df_lst*57.3);
        xlabel(xlbl);
        ylabel('\delta_f / ^{o}');
        hold on;
        grid on;
        
        subplot(324);
        plot(v_lst, de_lst*57.3);
        xlabel(xlbl);
        ylabel('\delta_e / ^{o}');
        hold on;
        grid on;
        
        subplot(325);
        plot(v_lst, thr_lst);
        xlabel(xlbl);
        ylabel('\delta_T');
        hold on;
        grid on;
        
        subplot(326);
        plot(v_lst, power_lst);
        xlabel(xlbl);
        ylabel('P_r / W');
        hold on;
        grid on;
        
    end
end


function trimGamma_df_report(rpt,opts,mdlWks,h0,v0,throttle)

    hf = figure();
    set(gcf,'Position',[0 0 700 600]);
    df_lst = [0,   10,  20,  30];
    lbl = {'df=0','df=10','df=20','df=30'};
    for idx = 1:length(df_lst)
        df = df_lst(idx);
%         lbl{idx} = ['df=',num2str(df),'^o'];
        assignin(mdlWks,'df',df./57.3); 
        plotTrimGammaDf(hf, ...
            opts, ...
            h0, v0, ...
            throttle, 'CAS/ m/s');
    end
    assignin(mdlWks,'df',0./57.3); 
    lgd = legend(lbl,'location','best');
    set(lgd,...
    'Position',[0.720333850982767 0.162299541026901 0.152500001464571 0.125950573061809]);
    fig = mlreportgen.report.Figure();
    fig.Snapshot.Caption = sprintf('定爬升配平结果_变襟翼(H=%.0f)',h0);
    add(rpt,fig);
    savefig(sprintf('trimGamma_df_h%.0f',h0));
    delete(gcf);
end


function plotTrimGammaDf(hf, opts, h0, v0, thr0, xlbl)
    if nargin < 6
        xlbl = 'v / m/s';
    end
    j=1;
    for v =v0
        [success,opts1] = trimGamma( opts, h0, v, thr0 );
        if success && opts1.gamma > -0.2
            opts = opts1;
            v_lst(j) = v;
            aoa_lst(j) = opts.alpha;
            de_lst(j) = opts.de;
            thr_lst(j) = opts.throttle;           
            gamma_lst(j) = opts.gamma;
            hdot_lst(j) = sin(opts.gamma)*opts.TAS;
            thrust_lst(j) = opts.thrust;
            RPM_lst(j) = opts.RPM;
            tht_lst(j) = opts.attitude(2);
            j=j+1;
        end
    end
  
    if j>1
        figure(hf);
        subplot(321);
        plot(v_lst, aoa_lst*57.3);
        xlabel(xlbl);
        ylabel('\alpha / ^{o}');
        hold on;
        grid on;
        
        subplot(322);
        plot(v_lst, de_lst*57.3);
        xlabel(xlbl);
        ylabel('\delta_e / ^{o}');
        hold on;
        grid on;
        
        subplot(323);
        plot(v_lst, gamma_lst*57.3);
        xlabel(xlbl);
        ylabel('\gamma / ^{o}');
        hold on;
        grid on;
        
        subplot(324);
        plot(v_lst, hdot_lst);
        xlabel(xlbl);
        ylabel('RoC / m/s');
        hold on;
        grid on;
        
        subplot(325);
        plot(v_lst, tht_lst*57.3);
        xlabel(xlbl);
        ylabel('\theta / ^{o}');
        hold on;
        grid on;
        
        subplot(326);
        plot(v_lst, thr_lst);
        xlabel(xlbl);
        ylabel('\delta_T');
        hold on;
        grid on;
    end
end