function H = histw_phase(sample, weight, action, hlist)
% Histogram for trajectories including phases
% INPUT: sample - what to hist; weight - which weight for sample; action -
% phase input (will be in exp(i*action/hbar); hlist - output xlist
% OUTPUT: H - ylist of histogram (complex amplitude)

load_const;

H = zeros(size(hlist));

hmin = min(hlist);
dh = hlist(2) - hlist(1);

for n = 1:numel(sample)
    
    index = 1 + floor(sample(n)/dh);

    gain_r = 1 + (sample(n)/dh) - index;
    gain_l = 1 - gain_r;

    if index < 1
        continue;
    end
    
    if index > numel(H)-1
        continue;
    end
    
    H(index) = H(index) + gain_l * sqrt(weight(n)) * exp(1i*action(n)/const.hbar) / dh;
    H(index+1) = H(index+1) + gain_r * sqrt(weight(n)) * exp(1i*action(n)/const.hbar) / dh;
end

% disp("returns complex spectrum");

end