-- Gkyl ------------------------------------------------------------------------
--
-- Rules for Gaussian quadrature in 1D. Multiple dimensional integrals
-- can be computed by tensor products.
--
--    _______     ___
-- + 6 @ |||| # P ||| +
--------------------------------------------------------------------------------

local ffi  = require "ffi"
local xsys = require "xsys"
local new, copy, fill, sizeof, typeof, metatype = xsys.from(ffi,
     "new, copy, fill, sizeof, typeof, metatype")

-- Ordinates for Gaussian quadrature. The first entry is just a buffer
-- to allow 1-based indexing
local ordinates = {
   -- N=1
   new("double[2]", 0, 0.0),
   -- N=2
   new("double[3]", 0, -0.5773502691896257645091, 0.5773502691896257645091),
   -- N=3
   new("double[4]", 0, -0.7745966692414833770359, 0, 0.7745966692414833770359),
   -- N=4
   new("double[5]", 0, -0.8611363115940525752239, -0.3399810435848562648027, 0.3399810435848562648027, 0.8611363115940525752239),
   -- N=5
   new("double[6]", 0, -0.9061798459386639927976, -0.5384693101056830910363, 0, 0.5384693101056830910363, 0.9061798459386639927976),
   -- N=6
   new("double[7]", 0, -0.9324695142031520278123, -0.6612093864662645136614, -0.2386191860831969086305, 0.2386191860831969086305, 0.6612093864662645136614, 0.9324695142031520278123),
   -- N=7
   new("double[8]", 0, -0.9491079123427585245262, -0.7415311855993944398639, -0.4058451513773971669066, 0, 0.4058451513773971669066, 0.7415311855993944398639, 0.9491079123427585245262),
   -- N=8
   new("double[9]", 0, -0.960289856497536231684, -0.7966664774136267395916, -0.5255324099163289858177, -0.1834346424956498049395, 0.1834346424956498049395, 0.525532409916328985818, 0.796666477413626739592, 0.9602898564975362316836),
}

-- Weights for Gaussian quadrature. The first entry is just a buffer
-- to allow 1-based indexing
local weights = {
   -- N=1
   new("double[2]", 0, 2.0),
   -- N=2
   new("double[3]", 0, 1.0, 1.0),
   -- N=3
   new("double[4]", 0, 0.5555555555555555555556, 0.888888888888888888889, 0.555555555555555555556),
   -- N=4
   new("double[5]", 0, 0.3478548451374538573731, 0.6521451548625461426269, 0.652145154862546142627, 0.3478548451374538573731),
   -- N=5
   new("double[6]", 0, 0.2369268850561890875143, 0.4786286704993664680413, 0.568888888888888888889, 0.478628670499366468041, 0.236926885056189087514),
   -- N=6
   new("double[7]", 0, 0.1713244923791703450403, 0.36076157304813860757, 0.46791393457269104739, 0.46791393457269104739, 0.36076157304813860757, 0.1713244923791703450403),
   -- N=7
   new("double[8]", 0, 0.1294849661688696932706, 0.279705391489276667901, 0.38183005050511894495, 0.4179591836734693877552, 0.38183005050511894495, 0.279705391489276667901, 0.1294849661688696932706),
   -- N=8
   new("double[9]", 0, 0.1012285362903762591525, 0.222381034453374470544, 0.313706645877887287338, 0.3626837833783619829652, 0.3626837833783619829652, 0.31370664587788728734, 0.222381034453374470544, 0.1012285362903762591525),
}

return {
   ordinates = ordinates,   
   weights = weights,
}

