PROGRAM isw_anomalies
  ! Modified Feb 2012 - cleaned version
  ! Written by Anais Rassat, June 2011
  ! Test for octopole planarity
  
  ! TO COMPILE 
  ! From 7/7/2011 Use the following: 
  ! gfortran anomalies_octplan.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_octplan.o
  
  
  use healpix_types
  use alm_tools
  use pix_tools
  use fitstools, ONLY: input_map, output_map, write_asctab
  use rngmod!, only: rand_init, rng_handle
  
!===============================================================
! Begin - numbers to change if you want to speed up calculation
! Number of GRF simulations to obtain t_stat probability distribution
  integer(i4b) :: nsimmax=1000
! Note for calculation to be included in a paper you should use nsimmax = 1d4
! However to speed up calculation can reduce this number
! End - numbers to change if you want to speed up calculation
!===============================================================

  integer(I4B), parameter::nside=512, npix=12*nside**2
  integer(i4b) :: lmax, nmaps=1  
  real(dp), allocatable, dimension(:,:) :: dw8
  real(dp), dimension(2) :: costheta
  real(dp) :: avg, sd
  
  real(sp) :: t_stat,  probt! de Oliveira-Costa Octopole planarity test
  real(sp), allocatable, dimension(:) :: stat
  
  real(sp), allocatable, dimension(:,:) :: map, doppler
  complex(spc), allocatable, dimension(:,:,:) :: alm, origalm
  
  integer, parameter::nlheader=80
  character(len=80), dimension (1:nlheader) ::header
  character(len=80), allocatable, dimension(:) :: cmbnames
  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(stat(0:nsimmax-1)) ! For t statistic
  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))
  allocate(cmbnames(1:nmaps))
  header(:) = ' '
  alm = 0.0
  origalm = 0.0
  
  write(*,*) '-------------------------'
  write(*,*) '-------------------------'
  write(*,*) '-------------------------'
  
  !=================================
  ! Read in maps
  !=================================
  !--------------------------------
  ! OPEN kinetic DOPPLER MAP
  !--------------------------------
  call input_map('../kineticDoppler/kDopplerquad_map_n512.fits', doppler, npix, 1) ! NESTED order
  call convert_nest2ring(nside, doppler) !dquad is now in ring order 
  !--------------------------------
  ! OPEN CMB MAP
  !--------------------------------
  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 / 1d3 / 2.725d0 ! Convert to unitless values
  write(*,*) 'CMB map read'
  !--------------------------------
  !SUBTRACT kinetic DOPPLER MAP
  !--------------------------------
  map(:,1:1) = map(:,1:1) - doppler(:,:) !Both maps should be in RING 
  !=================================
  ! Calculate Alm's of unrotated map
  !=================================
  alm = 0.d0
  call map2alm(nside, lmax, lmax, map, alm, costheta, dw8) !Calculte alm's
  origalm = alm
  !=================================
  ! Calculate t-stat
  !=================================
  call findt_stat(alm, lmax, t_stat)
  write(*,*) ' t_stat: ', t_stat                    
  !===================
  ! Calculate probability
  !===================
  call init_calctprob(alm,lmax,nsimmax,stat, input)            
  call calctprob(stat, nsimmax,t_stat, probt)
  write(*,*) 'Prob (%)', probt
  
  
END PROGRAM isw_anomalies


SUBROUTINE findt_stat(alm,lmax, test)
use healpix_types
use alm_tools
use pix_tools
use fitstools
! Written by Anais Rassat, June 2011
! Calculates t statistic defined Equation (3) of de Oliveira-Costa et al 2004.
! http://adsabs.harvard.edu/abs/2004PhRvD..69f3516D
! Define variables

integer(i4b) :: lmax, nsim, isim, imax
integer(I4B), parameter::nside=128, npix=12*nside**2

real(sp) :: test, test2
real(sp), allocatable, dimension(:) :: stat
real(sp), allocatable, dimension(:) :: maxstatloc

real(dp) :: psi, theta, phi

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(origalm(1:3, 0:lmax, 0:lmax))
allocate(maxstatloc(0:nsim-1))
maxstatloc = 0.0
origalm = 0.0 
origalm = alm ! save the value of the original alm
do isim=0, nsim-1
   test = 0.0
   call pix2ang_ring(nside,isim,theta,phi)!calculate theta,phi corresponding to pixel isim
   psi = 1.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
   test2 = abs(alm(1,3,0))**2 + 2.d0*abs(alm(1,3,1))**2 + 2.d0*abs(alm(1,3,2))**2 + 2.d0*abs(alm(1,3,3))**2
   test =  2.d0*abs(alm(1,3,3))**2/test2
   stat(isim) = test
enddo
maxstatloc = maxloc(stat) !find location of maximum quadrupole
imax = INT(maxstatloc(0)) !index of max quad
!write(*,*) imax 
!write(*,*) stat(imax)
test = stat(imax)
!write(*,*) INT(maxstatloc(0:3))
alm = origalm
END SUBROUTINE findt_stat

SUBROUTINE init_calctprob(alm,lmax,nsimmax, stat, input)            
use healpix_types
use alm_tools
use pix_tools
use fitstools
use rngmod!, only: rand_init, rng_handle
! Modified by Anais Rassat Feb 2012 - Fixed bug in m=0 alm for GRF (was an extra sqrt(2) that has now been removed)
! Written by Anais Rassat, June 2011
! Calculates probability t statistic defined Equation (3) of de Oliveira-Costa et al 2004 for GRF simulations
! http://adsabs.harvard.edu/abs/2004PhRvD..69f3516D
! Define variables
integer(i4b) :: lmax, isim!, nsim, isim, imax
integer(i4b), parameter::nlheader=80
integer(i4b) :: nsimmax ! Number of l=3 simulations
character(len=80), dimension (1:nlheader) ::header
real(sp), allocatable, dimension(:,:) :: clmap
type(planck_rng) :: rng_handle
integer(I4B), parameter::nside=512, npix=12*nside**2
real(sp) :: t_stat, countt, probt! de Oliveira-Costa Octopole planarity test
real(sp), dimension(0:nsimmax-1) :: stat
complex(spc), dimension(1:3, 0:lmax, 0:lmax) :: alm, alm_grf
!real(sp) :: x, y
character(len=80):: cmbname, cmbnameshort
integer(i4b) :: input,l,m
character(len=7) :: nsimmaxnumber
character(len=3) :: inputname
character(len=50) :: filename

! Allocate variables
!allocate(stat(0:nsimmax-1)) ! For t statistic
allocate(clmap(0:lmax, 1:3))
write(inputname, '(I3)') input
write(nsimmaxnumber, '(I7)') nsimmax    ! this converts the simulation number to a string
nsimmaxnumber = adjustl(nsimmaxnumber) ! removes spaces ?
inputname = adjustl(inputname)
cmbnameshort = trim(cmbname)
!filename = 'tstat_noISW'//trim(inputname)//'_'//trim(nsimmaxnumber)//'.dat'!cmbnameshort!//trim(nsimmaxnumber)//'.dat'
!write(*,*) filename
!!write(*,*) nsimmaxnumber
!!write(*,*) input, inputname
!!stop
!!write(*,*) cmbname
!open (4, file=filename, status='replace')
call alm2cl(lmax, lmax, alm, clmap)
header = ''
!call write_asctab(clmap, lmax, 1, header, nlheader, '!cl.fits')
call rand_init(rng_handle, -1)
!write(*,*) rng_handle
!stop
!=============================================
! Begin - GRFs simulations with only octopole
do isim = 0, nsimmax-1  
write(*,*) isim
   alm_grf = alm
   do l=3,3 !Only the m=3 terms are taken as GRF
      do m=0,l
         if (m .ne. 0) then 
            alm_grf(1, l, m) = CMPLX(rand_gauss(rng_handle)*sqrt(clmap(l,1)/2.d0), rand_gauss(rng_handle)*sqrt(clmap(l,1)/2.d0))
         endif
         if (m .eq. 0) then 
            alm_grf(1, l, m) = CMPLX(rand_gauss(rng_handle)*sqrt(clmap(l,1)), 0.d0)
         endif
      enddo
   enddo
   
   call findt_stat(alm_grf, lmax, t_stat) ! Calculate t_stat for each simulation
   stat(isim) = t_stat
!   write(4,*) t_stat
!   call flush(4)
enddo
! End - GRFs simulations with only octopole
!=============================================
!close(4)
END SUBROUTINE init_calctprob

SUBROUTINE calctprob(stat, nsimmax,t_stat, probt)
use healpix_types
use alm_tools
use pix_tools
use fitstools
use rngmod!, only: rand_init, rng_handle
! Written by Anais Rassat, June 2011
! For a given 't' value (see de Oliveira-Costa et al 2004), 
! And a pre-initialised stat, calculates the probability probt 
! init_calctprob must be run beforehand
integer(i4b) :: nsimmax ! Number of l=3 simulations
integer(i4b) :: isim
real(sp) :: t_stat, countt, probt! de Oliveira-Costa Octopole planarity test
real(sp), dimension(0:nsimmax-1) :: stat

!allocate(stat(0:nsimmax-1)) ! For t statistic

countt = 0.d0 
do isim = 0, nsimmax-1 
   if (stat(isim) .gt. t_stat) then 
      countt = countt + 1.d0
   endif
enddo
write(*,*), t_stat
            write(*,*) 'Count t stat' ,countt
probt = countt / size(stat)*100.d0 ! probability in %

END SUBROUTINE calctprob


! 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) 
!     write(*,*) sum, x(i)
  END DO
!  write(*,*) avg, sum/npts, sum, npts
  AVG = SUM / NPTS 
  
!Now calculate Standard deviation
  DO I = 1, NPTS
     SUMSQ = SUMSQ + (X(I)-AVG)**2 
!write(*,*) sumsq, (x(i)-avg)**2, npts, sqrt(sumsq/(npts-1))
  END DO


!  WRITE(*,*) SUM, NPTS
!  WRITE(*,*) SUMSQ, AVG, NPTS-1

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