% GETCDMAP
% 
% [dH,beta,CD]=getCDmap(show)
% 
% if show='plot' then the curve will be plotted
%
% Create data to populate  Simulink  Look-up Table (n-D) block to map
% $(\delta_H,\beta) \to C_D$
% The look-up table data is generated using the 'fit' surface fitting tool 
% from the MATLAB Curve Fitting Toolbox using the data from  
%   Jason Barden, “Active aerodynamic control of heavy goods vehicles,” 
%   PhD Thesis, Cranfield University, 2013
% and using piecewise linear interpolation. 
% Four additional points were introduced at the domain vertices to extend 
% the curve to the full domain. 
% Linear interpolation is used rather than higher order schemes to 
% avoid the introduction of local minima.
%
% 
% Code is used by truck_simulation_UKACC2018.m which is the 
% simulation codes for 
%   "Extremum Seeking Control for Truck Drag Reduction"
%   G Papageorgiou, J Barden, J Whidborne, K Garry
%  12th UKACC International Conference on Control 
%  Sheffield, UK, 5th - 7th September 2018

% (c) JF Whidborne, Cranfield University, 12 June 2018
 

function [dH,beta,CD]=getCDmap(show)

 
if(nargin==0), show='noplot'; end


%% scanned-in data from Jason Barden thesis Figure 7.4
% (replicated data from g1335 commented out)
Fig4_7_data=[0.8689080459770114, 0.6446884272997033
    0.8898275862068965, 0.62973293768546
    0.9119540229885057, 0.6108605341246292
    0.9147701149425287, 0.6119287833827893
    0.9167816091954022, 0.6048071216617211
    0.9228160919540229, 0.5998219584569734
    %  0.9360919540229884, 0.5980415430267063
    0.9377011494252874, 0.5930563798219585
    %  0.9417241379310344, 0.590919881305638
    0.9429310344827586, 0.5887833827893175
    %  0.9477586206896551, 0.5880712166172107
    %  0.9537931034482758, 0.5809495548961425
    %  0.9594252873563218, 0.5759643916913947
    0.9658620689655172, 0.570979228486647
    %  0.9654597701149424, 0.5699109792284867
    %  0.9710919540229884, 0.5588724035608309
    0.97203125, 0.5561206896551725
    %  0.9767241379310344, 0.5489020771513353
    %  0.9827586206896551, 0.5410682492581603
    %  0.9883908045977011, 0.54
    %  0.9940229885057471, 0.538219584569733
    0.99890625, 0.5449137931034483
    1.00203125, 0.5550431034482759
    1.008125, 0.5580603448275863
    1.0140625, 0.56;
    1.0197701149425287, 0.5610089020771514
    1.0270114942528736, 0.5638575667655786
    1.0330459770114941, 0.5649258160237389
    1.0386781609195401, 0.5670623145400594
    1.0447126436781609, 0.5699109792284867
    1.0507471264367816, 0.570979228486647
    1.1078735632183907, 0.6147774480712167
    1.1348275862068964, 0.6251038575667657
    1.1557471264367816, 0.6279525222551929
    1.1758620689655173, 0.6322255192878339];



dH0=Fig4_7_data(:,1)';
CD0=Fig4_7_data(:,2)';
beta0=0*dH0;


%% now compare with g1335 data
%dH1335= [0.9360 0.9420 0.9480 0.9540 0.9590    0.9650    0.9710    0.9770    0.9830    0.9880 0.9940];
%CD1335=[ 0.598	0.591	0.588	0.581	0.576	0.57	0.559	0.549	0.541	0.54	0.538];
%dH=[dH0,dH1335];
%CD=[CD0,CD1335];
%plot(dH0,CD0,'--*',dH1335,CD1335,'-^')

%%  include g1335 data
%deflector.CD=xlsread('CD_data.xlsx','g1335' );
deflector.CD=[
    1.032	1.027	1.028	1.023	1.018	1.013	1.013	1.01	1.007	1.004	1.005
    0.898	0.894	0.892	0.889	0.885	0.88	0.883	0.883	0.882	0.882	0.884
    0.709	0.704	0.701	0.695	0.693	0.688	0.687	0.691	0.691	0.691	0.693
    0.672	0.667	0.665	0.657	0.655	0.65	0.646	0.647	0.646	0.648	0.65
    0.64	0.634	0.631	0.625	0.62	0.615	0.608	0.605	0.604	0.605	0.605
    0.615	0.607	0.605	0.598	0.594	0.588	0.579	0.571	0.567	0.567	0.568
    0.6	    0.593	0.59	0.583	0.578	0.571	0.561	0.544	0.544	0.54	0.544
    0.598	0.591	0.588	0.581	0.576	0.57	0.559	0.549	0.541	0.54	0.538
    0.6 	0.593	0.59	0.583	0.578	0.571	0.561	0.544	0.544	0.54	0.544
    0.615	0.607	0.605	0.598	0.594	0.588	0.579	0.571	0.567	0.567	0.568
    0.64	0.634	0.631	0.625	0.62	0.615	0.608	0.605	0.604	0.605	0.605
    0.672	0.667	0.665	0.657	0.655	0.65	0.646	0.647	0.646	0.648	0.65
    0.709	0.704	0.701	0.695	0.693	0.688	0.687	0.691	0.691	0.691	0.693
    0.898	0.894	0.892	0.889	0.885	0.88	0.883	0.883	0.882	0.882	0.884
    1.032	1.027	1.028	1.023	1.018	1.013	1.013	1.01	1.007	1.004	1.005];
deflector.beta=[-15,-10,-5,-4,-3,-2,-1,0,1,2,3,4,5,10,15];
%deflector.deltaH=[484,487,490,493,496,499,502,505,508,511,514]/500;
deflector.deltaH= [0.9360 0.9420 0.9480 0.9540 0.9590  0.9650  0.9710  0.9770  0.9830  0.9880 0.9940];

for i=1:15 % add g1335 data
    CD0=[CD0,deflector.CD(i,:)];
    dH0=[dH0,deflector.deltaH];
    beta0=[beta0,deflector.beta(i)*ones(1,11)];
end

%% surface fit
% add points at corners to improve extrapolation
% CD0=[CD0,1.25,1.25,1.08,1.08,1.1,1.1];
% beta0=[beta0, -15,15, -15,15, -15,15] ;
% dH0=[dH0,1.1758620689655173,1.1758620689655173, 1.1,1.1,0.8689080459770114,0.8689080459770114];
% [CDmap, gof] = createFitCDmap(dH0, beta0, CD0);

CD0=[CD0,1.25,1.25, 1.1,1.1];
beta0=[beta0, -15,15,  -15,15] ;
dH0=[dH0,1.1758620689655173,1.1758620689655173,  0.8689080459770114,0.8689080459770114];
[CDmap] = createFitCDmap(dH0, beta0, CD0, show);

beta=[-15:0.2:15];
dH=0.87:0.001:1.17;
CD=zeros(length(dH),length(beta));
for i=1:length(dH)
    for j=1:length(beta)
        CD(i,j)=CDmap(dH(i),beta(j));
    end
end

%% check
if(0)
    figure(11)
    surfc(dH,beta,CD')
    xlabel('$\delta_H$','Interpreter' ,'Latex','fontsize',14)
    ylabel('$\beta$','Interpreter' ,'Latex','fontsize',14)
    zlabel('$C_D$','Interpreter' ,'Latex','fontsize',14)
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [CDmap, gof] = createFitCDmap(dH0, beta0, CD0, show)
%CREATEFITCDMAP(DH0,BETA0,CD0)
%  Create a fit.
%
%  Data for 'CDdeflectorMap' fit:
%      X Input : dH0
%      Y Input : beta0
%      Z Output: CD0
%  Output:
%      fitresult : a fit object representing the fit.
%      gof : structure with goodness-of fit info.
%
%  See also FIT, CFIT, SFIT.

%  Auto-generated by MATLAB on 20-Sep-2017 10:16:50


%% Fit: 'CDdeflectorMap'.
[xData, yData, zData] = prepareSurfaceData( dH0, beta0, CD0 );

% Set up fittype and options.
fittype = 'linearinterp';
%fittype = 'cubicinterp'; % gives local minima

% Fit model to data.
[CDmap, gof] = fit( [xData, yData], zData, fittype, 'Normalize', 'on' );

if strcmp(show,'plot')
    % Plot fit with data.
    figure( 'Name', 'CDdeflectorMap' );
    h = plot( CDmap, [xData, yData], zData );
    %legend( h, 'CDdeflectorMap', 'CD0 vs. dH0, beta0', 'Location', 'NorthEast' );
    % Label axes
    xlabel('$\delta_H$','Interpreter' ,'Latex','fontsize',14)
    ylabel('$\beta$','Interpreter' ,'Latex','fontsize',14)
    zlabel('$C_D$','Interpreter' ,'Latex','fontsize',14)
    grid on
    view( -33.5, 10);
    if(0) % print
        fig = gcf;
        fig.PaperPositionMode = 'auto';
        fig_pos = fig.PaperPosition;
        fig.PaperSize = [fig_pos(3) fig_pos(4)];
        saveas(gcf,'CDMap1', 'pdf')  
    end
end
return