s2let_demo6

 s2let_demo6

0001 % s2let_demo6
0002 mapTfile = '../../../../ebsep/data/wmap_T_rot.fits';
0003 mapQfile = '../../../../ebsep/data/wmap_Q_rot.fits';
0004 mapUfile = '../../../../ebsep/data/wmap_U_rot.fits';
0005 maskTfile = '../../../../ebsep/data/wmap_T_mask_rot.fits';
0006 maskQUfile = '../../../../ebsep/data/wmap_QU_mask_rot.fits';
0007 
0008 [mapT_hpx, ~] = s2let_hpx_read_real_map(mapTfile);
0009 [maskT_hpx, ~] = s2let_hpx_read_real_map(maskTfile);
0010 [mapQ_hpx, ~] = s2let_hpx_read_real_map(mapQfile);
0011 [mapU_hpx, ~] = s2let_hpx_read_real_map(mapUfile);
0012 [maskQU_hpx, ~] = s2let_hpx_read_real_map(maskQUfile);
0013 
0014 Spin = 2 
0015 L = 192
0016 B = 3
0017 N = 3
0018 J_min = 2
0019 downsampling = false;
0020 mode = 3;
0021 PhiRot = pi;
0022 threshold = -1;
0023 
0024 % Analyse Q/U maps with Healpix routines to produce full-sky E/B alms
0025 [E_lm, B_lm] = s2let_hpx_map2alm_spin(mapQ_hpx, mapU_hpx, 'L', L);
0026 
0027 % Compute full sky Q+iU and Q-iU maps from E/B alms
0028 [QpU_lm] = ebsep_compute_QpU_fullsky(E_lm, B_lm);
0029 [QmU_lm] = ebsep_compute_QmU_fullsky(E_lm, B_lm);
0030 QpU = ssht_inverse(QpU_lm, L, 'spin', 2);
0031 QmU = ssht_inverse(QmU_lm, L, 'spin', -2);
0032 maskQU = s2let_hpx2mw(maskQU_hpx, 'L', L);
0033 maskQU(maskQU > threshold) = 1;
0034 maskQU(maskQU <= threshold) = 0;
0035 maskQU_lm = ssht_forward(maskQU, L);
0036 
0037 [maskQU, mask_wav, mask_scal] = ebsep_sample_mask(maskQU_lm, B, L, J_min, 'threshold', threshold, 'Downsample', downsampling);
0038 %
0039 f = QpU .* maskQU;
0040 %f = QpU;
0041 %
0042 
0043 J = s2let_jmax(L, B);
0044 zoomfactor = 1.7;
0045 ns = ceil(sqrt(2+J-J_min+1)) ;
0046 ny = 3; % ns - 1 + rem(2+J-J_min+1 , ns) ;
0047 nx = 3; % ns;
0048 
0049 maxfigs = nx*ny;
0050 pltroot = '../../../figs'
0051 configstr = ['Spin',int2str(Spin),'_N',int2str(N),'_L',int2str(L),'_B',int2str(B),'_Jmin',int2str(J_min)]
0052 
0053 [f_wav, f_scal] = s2let_transform_analysis_mw(f, 'B', B, 'J_min', J_min, 'N', N, 'Upsample', ~downsampling);
0054 
0055 bl = L;
0056 figure('Position',[100 100 1300 1000])
0057 subplot(nx, ny, 1);
0058 ssht_plot_mollweide(f, L, 'Mode', mode);
0059 title('Initial data')
0060 campos([0 0 -1]); camup([0 1 0]); zoom(zoomfactor)
0061 v = caxis;
0062 temp = max(abs(v));
0063 caxis([-temp temp])
0064 subplot(nx, ny, 2);
0065 if downsampling
0066     bl = min([ s2let_bandlimit(J_min-1,J_min,B,L), L ]);
0067 end
0068 ssht_plot_mollweide(f_scal .* mask_scal, bl, 'Mode', mode);
0069 campos([0 0 -1]); camup([0 1 0]); zoom(zoomfactor)
0070 v = caxis;
0071 temp = max(abs(v));
0072 caxis([-temp temp])
0073 title('Scaling fct')
0074 ind = 3
0075 for j = J_min:J
0076     for en = 1:N
0077         ind = ind + 1
0078         if ind <= maxfigs
0079             subplot(nx, ny, ind);
0080             if downsampling
0081                 bl = min([s2let_bandlimit(j, J_min, B, L), L]);
0082             end
0083             ssht_plot_mollweide(f_wav{j-J_min+1,en}.*mask_wav{j-J_min+1,1}, bl, 'Mode', mode);
0084             campos([0 0 -1]); camup([0 1 0]); zoom(zoomfactor)
0085             v = caxis;
0086             temp = max(abs(v));
0087             caxis([-temp temp])
0088             title(['Wavelet scale j=',int2str(j)-J_min+1,', n=',int2str(en)])
0089          end
0090     end
0091 end
0092 
0093 fname = [pltroot,'/s2let_demo6_', configstr, '_wmap_unmasked.png']
0094 print('-r200', '-dpng', fname)
0095 
0096