ssht_sampling_mod.f90

Go to the documentation of this file.

! SSHT package to perform spin spherical harmonic transforms
! Copyright (C) 2011  Jason McEwen
! See LICENSE.txt for license details


!------------------------------------------------------------------------------
! ssht_sampling_mod  -- SSHT library sampling class
! 

module ssht_sampling_mod

  use ssht_types_mod
  use ssht_error_mod

  implicit none

  private


  !---------------------------------------
  ! Subroutine and function scope
  !---------------------------------------

  public :: &
       ssht_sampling_elm2ind, &
       ssht_sampling_ind2elm, &
       ssht_sampling_dh_t2theta, &
       ssht_sampling_dh_p2phi, &
       ssht_sampling_gl_p2phi, &
       ssht_sampling_mw_t2theta, &
       ssht_sampling_mw_p2phi, &
       ssht_sampling_mweo_t2theta, &
       ssht_sampling_mweo_p2phi, &
       ssht_sampling_weight_dh, &
       ssht_sampling_weight_mw, &
       ssht_sampling_gl_thetas_weights


  !---------------------------------------
  ! Global variables
  !---------------------------------------

  integer, public, parameter :: SSHT_SAMPLING_DH = 1

  integer, public, parameter :: SSHT_SAMPLING_GL = 2

  integer, public, parameter :: SSHT_SAMPLING_MW = 3

  integer, public, parameter :: SSHT_SAMPLING_MWEO = 4

  integer, public, parameter :: SSHT_SAMPLING_DEFAULT = SSHT_SAMPLING_MW


  !---------------------------------------
  ! Data types
  !---------------------------------------

  ! None.


  !----------------------------------------------------------------------------

contains


  !============================================================================
  ! Sampling weights
  !============================================================================


  !--------------------------------------------------------------------------
  ! ssht_sampling_weight_dh
  !

  function ssht_sampling_weight_dh(theta_t, L) result(w)

    real(dp), intent(in) :: theta_t
    integer, intent(in) :: L
    real(dp) :: w       

    integer :: k

    w = 0d0
    do k = 0,L-1
       w = w + sin((2d0*k+1d0)*theta_t) / real(2d0*k+1d0,dp)
    end do
    w = (2d0/real(L,dp)) * sin(theta_t) * w

  end function ssht_sampling_weight_dh


  !--------------------------------------------------------------------------
  ! ssht_sampling_weight_mw
  !

  function ssht_sampling_weight_mw(p) result(w)

    integer, intent(in) :: p
    complex(dpc) :: w

    if(p == 1) then
       w = I * PION2
    elseif(p == -1) then
       w = - I * PION2
    elseif(mod(p,2) == 0 ) then
       ! Even case
       w = 2d0 / (1d0 - p**2)
    else
       ! Odd case (|p| /= 1)
       w = 0d0
    end if

  end function ssht_sampling_weight_mw


  !--------------------------------------------------------------------------
  ! ssht_sampling_gl_thetas_weights
  !

  subroutine ssht_sampling_gl_thetas_weights(thetas, weights, L)

    integer, intent(in) :: L
    real(dp), intent(out) :: thetas(0:L-1)
    real(dp), intent(out) :: weights(0:L-1)
    
    real(dp) :: znodes(0:L-1)
    integer :: p

    call gauleg(-1d0, 1d0, znodes(0:L-1), weights(0:L-1), L)

    do p = 0, L-1
       thetas(p) = acos(znodes(p))
    end do

  end subroutine ssht_sampling_gl_thetas_weights


  !--------------------------------------------------------------------------
  ! gauleg
  !

  subroutine gauleg(x1, x2, x, w, n)

    integer, intent(in) :: n
    real(dp), intent(in) :: x1,x2
    real(dp), intent(out) :: x(1:n),w(1:n)

    real(dp), parameter :: EPS = 1d-14
    integer :: i,j,m
    real(dp) :: p1,p2,p3,pp,xl,xm,z,z1

    m=(n+1)/2
    xm=0.5d0*(x2+x1)
    xl=0.5d0*(x2-x1)
    do i=1,m
       z=cos(3.141592654d0*(i-.25d0)/(n+.5d0))
       do 
          p1=1.d0
          p2=0.d0
          do j=1,n
             p3=p2
             p2=p1
             p1=((2.d0*j-1.d0)*z*p2-(j-1.d0)*p3)/j
          end do
          pp=n*(z*p1-p2)/(z*z-1.d0)
          z1=z
          z=z1-p1/pp
          if(abs(z-z1) <= EPS) exit
       end do
       x(i)=xm-xl*z
       x(n+1-i)=xm+xl*z
       w(i)=2.d0*xl/((1.d0-z*z)*pp*pp)
       w(n+1-i)=w(i)
    end do

  end subroutine gauleg


  !============================================================================
  ! Sampling relations
  !============================================================================


  !----------------------------------------------------------------------------
  ! ssht_sampling_dh_t2theta
  !
  
  function ssht_sampling_dh_t2theta(t, L) result (theta)

    integer, intent(in) :: t
    integer, intent(in) :: L
    real(dp) :: theta

    theta = (2d0*t+1d0)*PI / (4d0*L)

  end function ssht_sampling_dh_t2theta


  !----------------------------------------------------------------------------
  ! ssht_sampling_dh_p2phi
  !
  
  function ssht_sampling_dh_p2phi(p, L) result (phi)

    integer, intent(in) :: p
    integer, intent(in) :: L
    real(dp) :: phi

    phi = 2d0*p*PI / (2d0*L - 1d0)      

  end function ssht_sampling_dh_p2phi


  !----------------------------------------------------------------------------
  ! ssht_sampling_gl_p2phi
  !
  
  function ssht_sampling_gl_p2phi(p, L) result (phi)

    integer, intent(in) :: p
    integer, intent(in) :: L
    real(dp) :: phi

    phi = 2d0*p*PI / (2d0*L - 1d0)      

  end function ssht_sampling_gl_p2phi


  !----------------------------------------------------------------------------
  ! ssht_sampling_mw_t2theta
  !
 
  function ssht_sampling_mw_t2theta(t, L) result (theta)

    integer, intent(in) :: t
    integer, intent(in) :: L
    real(dp) :: theta

    theta = (2d0*t+1d0)*PI / (2d0*L - 1d0)

  end function ssht_sampling_mw_t2theta


  !----------------------------------------------------------------------------
  ! ssht_sampling_mw_p2phi
  !

  function ssht_sampling_mw_p2phi(p, L) result (phi)

    integer, intent(in) :: p
    integer, intent(in) :: L
    real(dp) :: phi

    phi = 2d0*p*PI / (2d0*L - 1d0)      

  end function ssht_sampling_mw_p2phi


  !----------------------------------------------------------------------------
  ! ssht_sampling_mweo_t2theta
  !

  function ssht_sampling_mweo_t2theta(t, L) result (theta)

    integer, intent(in) :: t
    integer, intent(in) :: L
    real(dp) :: theta

    theta = (2d0*t+1d0)*PI / (2d0*L - 1d0)

  end function ssht_sampling_mweo_t2theta


  !----------------------------------------------------------------------------
  ! ssht_sampling_mweo_p2phi
  !

  function ssht_sampling_mweo_p2phi(p, L) result (phi)

    integer, intent(in) :: p
    integer, intent(in) :: L
    real(dp) :: phi

    phi = (2d0*p+1d0)*PI / (2d0*L - 1d0)

  end function ssht_sampling_mweo_p2phi


  !============================================================================
  ! Harmonic index relations
  !============================================================================

  
  !----------------------------------------------------------------------------
  ! ssht_sampling_elm2ind
  !

  subroutine ssht_sampling_elm2ind(ind, el, m)

    integer, intent(out) :: ind
    integer, intent(in) :: el, m

    ind = el**2 + el + m

  end subroutine ssht_sampling_elm2ind


  !----------------------------------------------------------------------------
  ! ssht_sampling_ind2elm
  !

  subroutine ssht_sampling_ind2elm(el, m, ind)

    integer, intent(out) :: el, m
    integer, intent(in) :: ind

    el = floor(sqrt(real(ind,dp)))
    m = ind - el**2 - el

  end subroutine ssht_sampling_ind2elm


end module ssht_sampling_mod