function [T_new,ET,EFl,EFi,varargout] = bht_energy(varargin) %#ok<STOUT>
% BHT_ENERGY
%
% Compute the balance of energy, the torques of the articulations and the
% power expended by the articulated bodies.
%
% Syntax
%
% [T,E_tot,E_fl,E_fi] = bht_energy('DataFilename',filename)
%
% [...,E_pot,E_def,E_net,torques,power] = bht_energy(..)
%
% [... ] = bht_energy(...,options)
%
% Description
%
% The function loads data from the MAT-File filename_results (generated
% by the function bht_traject_compute from the DAT-File filename) and
% returns:
%
% o A linearly spaced vector T of the time steps.
% o The total kinetic energy E_tot of the fluid-bodies system
% (a column vector, each row corresponding to one time step).
% o The kinetic energy of the fluid E_fl.
% o The total kinetic energy of the bodies E_fi (also a column
% vector).
%
% When all of the solids composing the bodies are not neutrally
% buoyant, the function also returns:
%
% o The potential energy E_pot of the system fluid-bodies (same
% format as above) related to the buoyant force,
%
% and when all of the bodies are not reduced to single solids, you can
% get:
%
% o The kinetic energy of bodies' inner deformations (i.e. due
% to the relative motion of the links in the body's frame) E_def
% (a column vector, each row corresponding to one time step).
% o The bodies' kinetic energy of rigid motion E_net (same format
% as above). Bodies' velocity can be decomposed into the velocity
% of deformations (due to the relative displacement of the links
% expressed in the local moving frame) and the rigid velocity
% which is the motion of the body considered as rigid in the
% reference frame. This decomposition allows to decompose the
% total kinetic energy of the bodies into the rigid kinetic
% energy and the kinetic energy of deformations. In particular,
% we have always E_fi = E_net + E_def.
% o The torques torques of the articulations (an array of size
% numbers of time step x numbers of controls).
% o The power power expended by the articulations (same format
% as torques) obtained by multiplying the torques by the angular
% velocities.
%
% Function's options
%
% The options are described in the following list:
%
% * TimeRange:
% - Expected value: array [Tmin Tmax]
% - Default value: Time range of the DAT-File
% - Description: Time interval over which the computations are
% done.
%
% * DrawGraph:
% - Expected value: 'on' or 'off'
% - Default value: 'on'
% - Description: Whether energy balances and power expended by
% the bodies are ploted at the end of the computations.
%
% * latex:
% - Expected value: 'on' or 'off'
% - Default value: 'off'
% - Description: Use latex interpreter for graphs.
%
% * SaveResults:
% - Expected value: 'on' or 'off'
% - Default value: 'on'
% - Description: Whether the results are saved in the MAT-File
% filename_energy_results at the end of the computations.
%
% * NoComputations
% - Expected value: 'on' or 'off'
% - Default value: 'off'
% - Description: switch this option to 'on' ifarticulations the
% computations have already be done and the results saved on
% disk. The file 'filename_energetic_results.mat' is loaded and
% used to build the graphs.
%
%
% Example
%
% The statement :
%
% [T,E_tot,E_fl,E_fi,E_pot,E_def,E_net,torques,power] =
% bht_energy('DataFilename',filename,'TimeRange',[50
% 55],'DrawGraph','no','SaveResults','off')
%
% will return T,E_tot,E_fl,E_fi,E_pot,E_def,E_net,torques and power
% related to the experiment described in the DAT-File filename, over
% the time interval [50 55]. No graph will be displayed and the results
% will not be saved on disk.
%
% See also BHT_TRAJECT_COMPUTE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% 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 %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ==================================
% default values for options
draw_graph = 1;
time_range = [];
save_results = 1;
latex = 0;
no_filename = 1;
computations = 1;
%==========================================================================
% reading options
%================
nb_varargin = numel(varargin);
kk = 1;
while kk <= nb_varargin
if ~ischar(varargin{kk})
error(['Unknown option ',num2str(varargin{kk})]) %#ok<ST2NM>
end
%------------------------------------------------------------------
switch lower(varargin{kk})
%--------------------------------------------------------------
case 'datafilename'
if ischar(varargin{kk+1})
filename = varargin{kk+1};
exten = regexp(filename,'\w*\.(\w*)', 'tokens');
no_filename = 0;
if isempty(exten)
filename = [filename,'.dat']; %#ok<AGROW>
end
else
error('Unexpected expression for option DataFilename.')
end
kk = kk+2;
%--------------------------------------------------------------
case 'drawgraph'
if ischar(varargin{kk+1})
if strcmpi(varargin{kk+1},'off')
draw_graph = 0;
elseif strcmpi(varargin{kk+1},'on')
draw_graph = 1;
else
error('Unexpected value for option DrawGraph. Must be ''on'' or ''off''.')
end
else
error('Unexpected value for option DrawGraph. Must be ''on'' or ''off''.')
end
kk = kk+2;
%-------------
case 'saveresults'
if ischar(varargin{kk+1})
if strcmpi(varargin{kk+1},'off')
save_results = 0;
elseif strcmpi(varargin{kk+1},'on')
save_results = 1;
else
error('Unexpected value for option SaveResults. Must be ''on'' or ''off''.')
end
else
error('Unexpected value for option SaveResults. Must be ''on'' or ''off''.')
end
kk = kk+2;
%----------------------------------------------------------
case 'timerange'
if isvector(varargin{kk+1})
if length(varargin{kk+1}) == 2
time_range = varargin{kk+1};
else
error('Unexpected value for option TimeRange. Must be a vector [Tmin Tmax].')
end
else
error('Unexpected value for option TimeRange. Must be a vector [Tmin Tmax].')
end
kk = kk+2;
%----------------------------------------------------------
case 'nocomputations'
if ischar(varargin{kk+1})
if strcmpi(varargin{kk+1},'off')
computations = 1;
elseif strcmpi(varargin{kk+1},'on')
computations = 0;
else
error('Unexpected value for option NoComputations. Must be ''on'' or ''off''.')
end
else
error('Unexpected value for option NoComputations. Must be ''on'' or ''off''.')
end
kk = kk+2;
%----------------------------------------------------------
case 'latex'
if ischar(varargin{kk+1})
if strcmpi(varargin{kk+1},'off')
latex = 0;
elseif strcmpi(varargin{kk+1},'on')
latex = 1;
else
error('Unexpected value for option Latex. Must be ''on'' or ''off''.')
end
else
error('Unexpected value for option Latex. Must be ''on'' or ''off''.')
end
kk = kk+2;
%----------------------------------------------------------
otherwise
error(['Unknown option ''',varargin{kk},''''])
end
end
% =========================================================================
% DataFilename is missing
if no_filename
filename = input('Enter data file''s name (hit Return to browse): ','s');
if isempty(filename);
filename = uigetfile('*.dat','Select data file');
end
end
%==========================================================================
base_filename = regexp(filename,'(\w*)', 'tokens');
base_filename = base_filename{1}{1};
filename1 = [base_filename,'_data'];
filename2 = [base_filename,'_results'];
filename3 = [base_filename,'_energetic_results'];
kine_filename = [base_filename,'_kine'];
% need datacompil =========================================================
load(filename1,'gene_data','bound_data','ode_data','phys_data','fish');
load(filename2,'cont_filename','cont_parameters','T','p','dp','q','dq',...
'Q','dQ','JJ');
% =========================================================================
nb_solids = gene_data.nb_solids;
nb_fishes = gene_data.nb_fishes;
nb_boundaries = gene_data.nb_boundaries;
array_fishes_solids = gene_data.array_fishes_solids;
buoyant = 1-min(phys_data.rhos == phys_data.rhof);
noonlysolids = nb_solids ~= nb_fishes;
% =========================================================================
if computations
% =========================================================================
if ~isempty(time_range)
ind_start = find(T<=time_range(1),1,'last'); %#ok<NODEF>
ind_end = find(T>=time_range(2),1,'first');
else
ind_start = 1;
ind_end = length(T);
end
nb_time_step = ind_end-ind_start+1;
%==========================================================================
rhos = phys_data.rhos;
rhof = phys_data.rhof;
Ms = phys_data.Ms;
mm = phys_data.fishes_mass;
%------------------------------------
% preallocating variables used later on ===================================
ET = zeros(nb_time_step,1);
EFi = cell(1,nb_fishes);
E_net = cell(1,nb_fishes);
Edef = cell(1,nb_fishes);
for kk = 1:nb_fishes;
EFi{kk} = zeros(nb_time_step,1);
E_net{kk} = zeros(nb_time_step,1);
Edef{kk} = zeros(nb_time_step,1);
end
E_net_tot = zeros(nb_time_step,1);
E_pot = zeros(nb_time_step,1);
EFl = zeros(nb_time_step,1);
T_new = zeros(nb_time_step,1);
torques = zeros(nb_time_step,nb_solids-nb_fishes);
power = zeros(nb_time_step,nb_solids-nb_fishes);
%-------------------------
g = cell(nb_solids,1);
gg = cell(nb_solids,1);
G = cell(nb_solids,3);
Mf = zeros(3*nb_solids,3*nb_solids);
%---------------------------
S1 = zeros(3*nb_solids,1);
S2 = zeros(3*nb_solids,1);
S3 = zeros(3*nb_solids,1);
S4 = zeros(3*nb_solids,1);
B1 = zeros(3*nb_solids,1);
%-------------------------
X = cell(1,nb_boundaries);
N = cell(1,nb_boundaries);
for kk = nb_solids+1:nb_boundaries
N{kk} = bound_data(kk).N;
X{kk} = bound_data(kk).X;
end
%--------------------------
n = zeros(1,nb_boundaries);
m = zeros(1,nb_boundaries);
dt = zeros(1,nb_boundaries);
for kk = 1:nb_boundaries
n(kk) = bound_data(kk).nbr_modes;
m(kk) = bound_data(kk).nbr_points;
dt(kk) = bound_data(kk).dt;
end
mtot = sum(m);
mc = [0,cumsum(m)];
A = zeros(mtot,mtot);
b = zeros(mtot,3*(nb_solids));
phi = cell(nb_solids,3);
toto = cell(nb_solids,3);
% =========================================================================
% recomputation of the added mass matrix ==================================
rel_time1 = 0;
compt_loc = 0;
time_range = T(ind_end)-T(ind_start);
disp('Computations in progress... ')
disp('0%------20%-------40%-------60%-------80%------100%')
for compt = ind_start:ind_end
rel_time = floor((T(compt)-T(ind_start))/time_range*50);
string = repmat('|',1,rel_time-rel_time1);
rel_time1 = rel_time;
fprintf(string)
compt_loc = compt_loc +1;
T_new(compt_loc) = T(compt);
%----------------------------------------------------------------------
MM = cell(1,nb_fishes);
MS = cell(1,nb_fishes);
for kk = 1:nb_fishes
MM{kk} = zeros(3*nb_fishes,3*nb_fishes);
MS{kk} = zeros(3*nb_solids,3*nb_solids);
MM{kk}(3*(kk-1)+1:3*kk,3*(kk-1)+1:3*kk) = diag([mm(kk),mm(kk),JJ(kk,compt)]);
indices = 3*[0, cumsum(array_fishes_solids)];
MS{kk}(indices(kk)+1:indices(kk+1),indices(kk)+1:indices(kk+1)) = ...
Ms(indices(kk)+1:indices(kk+1),indices(kk)+1:indices(kk+1));
end
%----------------------------------------------------------------------
% computing controls and kinematic ================================
pp = p(:,compt);%#ok
dpp = dp(:,compt);%#ok
[c,dc,dttc] = feval(cont_filename,T(compt),cont_parameters);
[qq,dpq,dcq,dppqdp,dpcqdc,dccqdc] = feval(kine_filename,pp,dpp,c,dc);
qq = reshape(q(:,compt),3,nb_solids);%#ok
% moving solids ===================================================
for k = 1:nb_solids
h = qq(1:2,k);
theta = qq(3,k);
R = [cos(theta),-sin(theta);...
sin(theta), cos(theta)];
X{k} = R*bound_data(k).X + h*ones(1,bound_data(k).nbr_points);
N{k} = R*bound_data(k).N;
end
% computing potentials using nystrom's method =====================
[uu,duu] = mnystrom(X,N);
% computing boundary conditions ===================================
for kk = 1:nb_solids
g{kk} = [N{kk}(1,:)./bound_data(kk).dl;...
N{kk}(2,:)./bound_data(kk).dl;...
bound_data(kk).flux];
gg{kk} = bound_data(kk).fflux;
end
% building added mass matrix ======================================
for s=1:nb_solids
for l=1:3
for jj=1:nb_solids
for kk=1:3
Mf(3*(s-1)+l,3*(jj-1)+kk) =...
0.5*rhof*(bound_data(s).dt*...
sum(bound_data(s).dl.*g{s}(l,:).*uu{s}(:,3*(jj-1)+kk)')...
+ bound_data(jj).dt*...
sum(bound_data(jj).dl.*g{jj}(kk,:).*uu{jj}(:,3*(s-1)+l)'));
end
end
end
if buoyant
E_pot(compt_loc) = E_pot(compt_loc) - ...
9.80665*(rhos(s)-rhof)*bound_data(s).dt...
*sum(bound_data(s).dl.*g{s}(1,:).*X{s}(2,:).*X{s}(1,:));
else
E_pot = [];
end
end
% =====================================================================
% Mass matrix of the system
M = Ms + Mf;
%########################################################################
ET(compt_loc) = 0.5.*dq(:,compt)'*M*dq(:,compt);%#ok
for kk = 1:nb_fishes;
EFi{kk}(compt_loc) = 0.5.*dq(:,compt)'*MS{kk}*dq(:,compt);%#ok
if noonlysolids
E_net{kk}(compt_loc) = 0.5.*dQ(:,compt)'*MM{kk}*dQ(:,compt);%#ok
end
end
EFl(compt_loc) = 0.5.*dq(:,compt)'*Mf*dq(:,compt);%#ok
if noonlysolids
% =====================================================================
% computation of the torques and powers of the fishes
for kk = 1:nb_solids
G{kk,1} = [bound_data(kk).K.*(g{kk}(2,:).^2-g{kk}(1,:).^2);...
-2*bound_data(kk).K.*g{kk}(1,:).*g{kk}(2,:);...
-bound_data(kk).K.*(g{kk}(2,:).*gg{kk}+g{kk}(3,:).*g{kk}(1,:))-g{kk}(2,:)];
G{kk,2} = [-2*bound_data(kk).K.*g{kk}(1,:).*g{kk}(2,:);...
bound_data(kk).K.*(g{kk}(1,:).^2-g{kk}(2,:).^2);...
bound_data(kk).K.*(gg{kk}.*g{kk}(1,:)-g{kk}(3,:).*g{kk}(2,:))+g{kk}(1,:)];
G{kk,3} = [-bound_data(kk).K.*(g{kk}(2,:).*gg{kk}+g{kk}(3,:).*g{kk}(1,:))-g{kk}(2,:);...
bound_data(kk).K.*(g{kk}(1,:).*gg{kk}-g{kk}(2,:).*g{kk}(3,:))+g{kk}(1,:);...
bound_data(kk).K.*(gg{kk}.^2-g{kk}(3,:).^2)+gg{kk}];
end
M1 = dpq'*M*dpq; % reduced mass matrix
% =================================================================
%pp = p(:,compt);
dpp = dp(:,compt);%#ok
% computing other terms of the ode ================================
W = dpq*dpp+dcq*dc;
W1 = reshape(W,3,nb_solids);
Z = dppqdp*dpp+2*dpcqdc*dpp+dccqdc*dc+dcq*dttc;
%Z2 = dppqdp*dpp+2*dpcqdc*dpp+dccqdc*dc+dpq*dttp;
%==================================================================
% term <dqM,W>W-1/2W^t<dqM,.>W of the ode
for s = 1:nb_solids
for l = 1:3
S1p = 0;
S2p = 0;
for k = 1:nb_solids
S1p = S1p+bound_data(k).dt*...
sum(bound_data(k).dl.*duu{k}(:,3*(s-1)+l)'.*(duu{k}*W)'.*(W1(:,k)'*g{k}));
S2p = S2p+bound_data(k).dt*...
sum(bound_data(k).dl.*uu{k}(:,3*(s-1)+l)'.*...
(W1(:,k)'*(G{k,1}*W1(1,k)+G{k,2}*W1(2,k)+G{k,3}*W1(3,k))));
end;
S1(3*(s-1)+l) = S1p;
S2(3*(s-1)+l) = S2p;
S3(3*(s-1)+l) = bound_data(s).dt*...
sum(bound_data(s).dl.*g{s}(l,:).*((duu{s}*W).^2)');
S4(3*(s-1)+l) = bound_data(s).dt*...
sum(bound_data(s).dl.*g{s}(l,:).*((W1(:,s)'*g{s}).^2));
% buoyant force ==========================================
B1(3*(s-1)+l) = 9.80665*(rhos(s)-rhof)*bound_data(s).dt*...
sum(bound_data(s).dl.*g{s}(l,:).*X{s}(2,:));
end;
end;
%==================================================================
% right hand side term
F1 = rhof*dpq'*(-S1+0.5*S3+S2+0.5*S4);
F2 = dpq'*M*Z;
B = dpq'*B1;
%+++++++++
F = F1+F2-B;
%+++++++++
%============
dttp = -M1\F;
%============
M1 = dcq'*M*dcq; % reduced mass matrix
Z = dppqdp*dpp+2*dpcqdc*dpp+dccqdc*dc+dpq*dttp;
F1 = rhof*dcq'*(-S1+0.5*S3+S2+0.5*S4);
F2 = dcq'*M*Z;
B = dcq'*B1;
F = F1+F2-B;
diaddc=diag(dc);
%+++++++++++++++++++++++++++++++++++++++++++++++++
torques(compt_loc,:) = (M1*dttc+F)';
power(compt_loc,:) = torques(compt_loc,:)*diaddc';
%+++++++++++++++++++++++++++++++++++++++++++++++++
end
end
%=================================================
if buoyant
E_pot = E_pot-min(E_pot)*ones(nb_time_step,1);
ET = ET + E_pot;
end
%=================================================
if noonlysolids
for kk = 1:nb_fishes
Edef{kk} = EFi{kk}-E_net{kk};
E_net_tot = E_net_tot + E_net{kk};
end
else
Edef = [];
E_net = [];
torques = [];
power = [];
end
%++++++++++++++++++++++
varargout{1} = E_pot;
varargout{2} = Edef;
varargout{3} = E_net;
varargout{4} = torques;
varargout{5} = power;
%++++++++++++++++++++++
% =========================================================================
% saving the results
if save_results
if noonlysolids
save(filename3,'cont_filename','cont_parameters','T_new','ET','EFi',...
'EFl','E_pot','Edef','E_net','torques','power');
else
save(filename3,'cont_filename','cont_parameters','T_new','ET','EFi',...
'EFl','E_pot');
end
fprintf('\n\n')
disp(['Solution saved in file ''',filename3,'.mat'''])
end
%##########################################################################
elseif ~exist([filename3,'.mat'],'file')
error(['Filename ''',filename3','.mat'' does not exist'])
else
load(filename3);
end
%##########################################################################
% figure
%##########################################################################
if draw_graph
% =========================================================================
% figure's properties
fig = figure(1);
set(fig,'visible','off','Position',[132 133 1156 676],'Name',...
'Energy balance, efficiency','numbertitle','off')
%==========================================================================
% plot 1
if noonlysolids
subplot(2,2,1)
end
%--------------------------------
% plotting
Y = [];
if noonlysolids
for kk = 1:nb_fishes;
Y = [ Y , E_net{kk},Edef{kk}]; %#ok<AGROW>
end
else
for kk = 1:nb_fishes;
Y = [ Y , EFi{kk}]; %#ok<AGROW>
end
end
Y = [Y , EFl,E_pot];
area(T_new,Y);
%--------------------------------
% color
colormap summer
%--------------------------------
%axis
set(gca,'Layer','top','xlim',[T_new(1),T_new(end)],'fontsize',12)
grid on
if latex
xlabel('Time','fontsize',16,'interpreter','latex')
ylabel('Energy (Joules)','fontsize',16,'interpreter','latex')
else
xlabel('Time','fontsize',16)
ylabel('Energy (Joules)','fontsize',16)
end
%----------------------------------
% title
if latex
title('Energy balance','fontsize',20,...
'interpreter','latex')
else
title('Energy balance','fontsize',20)
end
% ---------------------------------
% legend
if noonlysolids
if buoyant
legend_fish = cell(1,2*nb_fishes+2);
else
legend_fish = cell(1,2*nb_fishes+1);
end
elseif buoyant
legend_fish = cell(1,nb_fishes+2);
else
legend_fish = cell(1,nb_fishes+1);
end
if latex
for kk = 1:nb_fishes
if noonlysolids
legend_fish{2*kk-1} = ['fish ',num2str(kk),': global motion~~'];
legend_fish{2*kk} = ['fish ',num2str(kk),': deformations'];
else
legend_fish{kk} = ['Solid ',num2str(kk)];
end
end
if noonlysolids
legend_fish{2*nb_fishes+1} = 'fluid';
if buoyant
legend_fish{2*nb_fishes+2} = 'buoyant';
end
else
legend_fish{nb_fishes+1} = 'fluid';
if buoyant
legend_fish{nb_fishes+2} = 'buoyant';
end
end
h = legend(legend_fish);
set(h,'interpreter','latex')
else
for kk = 1:nb_fishes
if noonlysolids
legend_fish{2*kk-1} = ['fish ',num2str(kk),': global motion'];
legend_fish{2*kk} = ['fish ',num2str(kk),': deformations'];
else
legend_fish{kk} = ['Solid ',num2str(kk)];
end
end
if noonlysolids
legend_fish{2*nb_fishes+1} = 'fluid';
if buoyant
legend_fish{2*nb_fishes+2} = 'buoyant';
end
else
legend_fish{nb_fishes+1} = 'fluid';
if buoyant
legend_fish{nb_fishes+2} = 'buoyant';
end
end
legend(legend_fish);
end
%==========================================================================
if noonlysolids
%==========================================================================
% plot 2
subplot(2,2,2)
%--------------------------------
% plotting
graph_torques = cell(1,2*(nb_solids-nb_fishes));
for k = 1:nb_solids-nb_fishes
graph_torques{2*k+1} = T_new;
graph_torques{2*k+2} = torques(:,k);
end
plot(graph_torques{:})
% -------------------------------
% axis
set(gca,'fontsize',12,'Layer','top','xlim',[T_new(1),T_new(end)])
if latex
xlabel('Time','fontsize',16,'interpreter','latex')
ylabel('Power (Watts)','fontsize',16,'interpreter','latex')
else
xlabel('Time','fontsize',16)
ylabel('Torques (N.m)','fontsize',16)
end
grid on
%---------------------------------
% title
if latex
title('Torques of the articulations','fontsize',20,...
'interpreter','latex')
else
title('Torques of the articulations','fontsize',20)
end
%-----------------------------------
% legend
legend0 = cell(1,nb_solids-nb_fishes);
for k = 1:nb_solids-nb_fishes
if latex
legend0{k} = ['$c_',num2str(k),'$'];
else
legend0{k} = ['c(',num2str(k),')'];
end
end
h = legend(legend0);
set(h,'fontsize',14)
if latex
set(h,'interpreter','latex')
end
%==========================================================================
% plot 3
subplot(2,2,3)
%--------------------------------
% plotting
Y = [];
for kk = 1:nb_fishes;
Y = [Y,100*E_net{kk}./ET,100*Edef{kk}./ET]; %#ok<AGROW>
end
Y = [Y,100*EFl./ET];
if buoyant
Y = [Y,100*E_pot./ET];
end
area(T_new,Y);
% -------------------------------
% color
colormap summer
% -------------------------------
% axis
set(gca,'Layer','top','fontsize',12,'xlim',[T_new(1),T_new(end)])
grid on
if latex
xlabel('Time','fontsize',16,'interpreter','latex')
ylabel('Energy (\%)','interpreter','latex','fontsize',16)
else
xlabel('Time','fontsize',16)
ylabel('Energy (%)','fontsize',16)
end
%---------------------------------
% title
if latex
title ('Energy balance (percentages)','fontsize',20,...
'interpreter','latex')
else
title ('Energy balance (percentages)','fontsize',20)
end
%----------------------------------
% legend
if latex
for kk = 1:nb_fishes
legend_fish{2*kk-1} = ['fish ',num2str(kk),': global motion~~'];
legend_fish{2*kk} = ['fish ',num2str(kk),': deformations'];
end
legend_fish{2*nb_fishes+1} = 'fluid';
if buoyant
legend_fish{2*nb_fishes+2} = 'buoyant';
end
h = legend(legend_fish);
set(h,'interpreter','latex')
else
for kk = 1:nb_fishes
legend_fish{2*kk-1} = ['fish ',num2str(kk),': global motion'];
legend_fish{2*kk} = ['fish ',num2str(kk),': deformations'];
end
legend_fish{2*nb_fishes+1} = 'fluid';
if buoyant
legend_fish{2*nb_fishes+2} = 'buoyant';
end
legend(legend_fish);
end
%==========================================================================
% plot 4
subplot(2,2,4)
%--------------------------------
% plotting
% plotting
graph_power = cell(1,2*(nb_solids-nb_fishes));
for k = 1:nb_solids-nb_fishes
graph_power{2*k+1} = T_new;
graph_power{2*k+2} = power(:,k);
end
plot(graph_power{:})
% -------------------------------
% axis
set(gca,'fontsize',12,'Layer','top','xlim',[T_new(1),T_new(end)])
if latex
xlabel('Time','fontsize',16,'interpreter','latex')
ylabel('Power (Watts)','fontsize',16,'interpreter','latex')
else
xlabel('Time','fontsize',16)
ylabel('Power (Watts)','fontsize',16)
end
grid on
%---------------------------------
% title
if latex
title('Power expended by the articulations','fontsize',20,...
'interpreter','latex')
else
title('Power expended by the articulations','fontsize',20)
end
%-----------------------------------
% legend
h = legend(legend0);
set(h,'fontsize',14)
if latex
set(h,'interpreter','latex')
end
end
%==========================================================================
set(fig,'visible','on')
end
%##########################################################################
% Appendix II
% function mnistrom_new
%##########################################################################
function [uu,duu] = mnystrom(X,N)
% this function computes the potentials and its tangential
% derivative on the boundaries of the fluid
for j = 1:nb_solids
phi{j,1} = N{j}(1,:);
phi{j,2} = N{j}(2,:);
phi{j,3} = bound_data(j).flux.*bound_data(j).dl;
for k1 = 1:3
toto{j,k1} = real(ifft(fft(phi{j,k1})./[1,1:n(j),n(j):-1:1]));
end
end
for ki = 1:nb_boundaries
for kj = 1:nb_boundaries
% computing the matricial bloc A(ki,kj)
dx1 = repmat(X{kj}(1,:),m(ki),1) - repmat((X{ki}(1,:))',1,m(kj));
dx2 = repmat(X{kj}(2,:),m(ki),1) - repmat((X{ki}(2,:))',1,m(kj));
dist2 = dx1.^2 + dx2.^2;
A(mc(ki)+1:mc(ki+1),mc(kj)+1:mc(kj+1)) = dt(kj)*(dx1.*repmat(N{kj}(1,:),m(ki),1)...
+ dx2.*repmat(N{kj}(2,:),m(ki),1))./dist2;
if ki == kj
% the diagonal bloc A(ki,kj) must be recomputed as
% follows :
id_diag = sub2ind([mtot,mtot],mc(ki)+1:mc(ki+1),mc(ki)+1:mc(ki+1));
A(id_diag) = -pi - 0.5*dt(ki)*bound_data(ki).K.*bound_data(ki).dl;
end
% second right hand side
if kj <= nb_solids
if kj ~= ki
J0 = 3*(kj-1);
for j=1:3
b(mc(ki)+1:mc(ki+1),J0+j) = 0.5*dt(kj)*log(dist2)*phi{kj,j}';
end
else % ki == kj the kernel has to be splited into two parts
temp = abs(sqrt(dist2)./(2*sin(0.5*(repmat(bound_data(ki).t,m(ki),1)...
-repmat(bound_data(ki).t',1,m(ki))))));
id_diag = sub2ind([m(ki),m(ki)],1:m(ki),1:m(ki));
temp(id_diag) = bound_data(ki).dl;
temp = log( exp(0.5)*temp );
J0 = 3*(ki-1);
for j=1:3
b(mc(ki)+1:mc(ki+1),J0+j) = dt(ki)* (temp * phi{ki,j}') - pi*toto{ki,j}';
end
end
end
end
end
% computation of the potentials
if gene_data.pb_type == 1;
nA = size(A,1);
A(1:end-1,:) = A(1:end-1,:) - ones(nA-1,1)*A(end,:);
b(1:end-1,:) = b(1:end-1,:) - ones(nA-1,1)*b(end,:);
A(end,:) = discret_data.last_line;
b(end,:) = 0;
end;
u = A\b;
%------------------------------
uu = cell(nb_solids,1);
duu = cell(nb_solids,1);
%==============================
% computation of the tangiential derivative
for j=1:nb_solids
uu{j}=u(mc(j)+1:mc(j+1),:);
duu{j}=ifft(i*([0:n(j),-n(j):-1]'*...
ones(1,3*nb_solids)).*fft(uu{j}))./(bound_data(j).dl'*ones(1,3*nb_solids));
end
end
%##########################################################################
end
