function bht_controls_check(filename,range,tol,varargin)
% BHT_CONTROLS_CHECK
%
% Allow to check controls M-Files.
%
% Syntax
%
% BHT_CONTROLS_CHECK(filename,t_range,tol)
%
% BHT_CONTROLS_CHECK(...,cont_parameters)
%
% Description
%
% A valid controls M-File contains functions that give the time
% evolving joints' angles of the articulated bodies (called the
% controls functions) but it contains also the expressions of the first
% and second derivatives of these functions. Any error in these
% expressions will impact deeply the bodies' motions.
%
% The function BHT_CONTROLS_CHECK computes numeric first and second
% derivatives of the controls functions given in the M-File filename,
% over the time interval t_range = [T(1),T(2)] and with theoretical
% accuracy tol, using the central difference formula. A figure is next
% displayed with three subplots. The first one is the plot of the
% controls functions. The others are the plots of the errors between
% the numeric derivatives and user's expressions of these derivatives
% given in the M-File. We point out that an error larger than the
% theoretical one does not necessarily means that the formula is wrong.
% Conversely, there can be an error in the formula even if the computed
% error is within the theoretical one.
%
% Add the optional input 'cont_parameters' when the controls M-File
% requires a parameters array.
%
% See also BHT_DATA_COMPILE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% BIOHYDRODYNAMICS MATLAB TOOLBOX %
% %
% A. MUNNIER and B. PINCON %
% %
% alexandre.munnier@iecn.u-nancy.fr bruno.pincon@iecn.u-nancy.fr %
% http://www.iecn.u-nancy.fr/~munnier http://www.iecn.u-nancy.fr/~pincon %
% %
% %
% INSTITUT ELIE CARTAN, NANCY 1 %
% http://www.iecn.u-nancy.fr/ %
% %
% INRIA Lorraine, Projet CORIDA %
% http://www.iecn.u-nancy.fr/~corida/ %
% %
% %
% %
% %
% August 15th 2008 %
% %
% GNU GPL v3 licence %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%##########################################################################
% Reading options
if ~isempty(varargin)
cont_parameters = varargin{1};
else
cont_parameters = [];
end
%-----------------------
if numel(varargin)>1 && strcmpi(varargin{2},'latex')
latex = 1;
else
latex = 0;
end
%==========================================================================
% Computing number of points for the discretization
T0 = range(1);
T1 = range(2);
nb_points = floor((T1-T0)/sqrt(tol));
t = linspace(T0,T1,nb_points+1);
t2 = linspace(T0,T1,2*nb_points+1);
dt = (T1-T0)/nb_points;
dt2 = 0.5*(T1-T0)/nb_points;
%------------------------------------------
% Evaluating the controls
for k = 1:nb_points+1
[c(:,k),dc(:,k),dttc(:,k)] = feval(filename,t(k),cont_parameters); %#ok<AGROW>
end
for k =1:2*nb_points+1
[c2(:,k),dc2(:,k)] = feval(filename,t2(k),cont_parameters); %#ok<AGROW>
end
%------------------
nb_var = size(c,1);
%==========================================================================
% preallocating numeric derivatives
dc_num = zeros(nb_var,nb_points-1);
dc2_num = zeros(nb_var,nb_points-1);
dc1 = zeros(nb_var,nb_points-1);
dttc_num = zeros(nb_var,nb_points-1);
dttc2_num = zeros(nb_var,nb_points-1);
dttc1 = zeros(nb_var,nb_points-1);
%-----------------------------------
% Computing numeric deriatives
for k = 2:nb_points
dc_num(:,k-1) = 0.5*(c(:,k+1)-c(:,k-1))./dt;
dc2_num(:,k-1) = 0.5*(c2(:,2*k)-c2(:,2*k-2))./dt2;
dc1(:,k-1) = dc(:,k);
dttc_num(:,k-1) = 0.5*(dc(:,k+1)-dc(:,k-1))./dt;
dttc2_num(:,k-1) = 0.5*(dc2(:,2*k)-dc2(:,2*k-2))./dt2;
dttc1(:,k-1) = dttc(:,k);
end
%-----------------------------------
% Computing errors
dc_err = abs(dc1-dc_num);
dc2_err = abs(dc1-dc2_num);
dttc_err = abs(dttc1-dttc_num);
dttc2_err = abs(dttc1-dttc2_num);
%------------------------------------
% Displaying conclusions
if (max(max(4*dc2_err-dc_err))<tol) %#ok<BDSCI>
mess1 = ' -> Ok';
else
mess1 = ' -> Something''s wrong';
end
if (max(max(4*dttc2_err-dttc_err))<tol) %#ok<BDSCI>
mess2 = ' -> Ok';
else
mess2 = ' -> Something''s wrong';
end
disp(['Error order estimated: ',num2str(tol)])
disp(['Max error for dc/dt: ',num2str(max(dc_err,[],2)'),mess1])
disp(['Max error for d^2c/dt^2: ',num2str(max(dttc_err,[],2)'),mess2])
%##########################################################################
% figure
%##########################################################################
t_num = t;
t_num(1) = [];
t_num(end) = [];
% building legends
legend0 = cell(nb_var,1);
legend1 = cell(nb_var,1);
legend2 = cell(nb_var,1);
for k = 1:nb_var
if latex
legend0{k} = ['$c_',num2str(k),'$'];
legend1{k} = ['${dc_',num2str(k),'}\,/{dt}$'];
legend2{k} = ['${d^2c_',num2str(k),'}\,/{dt^2}$'];
else
legend0{k} = ['c(',num2str(k),')'];
legend1{k} = ['dc(',num2str(k),')/dt'];
legend2{k} = ['d^2c(',num2str(k),')/dt'];
end
end
%==========================================================================
% figure properties
figure(1)
set(1,'visible','off','Position',[132 61 1152 748],...
'Name',['Controls M-File: ',filename] ,'numbertitle','off')
%==========================================================================
% plot 1
subplot(3,1,1)
plot(t,c)
% -------------------
% title
if latex
title('The controls $c$','fontsize',16,'interpreter','latex')
else
title('The controls c','fontsize',16)
end
%-----------------------------------
% axis
set(gca,'fontsize',14,'layer','top')
grid on
% ----------------------------------
% legend
h = legend(legend0);
set(h,'fontsize',14,'location','eastoutside')
if latex
set(h,'interpreter','latex')
end
%==========================================================================
% plot 2
subplot(3,1,2)
plot(t_num,dc_err)
%----------------------------
% title
if latex
title('Error for ${dc}\,/{dt}$','fontsize',16,'interpreter','latex')
else
title('Error for dc/dt','fontsize',16)
end
%-----------------------------------
% axis
set(gca,'fontsize',14,'layer','top')
grid on
%-----------------------------------
% legend
h = legend(legend1);
set(h,'fontsize',14,'location','eastoutside')
if latex
set(h,'interpreter','latex')
end
%==========================================================================
% plot 3
subplot(3,1,3)
plot(t_num,dttc_err)
%----------------------------------
% title
if latex
title('Error for ${d^2c}\,/{dt^2}$','fontsize',16,'interpreter','latex')
else
title('Error for d^2c/dt^2','fontsize',16)
end
%-----------------------------------
% axis
set(gca,'fontsize',14,'layer','top')
grid on
%----------------------------------
% legend
h = legend(legend0);
set(h,'fontsize',14,'location','eastoutside')
if latex
set(h,'interpreter','latex')
end
%==========================================================================
set(1,'visible','on')
%==========================================================================
end
