!>
!> @file disrhs.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
SUBROUTINE disrhs(mbess, spl, rhs)
!
!   Assembly the RHS using 2d spline spl
!
  USE bsplines
  IMPLICIT NONE
  INTEGER, INTENT(in) :: mbess
  TYPE(spline2d), INTENT(in) :: spl
  DOUBLE PRECISION, INTENT(out) :: rhs(:)
  INTEGER :: n1, nidbas1, ndim1, ng1
  INTEGER :: n2, nidbas2, ndim2, ng2
  INTEGER :: i, j, ig1, ig2, k1, k2, i1, j2, ij, nrank
  DOUBLE PRECISION, ALLOCATABLE :: xg1(:), wg1(:), fun1(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: xg2(:), wg2(:), fun2(:,:)
  DOUBLE PRECISION:: contrib
!===========================================================================
!              1.0 Prologue
!
!   Properties of spline space
!
  CALL get_dim(spl%sp1, ndim1, n1, nidbas1)
  CALL get_dim(spl%sp2, ndim2, n2, nidbas2) 
!
  ALLOCATE(fun1(0:nidbas1,1)) ! needs only basis functions (no derivatives)
  ALLOCATE(fun2(0:nidbas2,1)) ! needs only basis functions (no derivatives)
!
!   Gauss quadature
!
  CALL get_gauss(spl%sp1, ng1)
  CALL get_gauss(spl%sp2, ng2)
  WRITE(*,'(/a, 2i3)') 'Gauss points and weights, ngauss =', ng1, ng2
  ALLOCATE(xg1(ng1), wg1(ng1))
  ALLOCATE(xg2(ng1), wg2(ng1))
!===========================================================================
!              2.0 Assembly loop
!
  nrank = SIZE(rhs)
  rhs(1:nrank) = 0.0d0
!
  DO i=1,n1
     CALL get_gauss(spl%sp1, ng1, i, xg1, wg1)
     DO ig1=1,ng1
        CALL basfun(xg1(ig1), spl%sp1, fun1, i)
        DO j=1,n2
           CALL get_gauss(spl%sp2, ng2, j, xg2, wg2)
           DO ig2=1,ng2
              CALL basfun(xg2(ig2), spl%sp2, fun2, j)
              contrib = wg1(ig1)*wg2(ig2) * &
                   &    rhseq(xg1(ig1),xg2(ig2), mbess)
              DO k1=0,nidbas1
                 i1 = i+k1
                 DO k2=0,nidbas2
                    j2 = MODULO(j+k2-1,n2) + 1
                    ij = j2 + (i1-1)*n2
                    rhs(ij) = rhs(ij) + contrib*fun1(k1,1)*fun2(k2,1)
                 END DO
              END DO
           END DO
        END DO
     END DO
  END DO
!===========================================================================
!              9.0  Epilogue
!
  DEALLOCATE(xg1, wg1, fun1)
  DEALLOCATE(xg2, wg2, fun2)
!
CONTAINS
  DOUBLE PRECISION FUNCTION rhseq(x1, x2, m)
    DOUBLE PRECISION, INTENT(in) :: x1, x2
    INTEGER, INTENT(in) :: m
    rhseq = REAL(4*(m+1),8)*x1**(m+1)*COS(REAL(m,8)*x2)
  END FUNCTION rhseq
END SUBROUTINE disrhs
!
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE disrhs3(mbess, npow, spl, rhs)
!
!   Assembly the RHS using 3d spline spl
!
  USE bsplines
  IMPLICIT NONE
  INTEGER, INTENT(in)           :: mbess, npow
  TYPE(spline2d1d), TARGET      :: spl
  DOUBLE PRECISION, INTENT(out) :: rhs(:,:)
!
  TYPE(spline1d), POINTER :: sp1, sp2, sp3
  INTEGER :: n1, nidbas1, ndim1, ng1
  INTEGER :: n2, nidbas2, ndim2, ng2
  INTEGER :: n3, nidbas3, ndim3, ng3
  INTEGER :: i, j, k, ig1, ig2, ig3, k1, k2, k3, i1, j2, ij, kk, nrank
  DOUBLE PRECISION, ALLOCATABLE :: xg1(:), wg1(:), fun1(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: xg2(:), wg2(:), fun2(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: xg3(:), wg3(:), fun3(:,:)
  DOUBLE PRECISION:: contrib
!===========================================================================
!              1.0 Prologue
!
!   Properties of spline space
!
  sp1 => spl%sp12%sp1
  sp2 => spl%sp12%sp2
  sp3 => spl%sp3
!
  CALL get_dim(sp1, ndim1, n1, nidbas1)
  CALL get_dim(sp2, ndim2, n2, nidbas2) 
  CALL get_dim(sp3, ndim3, n3, nidbas3) 
!
  ALLOCATE(fun1(0:nidbas1,1)) ! needs only basis functions (no derivatives)
  ALLOCATE(fun2(0:nidbas2,1)) ! needs only basis functions (no derivatives)
  ALLOCATE(fun3(0:nidbas3,1)) ! needs only basis functions (no derivatives)
!
!   Gauss quadature
!
  CALL get_gauss(sp1, ng1)
  CALL get_gauss(sp2, ng2)
  CALL get_gauss(sp3, ng3)
  WRITE(*,'(/a, 3i3)') 'Gauss points and weights, ngauss =', ng1, ng2, ng3
  ALLOCATE(xg1(ng1), wg1(ng1))
  ALLOCATE(xg2(ng2), wg2(ng2))
  ALLOCATE(xg3(ng3), wg3(ng3))
!===========================================================================
!              2.0 Assembly loop
!
  nrank = SIZE(rhs,1)
  rhs(1:nrank,1:n3) = 0.0d0
!
  DO i=1,n1
     CALL get_gauss(sp1, ng1, i, xg1, wg1)
     DO ig1=1,ng1
        CALL basfun(xg1(ig1), sp1, fun1, i)
        DO j=1,n2
           CALL get_gauss(sp2, ng2, j, xg2, wg2)
           DO ig2=1,ng2
              CALL basfun(xg2(ig2), sp2, fun2, j)
              DO k=1,n3
                 CALL get_gauss(sp3, ng3, k, xg3, wg3)
                 DO ig3=1,ng3
                    CALL basfun(xg3(ig3), sp3, fun3, k)
                    contrib = wg1(ig1)*wg2(ig2)*wg3(ig3) * &
                         &    rhseq(xg1(ig1), xg2(ig2), xg3(ig3), mbess, npow)
                    DO k1=0,nidbas1
                       i1 = i+k1
                       DO k2=0,nidbas2
                          j2 = MODULO(j+k2-1,n2) + 1
                          ij = j2 + (i1-1)*n2
                          DO k3=0,nidbas3
                             kk = MODULO(k+k3-1,n3) + 1
                             rhs(ij,kk) = rhs(ij, kk) + &
                                  &  contrib*fun1(k1,1)*fun2(k2,1)*fun3(k3,1)
                          END DO
                       END DO
                    END DO
                 END DO
              END DO
           END DO
        END DO
     END DO
  END DO
!===========================================================================
!              9.0  Epilogue
!
  DEALLOCATE(xg1, wg1, fun1)
  DEALLOCATE(xg2, wg2, fun2)
  DEALLOCATE(xg3, wg3, fun3)
!
CONTAINS
  DOUBLE PRECISION FUNCTION rhseq(x1, x2, x3, m, n)
    DOUBLE PRECISION, INTENT(in) :: x1, x2, x3
    INTEGER, INTENT(in) :: m, n
    rhseq = REAL(4*(m+1),8)*x1**(m+1)*COS(REAL(m,8)*x2)*COS(x3)**n
  END FUNCTION rhseq
END SUBROUTINE disrhs3