small efficiency improvements?

% Parameters for Sco X-1

 % s 

f_max = 250

f_min = 100

f_max = 1000

f_min = 950

f_step = 1/Tseg

f_gr = (f_min:f_step:f_max).';

% Porb = 68023.919 % s

% porb_gr = (Porb).'; porb_max = Porb; porb_min = Porb; % s

Porb = 68023.919 % s

porb_gr = (Porb).'; porb_max = Porb; porb_min = Porb; % s

% porb_unc = 2*0.017; % s

porb_max = 21600

porb_min = 18000

% porb_max = 21600

% porb_min = 18000

% porb_step = 60

% porb_gr = (porb_min:porb_step:porb_max);

% Tasc_old = 56722.6303715306713; % MJD

Tasc_old = 57508.3323; % MJD

Porb = 0.2415361*86400.0

Norb=round((MJDREF-Tasc_old)/(Porb/86400.0));

Tasc_old = 56722.6303715306713; % MJD

% Tasc_old = 57508.3323; % MJD

% Porb = 0.2415361*86400.0

Norb=round((MJDREF-Tasc_old+1)/(Porb/86400.0));

% Norb=round((MJDREF-Tasc_old)/(Porb/86400.0));

Tasc_propagato = Tasc_old+Norb*(Porb/86400.0)

porb_unc =  0.0000036*86400.0;

% Tasc_old_unc = 2*33 + 4; % in sec

Tasc_old_unc = 0.0053*86400.0; % in sec

% porb_unc =  0.0000036*86400.0;

porb_unc = 0.017*2;

Tasc_old_unc = 2*33 + 4; % in sec

% Tasc_old_unc = 0.0053*86400.0; % in sec

Tasc_new_unc = sqrt(Tasc_old_unc^2 + (Norb*porb_unc)^2); % s 

tasc_mid = (Tasc_propagato - MJDREF)*86400

tasc_min = tasc_mid - 2.0*Tasc_new_unc

tasc_max = tasc_mid + 2.0*Tasc_new_unc

if ((tasc_max - tasc_min)>porb_max)

    disp("Range of uncertainties on T_asc larger than P_orb!")

    tasc_min = tasc_mid -0.5*porb_max

    tasc_max = tasc_mid +0.5*porb_max

end

% tasc_min = ((59698.7231 - MJDREF).*86400)

% tasc_max = (((59698.9232) - MJDREF).*86400)

% tasc_gr = (((58106.9494:8e-4:58107.1595) - MJDREF).*86400).';

% tasc_min = min(tasc_gr)

% tasc_max = max(tasc_gr)

% tasc_mid = 0.5*(tasc_max+tasc_min)

a_min = 0.0008

a_max = 0.26

% a_min = 1.44

% a_max = 3.25

% a_min = 0.0008

% a_max = 0.26

a_min = 1.44

a_max = 3.25

nrot = ceil(f_max*Tseg);

toff = tstart - tasc_mid;

% porb_max = 21600

% porb_min = 18000

% porb_step = 60

tol = 0.9

delporb(a_max,porb_min,f_max,toff,Tseg,M,tol)

(porb_min^2)*sqrt(0.1*(1-tol))/(a_max*2*pi*Tseg*f_max)

stepparino = max(delporb(a_max,porb_min,f_max,toff,Tseg,M,tol),(porb_min^2)*sqrt(0.1*(1-tol))/(a_max*2*pi*Tseg*f_max))

porb_step = (porb_max-porb_min)/ceil((porb_max-porb_min)/stepparino)

porb_gr = (porb_min:porb_step:porb_max).';

% porb_gr = (porb_min:porb_step:porb_max).';

%%

stepparino = delasenc(porb_min,f_max,toff,Tseg,M,tol)

delasenc(porb_min,f_max,toff,Tseg,M,tol)

1.0/(2.0*f_max*tol)

%%

a_step = (a_max-a_min)/ceil((a_max-a_min)/stepparino)

% a_step = 1.0/(2.0*f_max*0.9)

%%

deltasc(a_max,porb_min,f_max,toff,Tseg,M,tol)

0.1025*porb_min/(pi*f_max*a_max*tol)

stepparino = max(deltasc(a_max,porb_min,f_max,toff,Tseg,M,tol),0.1025*porb_min/(pi*f_max*a_max*tol))

%%

tasc_step = (tasc_max-tasc_min)/ceil((tasc_max-tasc_min)/stepparino)

% tasc_max = max(tasc_gr)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nodesnum = length(porb_gr)*length(a_gr)*length(tasc_gr)*length(f_gr)

%% 

% Add step to generate par. extrema

avetime = 0.0;

guess = nino*M*3071.5583/(3100.0*28);

disp(['Rough estimate of the time needed for the first loop = ',num2str(guess),' s']);  disp(' ');

%%

for m=1:M

    segstart = tic;

    ttempdifl(N+1) = (tm(m,N)-tmid(m)+0.5*dt);

    % xtemp = (x(:,m).');

    xtemp = gpuArray(x(:,m).');

    sumx = sum(xtemp);

    avex = mean(xtemp);

    % ttemp = tmrado(m,:);

    % tedge(1:N) = ttemp(:) - 0.5*dt;  

    % tedge(N+1) = ttemp(end) + 0.5*dt;

    tedge(1:N) = tm(m,1:N) - 0.5*dt;  

    tedge(N+1) = tm(m,N) + 0.5*dt;

    tedge(1:N) = tm(m,1:N) - 0.5*dt -tmid(m);  

    tedge(N+1) = tm(m,N) + 0.5*dt -tmid(m);

    Y0 = fft(xtemp);

    norman = 2./mean((abs(Y0(F1>1000 & F1<Nyq)).^2)./sumx);

    % norman = 2./mean((abs(Y0(F1>1000 & F1<Nyq)).^2)./sumx);

    norman = 2./(var(xtemp)/mean(xtemp));

    % toc

        % Ricampionamento (Eq. 17 MP2015)

        % Per prima cosa, generiamo la nuova coordinata temporale per 

        % questa templ combini 

        tau(1:Nphot) = tmid(m) + sum((nibank(i,1:s_s)./(factorial(1:s_s))).*((ttempdif(1:Nphot)).^((1:s_s).').'),2);

        % tau(1:Nphot) = tmid(m) + sum((nibank(i,1:s_s)./(factorial(1:s_s))).*((ttempdif(1:Nphot)).^((1:s_s).').'),2);

        tau(1:Nphot) = sum((nibank(i,1:s_s)./(factorial(1:s_s))).*((ttempdif(1:Nphot)).^((1:s_s).').'),2);

        % taul(1:N+1) = tmid(m) + sum((nibank(i,1:s_s)./(factorial(1:s_s))).*((ttempdifl(1:N+1)).^((1:s_s).').'),2);

        % dtau(1:N) = (taul(2:N+1)-taul(1:N));

        % Y0 = fft(xtemp);

        % inzo = (interp1(tau(:),(xtemp(:)./dtau),ttemp(:),'linear',0)).*dt;

        Y1 = fft(mat).';

        sumo = real(Y1(1));

        % sumo = real(Y1(1));

        % plot(F1,abs(Y0)./sumx,F1,abs(Y1)./sumo)

        % pause

        % Careful! ni_0 still needs to be defined (it's the spin freq. 

        % currently being considered)

        Yenne(:)=Y1(f_ind(:)); 

        Lambda(:,i) = norman.*(abs(Yenne(:)).^2)./sumo;

        Lambda(:,i) = (abs(Yenne(:)).^2);

        % Lambda(:,i) = norman.*(abs(Yenne(:)).^2)./sumo;

        % 

        % for n = 1:length(f_gr)

        % 

        % end

        % disp('1ni')

    end

    Lambda(:,:) = Lambda(:,:).*norman;

    % disp(max(Lambdacp(:,:)-Lambda(:,:)))

    % disp(min(Lambdacp(:,:)-Lambda(:,:)))

    % toc

memthres = 44.0 - (length(f_gr)*length(nibank)*8.0)/1024^3;

if (nodesmem > memthres)

    disp('Add part to split par. space!!!')

end

end

deltaUTC2TT=32.184+leap;

end