Commit 6854ec69 authored by Riccardo La Placa's avatar Riccardo La Placa
Browse files

Committing only to have latest updates (tiny)

parent 3afa76ce

scsearch.m

+42 −38
Original line number Diff line number Diff line
@@ -19,6 +19,7 @@ t_raw = t_raw{1}; 
t0 = (t_raw(end)+t_raw(1))/2;



% Prendi da file il tspan

% Add step to generate par. extrema

tspan = t_raw(end)-t_raw(1);

param = [f_min,f_max,porb_min,porb_max,a_min,a_max,tasc_min,tasc_max];
@@ -28,6 +29,36 @@ Omega_min = 2*pi/porb_max; 
gamma_min = Omega_min*(t0-tasc_max);

gamma_max = Omega_max*(t0-tasc_min);



% Find largest singam

% M2015 eq. 15

gam = [gamma_max,gamma_min];

singal = zeros(4,1);

singah = zeros(4,1);

nismin = zeros(4,1);

nismax = zeros(4,1);

for s = 1:4

    singah(s) = max(sin(gam - s*pi/2));

    singal(s) = min(sin(gam - s*pi/2));

    % 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.

    nismin(s) = f_min*a_min*(Omega_min^s)*singal(s);

    nismax(s) = f_max*a_max*(Omega_max^s)*singah(s);

end



% Try s* and check \nu_s range

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

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

s_s=int16(4);

while(1)

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

        s_s=s_s-1;

    else

        break;

    end

end



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

Nyq=2000; %cambiare in caso
@@ -35,9 +66,9 @@ dt_psd=1/(2*Nyq); %risoluzione spettro di potenza, diverso dal dt dei segmenti 
nbin=round((t_raw(end)-t_raw(1))/dt_psd);

[C,t]=(histcounts(t_raw,nbin)); 

C=C.'; %conteggi

t=(t(2:end)).'; %vettore tempi rebbinato, assegno al tempo l'edge destro di ogni bin

%se serve il centro del bin devo fare: t=t(1:end-1)+(t(2)-t(1))/2; 

%perchè dt=t(2)-t(1)

t=(t(1:end-1)+(t(2)-t(1))/2).'; % vettore tempi rebinnato, prendo il centro del bin 

% t=(t(2:end)).'; %vettore tempi rebinnato, assegno al tempo l'edge destro di ogni bin



%VEDERE SE E' NECESSARIO FARE TUTTO IL CICLO SOTTO PER I tm

%Come fatto qui sopra, tutti i conteggi sono assegnati al centro di ogni

%bin che dura dt
@@ -67,38 +98,8 @@ x = zeros(M,N); 
for m = 1:M

    for j = 1:N

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

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

    end

end





% Find largest singam

% M2015 eq. 15

gam = [gamma_max,gamma_min];

singal = zeros(4,1);

singah = zeros(4,1);

nismin = zeros(4,1);

nismax = zeros(4,1);

for s = 1:4

    singah(s) = max(sin(gam - s*pi/2));

    singal(s) = min(sin(gam - s*pi/2));

    % 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.

    nismin(s) = f_min*a_min*(Omega_min^s)*singal(s);

    nismax(s) = f_max*a_max*(Omega_max^s)*singah(s);

end



% Try s* and check \nu_s range

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

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

s_s=int16(4);

while(1)

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

        s_s=s_s-1;

    else

        break;

        Cmj=C(pippo);

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

    end

end
@@ -115,8 +116,8 @@ end 
Nni = zeros(s_s,1);

for s = 1:s_s

    Nni(s) = ceil((nismax(s)-nismin(s))/delta_ni(mu_s,s,g_jj(s)));

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

    nis{s}=a;

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

    nis{s}=franco;

end



% Resampling the times over nism{m,s} templates
@@ -132,8 +133,11 @@ nibank = combinations(nis{:}).Variables; 


%for each segment (lavoro su tm)

for m=1:M

    [Cm,edges]=(histcounts(tm(m,:),round((tm(m,end)-tm(m,1))/dt_psd))); 

    [Cm,edges]=(histcounts(x(m,:),round((tm(m,end)-tm(m,1))/dt_psd))); % R - convincitene

    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 

    %% La riga 137 non può essere giusta, si perde sicuramente un pezzo 

    %% di tm(1) che è preso come primo edge da histcounts: va sostituito 

    %% edges(2) con edges(1)

    Y=fft(Cm).'; clear Cm

    F=((0:length(Y)-1)./edges).'; clear edges

    L=length(F); %lunghezza iniziale, servirà per lo zero-padding