Commit 2499b91e authored by Riccardo La Placa's avatar Riccardo La Placa
Browse files

Quanto siamo forti! 

Abbiamo sistemato i problemi con ni_zero, con tau, con Lambda e il resto; manca da ridurre il tempo di sversamento dati e da sistemare il posizionamento dei par. fisici
parent 60556dfb

J1023_B_2017_Bary.fits

0 → 100644
+111 MiB

File added.

No diff preview for this file type.

View file

scsearch.m

+30 −11
Original line number Diff line number Diff line
@@ -4,7 +4,10 @@ close all; 
% Define variables, M is segm. number

pathfi = './';



finame = 'franco.fits';

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

tic



finame = 'J1023_B_2017_Bary.fits';

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

t_raw = t_raw{1};

% info = fitsinfo(finame).BinaryTable.Keywords;
@@ -14,6 +17,9 @@ t_raw = t_raw{1}; 
%     end

% end



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

toc



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

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

% Lines added just to test what was done until now
@@ -27,6 +33,8 @@ porb_gr = zeros(5,1); 
a_gr = zeros(5,1);

tasc_gr = zeros(5,1);



tic



for j = 1:5

    f_gr(j) = f_tru*((j^3)/9);

    porb_gr(j) = porb_tru*((j^3)/9);
@@ -79,7 +87,9 @@ for s = 1:4 
    nismax(s) = f_max*a_max*(Omega_max^s)*singah(s);

end



toc



tic

%Rebin -------------------------------------------------------------------

Nyq=2000; %cambiare in caso
@@ -94,6 +104,8 @@ t=(t(1:end-1)+(t(2)-t(1))/2).'; % vettore tempi rebinnato, prendo il centro del 
%Come fatto qui sopra, tutti i conteggi sono assegnati al centro di ogni

%bin che dura dt



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
@@ -117,15 +129,20 @@ end 


% Time series for each segment

x = zeros(M,N);

toc



tic

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

%% Lasciarlo così per circa 13E+6 bin temporali richiede ~88 ore (su questo

%% computer): trovare modo più sensato

for m = 1:M

    for j = 1:N

        pippo = t>=tm(m,j)-dt/2 & t<tm(m,j)+dt/2;

        Cmj=C(pippo);

        Cmj=C(pippo); clear pippo

        x(m,j)=sum(Cmj);

    end

end



toc

% Try s* and check \nu_s range

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

mu_s=0.001; %massimo mismatch sulla griglia coerente da scegliere
@@ -168,7 +185,7 @@ nis{s_s+1}=f_gr; 
% vari ni1,ni2,...,nis_s

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

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

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

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



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

%per mantenere l'informazione di fase, non faccio il valore assoluto al quadrato della fft
@@ -201,14 +218,14 @@ for m=1:M 
        nizero = nibank(i,s_s+1);

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

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

        ta = zeros(N,1);

        tau = zeros(N,1);

        for j = 1:N

            ta(j) = ;

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

            for s = 1:s_s

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

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

            end

        end

        X1 = interp1(tm(m,:),X,ta,'linear',0);

        X1 = interp1(tm(m,:),X,tau,'linear',0);



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

        %zero-padding (metto gli zeri alla fine) -------------------------------
@@ -218,7 +235,8 @@ for m=1:M 
        %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?????

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

        Y1 = fft(X1).';

        clear X1

        F1=((0:length(Y1)-1)./(tm(m,end)-tm(m,1))).';

        F1=F1(1:round(length(F1)/2));
@@ -226,9 +244,10 @@ for m=1:M 


        %Calcolo della detection statistic

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

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

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

    end



    

end