ni_s only contains derivatives now, not the frequency (this shrinks nibank).

%%  Save this template bank, we'll need it later

% Nni+1 is the number of ni_s points in the grid (counting borders, 

% Nni is the number of steps) for each s couple.

% We make sure to oversample with respect to eq. 23 M2015 to ensure the max

% mismatch is not exceeded. -------------

% We make sure to oversample with respect to eq. 23 M2015 to ensure the 

% max mismatch is not exceeded. -------------

% nis is a structure which holds the various vectors of ni_s^(m) (our

% coherent template bank)

% each nis{s} is the array giving the respective grid, and its i-th 

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nis{s_s+1}=f_gr.';

% Changing nis to stop at s*(= s_s): this way nibank is length(f_gr)

% times smaller, which also helps in the NNsearch

%%%%%%%%%%%% -> commenting the next line and adding a cycle on f_gr

% nis{s_s+1}=f_gr.';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

toc

%%%[kung,fu,fight] = ndgrid(nis{1},nino{2},nino{3});

tm=gpuArray(tm);

nibank = gpuArray(combinations(nis{:}).Variables);

Lambda = zeros(length(nibank),M,'gpuArray');

Lambda = zeros(length(nibank),length(f_gr),M,'gpuArray');

toc 

%Fourier transform on original time-series --------------------------------

%per mantenere l'informazione di fase, non faccio il valore assoluto al quadrato della fft

    tic

    for i = 1:length(nibank)

        % Tento ricampionamento (Eq. 17 MP2015)

        % Per prima cosa, generiamo la nuova coordinata temporale per 

        % questa templ combini 

        % Careful! ni_0 still needs to be defined (it's the spin freq. 

        % currently being considered)

        for n = 1:length(f_gr)

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        nizero = nibank(i,s_s+1);

            nizero = f_gr(n);

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            % tau = zeros(N,1);

            % tic

            Y1=Y1(cond);

            %Calcolo della detection statistic

        %QUESTA SARA' DA INIZIALIZZARE ALL'ESTERNO tipo Lambda = zeros(length(template),1);

        Lambda(i,m)=sum(abs(Y1).^2)/sum(x(m,:)); %CREDO (oppure prendono la potenza massima?)

            Lambda(i,n,m)=sum(abs(Y1).^2)/sum(x(m,:)); %CREDO (oppure prendono la potenza massima?)

            toc

        end

        disp('1ni')

    end

    toc

nisearcher = KDTreeSearcher(nibank,'BucketSize',100);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

parbank = combinations(f_gr,porb_gr,a_gr,tasc_gr).Variables;

bestpar = zeros(5,1);

% parbank = combinations(f_gr,porb_gr,a_gr,tasc_gr).Variables;

parbank = combinations(porb_gr,a_gr,tasc_gr).Variables;

bestpar = zeros(1,5);

toc

%% Find out how to search for nearest neighbour while keeping one par. fixed

tic

for n=1:length(f_gr)

    for i=1:length(parbank)

    curpar = parbank(i,:);

        curpar = [f_gr(n),parbank(i,:)];

        % Move from ν,P,a,T_asc to ν,Ω,a,γ

        curpar(2) = 2*pi/curpar(2);

        totlam = 0;

        for m = 1:M

            curpar(4) = curpar(2)*(tmid(m) - curpar(4));

        curni=zeros(1,s_s+1);

            curni=zeros(1,s_s);

            for s=1:s_s

                curni(s) = (curpar(2)^s)*sin(curpar(4)-0.5*s*pi);

            end

            curni=curni.*curpar(1).*curpar(2);

        curni(s_s+1) = curpar(1);

            % curni(s_s+1) = curpar(1);

            [Idx,D] = knnsearch(nisearcher,curni);

        totlam = totlam + (Lambda(Idx,m));

            totlam = totlam + (Lambda(Idx,n,m));

            %        nibank(Idx,:)

            %        curpar(1)

        end

        if (totlam > bestpar(1))

        bestpar = [totlam;parbank(i,:).'];

            bestpar = [totlam,f_gr(n),parbank(i,:)];

        end

    end

end

toc