import numpy as np

def computeQuadrature(integrationRule,numIntegrPoints):
    #Compute integration points and weights for numerical integration between -1 and 1
    if integrationRule=='Gauss-Lobatto':
        integrPoints,integrWeights=computeQuadrature_GaussLobatto(numIntegrPoints)

    #Compute integration points and weights for numerical integration between 0 and 1
    integrPoints=(integrPoints+1)/2
    integrWeights=integrWeights/2

    return integrPoints,integrWeights


def computeQuadrature_GaussLobatto(numIntegrPoints):
    #integration between -1 and -1
    if numIntegrPoints==3:
        integrPoints=np.array([-1,0,1])
        integrWeights=np.array([1/3,4/3,1/3])
    elif numIntegrPoints==4:
        integrPoints=np.array([-1,-1/np.sqrt(5),1/np.sqrt(5),1])
        integrWeights=np.array([1/6,5/6,5/6,1/6])
    elif numIntegrPoints==5:
        integrPoints=np.array([-1,-np.sqrt(3/7),0,np.sqrt(3/7),1])
        integrWeights=np.array([1/10,49/90,32/45,49/90,1/10])
    elif numIntegrPoints==6:
        integrPoints=np.array([-1, -0.765055,-0.285232,0.285232, 0.765055, 1])
        integrWeights=np.array([0.066667, 0.378475, 0.554858, 0.554858, 0.378475, 0.066667])

    return integrPoints,integrWeights