function [h, histmat] = hist2dplot(xdata,ydata,xmin,ymin,nx,ny,dx,dy,contSwi)
    %xname,yname,xeinheit,yeinheit)
% Zeichnen eines zweidimensionalen Histogramms
% INPUT: 
%  x(y)data : Zeitreihe
%  x(y)min  : Untere Schrank des Histogramms
%  nx(y)    : Anzahl der Intervalle des Histogramms
%  dx(y)    : Intervallbreite
%  x(y)name : Bezeichnung der x(y)-Variable
%  x(y)einheit: Einheit der x(y)-Variable

%% Vorbereitungen
xmax=xmin+nx*dx; % obere Schranke im Histogramm
ymax=ymin+ny*dy; %         "

% Berechnung des 2d-Histogramms
ix = floor((xdata-xmin)/(xmax-xmin)*nx);
iy = floor((ydata-ymin)/(ymax-ymin)*ny);
index=~isnan(ix);ix(index) = max(0,min(nx-1,ix(index)));
index=~isnan(iy);iy(index) = max(0,min(ny-1,iy(index)));
iz= ix*ny+iy;
histmat=reshape(hist(iz,0:(nx*ny-1)),ny,nx);
histmat=histmat/sum(sum(histmat))/dx/dy; % Normierung
histmat=smooth2(histmat,5);
%% Plotten des 2d-Histogramms
hold off
colormap([ 1 1 1 ; jet(127)])
xE = [xmin+0.5*dx xmax-0.5*dx];
yE = [ymin+0.5*dy ymax-0.5*dy];
xs = xmin:dx:xmax-dx;    % center of x ticks
ys = ymin:dy:ymax-dy;    % center of x ticks
xc = (xmin+dx/2):dx:(xmax-dx/2);    % center of x ticks
yc = (ymin+dy/2):dy:(ymax-dy/2);    % center of x ticks

if(contSwi)
    [~, h] = contourf(xc,yc,histmat,9);
    h.LineColor = 'none';
else
    h = spcolor(xc,yc,histmat);
    axis('tight');
    axis([xmin xmax ymin ymax]);
    shading('flat');
end



xlim([xmin xmax]);
ylim([ymin ymax]);
set(gca,'YDir','normal')
%colorbar
%xlabel([xname ' [' xeinheit ']']);
%ylabel([yname ' [' yeinheit ']']);
%text(0.9,-0.1,['Haeufigkeit [1/(' xeinheit '*' yeinheit ')]'],...
%    'units','normalized');



% This guys is going to smooth it for 2 dimensions I hope;
end

function h = spcolor(x,y,A)

%%

[S1, S2] = size(A);

A2 = [A zeros(S1,1);
    zeros(1,S2) 0];

xx = [x 2*x(end)-x(end-1)]-(x(end)-x(end-1))/2;
yy = [y 2*y(end)-y(end-1)]-(y(end)-y(end-1))/2;

h = pcolor(xx,yy,A2);

end



function yout = smooth2(A,n)
if length(A)<3
    warning('Sorry bud, you can''t smooth less than 3 points, thats silly');
    yout = A;
    return;
end

if length(n)==1
    if n<2
        yout = A;
        return;
    end
    
%     forcing you to have an odd 'n'
    if double(~logical(round(n/2)-n/2))
        n = n+1;
    end
    
    bee = linspace(-1.96,1.96,n); % normal distribution with 95% confidence bounds
    [BX, BY] = meshgrid(bee);
    R2 = BX.^2+BY.^2;
    toConvolve = exp(-R2)/sum(exp(-R2(:)));
    
%     toConvolve = exp(-bee.^2)/sum(exp(-bee.^2));
%     old moving average computation: 
% toConvolve = ones(n,1)/n;
else
    toConvolve = n;
    if round(sum(toConvolve)*100000)/100000~=1
        warning('the sum here does not equal to one.');
    end
    n = length(toConvolve);
end

if min(size(A))<=n
    warning('Sorry bud, you can''t smooth that, pick a smaller n or use more points');
    yout = A;
    return;
end



% isVertical = size(y,1)>size(y,2);

% y=y(:)';

% padding on the left
padLeft = repmat(mean(A(:,1:n),2),1,n);
padRite = repmat(mean(A(:,end-n:end),2),1,n);
A = [padLeft A padRite];

padTop = repmat(mean(A(1:n,:),1),n,1);
padBot = repmat(mean(A(end-n:end,:),1),n,1);

A = [padTop; A; padBot];

% the main event
As=conv2(A,toConvolve);
% % % % % % % % % % % % % 
% clf;

% sanePColor(As);

% outputting a centered subset
isEven = double(~logical(round(n/2)-n/2));
if isEven
    yout = As(n+(n/2):end-n-(n/2),n+(n/2):end-n-(n/2));
else % it is odd then
    yout = As(n+(n/2)+.5:end-n-(n/2)+.5,n+(n/2)+.5:end-n-(n/2)+.5);
end

end
