Commit c20554ab authored by Riccardo La Placa's avatar Riccardo La Placa
Browse files

FINISH CHANGES TO SCO STUPIDO

parent 2515aec9

scsearch_SCO_multiobs_multifreq.m

+34 −12
Original line number Original line Diff line number Diff line
@@ -4,7 +4,7 @@ close all; 
import matlab.io.*

import matlab.io.*

addpath('functions/')

addpath('functions/')





maxNumCompThreads(72);

maxNumCompThreads(96);





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

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

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

pathfi = '/data/Sorgenti/SCORPIUS-X-1/2022-2023/';
@@ -65,9 +65,9 @@ porb_min = Porb - 3*porb_unc; % s 
porb_max = Porb + 3*porb_unc; % s

porb_max = Porb + 3*porb_unc; % s

porb_gr = (Porb).';

porb_gr = (Porb).';

Tasc_old = 56722.6303715306713

Tasc_old = 56722.6303715306713

Tasc_old_unc = 3*(33 + 4); % s (3sigma Killestein + 4 s GPS def. uncertainty)

Tasc_old_unc = 3*(33); % s (3sigma Killestein + 4 s GPS def. uncertainty)

a_fact = porb_max/(2.0*pi*299792458);

a_fact = porb_max/(2.0*pi*299792458);

a_min = 1.44

a_min = 1.45

a_max = 3.25

a_max = 3.25

% tasc_min = (58878.9911009 - MJDREF).*86400 - porb_max/2

% tasc_min = (58878.9911009 - MJDREF).*86400 - porb_max/2

% tasc_max = (58878.9911009 - MJDREF).*86400 + porb_max/2

% tasc_max = (58878.9911009 - MJDREF).*86400 + porb_max/2
@@ -136,16 +136,17 @@ for fifo = 1:n_fr_intervals 
        save([pathfi,countsname],"x","M","N","tm","tmid","dt","t_raw","tstart","nbin","mjdref","-v7.3","-nocompression");

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

    end

    end

    

    

    Nphotot(1) = sum(x,"all");

    Mtot = M;

    Mtot = M;

    MJDREF(1)=str2double(mjdref);

    MJDREF(1)=str2double(mjdref);

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

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

    Tasc_propagato = Tasc_old+Norb*(Porb/86400.0)

    Tasc_propagato = Tasc_old+Norb*(Porb/86400.0)

    % porb_unc =  0.0000036*86400.0;

    % porb_unc =  0.0000036*86400.0;

    % Tasc_old_unc = 0.0053*86400.0; % 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_new_unc = sqrt(Tasc_old_unc^2 + (Norb*3*porb_unc)^2); % s 

    tasc_mid = (Tasc_propagato - MJDREF(1))*86400

    tasc_mid = (Tasc_propagato - MJDREF(1))*86400

    tasc_min = tasc_mid - 2.0*Tasc_new_unc

    tasc_min = tasc_mid - Tasc_new_unc

    tasc_max = tasc_mid + 2.0*Tasc_new_unc

    tasc_max = tasc_mid + Tasc_new_unc

    % tasc_min = (57449.72579 - MJDREF(1)).*86400 - porb_max/2

    % tasc_min = (57449.72579 - MJDREF(1)).*86400 - porb_max/2

    % tasc_max = (57449.72579 - MJDREF(1)).*86400 + porb_max/2

    % tasc_max = (57449.72579 - MJDREF(1)).*86400 + porb_max/2

    % tasc_mid = (tasc_max + tasc_min)/2

    % tasc_mid = (tasc_max + tasc_min)/2
@@ -224,18 +225,18 @@ for fifo = 1:n_fr_intervals 
            clear C

            clear C

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

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

        end

        end

        Mtot = Mtot + M;



        MJDREF(o)=str2double(mjdref);

        MJDREF(o)=str2double(mjdref);

        Nphotot(o) = sum(x,"all");

        Mtot = Mtot + M;

        if (round(abs((MJDREF(o)-Tasc_old))/(Porb/86400.0)) > Norb)

        if (round(abs((MJDREF(o)-Tasc_old))/(Porb/86400.0)) > Norb)

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

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

            Tasc_propagato = Tasc_old+Norb*(Porb/86400.0)

            Tasc_propagato = Tasc_old+Norb*(Porb/86400.0)

            % porb_unc =  0.0000036*86400.0;

            % porb_unc =  0.0000036*86400.0;

            % Tasc_old_unc = 0.0053*86400.0; % 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_new_unc = sqrt(Tasc_old_unc^2 + (Norb*3*porb_unc)^2); % s

            tasc_mid = (Tasc_propagato - MJDREF(1))*86400

            tasc_mid = (Tasc_propagato - MJDREF(1))*86400

            tasc_min = tasc_mid - 2.0*Tasc_new_unc

            tasc_min = tasc_mid - Tasc_new_unc

            tasc_max = tasc_mid + 2.0*Tasc_new_unc

            tasc_max = tasc_mid + Tasc_new_unc

        end

        end

        toff = tstart - tasc_mid;

        toff = tstart - tasc_mid;

        stepparino = delporb(a_max,porb_min,f_max,toff,Tseg,M,tol);

        stepparino = delporb(a_max,porb_min,f_max,toff,Tseg,M,tol);
@@ -398,7 +399,7 @@ for fifo = 1:n_fr_intervals 
            graph1 = [pathfigu,'PS_',finame,'_Tseg',num2str(Tseg),'_m',num2str(m),'.pdf'];

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

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

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

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

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

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

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

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

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

            end

            end

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

            norman = gather(2.*abs(Y0(1))./mean((abs(Y0(F1>1000 & F1<Nyq)).^2)));
@@ -592,6 +593,27 @@ for fifo = 1:n_fr_intervals 




end

end





[maxoni,maxindi] = maxk(bruttone,20);

if (maxoni(1) > sigmastar2) 

    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')

    disp(['A signal was found (total summed power ',num2str(maxoni(1)),') at a frequency = ',num2str(totfgr(maxindi(1)),12),' Hz']);

    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')

else

    pup = power_upper_limit(maxoni,Mtot,0.0027);

    if (exist('back','var')==1)

        amplicost = 1.61*(sum(Nphotot(o)/(Tseg*((Nphot/Tseg) + 15000)))*(Nphot*51)^(-0.5);

        Auls = amplicost.*sqrt(pup)./(pi.*sinc(totfgr(maxindi)./2.*Nyq));

        [maxau,maxauind] = max(Auls);

        disp(['No signal was found, our UL on amplitude is ',num2str(maxau),' (p_ul at 3 sig) at a frequency = ',num2str(totfgr(maxindi(maxauind))),' Hz (no background)']);



    else

        amplicost = 1.61*(Nphotot)^(-0.5); %No background for now

        Auls = amplicost.*sqrt(pup)./(pi.*sinc(totfgr(maxindi)./2.*Nyq));

        [maxau,maxauind] = max(Auls);

        disp(['No signal was found, our UL on amplitude is ',num2str(maxau),' (p_ul at 3 sig) at a frequency = ',num2str(totfgr(maxindi(maxauind))),' Hz (no background)']);

    end

end







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

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

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

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