function JoinFragmNew()
% constructs vibrational model using fragment of ready models.
% Instructions for new model design are within
% LevInfinite\Config\ToJoin.txt.

global Type InfTable AtomTypes AtomCoords Bonds
global Alphas Gammas Ros His 

fragm_connect = zeros(2,4);
filename = 'C:\LevInfinite\Config\ToJoin.txt';
fid = fopen(filename,'rt');
s = fgetl(fid);
fragm1name = strtok(s);
s = fgetl(fid);
fragm_connect(1,1) = str2double(strtok(s));
s = fgetl(fid);
fragm_connect(1,2) = str2double(strtok(s));
s = fgetl(fid);
fragm_connect(1,3) = str2double(strtok(s));
s = fgetl(fid);
fragm_connect(1,4) = str2double(strtok(s));
s = fgetl(fid);
fragm2name = strtok(s);
s = fgetl(fid);
fragm_connect(2,1) = str2double(strtok(s));
s = fgetl(fid);
fragm_connect(2,2) = str2double(strtok(s));
s = fgetl(fid);
fragm_connect(2,3) = str2double(strtok(s));
s = fgetl(fid);
fragm_connect(2,4) = str2double(strtok(s));
s = fgetl(fid);
jointname = strtok(s);
fclose(fid);

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

fullname = strcat(modelfolder,'Bstruct.mat');
load(fullname);

% preserve e2 for turn 1 of fragment 2
% preserve e4 for turn 2 of fragment 2
e2 = atom_pair_vector(fragm_connect(2,2), fragm_connect(2,1));
e4 = atom_pair_vector(fragm_connect(2,3), fragm_connect(2,4));

% Shelter Fragm2-globals
InfTable_Fragm2 = InfTable;
Type_Fragm2 = Type;
AtomTypes_Fragm2 = AtomTypes;
AtomCoords_Fragm2 = AtomCoords;
Bonds_Fragm2 = Bonds;
Alphas_Fragm2 = Alphas;
Gammas_Fragm2 = Gammas;
Ros_Fragm2 = Ros;
His_Fragm2 = His;

modelfolder = strcat(molfolder,fragm1name);
modelfolder = strcat(modelfolder,'\');
fullname = strcat(modelfolder,'Bstruct.mat');
load(fullname);

% turn 1 of fragment 2
nBonds1 = InfTable(2);
e1 = atom_pair_vector(fragm_connect(1,1), fragm_connect(1,2));
t = cross(e1,e2);
normt = norm(t);
if normt < 1e-3
    for i = 1:nBonds1
        eaux = atom_pair_vector(Bonds(i,1), Bonds(i,2));
        t = cross(eaux,e2);
        normt = norm(t);
        if normt > 1e-3,break, end
    end % i
end % if
    t = t/normt;
    phi = acos(e1*e2');
    M = matrturn(t, phi);
% Try e2 turn
e2new = M*e2';
cos_tryphi = e1*e2new;
if (sign(cos_tryphi) ~= 1)
    phi = -phi;
    M = matrturn(t, phi);
end % if
nAtoms2 = InfTable_Fragm2(1);
for i = 1:nAtoms2
    r = M*AtomCoords_Fragm2(i,:)';
    AtomCoords_Fragm2(i,:) = r';
end % i
e4 = M*e4'; % turn e4 with fragment 2

% turn 2 of fragment 2
 e3 = atom_pair_vector(fragm_connect(1,3), fragm_connect(1,4));
a = (e3*e1')*(e1*e4) - (e3*e4);
b = -(e1*(cross(e3,e4))');
phi1 = atan(b/a);
phi2 = phi1 + pi;
if (a*cos(phi1)+b*sin(phi1)) < 0
    phi = phi1;
else
    phi = phi2;
end % if
M = matrturn(e1, -phi);
for i = 1:nAtoms2
    r = M*AtomCoords_Fragm2(i,:)';
    AtomCoords_Fragm2(i,:) = r';
end % i

% shift atoms of fragment 2 to new positions
r = AtomCoords(fragm_connect(1,1),:) ...
    - AtomCoords_Fragm2(fragm_connect(2,2),:);
for i = 1:nAtoms2
    AtomCoords_Fragm2(i,:) = r + AtomCoords_Fragm2(i,:);
end % i

% Join atoms. Build Atoms2Renum table
Atoms2Renum = zeros(nAtoms2,1);
nAtoms1 = InfTable(1);
inew = nAtoms1;
for i = 1:nAtoms2
    if i == fragm_connect(2,1)
        Atoms2Renum(i) = fragm_connect(1,2);
        continue
    end
    if i == fragm_connect(2,2)
        Atoms2Renum(i) = fragm_connect(1,1);
        continue
    end
    AtomCoords = [AtomCoords;AtomCoords_Fragm2(i,:)];
    inew = inew + 1;
    Atoms2Renum(i) = inew;
    Symbol = Type_Fragm2(AtomTypes_Fragm2(i)).Symbol;
    Mass = Type_Fragm2(AtomTypes_Fragm2(i)).Mass;
    nTypes = length(Type);
    coincide = 0;
    for j = 1:nTypes
        if (Symbol == Type(j).Symbol)&&(Mass == Type(j).Mass)
            coincide = 1;
            AtomTypes = [AtomTypes;j];
            break
        end % if
    end % j
    if ~coincide
        nTypes = nTypes + 1;
        Type(nTypes).Symbol = Symbol;
        Type(nTypes).Mass = Mass;
        AtomTypes = [AtomTypes;nTypes];
    end % if
end % i

% find Fragment2 Bond with fragm_connect(2,1)&& fragm_connect(2,2)
nBonds_Fragm2 = InfTable_Fragm2(2);
ac1 = fragm_connect(2,1);
ac2 = fragm_connect(2,2);
for i = 1:nBonds_Fragm2
    a1 = Bonds_Fragm2(i,1);
    a2 = Bonds_Fragm2(i,2);
    if ((a1 == ac1)&&(a2 == ac2))||((a1 == ac2)&&(a2 == ac1)),break,end
end % i
Bonds_Fragm2(i,:) = [];

% Renumber Fragment2 vibrational coordinates
InfTable(1) = InfTable(1) + nAtoms2 - 2;
for i = 1:nBonds_Fragm2-1
    a1 = Bonds_Fragm2(i,1);
    a2 = Bonds_Fragm2(i,2);
    Bonds_Fragm2(i,1) = Atoms2Renum(a1);
    Bonds_Fragm2(i,2) = Atoms2Renum(a2);
end % i
Bonds = [Bonds;Bonds_Fragm2];
InfTable(2) = InfTable(2) + nBonds_Fragm2-1;

nAlphas_Fragm2 = InfTable_Fragm2(3);
if nAlphas_Fragm2 > 0
    for i = 1:nAlphas_Fragm2
        a1 = Alphas_Fragm2(i,1);
        a2 = Alphas_Fragm2(i,2);
        a3 = Alphas_Fragm2(i,3);
        Alphas_Fragm2(i,1) = Atoms2Renum(a1);
        Alphas_Fragm2(i,2) = Atoms2Renum(a2);
        Alphas_Fragm2(i,3) = Atoms2Renum(a3);
    end % i
    Alphas = [Alphas;Alphas_Fragm2];
    InfTable(3) = InfTable(3) + nAlphas_Fragm2;
end % if nAlphas_Fragm2 > 0

nGammas_Fragm2 = InfTable_Fragm2(4);
if nGammas_Fragm2 > 0
    for i = 1:nGammas_Fragm2
        a1 = Gammas_Fragm2(i,1);
        a2 = Gammas_Fragm2(i,2);
        a3 = Gammas_Fragm2(i,3);
        Gammas_Fragm2(i,1) = Atoms2Renum(a1);
        Gammas_Fragm2(i,2) = Atoms2Renum(a2);
        Gammas_Fragm2(i,3) = Atoms2Renum(a1);
        if (~Gammas_Fragm2(i,4))
            a5 = Gammas_Fragm2(i,5);
            a6 = Gammas_Fragm2(i,6);
            Gammas_Fragm2(i,5) = Atoms2Renum(a5);
            Gammas_Fragm2(i,6) = Atoms2Renum(a5);
        end % if
    end % i
    Gammas = [Gammas;Gammas_Fragm2];
    InfTable(4) = InfTable(4) + nGammas_Fragm2;
end % if nGammas_Fragm2 > 0

nRos_Fragm2 = InfTable_Fragm2(5);
if nRos_Fragm2 > 0
    for i = 1:nRos_Fragm2
        a1 = Ros_Fragm2(i,1);
        a2 = Ros_Fragm2(i,2);
        a3 = Ros_Fragm2(i,3);
        a4 = Ros_Fragm2(i,4);
        Ros_Fragm2(i,1) = Atoms2Renum(a1);
        Ros_Fragm2(i,2) = Atoms2Renum(a2);
        Ros_Fragm2(i,3) = Atoms2Renum(a3);
        Ros_Fragm2(i,4) = Atoms2Renum(a4);
    end % i
    Ros = [Ros;Ros_Fragm2];
    InfTable(5) = InfTable(5) + nRos_Fragm2;
end % if nRos_Fragm2 > 0

nHis_Fragm2 = InfTable_Fragm2(6);
if nHis_Fragm2 > 0
    for i = 1:nHis_Fragm2
        a1 = His_Fragm2(i,1);
        a2 = His_Fragm2(i,2);
        a3 = His_Fragm2(i,3);
        His_Fragm2(i,1) = Atoms2Renum(a1);
        His_Fragm2(i,2) = Atoms2Renum(a2);
        His_Fragm2(i,3) = Atoms2Renum(a3);
        if (~His_Fragm2(i,4))
            a5 = His_Fragm2(i,5);
            a6 = His_Fragm2(i,6);
            His_Fragm2(i,5) = Atoms2Renum(a5);
            His_Fragm2(i,6) = Atoms2Renum(a6);
        end % if
        a7 = His_Fragm2(i,7);
        a8 = His_Fragm2(i,8);
        a9 = His_Fragm2(i,9);
        His_Fragm2(i,7) = Atoms2Renum(a7);
        His_Fragm2(i,8) = Atoms2Renum(a8);
        His_Fragm2(i,9) = Atoms2Renum(a9);
        if (~His_Fragm2(i,10))
            a11 = His_Fragm2(i,11);
            a12 = His_Fragm2(i,12);
            His_Fragm2(i,11) = Atoms2Renum(a11);
            His_Fragm2(i,12) = Atoms2Renum(a12);
        end % if
    end % i
    His = [His;His_Fragm2];
    InfTable(6) = InfTable(6) + nHis_Fragm2;
end % if nHis_Fragm2 > 0

InfTable(1) = size(AtomCoords,1);
shiftMC;

InfTable(7) = InfTable(2) + InfTable(3) + InfTable(4) + ...
    InfTable(5) + InfTable(6); 
modelfolder = strcat(molfolder,jointname);
modelfolder = strcat(modelfolder,'\');

fullname = strcat(modelfolder,'Bstruct.mat');
save(fullname,'Type', 'InfTable', 'AtomTypes', 'AtomCoords', 'Bonds', ...
           'Alphas', 'Gammas', 'Ros', 'His');

to_calc_tp(modelfolder);
view_bonds;
update_ini(modelfolder,jointname);


% My sub-functions ==============================================

function [r,varargout] = atom_pair_vector(Atom1, Atom2)
% construct unit vector from Atom1 to Atom2 using AtomCoords
global AtomCoords

r = AtomCoords(Atom2,:) - AtomCoords(Atom1,:);
s = norm(r);
r = r/s;
if nargout > 0, varargout{1} = s; end % if

% ---------------------------------------
function M = matrturn(ax, phirad)
% M is result - 3x3-matrix of rotation
% ax is unit vector of axis of turn
% phirad is signed angle given radians
if abs(norm(ax) - 1) > 1e-3
    error('axis must be of unit size')
end % if
cs = cos(phirad);
sn = sin(phirad);
M = zeros(3);
M(1,1) = ax(1)^2*(1 - cs) + cs;
M(1,2) = ax(1)*ax(2)*(1 - cs) - ax(3)*sn;
M(1,3) = ax(1)*ax(3)*(1 - cs) + ax(2)*sn; 
M(2,1) = ax(2)*ax(1)*(1 - cs) + ax(3)*sn;
M(2,2) = ax(2)^2*(1 - cs) + cs;
M(2,3) = ax(2)*ax(3)*(1 - cs) - ax(1)*sn;
M(3,1) = ax(3)*ax(1)*(1 - cs) - ax(2)*sn;
M(3,2) = ax(3)*ax(2)*(1 - cs) + ax(1)*sn;
M(3,3) = ax(3)^2*(1 - cs) + cs;

% ----------------------------------------
function to_calc_tp(modelfolder)
% To calculi and viewer save txt-files

global Type AtomTypes InfTable AtomCoords Bonds 
global Alphas Gammas Ros His

Vectors = zeros(6,3);
Vectors(1,1) = 1;   % ort x
Vectors(2,1) = -1;  % ort -x
Vectors(3,2) = 1;   % ort y
Vectors(4,2) = -1;  % ort -y
Vectors(5,3) = 1;   % ort z
Vectors(6,3) = -1;   % ort -z

nAtoms = InfTable(1);
nBonds = InfTable(2);
nAlphas = InfTable(3);
nGammas = InfTable(4);
nRos = InfTable(5);
nHis = InfTable(6);

fullname = strcat(modelfolder,'InfTable.txt');
fid = fopen(fullname,'w');
    fprintf(fid,'%12i%12i%12i%12i%12i%12i \r\n',InfTable(1:6));
    fprintf(fid,'%12i%12i \r\n',InfTable(7:8));
fclose(fid);

fullname = strcat(modelfolder,'AtomCoords.txt');
fid = fopen(fullname,'w');
for i = 1:nAtoms
    fprintf(fid,'%11.4f%11.4f%11.4f \r\n',AtomCoords(i,:));
end % i
fclose(fid);

fullname = strcat(modelfolder,'AtomMasses.txt');
fid = fopen(fullname,'w');
for i = 1:nAtoms
    mass = Type(AtomTypes(i)).Mass;
    fprintf(fid,'%11.4f \r\n', mass);
end % i
fclose(fid);

fullname = strcat(modelfolder,'AtomSymbols.txt');
fid = fopen(fullname,'w');
for i = 1:nAtoms
    symbol = Type(AtomTypes(i)).Symbol;
    fprintf(fid,'%s \r\n', symbol);
end % i
fclose(fid);

BondUnitVectors = zeros(nBonds,3);
BondLengthes = zeros(nBonds,1);
for i = 1:nBonds
    Atom1 = Bonds(i,1);
    Atom2 = Bonds(i,2);
    [r,s] = atom_pair_vector(Atom1, Atom2);
    BondUnitVectors(i,:) = r / norm(r);
    BondLengthes(i) = s;
end % i

fullname = strcat(modelfolder,'BondVectors.txt');
fid = fopen(fullname,'w');
for i = 1:nBonds
    fprintf(fid,'%11.4f%11.4f%11.4f \n', BondUnitVectors(i,:));
end % i
fclose(fid);

fullname = strcat(modelfolder,'BondLengthes.txt');
fid = fopen(fullname,'w');
for i = 1:nBonds
    fprintf(fid,'%11.4f \n', BondLengthes(i));
end % i
fclose(fid);

fullname = strcat(modelfolder,'BqStructures.txt');
fid = fopen(fullname,'w');
for i = 1:nBonds
    fprintf(fid,'%4i%4i%11.4f%11.4f%11.4f%11.4f%11.4f%11.4f \r\n', ...
        Bonds(i, :), -BondUnitVectors(i,:), BondUnitVectors(i,:));
end % i
fclose(fid);

if nAlphas > 0
    fullname = strcat(modelfolder,'BaStructures.txt');
    fid = fopen(fullname,'w');
    for ia = 1:nAlphas
        Atom1 = Alphas(ia, 1);
        Atom2 = Alphas(ia, 2);
        Atom3 = Alphas(ia, 3);
        [e1,s1] = atom_pair_vector(Atom2, Atom1);
        s1 = 1.09/s1; 
        [e2,s2] = atom_pair_vector(Atom2, Atom3);
        s2 = 1.09/s2; 
        cs = e1*e2';
        sn = sqrt(1 - cs^2);
        b1 = (e2 - e1*cs)*s1/sn;
        b3 = (e1 - e2*cs)*s2/sn;
        b2 = b1 + b3;
        b1 = -b1;
        b3 = -b3;
        fprintf(fid,'%4i%4i%4i', Atom1, Atom2, Atom3);
        fprintf(fid,'%11.4f%11.4f%11.4f', b1(:));
        fprintf(fid,'%11.4f%11.4f%11.4f', b2(:));
        fprintf(fid,'%11.4f%11.4f%11.4f \n', b3(:));
    end % ia
    fclose(fid);
end % if nAlphas > 0

if nGammas > 0
    fullname = strcat(modelfolder,'BgStructures.txt');
    fid = fopen(fullname,'w');
    for ig = 1:nGammas
        Atom1 = Gammas(ig, 1);
        Atom2 = Gammas(ig, 2);
        Atom3 = Gammas(ig, 3);
        [e1,s1] = atom_pair_vector(Atom2, Atom1);
        s1 = 1.09/s1; 
        [e2,s2] = atom_pair_vector(Atom2, Atom3);
        s2 = 1.09/s2;
        gVectorID = Gammas(ig, 4);
        if gVectorID > 0
            gVector = Vectors(gVectorID,:);
        else
            eaux = atom_pair_vector(Gammas(ig, 5), Gammas(ig, 6));
            normalaux = cross(e1, eaux);
            normalaux = normalaux/norm(normalaux);
            gVector = cross(normalaux, e1);
        end % if gVectorID > 0
        b1 = gVector*s1;
        b3 = gVector*s2;
        b2 = b1 + b3;
        b1 = -b1;
        b3 = -b3;
        fprintf(fid,'%4i%4i%4i', Atom1, Atom2, Atom3);
        fprintf(fid,'%11.4f%11.4f%11.4f', b1(:));
        fprintf(fid,'%11.4f%11.4f%11.4f', b2(:));
        fprintf(fid,'%11.4f%11.4f%11.4f \n', b3(:));
    end % ig
    fclose(fid);
end % if nGammas > 0

if nRos > 0
    fullname = strcat(modelfolder,'BrStructures.txt');
    fid = fopen(fullname,'w');
    for ir = 1:nRos
        Atom1 = Ros(ir, 1);
        Atom2 = Ros(ir, 2);
        Atom3 = Ros(ir, 3);
        Atom4 = Ros(ir, 4);
        [e1, s1] = atom_pair_vector(Atom2, Atom1);
        [e2, s2] = atom_pair_vector(Atom2, Atom3);
        [e4, s4] = atom_pair_vector(Atom2, Atom4);
        s1 = 1.09/s1;
        s2 = 1.09/s2;
        s4 = 1.09/s4;
        eAux = cross(e1,e2);
        s = norm(eAux);
        eAux = eAux/s; % it is normal to the Plane
        cs = e4*eAux';
        sn = sqrt(1 - cs^2);
        Auxr1 = (e4 - eAux*cs)/sn; % it is g-vector
        Auxr2 = (eAux - e4*cs)/sn; % it is f-vector
        temp = cross(e2, Auxr1);
        e3 = cross(e1, Auxr1);
        b1 = s1*temp/s;
        b3 = -s2 * e3/s;
        b4 = s4*Auxr2;
        b2 = -(b4 + b1 + b3);
        fprintf(fid,'%4i%4i%4i%4i', Atom1, Atom2, Atom3, Atom4);
        fprintf(fid,'%11.4f%11.4f%11.4f', b1(:));
        fprintf(fid,'%11.4f%11.4f%11.4f', b2(:));
        fprintf(fid,'%11.4f%11.4f%11.4f', b3(:));
        fprintf(fid,'%11.4f%11.4f%11.4f \n', b4(:));
    end % ir
    fclose(fid);
end % nRos > 0

if nHis > 0
    fullname = strcat(modelfolder,'BhStructures.txt');
    fid = fopen(fullname,'w');
    for ih = 1:nHis
        % Plane 1
        Atom1 = His(ih, 1);
        Atom2 = His(ih, 2);
        Atom3 = His(ih, 3);
        [e1, s1] = atom_pair_vector(Atom2, Atom1);
        [e2, s2] = atom_pair_vector(Atom2, Atom3);
        s1 = 1.09/s1;
        s2 = 1.09/s2;
            % Get eaux as the turn axis
        i = His(ih, 4);
        if i > 0    % axis is a ready unit vector
            eaux = Vectors(i,:);
        else
            eaux = atom_pair_vector(His(ih, 5), His(ih, 6));
        end % if > 0
            % Get normal
        normal = cross(e1, e2);
        normal = normal/norm(normal);
        ete1 = eaux*e1';
        ete2 = eaux*e2';
        e1e2 = e1*e2';
        A = sqrt(ete1^2 + ete2^2 -2*e1e2*ete1*ete2);
        b1 = s1*normal*ete2/A;
        b3 = -s2*normal*ete1/A;
        b2 = -(b1 + b3);
        % Plane 2
        Atom4 = His(ih, 1+6);
        Atom5 = His(ih, 2+6);
        Atom6 = His(ih, 3+6);
        [e1, s1] = atom_pair_vector(Atom5, Atom4);
        [e2, s2] = atom_pair_vector(Atom5, Atom6);
        s1 = 1.09/s1;
        s2 = 1.09/s2;
            % Get eaux as the turn axis
        i = His(ih, 4+6);
        if i > 0    % axis is a ready unit vector
            eaux = Vectors(i,:);
        else
            eaux = atom_pair_vector(His(ih, 5+6), His(ih, 6+6));
        end % if  > 0
            % Get normal
        normal = cross(e1, e2);
        normal = normal/norm(normal);
        ete1 = eaux*e1';
        ete2 = eaux*e2';
        e1e2 = e1*e2';
        A = sqrt(ete1^2 + ete2^2 -2*e1e2*ete1*ete2);
        b4 = s1*normal*ete2/A;
        b6 = -s2*normal*ete1/A;
        b5 = -(b4 + b6);
        fprintf(fid,'%4i%4i%4i%4i%4i%4i', ...
            Atom1, Atom2, Atom3, Atom4, Atom5, Atom6);
        fprintf(fid,'%11.4f%11.4f%11.4f', b1(:));
        fprintf(fid,'%11.4f%11.4f%11.4f', b2(:));
        fprintf(fid,'%11.4f%11.4f%11.4f', b3(:));
        fprintf(fid,'%11.4f%11.4f%11.4f', b4(:));
        fprintf(fid,'%11.4f%11.4f%11.4f', b5(:));
        fprintf(fid,'%11.4f%11.4f%11.4f \n', b6(:));
    end % ih
    fclose(fid);
end % if nHis > 0

% ----------------------------------
function view_bonds
% View Model's bonds within MatLab

global InfTable AtomCoords Bonds 
hold off
Coords = zeros(2,3);
nBonds = InfTable(2);
for iBond = 1:nBonds
    a1 = Bonds(iBond,1);
    a2 = Bonds(iBond,2);
    Coords(1,:) = AtomCoords(a1,:);
    Coords(2,:) = AtomCoords(a2,:);
    plot3(Coords(1:2,1),Coords(1:2,2),Coords(1:2,3),'b-o','MarkerSize',12)
    hold on
end % iBond
xlabel('x')
ylabel('y')
zlabel('z')
grid on

nAtoms = InfTable(1);
for i = 1:nAtoms
    num = int2str(i);
    r = AtomCoords(i,:);
    text('String',num, 'Position',r);
end % i Atoms

% ------------------------------

function update_ini(pname,jointname)
diskname = pname(1:2);
filename = strcat(diskname,'\LEVINFINITE\Config\ToCalculate.ini');
fid = fopen(filename,'w');
fprintf(fid,'%s\n ',jointname);
fclose(fid);


function shiftMC()
global Type InfTable AtomTypes AtomCoords 
nAtoms = InfTable(1);
rmc = zeros(1,3);
cmass = 0;
for i = 1:nAtoms
    r = AtomCoords(i,:);
    it = AtomTypes(i);
    m = Type(it).Mass;
    cmass = cmass + m;
    rmc = rmc + m*r; 
end % i Atoms
rmc = rmc/cmass;
rmc = - rmc;
 for i = 1:nAtoms
    AtomCoords(i,:) = AtomCoords(i,:) + rmc;
end % i Atoms

 
% End My functions-------------------------------------
