Commit 1a6051e0 authored by Riccardo La Placa's avatar Riccardo La Placa
Browse files

Clean and done, fixed exporting size to A4

parent 645a9abf

powdist_f.m

+13 −9
Original line number Diff line number Diff line 
function [extrVoverE,CrossTalkProb,HalfMeanPow] = powdist_f(NyqFreq,LC,Binned,PSname,Noisename)

%POWDIST_F PS, rebinning and crosstalk estimate; graphs for visual inspection 

%

%   Input arguments:

%   NyqFreq is the Nyquist frequency in Hz 

%   LC is the lightcurve
@@ -7,6 +8,7 @@ function [extrVoverE,CrossTalkProb,HalfMeanPow] = powdist_f(NyqFreq,LC,Binned,PS 
%   otherwise we assume LC are actually the times of arrival

%   PSname and Noisename are the filenames (or complete paths) of graphs to

%   be saved, resp. PSs+lightcurve and counts dist. with crosstalk estimate

%   

%   Output arguments:

%   extrVoverE are the three values of the expected Var/E assuming 8-pixel,

%   4-pixel, and 1-pixel crosstalk probability distributions
@@ -22,7 +24,7 @@ dt = 1/(2*NyqFreq); 
Hfreq = 1000;



if (~Binned)

    N = fix((LC(end)-LC(1))/dt);

    N = round((LC(end)-LC(1))/dt);

    Tseg = dt*N;

    ted = (0:dt:Tseg);

    tm = (0.5*dt:dt:(Tseg-0.5*dt));
@@ -36,7 +38,7 @@ else 
    xtemp = LC;

end

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

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

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

nblog = 100;

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

for i = 1:length(Flog)-1
@@ -49,7 +51,7 @@ 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);

Y0 = abs(fft(xtemp));

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

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

Yqlbin = zeros(size(Fqlbin));

for nino = 1:nbqlog

    Yqlbin(nino) = sum(Y0((Fqlog(nino)<=F0)&(F0<Fqlog(nino+1))));
@@ -57,16 +59,18 @@ end 
Yqlbin(:) = Yqlbin(:)./(Fqlbin(:).*Tseg);







t = tiledlayout(3,1);

figure('Units','centimeters','Position', [0.75 0.75 21 29.7])

t = tiledlayout(3,1,'Padding','tight');

% t.Units = 'centimeters';

% t.OuterPosition = [0.75 0.75 29.7 21];

nexttile

semilog(F0,Y0)

semilogx(F0,Y0)

title('Power spectrum (PS) up to Nyq. freq.')

xlabel('Frequency [Hz]')

ylabel('Leahy-normalized power (no corr.)')



nexttile

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

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

semilogx(Fqlmid,Yqlbin)

xlim([10 NyqFreq+1000])

ylim([yme*0.8 yme*1.2])
@@ -102,7 +106,7 @@ if (isempty(xtemp(xtemp ==0))) 
            N = 2*N;

        end

        Y0 = abs(fft(xtemp));

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

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

    else

        disp('Data were provided already binned, nothing we can do')

        disp('Try with times of arrival instead, and set Binned to false')
@@ -124,7 +128,7 @@ mi = -log(ptot0); 
% mivalues(i,m) = mi;

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

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

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

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

vsuedat = var(xtemp)/mean(xtemp);



ncr = 8;