Commit 698c369e authored by Riccardo La Placa's avatar Riccardo La Placa
Browse files

Scappo

parent 0a1d1edc

power_dist_test.m

+201 −176
Original line number Diff line number Diff line
@@ -5,8 +5,9 @@ cd /data/Sorgenti/power_dist_test/ 
diary(['log_pow_dist_',char(datetime('now','Format','dd_MM')),'.log'])

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

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

mario = sort(mario);

import matlab.io.*

Nyq = 4000

Nyq = 2000

Tseg = 512

dt = 1/(2*Nyq);

N = fix(Tseg/dt);
@@ -26,6 +27,10 @@ Fqlog = unique(Fqlog); 
nbqlog = length(Fqlog)-1;

Fqlmid(1:nbqlog) = 0.5*(Fqlog(2:end)+Fqlog(1:end-1));

Fqlbin(1:nbqlog) = Fqlog(2:end)-Fqlog(1:end-1);

epsmedio = 0;

mivalues = zeros(length(mario),100);

epsvalues = zeros(length(mario),100);

varvalues = zeros(length(mario),100,3);

for i = 1:length(mario)

    filoc = [pwd,'/',char(mario(i))]

    % file = fits.openFile(filoc)
@@ -50,6 +55,7 @@ for i = 1:length(mario) 
        end

        Yqlbin(:) = Yqlbin(:)./(Fqlbin(:).*Tseg);



        if (i == 4 || i == 6 || i == 7)

            tiledlayout(2,1)

            nexttile
@@ -65,6 +71,7 @@ for i = 1:length(mario) 
            plot(tm,gather(xtemp))



            pause

        end





        if (~isempty(xtemp(xtemp ==0)))
@@ -74,16 +81,19 @@ for i = 1:length(mario) 
                datadist(bu) = length(xtemp(xtemp == bu))/N;

            end

            ptot0 = length(xtemp(xtemp == 0))/N;

            mi = -log(ptot0)

            mi = -log(ptot0);

            mivalues(i,m) = mi;

            extrpoiss = poisspdf([1:maxcount],mi);

            ncr = 8;

            epscr = 1 - (length(xtemp(xtemp == 1))/N)/(mi*ptot0)

            epscr = 1 - (length(xtemp(xtemp == 1))/N)/(mi*ptot0);

            epsvalues(i,m) = epscr;

            epsmedio = epsmedio + epscr/M;

            pcr = 1- (1-epscr)^(1/ncr);

            expec1 = 1+ epscr + (3*(ncr-1)/(2*ncr))*epscr^2;

            var1 = epscr + (epscr^2)*(7-9/ncr)/2;

            enf = 1 + var1/(expec1^2);

            expectot = expec1*mi

            vartot = var1*mi + mi*expec1^2

            expectot = expec1*mi;

            vartot = var1*mi + mi*expec1^2;

            qcr = (1-epscr)^(1/ncr);

            p8_1 = zeros(1,maxcount);

            p8_1(1) = (1-epscr);
@@ -95,12 +105,12 @@ for i = 1:length(mario) 
            for k = 6:maxcount

                p8_1(k) = p8_1(5)*(1-r8_1)^(k-5);

            end

            vard = var(xtemp)

            vsuetot = vartot/expectot

            vsuedat = vard/mean(xtemp)

            meanpow = mean(Y0(100*Tseg+1:end))/2

            meanpow/vsuetot

            meanpow/vsuedat

            vard = var(xtemp);

            vsuetot = vartot/expectot;

            vsuedat = vard/mean(xtemp);

            meanpow = mean(Y0(100*Tseg+1:end))/2;

            varvalues(i,m,1) = meanpow/vsuetot;

            % meanpow/vsuedat

            extrtot8 = ptotk(p8_1,extrpoiss);
@@ -110,8 +120,8 @@ for i = 1:length(mario) 
            expec1 = 1+ epscr + (3*(ncr-1)/(2*ncr))*epscr^2;

            var1 = epscr + (epscr^2)*(7-9/ncr)/2;

            enf = 1 + var1/(expec1^2);

            expectot = expec1*mi

            vartot = var1*mi + mi*expec1^2

            expectot = expec1*mi;

            vartot = var1*mi + mi*expec1^2;

            qcr = (1-epscr)^(1/ncr);

            p4_1 = zeros(1,maxcount);

            p4_1(1) = (1-epscr);
@@ -123,12 +133,12 @@ for i = 1:length(mario) 
            for k = 6:maxcount

                p4_1(k) = p4_1(5)*(1-r4_1)^(k-5);

            end

            vard = var(xtemp)

            vsuetot = vartot/expectot

            vsuedat = vard/mean(xtemp)

            meanpow = mean(Y0(100*Tseg+1:end))/2

            meanpow/vsuetot

            meanpow/vsuedat

            vard = var(xtemp);

            vsuetot = vartot/expectot;

            vsuedat = vard/mean(xtemp);

            meanpow = mean(Y0(100*Tseg+1:end))/2;

            varvalues(i,m,2) = meanpow/vsuetot;

            % meanpow/vsuedat

            extrtot4 = ptotk(p4_1,extrpoiss);



            ncr = 1;
@@ -137,19 +147,21 @@ for i = 1:length(mario) 
            expec1 = 1+ epscr + (3*(ncr-1)/(2*ncr))*epscr^2;

            var1 = epscr + (epscr^2)*(7-9/ncr)/2;

            enf = 1 + var1/(expec1^2);

            expectot = expec1*mi

            vartot = var1*mi + mi*expec1^2

            expectot = expec1*mi;

            vartot = var1*mi + mi*expec1^2;

            qcr = (1-epscr)^(1/ncr);

            p1_1 = zeros(1,maxcount);

            p1_1(1:maxcount) = (pcr.^(0:maxcount-1)).*(1-pcr);

            vard = var(xtemp)

            vsuetot = vartot/expectot

            vsuedat = vard/mean(xtemp)

            meanpow = mean(Y0(100*Tseg+1:end))/2

            meanpow/vsuetot

            meanpow/vsuedat

            vard = var(xtemp);

            vsuetot = vartot/expectot;

            vsuedat = vard/mean(xtemp);

            meanpow = mean(Y0(100*Tseg+1:end))/2;

            varvalues(i,m,3) = meanpow/vsuetot;

            % meanpow/vsuedat

            extrtot1 = ptotk(p1_1,extrpoiss);





            if (i == 4 || i == 6 || i == 7)

            tiledlayout(2,1)

            nexttile

            semilogy([1:maxcount],extrtot1,[1:maxcount],extrtot4,[1:maxcount],extrtot8,[1:maxcount],datadist,'+')
@@ -159,150 +171,163 @@ for i = 1:length(mario) 
            ylim([0.0001 100])

            % legend('abs(extrtot1-datadist)/datadist','abs(extrtot4-datadist)/datadist','abs(extrtot8-datadist)/datadist')

            yline([0.01 0.1 0.5 1.0])

        else

            disp('Too many counts!')

            disp(['Avg. counts per bin :',num2str(mean(xtemp)),'.']);

        end

            pause

            end

end



%% 

filoc = ['/data/Sorgenti/power_dist_test/dark.mat']

load(filoc)

% file = fits.openFile(filoc)

% fits.closeFile(file);

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

t_raw = ToAs2;

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

for m = 1:M

    disp(['Segment ',num2str(m),' of ',num2str(M),' in file',filoc,' '])

    tsplit = t_raw((t_raw>(t_raw(1)+Tseg*(m-1)))&(t_raw<(t_raw(1)+Tseg*m)));

    % t_raw = t_raw(t_raw<(t_raw(1)+Tseg));

    tsplit = tsplit - tsplit(1);

    Nphot = length(tsplit);

    [xtemp,~]=(histcounts(tsplit,ted));

    xtemp = gpuArray(xtemp);

    Y0 = abs(fft(xtemp));

    Y0 = 2.0.*(Y0(1:N/2+1).^2)./Y0(1);

    % Y02 = Y0(end/2:end);

    Yqlbin = zeros(size(Fqlbin));

    for nino = 1:nbqlog

        Yqlbin(nino) = sum(Y0((Fqlog(nino)<=F0)&(F0<Fqlog(nino+1))));

    end

    Yqlbin(:) = Yqlbin(:)./(Fqlbin(:).*Tseg);



    tiledlayout(1,1)

    nexttile



    yme = mean(Y0(100*Tseg+1:end));

    semilogx(Fqlmid,Yqlbin)

    xlim([10 3000])

    ylim([yme*0.8 yme*1.2])

    yline(yme)

    xline(100)



    pause





    if (~isempty(xtemp(xtemp ==0)))

        maxcount = max(xtemp);

        datadist = zeros(1,maxcount);

        for bu = 1:maxcount

            datadist(bu) = length(xtemp(xtemp == bu))/N;

        end

        ptot0 = length(xtemp(xtemp == 0))/N;

        mi = -log(ptot0)

        extrpoiss = poisspdf([1:maxcount],mi);

        ncr = 8;

        epscr = 1 - (length(xtemp(xtemp == 1))/N)/(mi*ptot0)

        pcr = 1- (1-epscr)^(1/ncr);

        expec1 = 1+ epscr + (3*(ncr-1)/(2*ncr))*epscr^2;

        var1 = epscr + (epscr^2)*(7-9/ncr)/2;

        enf = 1 + var1/(expec1^2);

        expectot = expec1*mi

        vartot = var1*mi + mi*expec1^2

        qcr = (1-epscr)^(1/ncr);

        p8_1 = zeros(1,maxcount);

        p8_1(1) = (1-epscr);

        p8_1(2) = (8*pcr*qcr^14);

        p8_1(3) = (12*(pcr^2)*(qcr^18)*(1+2*qcr+4*qcr^2));

        p8_1(4) = 4*(pcr^3)*(qcr^20)*(1+3*qcr+14*qcr^2+30*qcr^3+61*qcr^4+59*qcr^5+72*qcr^6);

        p8_1(5) = 5*(pcr^4)*(qcr^24)*(9+36*qcr+98*qcr^2+188*qcr^3+310*qcr^4+372*qcr^5+520*qcr^6+396*qcr^7+341*qcr^8);

        r8_1 = p8_1(5)/(1-sum(p8_1(1:4)));

        for k = 6:maxcount

            p8_1(k) = p8_1(5)*(1-r8_1)^(k-5);

        end

        vard = var(xtemp)

        vsuetot = vartot/expectot

        vsuedat = vard/mean(xtemp)

        meanpow = mean(Y0(100*Tseg+1:end))/2

        meanpow/vsuetot

        meanpow/vsuedat

        extrtot8 = ptotk(p8_1,extrpoiss);





        ncr = 4;

        epscr = 1 - (length(xtemp(xtemp == 1))/N)/(mi*ptot0);

        pcr = 1- (1-epscr)^(1/ncr);

        expec1 = 1+ epscr + (3*(ncr-1)/(2*ncr))*epscr^2;

        var1 = epscr + (epscr^2)*(7-9/ncr)/2;

        enf = 1 + var1/(expec1^2);

        expectot = expec1*mi

        vartot = var1*mi + mi*expec1^2

        qcr = (1-epscr)^(1/ncr);

        p4_1 = zeros(1,maxcount);

        p4_1(1) = (1-epscr);

        p4_1(2) = (4*pcr*qcr^6);

        p4_1(3) = (18*(pcr^2)*(qcr^8));

        p4_1(4) = 4*(pcr^3)*(qcr^8)*(1+3*qcr+18*qcr^2);

        p4_1(5) = 5*(pcr^4)*(qcr^10)*(8+24*qcr+55*qcr^2);

        r4_1 = p4_1(5)/(1-sum(p4_1(1:4)));

        for k = 6:maxcount

            p4_1(k) = p4_1(5)*(1-r4_1)^(k-5);

        end

        vard = var(xtemp)

        vsuetot = vartot/expectot

        vsuedat = vard/mean(xtemp)

        meanpow = mean(Y0(100*Tseg+1:end))/2

        meanpow/vsuetot

        meanpow/vsuedat

        extrtot4 = ptotk(p4_1,extrpoiss);



        ncr = 1;

        epscr = 1 - (length(xtemp(xtemp == 1))/N)/(mi*ptot0);

        pcr = 1- (1-epscr)^(1/ncr);

        expec1 = 1+ epscr + (3*(ncr-1)/(2*ncr))*epscr^2;

        var1 = epscr + (epscr^2)*(7-9/ncr)/2;

        enf = 1 + var1/(expec1^2);

        expectot = expec1*mi

        vartot = var1*mi + mi*expec1^2

        qcr = (1-epscr)^(1/ncr);

        p1_1 = zeros(1,maxcount);

        p1_1(1:maxcount) = (pcr.^(0:maxcount-1)).*(1-pcr);

        vard = var(xtemp)

        vsuetot = vartot/expectot

        vsuedat = vard/mean(xtemp)

        meanpow = mean(Y0(100*Tseg+1:end))/2

        meanpow/vsuetot

        meanpow/vsuedat

        extrtot1 = ptotk(p1_1,extrpoiss);



        tiledlayout(2,1)

        nexttile

        semilogy([1:maxcount],extrtot1,[1:maxcount],extrtot4,[1:maxcount],extrtot8,[1:maxcount],datadist,'+')

        legend('extrtot1','extrtot4','extrtot8','datadist')

        nexttile

        semilogy([1:maxcount],abs(extrtot1-datadist)./datadist,[1:maxcount],abs(extrtot4-datadist)./datadist,[1:maxcount],abs(extrtot8-datadist)./datadist)

        ylim([0.0001 100])

        legend('abs(extrtot1-datadist)/datadist','abs(extrtot4-datadist)/datadist','abs(extrtot8-datadist)/datadist')

        yline([0.01 0.1 0.5 1.0])

        else

            disp('Too many counts!')

            disp(['Avg. counts per bin :',num2str(mean(xtemp)),'.']);

        end

    pause

        % pause

    end

end

epsmedio = epsmedio/length(mario)

epsnans = epsvalues;

epsnans(epsnans == 0) = NaN;

% controllato che tutte quelle sopra 0.1 sono dovute a buchi

% nell'osservazione

epsnans(epsnans > 0.1) = NaN;

minans = mivalues;

minans(minans == 0) = NaN;

%% 

% disp('Finished with the fits files, move to dark obs.?')

% pause

% filoc = ['/data/Sorgenti/power_dist_test/dark.mat']

% load(filoc)

% % file = fits.openFile(filoc)

% % fits.closeFile(file);

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

% %%

% t_raw = ToAs2;

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

% for m = 1:M

%     disp(['Segment ',num2str(m),' of ',num2str(M),' in file',filoc,' '])

%     tsplit = t_raw((t_raw>(t_raw(1)+Tseg*(m-1)))&(t_raw<(t_raw(1)+Tseg*m)));

%     % t_raw = t_raw(t_raw<(t_raw(1)+Tseg));

%     tsplit = tsplit - tsplit(1);

%     Nphot = length(tsplit);

%     [xtemp,~]=(histcounts(tsplit,ted));

%     xtemp = gpuArray(xtemp);

%     Y0 = abs(fft(xtemp));

%     Y0 = 2.0.*(Y0(1:N/2+1).^2)./Y0(1);

%     % Y02 = Y0(end/2:end);

%     Yqlbin = zeros( ...

%         size(Fqlbin));

%     for nino = 1:nbqlog

%         Yqlbin(nino) = sum(Y0((Fqlog(nino)<=F0)&(F0<Fqlog(nino+1))));

%     end

%     Yqlbin(:) = Yqlbin(:)./(Fqlbin(:).*Tseg);

% 

%     tiledlayout(1,1)

%     nexttile

% 

%     yme = mean(Y0(100*Tseg+1:end));

%     semilogx(Fqlmid,Yqlbin)

%     xlim([10 3000])

%     ylim([yme*0.8 yme*1.2])

%     yline(yme)

%     xline(100)

% 

%     pause

% 

% 

%     if (~isempty(xtemp(xtemp ==0)))

%         maxcount = max(xtemp);

%         datadist = zeros(1,maxcount);

%         for bu = 1:maxcount

%             datadist(bu) = length(xtemp(xtemp == bu))/N;

%         end

%         ptot0 = length(xtemp(xtemp == 0))/N;

%         mi = -log(ptot0)

%         extrpoiss = poisspdf([1:maxcount],mi);

%         ncr = 8;

%         epscr = 1 - (length(xtemp(xtemp == 1))/N)/(mi*ptot0)

%         pcr = 1- (1-epscr)^(1/ncr);

%         expec1 = 1+ epscr + (3*(ncr-1)/(2*ncr))*epscr^2;

%         var1 = epscr + (epscr^2)*(7-9/ncr)/2;

%         enf = 1 + var1/(expec1^2);

%         expectot = expec1*mi

%         vartot = var1*mi + mi*expec1^2

%         qcr = (1-epscr)^(1/ncr);

%         p8_1 = zeros(1,maxcount);

%         p8_1(1) = (1-epscr);

%         p8_1(2) = (8*pcr*qcr^14);

%         p8_1(3) = (12*(pcr^2)*(qcr^18)*(1+2*qcr+4*qcr^2));

%         p8_1(4) = 4*(pcr^3)*(qcr^20)*(1+3*qcr+14*qcr^2+30*qcr^3+61*qcr^4+59*qcr^5+72*qcr^6);

%         p8_1(5) = 5*(pcr^4)*(qcr^24)*(9+36*qcr+98*qcr^2+188*qcr^3+310*qcr^4+372*qcr^5+520*qcr^6+396*qcr^7+341*qcr^8);

%         r8_1 = p8_1(5)/(1-sum(p8_1(1:4)));

%         for k = 6:maxcount

%             p8_1(k) = p8_1(5)*(1-r8_1)^(k-5);

%         end

%         vard = var(xtemp)

%         vsuetot = vartot/expectot

%         vsuedat = vard/mean(xtemp)

%         meanpow = mean(Y0(100*Tseg+1:end))/2

%         meanpow/vsuetot

%         meanpow/vsuedat

%         extrtot8 = ptotk(p8_1,extrpoiss);

% 

% 

%         ncr = 4;

%         epscr = 1 - (length(xtemp(xtemp == 1))/N)/(mi*ptot0);

%         pcr = 1- (1-epscr)^(1/ncr);

%         expec1 = 1+ epscr + (3*(ncr-1)/(2*ncr))*epscr^2;

%         var1 = epscr + (epscr^2)*(7-9/ncr)/2;

%         enf = 1 + var1/(expec1^2);

%         expectot = expec1*mi

%         vartot = var1*mi + mi*expec1^2

%         qcr = (1-epscr)^(1/ncr);

%         p4_1 = zeros(1,maxcount);

%         p4_1(1) = (1-epscr);

%         p4_1(2) = (4*pcr*qcr^6);

%         p4_1(3) = (18*(pcr^2)*(qcr^8));

%         p4_1(4) = 4*(pcr^3)*(qcr^8)*(1+3*qcr+18*qcr^2);

%         p4_1(5) = 5*(pcr^4)*(qcr^10)*(8+24*qcr+55*qcr^2);

%         r4_1 = p4_1(5)/(1-sum(p4_1(1:4)));

%         for k = 6:maxcount

%             p4_1(k) = p4_1(5)*(1-r4_1)^(k-5);

%         end

%         vard = var(xtemp)

%         vsuetot = vartot/expectot

%         vsuedat = vard/mean(xtemp)

%         meanpow = mean(Y0(100*Tseg+1:end))/2

%         meanpow/vsuetot

%         meanpow/vsuedat

%         extrtot4 = ptotk(p4_1,extrpoiss);

% 

%         ncr = 1;

%         epscr = 1 - (length(xtemp(xtemp == 1))/N)/(mi*ptot0);

%         pcr = 1- (1-epscr)^(1/ncr);

%         expec1 = 1+ epscr + (3*(ncr-1)/(2*ncr))*epscr^2;

%         var1 = epscr + (epscr^2)*(7-9/ncr)/2;

%         enf = 1 + var1/(expec1^2);

%         expectot = expec1*mi

%         vartot = var1*mi + mi*expec1^2

%         qcr = (1-epscr)^(1/ncr);

%         p1_1 = zeros(1,maxcount);

%         p1_1(1:maxcount) = (pcr.^(0:maxcount-1)).*(1-pcr);

%         vard = var(xtemp)

%         vsuetot = vartot/expectot

%         vsuedat = vard/mean(xtemp)

%         meanpow = mean(Y0(100*Tseg+1:end))/2

%         meanpow/vsuetot

%         meanpow/vsuedat

%         extrtot1 = ptotk(p1_1,extrpoiss);

% 

%         tiledlayout(2,1)

%         nexttile

%         semilogy([1:maxcount],extrtot1,[1:maxcount],extrtot4,[1:maxcount],extrtot8,[1:maxcount],datadist,'+')

%         legend('extrtot1','extrtot4','extrtot8','datadist')

%         nexttile

%         semilogy([1:maxcount],abs(extrtot1-datadist)./datadist,[1:maxcount],abs(extrtot4-datadist)./datadist,[1:maxcount],abs(extrtot8-datadist)./datadist)

%         ylim([0.0001 100])

%         legend('abs(extrtot1-datadist)/datadist','abs(extrtot4-datadist)/datadist','abs(extrtot8-datadist)/datadist')

%         yline([0.01 0.1 0.5 1.0])

%     else

%         disp('Too many counts!')

%         disp(['Avg. counts per bin :',num2str(mean(xtemp)),'.']);

%     end

%     pause

% end

% % xtemp = x(:,1);

% % length(xtemp(xtemp == 0))

% % length(xtemp(xtemp == 0))/N

scsearch.m

+34 −29
Original line number Diff line number Diff line
@@ -5,8 +5,8 @@ import matlab.io.* 


maxNumCompThreads(48);



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

% 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)
@@ -31,38 +31,43 @@ tic 
% finame = 'ScoX1_20230423_Bary.fits'

% finame = 'J1023_B_2017_Bary.fits'

% finame = 'EPN_0744840201_bary.fits'

% finame = 'ScoX1_20220628_N2_Targ_Bary.fits'

finame = 'ofba01010_tag_bary_lambda_165_310_chan994_1006.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);

% fits.movAbsHDU(file,2);

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

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

% t_raw = t_raw{1};



%%%%%%% HST CASE

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

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

MJDREF=str2double(mjdref);

fits.movAbsHDU(file,2);

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

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

% 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



%%%%%%%



fits.closeFile(file);
@@ -85,7 +90,7 @@ pathfi = [folname,'/']; 
t0 = (t_raw(end)+t_raw(1))/2;

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

Nyq = 2000 %cambiare in caso

Tseg = 530

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