bianchi2_lut_mod.f90

Go to the documentation of this file.

!------------------------------------------------------------------------------
! bianchi2_lut_mod
!
!------------------------------------------------------------------------------

module bianchi2_lut_mod

  use s2_types_mod
  use bianchi2_error_mod

  implicit none

  private


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

  public :: bianchi2_lut_init, &
            bianchi2_lut_write, &
            bianchi2_lut_read, &
            bianchi2_lut_access, &
            bianchi2_lut_free, &
            bianchi2_lut_plgndr


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

  type, public :: bianchi2_lut
     private
     integer :: lmax_LUT = 0

     integer :: Ntheta_LUT = 0

     real(s2_dp), dimension(:,:), allocatable :: table_LUT
  end type bianchi2_lut


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

  contains


    !-------------------------------------------------------------------------
    ! bianchi2_lut_init
    !
    !-------------------------------------------------------------------------

    function bianchi2_lut_init(lmax,Ntheta) result(lut)
      
      integer, intent(in) :: lmax, ntheta
      type(bianchi2_lut) :: lut

      integer :: l, itheta
      real(s2_dp) :: theta,dtheta

      lut%lmax_LUT = lmax
      lut%Ntheta_LUT = ntheta

      allocate(lut%table_LUT(1:lmax,0:ntheta-1))

      dtheta = (pi - 0d0) / Real(ntheta, s2_dp)
      do l=1, lmax
         do itheta=0, ntheta-1
            theta = itheta*dtheta
            ! Compute the data with plgndr.
            lut%table_LUT(l,itheta) = bianchi2_lut_plgndr(l,1,cos(theta))
            theta = theta + dtheta
         end do
      end do

    end function bianchi2_lut_init


    !-------------------------------------------------------------------------
    ! bianchi2_lut_write
    !
    !-------------------------------------------------------------------------

    subroutine bianchi2_lut_write(lut,filename_LUT)

      type(bianchi2_lut), intent(in) :: lut
      character(len=S2_STRING_LEN), intent(in) :: filename_lut

      integer :: l,itheta

      open(10,file=trim(filename_lut),form='unformatted',&
           status='replace',action='write')

      ! Write the sizes of the table.
      write(10) lut%lmax_LUT
      write(10) lut%Ntheta_LUT

      ! Write the table
      do l=1, lut%lmax_LUT
         do itheta=0, (lut%Ntheta_LUT - 1)
            write(10) lut%table_LUT(l,itheta)
         end do
      end do

      close(10)

    end subroutine bianchi2_lut_write


    !-------------------------------------------------------------------------
    ! bianchi2_lut_read
    !
    !-------------------------------------------------------------------------

    function bianchi2_lut_read(lmax,Ntheta,filename_LUT) result(lut)
      
      integer, intent(in) :: lmax,ntheta
      character(len=S2_STRING_LEN) :: filename_lut
      type(bianchi2_lut) :: lut

      integer :: l,itheta

      open(10,file=trim(filename_lut),form='unformatted',&
              status='old',action='read')

      read(10) lut%lmax_LUT
      read(10) lut%Ntheta_LUT
      allocate(lut%table_LUT(1:lmax,0:ntheta-1))

      if(lmax .ne. lut%lmax_LUT) then
         call bianchi2_error(bianchi2_error_plm1table_l_invalid, &
                 'bianchi2_plm1table_getval')
  
      elseif(ntheta .ne. lut%Ntheta_LUT) then
         call bianchi2_error(bianchi2_error_plm1table_theta_invalid, &
              'bianchi2_plm1table_getval')

      else     
         do l=1, lmax      
            do itheta=0, ntheta-1
               read(10) lut%table_LUT(l,itheta)
            end do
         end do
      end if

      close(10)
      
    end function bianchi2_lut_read


    !--------------------------------------------------------------------------
    ! bianchi2_lut_access
    !
    !--------------------------------------------------------------------------

    function bianchi2_lut_access(lut,l,itheta) result(plm1)
      
      type(bianchi2_lut), intent(in) :: lut
      integer, intent(in) :: l, itheta
      real(s2_dp) :: plm1

      plm1 = lut%table_LUT(l,itheta)

    end function bianchi2_lut_access

    !--------------------------------------------------------------------------
    ! bianchi2_lut_free
    !
    !--------------------------------------------------------------------------

    subroutine bianchi2_lut_free(lut)
      
      type(bianchi2_lut), intent(inout) :: lut

      if(allocated(lut%table_LUT)) deallocate(lut%table_LUT)
      lut%lmax_LUT = 0
      lut%Ntheta_LUT = 0

    end subroutine bianchi2_lut_free


    !--------------------------------------------------------------------------
    ! bianchi2_lut_plgndr
    !
    !--------------------------------------------------------------------------

    function bianchi2_lut_plgndr(l,m,x) result(plgndr)

      integer :: l, m
      real(s2_dp) :: plgndr, x

      integer :: i, ll
      real(s2_dp) ::  fact, pll, pmm, pmmp1, somx2

      if(m<0 .or. m>l .or. abs(x)>1d0) stop 'bad arguments in plgndr'

      pmm=1.                     ! Compute Pmm.
      if(m.gt.0) then
         somx2=sqrt((1.-x)*(1.+x))
         fact=1.
         do i=1,m
            pmm=-pmm*fact*somx2
            fact=fact+2.
         enddo
      endif
      if(l.eq.m) then
         plgndr=pmm
      else
         pmmp1=x*(2*m+1)*pmm      ! Compute Pm m+1.
         if(l.eq.m+1) then
            plgndr=pmmp1
         else                    ! Compute Pm  l , l > m+ 1.
            do ll=m+2,l
               pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m)
               pmm=pmmp1
               pmmp1=pll
            enddo
            plgndr=pll
         endif
      endif
      return

    end function bianchi2_lut_plgndr


end module bianchi2_lut_mod