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, find the axis of evil

!Example of how to compile
! gfortran anomalies_aoe_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_aoe_public.o

use healpix_types
use alm_tools
use pix_tools
use fitstools
use rngmod

integer(i4b) :: lmax, nmaps=1, mmax_aoe
integer(i4b) :: lorder, lorder2, nsimisw, isimisw, nsimiswmax, lordermax=5
integer(i4b) :: test
integer(I4B), parameter::nside=512, npix=12*nside**2
integer(I4B)  :: input ! Choice of CMB map 
integer(I4B)  :: subtractisw, addnvss! Choice of ISW subtraction

real(dp), allocatable, dimension(:,:) :: dw8
real(dp), dimension(2) :: costheta
real(dp) :: glmax2, gbmax2, glmax3, gbmax3
real(dp), dimension(3) :: vector2, vector4
real(Dp), allocatable, dimension(:,:) :: vector3
real(dp) :: dotproduct, testangle
real(dp), dimension(6) :: angle
  Real(dp):: AVG, SD

real(sp) :: rl
real(sp), allocatable, dimension(:,:) :: map, tisw ,tisw2, doppler!,dquad
complex(spc), allocatable, dimension(:,:,:) :: alm, origalm

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

!type(planck_rng) :: rng_handle

!=================================
!Definitions and allocations
!=================================
subtractisw = 1 ! 1 - do not subtract ISW, 2 - subtract ISW
addnvss = 1 ! if =0 use only 2mass maps, if =1 use also nvss isw map
randomisw = 1 ! 1 - do not create ISW simulations 2 - Create ISW simulations
lmax = 15
nsimiswmax = 1 ! Number of ISW simulations
!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(tisw(0:nside2npix(nside)-1,1:1)) ! ISW map for 2MASS
allocate(tisw2(0:nside2npix(nside)-1,1:1)) ! ISW map for NVSS
allocate(doppler(0:nside2npix(nside)-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))
allocate(vector3(3,1:lordermax))

alm = 0.0
origalm = 0.0

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

!=================================
! Read in maps
!=================================
! First read in the kinetic Doppler quadrupole map
call input_map('dquad_512.fits', doppler, npix, 1) ! NESTED order
call convert_nest2ring(nside, doppler) !dquad is now in ring order 
!--------------------------------
!Now chose which of the CMB maps you would like to test
!--------------------------------
Do subtractisw = 1,2! 2
   input = 1
   call input_map('../WMAP/W9_inpainted_sparsity_jan13mono_pows_100.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
   !--------------------------------
   !Now substract the kinetic Doppler quadrupole
   !--------------------------------
   map(:,1:1) = map(:,1:1) - doppler(:,:) !Both maps should be in RING 
   write(*,*) 'CMB map: Kinetic Doppler map has been subtracted.'
   if (subtractisw == 2) then 
      tisw = tisw / 1d6/2.725 ! convert to unitless values
      tisw2 = tisw2 / 1d6 / 2.725 ! convert to unitless values
      map(:,1:1) = map(:,1:1) - tisw(:,:) !Both maps should be in RING 
      if (addnvss == 1) then
         map(:,1:1) = map(:,1:1) - tisw2(:,:) ! 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, lmax, lmax, map, alm, costheta, dw8) !Calculte alm's
   origalm = alm
   
   do lorder=2, lordermax 
      call findaoe(alm, lmax , lorder, vector2, glmax2, gbmax2, rl, mmax_aoe) 
      write(*,*), 'lorder :', lorder, 'axis: ', glmax2, gbmax2, 'rl: ', rl, ' mmax: ', mmax_aoe
      vector3(:,lorder) = vector2
   enddo
   test = 0
   do lorder=2, lordermax
      do lorder2= lorder+1, lordermax
         vector2 = vector3(:, lorder)
         vector4 = vector3(:,lorder2)
         dotproduct =vector2(1)*vector4(1) + vector2(2)*vector4(2) + vector2(3)*vector4(3)
         testangle = acos(dotproduct)*rad2deg
         if (testangle .gt. 90) then 
            testangle = 180.-testangle
         endif
         write(*,*), test, lorder, lorder2, testangle
         test = test+1
         angle(test) = testangle
      enddo
   enddo
   Call Meansd(angle, 6, avg, sd)
   write(*,*) ' Angle: ', avg, ' +/- ', sd
enddo
END PROGRAM isw_anomalies
    

! Following is for Axis of Evil
SUBROUTINE findaoe(alm, lmax, lorder, vector, glmax, gbmax, rl, mmax_aoe)
use healpix_types
use alm_tools
use pix_tools
use fitstools

! Written by Anais Rassat
! July 2012
! Cleaned July 2013 for public release
! Based on findaxis code for alignment of octopole and quadrupole written by Anais Rassat
! 
! Define variables
integer(i4b) :: lmax, nsim, isim, imax, mmax, mmax_aoe!, imax2
integer(I4B), parameter::nside=128, npix=12*nside**2
integer(I4B) :: lorder, morder, mvalue!, lorder2

real(sp) :: test,cl, rl !test2
real(sp), allocatable, dimension(:,:) :: stat!, stat2
real(sp), allocatable, dimension(:) :: maxstatloc, maxstatloc2


real(dp) :: psi, theta, phi
real(dp) :: thetamax, phimax!, thetamax2, phimax2
real(dp) :: glmax, gbmax!, glmax2, gbmax2
real(dp), dimension(3) :: vector!, vector2

complex(spc), allocatable, dimension(:,:,:) :: origalm
complex(spc), dimension(1:3, 0:lmax, 0:lmax) :: alm


! Allocate variables
nsim = nside2npix(nside)-1
mmax = 2*lorder+1 ! number of m-modes to search for preferred axis
allocate(stat(0:nsim-1,0:mmax))
allocate(origalm(1:3, 0:lmax, 0:lmax))
allocate(maxstatloc(0:nsim-1))
allocate(maxstatloc2(0:nsim-1))

maxstatloc = 0.0
maxstatloc2 = 0.0
origalm = 0.0 
origalm = alm ! save the value of the original alm

! calculate cl - which corresponds to what is called (2l+1)Cl in equation 2 of Land & Magueijo 2005 (Axis of Evil paper). 
cl = abs(alm(1,lorder,0))**2
do m = 1, lorder
   cl = cl + 2.d0*(abs(alm(1,lorder,m))**2) ! sum_m 2*|alm|^2
enddo

do mvalue=0,lorder ! loop over m values to calculate statistics
  ! write(*,*) 'Findaxis: Calculating axis for l= ', lorder, ' m= ', mvalue!, lorder2
   do isim=0, nsim-1
      !      test = 0.0
      !   test2 = 0.0
      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 = origalm
      call rotate_alm(lmax, alm, phi,theta ,psi)!Rotate the map / alm's
      if (mvalue.eq.0) then
         test = 0.0
         test = (abs(alm(1,lorder, 0))**2) !sum_m=0 |alm|^2
      else
         test = 0.0
         test =  2.d0*(abs(alm(1,lorder,mvalue)))**2 !sum_m 2*|alm|^2
      endif
      stat(isim, mvalue) = test / cl
   enddo
enddo

maxstatloc = maxloc(stat) !find location of maximum quadrupole
imax = INT(maxstatloc(0)) !index of pixel
mmax_aoe = INT(maxstatloc(1))-1 !index of max m
rl = stat(imax, mmax_aoe)

! Special case for lorder =2 and m=1 - then prefer m=2
if (lorder .eq. 2) then 
   if (mmax_aoe .eq. 1) then 
write(*,*) 'lorder = 2 and mpreferred = 1'
write(*,*), 'lorder', lorder, 'mpref', mmax_aoe, 'rl', rl
maxstatloc = maxloc(stat(:, lorder))
imax = INT(maxstatloc(0))
mmax_aoe = lorder
rl = stat(imax, mmax_aoe)
write(*,*), 'lorder', lorder, 'mpref', mmax_aoe, 'rl', rl
endif
endif
!---------------------------------------------------
!For this study
call pix2vec_ring(nside,imax, vector)
call pix2ang_ring(nside,imax, thetamax, phimax)
!---------------------------------------------------
glmax = 180.d0 - phimax*rad2deg  
gbmax = 90.d0 - thetamax*rad2deg
END SUBROUTINE findaoe


! Taken from http://programming.ccp14.ac.uk/fortran-resources/~cgp/prof77/node3.html
! And then translated into Fortran 90
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
!  Real(sp), allocatable, dimension(:) :: X 
  !Allocate
!  allocate(x(1:npts))
  
  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