clear

% Import the files
DELIMITER = ' ';
HEADERLINES = 0;
ExpSolCoulomb = importdata('./Solv/wwd.out.force.xyz', DELIMITER, HEADERLINES);
ImpSolCoulomb = importdata('./Impl/Coulomb/wwd.out.force.xyz', DELIMITER, HEADERLINES);
ImpSolCoulombDiElecEPS1D78S03 = importdata('./Impl/CoulombDiElecEPS1/wwd.out.force.xyz', DELIMITER, HEADERLINES);
ImpSolCoulombDiElecEPS1D53S02 = importdata('./Impl/CoulombDiElecEPS1D53S02/wwd.out.force.xyz', DELIMITER, HEADERLINES);
ImpSolCoulombDiElecEPS2D53S15 = importdata('./Impl/CoulombDiElecEPS2D53S15/wwd.out.force.xyz', DELIMITER, HEADERLINES);
masses = importdata('./mass.txt');
sqmasses = sqrt(masses);
sqmassesInv = 1 ./ sqmasses;

numatoms = size(ExpSolCoulomb,1);

% matrix of differences
difdata = zeros(numatoms,12);
difdataMnorm = zeros(numatoms,12);
difdataVnorm = zeros(numatoms,4);
difdataVnormMnorm = zeros(numatoms,4);

%difXYZ
difdata(:,1) = abs(ExpSolCoulomb(:,1) - ImpSolCoulomb(:,1));
difdata(:,2) = abs(ExpSolCoulomb(:,2) - ImpSolCoulomb(:,2));
difdata(:,3) = abs(ExpSolCoulomb(:,3) - ImpSolCoulomb(:,3));

%difXYZ
difdata(:,4) = abs(ExpSolCoulomb(:,1) - ImpSolCoulombDiElecEPS1D78S03(:,1));
difdata(:,5) = abs(ExpSolCoulomb(:,2) - ImpSolCoulombDiElecEPS1D78S03(:,2));
difdata(:,6) = abs(ExpSolCoulomb(:,3) - ImpSolCoulombDiElecEPS1D78S03(:,3));

%difXYZ
difdata(:,7) = abs(ExpSolCoulomb(:,1) - ImpSolCoulombDiElecEPS1D53S02(:,1));
difdata(:,8) = abs(ExpSolCoulomb(:,2) - ImpSolCoulombDiElecEPS1D53S02(:,2));
difdata(:,9) = abs(ExpSolCoulomb(:,3) - ImpSolCoulombDiElecEPS1D53S02(:,3));

%difXYZ
difdata(:,10) = abs(ExpSolCoulomb(:,1) - ImpSolCoulombDiElecEPS2D53S15(:,1));
difdata(:,11) = abs(ExpSolCoulomb(:,2) - ImpSolCoulombDiElecEPS2D53S15(:,2));
difdata(:,12) = abs(ExpSolCoulomb(:,3) - ImpSolCoulombDiElecEPS2D53S15(:,3));


% metrics for raw forces
difdataMEAN = mean(difdata);
difdataSTD = std(difdata);
difdataMAX = max(difdata);

% metrics for mass weighted forces by direction
% for i=1:12
%     difdataMnorm(:,i) = sqmassesInv .* difdata(:,i);
%     if i<4
%         ExpSolCoulombMnorm(:,i) = sqmassesInv .* abs(ExpSolCoulomb(:,i));
%     end
% end
% difdataMnormMEAN = mean(difdataMnorm);
% difdataMnormSTD = std(difdataMnorm);
% difdataMnormMAX = max(difdataMnorm);
% ExpSolCoulombMnormMEAN = [mean(ExpSolCoulombMnorm(:,1)),mean(ExpSolCoulombMnorm(:,2)),mean(ExpSolCoulombMnorm(:,3)),...
%                     mean(ExpSolCoulombMnorm(:,1)),mean(ExpSolCoulombMnorm(:,2)),mean(ExpSolCoulombMnorm(:,3)),...
%                     mean(ExpSolCoulombMnorm(:,1)),mean(ExpSolCoulombMnorm(:,2)),mean(ExpSolCoulombMnorm(:,3)),...
%                     mean(ExpSolCoulombMnorm(:,1)),mean(ExpSolCoulombMnorm(:,2)),mean(ExpSolCoulombMnorm(:,3))];
% avgFE = difdataMnormMEAN ./ ExpSolCoulombMnormMEAN                
% maxFE = difdataMnormMAX ./ ExpSolCoulombMnormMEAN

% metrics for mass weighted forces by vector norm
for i=1:4
    j=i-1;
    difdataVnorm(:,i) = sqrt(difdata(:,3*j+1) .^2 + difdata(:,3*j+2) .^2 + difdata(:,3*j+3) .^2);
    difdataVnormMnorm(:,i) = sqmassesInv .* difdataVnorm(:,i);
end

avgVnorm = mean(difdataVnorm)
maxVnorm = max(difdataVnorm)

ExpSolCoulombVnorm(:,1) = sqrt(ExpSolCoulomb(:,1) .^2 + ExpSolCoulomb(:,2) .^2 + ExpSolCoulomb(:,3) .^2);
ExpSolCoulombVnormMnorm(:,1) = sqmassesInv .* ExpSolCoulombVnorm(:,1);
ExpSolCoulombVnormMnormMEAN = [mean(ExpSolCoulombVnormMnorm(:,1)),...
                    mean(ExpSolCoulombVnormMnorm(:,1)),...
                    mean(ExpSolCoulombVnormMnorm(:,1)),...
                    mean(ExpSolCoulombVnormMnorm(:,1))];

avgVnormMnorm = mean(difdataVnormMnorm) ./ ExpSolCoulombVnormMnormMEAN              
maxVnormMnorm = max(difdataVnormMnorm) ./ ExpSolCoulombVnormMnormMEAN

%figure(2)
%plot(ExpSolCoulomb(:,1))
%hold all
%plot(ImpSolCoulomb(:,1))
%plot(ImpSolCoulombDiElecEPS1D78S03(:,1))
%plot(ImpSolCoulombDiElecEPS1D53S02(:,1))
%plot(ImpSolCoulombDiElecEPS2D53S15(:,1))
