diff --git a/fortran/sharp.f90 b/fortran/sharp.f90 index 7dc815d..fb80c3c 100644 --- a/fortran/sharp.f90 +++ b/fortran/sharp.f90 @@ -103,12 +103,32 @@ module sharp type(c_ptr), intent(in) :: alm(*), map(*) end subroutine c_sharp_execute_mpi + ! Legendre transforms + subroutine c_sharp_legendre_transform(bl, recfac, lmax, x, out, nx) & + bind(c, name='sharp_legendre_transform') + use iso_c_binding + integer(c_ptrdiff_t), value :: lmax, nx + real(c_double) :: bl(lmax + 1), x(nx), out(nx) + real(c_double), optional :: recfac(lmax + 1) + end subroutine c_sharp_legendre_transform + + subroutine c_sharp_legendre_transform_s(bl, recfac, lmax, x, out, nx) & + bind(c, name='sharp_legendre_transform_s') + use iso_c_binding + integer(c_ptrdiff_t), value :: lmax, nx + real(c_float) :: bl(lmax + 1), x(nx), out(nx) + real(c_float), optional :: recfac(lmax + 1) + end subroutine c_sharp_legendre_transform_s end interface interface sharp_execute module procedure sharp_execute_d end interface + interface sharp_legendre_transform + module procedure sharp_legendre_transform_d, sharp_legendre_transform_s + end interface sharp_legendre_transform + contains ! alm info @@ -240,6 +260,25 @@ contains end if end subroutine sharp_execute_d + subroutine sharp_legendre_transform_d(bl, x, out) + use iso_c_binding + real(c_double) :: bl(:) + real(c_double) :: x(:), out(size(x)) + !-- + integer(c_ptrdiff_t) :: lmax, nx + call c_sharp_legendre_transform(bl, lmax=int(size(bl) - 1, c_ptrdiff_t), & + x=x, out=out, nx=int(size(x), c_ptrdiff_t)) + end subroutine sharp_legendre_transform_d + + subroutine sharp_legendre_transform_s(bl, x, out) + use iso_c_binding + real(c_float) :: bl(:) + real(c_float) :: x(:), out(size(x)) + !-- + integer(c_ptrdiff_t) :: lmax, nx + call c_sharp_legendre_transform_s(bl, lmax=int(size(bl) - 1, c_ptrdiff_t), & + x=x, out=out, nx=int(size(x), c_ptrdiff_t)) + end subroutine sharp_legendre_transform_s end module diff --git a/fortran/test_sharp.f90 b/fortran/test_sharp.f90 index f42cf0c..0b7cce2 100644 --- a/fortran/test_sharp.f90 +++ b/fortran/test_sharp.f90 @@ -54,4 +54,31 @@ program test_sharp print *, 'DONE' call MPI_Finalize(ierr) + print *, 'LEGENDRE TRANSFORMS' + + call test_legendre_transforms() + +contains + subroutine test_legendre_transforms() + integer, parameter :: lmax = 20, nx=10 + real(c_double) :: bl(0:lmax) + real(c_double) :: x(nx), out(nx) + real(c_float) :: out_s(nx) + !-- + integer :: l, i + + do l = 0, lmax + bl(l) = 1.0 / real(l + 1, c_double) + end do + do i = 1, nx + x(i) = 1 / real(i, c_double) + end do + out = 0 + call sharp_legendre_transform(bl, x, out) + print *, out + call sharp_legendre_transform(real(bl, c_float), real(x, c_float), out_s) + print *, out_s + end subroutine test_legendre_transforms + + end program test_sharp