

program poker_mbb
  use healpix_types
  use extension
  use misc_utils

  use poker_types
  use poker_input
  use poker_tools

  implicit none

  character(len=200) :: parafile
  integer :: i_parfile
  type(poker_str) :: pars

  integer(I4B) :: i, j, ios, n_iter, npix, b, b1, m, n, m1, n1, np, ii, jj, i_alpha, j_alpha, itap, nbins
  integer(I4B) :: myid, nproc, ierr, npix_proc, i_start, i_end, naxis1, naxis2
  real(DP),     dimension(:,:), allocatable :: mtt_bb_loc, mtt_bb
  !real(DP),     dimension(:,:), allocatable :: mtt_bb_loc_1, mtt_bb_loc_x, mtt_bb_1, mtt_bb_x
  real(DP),     dimension(:), allocatable :: nb, nb_loc, bintab

  real(DP), dimension(0:1) :: kv, k1v, x, y
  real(DP) :: cosphi, sinphi, cos2phi, sin2phi, cos2phi2, sin2phi2, sin4phi, sigma, bk

  real(DP),     dimension(:,:),   allocatable :: pix_mn, map_k_binning, map_b
  complex(DPC), dimension(:,:), allocatable :: w_mn!, beam_mn, beam1_mn!, kernel
  real(DP),     dimension(:,:), allocatable :: beam_mn, beam1_mn
  real(DP) :: err, ww, kx, ky, k, k1, p, q, p1, kx1, ky1, kk, r, rm1, dummy_min, dummy_max
  complex(DPC) :: ww_c

  character(len=2) :: si, sj
  character(len=30) :: code
  character(len=200) :: file

  real(SP)     :: time0, time1

  !fftw parameter
  integer(I8B) :: plan
  !-------------------------------------------------------------------------------------------------

  INCLUDE 'fftw3.f'

#ifdef PARA
  INCLUDE 'mpif.h'
  CALL mpi_init( ierr)
  CALL mpi_comm_size( MPI_COMM_WORLD, nproc, ierr)
  CALL mpi_comm_rank( MPI_COMM_WORLD, myid,  ierr)
#else
  myid = 0
  nproc = 1
#endif

  code = "poker_mbb"

  if (myid == 0) call cpu_time(time0)

  i_parfile = 1
  call poker_read_params(i_parfile, myid, pars, parafile)
  if ((myid == 0) .and. (pars%verb /= 0)) print*, "input parameters done"

  !get map size parameters
  call get_naxis12( trim(pars%mask), naxis1, naxis2)
  npix = naxis1 * naxis2

  !get binning parameters
  call get_naxis1( trim(pars%input_bintab), nbins)
  nbins = nbins - 1 !bintab has nbins+1 elements

  allocate( nb_loc(0:nbins-1), stat=ios)
  call assert_alloc( ios, code, "nb_loc")

  allocate( bintab(0:nbins+1-1), stat=ios)
  call assert_alloc( ios, code, "bintab")

  allocate( mtt_bb_loc(0:nbins-1,0:nbins-1), stat=ios)
  call assert_alloc( ios, code, "mtt_bb_loc")
  mtt_bb_loc = 0.0d0 !init

  allocate( w_mn( 0:naxis1-1, 0:naxis2-1), stat=ios)
  call assert_alloc( ios, code, "w_mn")
  w_mn = 0.0d0

  allocate( beam_mn(0:naxis1-1, 0:naxis2-1), stat=ios)
  call assert_alloc( ios, code, "beam_mn")
  beam_mn = 0.0d0

  allocate( beam1_mn(0:naxis1-1, 0:naxis2-1), stat=ios)
  call assert_alloc( ios, code, "beam1_mn")
  beam1_mn = 0.0d0

  allocate( pix_mn(0:naxis1-1, 0:naxis2-1), stat=ios)
  call assert_alloc( ios, code, "pix_mn")
  pix_mn = 0.0d0

  allocate( map_k_binning( 0:naxis1-1, 0:naxis2-1), stat=ios)
  call assert_alloc( ios, code, "map_k_binning")
  map_k_binning = 0.0d0

  if( myid == 0) then
     
     ! Useful quantity to bin 2D power spectra
     allocate( map_b( 0:naxis1-1, 0:naxis2-1), stat=ios)
     call assert_alloc( ios, code, "map_b")
     map_b = 0.0d0

     !Read bin boundaries
     call read_fits_1d( trim(pars%input_bintab), nbins+1, bintab)

     !Read binning weights
     call readimage( trim(pars%map_k_binning), naxis1, naxis2, map_k_binning)
     if (pars%verb /= 0) print*, "read "//trim(pars%map_k_binning)//" done."

     !Read mask (use pix_mn array to save memory)
     call readimage( trim(pars%mask), naxis1, naxis2, pix_mn)
     if (pars%verb /= 0) print*, "read "//trim(pars%mask)//" done."
     w_mn = pix_mn !passing R8 values to complex array

     ! Fourier transform the mask
     call dfftw_plan_dft_2d( plan, naxis1, naxis2, w_mn, w_mn, -1, fftw_estimate)
     call dfftw_execute_dft( plan, w_mn, w_mn)
     call dfftw_destroy_plan(plan) 

     !Correct for FFTW normalization
     w_mn = w_mn/npix

     !read beam transfer function (use pix_mn array to save memory)
     call readimage( trim(pars%beam), naxis1, naxis2, pix_mn)
     if (pars%verb /= 0) print*, "read "//trim(pars%beam)//" done."
     beam_mn = pix_mn !passing R8 values to complex array

     call readimage( trim(pars%beam1), naxis1, naxis2, pix_mn)
     if (pars%verb /= 0) print*, "read "//trim(pars%beam1)//" done."
     beam1_mn = pix_mn !passing R8 values to complex array

  endif

#ifdef PARA
  call MPI_Barrier( MPI_COMM_WORLD, ierr)
  call MPI_Bcast( bintab,        nbins+1, MPI_Double_Precision, 0, MPI_COMM_WORLD, ierr)
  call MPI_BCAST( w_mn,          npix,      MPI_double_complex,   0, MPI_COMM_WORLD, ierr)
  call MPI_BCAST( beam_mn,       npix,      MPI_double_complex,   0, MPI_COMM_WORLD, ierr)
  call MPI_BCAST( beam1_mn,   npix,      MPI_double_complex,   0, MPI_COMM_WORLD, ierr)
  call MPI_BCAST( map_k_binning, npix,      MPI_double_precision, 0, MPI_COMM_WORLD, ierr)
  call MPI_Barrier( MPI_COMM_WORLD, ierr)
#endif

  if ((myid == 0) .and. (pars%verb /= 0)) print*, "MPI_Bcast done"

  ! Square pixel window function
  pix_mn = 1.0d0 !init to default
  if (pars%include_pix_window_function /= 0) then
     do i=0, npix-1
        call pix2mn( i, naxis1, naxis2, m, n)
        call mn2k(  m, n, naxis1, naxis2, kx, ky, k)
        if ( (kx /= 0.) .AND. (ky /= 0.)) then
           pix_mn(m,n) = dsin( PI*kx)*dsin( PI*ky) / (PI*kx * PI*ky)
        endif
     enddo
  endif

  ! Distribute pixels amond CPUs
  npix_proc  = npix/nproc !integer
  i_start = myid*npix_proc
  if (myid /= nproc-1) then
     i_end = i_start + npix_proc -1
  else
     i_end = npix-1
  endif

  !----------------------------------------- Main loop -----------------------------------------
  nb_loc       = 0.0d0 ! Init normalization factor in binning matrix P

  if ((myid == 0) .and. (pars%verb /= 0)) print*, "Starting main loop"
  if (pars%simul_type == 1) then

     do i=i_start, i_end
        call pix2mn( i, naxis1, naxis2, m, n)
        call mn2k( m, n, naxis1, naxis2, kx, ky, k)
        call k2bin( k, nbins, bintab, b)

        if (b >= 0) then
           p = map_k_binning( m, n)
           
           nb_loc(b) = nb_loc(b) + 1.0d0

           do j=0, npix-1
              call pix2mn( j, naxis1, naxis2, m1, n1)
              call mn2k( m1, n1, naxis1, naxis2, kx, ky, k)
              call k2bin( k, nbins, bintab, b1)
              
              if (b1 >= 0) then
                 
                 q = 1.0d0/map_k_binning(m1,n1)

                 if ( (m1 <= m) .and. (n1 <= n)) ww = abs( pix_mn(m1,n1) * sqrt( beam_mn( m1, n1) * beam1_mn( m1, n1)) * w_mn( m-m1,              n-n1))**2
                 if ( (m1 <= m) .and. (n1 >  n)) ww = abs( pix_mn(m1,n1) * sqrt( beam_mn( m1, n1) * beam1_mn( m1, n1)) * w_mn( m-m1,       naxis2+n-n1))**2
                 if ( (m1 >  m) .and. (n1 <= n)) ww = abs( pix_mn(m1,n1) * sqrt( beam_mn( m1, n1) * beam1_mn( m1, n1)) * w_mn( naxis1+m-m1,       n-n1))**2
                 if ( (m1 >  m) .and. (n1 >  n)) ww = abs( pix_mn(m1,n1) * sqrt( beam_mn( m1, n1) * beam1_mn( m1, n1)) * w_mn( naxis1+m-m1,naxis2+n-n1))**2
                 mtt_bb_loc(b1,b) = mtt_bb_loc(b1,b) + p*ww*q

              endif !b1>=0
           enddo!j=0, npix-1
        endif !b>=0
     enddo !i=i_start, i_end

  else
     print*, "poker_mbb : simul_type /= 1"
     print*, "polarization not implemented yet"
     stop
  endif!simul_type==1

  if ((myid == 0) .and. (pars%verb /= 0)) print*, "Main loop done."

  if (myid == 0) then
     allocate( mtt_bb(0:nbins-1,0:nbins-1), stat=ios)
     call assert_alloc( ios, code, "mtt_bb")
     mtt_bb = 0.0d0

     allocate( nb(0:nbins-1), stat=ios)
     call assert_alloc( ios, code, "nb")
     nb = 0.0d0
  endif

#ifdef PARA
  call MPI_Barrier( MPI_COMM_WORLD, ierr)
  call MPI_Reduce( mtt_bb_loc, mtt_bb, nbins**2,     MPI_Double_Precision, MPI_SUM, 0, MPI_COMM_WORLD, ierr)
  call MPI_Reduce( nb_loc,     nb,     nbins,        MPI_Double_Precision, MPI_SUM, 0, MPI_COMM_WORLD, ierr)
  call MPI_Barrier( MPI_COMM_WORLD, ierr)
#else
  mtt_bb    = mtt_bb_loc
  nb        = nb_loc
#endif

  if ((myid == 0) .and. (pars%verb /= 0)) print*, "Reduce done."

  if (myid == 0) then

     !Compute map_b
     do i=0, npix-1
        call pix2mn( i, naxis1, naxis2, m, n)
        call mn2k( m, n, naxis1, naxis2, kx, ky, k)
        call k2bin( k, nbins, bintab, b)
        map_b(m,n) = b
     enddo

     !Output
     call writeimage( trim(pars%file_map_b), naxis1, naxis2, map_b)

     ! Account for normalization factor in P (number of modes per bin, see Hivon et al, Eq. (20))
     do b=0, nbins-1
        if (nb(b) /= 0.) then
           mtt_bb(:,b) = mtt_bb(:,b)/nb(b)
        endif
     enddo

     !Output
     call writeimage( trim(pars%file_mtt_bb_x), nbins, nbins, mtt_bb)
     if ((myid == 0) .and. (pars%verb /= 0)) print*, "wrote mtt_bb (x)"

     call cpu_time(time1)
     open(  2, file=trim(pars%report_card), status='unknown')
     write( 2, *) "npix  = ", npix
     write( 2, *) "nproc = ", nproc
     write( 2, *) "CPU time (sec): ", time1 - time0
     close( 2)

     if (myid ==0) print*, "CPU time (sec): ", time1 - time0

  endif

#ifdef PARA
  CALL mpi_finalize(ierr)
#else

  deallocate( nb_loc)
  deallocate( bintab)
  deallocate( mtt_bb_loc)
  deallocate( w_mn)
  deallocate( beam_mn)
  deallocate( pix_mn)
  deallocate( map_k_binning)
  deallocate( mtt_bb)
  deallocate( nb)

#endif

end program poker_mbb
