PROGRAM isw_anomalies
! Modified by Anais Rassat, Jan 2013: updated for public use
! Modified by Anais Rassat, Feb 2012: adapted code for general map as input
! Assume you input a 'clean' CMB map (e.g. Doppler already subtracted, isw map already subtracted or not...)
!=================================================
! Written by Anais Rassat, Feb 2011
! Purpose, find the preferred axis for l=2 and l=3
!=================================================
! gfortran anomalies_l2l3_axis.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_l2l3_axis.o

use healpix_types
use alm_tools
use pix_tools
use fitstools
use rngmod

integer(i4b) :: lmax
integer(i4b) :: lorder, lorder2
integer(I4B), parameter::nside=512, npix=12*nside**2

real(dp), allocatable, dimension(:,:) :: dw8
real(dp), dimension(2) :: costheta
real(dp) :: glmax2, gbmax2, glmax3, gbmax3
real(dp), dimension(3) :: vector2, vector3
real(dp) :: dotproduct

real(sp), allocatable, dimension(:,:) :: map, doppler
complex(spc), allocatable, dimension(:,:,:) :: alm, origalm

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

!type(planck_rng) :: rng_handle

!=================================
!Definitions and allocations
!=================================
lmax = 4
!DP
allocate(dw8(1:2*nside, 1:3))
dw8 = 1.0_dp
costheta =(-1.d0,1.d0)

!SP
allocate(doppler(0:nside2npix(nside)-1,1:1)) ! kinetic Doppler quadrupole map
allocate(map(0:nside2npix(nside)-1,1:3)) ! CMB map
allocate(alm(1:3, 0:lmax, 0:lmax))
allocate(origalm(1:3, 0:lmax, 0:lmax))

alm = 0.0
origalm = 0.0

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

!=================================
! Read in maps
!=================================
!--------------------------------
!Now chose which of the CMB maps you would like to test
!--------------------------------
call input_map('../kineticDoppler/kDopplerquad_map_n512.fits', doppler, npix, 1) ! NESTED order
call convert_nest2ring(nside, doppler) !dquad is now in ring order 

   call input_map('../CMB_LSS_products/W9_inpainted_sparsity_jan13mono_pows_100_l2l3.fits', map, npix, 1) ! NESTED
   call convert_nest2ring(nside, map)
   map = map / 1d6 / 2.725d0 ! Convert to unitless values
 map(:,1:1) = map(:,1:1) - doppler(:,:) !Both maps should be in RING 
!=================================
! Calculate Alm's of unrotated map
!=================================
call map2alm(nside, lmax, lmax, map, alm, costheta, dw8) !Calculte alm's
origalm = alm
!=================================
! Find the axis for l=2 and l=3 - output preferred axes
!=================================
lorder = 2 ! Quadrupole
lorder2 = 3 ! Octopole
call findaxis(alm, lmax , lorder, vector2, glmax2, gbmax2, lorder2, vector3, glmax3, gbmax3)
!=================================
! Calculate corresponding dotproduct
!=================================
dotproduct =vector2(1)*vector3(1) + vector2(2)*vector3(2) + vector2(3)*vector3(3)
!=================================
! PRINT OUT RESULTS
!=================================
write(*,*) '-------------------------'
write(*,*) '-------------------------'
write(*,*) 'Preferred axis for l=2'
write(*,*) glmax2  , gbmax2
write(*,*) '-------------------------'
write(*,*) '-------------------------'
write(*,*) 'Preferred axis for l=3'
write(*,*) glmax3 , gbmax3
write(*,*) '-------------------------'
write(*,*) '-------------------------'
write(*,*) 'n_2 . n_3,  angle (deg), Probability (%)'
write(*,*)  dotproduct, acos(dotproduct)*rad2deg, (1.d0 - dotproduct)*100.d0
write(*,*) '-------------------------'
write(*,*) '-------------------------'

END PROGRAM isw_anomalies

SUBROUTINE findaxis(alm, lmax, lorder, vector, glmax, gbmax, lorder2, vector2, glmax2, gbmax2)
use healpix_types
use alm_tools
use pix_tools
use fitstools
!------------------------------------------------------
! Written by Anais Rassat February 2011
! PURPOSE = find preferred axis for l=lorder and l=lorder2
!------------------------------------------------------
! INPUT
! alm              - alm of map to test
! lmax             - lmax must be greater than lorder and lorder2
! lorder           - first multipole for which to search preferred axis
! lorder2          - second multipole for which to search preferred axis
!------------------------------------------------------
! OUTPUT
! vector           - preferred axis for l=lorder
! (glmax, gbmax)   - corresponding Galactic coordinates
! vector2          - preferred axis for l=lorder2
! (glmax2, gbmax2) - corresponding Galactic coordinates for l=lorder2
!------------------------------------------------------

! Define variables
integer(i4b) :: lmax, nsim, isim, imax, imax2
integer(I4B), parameter::nside=512, npix=12*nside**2
! to speed up can put nside < 512, nside=512 gives you a precision of < 1' for the preferred axis (check in Oliveira paper)
integer(I4B) :: lorder, morder, lorder2

real(sp) :: test, 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
allocate(stat(0:nsim-1))
allocate(stat2(0:nsim-1))
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
!write(*,*) 'Findaxis: Calculating axis for l= ', lorder, 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
!   call rotate_alm(lmax, alm, phi,theta ,psi)!Rotate the map / alm's
!   write(*,*) 'Calculating for ', lorder
   do morder=1,lorder
      test =  test + 2.d0* morder**2*(abs(alm(1,lorder,morder)))**2 !sum_m m^2 |alm|^2
   enddo
!   write(*,*) 'Calculating for ', lorder2
   do morder=1,lorder2
      test2 =  test2 + 2.d0* morder**2*(abs(alm(1,lorder2,morder)))**2 !sum_m m^2 |alm|^2
   enddo
   stat(isim) = test
   stat2(isim) = test2
enddo
maxstatloc = maxloc(stat) !find location of maximum quadrupole
maxstatloc2 = maxloc(stat2) !find location of maximum octopole
imax = INT(maxstatloc(0)) !index of max quad
imax2 = INT(maxstatloc2(0)) !index of max quad
write(*,*) INT(maxstatloc(0:3))
write(*,*) INT(maxstatloc2(0:3))
!---------------------------------------------------
!For this study
call pix2vec_ring(nside,imax, vector)
call pix2ang_ring(nside,imax, thetamax, phimax)
call pix2vec_ring(nside,imax2, vector2)
call pix2ang_ring(nside,imax2, thetamax2, phimax2)
!---------------------------------------------------

glmax = 180.d0 - phimax*rad2deg  
gbmax = 90.d0 - thetamax*rad2deg


glmax2 = 180.d0 - phimax2*rad2deg  
gbmax2 = 90.d0 - thetamax2*rad2deg


END SUBROUTINE findaxis