Commit 5cad424f authored by Riccardo La Placa's avatar Riccardo La Placa
Browse files

add multiobs to be sure

parent 89846488

scsearch_1023_multiobs.m

0 → 100644
+270 −0
Original line number Diff line number Diff line 
clear;

close all;



import matlab.io.*

addpath('functions/')



maxNumCompThreads(48);



pathfi = '/data/Sorgenti/J1023/30-31_01_2020/';

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



diary([pathfi,'loggerino_metrica_',char(datetime('now','Format','dd_MM')),'.log'])

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

daterella = ['Date and time now ',char(datetime('now')),' (UTC), maybe'];

disp(daterella); clear daterella



tic



mario = split(ls('*.fits'));

mario = sort(mario(~cellfun('isempty',mario)));

numfiles = length(mario);



toc



pathfi = [folname,'/'];



Nyq = 2000 %cambiare in caso

Tseg = 256

dt_psd=1/(2*Nyq);

dt = dt_psd;

N = fix(Tseg/dt);

f_min = 550

f_max = 600

porb_min = 17115.3216 % s

porb_max = 17115.3216; porb_gr = (porb_max).'; % s

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

a_min = a_fact*256000*0.04

a_max = a_fact*280000*0.72

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

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

f_step = 1/Tseg

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

% tasc_mid = (tasc_max + tasc_min)/2

% tasc_step = tasc_max-tasc_min;

% porb_step = porb_max;

% a_step = a_max;

tol = 0.9

MJDREF = zeros(numfiles,1);



tic

%% First step in par., def. also max min tasc

finame = char(mario(1))

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

file=fits.openFile(filoc);

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

MJDREF(1)=str2double(mjdref);

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

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

tasc_mid = (tasc_max + tasc_min)/2

fits.movAbsHDU(file,2);

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

t_raw = fitsread(filoc,"binarytable");

t_raw = t_raw{1};

fits.closeFile(file);

nbin=ceil((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

M=fix((t_raw(end)-t_raw(1))/Tseg);

tm = zeros(M,N);

tmid = zeros(M,1);

for m = 1:M

    tm(m,1) = t(1)+(m-1)*Tseg; %t(j) è già centrato

    for j = 2:N

        tm(m,j)= tm(m,1)+(j-1)*dt;

    end

    tmid(m) = (tm(m,N)+tm(m,1))/2;

end



% Time series for each segment

x = zeros(N,M);

aux1 = M*N;

aux2 = C(1:aux1); clear aux1

x = reshape(aux2,[N,M]); clear aux2

clear C

countsname = ['counts_',finame,'_Tseg_',num2str(Tseg),'.mat'];

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



toff = tstart - tasc_mid;

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

%%

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

% 1.0/(2.0*f_max*tol);

%%

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

a_step = stepparino

%%

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

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

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

%%

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

tasc_step = stepparino



%% Multiple obs. (that's what w'ere here for)

if (numfiles>1)

for i = 2:numfiles

    finame = char(mario(i))

    %% WILL have to choose wether to do only the parameter space check here

    

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

    file=fits.openFile(filoc);

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

    MJDREF(i)=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=ceil((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

    M=fix((t_raw(end)-t_raw(1))/Tseg);

    tm = zeros(M,N);

    tmid = zeros(M,1);

    for m = 1:M

        tm(m,1) = t(1)+(m-1)*Tseg; %t(j) è già centrato

        for j = 2:N

            tm(m,j)= tm(m,1)+(j-1)*dt;

        end

        tmid(m) = (tm(m,N)+tm(m,1))/2;

    end



    % Time series for each segment

    x = zeros(N,M);

    aux1 = M*N;

    aux2 = C(1:aux1); clear aux1

    x = reshape(aux2,[N,M]); clear aux2

    clear C

    countsname = ['counts_',finame,'_Tseg_',num2str(Tseg),'.mat'];

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



    toff = tstart - tasc_mid;

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

    %%

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

    % 1.0/(2.0*f_max*tol)

    %%

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

    if (stepparino < a_step)

        a_step = stepparino

    end

    % % 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 = deltasc(a_max,porb_min,f_max,toff,Tseg,M,tol);

    %%

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

    if (stepparino < tasc_step)

        tasc_step = stepparino

    end

end

end

%% Finally make grid of orb. parameters with smallest steps (and define trupars)

a_gr = (a_min:a_step:a_max).';

tasc_gr = (tasc_min:tasc_step:tasc_max).';

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

trupar = [5.924214682724856e+02,1.71155216592e+04,0.343356,58878.9911009];

toc



%% Just getting an idea

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



%% Aux. param. and coherent grid generation

% Aux par

Omega_max = 2*pi/porb_min;

Omega_min = 2*pi/porb_max;



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)

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

% Adding special case for s = 1

s = 1;

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

nismax(s) = (1 + a_max*Omega_max);

for s = 2:4

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

    nismax(s) =  a_max*(Omega_max^s);

end



% Try s* and check ν^s range

tic

s_s = 4;

while(1)

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

        s_s=s_s-1;

    else

        break;

    end

end



%%  Create a coherent template bank in the ν_s space

% Nni+1 is the number of ni_s points in the grid (counting borders, 

% Nni is the number of steps) for each s couple.

% We do NOT use eq. 23 M2015, but instead use delta_ni_new so that our grid

% never exceeds half a freq. step in observed freq. from each contr.

% nis is a structure which holds the various vectors of ni_s^(m) (our

% coherent template bank)

% each nis{s} is the array giving the respective grid, and its i-th 

% element is accessed as nis{s}(i)

Nni = zeros(s_s,1);

for s = 1:s_s

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

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

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

    end

    franco=nismin(s):((nismax(s)-nismin(s))/(Nni(s))):nismax(s);

    nis{s}=franco;

end

infax = gpuArray(1./factorial(1:s_s)); % MATLAB's factorial is slow



toc



%% ADD CYCLE ON OBS.



tm=gpuArray(tm);

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

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

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

for n = 1:length(f_gr)

    [~,f_ind(n)] = min(abs(F1-f_gr(n)));

end
 
 No newline at end of file