*** Deconvolve C_l for the Mask with the MASTER algorithm ***
June 27, 2010
binned cl version by Raphael Flauger: April 27, 2013
beam transfer function added: April 30, 2014
Calculation of the coupling matrix is modified to use symmetry of the square of Wigner 3j: July 26, 2018 by R. Makiya
fixed the bug in treatment of b_l (since ver. July 26): February 25, 2019 by R. Makiya

Here we provide a simple program for implementing the so-called 
MASTER algorithm [Hivon et al., ApJ, 567, 2 (2002)] to deconvolve 
the power spectrum for the effect of a mask applied to a map.

1. Untar "deconvolve_mask.tar"
2. Edit "Makefile"
   - FC: To use the correct fortran90 compiler
   - LDFLAGS: To use the correct link to LAPACK and BLAS library
3. Then, "make"
4. To run it, type "deconvolve_mask"

This code will ask you to enter:
1. the maximum multipole,
2. the bin size of the output cl,
3. the filename of the raw power spectrum measured from a masked sky map, and 
4. the filename of the power spectrum of the mask. 
5. the filename of the beam transfer function, including the effect of the pixel
   window function 

NOTE: The beam transfer function, b_l, is defined such that the beam-smeared
      power spectrum is given by C_l[smeared]=C_l[true]*b_l^2. This should
      also include the appropriate pixel-window function.

#3 and #4 can be computed by Healpix's "anafast", but you will need to convert 
the output of anafast (which is in the fits format) to the ASCII format. 

The code assumes that the input spectra are unbinned, and are given at all multipoles
from l=0 to l=[maximum multipole]

The format of the ASCII file must be:

cl at l=0
cl at l=1
cl at l=2
.
.

I.e., a single-column file, and cl is WITHOUT the "l(l+1)/twopi" factor.

The code will generate "multipole_deconvolved_cl.txt", which contains:

multipole, binned cl

I.e., a two-column file.

As an example, a cut-sky power spectrum file "cl_cutsky.txt" and a mask power
spectrum file "cl_of_mask.txt" are included (for lmax=1024). You may compare the output,
"multipole_deconvolved_cl.txt" with the full-sky power spectrum, "cl_fullsky.txt", to
see how well the full-sky power spectrum can be reconstructed by this algorithm.
The beam transfer function file included in this package "bl.txt" contains just unity.