gauss_legendre.F90

Go to the documentation of this file.

!! Copyright (C) 2006 Hyllios
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program 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 General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module gauss_legendre_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
 
  implicit none
 
  private
  public :: &
    gauss_legendre_points
 
 
  ! 2 point
  real(real64), target :: GL_points_2(1) = (/                                   &
    -0.57735026919_real64 /)
  real(real64), target :: GL_weights_2(1) = (/                                  &
    m_one /)
 
  ! 3 point
  real(real64), target :: GL_points_3(2) = (/                                   &
    -0.774596669241_real64, m_zero /)
  real(real64), target :: GL_weights_3(2) = (/                                  &
    0.555555555556_real64, 0.888888888889_real64 /)
 
  ! 4 point
  real(real64), target :: GL_points_4(2) = (/                                   &
    -0.861136311594_real64, -0.339981043585_real64 /)
  real(real64), target :: GL_weights_4(2) = (/                                  &
    0.347854845137_real64, 0.652145154863_real64 /)
 
  ! 5 point
  real(real64), target :: GL_points_5(3) = (/                                   &
    -0.906179845939_real64, -0.538469310106_real64, m_zero /)
  real(real64), target :: GL_weights_5(3) = (/                                  &
    0.236926885056_real64, 0.478628670499_real64, 0.568888888889_real64 /)
 
  ! 6 point
  real(real64), target :: GL_points_6(3) = (/                                   &
    -0.932469514203_real64, -0.661209386466_real64, -0.238619186083_real64 /)
  real(real64), target :: GL_weights_6(3) = (/                                  &
    0.171324492379_real64, 0.360761573048_real64, 0.467913934573_real64 /)
 
  ! 7 point
  real(real64), target :: GL_points_7(4) = (/                                   &
    -0.949107912343_real64, -0.741531185599_real64, -0.405845151377_real64, &
    m_zero /)
  real(real64), target :: GL_weights_7(4) = (/                                  &
    0.129484966169_real64, 0.279705391489_real64, 0.381830050505_real64,    &
    0.417959183673_real64 /)
 
  ! 8 point
  real(real64), target :: GL_points_8(4) = (/                                   &
    -0.960289856498_real64, -0.796666477414_real64, -0.525532409916_real64, &
    -0.183434642496_real64 /)
  real(real64), target :: GL_weights_8(4) = (/                                  &
    0.10122853629_real64, 0.222381034453_real64, 0.313706645878_real64,     &
    0.362683783378_real64 /)
 
  ! 9 point
  real(real64), target :: GL_points_9(5) = (/                                   &
    -0.968160239508_real64, -0.836031107327_real64, -0.613371432701_real64, &
    -0.324253423404_real64, m_zero /)
  real(real64), target :: GL_weights_9(5) = (/                                  &
    0.0812743883616_real64, 0.180648160695_real64, 0.260610696403_real64,   &
    0.31234707704_real64, 0.330239355001_real64 /)
 
  ! 10 point
  real(real64), target :: GL_points_10(5) = (/                                  &
    -0.973906528517_real64, -0.865063366689_real64, -0.679409568299_real64, &
    -0.433395394129_real64, -0.148874338982_real64 /)
  real(real64), target :: GL_weights_10(5) = (/                                 &
    0.0666713443087_real64, 0.149451349151_real64, 0.219086362516_real64,   &
    0.26926671931_real64,   0.295524224715_real64 /)
 
  ! 11 point
  real(real64), target :: GL_points_11(6) = (/                                  &
    -0.978228658146_real64, -0.887062599768_real64, -0.730152005574_real64, &
    -0.519096129207_real64, -0.269543155952_real64, m_zero /)
  real(real64), target :: GL_weights_11(6) = (/                                 &
    0.0556685671162_real64, 0.125580369465_real64, 0.186290210928_real64,   &
    0.233193764592_real64,  0.26280454451_real64,  0.272925086778_real64 /)
 
  ! 12 point
  real(real64), target :: GL_points_12(6) = (/                                  &
    -0.981560634247_real64, -0.90411725637_real64,  -0.769902674194_real64, &
    -0.587317954287_real64, -0.367831498998_real64, -0.125233408511_real64 /)
  real(real64), target :: GL_weights_12(6) = (/                                 &
    0.0471753363866_real64, 0.106939325995_real64, 0.160078328543_real64,   &
    0.203167426723_real64,  0.233492536538_real64, 0.249147045813_real64 /)
 
contains
 
  ! --------------------------------------------------------------------
  subroutine gauss_legendre_points(n, points, weights)
    integer, intent(in)  :: n
    real(real64),   intent(out) :: points(:)
    real(real64),   intent(out) :: weights(:)
 
    integer :: i
 
    real(real64), pointer :: points_ref(:), weights_ref(:)
 
    nullify(points_ref)
    nullify(weights_ref)
 
    assert(n >= 2.and.n <= 12)
    select case (n)
    case (2)
      points_ref  => gl_points_2
      weights_ref => gl_weights_2
    case (3)
      points_ref  => gl_points_3
      weights_ref => gl_weights_3
    case (4)
      points_ref  => gl_points_4
      weights_ref => gl_weights_4
    case (5)
      points_ref  => gl_points_5
      weights_ref => gl_weights_5
    case (6)
      points_ref  => gl_points_6
      weights_ref => gl_weights_6
    case (7)
      points_ref  => gl_points_7
      weights_ref => gl_weights_7
    case (8)
      points_ref  => gl_points_8
      weights_ref => gl_weights_8
    case (9)
      points_ref  => gl_points_9
      weights_ref => gl_weights_9
    case (10)
      points_ref  => gl_points_10
      weights_ref => gl_weights_10
    case (11)
      points_ref  => gl_points_11
      weights_ref => gl_weights_11
    case (12)
      points_ref  => gl_points_12
      weights_ref => gl_weights_12
    end select
 
    do i = 1, n/2
      points(i)     =  points_ref(i)
      points(n-i+1) = -points_ref(i)
      weights(i)     =  weights_ref(i)
      weights(n-i+1) =  weights_ref(i)
    end do
    if (mod(n, 2) == 1) then
      points(n/2 + 1) = points_ref(n/2 + 1)
      weights(n/2 + 1) = weights_ref(n/2 + 1)
    end if
 
  end subroutine gauss_legendre_points
 
end module gauss_legendre_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: