Un sacco di tempo buttato probabilmente

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

%NOTE: t0 used to define the orbital phase γ is equal to the mid point

%of the observation span for dataset; M2015, Sec.7

%of the observation span for dataset; M2015, Sec.7 FAAAAAAAAAAAAAAAAALSE

% FAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAALSE

t0 = (t_raw(end)+t_raw(1))/2;

% Prendi da file il tspan

tic

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

Nyq=2000; %cambiare in caso

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(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

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)

% Time matrix with bin midpoints for each segment tm(m,j) ---------------

% tmid(m) is the midpoint in time for the m-th segment

% dt = Tseg/N; %

dt = dt_psd;

N = fix(Tseg/dt);

tm = zeros(M,N);

tmid = zeros(M,1);

for m = 1:M

    tm(m,1) = t(1)+(m-1)*Tseg; %t(j) è già centrato 

    for j = 2:N

        tm(m,j)= tm(m,1)+(j-1)*dt; 

    end

    tmid(m) = (tm(m,N)+tm(m,1))/2;

end

%CONTROLLARE SE CON LA FUNZIONE RESHAPE PUOI EVITARE IL CICLO FOR

% Time series for each segment

x = zeros(M,N);

toc

tic

aux1 = M*N;

aux2 = C(1:aux1); clear aux1

x = reshape(aux2,[M,N]); clear aux2

toc

% 

% 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];

% Aux par

Omega_max = 2*pi/porb_min;

Omega_min = 2*pi/porb_max;

gamma_min = Omega_min*(t0-tasc_max);

gamma_max = Omega_max*(t0-tasc_min);

tdiff = tmid - [tasc_max,tasc_min];

gam = tdiff(:).*[Omega_min,Omega_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),[],'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;

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

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

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;

    end

end

toc

tic

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

Nyq=2000; %cambiare in caso

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(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

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)

% Time matrix with bin midpoints for each segment tm(m,j) ---------------

% tmid(m) is the midpoint in time for the m-th segment

% dt = Tseg/N; %

dt = dt_psd;

N = fix(Tseg/dt);

tm = zeros(M,N);

tmid = zeros(M,1);

for m = 1:M

    tm(m,1) = t(1)+(m-1)*Tseg; %t(j) è già centrato 

    for j = 2:N

        tm(m,j)= tm(m,1)+(j-1)*dt; 

% Adding special case for s = 1

s = 1;

if a_max*Omega_max>1 %%%%%%%%%%%%%%%%%%%%%%%%%% FINIREEEEEEEEEEEEEEE %%%

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

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

for s = 2: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.

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

    tmid(m) = (tm(m,N)+tm(m,1))/2;

end

%CONTROLLARE SE CON LA FUNZIONE RESHAPE PUOI EVITARE IL CICLO FOR

% Time series for each segment

x = zeros(M,N);

toc

tic

aux1 = M*N;

aux2 = C(1:aux1); clear aux1

x = reshape(aux2,[M,N]); clear aux2

toc

tic

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

%for m = 1:M