
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