%MyFont = {'FontName', 'Times New Roman', 'FontSize', 16};
%Number of pairs in B3
a = 1:1240;
stream0 = RandStream('mlfg6331_64','Seed',0);
RandStream.setDefaultStream(stream0);

%Shuffling of the array
a1 = a(randperm(length(a)));
%Divide into 10 groups of equal len
tile_combos = mat2cell(a1,1,[155,155,155,155,155,155,155,155]);
load('maps/5901/real_g_filtp3_weight3_gs_dp0000_jack02.mat'); uk = get_ukpervolt('2017','bicep3');
acc= cal_coadd_ac(ac,uk);
%setwinsize(gcf,600,1000)
b = tile_combos{1}
apsopt.coaddrx=1;
ukpv = get_ukpervolt('2017','bicep3');
for i = 1:8
    if i>1
        b = cat(2,b,tile_combos{i});
    end
    b1 = reshape(b,1,[])
    length(b1)
    acc1=coadd_ac_overrx(acc(b1));
    map(i)=make_map(acc1,m);
    it(i)=nansum(map(i).Pwitime(:));
    % the Pw field is generated in add_masks as 1/(Qvar+Uvar)
    % for total weight it might make more sense to use 1/Qvar+1/Uvar
    tw(i)=nansum(1./map(i).Qvar(:)+1./map(i).Uvar(:))
%   t_name = sprintf('Pol wt of %drx',2*i)
%    t_name = sprintf('Pol wt map of %d coadded pairs', length(b1));
%    subplot(4,2,i); plot_map(m,1./map(i).Qvar+1./map(i).Uvar); colorbar;
%    title(t_name);
end
%set(gca, MyFont{:});
%saveas(gcf,'fig1_b3_2017_jack02.png');
%clf;	

%[p,ind] = get_array_info('20170505');
%% account for non-zero drum angle
%p = rotarray(p, -median(p.drumangle(ind.la)));
%% do approximate transformation from local to celestial coordinates;
%% namely, reversal between postive el and postive dec.
%p = rotarray(p, -180);
%% r,theta to x,y coordinates
%[y,x] = reckon(0, 0, p.r, p.theta-90);
%for i=1:8
%	figure(); setwinsize(gcf(), 300, 300);
%	% Insert whatever pair selection you're really doing here.
%	b = tile_combos{1};
%    	if i>1
%        	b = cat(1,b,tile_combos{2:i})
%    	end
%    	b1 = reshape(b,1,[]);
%    	length(b1)
%	
%	mask = false(length(ind.a),1);
%	sel = b1;
%	mask(sel) = true;
%	% Then just plot as scatter points on the plot, masking unkept detectors by
%	% NaN-ing their coordinates.
%	xm = x(ind.a); xm(~mask) = NaN;
%	ym = y(ind.a); ym(~mask) = NaN;
%	scatter(xm, ym, '.');
%	set(gca(), 'XDir', 'reverse');
%        t_name = sprintf('%d pairs', length(b1));
%	map_name = sprintf('%d.png',i);title(t_name);
%	saveas(gcf,map_name);clf;
%end

%Loading the xxx6 realizations of the different subsets of pairs
load('aps/5901/xxx6_g_filtp3_weight3_gs_dp1100_jack02_overrx_155pairs.mat'); aps1=aps;
load('aps/5901/xxx6_g_filtp3_weight3_gs_dp1100_jack02_overrx_310pairs.mat'); aps2=aps;
load('aps/5901/xxx6_g_filtp3_weight3_gs_dp1100_jack02_overrx_465pairs.mat'); aps3=aps;
load('aps/5901/xxx6_g_filtp3_weight3_gs_dp1100_jack02_overrx_620pairs.mat'); aps4=aps;
load('aps/5901/xxx6_g_filtp3_weight3_gs_dp1100_jack02_overrx_775pairs.mat'); aps5=aps;
load('aps/5901/xxx6_g_filtp3_weight3_gs_dp1100_jack02_overrx_930pairs.mat'); aps6=aps;
load('aps/5901/xxx6_g_filtp3_weight3_gs_dp1100_jack02_overrx_1085pairs.mat'); aps7=aps;
load('aps/5901/xxx6_g_filtp3_weight3_gs_dp1100_jack02_overrx_1240pairs.mat'); aps8=aps;

sf=2*pi./(aps.l.^2);

s=1;

i1=aps.l>10&aps.l<50;
i2=aps.l>330&aps.l<400;
i = aps.l>10

setwinsize(gcf,450,450); clf
h1=loglog(aps.l(i),squeeze(repmat(sf(i),[1,1,size(aps1.Cs_l,3)]).*aps1.Cs_l(i,s,:)), 'r-');
hold on
h2=loglog(aps.l(i),squeeze(repmat(sf(i),[1,1,size(aps1.Cs_l,3)]).*aps2.Cs_l(i,s,:)), 'g-');
h3=loglog(aps.l(i),squeeze(repmat(sf(i),[1,1,size(aps1.Cs_l,3)]).*aps3.Cs_l(i,s,:)), 'b-');
h4=loglog(aps.l(i),squeeze(repmat(sf(i),[1,1,size(aps1.Cs_l,3)]).*aps4.Cs_l(i,s,:)), 'c-');
h5=loglog(aps.l(i),squeeze(repmat(sf(i),[1,1,size(aps1.Cs_l,3)]).*aps5.Cs_l(i,s,:)), 'y-');
h6=loglog(aps.l(i),squeeze(repmat(sf(i),[1,1,size(aps1.Cs_l,3)]).*aps6.Cs_l(i,s,:)), 'k-');
h7=loglog(aps.l(i),squeeze(repmat(sf(i),[1,1,size(aps1.Cs_l,3)]).*aps7.Cs_l(i,s,:)), 'r--');
h8=loglog(aps.l(i),squeeze(repmat(sf(i),[1,1,size(aps1.Cs_l,3)]).*aps8.Cs_l(i,s,:)), 'g--');
hold off
grid
%
xlabel('multipole'); ylabel('TT N_l (\muK^2)');
legend([h1(1),h2(1),h3(1),h4(1),h5(1),h6(1),h7(1),h8(1)],{'155pairs','310pairs','465pairs','620pairs','785pairs','930pairs','1085pairs','1240pairs'});
title('First 100 noise realizations');
%set(gca, MyFont{:});
saveas(gcf,'fig2a_B3_2017_TT_jack02.png');clf;


%Taking the average of the 100 N\ell
aps1.Cs_l=mean(aps1.Cs_l,3);
aps2.Cs_l=mean(aps2.Cs_l,3);
aps3.Cs_l=mean(aps3.Cs_l,3);
aps4.Cs_l=mean(aps4.Cs_l,3);
aps5.Cs_l=mean(aps5.Cs_l,3);
aps6.Cs_l=mean(aps6.Cs_l,3);
aps7.Cs_l=mean(aps7.Cs_l,3);
aps8.Cs_l=mean(aps8.Cs_l,3);

setwinsize(gcf,450,450); 
h1=loglog(aps.l(i),sf(i).*aps1.Cs_l(i,s),'r-');
hold on
h2=loglog(aps.l(i),sf(i).*aps2.Cs_l(i,s),'g-');
h3=loglog(aps.l(i),sf(i).*aps3.Cs_l(i,s),'b-');
h4=loglog(aps.l(i),sf(i).*aps4.Cs_l(i,s),'c-');
h5=loglog(aps.l(i),sf(i).*aps5.Cs_l(i,s),'y-');
h6=loglog(aps.l(i),sf(i).*aps6.Cs_l(i,s),'k-');
h7=loglog(aps.l(i),sf(i).*aps7.Cs_l(i,s),'r--');
h8=loglog(aps.l(i),sf(i).*aps8.Cs_l(i,s),'g--');
hold off
grid

xlabel('multipole'); ylabel('TT N_l (\muK^2)');
legend([h1(1),h2(1),h3(1),h4(1),h5(1),h6(1),h7(1),h8(1)],{'155pairs','310pairs','465pairs','620pairs','785pairs','930pairs','1085pairs','1240pairs'});
title('Mean of first 100 noise realizations');
%set(gca, MyFont{:});
saveas(gcf,'fig2b_B3_2017_TT_jack02.png');clf;




tot(1)=mean(sf(i).*aps1.Cs_l(i,s));
tot(2)=mean(sf(i).*aps2.Cs_l(i,s));
tot(3)=mean(sf(i).*aps3.Cs_l(i,s));
tot(4)=mean(sf(i).*aps4.Cs_l(i,s));
tot(5)=mean(sf(i).*aps5.Cs_l(i,s));
tot(6)=mean(sf(i).*aps6.Cs_l(i,s));
tot(7)=mean(sf(i).*aps7.Cs_l(i,s));
tot(8)=mean(sf(i).*aps8.Cs_l(i,s));

of(1)=mean(sf(i1).*aps1.Cs_l(i1,s));
of(2)=mean(sf(i1).*aps2.Cs_l(i1,s));
of(3)=mean(sf(i1).*aps3.Cs_l(i1,s));
of(4)=mean(sf(i1).*aps4.Cs_l(i1,s));
of(5)=mean(sf(i1).*aps5.Cs_l(i1,s));
of(6)=mean(sf(i1).*aps6.Cs_l(i1,s));
of(7)=mean(sf(i1).*aps7.Cs_l(i1,s));
of(8)=mean(sf(i1).*aps8.Cs_l(i1,s));

wn(1)=mean(sf(i2).*aps1.Cs_l(i2,s));
wn(2)=mean(sf(i2).*aps2.Cs_l(i2,s));
wn(3)=mean(sf(i2).*aps3.Cs_l(i2,s));
wn(4)=mean(sf(i2).*aps4.Cs_l(i2,s));
wn(5)=mean(sf(i2).*aps5.Cs_l(i2,s));
wn(6)=mean(sf(i2).*aps6.Cs_l(i2,s));
wn(7)=mean(sf(i2).*aps7.Cs_l(i2,s));
wn(8)=mean(sf(i2).*aps8.Cs_l(i2,s));

setwinsize(gcf,450,450)
loglog(tw,of,'-ob');
hold on
loglog(tw,wn,'-or');
xlim([5e4,1e6]);
%extrapolate point 1 based on map weight
ofp(1)=of(1); ofp(2:8)=of(1)*tw(1)./tw(2:8);
wnp(1)=wn(1); wnp(2:8)=wn(1)*tw(1)./tw(2:8);
loglog(tw,ofp,'b--')
loglog(tw,wnp,'r--')
hold off
grid


xlabel('total polarized map weight (155, 310, 465, 620, 785, 930, 1085, 1240) pairs respectively'); ylabel('mean TT N_l (\muK^2)');
legend({'ell ~ 40 bandpower','ell ~  350 bandpower'})
%set(gca, MyFont{:});
saveas(gcf,'fig3_B32017_TT_jack02.png');clf

setwinsize(gcf,450,450)
loglog(tw,tot,'-ok');
hold on
%loglog(tw,tot,'-or');
xlim([5e4,1e6]);
%extrapolate point 1 based on map weight
%ofp(1)=of(1); ofp(2:8)=mean(of(2))*mean(tw(2))./tw(2:8);
%wnp(1)=wn(1); wnp(2:8)=mean(wn(2))*mean(tw(2))./tw(2:8);
totp(1) = tot(1); totp(2:8) = tot(1)*tw(1)./tw(2:8)
%loglog(tw,ofp,'b--')
%loglog(tw,wnp,'r--')
loglog(tw,totp,'k--')
hold off
grid


xlabel('total polarized map weight (155, 310, 465, 620, 785, 930, 1085, 1240) pairs respectively'); ylabel('mean TT N_l (\muK^2)');
legend({'l > 10'})
%set(gca, MyFont{:});
saveas(gcf,'fig3_B32017_TT_combines_jack02.png');clf

