function [hfdf_map,epoch,dx]=LongT_GENERALIZED_fasthough(peaks,hm_job)

%LongT_GENERALIZED_fasthough  creates a x/k Hough map (based on Andrew Miller 2019 version)
%
% INPUT
%    peaks(3,n)      peaks of the peakmap as [t,fr,amp] (corrected for the Doppler effect)
%                      time is in seconds, frequency in Hz, amp not used
%
%    hm_job          hough map structure 
%        .minf       minimum frequency
%        .maxf       maximum frequency
%        .df         frequency resolution
%        .frenh      over-resolution factor (default as 1)
%        .dur        map duration, in seconds
%        .ref_perc_time  percentile of reference time, from 0 to 1
%        .n          braking index of the signal evolution power law
%        .gridk      vector containing the grid of k values to be seached
%
% OUTPUT
%    hfdf            hough map (gd2)
%    job_info        job info structure
%    checkE          service structure for test and debug
%
% Version 1.0 - February 2024
% Developed by Lorenzo Pierini within PNRR ICSC, Spoke 2 project
% Progetto ICSC - Spoke 2 - Codice CN00000013 - CUP I53C21000340006 - 
% Missione 4 Istruzione e ricerca - 
% Componente 2 Dalla ricerca all'impresa – Investimento 1.4.
%
% Part of Snag toolbox - Signal and Noise for Gravitational Antennas

pow=hm_job.n-1;

if ~isfield(hm_job,'ref_perc_time')
    ref_perc_time=0.5;
else
    ref_perc_time=hm_job.ref_perc_time;
end

[~,n_peaks]=size(peaks);
weights=ones(1,n_peaks);
weights(peaks(2,:)<0)=0;

try
    epoch=quantile(peaks(1,:),ref_perc_time); %beginning, middle, or end of observing run
catch
    epoch=min(peaks(1,:)); %if the above fails, place t0 at the beginning of the observing run
end

peaks(1,:)=(peaks(1,:)-epoch)*86400;
Tmax=max(abs(peaks(1,:)));


minf0=hm_job.minf;
maxf0=hm_job.maxf;
df=hm_job.df;
enh=hm_job.frenh;


if hm_job.n==1 % case of pulsar winds
    xx=log(peaks(2,:));
    dx=df/maxf0;
    ddx=dx/enh;
    minx0=floor(log(minf0)/ddx)*ddx;
    maxx0=ceil(log(maxf0)/ddx)*ddx;
    pow=1; % not physical
else
    xx=peaks(2,:).^-pow;
    dx=pow*df/(maxf0^hm_job.n);
    ddx=dx/enh;
    minx0=1/maxf0^pow;%floor((1/maxf0^pow)/ddx)*ddx;
    maxx0=1/minf0^pow;%ceil((1/minf0^pow)/ddx)*ddx;
    
    if maxx0<minx0
        temp1=minx0;
        temp2=maxx0;
        maxx0=temp1;
        minx0=temp2;
    end
    
    % smallest possible step
    dx=abs(dx);
end

gridk=hm_job.gridk;
if hm_job.n==11/3
%     disp('chirp, flipping spindowns to spinups')
    gridk=-gridk;
end

sd_belt=ceil(Tmax*pow*max(abs(gridk))/ddx)+10;
nbin_fbandx=ceil(abs(maxx0-minx0)/ddx);

inix=minx0-sd_belt*ddx;


I1000=0;
I500=I1000;

nbin_x=nbin_fbandx+2*sd_belt+I1000; 

% INDEXES NO MORE USED
% n_of_peaks=length(peaks);
% ii=find(diff(peaks(1,:)));  % find the different times
% ii=[ii n_of_peaks]; 
% nTimeSteps=length(ii); % number of times of the peakmap

% dk is used to normalize the gridk values (numerical errors less likely)
dk=min(abs(diff(gridk)));
% normalized
xx_norm=(xx-inix)/ddx-I500;
tt_norm=peaks(1,:)*dk*pow/ddx;
slopes=gridk/(dk);

% CONTROLS ON INDEXES IN CASE OF ERRORS
% error_n_x_bins=max(xx_norm)-min(xx_norm)-nbin_fbandx
% lower_expect_hm_x=min(xx_norm)+min([-max(tt_norm)*max(slopes),min(tt_norm)*max(slopes)])
% expected_hm_x_exceed=nbin_x-max(xx_norm)-max([max(tt_norm)*max(slopes),-min(tt_norm)*max(slopes)])

% The following line substitutes the whole hough loop
hfdf_cell=arrayfun(@(x) accumarray(1+round(xx_norm-tt_norm.*x).',weights.',[nbin_x,1]),slopes,'UniformOutput',false);

% Here I cut the sd-belts added to do the histogram
for i=1:length(slopes)
    hfdf_cell{i}=hfdf_cell{i}(sd_belt+1:sd_belt+nbin_fbandx);
end
binh_df0=cell2mat(hfdf_cell).';

% half to full diff. map - Carl Sabottke idea - STILL TO BE TESTED!!!
halfstrip=0;
if halfstrip >0
    binh_df0(:,halfstrip*2+1:nbin_x)=binh_df0(:,halfstrip*2+1:nbin_x)-binh_df0(:,1:nbin_x-halfstrip*2); 
    binh_df0=cumsum(binh_df0,2);
end


inix=minx0;
hfdf_map=gd2(binh_df0.'); 

hfdf_map=edit_gd2(hfdf_map,'dx',ddx,'ini',inix-ddx*I500,'capt','Histogram of y0-x0'); %inix-I500*ddf

end