% toy-toy sims of mask diffusivity

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

% make circular and apodized hit maps with the same sum
h1=real(ad.t_r_deg<5.4);
h2=cosd(ad.t_r_deg*9).^2; h2(ad.t_r_deg*9>90)=0;
[sum(h1(:)) sum(h2(:))]

for i=1:10000
  % signal is uniform
  s=randn(np);
  
  % noise goes up around the edge
  n1=randn(np)./sqrt(h1); n1(~isfinite(n1))=0;
  n2=randn(np)./sqrt(h2); n2(~isfinite(n2))=0;
  
  % calc var of map*hits
  sv1(i)=mean((s(:).*h1(:)).^2);
  sv2(i)=mean((s(:).*h2(:)).^2);  
  nv1(i)=mean((n1(:).*h1(:)).^2);
  nv2(i)=mean((n2(:).*h2(:)).^2);
end

% 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)
sv1=sv1./mean(h1(:).^2);
sv2=sv2./mean(h2(:).^2);
nv1=nv1./mean(h1(:).^2);
nv2=nv2./mean(h2(:).^2);

% plot maps for last realization
setwinsize(gcf,400,600);
t=ad.t_val_deg;
subplot(3,2,1); imagesc(t,t,s); title('s'); caxis([-3,3]);
subplot(3,2,3); imagesc(t,t,h1); title('h_1'); caxis([0,1]);
subplot(3,2,4); imagesc(t,t,h2); title('h_2'); caxis([0,1]);
subplot(3,2,5); imagesc(t,t,n1); title('n_1'); caxis([-3,3]);
subplot(3,2,6); imagesc(t,t,n2); title('n_2'); caxis([-3,3]);

setwinsize(gcf,400,400);
subplot(2,2,1); imagesc(t,t,s.*h1); title('sh_1'); caxis([-3,3]);
subplot(2,2,2); imagesc(t,t,s.*h2); title('sh_2'); caxis([-3,3]);
subplot(2,2,3); imagesc(t,t,n1.*h1); title('n_1h_1'); caxis([-3,3]);
subplot(2,2,4); imagesc(t,t,n2.*h2); title('n_2h_2'); caxis([-3,3]);

setwinsize(gcf,400,400);
subplot(2,2,1); hfill(sv1,30);
title(sprintf('V_{s1}: m=%.2f, \\sigma=%.3f, \\sigma/m=%.3f',mean(sv1),std(sv1),std(sv1)/mean(sv1)));
subplot(2,2,2); hfill(sv2,30);
title(sprintf('V_{s2}: m=%.2f, \\sigma=%.3f, \\sigma/m=%.3f',mean(sv2),std(sv2),std(sv2)/mean(sv2)));
subplot(2,2,3); hfill(nv1,30);
title(sprintf('V_{n1}: m=%.2f, \\sigma=%.3f, \\sigma/m=%.3f',mean(nv1),std(nv1),std(nv1)/mean(nv1)));
subplot(2,2,4); hfill(nv2,30);
title(sprintf('V_{n2}: m=%.2f, \\sigma=%.3f, \\sigma/m=%.3f',mean(nv2),std(nv2),std(nv2)/mean(nv2)));


% do the above again using the actual masks
h1=read_fits_map('cmbs4_dat/expt_xx/04/rhits/n0512.fits'); h1=ud_grade(h1,64);
h2=read_fits_map('cmbs4_dat/expt_xx/04b/rhits/n0512.fits'); h2=ud_grade(h2,64);
h3=read_fits_map('cmbs4_dat/expt_xx/04c/rhits/n0512.fits'); h3=ud_grade(h3,64);
h1=h1.map; h2=h2.map; h3=h3.map;
h1(isnan(h1))=0;
h2(isnan(h2))=0;
h3(isnan(h3))=0;
[nansum(h1) nansum(h2) nansum(h3)]

np=size(h1);

for i=1:1000
  % signal is uniform
  s=randn(np);
  
  % noise goes up around the edge
  n1=randn(np)./sqrt(h1); n1(~isfinite(n1))=0;
  n2=randn(np)./sqrt(h2); n2(~isfinite(n2))=0;
  n3=randn(np)./sqrt(h3); n3(~isfinite(n3))=0;
  
  % calc var of map*hits
  sv1(i)=mean((s(:).*h1(:)).^2);
  sv2(i)=mean((s(:).*h2(:)).^2);  
  sv3(i)=mean((s(:).*h3(:)).^2);  
  nv1(i)=mean((n1(:).*h1(:)).^2);
  nv2(i)=mean((n2(:).*h2(:)).^2);
  nv3(i)=mean((n3(:).*h3(:)).^2);
end

% 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)
sv1=sv1./mean(h1(:).^2);
sv2=sv2./mean(h2(:).^2);
sv3=sv3./mean(h3(:).^2);
nv1=nv1./mean(h1(:).^2);
nv2=nv2./mean(h2(:).^2);
nv3=nv3./mean(h3(:).^2);

setwinsize(gcf,800,500);
subplot(2,3,1); hfill(sv1,30);
title(sprintf('04 sig: m=%.2f, \\sigma=%.3f, \\sigma/m=%.3f',mean(sv1),std(sv1),std(sv1)/mean(sv1)));
subplot(2,3,2); hfill(sv2,30);
title(sprintf('04b sig: m=%.2f, \\sigma=%.3f, \\sigma/m=%.3f',mean(sv2),std(sv2),std(sv2)/mean(sv2)));
subplot(2,3,3); hfill(sv3,30);
title(sprintf('04c sig: m=%.2f, \\sigma=%.3f, \\sigma/m=%.3f',mean(sv3),std(sv3),std(sv3)/mean(sv3)));

subplot(2,3,4); hfill(nv1,30);
title(sprintf('04 noi: m=%.2f, \\sigma=%.3f, \\sigma/m=%.3f',mean(nv1),std(nv1),std(nv1)/mean(nv1)));
subplot(2,3,5); hfill(nv2,30);
title(sprintf('04b noi: m=%.2f, \\sigma=%.3f, \\sigma/m=%.3f',mean(nv2),std(nv2),std(nv2)/mean(nv2)));
subplot(2,3,6); hfill(nv3,30);
title(sprintf('04c noi: m=%.2f, \\sigma=%.3f, \\sigma/m=%.3f',mean(nv3),std(nv3),std(nv3)/mean(nv3)));


% plot signal and noise curves to look at weirdness
load aps/1699/xxx5_04p00_jack0.mat
s04=aps;
load aps/4321/xxx5_04bp00_jack0.mat
s04b=aps;
load aps/4322/xxx5_04cp00_jack0.mat
s04c=aps;

load aps/1699/xxx6_04p00_jack0.mat
n04=aps;
load aps/4321/xxx6_04bp00_jack0.mat
n04b=aps;
load aps/4322/xxx6_04cp00_jack0.mat
n04c=aps;

b=7
setwinsize(gcf,600,600); clf
semilogy(aps(b).l(2:10),mean(s04c(b).Cs_l(2:10,4,:),3),'g'); ylabel('BB D_l')
hold on
plot(aps(b).l(2:10),mean(s04(b).Cs_l(2:10,4,:),3),'b');
plot(aps(b).l(2:10),mean(s04b(b).Cs_l(2:10,4,:),3),'r');
plot(aps(b).l(2:10),mean(n04c(b).Cs_l(2:10,4,:),3),'g--'); ylabel('noise BB D_l')
plot(aps(b).l(2:10),mean(n04(b).Cs_l(2:10,4,:),3),'b--');
plot(aps(b).l(2:10),mean(n04b(b).Cs_l(2:10,4,:),3),'r--');
hold off
title({'155GHz - masks: blue=circular 04, red=Chile 04b, green=Pole 04c',...
       'solid=LCDM signal (A_L=1), dashed=noise'});
