function intir()
% first, you should prepare
% freqs-shifts,
% miu and dmdq in a folder
% mentioned within ToCalculate.ini

% define foldername
modelfolder = getmodel;

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

Miu = zeros(1,nBonds);
filename = strcat(modelfolder,'Miu.txt');
fid = fopen(filename,'r');
Miu = fscanf(fid,'%f',nBonds);
fclose(fid);

filename = strcat(modelfolder,'dmdqSparse.txt');
load(filename);
dMdq = zeros(nBonds,nVibrCoords);
mdmdq = size(dmdqSparse, 1);
for i   = 1:mdmdq
    idmdq = dmdqSparse(i,1);
    jdmdq = dmdqSparse(i,2);
    r = dmdqSparse(i,3);
    dMdq(idmdq, jdmdq) = r;
end % i
dmdq = dMdq;
clear dMdq dmdqSparse

numba = zeros(nBonds,2);
filename = strcat(modelfolder,'BqStructures.txt');
load(filename);
for i = 1:nBonds
    numba(i,1) = BqStructures(i,1);
    numba(i,2) = BqStructures(i,2);
end % i
clear BqStructures

filename = strcat(modelfolder,'BondLengthes.txt');
load(filename);
sigma = 1.00./BondLengthes;
clear BondLengthes

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

IntegrIntIR = zeros(1,nFreqs);
AtomShifts = zeros(nAtoms,3);
filename = strcat(modelfolder,'Shifts.txt');
fidShifts = fopen(filename,'rt');
Lq = zeros(1,nVibrCoords);
filename = strcat(modelfolder,'Lq.txt');;
fidLq = fopen(filename,'rt');
filename = strcat(modelfolder,'IRSabsint.txt');
fidIRS = fopen(filename,'w');
for i = 1:nFreqs
    for iatom = 1:nAtoms
        AtomShifts(iatom,1:3) = fscanf(fidShifts,'%f',3);
    end % iatom
    Lq = fscanf(fidLq,'%f',nVibrCoords);
    Ix = 0;
    Iy = 0;
    Iz = 0;
    for iq = 1:nBonds
        Ix = Ix + Miu(iq)*sigma(iq)* ...
            ((AtomShifts(numba(iq,2),1)-AtomShifts(numba(iq,1),1))- ...
            BondVectors(iq,1)*Lq(iq));
        Iy = Iy + Miu(iq)*sigma(iq)* ...
            ((AtomShifts(numba(iq,2),2)-AtomShifts(numba(iq,1),2))- ...
            BondVectors(iq,2)*Lq(iq));
        Iz = Iz + Miu(iq)*sigma(iq)* ...
            ((AtomShifts(numba(iq,2),3)-AtomShifts(numba(iq,1),3))- ...
            BondVectors(iq,3)*Lq(iq));
    end % iq
    Ix = Ix + BondVectors(:,1)'*dmdq*Lq;
    Iy = Iy + BondVectors(:,2)'*dmdq*Lq;
    Iz = Iz + BondVectors(:,3)'*dmdq*Lq;
    IntegrIntIR(i) = 20*(Ix^2 + Iy^2 + Iz^2);
    fprintf(fidIRS,'%7.2f  %7.2f  %7.2f  %7.2f \n',IntegrIntIR(i), ...
            Ix, Iy, Iz);
end % i
fclose(fidShifts);
fclose(fidLq);
fclose(fidIRS);

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,'\');
