Commit 45988d17 authored by Riccardo La Placa's avatar Riccardo La Placa
Browse files

Merge branch 'devel' into 'main' 

A lot of changes, now multiobs and multifreq automatically

See merge request heag-oar/scsearch!1
parents e67a54cd cf97617a

functions/fastcounts.m

0 → 100644
+44 −0
Original line number Diff line number Diff line 
function counts = fastcounts(times,tedges)

%FASTCOUNTS Histogram of array of times in bins with tedges (use on RAM!)

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

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

%   but it still is faster in  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



if (~issorted(times))

    times = sort(times);

end

Nedges = length(tedges);

Nbins = length(tedges)-1;

Nphot = length(times);

counts = zeros(size(tedges));

counts = counts(1:Nbins);

% It appears matlab is faster with doubles anyway, we can take the

% negligible hit in memory

% counts = zeros(size(tedges),'int16');

% counts = counts(1:Nbins);

% if (Nphot<intmax('int16'))

%     counts = zeros(size(counts),'int16');

% elseif (Nphot<intmax('int32'))

%     counts = zeros(size(counts),'int32');

% elseif (Nphot<intmax('int64'))

%     counts = zeros(size(counts),'int64');

% else

%     disp(['The length of times is greater than ',num2str(intmax('int64')),', nothing will work!'])

% end



j = 1;

i = 1;

while times(j) < tedges(1)

    j=j+1;

end

while i < Nedges

    while (j<=Nphot && times(j)<tedges(i+1))

        counts(i) = counts(i)+1;

        j = j+1;

    end

    i = i+1;

end



end
 
 No newline at end of file

functions/lambdasum.m

0 → 100644
+64 −0
Original line number Diff line number Diff line 
function [outputArg1,outputArg2] = lambdasum(inputArg1,inputArg2,ParBank,NiSearcher)

%lambdasum Sums the Lambdas for all parameter combinations over obs.s/seg.s

%   Detailed explanation goes here





% 

% ParBank = combinations(porb_gr,a_gr,tasc_gr).Variables;



parno = length(ParBank);

freqno = length(f_gr);

nino = length(NiSearcher.X);

if (nino<intmax('uint16'))

    niwalk = ones(M,parno,'uint16');

elseif (nino<intmax('uint32'))

    niwalk = ones(M,parno,'uint32');

elseif (nino<intmax('uint64'))

    niwalk = ones(M,parno,'uint64');

else

    disp(['The length of nibank in greater than ',num2str(intmax('uint64')),', something definitely makes no sense here'])

end



disp(['length(parbank) = ',num2str(parno),'']);

% Writing on matrix of ones is faster in MATLAB than on zeroes, we subtract

% one after all calculations are done

totlam = ones(freqno,parno,'single');

guess = parno*M*16192.4539/(4.0*7546032.0);

disp(['(Probably terrible) Estimate of the time needed for the second loop = ',num2str(ceil(guess)),' s']);  disp(' ');



avetime = 0.0;

for m = 1:M

    segstart = tic;

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

    load([pathfi,lamfiname],"Lambda","");

    for i=1:parno

        

        curpar = ParBank(i,:);

        % Move from ν,P,a,T_asc to ν,Ω,a,γ

        % Move from P,a,T_asc to Ω,a,γ

        curpar(1) = 2*pi/curpar(1);

        curpar(3) = curpar(1)*(tmid(m) - curpar(3));

        curni=zeros(1,s_s);

        for s=1:s_s

            curni(s) = (curpar(1)^s)*sin(curpar(3)+0.5*s*pi);

        end

        curni = -curni.*curpar(2);

        curni(1) = curni(1) + 1;



        [Idx,~] = knnsearch(nisearcher,curni);



        niwalk(m,i) = Idx;



        totlam(:,i) = totlam(:,i)+Lambda(:,Idx);



    end

    segend = toc(segstart);

    avetime = avetime + segend;

    disp(['Just finished segment num.: ',num2str(m),'.']);

    disp(['Elapsed time since loop began = ',num2str(avetime),' s, estimated time left = ',num2str(((M-m)*avetime/m)),' s.']); disp(' ');

    clear Lambda

end



outputArg1 = inputArg1;

outputArg2 = inputArg2;

end
 
 No newline at end of file

functions/powdist.m

+1 −0
Original line number Diff line number Diff line
@@ -191,6 +191,7 @@ t = tiledlayout(2,1); 
nexttile

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

legend({'1-pixel crosstalk','4-pixel crosstalk','8-pixel crosstalk','Pure Poissonian','Data distribution'},'Location','southwest')

text(0.8,0.870,['$\epsilon_{ct} = $',num2str(epscr)])

xlabel('Counts')

ylabel('Number of bins with x counts')

title('Extrapolated generalized Poisson distributions vs data')

functions/power_upper_limit.m

0 → 100644
+24 −0
Original line number Diff line number Diff line 
function p_ul = power_upper_limit(max_power,n_ps,prob)

%POWER_UPPER_LIMIT Upper limit on signal at given prob. for measured power

%   This function calculates "the [power(s)] of the signal that, if it were 

%   present in the bin(s) with the lagest power P_max would, with 

%   probability [1-prob], have produced a power greater than P_max in that

%   bin" (Vaughan+1994, Sec. 3.3.). It works by inverting numerically the

%   noncentral chi^2 CDF which describes the probability of measuring P_max

%   based on the number of power spectra and underlying signal power.

%   Inputs:

%   - max_power (scalar or array): measured power(s) to which to associate

%     the underlying signal power upper limit(s) 

%   - n_ps (scalar): number of power spectra summed to obtain max_power

%   - prob (scalar): probability that a signal with P_ul would result in a

%     power <= P_max 

%   Output: 

%   - p_ul (scalar or array): signal power upper limit(s)

    parser=inputParser;

    parser.KeepUnmatched=true;

    % Simple check on 0<prob<=1 

    addRequired(parser,'prob',@(x) isreal(x) && x<=1 && x>0);

    %p_ul = fzero(@(x) ncx2cdf(max_power,2*n_ps,x)-prob,max_power);

    p_ul = arrayfun(@(max_power) fzero(@(x) ncx2cdf(max_power,2*n_ps,x)-prob,max_power),max_power);

end

scsearch.m

+116 −130
Original line number Diff line number Diff line
@@ -2,20 +2,32 @@ clear; 
close all;



import matlab.io.*

addpath('functions/')



maxNumCompThreads(48);

maxNumCompThreads(24);



% pathfi = '/data/Sorgenti/3FGLJ1544.6-1125/HST/';

pathfi = '/data/Sorgenti/SCORPIUS-X-1/';

% pathfi = '/data/Sorgenti/J1023/';

folname = [pathfi,'SP_search_',char(datetime('now','Format','dd_MM'))];

figfolname = [folname,'/figures'];

if (exist(folname,"dir")==7)

    disp(['Folder ',folname,' already exists, previous Lambda files (if present) may be lost']);

    if (exist(figfolname,"dir")==7)

        disp(['Folder ',figfolname,' already exists, previous graphs (if present) may be lost']);

    else

        mkdir(figfolname);

        disp(['Folder ',figfolname,' created']);

    end

else

    mkdir(folname);

    disp(['Folder ',folname,' created']);

    mkdir(figfolname);

    disp(['Folder ',figfolname,' created']);

end

pathfi = [folname,'/'];

pathfigu = [figfolname,'/'];





reset(gpuDevice(1));

gpuDevice(1);
@@ -34,17 +46,17 @@ tic 
finame = 'ScoX1_20220628_N2_Targ_Bary.fits'

% finame = 'ofba01010_tag_bary_lambda_165_310_chan994_1006.fits';



filoc = [pathfi,'../',finame];

file=fits.openFile(filoc);

% % filoc = [pathfi,'../',finame];

% % file=fits.openFile(filoc);



%%%%%% SiFAP CASE



mjdref=fits.readKey(file,'MJDREF');

MJDREF=str2double(mjdref);

fits.movAbsHDU(file,2);

tstart = str2double(fits.readKey(file,'TSTART'));

t_raw = fitsread(filoc,"binarytable");

t_raw = t_raw{1};

% % mjdref=fits.readKey(file,'MJDREF');

% % MJDREF=str2double(mjdref);

% % fits.movAbsHDU(file,2);

% % tstart = str2double(fits.readKey(file,'TSTART'));

% % t_raw = fitsread(filoc,"binarytable");

% % t_raw = t_raw{1};



%%%%%%% HST CASE
@@ -69,7 +81,7 @@ t_raw = t_raw{1}; 


%%%%%%%



fits.closeFile(file);

% % fits.closeFile(file);





% info = fitsinfo(finame).BinaryTable.Keywords;
@@ -84,24 +96,30 @@ toc 


pathfi = [folname,'/'];



%NOTE: t0 used to define the orbital phase γ is equal to the mid point

%of the observation span for dataset; M2015, Sec.7 FAAAAAAAAAAAAAAAAALSE

% FAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAALSE

t0 = (t_raw(end)+t_raw(1))/2;

%Rebin -------------------------------------------------------------------

Nyq = 2000 %cambiare in caso

Tseg = 256



dt_psd=1/(2*Nyq); %risoluzione spettro di potenza, diverso dal dt dei segmenti

nbin=round((t_raw(end)-t_raw(1))/dt_psd);

% nbin=round((t_raw(end)-t_raw(1))/dt_psd);



countsname = sprintf('counts_check_Tseg_%d.mat',Tseg);

%%

if((Nyq == 2000) & (exist([pathfi,countsname],"file")==2))

    load([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw");

    load([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw","mjdref","MJDREF","tstart","nbin");

    % clear t_raw

else

    tic

    filoc = [pathfi,'../',finame];

    file=fits.openFile(filoc);

    mjdref=fits.readKey(file,'MJDREF');

    MJDREF=str2double(mjdref);

    fits.movAbsHDU(file,2);

    tstart = str2double(fits.readKey(file,'TSTART'));

    t_raw = fitsread(filoc,"binarytable");

    t_raw = t_raw{1};

    fits.closeFile(file);

    nbin=round((t_raw(end)-t_raw(1))/dt_psd);

    [C,t]=(histcounts(t_raw,nbin));

    C=C.'; %conteggi

    t=(t(1:end-1)+(t(2)-t(1))/2).'; % vettore tempi rebinnato, prendo il centro del bin
@@ -144,7 +162,7 @@ else 
    x = reshape(aux2,[N,M]); clear aux2

    toc

    clear C

    save([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw","-v7.3","-nocompression");

    save([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw","mjdref","MJDREF","tstart","nbin","-v7.3","-nocompression");

end



%%
@@ -209,8 +227,7 @@ porb_gr = (Porb).'; porb_max = Porb; porb_min = Porb; % s 
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));

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

Tasc_propagato = Tasc_old+Norb*(Porb/86400.0)

% porb_unc =  0.0000036*86400.0;

porb_unc = 0.017*2;
@@ -286,36 +303,13 @@ nodesnum = length(porb_gr)*length(a_gr)*length(tasc_gr)*length(f_gr) 
Omega_max = 2*pi/porb_min;

Omega_min = 2*pi/porb_max;



% tdiff = tmid - [tasc_max,tasc_min];

% gam = tdiff(:).*[Omega_min,Omega_max];

% 

% % gamma_min = Omega_min*(t0-tasc_max);

% % gamma_max = Omega_max*(t0-tasc_min);

% 

% % Find largest singam

% % M2015 eq. 15

% singal = zeros(4,1);

% singah = zeros(4,1);

nismin = zeros(4,1);

nismax = zeros(4,1);



% Find the range in ν^s (i.e. the maximum span of Eq. (15) covered by 

% combining the search parameters over their respective ranges)

% Sin(gamma) going from -1 to 1

% We take sin(gamma) going from -1 to 1

% Adding special case for s = 1

% s = 1;

% if a_max*Omega_max>1 

%     nismin(s) = -f_max*a_max*(Omega_max) + f_max;

%     if a_min*Omega_min>1

%         nismax(s) = -f_min*a_min*(Omega_min) + f_min;

%     else

%         nismax(s) = -f_max*a_min*(Omega_min) + f_max;

%     end

% else

%     nismin(s) = -f_min*a_max*(Omega_max) + f_min;

%     nismax(s) = -f_max*a_min*(Omega_min) + f_max;

% end

% Is (1-Ome*a) more than fmin-fmax? to ask myself

s = 1;

nismin(s) = (1 - a_max*Omega_max);

nismax(s) = (1 + a_max*Omega_max);
@@ -387,11 +381,11 @@ tic 
% vari ni1,ni2,...,nis_s

%%%[kung,fu,fight] = ndgrid(nis{1},nino{2},nino{3});

tm=gpuArray(tm);

nibank = gpuArray(combinations(nis{:}).Variables);

nifax = nibank(:,1:s_s).*infax(1:s_s);

nibank = combinations(nis{:}).Variables;

nifax = gpuArray(nibank(:,1:s_s).*infax(1:s_s));

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

% stores array elements in a column-major layout

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

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

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

f_ind = zeros(1,length(f_gr),'uint32');

for n = 1:length(f_gr)
@@ -451,7 +445,13 @@ for m=1:M 
    tedge(N+1) = tm(m,N) + 0.5*dt -tmid(m);

    % tedge = (tm(m,1)-1.5*dt-tmid)

    Y0 = fft(xtemp);

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



    graph1 = [pathfigu,'PS_',finame,'_Tseg',num2str(Tseg),'_m',num2str(m),'.pdf'];

    graph2 = [pathfigu,'CT_',finame,'_Tseg',num2str(Tseg),'_m',num2str(m),'.pdf'];

    if (~(isfile(graph1)) || ~(isfile(graph2)))

        [VsuEtot,epsct,hmp] = powdist(Nyq,ttemp,false,graph1,graph2);

    end

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

    % norman = 2./(Nphot.*var(xtemp)./mean(xtemp));

    tau = ones(size(ttemp),"like",ttemp);
@@ -478,14 +478,14 @@ for m=1:M 


        % Y0 = fft(xtemp);

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

        Y1 = fft(mat).';

        Y1 = abs(fft(mat).');

        % 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) = Y1(f_ind(:)); 

        Lambda(:,i) = gather(Y1(f_ind(:))./sqrt(Y1(1))); 

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

        % 

        % for n = 1:length(f_gr)
@@ -496,7 +496,7 @@ for m=1:M 
        % end

        % disp('1ni')

    end

    Lambda = (abs(Lambda).^2).*norman;

    Lambda = (Lambda.^2).*norman;

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

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

    % toc
@@ -559,11 +559,12 @@ else 
    disp(['The length of nibank in greater than ',num2str(intmax('uint64')),', something definitely makes no sense here'])

end



disp(['length(parbank) = ',num2str(length(parbank)),'']);

nimask = false(length(nibank),1);

% totlam = zeros(length(f_gr),length(parbank),'single');

totlam = ones(length(f_gr),length(parbank),'single');

guess = length(parbank)*M*16192.4539/(4.0*7546032.0);

disp(['length(parbank) = ',num2str(parno),'']);

% nimask = false(length(nibank),1);

% Writing on matrix of ones is faster in MATLAB than on zeroes, we subtract

% one after all calculations are done

totlam = ones(length(f_gr),parno,'single');

guess = parno*M*16192.4539/(4.0*7546032.0);

disp(['(Probably terrible) Estimate of the time needed for the second loop = ',num2str(ceil(guess)),' s']);  disp(' ');



bles = cospi(1/2);
@@ -572,10 +573,7 @@ for m = 1:M 
    segstart = tic;

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

    load([pathfi,lamfiname]);

    % toc

    % disp(m);

    % tic

    for i=1:length(parbank)

    for i=1:parno

        

        curpar = [bles,parbank(i,:)];

        % Move from ν,P,a,T_asc to ν,Ω,a,γ
@@ -588,70 +586,43 @@ for m = 1:M 
        curni = -curni.*curpar(3);

        curni(1) = curni(1) + 1;



        [Idx,D] = knnsearch(nisearcher,curni);

        % if (~nimask(Idx))

        nimask(Idx) = 1;

        % end

        % [Idx,D] = knnsearch(nisearcher,curni);

        [schi3,fo3] = min(abs(nibank(:,:)-curni),[],1);

        Idx = sum(fo3)-length(fo3)+1;



        % nimask(Idx) = 1;



        niwalk(m,i) = Idx;



        totlam(:,i) = totlam(:,i)+Lambda(:,Idx);



    end

    % toc 

    segend = toc(segstart);

    avetime = avetime + segend;

    disp(['Just finished segment num.: ',num2str(m),'.']);

    disp(['Elapsed time since loop began = ',num2str(avetime),' s, estimated time left = ',num2str(((M-m)*avetime/m)),' s.']); disp(' ');

    clear Lambda

end

% % for n=1:length(f_gr)

% %     for i=1:length(parbank)

% %         curpar = [f_gr(n),parbank(i,:)];

% %         % Move from ν,P,a,T_asc to ν,Ω,a,γ

% %         curpar(2) = 2*pi/curpar(2);

% %         for m = 1:M

% %             curpar(4) = curpar(2)*(tmid(m) - parbank(i,3));

% %             curni=zeros(1,s_s);

% %             for s=1:s_s

% %                 curni(s) = (curpar(2)^s)*sin(curpar(4)+0.5*s*pi);

% %             end

% %             curni = -curni.*curpar(1).*curpar(3);

% %             curni(1) = curni(1) + curpar(1);

% % 

% %             lamfiname = sprintf('Lambda_seg_%d.mat', m);

% %             % load(lamfiname);

% %             load([pathfi,lamfiname]);

% % 

% %             [Idx,D] = knnsearch(nisearcher,curni);

% %             totlam(i,n) = totlam(i,n)+Lambda(Idx,n);

% % 

% %             clear Lambda

% % 

% %         end

% %         % % if (totlam > bestpar(1))

% %         % %     bestpar = [totlam,f_gr(n),parbank(i,:)];

% %         % % end

% %     end

% % end

% toc

%% 



% totlam was initialized as ones

totlam = totlam - 1;



%% Calculate detection threshold for a multi-trial false alarm prob. = pfa

pfa = 0.01;

sigmastar = 2*gammaincinv((1-(1-pfa)^(1/(length(parbank)*length(f_gr)))),M,'upper');

sigmastar = 2*gammaincinv((1-(1-pfa)^(1/(parno*length(f_gr)))),M,'upper');

%sigmastar = 2*gammaincinv((1-(1-pfa)^(1/(sum(nimask)*length(f_gr)))),M,'upper');

[bru,tto] = max(totlam,[],2);

[maxtotlam,maxinde] = max(bru);

bestpar = [f_gr(maxinde), parbank(tto(maxinde),:), double(maxtotlam)]; 

bparch = string(bestpar);

tascmjd = MJDREF+tasc_gr./86400;

pfamax = 1.0-( (1.0-gammainc(0.5*double(maxtotlam),M,'upper'))^(length(parbank)*length(f_gr)));

pfamax = 1.0-( (1.0-gammainc(0.5*double(maxtotlam),M,'upper'))^(parno*length(f_gr)));

disp('Our most significant peak has parameters:');

disp(['f_spin = '+bparch(1)+' Hz, p_orb = '+bparch(2)+' s, a*sin(i)/c = '+bparch(3)+' s, t_asc = '+bparch(4)+' s.']);

disp(['t_asc in days is = ',num2str((MJDREF+bestpar(4)/86400)),'.']);

disp(['Its detection statistics is Σ_max = ',num2str(maxtotlam),',']);

disp(['which considering all parameters, corresponds to a false alarm probability of ',num2str(pfamax),','])

disp(['which considering all parameters (',num2str(parno),') corresponds to a false alarm probability of ',num2str(pfamax),','])

disp(['while the multi-trial threshold for a false alarm probability = ',num2str(pfa),' is eq. to ',num2str(sigmastar),'.']);

truwalk = size(unique(niwalk.', 'rows'),1);

sigmastar2 = 2*gammaincinv((1-(1-pfa)^(1/(truwalk*length(f_gr)))),M,'upper');
@@ -691,8 +662,13 @@ diary off 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%% Fai il grafico delle totlam a posteriori %%%%%%%%%%%%%%%%



figu = figure('visible','on');

% figu = figure('visible','on');

rina = [pathfigu,'f_maxtotlam.m'];

plot(f_gr,bru)

yline(sigmastar,'--','Σ*')

yline(sigmastar2,'--','Σ*2')

saveas(figu,rina);



% for n = 1:length(f_gr)

%     hold on
@@ -707,41 +683,51 @@ yline(sigmastar,'--','Σ*') 
% contour(a_gr,tasc_gr,squeeze(max(provaccia,[],3)),[20:20:120])

% contour(a_gr,porb_gr,squeeze(max(provaccia,[],1)).',[20:20:120])



% provaccia = reshape(totlam.',length(parbank)*(length(f_gr)),1);

% provaccia = reshape(provaccia,length(tasc_gr),length(a_gr),length(porb_gr),length(f_gr));

% figstep = (10);

if (inGB<350)

    provaccia = reshape(totlam.',length(parbank)*(length(f_gr)),1);

    provaccia = reshape(provaccia,length(tasc_gr),length(a_gr),length(porb_gr),length(f_gr));

    minfiggrid = min(maxtotlam*0.9,sigmastar2*0.9);

    maxfiggrid = max(maxtotlam*1.1,sigmastar2*1.1);

    figstep = (maxfiggrid-minfiggrid)/4;

    fig_grid = [minfiggrid:figstep:maxfiggrid];

    figu = figure('visible','on');

    % figu = figure('visible','on');

    rina = [pathfigu,'porb_tasc.m'];

    contourf(porb_gr,tascmjd,squeeze(max(provaccia,[],[2 4])),fig_grid,'ShowText','on')

    hold on

    scatter(bestpar(2),(MJDREF+bestpar(4)/86400))

    if (exist('trupar') == 1)

        scatter(trupar(2),trupar(4),'+')

    end

    hold off

    saveas(figu,rina);

    % figu = figure('visible','off');

% rina = "porb_tasc.m";

% contourf(porb_gr,tascmjd,squeeze(max(provaccia,[],[2 4])),[sigmastar*0.9:figstep:1.1*maxtotlam],'ShowText','on')

% hold on

% scatter(bestpar(2),(MJDREF+bestpar(4)/86400))

% if (exist('trupar') == 1)

%     scatter(trupar(2),trupar(4),'+')

% end

% hold off

% saveas(figu,rina);

    figu = figure('visible','on');

    rina = [pathfigu,'porb_a.m'];

    contourf(porb_gr,a_gr,squeeze(max(provaccia,[],[1 4])),fig_grid,'ShowText','on')

    hold on

    scatter(bestpar(2),bestpar(3))

    if (exist('trupar') == 1)

        scatter(trupar(2),trupar(3),'+')

    end

    hold off

    saveas(figu,rina);

    % figu = figure('visible','off');

% rina = "porb_a.m";

% contourf(porb_gr,a_gr,squeeze(max(provaccia,[],[1 4])),[sigmastar:figstep:1.1*maxtotlam],'ShowText','on')

% hold on

% scatter(bestpar(2),bestpar(3))

% if (exist('trupar') == 1)

%     scatter(trupar(2),trupar(3),'+')

% end

% hold off

% saveas(figu,rina);

% figu = figure('visible','off');

% rina = "f_porb.m";

% contourf(f_gr,porb_gr,squeeze(max(provaccia,[],[1 2])),[sigmastar:figstep:1.1*maxtotlam],'ShowText','on')

% hold on

% scatter(bestpar(1),bestpar(2))

% if (exist('trupar') == 1)

    % scatter(trupar(1),trupar(2),'+')

% end

% hold off

% saveas(figu,rina);

    figu = figure('visible','on');

    rina = [pathfigu,'f_porb.m'];

    contourf(f_gr,porb_gr,squeeze(max(provaccia,[],[1 2])),fig_grid,'ShowText','off')

    hold on

    scatter(bestpar(1),bestpar(2))

    if (exist('trupar') == 1)

        scatter(trupar(1),trupar(2),'+')

    end

    hold off

    saveas(figu,rina);

end



% if (ceil(timerapp)<11)

% exit

% end



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

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