function [h2, T, res_var]=hist_eq2(h1,N)
%h1=[10 30 60 20];
K = length(h1);
%N = 2;
h2 = zeros([1 N]);
q = sum(h1)/N;
t=1;
T = zeros([1 K]);
for n=1:N,
  if t > K, 
      break;
  end
  flag = 1;
  h2(n) = h2(n) + h1(t);
  T(t) = n;
  t = t+1;

  while t <= K & flag ==1,
    flag = 0;
    if (h2(n)+h1(t)-q).^2 <= (h2(n)-q).^2,
       h2(n) = h2(n) + h1(t);
       T(t) = n;
       t = t+1;
       flag = 1;
    end
  end
end
while t <= K,
  T(t) = N;
  h2(N) = h2(N)+h1(t);
  t= t+1;
end

disp([sprintf('Var = %12.10f', sum((h2-q).^2))]);
res_var = sum((h2-q).^2);

return;
