! Compute Power Spectrum and Power Spectrum Density
!                        of 
! Spherical Harmonic Function
! by Y.Ishihara @ RISE Project, NAOJ
program SHPOW
      use SHTOOLS

      implicit none
      integer, parameter ::   degmax = 3600
      real*8 ::               cilm(2,degmax+1, degmax+1), spectrum1(degmax+1), spectrum2(degmax+1)
      integer ::              lmax,  i
      character*80 ::         shfile
      
      write(*,*) 'Input file name of Spherical Harmonics'
      read(*,*) shfile
      call SHRead(shfile, cilm, lmax)
      write(*,*) 'Max. Degree of Input Spherical Harmonics is ...', lmax
      call SHPowerSpectrum(cilm, lmax, spectrum1)
      call SHPowerSpectrumDensity(cilm, lmax, spectrum2)
      open (100, file='shpow.out',status='unknown')
      do i=1, lmax+1
        write(100,*) i-1, spectrum1(i), spectrum2(i)
      enddo
      close (100)
      stop
end program  SHPOW
