function [ H,S] = calcHamiltonian( orb,Z,xAtom,yAtom,zAtom )
%CALCHAMILTONIAN calculates Hamiltonian and overlap matrix
% Uses molecular orbitals and doesn't take nuclei-nuclei interaction into
% account. 
%
% Z ... a [m x 1] matrix containing the atomic numbers of all atoms in the
%       molecule. Atoms with n orbitals have to mentioned n times
% xAtom, yAtom and zAtom are [m x 1] matrizes containing relative positions
%       of all molecules (in Meter). 
% orb ... {m x 1} cell containing {1 x n} cells containing strings
%       describing the n orbitals which
%       should be used for calculation (1s, 2s, 2px, 2py, 2pz ...).
%
%H is a matrix containing energy expectation values, S is the overal matrix

a0 = 5.2917721E-11; % Bohr radius in meters
%  
%  x_lim = a0*[-300 300];
%  y_lim = a0*[-300 300];
%  z_lim = a0*[-300 300];

universallimit = 50*a0;

x_lim_l = -universallimit;
x_lim_u = universallimit;

y_lim_l = -universallimit;
y_lim_u = universallimit;

z_lim_l = -universallimit;
z_lim_u = universallimit;

hbar = 1.05457266E-34;
me = 9.1093897E-31;
eps0 = 8.854187817E-12;
ec = 1.60217733E-19;
k1 = 1/2*hbar^2/me;
k2 = ec^2/(4*pi*eps0);

%---------------------Set tolerance of integrals---------------------------
ATol = ec/1E2;
RTol = 1E-5;

num_orbs_per_atom = cellfun('size',orb,2);
nopa_cumsum = cumsum(num_orbs_per_atom);
N_atom = length(Z);
N = sum(num_orbs_per_atom);
orblist = {};
for k = 1:length(orb)
    orblist = [orblist orb{k}];
end

H = NaN(N);
S = H;



for k = 1:N
    for l = k:N
        atom_k = find(k >= nopa_cumsum,1,'last');
        atom_l = find(l >= nopa_cumsum,1,'last');
   
        if(k==l) S(k,l) = 1;
        else
         S(k,l) = integral3(@(x,y,z) atorbit(x-xAtom(atom_k),y- yAtom(atom_k), z-zAtom(atom_k),Z(atom_k),false,orblist{k}).* ...
            atorbit(x-xAtom(atom_l),y- yAtom(atom_l),z- zAtom(atom_l),Z(atom_l),false,orblist{l}), ...
            x_lim_l,x_lim_u,y_lim_l,y_lim_u,z_lim_l,z_lim_u,'AbsTol',ATol,'RelTol',RTol);
        end
        fprintf('Done S(%i,%i)\n',k,l);
    end
end
S(tril(true(N),-1)) = S(triu(true(N),1));
 
for k = 1:N
    for l = 1:N
        atom_k = find(k >= nopa_cumsum,1,'last');
        atom_l = find(l >= nopa_cumsum,1,'last');



        orb1 = orblist{k};
        orb2 = orblist{l};

        pos1 = [xAtom(atom_k) yAtom(atom_k) zAtom(atom_k)];
        pos2 =  [xAtom(atom_l) yAtom(atom_l) zAtom(atom_l)];


        %n1 = str2double(orb1(1));
        n2 = str2double(orb2(1));

        %eigenE1 = -me*ec^4/(2*(4*pi*eps0)^2*hbar^2)/n1^2*Z1^2;
        eigenE2 = -me*ec^4/(2*(4*pi*eps0)^2*hbar^2)/n2^2*Z(atom_l)^2;

        H(k,l) = eigenE2*S(k,l);
        
        % Exclude atom corresponding to right eigenfunction (atom_l)
        for j = [1:atom_l-1 atom_l+1:N_atom]
            H(k,l) = H(k,l) - k2*Z(j)*...
            integral3(@(x,y,z) atorbit(x-pos1(1),y-pos1(2),z-pos1(3),Z(atom_k),false,orb1)...
                ./sqrt((x-xAtom(j)).^2+(y-yAtom(j)).^2+(z-zAtom(j)).^2).* ...
                atorbit(x-pos2(1),y-pos2(2),z-pos2(3),Z(atom_l),false,orb2),...
                x_lim_l,x_lim_u,y_lim_l,y_lim_u,z_lim_l,z_lim_u,'AbsTol',ATol,'RelTol',RTol);
            fprintf('Done H(%i,%i) - Atom %i\n',k,l,j);
        end
        
        %
        % mask = 1:N_atom ~= atom_l;
        % H(k,l) = H(k,l) - k2*...
        %         integral3(@(x,y,z) expecval(x,y,z,xAtom(mask),yAtom(mask),zAtom(mask),Z(mask),orb1,orb2,pos1,pos2,Z(atom_k),Z(atom_l)),...
        %         x_lim_l,x_lim_u,y_lim_l,y_lim_u,z_lim_l,z_lim_u,'AbsTol',1e-20,'RelTol',1e-5);
    end
end


 % integral3(@(x,y,z) atorbit(x,y,z,0,0,0,1,false,'s2'),-1,1,-1,1,-1,1);


end

function H = expecval(x,y,z,xAtom,yAtom,zAtom,Z,orb1,orb2,pos1,pos2,Z1,Z2)
%atom at pos2 has to be excluded from atommatrix xAtom,yAtom etc.
% 
%   hbar = 1.05457266E-34;
%  me = 9.1093897E-31;
%  eps0 = 8.854187817E-12;
%  ec = 1.60217733E-19;
%  %k1 = 1/2*hbar^2/me;
%  k2 = ec^2/(4*pi*eps0);
%  
%  %n1 = str2double(orb1(1));
%  n2 = str2double(orb2(1));
%  
%  %eigenE1 = -me*ec^4/(2*(4*pi*eps0)^2*hbar^2)/n1^2*Z1^2;
%  eigenE2 = -me*ec^4/(2*(4*pi*eps0)^2*hbar^2)/n2^2*Z2^2;
%  
% H = eigenE2*atorbit(x-pos1(1),y-pos1(2),z-pos1(3),Z1,false,orb1).*atorbit(x-pos2(1),y-pos2(2),z-pos2(3),Z2,false,orb2);

H=0;
for j = 1:length(Z)
    H = H - atorbit(x-pos1(1),y-pos1(2),z-pos1(3),Z1,false,orb1).* ...
        Z(j)./sqrt((x-xAtom(j)).^2+(y-yAtom(j)).^2+(z-zAtom(j)).^2).* ...
        atorbit(x-pos2(1),y-pos2(2),z-pos2(3),Z2,false,orb2);
end

end