% This Matlab Script Calculates a Ramachandran Distribution 
function [Bins]=RamachandranFe(fn,maxfe)
%clear
numTestSets = 1; % adjust the terminal iteration to match number of test sets

for i=1:numTestSets
    index = num2str(i);
    dihedraldata = load(fn);
    % WW ARG12 are dihedrals 32 and 33 for matlab +2 to column index
    Phis(:,i) = dihedraldata(:,1); % dihedral index 10 in ProtoMol
    Psis(:,i) = dihedraldata(:,2); % dihedral index 17 in ProtoMol
end

dataSize = size(Phis(:,1));
binRes =  5; % Bin resolution in degrees

PhiBins = -180:binRes:180;
PsiBins = -180:binRes:180;
tbinCount = size(PhiBins);
binCount = tbinCount(2);
Bins = zeros(binCount,binCount);

confFoundA = 0;

for j=1:dataSize
    for k=1:numTestSets
        
        % added to shift axis to match published results
        while Phis(j,k) >= 180
            Phis(j,k) = Phis(j,k) - 360;
        end
        while Psis(j,k) >= 180
            Psis(j,k) = Psis(j,k) - 360;
        end
        while Phis(j,k) < -180
            Phis(j,k) = Phis(j,k) + 360;
        end
        while Psis(j,k) < -180
            Psis(j,k) = Psis(j,k) + 360;
        end
        % end shift
        
        curPhiBin = floor((Phis(j,k)+180)/binRes)+1;
        curPsiBin = floor((Psis(j,k)+180)/binRes)+1;
        Bins(curPsiBin,curPhiBin) = Bins(curPsiBin,curPhiBin) + 1;
    end
end

% confAFoundat;
% confBFoundat; 

Bins;
Bins = Bins ./ sum(sum(Bins));
for i=1:binCount
    for j=1:binCount
        if Bins(i,j) > 0
            Bins(i,j) = -log(Bins(i,j));%0.6*
        end
    end
end

maxb = max(max(Bins));
if exist('maxfe') == 1
    maxb = maxfe;
end

for i=1:binCount
    for j=1:binCount
        if Bins(i,j) == 0
            Bins(i,j) = maxb;
        end
        if Bins(i,j) > maxb;
            Bins(i,j) = maxb;
        end
    end
end
%flip y data if imagesc
Binsout=Bins;
for i=1:(binCount/2)
    for j=1:binCount
       temp = Binsout(i,j);
       Binsout(i,j) = Binsout(binCount-i+1,j);
       Binsout(binCount-i+1,j) = temp;
    end
end

figure(3);
load('MyColormaps','mycmap');
set(gcf,'Colormap',mycmap);
imagesc(PhiBins,PsiBins,Binsout)
%mesh(PhiBins,PsiBins,Bins);
hold all
set(gca,'FontSize',17,...
  'Layer','top',...
  'XTick',[-150 -100 -50 0 50 100 150 ],...
  'XTickLabel',{'-150','-100','-50','0','50','100','150'},...
  'YTick',[-150 -100 -50 0 50 100 150 ],...
  'YTickLabel',{'150','100','50','0','-50','-100','-150'})
% annotation1 = annotation(...
%   gcf,'line',...
%   [0.0833 0.0],[0.0833 0.3333],[0.5 0.0],[0.5 0.3333],...
%   'LineWidth',2);
colorbar1 = colorbar('peer',...
  gca,[0.9225 0.131 0.04771 0.794],...
  'Box','on',...
  'Location','manual');xlabel('\Phi')
xlim([-182.5 177.5])
ylim([-177.5 182.5])
ylabel('\Psi')
plot([-185 -60 -60],[150 150 185],[-185 -60 -60],[-110 -110 -180],...
    [-110 -110 -60 -60],[-110 -45 -45 -110],...
    [-110 -110 -60 -60 -110],[35 70 70 35 35],...
    [160 160 180],[185 150 150],...
    [160 160 180],[-180 -110 -110],...
    'Color',[0 0 0],...
    'LineWidth',3)
annotation1 = annotation(...
  gcf,'textarrow',...
  [0.4857 0.4214],[0.5 0.4595],...
  'LineWidth',2,...
  'String',{'\alpha_R'},...
  'FontName','Times New Roman',...
  'FontSize',17,...
  'TextLineWidth',2);
annotation2 = annotation(...
  gcf,'textarrow',...
  [0.4821 0.4268],[0.619 0.631],...
  'LineWidth',2,...
  'String',{'C7_{eq}'},...
  'FontName','Times New Roman',...
  'FontSize',17,...
  'TextLineWidth',2);
annotation3 = annotation(...
  gcf,'textarrow',...
  [0.4911 0.425],[0.8357 0.8571],...
  'LineWidth',2,...
  'String',{'C5_{ax}'},...
  'FontName','Times New Roman',...
  'FontSize',17,...
  'TextLineWidth',2);
annotation4 = annotation(...
  gcf,'textarrow',...
  [0.4857 0.4232],[0.231 0.2024],...
  'LineWidth',2,...
  'String',{'C5_{ax}'},...
  'FontName','Times New Roman',...
  'FontSize',17,...
  'TextLineWidth',2);
annotation5 = annotation(...
  gcf,'arrow',...
  [0.5714 0.8589],[0.8429 0.8452],...
  'LineWidth',2);
hold off

% figure(4);
% mb = max(max(Bins))
% minb = min(min(Bins))
% %mb=0.0091
% contVals = minb:(mb-minb)/500:mb;
% %Bins = Bins / 6.15385
% contour(PhiBins,PsiBins,Bins,contVals);
% hold all
% xlabel('\Phi')
% xlim([-180 180])
% ylim([-180 180])
% ylabel('\Psi')
% title('Implicitly Solvated Alanine - 300K')
% %print -depsc '2Dalanine_ramaNM.eps'
% %print -depsc '2DimpSolWW_rama.png'
% hold off