https://github.com/Data2Dynamics/d2d
Tip revision: 4c697ee160f4d91719a65cfcfe3e3b50395b6656 authored by Helge Hass on 29 July 2015, 07:45:33 UTC
new C code version
new C code version
Tip revision: 4c697ee
arCalcUnits.m
function arCalcUnits(m)
clc
global ar
if(~exist('m','var'))
m = 1;
end
%% reset all entries
for jx = 1:length(ar.model(m).x)
ar.model(m).xUnits{jx,1} = '';
end
%% species units
for jv = 1:length(ar.model(m).fv)
% [ar.model(m).x' ar.model(m).xUnits]
xs = union(ar.model(m).fv_source{jv}, ar.model(m).fv_target{jv});
q = ismember(ar.model(m).x, xs);
xu = ar.model(m).xUnits(q,1);
xul = ar.model(m).xUnits(q,2);
xuu = unique(xu);
xuul = unique(xul);
if(length(xuul)>1)
fprintf(2, 'Inconsistent species units:\n');
for jj=find(q)
fprintf(2, '\t%s [%s]\n', ar.model(m).x{jj}, ar.model(m).xUnits{jj,2});
end
error('Inconsistent species units');
end
xuu = setdiff(xuu, {''});
if(isempty(xuu))
ar.model(m).xUnits(q,1) = repmat(xs(1), sum(q),1);
elseif(length(xuu)==1)
ar.model(m).xUnits(q,1) = repmat(xuu(1), sum(q),1);
else
ar.model(m).xUnits(q,1) = repmat(xuu(1), sum(q),1);
for jxuu=2:length(xuu)
ar.model(m).xUnits(:,1) = strrep(ar.model(m).xUnits(:,1), ...
xuu(jxuu), xuu(1));
end
end
end
known_vars = ar.model(m).x;
known_units = ar.model(m).xUnits(:,1)';
%% parameter units for volumes and initial conditions
ar.model(m).pUnits = cell(size(ar.model(m).p));
for jp = 1:length(ar.model(m).p)
% volumes
q = ismember(ar.model(m).pc, ar.model(m).p{jp});
if(sum(q) == 1)
ar.model(m).pUnits(jp) = ar.model(m).c(q);
known_vars{end+1} = ar.model(m).p{jp}; %#ok<AGROW>
known_units(end+1) = ar.model(m).pUnits(jp); %#ok<AGROW>
continue
elseif(sum(q) > 1)
error('double matches');
end
% initial conditions
q = ismember(ar.model(m).px0, ar.model(m).p{jp});
if(sum(q) == 1)
ar.model(m).pUnits(jp) = ar.model(m).xUnits(q,1);
known_vars{end+1} = ar.model(m).p{jp}; %#ok<AGROW>
known_units(end+1) = ar.model(m).pUnits(jp); %#ok<AGROW>
continue
elseif(sum(q) > 1)
error('double matches');
end
end
%% remaining dynamic parameters
px = {};
for jp = 1:length(ar.model(m).p)
q = ismember(ar.model(m).pv, ar.model(m).p{jp});
if(sum(q) == 1)
px{end+1} = ar.model(m).p{jp}; %#ok<AGROW>
elseif(sum(q) > 1)
error('double matches');
end
end
%% setup flux equations for unit solver
fv = ar.model(m).fv;
fv = sym(fv);
fv = subs(fv, sym(known_vars), sym(known_units));
for jv = 1:length(fv)
xs = union(ar.model(m).fv_source{jv}, ar.model(m).fv_target{jv});
fv(jv) = fv(jv) * sym(ar.model(m).tUnits{1}) / sym(xs{1}) - 1;
end
punits = unit_solver(fv, px);
% %% remaining dynamic parameters
%
% for jp = 1:length(ar.model(m).p)
% q = ismember(ar.model(m).pv, ar.model(m).p{jp});
% if(sum(q) == 1)
% ar.model(m).p{jp}
% q = ismembern(ar.model(m).pvs,ar.model(m).p{jp});
% vs = ar.model(m).fv(q);
% vsi = find(q);
% xs = union(ar.model(m).fv_source{vsi(1)}, ar.model(m).fv_target{vsi(1)});
% xsu = ar.model(m).xUnits{ismember(ar.model(m).x, xs{1}),1};
% vs = sym(vs);
% vs = subs(vs, sym(ar.model(m).x), sym(ar.model(m).xUnits(:,1)'))
% vs = vs * sym(ar.model(m).tUnits{1}) / sym(xsu) - 1;
% solve(vs(:), ar.model(m).p{jp})
%
% ptmp = '';
%
% ar.model(m).pUnits{jp} = ptmp;
% continue
% elseif(sum(q) > 1)
% error('double matches');
% end
% end
function q = ismembern(c,l)
q = logical(size(c));
for j=1:length(c)
q(j) = sum(ismember(c{j},l))>0;
end
function punits = unit_solver(fv, px)
psolved = {};
punits = {};
px = sym(px);
for j=1:length(px)
fvtmp = findasym(fv, px(j));
px(j)
fvtmp
S = solve(fvtmp, px(j))
end
function fxout = findasym(fx, x)
fxout = sym([]);
for j=1:length(fx)
q = sum(ismember(symvar(fx(j)),x))>0;
if(q)
fxout(end+1) = fx(j); %#ok<AGROW>
end
end