Added HST files support

maxNumCompThreads(48);

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

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

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

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

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

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

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

% finame = 'ScoX1_20230423_Bary.fits'

% finame = 'J1023_B_2017_Bary.fits'

% finame = 'EPN_0744840201_bary.fits'

finame = 'ScoX1_20220628_N2_Targ_Bary.fits'

% finame = 'ScoX1_20220628_N2_Targ_Bary.fits'

finame = 'ofba01010_tag_bary_lambda_165_310_chan994_1006.fits';

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

file=fits.openFile(filoc);

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

MJDREF=str2double(mjdref);

% 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

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

MJDREF=str2double(mjdref)+UTC2TT_HST(str2double(mjdref))/86400;

mjdref=num2str(MJDREF);

RA=fits.readKey(file,'RA_TARG');

DEC=fits.readKey(file,'DEC_TARG');

obj_name=fits.readKey(file,'TARGNAME');

datamode=fits.readKey(file,'OBSMODE');

fits.movAbsHDU(file,2);

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

fits.closeFile(file);

tunit1=fits.readKey(file,'TUNIT1');

timedel=str2double(fits.readKey(file,'TSCAL1'));

tcn=fits.getColName(file,'TIME');

ToAs=fits.readCol(file,tcn); clear tcn  

timesys='TDB';

ephfile='DE200';

tstart=(ToAs(1));

tstop=(ToAs(end));

n_events=length(ToAs);

t_raw = ToAs; clear ToAs

t_raw = fitsread(filoc,"binarytable");

t_raw = t_raw{1};

fits.closeFile(file);

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

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

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

Nyq = 2000 %cambiare in caso

Tseg = 512

Tseg = 530

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

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

% Parameters for Sco X-1

 % s 

f_max = 750

f_min = 700

f_max = 600

f_min = 550

f_step = 1/Tseg

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

Porb = 68023.919 % s

porb_gr = (Porb).'; porb_max = Porb; porb_min = Porb; % s

porb_unc = 2*0.017; % s

% porb_max = 28800

% porb_min = 14400

% Porb = 68023.919 % s

% porb_gr = (Porb).'; porb_max = Porb; porb_min = Porb; % s

% porb_unc = 2*0.017; % s

porb_max = 21600

porb_min = 18000

% porb_step = 60

% porb_gr = (porb_min:porb_step:porb_max);

Tasc_old = 56722.6303715306713; % MJD

% Tasc_old = 57508.3323; % MJD

% Porb = 0.2415361*86400.0

% Tasc_old = 56722.6303715306713; % 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 = 2*33 + 4; % in sec

% Tasc_old_unc = 0.0053*86400.0; % in sec

% Tasc_old_unc = 2*33 + 4; % 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_min = min(tasc_gr)

% tasc_max = max(tasc_gr)

% tasc_mid = 0.5*(tasc_max+tasc_min)

a_min = 1.44

a_max = 3.25

a_min = 0.0008

a_max = 0.26

% a_min = 1.44

% a_max = 3.25

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

M

avetime = 0.0;

guess = nino*M*3071.5583/(3100.0*28);

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

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

for m=1:M

    segstart = tic;

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

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

guess = length(parbank)*M*16192.4539/(4.0*7546032.0);

disp(['(Probably terrible) Estimate of the time needed for the second loop = ',num2str(guess),' s']);  disp(' ');

bles = cospi(1/2);

avetime = 0.0;

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

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

pfamax2 = 1.0-( (1.0-gammainc(0.5*double(maxtotlam),M,'upper'))^(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(['When considering only the nis walked instead, Σ_max corresponds to a false alarm probability of ',num2str(pfamax2),','])

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('totlam').bytes;

%inKB = memused/1024;

%inMB = memused/1024^2;

inGB = memused/1024^3;

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

timerapp = 44.0*memused/670891259696.0;

% disp(bestpar);

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

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

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

%% 

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;

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

        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