! 
! Copyright (c) 2008 Kristopher L. Kuhlman (klkuhlm at sandia dot gov)
! 
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
! 
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
! 
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! 
! Malama, B., K.L. Kuhlman, and W. Barrash, 2007. Semi-analytical solution
!  for flow in leaky unconfined aquifer-aquitard systems, Journal of
!  Hydrology, 346(1-2), 59–68.
! http://dx.doi.org/10.1016/j.jhydrol.2007.08.018
!
! Malama, B., K.L. Kuhlman, and W. Barrash, 2008. Semi-analytical solution
!  for flow in a leaky unconfined aquifer toward a partially penetrating
!  pumping well, Journal of Hydrology, 356(1-2), 234–244.
! http://dx.doi.org/10.1016/j.jhydrol.2008.03.029
!

module gauss_points

  implicit none
  private
  public :: Gauss_Quad_Setup

  contains

  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  subroutine Gauss_Quad_Setup(intOrd,coord,weight)
    use constants, only : I4B, EP

    integer(I4B), intent(in) :: intOrd    ! Gauss integration order
    real(EP), intent(out), dimension(:) :: coord,weight

    ! from Stroud & Secrest "Gaussian Quadrature Formulas", 1966: Table 1
    ! compared to quad-precision output from GAULEG() subroutine
    ! in Antia numerical methods book (weights and abscissa agree to >= 30 places)
    select case (intOrd)
    case(1)
       coord(1) = 0.0_EP
       weight(1) = 2.0_EP
    case(2)
       coord(1) = 0.577350269189625764509148780502_EP
       coord(2) = -coord(1)
       weight(1:2) = 1.0_EP
    case(3)
       coord(1:2) = (/ 0.774596669241483377035853079956_EP, 0.0_EP /)
       coord(3) = -coord(1)
       weight(1:2) = (/ 0.55555555555555555555555555556_EP, &
                      & 0.88888888888888888888888888889_EP /)
       weight(3) = weight(1)
    case(4)
       coord(1:2) = (/ 0.861136311594052575223946488893_EP, &
                     & 0.339981043584856264802665759103_EP /)
       coord(3:4) = -coord(2:1:-1)
       weight(1:2) = (/ 0.347854845137453857373063949222_EP, &
                      & 0.652145154862546142626936050778_EP /)
       weight(3:4) = weight(2:1:-1)
    case(5)
       coord(1:3) = (/ 0.906179845938663992797626878299_EP, &
                     & 0.538469310105683091036314420700_EP, 0.0_EP /)
       coord(4:5) = -coord(2:1:-1)
       weight(1:3) = (/ 0.236926885056189087514264040720_EP, &
                      & 0.478628670499366468041291514836_EP, &
                     &  0.568888888888888888888888888889_EP /)
       weight(4:5) = weight(2:1:-1)
    case(6)
       coord(1:3) = (/ 0.932469514203152027812301554494_EP, &
                     & 0.661209386466264513661399595020_EP, &
                     & 0.238619186083196908630501721681_EP /)
       coord(4:6) = -coord(3:1:-1)
       weight(1:3) = (/ 0.171324492379170345040296142173_EP, &
                      & 0.360761573048138607569833513838_EP, &
                     &  0.467913934572691047389870343990_EP /)
       weight(4:6) = weight(3:1:-1)
    case(7)
       coord(1:4) = (/ 0.949107912342758524526189684048_EP, &
                     & 0.741531185599394439863864773281_EP, &
                     & 0.405845151377397166906606412077_EP, 0.0_EP /)
       coord(5:7) = -coord(3:1:-1)
       weight(1:4) = (/ 0.129484966168869693270611432679_EP, &
                      & 0.279705391489276667901467771424_EP, &
                      & 0.381830050505118944950369775489_EP, &
                      & 0.417959183673469387755102040816_EP /)
       weight(5:7) = weight(3:1:-1)
    case(8)
       coord(1:4) = (/ 0.960289856497536231683560868569_EP, &
                     & 0.796666477413626739591553936476_EP, &
                     & 0.525532409916328985817739049189_EP, &
                     & 0.183434642495649804939476142360_EP /)
       coord(5:8) = -coord(4:1:-1)
       weight(1:4) = (/ 0.101228536290376259152531354310_EP, &
                      & 0.222381034453374470544355994426_EP, &
                      & 0.313706645877887287337962201987_EP, &
                      & 0.362683783378361982965150449277_EP /)
       weight(5:8) = weight(4:1:-1)
    case(9)
       coord(1:5) = (/ 0.968160239507626089835576202904_EP, &
                     & 0.836031107326635794299429788070_EP, &
                     & 0.613371432700590397308702039341_EP, &
                     & 0.324253423403808929038538014643_EP, 0.0_EP /)
       coord(6:9) = -coord(4:1:-1)
       weight(1:5) = (/ 0.0812743883615744119718921581105_EP, &
                      & 0.180648160694857404058472031243_EP, &
                      & 0.260610696402935462318742869419_EP, &
                      & 0.312347077040002840068630406584_EP, &
                      & 0.330239355001259763164525069287_EP /)
       weight(6:9) = weight(4:1:-1)
    case(10)
       coord(1:5) = (/ 0.973906528517171720077964012084_EP, &
                     & 0.865063366688984510732096688423_EP, &
                     & 0.679409568299024406234327365115_EP, &
                     & 0.433395394129247190799265943166_EP, &
                     & 0.148874338981631210884826001130_EP /)
       coord(6:10) = -coord(5:1:-1)
       weight(1:5) = (/ 0.0666713443086881375935688098933_EP, &
                      & 0.149451349150580593145776339658_EP, &
                      & 0.219086362515982043995534934228_EP, &
                      & 0.269266719309996355091226921569_EP, &
                      & 0.295524224714752870173892994651_EP /)
       weight(6:10) = weight(5:1:-1)
    case(12) 
       coord(1:6) = (/ 0.981560634246719250690549090149_EP, &
                     & 0.904117256370474856678465866119_EP, &
                     & 0.769902674194304687036893833213_EP, &
                     & 0.587317954286617447296702418941_EP, &
                     & 0.367831498998180193752691536644_EP, &
                     & 0.125233408511468915472441369464_EP /)
       coord(7:12) = -coord(6:1:-1)
       weight(1:6) = (/ 0.0471753363865118271946159614850_EP, &
                      & 0.106939325995318430960254718194_EP, &
                      & 0.160078328543346226334652529543_EP, &
                      & 0.203167426723065921749064455810_EP, &
                      & 0.233492536538354808760849898925_EP, &
                      & 0.249147045813402785000562436043_EP /)
       weight(7:12) = weight(6:1:-1)
    case(20) 
       coord(1:10) = (/ 0.993128599185094924786122388471_EP, &
                      & 0.963971927277913791267666131197_EP, &
                      & 0.912234428251325905867752441203_EP, &
                      & 0.839116971822218823394529061702_EP, &
                      & 0.746331906460150792614305070356_EP, &
                      & 0.636053680726515025452836696226_EP, &
                      & 0.510867001950827098004364050955_EP, &
                      & 0.373706088715419560672548177025_EP, &
                      & 0.227785851141645078080496195369_EP, &
                      & 0.0765265211334973337546404093988_EP /)
       coord(11:20) = -coord(10:1:-1)
       weight(1:10) = (/ 0.0176140071391521183118619623519_EP, &
                       & 0.0406014298003869413310399522749_EP, &
                       & 0.0626720483341090635695065351870_EP, &
                       & 0.0832767415767047487247581432220_EP, &
                       & 0.101930119817240435036750135480_EP, &
                       & 0.118194531961518417312377377711_EP, &
                       & 0.131688638449176626898494499748_EP, &
                       & 0.142096109318382051329298325067_EP, &
                       & 0.149172986472603746787828737002_EP, &
                       & 0.152753387130725850698084331955_EP /)
       weight(11:20) = weight(10:1:-1)

    case default
       print *, 'unsupported Gaussian Quadrature integration order:', intOrd
       stop
    end select

  end subroutine Gauss_Quad_Setup

end module gauss_points
