 figure(101)


 subplot(4,2,6) % histogram of all pitches mod 1
 hold off
 histogram(mod(p(:),1)*12,(0:.005:1)*12)
 xlabel('pitch (x mod 1)')

 subplot(4,2,7) %plot max of pitch histogram (useful for tracking thermalization of a perturbation)
 hold off
 hp=histcounts(mod(p(:),1)*12,(0:.005:1)*12);
 mp=max(hp);
 mps=[mps mp];
 plot(mps)

%  subplot(4,2,7) %pitch-pitch correlations along x and y and z directions
%  (interesting to track emergence of order, though not really used for
%  anything)
% corrs=[];
% for shft=1:n
%     pxshft=circshift(p,[shft, 0, 0]);
%     pyshft=circshift(p,[0, shft, 0]);
%     pzshft=circshift(p,[0, 0, shft]);
%     dx=mod(p-pxshft,1);
%     dy=mod(p-pyshft,1);
%     dz=mod(p-pzshft,1);
%     hax=0:.002:1;
%     %corr=histcounts([dy(:)],hax);
%     corr=histcounts([dz(:) dy(:) dx(:)],hax);
%     corrs=[corrs; corr];
% end
% mx=1*max(max(corrs(2:(end-1),:)));
% mn=1*min(min(corrs(2:(end-1),:)));
% 
% imagesc(hax,1:n,corrs,[mn mx])
% ylabel('dx, dy')
% xlabel('interval (x)')
% colormap('jet') 
% 
% subplot(4,2,8) % histogram of all pitches
% plot(hax(2:end),corrs(1,:))
% ylim([0 600])
% xlabel('Interval (x)')

subplot(4,2,[1 2]) % total energy
if mod(coun1,1)==0

up=0;
for jj=1:length(nnx)/2
    pshft=circshift(p,[nnx(jj) nny(jj) nnz(jj)]);
    dx=mod(p-pshft,1);
    u1=interp1(faxn,dn,dx,'linear',0);
    up=up+sum(u1(:));
end
tp=0.5*m*v.^2;

 utot(coun1)=up+sum(tp(:)); 
   
hold off
plot(utot)
xlabel('iteration')
ylabel('energy')
xlim([0 coun1])
end

%calculate order parameters
i=sqrt(-1);

for jj=[5 7 12] %which ones to calculate
    pk=exp(jj*i*(2*pi*p-pi/2));
    mk=sum(pk(:))/64000;
    mks(jj,coun1)=mk; 
end

coun1=coun1+1;
  
%plot order parameters  
subplot(4,2,[3 4])
hold off
for jj=[5 7 12]   
    plot(abs(mks(jj,:)))
    hold on
end
legend({'5' '7' '12'})
xlim([0 coun1])

subplot(4,2,5) %image plot of a slice of pitches with mod 1 colorscale
c=colormap('jet');

sclfac=2.5;
p0=squeeze(p(:,:,20));
fac=(floor(p0)/sclfac);
pr=interp1(1:256,c(:,1),ceil(mod(p0,1)*256));
pg=interp1(1:256,c(:,2),ceil(mod(p0,1)*256));
pb=interp1(1:256,c(:,3),ceil(mod(p0,1)*256));
phsv=rgb2hsv(cat(3,pr, pg, pb));
phue=squeeze(phsv(:,:,1));
psat=squeeze(phsv(:,:,2)).*min(2.^fac,1);
pval=squeeze(phsv(:,:,3)).*min(2.^(-fac),1);
pim=hsv2rgb(cat(3,phue, pval, psat));

image(pim)
daspect([1 1 1])

drawnow
