Commit 67f5a6a3 authored by riccardo's avatar riccardo
Browse files

Greatest strides ever

parent e5e8a5fc

scsearch.m

+179 −181
Original line number Diff line number Diff line
@@ -2,20 +2,23 @@ clear; 
close all;



% Define variables, M is segm. number

pathfi = './';



pathfi = '/data/Sorgenti/3FGLJ1544.6-1125/';

reset(gpuDevice(2));

gpuDevice(2);

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

tic



finame = 'SiFAP2_20220429_3FGLJ1544_N1_Bary_N3_NOTCORRECTED.fits';

% finame = 'J1023_B_2017_Bary.fits';

finame = 'EPN_0744840201_bary.fits';

% finame = 'EPN_0744840201_bary.fits';

t_raw = fitsread([pathfi,finame],"binarytable");

% t_raw = fitsread('C:\Users\Filippo\Desktop\EPN_0744840201_bary.fits','binarytable');

%t_raw = fitsread('C:\Users\Filippo\Desktop\J1544_2022\SiFAP2_20220429_3FGLJ1544_N1_Bary.fits','binarytable');

t_raw = t_raw{1};

t_raw=t_raw./86400+50814;

MJDREF=fix(t_raw(1));

t_raw=(t_raw-MJDREF).*86400;

t_raw=t_raw(t_raw>=t_raw(1) & t_raw<(t_raw(1)+ 6*256 + 1));

MJDREF=0.943055555602768;

% t_raw=t_raw./86400+50814;

% MJDREF=fix(t_raw(1));

% t_raw=(t_raw-MJDREF).*86400;

% t_raw=t_raw(t_raw>=t_raw(1) & t_raw<(t_raw(1)+ 6*256 + 1));



% info = fitsinfo(finame).BinaryTable.Keywords;

% for i = 1:length(info)
@@ -33,24 +36,29 @@ toc 
% f_tru = 592.42146827248556; %Hz

% porb_tru = 17115.5216592; %s

% a_tru = 0.343356; %lt-s

%tasc_tru = 58107.009477; %MJD?

% tasc_tru = (.009477-.2326620370367891)*86400; %MJD?



% Parameters for EPN_0744840201_bary.fits

f_tru = 598.8921309; %Hz

porb_tru = 8844.08; %s

a_tru = 0.0649905; %lt-s

% tasc_tru = 57231.437581; %MJD?

tasc_tru = (57231.437581-MJDREF)*86400; %in secondi

% f_tru = 598.8921309; %Hz

% porb_tru = 8844.08; %s

% a_tru = 0.0649905; %lt-s

% % tasc_tru = 57231.437581; %MJD?

% tasc_tru = (57231.437581-MJDREF)*86400; %in secondi



% f_gr = zeros(5,1);

% porb_gr = zeros(5,1);

% a_gr = zeros(5,1);

% tasc_gr = zeros(5,1);



f_gr=f_tru+(-2:2).';

porb_gr = [porb_tru-100.0,porb_tru-50.0,porb_tru,porb_tru+50.0,porb_tru+100.0].';

a_gr = [0.01, 0.03, a_tru, 0.12, 0.15].';

tasc_gr = tasc_tru+[-500.0,-250.0,0.0,250.0,500.0].';

Tseg=490; %segments' length in seconds



%f_gr=f_tru+(-2:2).';

f_gr = (768:1/Tseg:780).';

%f_gr=(772:1/Tseg:774).';

% porb_gr = (18000:10:21600).';

porb_gr=((20868.72-3*0.31):0.1:(20868.72+3*0.31));

a_gr = (0.01:8.0e-4:0.26).';

tasc_gr = (((0.724:8e-4:0.922)-MJDREF).*86400).';

% for j = 1:5

%     %f_gr(j) = f_tru*((j^2)/9);

%     porb_gr(j) = porb_tru*((j^2)/9);
@@ -91,7 +99,7 @@ t=(t(1:end-1)+(t(2)-t(1))/2).'; % vettore tempi rebinnato, prendo il centro del 


toc

tic

Tseg=256; %segments' length in seconds



M=fix((t(end)-t(1))/Tseg); %number of segments

%con fix prendo la parte intera, scarto l'ultimo segmento che tanto non

%sarà mai di lunghezza Tseg (molto improbabile)
@@ -142,61 +150,27 @@ Omega_min = 2*pi/porb_max; 
% singah = zeros(4,1);

nismin = zeros(4,1);

nismax = zeros(4,1);

% for s = 1:4

%     singah(s) = max(sin(gam + s*pi/2),[],'all');

%     singal(s) = min(sin(gam + s*pi/2),[],'all');

% end

% 

% %% Tutto sto macello e poi ora mi viene da pensare che valori intermedi di

% %% tasc possono portare a sin(gam) anche più alti o più bassi...

% %% Guarda come ora avrei fatto molto prima a mettere semplicemente +1 e -1

% % Adding special case for s = 1

% s = 1;

% if singal(s)>0

%     nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s) + f_max;

%     nismax(s) = -f_min*a_min*(Omega_min^s)*singal(s) + f_min;

% elseif singah(s)>0

%     nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s) + f_max;

%     nismax(s) = -f_max*a_max*(Omega_max^s)*singal(s) + f_max;

% else 

%     nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s) + f_max;

%     nismax(s) = -f_min*a_min*(Omega_min^s)*singal(s) + f_min;

% end

% for s = 2:4

%     % This range is computed by finding the maximum span of Equation (15) after varying the search

%     % parameters over their respective ranges (given in Table 2). This

%     % is done with the exception of ν which is held fixed at its

%     % maximum value within sub-bands over the frequency search space.

%     %% RECHECK EVERY COMBINATION

%     if singal(s)>0

%         nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s);

%         nismax(s) = -f_min*a_min*(Omega_min^s)*singal(s);

%     elseif singah(s)>0

%         nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s);

%         nismax(s) = -f_max*a_max*(Omega_max^s)*singal(s);

%     else

%         nismin(s) = -f_max*a_max*(Omega_max^s)*singah(s);

%         nismax(s) = -f_min*a_min*(Omega_min^s)*singal(s);

%     end

% end



% Find the range in ν^s (i.e. the maximum span of Eq. (15) covered by 

% combining the search parameters over their respective ranges)

% Sin(gamma) going from -1 to 1

% Adding special case for s = 1

% s = 1;

% if a_max*Omega_max>1 

%     nismin(s) = -f_max*a_max*(Omega_max) + f_max;

%     if a_min*Omega_min>1

%         nismax(s) = -f_min*a_min*(Omega_min) + f_min;

%     else

%         nismax(s) = -f_max*a_min*(Omega_min) + f_max;

%     end

% else

%     nismin(s) = -f_min*a_max*(Omega_max) + f_min;

%     nismax(s) = -f_max*a_min*(Omega_min) + f_max;

% end

% Is (1-Ome*a) more than fmin-fmax? to ask myself

s = 1;

if a_max*Omega_max>1 

    nismin(s) = -f_max*a_max*(Omega_max) + f_max;

    if a_min*Omega_min>1

        nismax(s) = -f_min*a_min*(Omega_min) + f_min;

    else

        nismax(s) = -f_max*a_min*(Omega_min) + f_max;

    end

else

    nismin(s) = -f_min*a_max*(Omega_max) + f_min;

    nismax(s) = -f_max*a_min*(Omega_min) + f_max;

end

% The other s's are easier

nismin(s) = f_max*(1 - a_max*Omega_max);

nismax(s) = f_max*(1 + a_max*Omega_max);

for s = 2:4

    nismin(s) = -f_max*a_max*(Omega_max^s);

    nismax(s) =  f_max*a_max*(Omega_max^s);
@@ -217,7 +191,7 @@ tic 


% Try s* and check \nu_s range

g_jj=(((pi*Tseg)^2)/3).*[1; (Tseg^2)/60; (Tseg^4)/1344; (Tseg^6)/172800]; %eq. 22 M2015 + calcoli da eq. 21 M2015

mu_s=0.005; %massimo mismatch sulla griglia coerente da scegliere

mu_s=0.05; %massimo mismatch sulla griglia coerente da scegliere

%s_s=uint8(4);

s_s = 4;

while(1)
@@ -245,7 +219,7 @@ for s = 1:s_s 
        Nni(s) = Nni(s)+1;

    end

    franco=nismin(s):((nismax(s)-nismin(s))/(Nni(s))):nismax(s);

    nis{s}=franco;

    nis{s}=franco/f_max;

end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -263,11 +237,17 @@ tic 
% Cercheremo un modo per combinare in tutti i possibili modi vettori di 

% vari ni1,ni2,...,nis_s

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

% tm=gpuArray(tm);

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

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

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

Lambda = zeros(length(nibank),length(f_gr));

tm=gpuArray(tm);

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

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

F1=gpuArray((0:N-1)./(N*dt_psd)).';

f_ind = zeros(1,length(f_gr),'uint32');

for n = 1:length(f_gr)

    [bles,f_ind(n)] = min(abs(F1-f_gr(n)));

end

clear bles

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

% Lambda = zeros(length(nibank),length(f_gr));

toc 

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

%per mantenere l'informazione di fase, non faccio il valore assoluto al quadrato della fft
@@ -296,7 +276,9 @@ for m=1:M 
    ttemp = tm(m,:);

    ttempdif = (ttemp(:)-tmid(m)).';

    ttempdifl = (ttemp(:)-tmid(m)-0.5*dt).';

    xtemp = x(:,m).';

    xtemp = gpuArray(x(:,m).');

    sumx = sum(xtemp);

    



    % toc
@@ -306,95 +288,45 @@ for m=1:M 
    tic

    for i = 1:length(nibank)



        % Tento ricampionamento (Eq. 17 MP2015)

        % 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 = f_gr(n);

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

            % tau = zeros(N,1);

            % tic

            % for j = 1:N

            %     % ta(j) = ; No Tau_zero: indifferente per il calcolo della FFT

            %     for s = 1:s_s

            %         tau(j) = tau(j) + (nibank(i,s)/(nizero*factorial(s)))*(tm(m,j)-tmid(m))^s;

            %     end

            % end

            % tic



            tau = tmid(m) + sum((nibank(i,1:s_s)./(nizero*factorial(1:s_s))).*((ttempdif(1:N)).^((1:s_s).').'),2);

            % taul(1:N) = tmid(m) + sum((nibank(i,1:s_s)./(nizero*factorial(1:s_s))).*((ttempdifl(1:N)).^((1:s_s).').'),2);

            % taul(N+1) = tmid(m) + sum((nibank(i,1:s_s)./(nizero*factorial(1:s_s))).*((ttemp(N)-tmid(m)+0.5*dt).^((1:s_s).').'));

            % for j=1:N

            %     dtau(j) = ((taul(j+1)-taul(j))); 

            % end

            taul(1:2) = tmid(m) + sum((nibank(i,1:s_s)./(nizero*factorial(1:s_s))).*((ttempdifl(1:2)).^((1:s_s).').'),2);

        tau = tmid(m) + sum((nibank(i,1:s_s)./(factorial(1:s_s))).*((ttempdif(1:N)).^((1:s_s).').'),2);

        taul(1:2) = tmid(m) + sum((nibank(i,1:s_s)./(factorial(1:s_s))).*((ttempdifl(1:2)).^((1:s_s).').'),2);

        dtau = ((taul(2)-taul(1)));





            % toc



            % tic

            % X1 = interp1(ttemp(:),xtemp(:),tau,'linear',0);

            % toc



            %X1 è la timeseries ricampionata (controllare che sia un vettore colonna)

            %zero-padding (metto gli zeri alla fine) -------------------------------

            % X1 = [X1;zeros(L-length(F),1)]; %concateno alla timeseries downsampled and resampled gli zeri

            %il secondo è un vettore colonna di zeri con lunghezza pari al numero di elementi che mi serve per tornare alla risoluzione originale



            %Fourier transform -----------------------------------------------------

            %[Cm,edges]=(histcounts(X1,round((X1(end)-X1(1))/dt_psd)));

            %edges=edges(end)-edges(2); %mi dà il tempo preciso di tutta la TdF, che sarà leggermente diversa da length(C)*dt per come è definito histcounts

            % Y1=(2./sum(X1).*abs(fft(X1)).^2).'; %normalizzazione Leahy, giusto?????

            % tic

            % Y1 = fft(interp1(ttemp(:),xtemp(:),tau,'linear',0)).';

            % Y1 = fft((interp1(tau(:),(xtemp(:)./dtau(:)),ttemp(:),'linear',0)).*dt).';

        % inzo = (interp1(tau(:),(xtemp(:)./dtau),ttemp(:),'linear',0)).*dt;

        Y1 = fft((interp1(tau(:),(xtemp(:)./dtau),ttemp(:),'linear',0)).*dt).';

            % clear X1

            %F1=((0:length(Y1)-1)./(ttemp(end)-ttemp(1))).'; 

            F1=((0:length(Y1)-1)./(length(Y1)*dt_psd)).';

            % F1=F1(1:round(length(F1)/2));

            % Y1=Y1(1:round(length(Y1)/2));

            % toc

            % tic



            % if n ~= 1 && n ~= length(f_gr)

            %     cond = F1>=0.5*(nizero+f_gr(n-1)) & F1<0.5*(nizero+f_gr(n+1));

            % elseif n == 1

            %     cond = F1>=(nizero) & F1<0.5*(nizero+f_gr(n+1));

            % else

            %     cond = F1>=0.5*(nizero+f_gr(n-1)) & F1<=nizero;

            % end

            [nnfreq,nnind] = min(abs(F1-nizero));

            % cond = F1>=F1(nnind-1) & F1<=F1(nnind+1);

            Y1=Y1(nnind); clear F1

        

            %Calcolo della detection statistic

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

            % [nnfreq,nnind] = min(abs(F1-curpar(1)));

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

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

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

        % currently being considered)

        for n = 1:length(f_gr)



            Yenne=Y1(f_ind(n)); 

            Lambda(i,n) = 2*(abs(Yenne)^2)/sumx; %CREDO (oppure prendono la potenza massima?)

            

            % toc

        end

        % disp('1ni')

    end

    toc

    

    disp(length(f_gr)*i);

    disp(m);

    lamfiname = sprintf('Lambda_seg_%d.mat', m);

    save(lamfiname,"Lambda");

    % save(lamfiname,"Lambda");

    save([pathfi,lamfiname],"Lambda","-v7.3");

    Lambda = zeros(length(nibank),length(f_gr));

    disp(m);



end





%% 





tic

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

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

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

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

%% The line above will have to be substituted with something along the lines

%% of the ones below to account for the metric g_jj in the phase derivatives

% gdistance = @(a,b)sqrt(((a-b).^2)*(g_jj(1:s_s)));
@@ -408,46 +340,112 @@ toc 


totlam = zeros(length(parbank),length(f_gr));

tic

for n=1:length(f_gr)

% for n=1:length(f_gr)

%     for i=1:length(parbank)

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

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

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

%         for m = 1:M

%             curpar(4) = curpar(2)*(tmid(m) - parbank(i,3));

%             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(3);

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

% 

%             lamfiname = sprintf('Lambda_seg_%d.mat', m);

%             % load(lamfiname);

%             load([pathfi,lamfiname]);

% 

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

%             totlam(i,n) = totlam(i,n)+Lambda(Idx,n);

% 

%             clear Lambda

% 

%         end

%         % % if (totlam > bestpar(1))

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

%         % % end

%     end

% end

bles = cospi(1/2);

for m = 1:M

    tic

    lamfiname = sprintf('Lambda_seg_%d.mat', m);

    % load(lamfiname);

    load([pathfi,lamfiname]);

    toc

    disp(m);

    tic

    for i=1:length(parbank)

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

        

        curpar = [bles,parbank(i,:)];

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

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

        for m = 1:M

        curpar(4) = curpar(2)*(tmid(m) - parbank(i,3));

        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(3);

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



            lamfiname = sprintf('Lambda_seg_%d.mat', m);

            load(lamfiname);

            

        curni = -curni.*curpar(3);

        curni(1) = curni(1) + 1;

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

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



            clear Lambda

        totlam(i,:) = totlam(i,:)+Lambda(Idx,:);



    end

        % % if (totlam > bestpar(1))

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

        % % end

    end

    toc 

    disp(m);

    clear Lambda

end

% % for n=1:length(f_gr)

% %     for i=1:length(parbank)

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

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

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

% %         for m = 1:M

% %             curpar(4) = curpar(2)*(tmid(m) - parbank(i,3));

% %             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(3);

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

% % 

% %             lamfiname = sprintf('Lambda_seg_%d.mat', m);

% %             % load(lamfiname);

% %             load([pathfi,lamfiname]);

% % 

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

% %             totlam(i,n) = totlam(i,n)+Lambda(Idx,n);

% % 

% %             clear Lambda

% % 

% %         end

% %         % % if (totlam > bestpar(1))

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

% %         % % end

% %     end

% % end

toc

%% 



% disp(bestpar);

% save('C:\Users\Filippo\Desktop\XMM_Jxxx\risultelli.mat');

save('risultelli.mat');

save([pathfi,'risultelli_07_03_',num2str(Tseg),'s.mat'],"-v7.3");

%% 



for n = 1:length(f_gr)

    hold on

    txt = ['f_{gr} = ',num2str(f_gr(n))];

    plot(totlam(:,n),'DisplayName',txt)

end

legend



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

%%%%%%%%%%%%% Fai il grafico delle totlam a posteriori %%%%%%%%%%%%%%%%



% for n = 1:length(f_gr)

%     hold on

%     txt = ['f_{gr} = ',num2str(f_gr(n))];

%     plot(totlam(:,n),'DisplayName',txt)

% end

% legend



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

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

test_efficiency.m

+111 −57
Original line number Diff line number Diff line
@@ -155,7 +155,7 @@ nismax = zeros(4,1); 
% combining the search parameters over their respective ranges)

% Sin(gamma) going from -1 to 1

% Adding special case for s = 1

s = 1;

% s = 1;

% if a_max*Omega_max>1 

%     nismin(s) = -f_max*a_max*(Omega_max) + f_max;

%     if a_min*Omega_min>1
@@ -168,8 +168,9 @@ s = 1; 
%     nismax(s) = -f_max*a_min*(Omega_min) + f_max;

% end

% Is (1-Ome*a) more than fmin-fmax? to ask myself

nismin(s) = -f_min*a_max*(Omega_max) + f_min;

nismax(s) = -f_max*a_min*(Omega_min) + f_max;

s = 1;

nismin(s) = f_max*(1 - a_max*Omega_max);

nismax(s) = f_max*(1 + a_max*Omega_max);

for s = 2:4

    nismin(s) = -f_max*a_max*(Omega_max^s);

    nismax(s) =  f_max*a_max*(Omega_max^s);
@@ -218,7 +219,7 @@ for s = 1:s_s 
        Nni(s) = Nni(s)+1;

    end

    franco=nismin(s):((nismax(s)-nismin(s))/(Nni(s))):nismax(s);

    nis{s}=franco;

    nis{s}=franco/f_max;

end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -242,9 +243,9 @@ Lambda = zeros(length(nibank),length(f_gr),'gpuArray'); 
F1=gpuArray((0:N-1)./(N*dt_psd)).';

f_ind = zeros(1,length(f_gr),'uint32');

for n = 1:length(f_gr)

    [ble,f_ind(n)] = min(abs(F1-f_gr(n)));

    [bles,f_ind(n)] = min(abs(F1-f_gr(n)));

end

clear ble

clear bles

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

% Lambda = zeros(length(nibank),length(f_gr));

toc
@@ -285,39 +286,27 @@ for m=1:M 
    % Dev'essere fatto per ciascun template, ergo mettiamo un bel ciclo 

    % sulle possibili combinazioni di ni_s

    tic

    % for i = 1:length(nibank)

    for i = 1:1

    for i = 1:length(nibank)

    % for i = 1:1



        % Tento ricampionamento (Eq. 17 MP2015)

        % 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 = f_gr(n);

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

            % tau = zeros(N,1);

            % tic

            % for j = 1:N

            %     % ta(j) = ; No Tau_zero: indifferente per il calcolo della FFT

            %     for s = 1:s_s

            %         tau(j) = tau(j) + (nibank(i,s)/(nizero*factorial(s)))*(tm(m,j)-tmid(m))^s;

            %     end

            % end

            % tic



            tau = tmid(m) + sum((nibank(i,1:s_s)./(nizero*factorial(1:s_s))).*((ttempdif(1:N)).^((1:s_s).').'),2);

            taul(1:2) = tmid(m) + sum((nibank(i,1:s_s)./(nizero*factorial(1:s_s))).*((ttempdifl(1:2)).^((1:s_s).').'),2);

        tau = tmid(m) + sum((nibank(i,1:s_s)./(factorial(1:s_s))).*((ttempdif(1:N)).^((1:s_s).').'),2);

        taul(1:2) = tmid(m) + sum((nibank(i,1:s_s)./(factorial(1:s_s))).*((ttempdifl(1:2)).^((1:s_s).').'),2);

        dtau = ((taul(2)-taul(1)));



        % inzo = (interp1(tau(:),(xtemp(:)./dtau),ttemp(:),'linear',0)).*dt;

        Y1 = fft((interp1(tau(:),(xtemp(:)./dtau),ttemp(:),'linear',0)).*dt).';

            Y1=Y1(f_ind(n)); 



            Lambda(i,n) = 2*(abs(Y1)^2)/sumx; %CREDO (oppure prendono la potenza massima?)

        

            % toc

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

        % currently being considered)

        for n = 1:length(f_gr)



            Yenne=Y1(f_ind(n)); 

            Lambda(i,n) = 2*(abs(Yenne)^2)/sumx; %CREDO (oppure prendono la potenza massima?)

            

        end

        % disp('1ni')

    end
@@ -332,6 +321,10 @@ for m=1:M 


end





%% 





tic

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

nisearcher = KDTreeSearcher(gather(nibank),'BucketSize',100);
@@ -348,40 +341,101 @@ toc 


totlam = zeros(length(parbank),length(f_gr));

tic

for n=1:length(f_gr)

% for n=1:length(f_gr)

%     for i=1:length(parbank)

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

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

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

%         for m = 1:M

%             curpar(4) = curpar(2)*(tmid(m) - parbank(i,3));

%             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(3);

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

% 

%             lamfiname = sprintf('Lambda_seg_%d.mat', m);

%             % load(lamfiname);

%             load([pathfi,lamfiname]);

% 

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

%             totlam(i,n) = totlam(i,n)+Lambda(Idx,n);

% 

%             clear Lambda

% 

%         end

%         % % if (totlam > bestpar(1))

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

%         % % end

%     end

% end

bles = cospi(1/2);

for m = 1:M

    tic

    lamfiname = sprintf('Lambda_seg_%d.mat', m);

    % load(lamfiname);

    load([pathfi,lamfiname]);

    toc

    disp(m);

    tic

    for i=1:length(parbank)

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

        

        curpar = [bles,parbank(i,:)];

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

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

        for m = 1:M

        curpar(4) = curpar(2)*(tmid(m) - parbank(i,3));

        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(3);

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



            lamfiname = sprintf('Lambda_seg_%d.mat', m);

            % load(lamfiname);

            load([pathfi,lamfiname]);



        curni = -curni.*curpar(3);

        curni(1) = curni(1) + 1;

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

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



            clear Lambda

        totlam(i,:) = totlam(i,:)+Lambda(Idx,:);



    end

        % % if (totlam > bestpar(1))

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

        % % end

    end

    toc 

    disp(m);

    clear Lambda

end

% % for n=1:length(f_gr)

% %     for i=1:length(parbank)

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

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

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

% %         for m = 1:M

% %             curpar(4) = curpar(2)*(tmid(m) - parbank(i,3));

% %             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(3);

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

% % 

% %             lamfiname = sprintf('Lambda_seg_%d.mat', m);

% %             % load(lamfiname);

% %             load([pathfi,lamfiname]);

% % 

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

% %             totlam(i,n) = totlam(i,n)+Lambda(Idx,n);

% % 

% %             clear Lambda

% % 

% %         end

% %         % % if (totlam > bestpar(1))

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

% %         % % end

% %     end

% % end

toc

%% 



% disp(bestpar);

% save('C:\Users\Filippo\Desktop\XMM_Jxxx\risultelli.mat');

save([pathfi,'risultelli_',num2str(Tseg),'s.mat'],"-v7.3");

save([pathfi,'risultelli_07_03_',num2str(Tseg),'s.mat'],"-v7.3");

%% 





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