Commit b96b5960 authored by Riccardo La Placa's avatar Riccardo La Placa
Browse files

LE RICERCHE SEMICOERENTI FUNZIONANO

parent 3a255b52

scsearch.m

+50 −23
Original line number Diff line number Diff line
@@ -13,9 +13,9 @@ t_raw = fitsread([pathfi,finame],"binarytable"); 
% t_raw = fitsread('C:\Users\Filippo\Desktop\EPN_0744840201_bary.fits','binarytable');

t_raw = t_raw{1};

t_raw=t_raw./86400+50814;

MJDREF=t_raw(1);

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)+ 10*256 + 1));

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)
@@ -48,8 +48,8 @@ tasc_tru = (57231.437581-MJDREF)*86400; %in secondi 
% tasc_gr = zeros(5,1);



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

porb_gr = [porb_tru-10.0,porb_tru-5.0,porb_tru,porb_tru+5.0,porb_tru+10.0].';

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

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].';

% for j = 1:5

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


toc

tic

Tseg=128; %segments' length in seconds

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)
@@ -216,12 +216,12 @@ tic 
%end



% 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

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

%s_s=uint8(4);

s_s = 4;

while(1)

    if((nismax(s_s)-nismin(s_s))<delta_ni(mu_s,s_s,g_jj(s_s)))

    if((nismax(s_s)-nismin(s_s))<0.5*delta_ni(mu_s,s_s,g_jj(s_s)))

        s_s=s_s-1;

    else

        break;
@@ -272,6 +272,8 @@ toc 
%Fourier transform on original time-series --------------------------------

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

%for each segment (lavoro su tm)

dtau = zeros(1,N);

taul = zeros(N+1,1);

disp(M);

for m=1:M
@@ -291,6 +293,8 @@ for m=1:M 
    % ttemp = gpuArray(tm(m,:));

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

    ttemp = tm(m,:);

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

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

    xtemp = x(:,m).';



    % toc
@@ -320,7 +324,13 @@ for m=1:M 
            % end

            % tic



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

            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



            % toc



            % tic
@@ -337,7 +347,8 @@ for m=1:M 
            %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(ttemp(:),xtemp(:),tau,'linear',0)).';

            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)).';
@@ -345,12 +356,19 @@ for m=1:M 
            % Y1=Y1(1:round(length(Y1)/2));

            % toc

            % tic

            cond = F1>=f_min & F1<=f_max;

            F1=F1(cond);

            Y1=Y1(cond);

            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

            Y1=Y1(cond); clear F1



            %Calcolo della detection statistic

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

            % 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?)

            

            % toc

        end
@@ -376,13 +394,18 @@ parbank = combinations(porb_gr,a_gr,tasc_gr).Variables; 
bestpar = zeros(1,5);

toc



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

% for i = 1:length(parbank)

%     for n=1:length(f_gr)

%         totlam(i,n) = totlam(i,n)+sum(Lambda(i,n,:));

%     end

% end

tic

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

        totlam = 0;

        for m = 1:M

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

            curni=zeros(1,s_s);
@@ -391,24 +414,28 @@ for n=1:length(f_gr) 
            end

            curni = -curni.*curpar(1).*curpar(3);

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

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



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

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

            %        nibank(Idx,:)

            %        curpar(1)

        end

        if (totlam > bestpar(1))

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

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

        end

        % % if (totlam > bestpar(1))

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

        % % end

    end

end

toc



disp(bestpar);

% disp(bestpar);

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

save('risultelli.mat');



for n = 1:length(f_gr)

    hold on

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

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

end

legend



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

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