% matlab script for computing the upper-limit on the signal amplitude h0
% from a radiometer search result
%
%
clear all;
close all;

% define observation params
Tobs = 1*86400;          % define the complete observation length (converges after >>1 days, but result is identical for integer days)
T = 4;                   % define the coherent length
Nseg = Tobs/T;           % define the number of segments
Nsky = 200;               % number of declination values to try
conf = 0.9;              % define the UL confidence

% define Sco X-1 sky location
alpha_sco = 4.2757;   % sco X-1
delta_sco = -0.2730;  % sco X-1

% vector of delination values
%skyvec(1,:) = [alpha_sco,delta_sco];
skyvec(:,1) = 2*pi*rand(1,Nsky);
skyvec(:,2) = asin(-1 + 2*rand(1,Nsky));

% define coarse time for computing antenna pattern integrals
dtc = 120;
tc = dtc/2:dtc:Tobs-dtc/2;

% define frequency vector
df = 1/T;
fmax = 2048;
f = 40:df:fmax-df;
Nf = length(f);

% define detector params from JKS 
Wr = 2*pi/86400;                      % earth angular rotation frequency 
phir0 = 1.0;                           % some initial earth rotation phase 
Ndet = 2;                             % number of detectors
lambda(1) = pi*(46.45/180);           % H1
gamma(1) = pi*(171.8/180);            % H1
phir(1) = phir0 + pi*(119.41/180);    % H1
lambda(2) = pi*(30.56/180);           % L1
gamma(2) = pi*(243.0/180);            % L1
phir(2) = phir0 + pi*(90.77/180);     % L1
lambda(3) = pi*(43.63/180);           % V1
gamma(3) = pi*(116.5/180);            % V1
phir(3) = phir0 + pi*(-10.5/180);     % V1

% define vectors of cosi and psi (we will integrate over these)
Ngrid = 100;                     % the number of grid points in each dimension
mustar = zeros(Ngrid,Ngrid);     
deta = 2/Ngrid;                  % the spacing on eta = cosi
etavec = -1:deta:1-deta;         % the vector of eta values ranging -1 -> 1
dpsi = (pi/2)/Ngrid;             % the spacing on psi
psivec = -pi/4:dpsi:pi/4-dpsi;   % the vector of psi values ranging -pi/4 -> pi/4

% define similar values for circularly polarised case cosi = 1
mustar_circ = zeros(1,Ngrid);
etac = 1;

% define alpha vector - alpha = h0*T/(P1*P2)^1/4
% we will integrate over this vector later to get Bayesian posteriors on h0
alphamax = 5;
Nalpha = 250;
alphavec = linspace(0,alphamax,Nalpha); 
dalpha = alphavec(2)-alphavec(1);

% loop over declination
for q=1:Nsky

    % current sky value
    alpha = skyvec(q,1);
    delta = skyvec(q,2);
    
    % initialise antenna functions
    at = zeros(Ndet,length(tc));
    bt = zeros(Ndet,length(tc));
    A = 0;
    B = 0;
    C = 0;
    
    % loop over detectors
    for j=1:Ndet

        % generate antenna pattern functions from JKS on coarse time vector
        % these are clowly varying functions so we don't need fine time
        % resolution
        at(j,:) = (1/16)*sin(2*gamma(j))*(3-cos(2*lambda(j)))*(3-cos(2*delta))*cos(2*(alpha-phir(j)-Wr*tc)) ...
            - (1/4)*cos(2*gamma(j))*sin(lambda(j))*(3-cos(2*delta))*sin(2*(alpha-phir(j)-Wr*tc)) ...
            + (1/4)*sin(2*gamma(j))*sin(2*lambda(j))*sin(2*delta)*cos(alpha-phir(j)-Wr*tc) ...
            - (1/2)*cos(2*gamma(j))*cos(lambda(j))*sin(2*delta)*sin(alpha-phir(j)-Wr*tc) ...
            + (3/4)*sin(2*gamma(j))*(cos(lambda(j))*cos(delta))^2;

        bt(j,:) = cos(2*gamma(j))*sin(lambda(j))*sin(delta)*cos(2*(alpha-phir(j)-Wr*tc)) ...
            + (1/4)*sin(2*gamma(j))*(3-cos(2*lambda(j)))*sin(delta)*sin(2*(alpha-phir(j)-Wr*tc)) ...
            + cos(2*gamma(j))*cos(lambda(j))*cos(delta)*cos(alpha-phir(j)-Wr*tc) ...
            + (1/2)*sin(2*gamma(j))*sin(2*lambda(j))*cos(delta)*sin(alpha-phir(j)-Wr*tc);

        % compute integrals as in JKS
        A = A + (1/Tobs)*sum(at(j,:).*at(j,:))*dtc;
        B = B + (1/Tobs)*sum(bt(j,:).*bt(j,:))*dtc;
        C = C + (1/Tobs)*sum(at(j,:).*bt(j,:))*dtc;

    end

    % compute expected variance of the cross correlation statistic
    % this multiplied by P1*P2/T^2 is the variance of the cross-corr stat
    % Eq 18.  It is the same for circ and non-circ
    sigsq = 2./sum((at(1,:).*at(2,:) + bt(1,:).*bt(2,:)).^2);

    % compute expected mean SNR value - here we average over eta and psi
    % loop over eta
    for j=1:Ngrid

        % compute the Fplus and Fcross for each detector as function of
        % time for each psi value
        Fp1 = at(1,:)'*cos(2*psivec) + bt(1,:)'*sin(2*psivec);
        Fp2 = at(2,:)'*cos(2*psivec) + bt(2,:)'*sin(2*psivec);
        Fc1 = bt(1,:)'*cos(2*psivec) - at(1,:)'*sin(2*psivec);
        Fc2 = bt(2,:)'*cos(2*psivec) - at(2,:)'*sin(2*psivec);

        % compute the components of Eq 17.  
        X = 0.25*((1+etavec(j)^2)^2)*Fp1.*Fp2 + Fc1.*Fc2*etavec(j)^2;
        Y = Fp1.*Fp2 + Fc1.*Fc2;
        V = sum(Y.^2);

        % the following when multiplied by h0^2 is the mean of the cross-corr
        % statistic Eq 17.  mustar is a 2D grid of values for each eta and
        % psi combination
        mustar(:,j) = (sum(X.*Y)./V);

    end

    % compute expected mean SNR value for circularly polarised signals   
    % here we simply fix eta = 1 and vary psi
    Fp1c = at(1,:)'*cos(2*psivec(1)) + bt(1,:)'*sin(2*psivec(1));
    Fp2c = at(2,:)'*cos(2*psivec(1)) + bt(2,:)'*sin(2*psivec(1));
    Fc1c = bt(1,:)'*cos(2*psivec(1)) - at(1,:)'*sin(2*psivec(1));
    Fc2c = bt(2,:)'*cos(2*psivec(1)) - at(2,:)'*sin(2*psivec(1));

    % compute the components of Eq 17 for circularly polarised signals
    Xc = 0.25*((1+etac^2)^2)*Fp1c.*Fp2c + Fc1c.*Fc2c*etac^2;
    Yc = Fp1c.*Fp2c + Fc1c.*Fc2c;
    Vc = sum(Yc.^2);

    % the following when multiplied by h0^2 is the mean of the cross-corr
    % statistic
    mustar_circ = (sum(Xc.*Yc)./Vc);

    % compute pdfs on h0

    % reshape mustar from a 2D square into a 1D vector for circ and
    % non-circ
    newmustar = reshape(mustar,1,[]);
    newmustar_circ = reshape(mustar_circ,1,[]);

    % initialise marginal likelihoods
    palpha = zeros(1,Nalpha);
    palpha_circ = zeros(1,Nalpha);
    
    % loop over h0 (alpha is propto h0) and compute the posterior on h0 for
    % each case
    for j=1:Nalpha

        % define true mean of the cross-corr statistic
        % the variance is independent of h0 (in the weak signal case)
        mu = newmustar*alphavec(j)^2;

        % compute likelihood for non-circular polarisation
        % here we simply intgerate the Gaussian likelihood over the
        % nuisance parameters eta and psi with flat priors
        % we ignore normalisation because all will get normalised later
        palpha(j) = sum(exp(-0.5*(mu.^2)./sigsq));

        % define mean cross-corr statistic for the circular polarisation
        % case
        mu_circ = newmustar_circ*alphavec(j)^2;

        % compute marginal likelihood for the circular polarisation case
        palpha_circ(j) = sum(exp(-0.5*(mu_circ.^2)./sigsq));

    end

    % normalise the posteriors so that the integrate to 1 
    palpha = palpha./(sum(palpha)*dalpha);
    palpha_circ = palpha_circ./(sum(palpha_circ)*dalpha);

    % compute 90% confidence UL on alpha
    [temp s] = unique(cumsum(palpha)*dalpha);
    alpha90 = interp1(temp,alphavec(s),conf);

    % compute 90% confidence UL on alpha
    [tempc sc] = unique(cumsum(palpha_circ)*dalpha);
    alphacirc90 = interp1(tempc,alphavec(sc),conf);

    % record the result as a function of declination
    result(q) = alpha90/alphacirc90;
    fprintf(1,'done %d/%d (%f %f) -> %f\n',q,Nsky,skyvec(q,1),skyvec(q,2),result(q));
    
end

