% look at sigma(r) as a function of mask

clear

% setup the grid
np=32;
ad=calc_ad(20,np);
ad.t_r_deg=ad.t_r*180/pi;

if(0)
  % make a series of binary hit maps with varying radii but same total hits
  for i=1:10
    h(:,:,i)=real(ad.t_r_deg<i);
  end
  
  % find the initial total of each mask
  s=squeeze(sum(sum(h,1),2));

  % normalize so all hit maps have approx same total as in next case
  % - i.e. impose equal observing effort
  for i=1:10
    h(:,:,i)=241*h(:,:,i)/s(i);
  end
else
  % make a series of hit maps with varying taper but same total
  tf=[0,0.16,0.295,0.41,0.517,0.61,0.70,0.78,0.855,0.925,0.995];
  for i=1:11
    x=tukeywin(101,tf(i));
    t{i}=x(51:end);
  end
  r=linspace(5.5,10,11);
  for i=1:11
    h(:,:,i)=interp1(linspace(0,r(i),length(t{1})),t{i},ad.t_r_deg);
  end
  h(~isfinite(h))=0;
end
  
% check the total hits in each map - in both cases forced equal across the set
th=squeeze(sum(sum(h,1),2));

% plot the hit maps
figure(1); setwinsize(gcf,900,500); clf
t=ad.t_val_deg;
for i=1:size(h,3)
  subplot(3,4,i); imagesc(t,t,h(:,:,i)); colorbar
  title(sprintf('hit map %d - total %.1f',i,th(i)));
end

% take the weights as hits as is conventionally dones
% (inverse noise variance weighting)
w=h;

% for some realizations
for j=1:10000
  % make noise for each
  for i=1:size(h,3)
    n(:,:,i)=randn(np)./sqrt(h(:,:,i));
  end
  n(~isfinite(n))=0;

  for i=1:size(h,3)
    % take the weighted var of noise
    x=mean(rvec(((n(:,:,i)).*w(:,:,i)).^2));
    % apply scale factor such that signal level is preserved -
    % as we must do when measuring the absolute level of signal
    % on the sky
    % (do things in this order as close to what we do when also
    % including FT/SHT)
    nv(j,i)=x./mean(rvec(w(:,:,i).^2));
  end
end

% the std of the variance is the "bandpower error bar"
% in the noise dominated regime this maps directly to sigma(r)

% Crawford/Hivon definition of "the area" - doesn't seem to make sense
% - see https://arxiv.org/pdf/astro-ph/0105302.pdf eqn 9 and 17
%for i=1:size(h,3)
%  i1(i)=sum(rvec(h(:,:,i).^2));
%  i2(i)=sum(rvec(h(:,:,i).^4));
%  nz(i)=sum(rvec(h(:,:,i)>0));
%end
%ta=nz.*(i1.^2)./i2;

% SO paper definition of "the area" - this does make sense
% - see https://arxiv.org/pdf/1808.07445.pdf page 16
for i=1:size(h,3)
  i1(i)=sum(rvec(h(:,:,i))).^2;
  i2(i)=sum(rvec(h(:,:,i).^2));
end
so=i1./i2;

% plot results
figure(2); clf
plot(sqrt(so),std(nv),'.-')
c=regress(std(nv)',sqrt(so)')
xl=[0,30]; yl=[0,0.18];
xlim(xl); ylim(yl);
hold on
plot(xl,xl*c,'r');
hold off
xlabel('sqrt(calc. effective area) (arb units)')
ylabel('toy sim bandpower uncertainty (\sigma(r) proxy)');

return

% get the f_sky numbers for the 3 hit maps in question
hm1=read_fits_map('cmbs4_dat/expt_xx/04/rhits/n0512.fits');
hm2=read_fits_map('cmbs4_dat/expt_xx/04b/rhits/n0512.fits');
hm3=read_fits_map('cmbs4_dat/expt_xx/04c/rhits/n0512.fits');

[nansum(hm1.map)./length(hm1.map) nansum(hm2.map)./length(hm2.map) nansum(hm3.map)./length(hm3.map)]

nansum(hm1.map).^2/nansum(hm1.map.^2)/length(hm1.map)
