- Author:
- WeiweiAi <wai484@aucklanduni.ac.nz>
- Date:
- 2023-04-04 15:51:18+12:00
- Desc:
- Add a steady state example
- Permanent Source URI:
- https://models.fieldml.org/workspace/6bc/rawfile/a7940fcacd455a09bdb7beac0da91eca79cf1adb/BG_fit/buildBG_matrix.m
function buildBG_matrix(N_f,N_r,Kname,kappa_name,importComp,newComp,modelName,modelPath)
%% Write to CellML
blank='';
% Model
modelfile=[modelPath modelName '.txt'];
fileID = fopen(modelfile,'w');
fprintf(fileID,'def model %s as\n',modelName);
% Unit
unitimport =["J_per_mol","fmol","fmol_per_sec","per_fmol"];
blanks=strcat('',strings(1,length(unitimport)));
unit_all=[blanks;unitimport;unitimport];
fprintf(fileID,'%3s def import using "../cellLib/BG/units_BG.cellml" for\n',blank);
fprintf(fileID,'%7s unit %s using unit %s;\n',unit_all(:));
fprintf(fileID,'%3s enddef;\n',blank);
fprintf(fileID,'\n');
% Import predefined models
for i=1:length(importComp)
if ~contains(importComp(i),"port")
fprintf(fileID,'%3s def import using "../cellLib/BG/%s_BG.cellml" for\n',blank,importComp(i));
fprintf(fileID,'%7s comp %s using comp %s_BG;\n',blank,newComp(i),importComp(i));
fprintf(fileID,'%3s enddef;\n',blank);
fprintf(fileID,'\n');
end
end
% Define the component
fprintf(fileID,'%3s def comp %s as\n',blank,modelName);
fprintf(fileID,'%7s var T: kelvin {pub: in, priv: out};\n',blank);
fprintf(fileID,'%7s var t: second {pub: in, priv: out};\n',blank);
% paras
for i=1:length(newComp)
ctype=importComp(i);
switch ctype
case 'Ce'
fprintf(fileID,'%7s var q_%s_init: fmol {pub: in, priv: out};\n',blank,newComp(i));
fprintf(fileID,'%7s var K_%s: per_fmol {pub: in, priv: out};\n',blank,newComp(i));
fprintf(fileID,'%7s var q_%s: fmol {pub: out, priv: in};\n',blank,newComp(i));
fprintf(fileID,'%7s var mu_%s: J_per_mol {pub: out, priv: in};\n',blank,newComp(i));
fprintf(fileID,'%7s var v_%s: fmol_per_sec {pub: out,priv: out};\n',blank,newComp(i));
case 'Se'
fprintf(fileID,'%7s var q_%s_init: fmol {pub: in, priv: out};\n',blank,newComp(i));
fprintf(fileID,'%7s var K_%s: per_fmol {pub: in, priv: out};\n',blank,newComp(i));
% fprintf(fileID,'%7s var q_%s: fmol {pub: out, priv: in};\n',blank,newComp(i));
fprintf(fileID,'%7s var mu_%s: J_per_mol {pub: out, priv: in};\n',blank,newComp(i));
fprintf(fileID,'%7s var v_%s: fmol_per_sec {pub: out,priv: out};\n',blank,newComp(i));
case 'port_Ce'
fprintf(fileID,'%7s var mu_%s: J_per_mol {pub: in};\n',blank,newComp(i));
fprintf(fileID,'%7s var v_%s: fmol_per_sec {pub: out};\n',blank,newComp(i));
case 'Re'
fprintf(fileID,'%7s var kappa_%s: fmol_per_sec {pub: in, priv: out};\n',blank,newComp(i));
fprintf(fileID,'%7s var v_%s: fmol_per_sec {pub: out,priv: in};\n',blank,newComp(i));
fprintf(fileID,'%7s var mu_%s_f: J_per_mol {pub: out,priv: out};\n',blank,newComp(i));
fprintf(fileID,'%7s var mu_%s_r: J_per_mol {pub: out,priv: out};\n',blank,newComp(i));
case 'port_Re'
fprintf(fileID,'%7s var mu_%s: J_per_mol {pub: out};\n',blank,newComp(i));
fprintf(fileID,'%7s var v_%s: fmol_per_sec {pub: in};\n',blank,newComp(i));
otherwise
fprintf('The component %s is not defined\n',ctype);
return
end
end
fprintf(fileID,'\n');
% Equations for conservation laws
for i=1:length(Kname)
findex=find(N_f(i,:));
if ~isempty(findex)
fRe=kappa_name(findex);
formatSpec = "-%f{dimensionless}*v_%s";
sfRe = strjoin(compose(formatSpec,N_f(i,findex)',fRe)');
else
sfRe='';
end
rindex=find(N_r(i,:));
if ~isempty(rindex)
rRe=kappa_name(rindex);
formatSpec = "+%f{dimensionless}*v_%s";
srRe = strjoin(compose(formatSpec,N_r(i,rindex)',rRe)');
else
srRe='';
end
fprintf(fileID,'%7s v_%s = %s%s;\n',blank,Kname(i),sfRe,srRe);
end
for i=1:length(kappa_name)
findex=find(N_f(:,i));
if ~isempty(findex)
fSpecie=Kname(findex);
formatSpec = "+%f{dimensionless}*mu_%s";
sfSpecie = strjoin(compose(formatSpec,N_f(findex,i),fSpecie)');
fprintf(fileID,'%7s mu_%s_f = %s;\n',blank,kappa_name(i),sfSpecie);
end
rindex=find(N_r(:,i));
if ~isempty(rindex)
rSpecie=Kname(rindex);
formatSpec = "+%f{dimensionless}*mu_%s";
srSpecie = strjoin(compose(formatSpec,N_r(rindex,i),rSpecie)');
fprintf(fileID,'%7s mu_%s_r = %s;\n',blank,kappa_name(i),srSpecie);
end
end
fprintf(fileID,'%3s enddef;\n',blank);
fprintf(fileID,'\n');
% Group
fprintf(fileID,'%3s def group as encapsulation for\n',blank);
fprintf(fileID,'%7s comp %s incl\n',blank,modelName);
nnewComp=length(newComp);
blanks=strcat('',strings(1,nnewComp));
data_all=[blanks;newComp'];
fprintf(fileID,'%11s comp %s;\n',data_all);
fprintf(fileID,'%7s endcomp;\n',blank);
fprintf(fileID,'%3s enddef;\n',blank);
fprintf(fileID,'\n');
% Mapping
for i=1:length(newComp)
fprintf(fileID,'%3s def map between %s and %s for\n',blank,modelName,newComp(i));
ctype=importComp(i);
switch ctype
case 'Ce'
varsc=["q","K","v","mu"];
formatSpec = "%s_%s";
names=strcat(newComp(i),strings(length(varsc),1));
varsn=compose(formatSpec,varsc',names);
vars2=["T","t","q_init",varsc];
newq=sprintf("q_%s_init",newComp(i));
vars1=["T","t",newq,varsn'];
case 'Se'
varsc=["K","v","mu"];
formatSpec = "%s_%s";
names=strcat(newComp(i),strings(length(varsc),1));
varsn=compose(formatSpec,varsc',names);
vars2=["T","t","q_init",varsc];
newq=sprintf("q_%s_init",newComp(i));
vars1=["T","t",newq,varsn'];
case 'Re'
vars2=["T","kappa","A_f","A_r","v"];
nkappa=sprintf("kappa_%s",newComp(i));
nmuf=sprintf("mu_%s_f",newComp(i));
nmur=sprintf("mu_%s_r",newComp(i));
nv=sprintf("v_%s",newComp(i));
vars1=["T",nkappa,nmuf,nmur,nv];
case {'port_Ce','port_Re'}
% do nothing
otherwise
fprintf('The component %s is not defined\n',ctype);
return
end
nvar=length(vars1);
blanks=strcat('',strings(1,nvar));
data_all=[blanks;vars1;vars2];
fprintf(fileID,'%7s vars %s and %s;\n',data_all(:));
fprintf(fileID,'%3s enddef;\n',blank);
fprintf(fileID,'\n');
end
fprintf(fileID,'enddef;\n');
fprintf(fileID,'\n');
fclose(fileID);
end
