Commit 573d299b authored by Riccardo La Placa's avatar Riccardo La Placa
Browse files

Merge commit 'dc668865'

parents 8291ff04 dc668865

functions/fastcounts.m

+1 −1
Original line number Diff line number Diff line
@@ -2,7 +2,7 @@ function counts = fastcounts(times,tedges) 
%FASTCOUNTS Histogram of array of times in bins with tedges (use on CPU!)

%   For longish arrays it's much faster than hisctounts and uses only a

%   fraction of its memory; operating on sorted times is where it shines

%   but it still is faster in  unsorted ones. 

%   but it still is faster on unsorted ones. 

%   The tedges can be non-uniformly spaced but the algorithm doesn't

%   rescale counts based on width: if you want that use histcounts

power_dist_test.m

+19 −9
Original line number Diff line number Diff line 
clear

close all



cd /data/Sorgenti/power_dist_test/

% cd /data/Sorgenti/power_dist_test/

diary(['log_pow_dist_',char(datetime('now','Format','dd_MM')),'.log'])

%%

mario = split(ls('*.fits'));
@@ -18,7 +18,7 @@ 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 = 100;

nblog = 200;

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/Tseg-(Flog(i))));
@@ -62,10 +62,14 @@ for i = 1:length(mario) 
            end

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



            if (i == 1)

                tiledlayout(2,1)

            % if (i == 1)

                tiledlayout(3,1)

                nexttile





                semilogx(F0,Y0)

                xlim([10 3000])

                nexttile 

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

                semilogx(Fqlmid,Yqlbin)

                xlim([10 3000])
@@ -77,8 +81,8 @@ for i = 1:length(mario) 
                xlim([ted(1) ted(end)])

                plot(tm,gather(xtemp))



                % pause

            end

                pause

            % end





            if (~isempty(xtemp(xtemp ==0)))
@@ -109,6 +113,9 @@ for i = 1:length(mario) 
                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)));

                if maxcount<5

                    p8_1 = p8_1(1:maxcount);

                end

                for k = 6:maxcount

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

                end
@@ -137,6 +144,9 @@ for i = 1:length(mario) 
                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)));

                if maxcount<5

                    p4_1 = p4_1(1:maxcount);

                end

                for k = 6:maxcount

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

                end
@@ -172,14 +182,14 @@ for i = 1:length(mario) 
                % if (i == 7)

                tiledlayout(2,1)

                nexttile

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

                semilogy([1:maxcount],extrtot1,[1:maxcount],extrtot4,[1:maxcount],extrtot8,[1:maxcount],extrpoiss,[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,'+')

                semilogy([1:maxcount],abs(extrtot1-datadist)./datadist,'+',[1:maxcount],abs(extrtot4-datadist)./datadist,'+',[1:maxcount],abs(extrtot8-datadist)./datadist,'+',[1:maxcount],abs(extrpoiss-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])

                % pause

                pause

                % end

            else

                disp('Too many counts!')