Commit 4528c9b9 authored by Riccardo La Placa's avatar Riccardo La Placa
Browse files

Steps forward, var/exp check

parent b19afa1d

power_dist_test.m

+122 −85
Original line number Diff line number Diff line 
clear

close all



cd /data/Sorgenti/SCORPIUS-X-1/

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));
@@ -15,16 +15,17 @@ tm = (0.5*dt:dt:(Tseg-0.5*dt)); 
F1=gpuArray((0:N-1)./(N*dt)).';

F0 = F1(1:N/2+1);

% F02 = F0(end/2:end);

nblog = 200;

nblog = 100;

Flog = logspace(log10(F1(2)/2),log10(F0(end)+0.5/Tseg),nblog);

for i = 1:length(Flog)-1

[~,minind] = min(abs(F0(:)-0.5/512-(Flog(i))));

[~,minind] = min(abs(F0(:)-0.5/Tseg-(Flog(i))));

Fqlog(i) = F0(minind)-0.5/Tseg;

end

Fqlog(nblog) = F0(end)+0.5/Tseg;

Fqlog = unique(Fqlog);

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

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

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

for i = 1:length(mario)

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

    % file = fits.openFile(filoc)
@@ -32,6 +33,7 @@ for i = 1:length(mario) 
    t_raw = fitsread(filoc,"binarytable");

    t_raw = t_raw{1};

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

    t_raw = t_raw - t_raw(1);

    Nphot = length(t_raw);

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

    xtemp = gpuArray(xtemp);
@@ -39,89 +41,124 @@ for i = 1:length(mario) 
    Y0 = 2.0.*(Y0(1:N/2+1).^2)./Y0(1);

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

    Yqlbin = zeros(size(Fqlbin));

    for nino = 1:length(Fqlbin)

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

    for nino = 1:nbqlog

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

    end

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

    semilogx(Flmid,Yqlbin,''+'')

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

    xlim([100 3000])

    semilogx(Fqlmid,Yqlbin)

    yline(mean(Y0(100*Tseg+1:end)))

    pause

end





% xtemp = x(:,1);

% length(xtemp(xtemp == 0))

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

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

% mean(xtemp)

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

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

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

        mi = -log(ptot0);

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

% ncross = 8;

% crn = 8;

        ncr = 8;

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

% clear creps ncross crn

        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/8)

% expec8L = 1*(1-epscr) + 2*(8*pcr*qcr^14) + 3*(84*(pcr^2)*(qcr^20)) + 4*24*(pcr^3)*(qcr^24)*(1+3*qcr+38*qcr^2)

% expec8L = 1*(1-epscr) + 2*(8*pcr*qcr^14) + 3*(84*(pcr^2)*(qcr^20))

% expec8L = 1*(1-epscr) + 2*(8*pcr*qcr^14) + 3*(84*(pcr^2)*(qcr^20)) + 4*24*(pcr^3)*(qcr^24)*(1+3*qcr+38*qcr^2) + 4*(pcr^4)*(qcr^30)*(180+540*qcr+2521*qcr2)

% expec8L = 1*(1-epscr) + 2*(8*pcr*qcr^14) + 3*(84*(pcr^2)*(qcr^20)) + 4*24*(pcr^3)*(qcr^24)*(1+3*qcr+38*qcr^2) + 4*(pcr^4)*(qcr^30)*(180+540*qcr+2521*qcr^2)

% p8l1 = (1-epscr);

% p8l2 = (8*pcr*qcr^14);

% p8l3 = (84*(pcr^2)*(qcr^20));

        expectot = expec1*mi

        vartot = var1*mi + mi*expec1^2

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

        p8l(1) = (1-epscr);

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

        p8l(3) = (84*(pcr^2)*(qcr^20));

        p8l(4) = 24*(pcr^3)*(qcr^24)*(1+3*qcr+38*qcr^2);

        p8l(5) = 4*(pcr^4)*(qcr^30)*(180+540*qcr+2521*qcr^2);

r8l = p8l(5)/(1-sum(p8l(1:4)))

% expec8L = 1*(1-epscr) + 2*(8*pcr*qcr^14) + 3*(84*(pcr^2)*(qcr^20)) + 4*24*(pcr^3)*(qcr^24)*(1+3*qcr+38*qcr^2) + p8l(5)*(1+4*r8l)/r8l^2

% sum(p8l(:))

        r8l = p8l(5)/(1-sum(p8l(1:4)));

        for k = 6:10

            p8l(k) = p8l(5)*(1-r8l)^(k-5);

        end

        vard = var(xtemp)

        vartot/expectot

        vard/mean(xtemp)

        mean(Y0(9*end/10:end))/2

    else

        disp('Too many counts!')

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

    end

    pause

end



% load('counts_check_Tseg_490_.mat')

% Y0 = 2.0.*(abs(fft(x(:,1))).^2)./sum(x(:,1));

% F1=((0:N-1)./(N*dt)).';

% mean(Y0)

% mean(Y0(9*N/20:N/2))

% Y0sub = Y0 - 2.300035862611959

% Y0sub = Y0 - 0.300035862611959;

% Y0div = Y0./(ans);

% mean(Y0div(9*N/20:N/2))

% Y0div = 2.0*Y0./(ans);

% Y0div = 2.0*Y0./(mean(Y0(9*N/20:N/2)));

% F0(9*N/20)

% F1(9*N/20)

% F1(N/2)

% F1(N/2+1)

% F0 = F1(2:N/2+1);

% Y0 = Y0(2:N/2+1);

% Y0div = Y0div(2:N/2+1);

% Y0sub = Y0sub(2:N/2+1);

% loglog(F0,Y0,F0,Y0sub,F0,Y0div)

% semilogy(F0,Y0,F0,Y0sub,F0,Y0div)

% yline(2)

% max(Y0)

% max(Y0div)

% max(Y0sub)

% powedges = (0.0:1.0:52.0);

% [histno,~] = histcounts(Y0,powedges,'normalization','pdf');

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

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

% Y0div2 = Y0div(end/2:end);

% Y0sub2 = Y0sub(end/2:end);

% F02 = F0(end/2:end);

% [histno,~] = histcounts(Y02,powedges,'normalization','pdf');

% [histnodiv,~] = histcounts(Y0div2,powedges,'normalization','pdf');

% [histnosub,~] = histcounts(Y0sub2,powedges,'normalization','pdf');

% powmid(1:52) = (powedges(2:53)+powedges(1:52))./2;

% semilogy(powmid,histno,'o',powmid,histnodiv,'+',powmid,chi2pdf(powmid,(2.0)),powmid,histnosub,'*')
 
 No newline at end of file
 


% % xtemp = x(:,1);

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

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

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

% % mean(xtemp)

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

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

% mi = -log(ptot0);

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

% % ncross = 8;

% % crn = 8;

% ncr = 8;

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

% % clear creps ncross crn

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

% pcr = 1-qcr;

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

% 

% % expec8L = 1*(1-epscr) + 2*(8*pcr*qcr^14) + 3*(84*(pcr^2)*(qcr^20)) + 4*24*(pcr^3)*(qcr^24)*(1+3*qcr+38*qcr^2)

% % expec8L = 1*(1-epscr) + 2*(8*pcr*qcr^14) + 3*(84*(pcr^2)*(qcr^20))

% % expec8L = 1*(1-epscr) + 2*(8*pcr*qcr^14) + 3*(84*(pcr^2)*(qcr^20)) + 4*24*(pcr^3)*(qcr^24)*(1+3*qcr+38*qcr^2) + 4*(pcr^4)*(qcr^30)*(180+540*qcr+2521*qcr2)

% % expec8L = 1*(1-epscr) + 2*(8*pcr*qcr^14) + 3*(84*(pcr^2)*(qcr^20)) + 4*24*(pcr^3)*(qcr^24)*(1+3*qcr+38*qcr^2) + 4*(pcr^4)*(qcr^30)*(180+540*qcr+2521*qcr^2)

% % p8l1 = (1-epscr);

% % p8l2 = (8*pcr*qcr^14);

% % p8l3 = (84*(pcr^2)*(qcr^20));

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

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

% p8l(3) = (84*(pcr^2)*(qcr^20));

% p8l(4) = 24*(pcr^3)*(qcr^24)*(1+3*qcr+38*qcr^2);

% p8l(5) = 4*(pcr^4)*(qcr^30)*(180+540*qcr+2521*qcr^2);

% r8l = p8l(5)/(1-sum(p8l(1:4)))

% % expec8L = 1*(1-epscr) + 2*(8*pcr*qcr^14) + 3*(84*(pcr^2)*(qcr^20)) + 4*24*(pcr^3)*(qcr^24)*(1+3*qcr+38*qcr^2) + p8l(5)*(1+4*r8l)/r8l^2

% % sum(p8l(:))

% for k = 6:10

% p8l(k) = p8l(5)*(1-r8l)^(k-5);

% end

% 

% % load('counts_check_Tseg_490_.mat')

% % Y0 = 2.0.*(abs(fft(x(:,1))).^2)./sum(x(:,1));

% % F1=((0:N-1)./(N*dt)).';

% % mean(Y0)

% % mean(Y0(9*N/20:N/2))

% % Y0sub = Y0 - 2.300035862611959

% % Y0sub = Y0 - 0.300035862611959;

% % Y0div = Y0./(ans);

% % mean(Y0div(9*N/20:N/2))

% % Y0div = 2.0*Y0./(ans);

% % Y0div = 2.0*Y0./(mean(Y0(9*N/20:N/2)));

% % F0(9*N/20)

% % F1(9*N/20)

% % F1(N/2)

% % F1(N/2+1)

% % F0 = F1(2:N/2+1);

% % Y0 = Y0(2:N/2+1);

% % Y0div = Y0div(2:N/2+1);

% % Y0sub = Y0sub(2:N/2+1);

% % loglog(F0,Y0,F0,Y0sub,F0,Y0div)

% % semilogy(F0,Y0,F0,Y0sub,F0,Y0div)

% % yline(2)

% % max(Y0)

% % max(Y0div)

% % max(Y0sub)

% % powedges = (0.0:1.0:52.0);

% % [histno,~] = histcounts(Y0,powedges,'normalization','pdf');

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

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

% % Y0div2 = Y0div(end/2:end);

% % Y0sub2 = Y0sub(end/2:end);

% % F02 = F0(end/2:end);

% % [histno,~] = histcounts(Y02,powedges,'normalization','pdf');

% % [histnodiv,~] = histcounts(Y0div2,powedges,'normalization','pdf');

% % [histnosub,~] = histcounts(Y0sub2,powedges,'normalization','pdf');

% % powmid(1:52) = (powedges(2:53)+powedges(1:52))./2;

% % semilogy(powmid,histno,'o',powmid,histnodiv,'+',powmid,chi2pdf(powmid,(2.0)),powmid,histnosub,'*')
 
 No newline at end of file