Post merge; primo accenno di nearest neighbour e derivate sistemate

close all;

% Define variables, M is segm. number

pathfi = './';

%pathfi = './';

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

tic

finame = 'J1023_B_2017_Bary.fits';

% finame = 'J1023_B_2017_Bary.fits';

%finame = 'EPN_0744840201_bary.fits';

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

%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(t_raw>=t_raw(1) & t_raw<t_raw(1)+900);

t_raw=t_raw./86400+50814;

MJDREF=t_raw(1);

t_raw=(t_raw-MJDREF).*86400;

t_raw=t_raw(t_raw>=t_raw(1) & t_raw<t_raw(1)+1200);

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

% for i = 1:length(info)

%     if (isequal(info{i},'TIMEDEL'))

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

% Lines added just to test what was done until now

% Parameters for J1023_B_2017_Bary.fits

f_tru = 592.42146827248556; %Hz

porb_tru = 17115.5216592; %s

a_tru = 0.343356; %lt-s

tasc_tru = 58107.009477; %MJD?

%f_tru = 592.42146827248556; %Hz

%porb_tru = 17115.5216592; %s

%a_tru = 0.343356; %lt-s

%tasc_tru = 58107.009477; %MJD?

% Parameters for EPN_0744840201_bary.fits

% f_tru = 598.8921309; %Hz

% porb_tru = 8844.08; %s

% a_tru = 0.0649905; %lt-s

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

Omega_max = 2*pi/porb_min;

Omega_min = 2*pi/porb_max;

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

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;

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

% 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

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

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

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

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

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

tm=gpuArray(tm);

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

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

toc 

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

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

tic

%for each segment (lavoro su tm)

for m=1:M

end

tic

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

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

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

parbank = combinations(f_gr,porb_gr,a_gr,tasc_gr).Variables;

bestpar = zeros(5,1);

toc

%% Find out how to search for nearest neighbour while keeping one par. fixed

tic

for i=1:length(parbank)

    curpar = parbank(i,:);

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

        bestpar = [totlam;parbank(i,:).'];

    end

end

toc

disp(bestpar);

save('risultelli.mat');