text: from casadi import *
  
  
  def construct_polynomials_basis(d, poly_type):
  
      # Get collocation points
      tau_root = np.append(0, collocation_points(d, poly_type))
  
      # Coefficients of the collocation equation
      C = np.zeros((d + 1, d + 1))
  
      # Coefficients of the continuity equation
      D = np.zeros(d + 1)
  
      # Coefficients of the quadrature function
      B = np.zeros(d + 1)
  
      # Construct polynomial basis
      for j in range(d + 1):
          # Construct Lagrange polynomials to get the polynomial basis at the collocation point
          p = np.poly1d([1])
          for r in range(d + 1):
              if r != j:
                  p *= np.poly1d([1, -tau_root[r]]) / (tau_root[j] - tau_root[r])
  
          # Evaluate the polynomial at the final time to get the coefficients of the continuity equation
          D[j] = p(1.0)
          # Evaluate the time derivative of the polynomial at all collocation points to get the coefficients of the continuity
          # equation
          pder = np.polyder(p)
          for r in range(d + 1):
              C[j, r] = pder(tau_root[r])
  
          # Evaluate the integral of the polynomial to get the coefficients of the quadrature function
          pint = np.polyint(p)
          B[j] = pint(1.0)
  
      return C, D, B
  
  