function diag2()
% Uq or Us diagonalization
% Run after diag1.m and MatrixU.m 

global  InfTable Bx By Bz

% define foldername
modelfolder = getmodel;

% read-write model files
IT_filename = strcat(modelfolder,'InfTable.txt');
fid = fopen(IT_filename,'r');
InfTable = fscanf(fid,'%i',9);
fclose(fid);
nAtoms = InfTable(1);
nVibrCoords = InfTable(7);
nBlocks = InfTable(8);
mVibrs = InfTable(9);
ns = zeros(1,nBlocks);
ms = zeros(1,nBlocks);
filename = strcat(modelfolder,'ms.txt');
fid = fopen(filename,'rt');
ms = fscanf(fid,'%i',nBlocks);
fclose(fid);

if nBlocks == 1
    Ufilename = strcat(modelfolder,'Uq.txt');
    ns(1) = nVibrCoords;
    Lqfilename = strcat(modelfolder,'Lq.txt');
    Lpfilename = strcat(modelfolder,'Lp.txt');
else
    Ufilename = strcat(modelfolder,'Us.txt');
    C = zeros(nVibrCoords,nVibrCoords);
    filename = strcat(modelfolder,'Cmatrix.txt');
    Cfid = fopen(filename,'rt');
    s = fgetl(Cfid);
    ns = fscanf(Cfid,'%i',nBlocks);
    irow = 0;
    for iblock = 1:nBlocks
        nrows = ns(iblock);
        for i = 1:nrows
            irow = irow + 1;
            ncols = fscanf(Cfid,'%i',1);
            for jcol = 1:ncols
                j = fscanf(Cfid,'%i',1);
                c = fscanf(Cfid,'%f',1);
                C(irow, j) = c;
            end % jcol
        end % i
    end % iblock
    fclose(Cfid);
    Lqfilename = strcat(modelfolder,'Ls.txt');
    Lpfilename = strcat(modelfolder,'Ps.txt');
end % nBlocks == 1

Ufid = fopen(Ufilename,'rt');
filename = strcat(modelfolder,'Bdr.txt');
Bdrfid = fopen(filename,'rt');
filename = strcat(modelfolder,'Tau.txt');
Taufid = fopen(filename,'rt');
Lqfid = fopen(Lqfilename,'w');
Lpfid = fopen(Lpfilename,'w');
freqs = [];

rest = nVibrCoords - 4*floor(nVibrCoords/4);
for iblock = 1:nBlocks
    nrows = ns(iblock); % rows in Us 
    nUs3 = nrows*(nrows+1)/2; % elements in Us upper triangle
    Us3 = zeros(1,nUs3); % Us upper triangle
    Us3 = fscanf(Ufid,'%f',nUs3);
    U  = zeros(nrows); % Us3 -> U square
    k = 0;
    for i = 1:nrows
        for j = i:nrows
            k = k + 1;
            t = Us3(k);
            U(i,j) = t;
            U(j,i) = t;
        end % j
    end % i
    m = ms(iblock);
    Ltau = zeros(nrows,m);
    tau = zeros(m,1);
    tau = fscanf(Taufid,'%e',m);
    for i = 1:m
        Ltau(:,i) = fscanf(Bdrfid,'%e',nrows);
    end % i  
    Utau = Ltau'*U*Ltau;
    [Lu, d] = eig(Utau);
    d = diag(d);
    [dsort,Ind] = sort(d);
    Lusort = Lu(:,Ind);
    niu = sqrt(dsort)*1000; % niu is alias for freqs within this block
    freqs = [freqs,niu'];
    Lq = Ltau*Lusort;
    for i = 1:m
        fprintf(Lqfid,'   %13.6f   %13.6f   %13.6f   %13.6f\n',Lq(:,i));
        if (rest~=0),fprintf(Lqfid,'\n');end
    end % i
    Ltau2 = zeros(nrows,m);
    for i = 1:m
        Ltau2(:,i) = Ltau(:,i)/(tau(i)^2);
    end % i    
    Lp = Ltau2*Lusort;
    for i = 1:m
        fprintf(Lpfid,'   %13.6f   %13.6f   %13.6f   %13.6f\n',Lp(:,i));
        if (rest~=0),fprintf(Lpfid,'\n');end
    end % i
end % iblock
fclose(Ufid);
fclose(Bdrfid);
fclose(Taufid);
fclose(Lqfid);
fclose(Lpfid);

if nBlocks > 1
    Lsfilename = strcat(modelfolder,'Ls.txt');
    Lsfid = fopen(Lsfilename,'rt');
    Psfilename = strcat(modelfolder,'Ps.txt');
    Psfid = fopen(Psfilename,'rt');
    Lqfilename = strcat(modelfolder,'Lq.txt');
    Lqfid = fopen(Lqfilename,'w');
    Lpfilename = strcat(modelfolder,'Lp.txt');
    Lpfid = fopen(Lpfilename,'w');
    ic1 = 1;
    for iblock = 1:nBlocks
        nrows = ns(iblock); % rows in Ls or Ps
        ic2 = ic1 + nrows - 1;
        m = ms(iblock);
        Lqblock = zeros(nVibrCoords,m);
        Lpblock = zeros(nVibrCoords,m);
        Lsblock = zeros(nrows,m);
        Psblock = zeros(nrows,m);
        for i = 1:m
            Lsblock(:,i) = fscanf(Lsfid,'%f',nrows);
            Psblock(:,i) = fscanf(Psfid,'%f',nrows);
        end % i
        Lqblock = C(ic1:ic2,:)'*Lsblock;
        Lpblock = C(ic1:ic2,:)'*Psblock;
        for i = 1:m
            fprintf(Lqfid,'   %13.6f   %13.6f   %13.6f   %13.6f\n',Lqblock(:,i));
            if (rest~=0),fprintf(Lqfid,'\n');end
            fprintf(Lpfid,'   %13.6f   %13.6f   %13.6f   %13.6f\n',Lpblock(:,i));
            if (rest~=0),fprintf(Lpfid,'\n');end
        end % i
        ic1 = ic1+ nrows;
    end % iblock
    fclose(Lqfid);
    fclose(Lpfid);
    fclose(Lsfid);
    fclose(Psfid);
    clear C
end % if nBlocks > 1

Lpfilename = strcat(modelfolder,'Lp.txt');
Lpfid = fopen(Lpfilename,'rt');
Shiftsfilename = strcat(modelfolder,'Shifts.txt');
Shiftsfid = fopen(Shiftsfilename,'w');

filename = strcat(modelfolder,'AtomMasses.txt');
load(filename);
E = eye(nAtoms);
for i = 1:nAtoms
    E(i,i) = 1.088/AtomMasses(i);
end
clear AtomMasses

matrb(modelfolder);

rx = zeros(nAtoms,nVibrCoords);
ry = zeros(nAtoms,nVibrCoords);
rz = zeros(nAtoms,nVibrCoords);
Lpcolumn = zeros(nVibrCoords,1);

for i = 1:mVibrs
    Lpcolumn = fscanf(Lpfid,'%f',nVibrCoords);
    rx = E*Bx'*Lpcolumn;
    ry = E*By'*Lpcolumn;
    rz = E*Bz'*Lpcolumn;
    for j = 1:nAtoms
     fprintf(Shiftsfid,'   %13.6f   %13.6f   %13.6f \n',rx(j),ry(j),rz(j));
    end % j
end % i

fclose(Shiftsfid);
fclose(Lpfid);

filename = strcat(modelfolder,'freqs.txt');
fid = fopen(filename,'w');
fprintf(fid, ...
'     %11.6f     %11.6f     %11.6f     %11.6f \n',freqs);
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
