Another wide variety of changes, first among them the photon-by-photon correction

clear;

close all;

% Define variables, M is segm. number

import matlab.io.*

pathfi = '/data/Sorgenti/3FGLJ1544.6-1125/';

% pathfi = '/data/Sorgenti/SCORPIUS-X-1/';

% pathfi = '/data/Sorgenti/J1023/';

folname = [pathfi,'single_photon_test_',char(datetime('now','Format','dd_MM'))];

if (exist(folname,"dir")==7)

    disp(['Folder ',folname,' already exists, previous Lambda files (if present) may be lost']);

else

    mkdir(folname);

    disp(['Folder ',folname,' created']);

end

pathfi = [folname,'/'];

reset(gpuDevice(2));

gpuDevice(2);

diary([pathfi,'loggerino',char(datetime('now','Format','dd_MM')),'.log'])

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

daterella = ['Date and time now ',char(datetime('now')),' (UTC), maybe'];

disp(daterella); clear daterella

tic

finame = 'SiFAP2_20220429_3FGLJ1544_N1_Bary.fits';

% finame = 'ScoX1_20230423_Bary.fits';

% finame = 'J1023_B_2017_Bary.fits';

% finame = 'EPN_0744840201_bary.fits';

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

filoc = [pathfi,'../',finame];

file=fits.openFile(filoc);

mjdref=fits.readKey(file,'MJDREF');

MJDREF=str2double(mjdref);

fits.movAbsHDU(file,2);

tstart = str2double(fits.readKey(file,'TSTART'));

fits.closeFile(file);

t_raw = fitsread(filoc,"binarytable");

%t_raw = fitsread('C:\Users\Filippo\Desktop\J1544_2022\SiFAP2_20220429_3FGLJ1544_N1_Bary.fits','binarytable');

t_raw = t_raw{1};

% MJDREF 1544

% MJDREF=59698.943055555602768;

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

toc

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

% 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 = (.009477-.2326620370367891)*86400; %MJD?

% Parameters for EPN_0744840201_bary.fits

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

% a_gr = zeros(5,1);

% tasc_gr = zeros(5,1);

% Parameters for 1544

Tseg=490; %segment length in seconds

MJDREF=59698.943055555602768;

f_gr = (100:1/Tseg:800).';

%f_gr=(772:1/Tseg:774).';

% porb_gr = (18000:10:21600).';

porb_gr=(18000:150:21600).';

a_gr = (0.01:4.0e-3:0.2).';

tasc_gr = (((59698.7231:8e-4:59698.9232) - MJDREF).*86400).';

% for j = 1:5

%     %f_gr(j) = f_tru*((j^2)/9);

%     porb_gr(j) = porb_tru*((j^2)/9);

%     a_gr(j) = a_tru*((j^2)/9);

%     tasc_gr(j) = tasc_tru*((j^2)/9);

% end

% Tseg=256; %segment length in seconds

% f_gr = (500:1/Tseg:600).';

% porb_gr=(17000:5:17200).';

% a_gr = (0.3:8.0e-4:0.4).';

% MJDREF = 58107.1236342593;

% tasc_gr = (((58106.9994:8e-4:58107.1095) - MJDREF).*86400).';

f_min = min(f_gr)

f_max = max(f_gr)

porb_min = min(porb_gr)

porb_max = max(porb_gr)

a_min = min(a_gr)

a_max = max(a_gr)

tasc_min = min(tasc_gr)

tasc_max = max(tasc_gr)

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

pathfi = [folname,'/'];

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

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

% FAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAALSE

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

tic

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

Nyq=2000; %cambiare in caso

Nyq=2000 %cambiare in caso

Tseg = 490

dt_psd=1/(2*Nyq); %risoluzione spettro di potenza, diverso dal dt dei segmenti

nbin=round((t_raw(end)-t_raw(1))/dt_psd);

countsname = sprintf('counts_check_Tseg_%d_.mat',Tseg);

%%

if((Nyq == 2000) & (exist([pathfi,countsname],"file")==2))

    load([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw");

    % clear t_raw

else

    tic

    [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

    % clear t_raw

    %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

    M=fix((t(end)-t(1))/Tseg); %number of segments

        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(N,M);

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

    toc

    clear C

% 

    save([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw","-v7.3","-nocompression");

end

%%

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

% 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 = (.009477-.2326620370367891)*86400; %MJD?

% Parameters for EPN_0744840201_bary.fits

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

% a_gr = zeros(5,1);

% tasc_gr = zeros(5,1);

% Parameters for 1544

% Tseg=490; %segment length in seconds

% MJDREF=59698.943055555602768;

% f_gr = (750:1/Tseg:800).';

% %f_gr=(772:1/Tseg:774).';

% % porb_gr = (18000:10:21600).';

% porb_gr=(18000:80:21600).';

% a_gr = (0.01:4.9e-3:0.5).';

% % porb_gr=(18000:100:21600).'; %f_max = 750, a_max = 0.5

% % a_gr = (0.01:6.0e-3:0.5).';

% tasc_gr = (((59698.7231:8e-4:59698.9232) - MJDREF).*86400).';

% for j = 1:5

%     %f_gr(j) = f_tru*((j^2)/9);

%     porb_gr(j) = porb_tru*((j^2)/9);

%     a_gr(j) = a_tru*((j^2)/9);

%     tasc_gr(j) = tasc_tru*((j^2)/9);

% end

% Tseg=256; %segment length in seconds

% f_gr = (500:1/Tseg:600).';

% porb_gr=(17000:5:17200).';

% a_gr = (0.3:8.0e-4:0.4).';

% MJDREF = 58107.1236342593;

% tasc_gr = (((58106.9994:8e-4:58107.1095) - MJDREF).*86400).';

% Parameters for Sco X-1

 % s 

f_max = 800

f_min = 750

f_step = 1/Tseg

f_gr = (f_min:f_step:f_max).';

% Porb = 68023.919 % s

% porb_gr = (Porb).'; % s

% porb_unc = 0.017; % s

% porb_max = 28800

% porb_min = 14400

% porb_step = 60

% porb_gr = (porb_min:porb_step:porb_max);

% Tasc_old = 56722.63037153067129629630; % MJD

Tasc_old = 57508.3323; % MJD

Porb = 0.2415361*86400.0

Norb=round((MJDREF-Tasc_old)/(Porb/86400.0));

% Norb=round((MJDREF-Tasc_old)/(Porb/86400.0));

Tasc_propagato = Tasc_old+Norb*(Porb/86400.0)

porb_unc =  0.0000036*86400.0;

% Tasc_old_unc = 37; % in sec

Tasc_old_unc = 0.0053*86400.0; % in sec

Tasc_new_unc = sqrt(Tasc_old_unc^2 + (Norb*porb_unc)^2); % s 

tasc_mid = (Tasc_propagato - MJDREF)*86400

tasc_min = tasc_mid - 2.0*Tasc_new_unc

tasc_max = tasc_mid + 2.0*Tasc_new_unc

% tasc_min = ((59698.7231 - MJDREF).*86400)

% tasc_max = (((59698.9232) - MJDREF).*86400)

% tasc_gr = (((58106.9494:8e-4:58107.1595) - MJDREF).*86400).';

% tasc_min = min(tasc_gr)

% tasc_max = max(tasc_gr)

tasc_mid = 0.5*(tasc_max+tasc_min)

a_min = 0.008

a_max = 0.15

nrot = ceil(f_max*Tseg);

toff = tstart - tasc_mid;

porb_max = 21600

porb_min = 18000

% porb_step = 60

tol = 0.9

stepparino = max(delporb(a_max,porb_min,f_max,toff,Tseg,M,tol),(porb_min^2)*sqrt(0.1*(1-tol))/(a_max*2*pi*Tseg*f_max))

porb_step = (porb_max-porb_min)/ceil((porb_max-porb_min)/stepparino)

porb_gr = (porb_min:porb_step:porb_max).';

stepparino = delasenc(porb_min,f_max,toff,Tseg,M,tol)

a_step = (a_max-a_min)/ceil((a_max-a_min)/stepparino)

% a_step = 1.0/(2.0*f_max*0.9)

%%

stepparino = max(deltasc(a_max,porb_min,f_max,toff,Tseg,M,tol),0.1025*porb_min/(pi*f_max*a_max*tol))

%%

tasc_step = (tasc_max-tasc_min)/ceil((tasc_max-tasc_min)/stepparino)

a_gr = (a_min:a_step:a_max).';

tasc_gr = (tasc_min:tasc_step:tasc_max).';

% trupar(1) = 592.42146827248556;

% trupar(2) = 17115.5216592;

% trupar(3) = 0.343356;

% trupar(4) = 58107.009477;

% f_min = min(f_gr)

% f_max = max(f_gr)

% porb_min = min(porb_gr)

% porb_max = max(porb_gr)

% a_min = min(a_gr)

% a_max = max(a_gr)

% tasc_min = min(tasc_gr)

% tasc_max = max(tasc_gr)

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

%% 

% Add step to generate par. extrema

param = [f_min,f_max,porb_min,porb_max,a_min,a_max,tasc_min,tasc_max];

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

% end

% Is (1-Ome*a) more than fmin-fmax? to ask myself

s = 1;

nismin(s) = f_max*(1 - a_max*Omega_max);

nismax(s) = f_max*(1 + a_max*Omega_max);

nismin(s) = (1 - a_max*Omega_max);

nismax(s) = (1 + a_max*Omega_max);

for s = 2:4

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

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

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

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

end

toc

%s_s=uint8(4);

s_s = 4;

while(1)

    if((nismax(s_s)-nismin(s_s))<delta_ni_new(s_s,Tseg))

    if((nismax(s_s)-nismin(s_s))<delta_ni_new(s_s,Tseg)/f_max)

        s_s=s_s-1;

    else

        break;

% element is accessed as nis{s}(i)

Nni = zeros(s_s,1);

for s = 1:s_s

    Nni(s) = ceil((nismax(s)-nismin(s))/delta_ni_new(s,Tseg));

    Nni(s) = ceil((nismax(s)-nismin(s))*f_max/delta_ni_new(s,Tseg));

    if (mod(Nni(s),2)==0)

        Nni(s) = Nni(s)+1;

    end

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

    nis{s}=franco/f_max;

    nis{s}=franco;

end

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

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

tm=gpuArray(tm);

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

Lambda = zeros(length(f_gr),length(nibank),'single','gpuArray');

% Lambda = zeros(length(f_gr),length(nibank),'single','gpuArray'); % MATLAB

% stores array elements in a column-major layout

Lambda = zeros(length(f_gr),length(nibank),'single');

F1=gpuArray((0:N-1)./(N*dt_psd)).';

f_ind = zeros(1,length(f_gr),'uint32');

for n = 1:length(f_gr)

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

taul = zeros(2,1);

% taul = zeros(N+1,1,'gpuArray');

ttempdifl = zeros(1,N+1,'gpuArray'); 

% taul = zeros(2,1);

disp(['length(nibank) = ',num2str(length(nibank)),'']);

nino = length(nibank);

s_s

M

avetime = 0.0;

guess = nino*M*5614.2461/(5376.0*48);

disp(['Rough estimate of the time necessary for the first loop = ',num2str(guess),'']);  disp(' ');

for m=1:M

    segstart = tic;

    % tic

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

 %   edges=edges(end)-edges(1); %mi dà il tempo preciso di tutta la TdF, che sarà leggermente diversa da length(C)*dt per come è definito histcounts 

    % ttemp = gpuArray(tm(m,:));

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

    ttemp = tm(m,:);

    % ttemp = tm(m,:);

    ttemp = gpuArray(t_raw((t_raw>=tm(m,1)-0.5*dt) & t_raw<=tm(m,N)+0.5*dt));

    Nphot = length(ttemp);

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

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

    ttempdifl(1:N) = (tm(m,1:N)-tmid(m)-0.5*dt).';

    ttempdifl(N+1) = (tm(m,N)-tmid(m)+0.5*dt);

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

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

    sumx = sum(xtemp);

    avex = mean(xtemp);

    % ttemp = tmrado(m,:);

    % tedge(1:N) = ttemp(:) - 0.5*dt;  

    % tedge(N+1) = ttemp(end) + 0.5*dt;

    tedge(1:N) = tm(m,1:N) - 0.5*dt;  

    tedge(N+1) = tm(m,N) + 0.5*dt;

    Y0 = fft(xtemp);

    norman = 2./mean((abs(Y0(F1>1000 & F1<Nyq)).^2)./sumx);

    % toc

    %FARE IL RICAMPIONAMENTO!!!!

    % Dev'essere fatto per ciascun template, ergo mettiamo un bel ciclo 

    % sulle possibili combinazioni di ni_s

    tic

    for i = 1:length(nibank)

    % tic

    for i = 1:nino

        % Ricampionamento (Eq. 17 MP2015)

        % Per prima cosa, generiamo la nuova coordinata temporale per 

        % questa templ combini 

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

        taul(1:2) = tmid(m) + sum((nibank(i,1:s_s)./(factorial(1:s_s))).*((ttempdifl(1:2)).^((1:s_s).').'),2);

        dtau = ((taul(2)-taul(1)));

        % inzo = (interp1(tau(:),(xtemp(:)./dtau),ttemp(:),'linear',0)).*dt;

        Y1 = fft((interp1(tau(:),(xtemp(:)./dtau),ttemp(:),'linear',0)).*dt).';

        tau(1:Nphot) = tmid(m) + sum((nibank(i,1:s_s)./(factorial(1:s_s))).*((ttempdif(1:Nphot)).^((1:s_s).').'),2);

        % taul(1:N+1) = tmid(m) + sum((nibank(i,1:s_s)./(factorial(1:s_s))).*((ttempdifl(1:N+1)).^((1:s_s).').'),2);

        % dtau(1:N) = (taul(2:N+1)-taul(1:N));

        % counvec = repelem(tau,xtemp);

        % [mat,~] = histcounts(counvec,tedge);

        [mat,~] = histcounts(tau,tedge);

        % Y0 = fft(xtemp);

        % inzo = (interp1(tau(:),(xtemp(:)./dtau),ttemp(:),'linear',0)).*dt;

        Y1 = fft(mat).';

        sumo = real(Y1(1));

        % plot(F1,abs(Y0)./sumx,F1,abs(Y1)./sumo)

        % pause

        % Careful! ni_0 still needs to be defined (it's the spin freq. 

        % currently being considered)

        Yenne(:)=Y1(f_ind(:)); 

        Lambda(:,i) = 2.*(abs(Yenne(:)).^2)./sumx;

        Lambda(:,i) = norman.*(abs(Yenne(:)).^2)./sumo;

        % 

        % for n = 1:length(f_gr)

        % 

    end

    % disp(max(Lambdacp(:,:)-Lambda(:,:)))

    % disp(min(Lambdacp(:,:)-Lambda(:,:)))

    toc

    disp(length(f_gr)*i);

    disp(m);

    % toc

    % disp(length(f_gr)*i);

    % disp(m);

    lamfiname = sprintf('Lambda_seg_%d.mat', m);

    % save(lamfiname,"Lambda");

    save([pathfi,lamfiname],"Lambda","-v7.3");

    Lambda = zeros(length(f_gr),length(nibank),'single','gpuArray');

    disp(m);

    save([pathfi,lamfiname],"Lambda","-v7.3","-nocompression");

    Lambda = zeros(length(f_gr),length(nibank),'single');

    % Lambda = zeros(length(f_gr),length(nibank),'single','gpuArray');

    segend = toc(segstart);

    avetime = avetime + segend;

    disp(['Just finished segment num.: ',num2str(m),'.']);

    disp(['Time elapsed in this loop = ',num2str(avetime),' s, estimated time left = ',num2str(((M-m)*avetime/m)),' s.']);  disp(' ');

end

%% 

%%

tic

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

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

%bestpar = zeros(1,5);

toc

disp(['length(parbank) = ',num2str(length(parbank)),'']);

nimask = false(length(nibank),1);

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

tic

% totlam = zeros(length(f_gr),length(parbank),'single');

totlam = ones(length(f_gr),length(parbank),'single');

guess = length(parbank)*M*29761.2957/(48.0*2180530.0);

disp(['Probably terrible estimate of the time necessary for the second loop = ',num2str(guess),'']);  disp(' ');

% tic

% for n=1:length(f_gr)

%     for i=1:length(parbank)

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

%     end

% end

bles = cospi(1/2);

avetime = 0.0;

for m = 1:M

    tic

    segstart = tic;

    lamfiname = sprintf('Lambda_seg_%d.mat', m);

    % load(lamfiname);

    load([pathfi,lamfiname]);

    toc

    disp(m);

    tic

    % toc

    % disp(m);

    % tic

    for i=1:length(parbank)

        curpar = [bles,parbank(i,:)];

        curni = -curni.*curpar(3);

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

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

        if (~nimask(Idx))

        % if (~nimask(Idx))

        nimask(Idx) = 1;

        end

        % end

        totlam(:,i) = totlam(:,i)+Lambda(:,Idx);

    end

    toc 

    disp(m);

    % toc 

    segend = toc(segstart);

    avetime = avetime + segend;

    disp(['Just finished segment num.: ',num2str(m),'.']);

    disp(['Elapsed time since loop began = ',num2str(avetime),' s, estimated time left = ',num2str(((M-m)*avetime/m)),' s.']); disp(' ');

    clear Lambda

end

% % for n=1:length(f_gr)

% %         % % end

% %     end

% % end

toc

% toc

%% 

totlam = totlam - 1;

%% Calculate detection threshold for a multi-trial false alarm prob. = pfa

pfa = 0.01;

%sigmastar = 2*gammaincinv((1-(1-pfa)^(1/(length(parbank)*length(f_gr)))),M,'upper');

sigmastar = 2*gammaincinv((1-(1-pfa)^(1/(sum(nimask)*length(f_gr)))),M,'upper');

sigmastar = 2*gammaincinv((1-(1-pfa)^(1/(length(parbank)*length(f_gr)))),M,'upper');

%sigmastar = 2*gammaincinv((1-(1-pfa)^(1/(sum(nimask)*length(f_gr)))),M,'upper');

[bru,tto] = max(totlam,[],2);

[maxtotlam,maxinde] = max(bru);

bestpar = [f_gr(maxinde), parbank(tto(maxinde),:), double(maxtotlam)]; 

bparch = string(bestpar);

tascmjd = MJDREF+tasc_gr./86400;

pfamax = 1.0-( (1.0-gammainc(0.5*double(maxtotlam),M,'upper'))^(length(parbank)*length(f_gr)));

disp('Our most significant peak has parameters:');

disp(['f_spin = '+bparch(1)+' Hz, p_orb = '+bparch(2)+' s, a*sin(i)/c = '+bparch(3)+' s, t_asc = '+bparch(4)+' s.']);

disp(['t_asc in days is = ',num2str((MJDREF+bestpar(4)/86400)),'.']);

disp(['Its detection statistics is Σ_max = ',num2str(maxtotlam),',']);

disp(['which considering all parameters, corresponds to a false alarm probability of ',num2str(pfamax),','])

disp(['while the multi-trial threshold for a false alarm probability = ',num2str(pfa),' is eq. to ',num2str(sigmastar),'.']);

sigmastar2 = 2*gammaincinv((1-(1-pfa)^(1/(sum(nimask)*length(f_gr)))),M,'upper');

pfamax2 = 1.0-( (1.0-gammainc(0.5*double(maxtotlam),M,'upper'))^(sum(nimask)*length(f_gr)));

disp(['When considering only the nis walked instead, Σ_max corresponds to a false alarm probability of ',num2str(pfamax2),','])

disp(['while the multi-trial threshold for a false alarm probability = ',num2str(pfa),' is eq. to ',num2str(sigmastar2),'.']);

% disp(bestpar);

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

save([pathfi,'risultelli_',char(datetime('now','Format','dd_MM')),'_',num2str(Tseg),'s.mat'],"-v7.3");

save([pathfi,'risultelli_test_',char(datetime('now','Format','dd_MM')),'_',num2str(Tseg),'s_amax',num2str(a_max),'_fmax',num2str(f_max),'.mat'],"-v7.3","-nocompression");

disp(['Results saved in ',pathfi,'risultelli_test_',char(datetime('now','Format','dd_MM')),'_',num2str(Tseg),'s_amax',num2str(a_max),'_fmax',num2str(f_max),'.mat']);

%% 

daterella = ['Date and time now ',char(datetime('now')),' (UTC), maybe'];

disp(daterella); clear daterella

diary off

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

%%%%%%%%%%%%% Fai il grafico delle totlam a posteriori %%%%%%%%%%%%%%%%

% contour(a_gr,tasc_gr,squeeze(max(provaccia,[],3)),[20:20:120])

% contour(a_gr,porb_gr,squeeze(max(provaccia,[],1)).',[20:20:120])

% provaccia = reshape(totlam.',length(parbank)*(length(f_gr)),1);

% provaccia = reshape(provaccia,length(tasc_gr),length(a_gr),length(porb_gr),length(f_gr));

% figstep = (10);

% figu = figure('visible','off');

% rina = "porb_tasc.m";

% contourf(porb_gr,tascmjd,squeeze(max(provaccia,[],[2 4])),[sigmastar*0.9:figstep:1.1*maxtotlam],'ShowText','on')

% hold on

% scatter(bestpar(2),(MJDREF+bestpar(4)/86400))

% if (exist('trupar') == 1)

%     scatter(trupar(2),trupar(4),'+')

% end

% hold off

% saveas(figu,rina);

% figu = figure('visible','off');

% rina = "porb_a.m";

% contourf(porb_gr,a_gr,squeeze(max(provaccia,[],[1 4])),[sigmastar:figstep:1.1*maxtotlam],'ShowText','on')

% hold on

% scatter(bestpar(2),bestpar(3))

% if (exist('trupar') == 1)

%     scatter(trupar(2),trupar(3),'+')

% end

% hold off

% saveas(figu,rina);

% figu = figure('visible','off');

% rina = "f_porb.m";

% contourf(f_gr,porb_gr,squeeze(max(provaccia,[],[1 2])),[sigmastar:figstep:1.1*maxtotlam],'ShowText','on')

% hold on

% scatter(bestpar(1),bestpar(2))

% if (exist('trupar') == 1)

    % scatter(trupar(1),trupar(2),'+')

% end

% hold off

% saveas(figu,rina);

exit

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

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

end

function res = delta_ni_new(s,Tcoh)

    %Compute eq. 23 M2015

    res = 0.5*(factorial(s-1))/(Tcoh^s);

end

function res = deltasc(asenc,porb,fspin,toff,tseg,M,eqtol)

    %Compute eq. 36 from Caliandro+2012

    Ome_orb = 2*pi/porb;

    omespin = 2*pi*fspin;

    nrot = ceil(fspin*tseg);

    pert = @(ind,step,m) (asenc*step*Ome_orb*cos(Ome_orb*(ind/fspin + toff + (m-1)*tseg)));

    % fun = @(step,m) ((sum(sin(omespin.*pert((0:nrot-1),step,m))).^2 + sum(cos(omespin.*pert((0:nrot-1),step,m))).^2)./(nrot).^2.0 - eqtol);

    res = 10.0^10;

    for m = 1:M

        fun = @(step) ((sum(sin(omespin.*pert((0:nrot-1),step,m))).^2 + sum(cos(omespin.*pert((0:nrot-1),step,m))).^2)./(nrot).^2.0 - eqtol);

        root = fzero(fun,[1 porb]);

        if (root<res)