PROGRAM isw_anomalies
!Modified by Anais Rassat, July 2013. Cleaned for public release.
!modified by Anais Rassat, June 2012
!Written by Anais Rassat, Feb 2011
!Purpose: calculate mirror parity (odd and even)

! Example of how to compile
! gfortran anomalies_parity_public.f90 -O3 -L/Users/arassat/Work/Healpix_2.15a/libf90 -I/Users/arassat/Work/Healpix_2.15a/includef90 -L/Users/arassat/Work/Libraries/cfitsio3230/lib -fopenmp -lhealpix -lcfitsio -o anomalies_parity_public.o

use healpix_types
use alm_tools
use pix_tools
use udgrade_nr
use fitstools
use rngmod

integer(i4b) :: lmax, nmaps=1
integer(i4b) :: lorder, lorder2, lordermax=5
integer(I4B), parameter::nside=512, npix=12*nside**2
integer(I4B), parameter::nside_parity=64, npix_parity=12*nside_parity**2
integer(I4B), parameter::nside_degrade=8, npix_degrade=12*nside_degrade**2
integer(I4B) :: lmax_parity
integer(I4B)  :: input! Choice of CMB map 
integer(I4B)  :: subtractisw, addnvss! Choice of ISW subtraction

real(sp) :: splus, smin

real(dp), allocatable, dimension(:,:) :: dw8
real(dp), dimension(2) :: costheta
  Real(dp):: AVG, SD

real(dp), allocatable, dimension(:,:) :: map, tisw ,tisw2, doppler
real(dp), allocatable, dimension(:,:) :: mapdegrade, tiswdegrade, tisw2degrade
real(dp), allocatable, dimension(:) :: smap
complex(dpc), allocatable, dimension(:,:,:) :: alm, origalm

character(len=80), dimension(60)::header
character(len=15), allocatable, dimension(:) :: cmbnames
header(:) = ''

!=================================
!Definitions and allocations
!=================================
addnvss = 1 ! if =0 use only 2mass maps, if =1 use also nvss isw map
lmax = 15
lmax_parity = 5
!DP
allocate(dw8(1:2*nside, 1:3))
dw8 = 1.0_dp
costheta =(-1.d0,1.d0)

!SP
allocate(map(0:nside2npix(nside)-1,1:3)) ! CMB map
allocate(mapdegrade(0:nside2npix(nside_degrade)-1,1:3)) ! CMB map
allocate(smap(0:npix_parity-1)) ! parity map
allocate(tisw(0:nside2npix(nside)-1,1:1)) ! ISW map for 2MASS
allocate(tisw2(0:nside2npix(nside)-1,1:1)) ! ISW map for NVSS
allocate(tiswdegrade(0:nside2npix(nside_degrade)-1,1:1)) ! CMB map
allocate(tisw2degrade(0:nside2npix(nside_degrade)-1,1:1)) ! CMB map
allocate(doppler(0:nside2npix(nside_degrade)-1,1:1)) ! kinetic Doppler quadrupole map
allocate(alm(1:3, 0:lmax, 0:lmax))
allocate(origalm(1:3, 0:lmax, 0:lmax))
allocate(cmbnames(1:nmaps))

alm = 0.0
origalm = 0.0

write(*,*) '-------------------------'
write(*,*) '-------------------------'
write(*,*) '-------------------------'

!=================================
! Read in maps
!=================================
! First read in the kinetic Doppler quadrupole map
call input_map('dquad_8.fits', doppler, npix_degrade,1)!NESTED This is the map in nside=8
call convert_nest2ring(nside_degrade, doppler) !dquad is now in ring order 
!--------------------------------
!Now chose which of the CMB maps you would like to test
!--------------------------------
Do subtractisw = 1, 2
   input = 1
   call input_map('../WMAP/W9_inpainted_sparsity_jan13mono_pows_100_l2_to_l5.fits', map, npix, 1) ! NESTED
   call convert_nest2ring(nside, map) ! Convert to Ring
   map = map / 1d6 / 2.725d0 ! Convert to unitless values
   write(*,*) 'CMB map: WMAP ILC 9yr (inpainted).'
   cmbnames(input) = 'W9 (inp)'
   call input_map('../2MASS/2mass_inpainted_theory0wmap58july13.fits', tisw, npix, 1) ! NESTED
   call convert_nest2ring(nside, tisw) !Convert ISW map to RING
   call input_map('../NVSS/nvss_inpainted_theory0wmap58july13.fits', tisw2, npix, 1) ! NESTED
   call convert_nest2ring(nside, tisw2) !Convert ISW map to RING
   ! For parity calculation - downgrade map to nside8=8
   call udgrade_ring(map, nside, mapdegrade, nside_degrade)         
   call udgrade_ring(tisw, nside, tiswdegrade, nside_degrade)         
   call udgrade_ring(tisw2, nside, tisw2degrade, nside_degrade)         
   !--------------------------------
   !Now substract the kinetic Doppler quadrupole
   !--------------------------------
   !Note Doppler quadrupole map is already in nside=8
   mapdegrade(:,1:1) = mapdegrade(:,1:1) - doppler(:,:) !Both maps should be in RING 
   write(*,*) 'CMB map: Kinetic Doppler map has been subtracted.'
   !--------------------------------
   
   if (subtractisw == 2) then 
      tiswdegrade = tiswdegrade / 1d6/2.725 ! convert to unitless values
      tisw2degrade = tisw2degrade / 1d6 / 2.725 ! convert to unitless values
      mapdegrade(:,1:1) = mapdegrade(:,1:1) - tiswdegrade(:,:) !Both maps should be in RING 
      if (addnvss == 1) then
         write(*,*) 'Subtracting both 2MASS and NVSS map'
         mapdegrade(:,1:1) = mapdegrade(:,1:1) - tisw2degrade(:,:) ! remove NVSS data too if required
      endif
      write(*,*) 'CMB map: ISW map has been subtracted.'
   else if (subtractisw == 1) then 
      write(*,*) 'CMB map: ISW map has NOT been subtracted.'
   endif
   !=================================
   ! Calculate Alm's of unrotated map
   !=================================
   call map2alm(nside_degrade, lmax, lmax, mapdegrade, alm, costheta, dw8) !Calculte alm's
   origalm = alm
   
   !=========================
   ! Calculate parity map
   !=========================
   call make_parity_map(smap, alm,nside_parity, npix_parity, lmax, lmax_parity, smin, splus) 
   write(*,*) 'Parity map calculated'
   write(*,*) 'S+ = ', splus
   write(*,*) 'S- = ', smin
enddo
END PROGRAM isw_anomalies

!======================================================================
!Following is to make new parity map from input map
SUBROUTINE make_parity_map(parity_map, orig_alm,nside,npix, lmax_alm, lmax, smin, splus)
! Updated by Anais Rassat, July 2013. Cleaned for public release.
! Written by Anais Rassat, July 2012
!INPUT: orig_alm: alm's of original CMB or other map
!       lmax: lmax for calculation of S statistic
!OUTPUT: parity_map: parity map for whole sky
use healpix_types
use alm_tools
use pix_tools
use fitstools

integer(i4b) :: lmax, nsim, isim, lmax_alm
integer(I4B) ::nside, npix
integer(I4B) :: lorder, m, mtest
real(sp), allocatable, dimension(:) :: splusind, sminind
real(sp) :: splus, smin

real(dp) :: cl
real(dp), dimension(0:npix-1) :: parity_map
real(dp) :: psi, theta, phi
real(dp) :: sstat, avg, sd

complex(dpc), dimension(1:3, 0:lmax_alm, 0:lmax_alm) :: orig_alm, alm

allocate(splusind(0:npix-1))
allocate(sminind(0:npix-1))

parity_map(:) = 0.d0 ! first create a parity map that is just zero everywhere

do isim=0, npix-1 ! run over different pixel values, each pixel corresponds to a different direction \hat{n})
   call pix2ang_ring(nside,isim,theta,phi)!calculate theta,phi corresponding to pixel isim
   psi = 0.d0 ! our rotation does not depend on the value of psi
   alm = orig_alm
   call rotate_alm(lmax_alm, alm, phi,theta ,psi)!Rotate the map / alm's
   sstat = 0.d0 ! one sstat value for each pixel direction
   do lorder = 2, lmax 
      cl = 0.d0 ! one cl value per lorder
      ! calculate cl, cl is invariant under rotation so only need to calculate once for each lorder
      cl = abs(alm(1,lorder,0))**2 ! start with value for m = 0
      do m = 1, lorder
         cl = cl + 2.d0*(abs(alm(1,lorder,m))**2) ! sum_m 2*|alm|^2 ; value for m>0
      enddo
      cl = cl / (2.d0 * lorder + 1.d0)
      do m = -lorder, lorder
         mtest = abs(m)
         sstat = sstat + (-1)**(lorder+m)*abs(alm(1,lorder,mtest))**2/cl !alm here are rotated ones
      enddo
   enddo
   parity_map(isim) = sstat - (lmax - 1.d0) ! one sstat value for each pixel direction
enddo
!now calculate S+ and S-
splusind = 0.d0
splusind = maxloc(parity_map) !find location of maximum quadrupole
splus = parity_map(int(splusind(0)))
sminind = 0.d0
sminind = minloc(parity_map) !find location of maximum quadrupole
smin = parity_map(int(sminind(0)))

call meansd_map(parity_map, npix, avg, sd)

splus = abs(splus - avg) / sd
smin = abs(smin - avg) / sd
write(*,*) splus, smin
END SUBROUTINE !make_parity_map
!======================================================================

! Taken from http://programming.ccp14.ac.uk/fortran-resources/~cgp/prof77/node3.html
! And then translated into Fortran 90 by Anais Rassat.
SUBROUTINE MEANSD(X, NPTS, AVG, SD) 

  use healpix_types
  use alm_tools
  use pix_tools
  use fitstools
  
  !Define variables
  Implicit None
  Integer(i4b) :: NPTS, I
  Real(dp):: AVG, SD, Sum, Sumsq
  Real(dp), dimension(1:npts) :: X
  
  SUM   = 0.0 
  SUMSQ = 0.0 
!First calculate Average
  DO  I = 1,NPTS 
     SUM   = SUM   + X(I) 
  END DO
  AVG = SUM / NPTS 
  
!Now calculate Standard deviation
  DO I = 1, NPTS
     SUMSQ = SUMSQ + (X(I)-AVG)**2 
  END DO

  SD  = SQRT(SUMSQ/(NPTS-1))
END SUBROUTINE MEANSD



! Taken from http://programming.ccp14.ac.uk/fortran-resources/~cgp/prof77/node3.html
! And then translated into Fortran 90 by Anais Rassat. 
SUBROUTINE MEANSD_MAP(X, NPTS, AVG, SD)  !Same as meansd but array start from zero

  use healpix_types
  use alm_tools
  use pix_tools
  use fitstools
  
  !Define variables
  Implicit None
  Integer(i4b) :: NPTS, I
  Real(dp):: AVG, SD, Sum, Sumsq
  Real(dp), dimension(0:npts-1) :: X
  
  SUM   = 0.0 
  SUMSQ = 0.0 
!First calculate Average
  DO  I = 0,NPTS-1 
     SUM   = SUM   + X(I) 
  END DO
  AVG = SUM / NPTS 
  
!Now calculate Standard deviation
  DO I = 0, NPTS-1
     SUMSQ = SUMSQ + (X(I)-AVG)**2 
  END DO

  SD  = SQRT(SUMSQ/(NPTS-1))
END SUBROUTINE MEANSD_MAP