New branch (to test MATLAB's commit too)

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

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

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

figfolname = [folname,'/figures'];

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

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

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

        disp(['Folder ',figfolname,' already exists, previous graphs (if present) may be lost']);

    else

        mkdir(figfolname);

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

    end

else

    mkdir(folname);

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

    mkdir(figfolname);

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

end

pathfi = [folname,'/'];

pathfigu = [figfolname,'/'];

reset(gpuDevice(1));

gpuDevice(1);

finame = 'ScoX1_20220628_N2_Targ_Bary.fits'

% finame = 'ofba01010_tag_bary_lambda_165_310_chan994_1006.fits';

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

file=fits.openFile(filoc);

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

% % file=fits.openFile(filoc);

%%%%%% SiFAP CASE

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

MJDREF=str2double(mjdref);

fits.movAbsHDU(file,2);

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

t_raw = fitsread(filoc,"binarytable");

t_raw = t_raw{1};

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

% % MJDREF=str2double(mjdref);

% % fits.movAbsHDU(file,2);

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

% % t_raw = fitsread(filoc,"binarytable");

% % t_raw = t_raw{1};

%%%%%%% HST CASE

%%%%%%%

fits.closeFile(file);

% % fits.closeFile(file);

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

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;

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

Nyq = 2000 %cambiare in caso

Tseg = 256

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

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

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

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

    % clear t_raw

else

    tic

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

    t_raw = fitsread(filoc,"binarytable");

    t_raw = t_raw{1};

    fits.closeFile(file);

    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

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

    toc

    clear C

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

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

end

%%

Tasc_old = 56722.6303715306713; % MJD

% Tasc_old = 57508.3323; % MJD

% Porb = 0.2415361*86400.0

Norb=round((MJDREF-Tasc_old+1)/(Porb/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;

porb_unc = 0.017*2;

nifax = nibank(:,1:s_s).*infax(1: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','gpuArray');

Lambda = ones(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)

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

    % tedge = (tm(m,1)-1.5*dt-tmid)

    Y0 = fft(xtemp);

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

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

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

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

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

    end

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

    % norman = 2./(Nphot.*var(xtemp)./mean(xtemp));

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

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

        % currently being considered)

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

        Lambda(:,i) = Y1(f_ind(:)); 

        Lambda(:,i) = gather(Y1(f_ind(:))); 

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

        % 

        % for n = 1:length(f_gr)