function MatrixT()
% forms Tp or Ts

global  InfTable Bx By Bz

% define foldername
modelfolder = getmodel;

% read model files
InfTable = zeros(1,9);
filename = strcat(modelfolder,'InfTable.txt');
fid = fopen(filename,'r');
InfTable = fscanf(fid,'%i',8);
fclose(fid);
nAtoms = InfTable(1);
nBonds = InfTable(2);
nAlphas = InfTable(3);
nGammas = InfTable(4);
nRos = InfTable(5);
nHis = InfTable(6);
nVibrCoords = InfTable(7);
nBlocks = InfTable(8);

filename = strcat(modelfolder,'AtomMasses.txt');
load(filename);

matrb(modelfolder);

Tp = zeros(nVibrCoords,nVibrCoords);
E = eye(nAtoms);
for i = 1:nAtoms
    E(i,i) = 1.088/AtomMasses(i);
end
Tp = Bx*E*Bx' + By*E*By' + Bz*E*Bz';
% clear Bx By Bz E AtomMasses

if nBlocks == 1
    filename = strcat(modelfolder,'Tp.txt');
    fid = fopen(filename,'w');
    nTp3 = nVibrCoords*(nVibrCoords+1)/2; % elements in Tp upper triangle
    Tp3 = zeros(1,nTp3); % Tp upper triangle
    k = 0;
    for i = 1:nVibrCoords
        for j = i:nVibrCoords
            k = k + 1;
            Tp3(k) = Tp(i,j);
        end % j    
    end % i
    fprintf(fid,'  %e  %e  %e  %e \n',Tp3);
    fclose(fid);
    return
end % nBlocks == 1

% Then form Ts
C = zeros(nVibrCoords,nVibrCoords);
nC = zeros(1,nBlocks);
filename = strcat(modelfolder,'Cmatrix.txt');
fid = fopen(filename,'rt');
s = fgetl(fid);
nC = fscanf(fid,'%i',nBlocks);
irow = 0;
for iblock = 1:nBlocks
    nrows = nC(iblock);
    for i = 1:nrows
        irow = irow + 1;
        ncols = fscanf(fid,'%i',1);
        for jcol = 1:ncols
            j = fscanf(fid,'%i',1);
            c = fscanf(fid,'%f',1);
            C(irow, j) = c;
        end % jcol
    end % i
end % iblock
fclose(fid);
Ts = C*Tp*C';

filename = strcat(modelfolder,'Ts.txt');
fid = fopen(filename,'w');
i1 = 0;
i2 = 0;
for iblock = 1:nBlocks
    nTs3 = nC(iblock)*(nC(iblock)+1)/2; % elements in Ts upper triangle
    Ts3 = zeros(1,nTs3); % Tp upper triangle
    k = 0;
    i1 = i2 + 1;
    i2 = i2 + nC(iblock);
    for i = i1:i2
        for j= i:nVibrCoords
            if j <= i2
                k = k + 1;
                Ts3(k) = Ts(i,j);
            else
                tsij = Ts(i,j);
                if abs(tsij) > 0.0005
                  s = sprintf('Error: i=%i j=%i Ts(i,j)=%f',i,j,tsij);
                    disp(s)
                    return
               end % if
            end % if
        end % j
    end % i
    fprintf(fid,'  %e  %e  %e  %e \n',Ts3);
    fprintf(fid,' \n');
end % iblock
fclose(fid);

% subfunctions ---------------------
function modelfolder = getmodel()
% gets model name in LevInfinite

filename = 'C:\LevInfinite\Config\Molecules.ini';
fid = fopen(filename,'rt');
molfolder = fgetl(fid);
fclose(fid);

filename = 'C:\LevInfinite\Config\ToCalculate.ini';
fid = fopen(filename,'rt');
modelname = fgetl(fid);
fclose(fid);

modelfolder = strcat(molfolder,modelname);
modelfolder = strcat(modelfolder,'\');

function matrb(modelfolder)
% Before run create txt files within THIS folder.

global  InfTable Bx By Bz

nAtoms = InfTable(1);
nBonds = InfTable(2);
nAlphas = InfTable(3);
nGammas = InfTable(4);
nRos = InfTable(5);
nHis = InfTable(6);
n = InfTable(7);
Bx = zeros(n,nAtoms);
By = zeros(n,nAtoms);
Bz = zeros(n,nAtoms);

ib = 0;
fullname = strcat(modelfolder,'BqStructures.txt');
load(fullname);
for i = 1:nBonds
    ib = ib + 1;
    ia = BqStructures(i,1);
    Bx(ib,ia) = BqStructures(i,3);
    By(ib,ia) = BqStructures(i,4);
    Bz(ib,ia) = BqStructures(i,5);
    ia = BqStructures(i,2);
    Bx(ib,ia) = BqStructures(i,6);
    By(ib,ia) = BqStructures(i,7);
    Bz(ib,ia) = BqStructures(i,8);
end % i
clear BqStructures
if nAlphas > 0
    fullname = strcat(modelfolder,'BaStructures.txt');
    load(fullname);
	for i = 1:nAlphas
        ib = ib + 1;
        ia = BaStructures(i,1);
        Bx(ib,ia) = BaStructures(i,4);
        By(ib,ia) = BaStructures(i,5);
        Bz(ib,ia) = BaStructures(i,6);
        ia = BaStructures(i,2);
        Bx(ib,ia) = BaStructures(i,7);
        By(ib,ia) = BaStructures(i,8);
        Bz(ib,ia) = BaStructures(i,9);
        ia = BaStructures(i,3);
        Bx(ib,ia) = BaStructures(i,10);
        By(ib,ia) = BaStructures(i,11);
        Bz(ib,ia) = BaStructures(i,12);
	end % i = 1:nAlphas
    clear BaStructures
end % if nAlphas > 0
if nGammas > 0
    fullname = strcat(modelfolder,'BgStructures.txt');
    load(fullname);
    for i = 1:nGammas
        ib = ib + 1;
        ia = BgStructures(i,1);
        Bx(ib,ia) = BgStructures(i,4);
        By(ib,ia) = BgStructures(i,5);
        Bz(ib,ia) = BgStructures(i,6);
        ia = BgStructures(i,2);
        Bx(ib,ia) = BgStructures(i,7);
        By(ib,ia) = BgStructures(i,8);
        Bz(ib,ia) = BgStructures(i,9);
        ia = BgStructures(i,3);
        Bx(ib,ia) = BgStructures(i,10);
        By(ib,ia) = BgStructures(i,11);
        Bz(ib,ia) = BgStructures(i,12);
	end % i = 1:nGammas
    clear BgStructures
end % if nGammas > 0
if nRos > 0
    fullname = strcat(modelfolder,'BrStructures.txt');
    load(fullname);
    for i = 1:nRos
        ib = ib + 1;
        ia = BrStructures(i,1);
        Bx(ib,ia) = BrStructures(i,5);
        By(ib,ia) = BrStructures(i,6);
        Bz(ib,ia) = BrStructures(i,7);
        ia = BrStructures(i,2);
        Bx(ib,ia) = BrStructures(i,8);
        By(ib,ia) = BrStructures(i,9);
        Bz(ib,ia) = BrStructures(i,10);
        ia = BrStructures(i,3);
        Bx(ib,ia) = BrStructures(i,11);
        By(ib,ia) = BrStructures(i,12);
        Bz(ib,ia) = BrStructures(i,13);
        ia = BrStructures(i,4);
        Bx(ib,ia) = BrStructures(i,14);
        By(ib,ia) = BrStructures(i,15);
        Bz(ib,ia) = BrStructures(i,16);
    end % i = 1:nRos
    clear BrStructures
end % if nRos > 0
if nHis > 0
    fullname = strcat(modelfolder,'BhStructures.txt');
    load(fullname);
	for i = 1:nHis
        ib = ib + 1;
        ia = BhStructures(i,1);
        Bx(ib,ia) = Bx(ib,ia)+BhStructures(i,7);
        By(ib,ia) = By(ib,ia)+BhStructures(i,8);
        Bz(ib,ia) = Bz(ib,ia)+BhStructures(i,9);
        ia = BhStructures(i,2);
        Bx(ib,ia) = Bx(ib,ia)+BhStructures(i,10);
        By(ib,ia) = By(ib,ia)+BhStructures(i,11);
        Bz(ib,ia) = Bz(ib,ia)+BhStructures(i,12);
        ia = BhStructures(i,3);
        Bx(ib,ia) = Bx(ib,ia)+BhStructures(i,13);
        By(ib,ia) = By(ib,ia)+BhStructures(i,14);
        Bz(ib,ia) = Bz(ib,ia)+BhStructures(i,15);
        ia = BhStructures(i,4);
        Bx(ib,ia) = Bx(ib,ia)+BhStructures(i,16);
        By(ib,ia) = By(ib,ia)+BhStructures(i,17);
        Bz(ib,ia) = Bz(ib,ia)+BhStructures(i,18);
        ia = BhStructures(i,5);
        Bx(ib,ia) = Bx(ib,ia)+BhStructures(i,19);
        By(ib,ia) = By(ib,ia)+BhStructures(i,20);
        Bz(ib,ia) = Bz(ib,ia)+BhStructures(i,21);
        ia = BhStructures(i,6);
        Bx(ib,ia) = Bx(ib,ia)+BhStructures(i,22);
        By(ib,ia) = By(ib,ia)+BhStructures(i,23);
        Bz(ib,ia) = Bz(ib,ia)+BhStructures(i,24);
    end % i = 1:nHis
    clear BhStructures
end % if nHis > 0
