function diag1()
% Tp or Ts diagonalization
% Run after MatrixT.m 

% define foldername
modelfolder = getmodel;

% read model files
IT_filename = strcat(modelfolder,'InfTable.txt');
fid = fopen(IT_filename,'r');
InfTable = fscanf(fid,'%i',8);
fclose(fid);
nVibrCoords = InfTable(7);
nBlocks = InfTable(8);
ns = zeros(1,nBlocks);
ms = zeros(1,nBlocks);
if nBlocks == 1
    Tfilename = strcat(modelfolder,'Tp.txt');
    ns(1) = nVibrCoords;
else
    Tfilename = strcat(modelfolder,'Ts.txt');
    filename = strcat(modelfolder,'Cmatrix.txt');
    fid = fopen(filename,'rt');
    s = fgetl(fid);
    ns = fscanf(fid,'%i',nBlocks);
    fclose(fid);
end % nBlocks == 1
Tfid = fopen(Tfilename,'rt');
filename = strcat(modelfolder,'Bdr.txt');
Bdrfid = fopen(filename,'w');
Tau = [];

for iblock = 1:nBlocks
    nrows = ns(iblock); % rows in Ts 
    nTs3 = nrows*(nrows+1)/2; % elements in Ts upper triangle
    Ts3 = zeros(1,nTs3); % Ts upper triangle
    Ts3 = fscanf(Tfid,'%f',nTs3);
    Tp  = zeros(nrows); % Ts -> Tp
    k = 0;
    for i = 1:nrows
        for j = i:nrows
            k = k + 1;
            t = Ts3(k);
            Tp(i,j) = t;
            Tp(j,i) = t;
        end % j
    end % i
    [Ltau, d] = eig(Tp);
    d = diag(d);
    indzero = [];
    for i = 1:nrows
       if abs(d(i)) < 0.0001
           indzero = [indzero, i];
       end % if
    end % i
    d(indzero) = [];
    m = length(d);
    Ltau(:,indzero) = [];
    tau  = sqrt(d);
    for i = 1:m
        Ltau(:,i) = Ltau(:,i)*tau(i);
        fprintf(Bdrfid,'%e \n',Ltau(:,i));
    end % i    
    ms(iblock) = m;
    Tau = [Tau; tau];
end % iblock
fclose(Tfid);
fclose(Bdrfid);

filename = strcat(modelfolder,'Tau.txt');
Taufid = fopen(filename,'w');
fprintf(Taufid,'  %e  %e  %e  %e  \n',Tau);
fclose(Taufid);
filename = strcat(modelfolder,'ms.txt');
fid = fopen(filename,'w');
fprintf(fid,'  %i',ms);
fclose(fid);

InfTable(9) = sum(ms);
fid = fopen(IT_filename,'w');
fprintf(fid,' %11i %11i %11i %11i %11i %11i \n',InfTable);
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,'\');

