First complete version of multiobs! Now running to check for errors

    mkdir(figfolname);

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

end

oldpathfi = pathfi;

pathfi = [folname,'/'];

pathfigu = [figfolname,'/'];

tic

cd(oldpathfi)

mario = split(ls('*.fits'));

mario = sort(mario(~cellfun('isempty',mario)));

numfiles = length(mario);

toc

pathfi = [folname,'/'];

pathfi = [folname,'/']

Nyq = 2000 %cambiare in caso

Tseg = 256

C=C.'; %conteggi

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

M=fix((t_raw(end)-t_raw(1))/Tseg);

Mtot = M;

Mold = 0;

tm = zeros(M,N);

tmid = zeros(M,1);

for m = 1:M

%% Multiple obs. (that's what w'ere here for)

if (numfiles>1)

for i = 2:numfiles

    finame = char(mario(i))

    %% WILL have to choose wether to do only the parameter space check here

for o = 2:numfiles

    finame = char(mario(o))

    %% WILL have to choose whether to do only the parameter space check here

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

    file=fits.openFile(filoc);

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

    MJDREF(i)=str2double(mjdref);

    MJDREF(o)=str2double(mjdref);

    fits.movAbsHDU(file,2);

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

    t_raw = fitsread(filoc,"binarytable");

    C=C.'; %conteggi

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

    M=fix((t_raw(end)-t_raw(1))/Tseg);

    Mtot = Mtot + M;

    tm = zeros(M,N);

    tmid = zeros(M,1);

    for m = 1:M

%% ADD CYCLE ON OBS.

tm=gpuArray(tm);

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

nifax = gpuArray(nibank(:,1:s_s).*infax(1:s_s));

nino = length(nibank);

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

s_s

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

% stores array elements in a column-major layout

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

for n = 1:length(f_gr)

    [~,f_ind(n)] = min(abs(F1-f_gr(n)));

end

for o = 1:numfiles

    if (numfiles ~= 1)

        finame = char(mario(o))

        countsname = ['counts_',finame,'_Tseg_',num2str(Tseg),'.mat'];

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

    end

    tm=gpuArray(tm);

    tedge = ones(1,N+1,'gpuArray');

    M

    avetime = 0.0;

    guess = nino*M*Tseg*450/(440.0*31*256);

    disp(['A rough estimate of the time needed for the first loop in the current file (',num2str(o),' of ',num2str(numfiles),') is'])

    if (guess > 60)

        guess = guess/60;

        if (guess > 60)

            guess = guess/60;

            disp(['',num2str(ceil(guess)),' hours']);  disp(' ');

        else

            disp(['',num2str(ceil(guess)),' minutes']);  disp(' ');

        end

    else

        disp(['',num2str(ceil(guess)),' s']);  disp(' ');

    end

    %% The first loop involves calculating Lambda for each point in nibank

    for m=1:M

        segstart = tic;

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

        Nphot = length(ttemp);

        ttempdifs = zeros(Nphot,s_s,'gpuArray');

        ttempdifs(:,1) = (ttemp(:)-tmid(m)).';

        for s = 2:s_s

            ttempdifs(:,s) = ttempdifs(:,1).^s;

        end

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

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

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

        Y0 = fft(xtemp);

        graph1 = [pathfigu,'PS_',finame,'_Tseg',num2str(Tseg),'_m',num2str(m),'.pdf'];

        graph2 = [pathfigu,'CT_',finame,'_Tseg',num2str(Tseg),'_m',num2str(m),'.pdf'];

        if (~(isfile(graph1)) || ~(isfile(graph2)))

            [VsuEtot,epsct,hmp] = powdist(Nyq,ttemp,false,graph1,graph2);

        end

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

        tau = ones(size(ttemp),"like",ttemp);

        for i = 1:nino

            % Ricampionamento (Eq. 17 MP2015)

            % Per prima cosa, generiamo la nuova coordinata temporale per

            % questa templ combini

            tau(1:Nphot) = sum((nifax(i,1:s_s)).*(ttempdifs(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));

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

            Y1 = abs(fft(mat).');

            Lambda(:,i) = gather(Y1(f_ind(:))./sqrt(Y1(1)));

        end

        Lambda = (Lambda.^2).*norman;

        lamfiname = ['Lambda_',finame,'_fmax_',num2str(f),'_seg_',num2str(m),'.mat'];

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

        Lambda = ones(size(Lambda),"like",Lambda);

        % 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

    %% The first time in this loop we need to set some general variables such as niwalk and totlam 

    if (o == 1)

        nodesnum = length(porb_gr)*length(a_gr)*length(tasc_gr)*length(f_gr);

        nodesmem = (nodesnum*8.0)/1024^3; %->to bytes -> to GBs

        memthres = 600.0 - (length(f_gr)*parno*8.0)/1024^3;

        if (nodesmem > memthres)

            % THIS DISCLAIMER IS STILL HERE AFTER MONTHS

            % FIX IIIIT

            disp('Add part to split par. space!!!')

        end

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

        parno = length(parbank);

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

        if (nino<intmax('uint16'))

            niwalk = ones(Mtot,parno,'uint16');

        elseif (nino<intmax('uint32'))

            niwalk = ones(Mtot,parno,'uint32');

        elseif (nino<intmax('uint64'))

            niwalk = ones(Mtot,parno,'uint64');

        else

            disp(['The length of nibank in greater than ',num2str(intmax('uint64')),', something definitely makes no sense here'])

        end

        % Writing on a matrix of ones is faster in MATLAB than on one made 

        % of zeroes, we shall subtract one after all calculations are done

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

    end

    if (o == 1)

        mjdiff = 0;

    else

        mjdiff = (MJDREF(o) - MJDREF(1))*86400.0;

    end

    guess = parno*M*16192.4539/(4.0*7546032.0);

    disp(['(Probably terrible) Estimate of the time needed for the second loop in the current file (',num2str(o),' of ',num2str(numfiles),') =']);

    if (guess > 60)

        guess = guess/60;

        if (guess > 60)

            guess = guess/60;

            disp(['',num2str(ceil(guess)),' hours']);  disp(' ');

        else

            disp(['',num2str(ceil(guess)),' minutes']);  disp(' ');

        end

    else

        disp(['',num2str(ceil(guess)),' s']);  disp(' ');

    end

    bles = cospi(1/2);

    avetime = 0.0;

    for m = 1:M

        segstart = tic;

        lamfiname = ['Lambda_',finame,'_fmax_',num2str(f),'_seg_',num2str(m),'.mat'];

        load([pathfi,lamfiname]);

        % toc

        % disp(m);

        % tic

        for i=1:length(parbank)

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

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

            curpar(2) = 2*pi/curpar(2);

            curpar(4) = curpar(2)*(tmid(m) - parbank(i,3) + mjdiff);

            curni=zeros(1,s_s);

            for s=1:s_s

                curni(s) = (curpar(2)^s)*sin(curpar(4)+0.5*s*pi);

            end

            curni = -curni.*curpar(3);

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

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

            [schi3,fo3] = min(abs(nibank(:,:)-curni),[],1);

            Idx = sum(fo3)-length(fo3)+1;

            % if (~nimask(Idx))

            % end

            niwalk(m+Mold,i) = Idx;

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

        end

        % 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

    Mold = Mold + M;

end

% totlam was initialized as ones

totlam = totlam - 1;

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

pfa = 0.01;

sigmastar = 2*gammaincinv((1-(1-pfa)^(1/(parno*length(f_gr)))),Mtot,'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-( (gammainc(0.5*double(maxtotlam),Mtot,'lower'))^(parno*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 (',num2str(parno),') 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),'.']);

truwalk = size(unique(niwalk.', 'rows'),1);

sigmastar2 = 2*gammaincinv(((1-pfa)^(1/(truwalk*length(f_gr)))),Mtot,'lower');

pfamax2 = 1.0-( (gammainc(0.5*double(maxtotlam),Mtot,'lower'))^(truwalk*length(f_gr)));

disp(['The number of combinations of ni combinations actually walked is ',num2str(truwalk),',']);

% 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(['therefore, when considering just them 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),'.']);

memused = whos;

summina = 0;

for k = 1:length(memused)

summina = summina + memused(k).bytes;

end

clear memused

%inKB = memused/1024;

%inMB = memused/1024^2;

inGB = summina/1024^3;

%fprintf("The workspace will need %.5f KB or %.5f MB or %.5f GB of memory", inKB, inMB, inGB)

timerapp = 24.0*summina/518260001575;

% disp(bestpar);

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

disp(['I am about to save around ',num2str(round(inGB)),' GB of data (+ headers), it should take roughly ',num2str(ceil(timerapp)),' minutes.']);

save([pathfi,'risultelli_tot_',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_tot_',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

%--------------------------------------------------------------------------

%FUNCTIONS

%--------------------------------------------------------------------------

function res = delta_ni(mu_s,s,g)

    %Compute eq. 23 M2015

    res= 2.*sqrt(mu_s./(s.*g));

end

function res = delta_ni_new(s,Tcoh)

    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)

            res = root;

        end

    end

end

function res = delasenc(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) (step*sin(Ome_orb*(ind/fspin + toff + (m-1)*tseg)));

    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,[0.00001 100]);

        if (root<res)

            res = root;

        end

    end

end

function res = delporb(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.*(ind./fspin + toff + (m-1).*tseg).*(cos(Ome_orb.*(ind./fspin + toff + (m-1).*tseg)))./(porb)));

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

            res = root;

        end

    end

end

function deltaUTC2TT=UTC2TT_HST(MJD_REF)

MJDleap=[41499;41683;42048;42413;42778;43144;43509;43874;44239;44786;45151;45516;46247;47161;47892;48257;48804;49169;49534;50083;50630;51179;53736;54832;56109;57204;57754];

nMJDleap=length(MJDleap);

a=find(MJDleap<=MJD_REF);

if(isempty(a))

    leap=nMJDleap+10;    

else

    leap=length(a)+10;

end

deltaUTC2TT=32.184+leap;

end