Great strides but I still have to finish

function [outputArg1,outputArg2] = powdist_f(NyqFreq,BinLC)

%POWDIST_F Summary of this function goes here

function [outputArg1,outputArg2] = powdist_f(NyqFreq,LC,Binned,PSname,Noisename)

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

%   Detailed explanation goes here

%   NyqFreq is the Nyquist frequency 

%   BinLC is the binned lightcurve (with dt assumed = 1/2*NyqFreq)

import matlab.io.*

%   NyqFreq is the Nyquist frequency in Hz 

%   LC is the lightcurve 

%   If Binned is true, we assume uniform binning and dt = 1/2*NyqFreq,

%   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

% import matlab.io.* prob./ unneeded

%% Part I: Power spectrum, quasi-log. rebinning and light-curve printing

dt = 1/(2*NyqFreq);

N = length(BinLC);

if (~Binned)

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

    Tseg = dt*N;

    ted = (0:dt:Tseg);

    tm = (0.5*dt:dt:(Tseg-0.5*dt));

    LC = LC-LC(1);

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

else

    N = length(LC);

    Tseg = dt*N;

    ted = (0:dt:Tseg);

    tm = (0.5*dt:dt:(Tseg-0.5*dt));

    xtemp = LC;

end

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

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

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/Tseg-(Flog(i))));

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

end

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

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

Y0 = abs(fft(xtemp));

Y0 = 2.0.*(Y0(1: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))));

end

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

t = tiledlayout(3,1);

nexttile

semilog(F0,Y0)

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

xlabel('Frequency [Hz]')

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

nexttile

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

semilogx(Fqlmid,Yqlbin)

xlim([10 3000])

ylim([yme*0.8 yme*1.2])

yline(yme,'--','Avg. pow. for f > 100 Hz')

xline(100,':')

xlabel('Frequency [Hz]')

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

title('Quasi-logarithmically-rebinned PS')

nexttile

xlim([ted(1) ted(end)])

plot(tm,xtemp)

xlabel('Time [s]')

ylabel('Counts per bin')

title('Lightcurve (binned as provided)')

exportgraphics(t,PSname)

%% Part II: Noise and crosstalk evaluation

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

    disp(['Too many counts per bin to estimate the underlying '...

        'distribution, which we assume is Poissonian'])

    disp('(is this you, Sco X-1?)')

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

        ' and no bin with 0 counts.']);

    if(~Binned)

        disp('Will rebin the times of arrival on a finer grid')

        while (length(xtemp(xtemp ==0))/N<exp(-1))

            dt = dt/2;

            ted = (0:dt:Tseg);

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

            N = 2*N;

        end

        Y0 = abs(fft(xtemp));

        Y0 = 2.0.*(Y0(1: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')

    end

end

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

    % mivalues(i,m) = mi;

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

    ncr = 8;

    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;

    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

    expec1 = p8_1(5)*(1+4*r8_1)/r8_1/r8_1 + sum((1:4).*p8_1(1:4));

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

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

    expectot = expec1*mi;

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

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

    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;

    % varvalues(i,m,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;

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

    % meanpow/vsuedat

    extrtot1 = ptotk(p1_1,extrpoiss);

    if (i == 7)

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

        pause

    end

outputArg1 = inputArg1;

outputArg2 = inputArg2;

end

function res = ptotk(pncr_1,extrprim)

    %Calculate crosstalk prob. as in Gallego+13

    % if (length(pncr_1)<maxcount)

    %     disp('p8_1 was not calculated for the whole distribution of observed counts')

    %     disp(['Setting the size of p8_m to that (',num2str(length(pncr_1)),') instead of ',num2str(maxcount),'.'])

    %     maxcount = length(pncr_1);

    % end

    maxcount = length(pncr_1);

    pncr_m = zeros(maxcount);

    pncr_m(1,:) = pncr_1(:);

    for m = 1:maxcount-1

        for k = 1:maxcount

            for i = 1:k-m

                pncr_m(m+1,k) = pncr_m(m+1,k) + pncr_m(m,k-i)*pncr_1(i);

            end

            % pncr_m(m+1,k) = sum(pncr_m(m,k-(1:k-m)).*pncr_1(1:k-m));

        end

    end

    res = zeros(size(pncr_1));

    for k = 1:maxcount

        for m = 1:k

            res(k) = res(k) + extrprim(m).*pncr_m(m,k);

        end

    end

end

 No newline at end of file

    lamfiname = sprintf('Lambda_fmax_%d_seg_%d.mat', f_max, m);

    save([pathfi,lamfiname],"Lambda","-v7.3","-nocompression");

    Lambda = ones(length(f_gr),length(nibank),'single','gpuArray');

    Lambda = ones(size(Lambda),"like",Lambda);

    % Lambda = zeros(length(f_gr),length(nibank),'single','gpuArray');

    segend = toc(segstart);

    avetime = avetime + segend;