Commit 9a2e2a84 authored by Riccardo La Placa's avatar Riccardo La Placa
Browse files

Using new Delta ni and fixed sign of T_asc, which worked splendidly on J1023

parent 51f68fc9

scsearch.m

+57 −40
Original line number Diff line number Diff line
@@ -2,8 +2,8 @@ clear; 
close all;



% Define variables, M is segm. number

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

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

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

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

reset(gpuDevice(2));

gpuDevice(2);
@@ -13,15 +13,15 @@ disp(daterella); clear daterella 


tic



% finame = 'SiFAP2_20220429_3FGLJ1544_N1_Bary.fits';

finame = 'J1023_B_2017_Bary.fits';

finame = 'SiFAP2_20220429_3FGLJ1544_N1_Bary.fits';

% finame = 'J1023_B_2017_Bary.fits';

% finame = 'EPN_0744840201_bary.fits';

t_raw = fitsread([pathfi,finame],"binarytable");

%t_raw = fitsread('C:\Users\Filippo\Desktop\J1544_2022\SiFAP2_20220429_3FGLJ1544_N1_Bary.fits','binarytable');

t_raw = t_raw{1};



% MJDREF 1544

MJDREF=59698.943055555602768;

% MJDREF=59698.943055555602768;





% info = fitsinfo(finame).BinaryTable.Keywords;
@@ -55,13 +55,14 @@ toc 
% tasc_gr = zeros(5,1);



% Parameters for 1544

% Tseg=490; %segment length in seconds

% f_gr = (100:1/Tseg:800).';

% %f_gr=(772:1/Tseg:774).';

% % porb_gr = (18000:10:21600).';

% porb_gr=(20868.72).';

% a_gr = (0.01:8.0e-4:0.26).';

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

Tseg=490; %segment length in seconds

MJDREF=59698.943055555602768;

f_gr = (100:1/Tseg:800).';

%f_gr=(772:1/Tseg:774).';

% porb_gr = (18000:10:21600).';

porb_gr=(18000:150:21600).';

a_gr = (0.01:4.0e-3:0.2).';

tasc_gr = (((59698.7231:8e-4:59698.9232) - MJDREF).*86400).';

% for j = 1:5

%     %f_gr(j) = f_tru*((j^2)/9);

%     porb_gr(j) = porb_tru*((j^2)/9);
@@ -69,21 +70,21 @@ toc 
%     tasc_gr(j) = tasc_tru*((j^2)/9);

% end



Tseg=256; %segment length in seconds

f_gr = (500:1/Tseg:600).';

porb_gr=(17000:5:17200).';

a_gr = (0.3:8.0e-4:0.4).';

MJDREF = 58107.1236342593;

tasc_gr = (MJDREF - ((58106.9994:8e-4:58107.1095)).*86400).';



f_min = min(f_gr);

f_max = max(f_gr);

porb_min = min(porb_gr);

porb_max = max(porb_gr);

a_min = min(a_gr);

a_max = max(a_gr);

tasc_min = min(tasc_gr);

tasc_max = max(tasc_gr);

% Tseg=256; %segment length in seconds

% f_gr = (500:1/Tseg:600).';

% porb_gr=(17000:5:17200).';

% a_gr = (0.3:8.0e-4:0.4).';

% MJDREF = 58107.1236342593;

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



f_min = min(f_gr)

f_max = max(f_gr)

porb_min = min(porb_gr)

porb_max = max(porb_gr)

a_min = min(a_gr)

a_max = max(a_gr)

tasc_min = min(tasc_gr)

tasc_max = max(tasc_gr)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -205,7 +206,7 @@ mu_s=0.05; %massimo mismatch sulla griglia coerente da scegliere 
%s_s=uint8(4);

s_s = 4;

while(1)

    if((nismax(s_s)-nismin(s_s))<0.5*delta_ni(mu_s,s_s,g_jj(s_s)))

    if((nismax(s_s)-nismin(s_s))<delta_ni_new(s_s,Tseg))

        s_s=s_s-1;

    else

        break;
@@ -224,7 +225,7 @@ end 
% element is accessed as nis{s}(i)

Nni = zeros(s_s,1);

for s = 1:s_s

    Nni(s) = ceil((nismax(s)-nismin(s))/delta_ni(mu_s,s,g_jj(s)));

    Nni(s) = ceil((nismax(s)-nismin(s))/delta_ni_new(s,Tseg));

    if (mod(Nni(s),2)==0)

        Nni(s) = Nni(s)+1;

    end
@@ -249,7 +250,7 @@ tic 
%%%[kung,fu,fight] = ndgrid(nis{1},nino{2},nino{3});

tm=gpuArray(tm);

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

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

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

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

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

for n = 1:length(f_gr)
@@ -265,9 +266,9 @@ toc 
% dtau = zeros(1,N);

% taul = zeros(N+1,1);

taul = zeros(2,1);

disp(length(nibank));

disp(s_s);

disp(M);

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

s_s

M

for m=1:M



    % tic
@@ -314,7 +315,7 @@ for m=1:M 
        % Careful! ni_0 still needs to be defined (it's the spin freq. 

        % currently being considered)

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

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

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

        % 

        % for n = 1:length(f_gr)

        %
@@ -332,7 +333,7 @@ for m=1:M 
    lamfiname = sprintf('Lambda_seg_%d.mat', m);

    % save(lamfiname,"Lambda");

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

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

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

    disp(m);



end
@@ -357,7 +358,7 @@ toc 




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

totlam = zeros(length(parbank),length(f_gr));

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

tic

% for n=1:length(f_gr)

%     for i=1:length(parbank)
@@ -413,7 +414,7 @@ for m = 1:M 
        if (~nimask(Idx))

            nimask(Idx) = 1;

        end

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

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



    end

    toc
@@ -456,12 +457,14 @@ toc 
pfa = 0.01;

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

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

[bru,tto] = max(totlam);

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

[maxtotlam,maxinde] = max(bru);

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

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

bparch = string(bestpar);

tascmjd = MJDREF+tasc_gr./86400;

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(['while the multi-trial threshold for a false alarm probability = ',num2str(pfa),' is eq. to ',num2str(sigmastar),'.']);
@@ -476,13 +479,22 @@ disp(daterella); clear daterella 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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



% plot(f_gr,max(totlam))

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



% for n = 1:length(f_gr)

%     hold on

%     txt = ['f_{gr} = ',num2str(f_gr(n))];

%     plot(totlam(:,n),'DisplayName',txt)

%     plot(totlam(n,:),'DisplayName',txt)

% end

% legend



% provina = [parbank(:,:),totlam(23661,:).'];

% provaccia = reshape(provina(:,4),length(tasc_gr),length(a_gr),length(porb_gr));

% contour(porb_gr,tasc_gr,squeeze(max(provaccia,[],2)),[20:20:120])

% contour(a_gr,tasc_gr,squeeze(max(provaccia,[],3)),[20:20:120])

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



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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -496,3 +508,8 @@ function res = delta_ni(mu_s,s,g) 
    %Compute eq. 23 M2015

    res= 2.*sqrt(mu_s./(s.*g));

end



function res = delta_ni_new(s,Tcoh)

    %Compute eq. 23 M2015

    res = 0.5*(factorial(s-1))/(Tcoh^s);

end