def get_parameter(parametrization_type):
    if parametrization_type == 0:  # normal parametrization
        # Get the current slider value
        t0 = t_value_input.value
    elif parametrization_type == 1:  # arc length parametrization
        f_x_str = x_component_input.value
        f_y_str = y_component_input.value
        f_x_sym = arc_functions.sym_parser(f_x_str)
        f_y_sym = arc_functions.sym_parser(f_y_str)
        from sympy.core import diff
        df_x_sym = diff(f_x_sym)
        df_y_sym = diff(f_y_sym)
        from sympy.abc import t
        from sympy import lambdify
        df_x = lambdify(t, df_x_sym, ['numpy'])
        df_y = lambdify(t, df_y_sym, ['numpy'])
        # compute arc length
        arc_length = arc_functions.arclength(df_x, df_y, arc_settings.t_value_max)
        # map input interval [t_value_min,t_value_max] to [0,arc_length]
        width_t = (arc_settings.t_value_max - arc_settings.t_value_min)
        t_fraction = (t_value_input.value - arc_settings.t_value_min) / width_t
        t_arc_length = t_fraction * arc_length
        # compute corresponding value on original parametrization
        t0 = arc_functions.s_inverse(df_x, df_y, t_arc_length)

    return t0

def lineseg_integrate(polygon, index, line_seg, expr, degree):
    """Helper function to compute the line integral of `expr` over `line_seg`

    Parameters
    ===========
    polygon : Face of a 3-Polytope
    index : index of line_seg in polygon
    line_seg : Line Segment

    >>> from sympy.integrals.intpoly import lineseg_integrate
    >>> polygon = [(0, 5, 0), (5, 5, 0), (5, 5, 5), (0, 5, 5)]
    >>> line_seg = [(0, 5, 0), (5, 5, 0)]
    >>> lineseg_integrate(polygon, 0, line_seg, 1, 0)
    5
    """
    if expr == S.Zero:
        return S.Zero
    result = S.Zero
    x0 = line_seg[0]
    distance = norm(tuple([line_seg[1][i] - line_seg[0][i] for i in
                           range(3)]))
    if isinstance(expr, Expr):
        expr_dict = {x: line_seg[1][0],
                     y: line_seg[1][1],
                     z: line_seg[1][2]}
        result += distance * expr.subs(expr_dict)
    else:
        result += distance * expr
    expr = diff(expr, x) * x0[0] + diff(expr, y) * x0[1] +\
        diff(expr, z) * x0[2]
    result += lineseg_integrate(polygon, index, line_seg, expr, degree - 1)
    result /= (degree + 1)
    return result

def calc_area(f_x_sym, f_y_sym, t_val):
    from sympy.abc import t

    f_x = lambdify(t, f_x_sym,['numpy'])
    f_y = lambdify(t, f_y_sym, ['numpy'])
    df_x = lambdify(t, diff(f_x_sym), ['numpy'])
    df_y = lambdify(t, diff(f_y_sym), ['numpy'])

    integrand = lambda tau: f_x(tau)*df_y(tau)-df_x(tau)*f_y(tau)
    return abs(0.5 * quad(integrand, 0, t_val)[0])

    def _eval_rewrite_as_Heaviside(self, *args, **kwargs):
        '''
        Rewrites a Singularity Function expression using Heavisides and DiracDeltas.

        '''
        x = self.args[0]
        a = self.args[1]
        n = sympify(self.args[2])

        if n == -2:
            return diff(Heaviside(x - a), x.free_symbols.pop(), 2)
        if n == -1:
            return diff(Heaviside(x - a), x.free_symbols.pop(), 1)
        if n.is_nonnegative:
            return (x - a)**n*Heaviside(x - a)

def integration_reduction(facets, index, a, b, expr, dims, degree):
    """Helper method for main_integrate. Returns the value of the input
    expression evaluated over the polytope facet referenced by a given index.

    Parameters
    ===========

    facets : List of facets of the polytope.
    index : Index referencing the facet to integrate the expression over.
    a : Hyperplane parameter denoting direction.
    b : Hyperplane parameter denoting distance.
    expr : The expression to integrate over the facet.
    dims : List of symbols denoting axes.
    degree : Degree of the homogeneous polynomial.

    Examples
    ========

    >>> from sympy.abc import x, y
    >>> from sympy.integrals.intpoly import integration_reduction,\
    hyperplane_parameters
    >>> from sympy.geometry.point import Point
    >>> from sympy.geometry.polygon import Polygon
    >>> triangle = Polygon(Point(0, 3), Point(5, 3), Point(1, 1))
    >>> facets = triangle.sides
    >>> a, b = hyperplane_parameters(triangle)[0]
    >>> integration_reduction(facets, 0, a, b, 1, (x, y), 0)
    5
    """
    if expr == S.Zero:
        return expr

    value = S.Zero
    x0 = facets[index].points[0]
    m = len(facets)
    gens = (x, y)

    inner_product = diff(expr, gens[0]) * x0[0] + diff(expr, gens[1]) * x0[1]

    if inner_product != 0:
        value += integration_reduction(facets, index, a, b,
                                       inner_product, dims, degree - 1)

    value += left_integral2D(m, index, facets, x0, expr, gens)

    return value/(len(dims) + degree - 1)

def polygon_integrate(facet, hp_param, index, facets, vertices, expr, degree):
    """Helper function to integrate the input uni/bi/trivariate polynomial
    over a certain face of the 3-Polytope.

    Parameters
    ===========

    facet : Particular face of the 3-Polytope over which `expr` is integrated
    index : The index of `facet` in `facets`
    facets : Faces of the 3-Polytope(expressed as indices of `vertices`)
    vertices : Vertices that constitute the facet
    expr : The input polynomial
    degree : Degree of `expr`

    Examples
    ========

    >>> from sympy.abc import x, y
    >>> from sympy.integrals.intpoly import polygon_integrate
    >>> cube = [[(0, 0, 0), (0, 0, 5), (0, 5, 0), (0, 5, 5), (5, 0, 0),\
                 (5, 0, 5), (5, 5, 0), (5, 5, 5)],\
                 [2, 6, 7, 3], [3, 7, 5, 1], [7, 6, 4, 5], [1, 5, 4, 0],\
                 [3, 1, 0, 2], [0, 4, 6, 2]]
    >>> facet = cube[1]
    >>> facets = cube[1:]
    >>> vertices = cube[0]
    >>> polygon_integrate(facet, [(0, 1, 0), 5], 0, facets, vertices, 1, 0)
    -25
    """
    expr = S(expr)
    if expr == S.Zero:
        return S.Zero
    result = S.Zero
    x0 = vertices[facet[0]]
    for i in range(len(facet)):
        side = (vertices[facet[i]], vertices[facet[(i + 1) % len(facet)]])
        result += distance_to_side(x0, side, hp_param[0]) *\
            lineseg_integrate(facet, i, side, expr, degree)
    if not expr.is_number:
        expr = diff(expr, x) * x0[0] + diff(expr, y) * x0[1] +\
            diff(expr, z) * x0[2]
        result += polygon_integrate(facet, hp_param, index, facets, vertices,
                                    expr, degree - 1)
    result /= (degree + 2)
    return result

 def _do(f, ab):
     dab_dsym = diff(ab, sym)
     if not dab_dsym:
         return S.Zero
     if isinstance(f, Integral):
         limits = [(x, x) if (len(l) == 1 and l[0] == x) else l
                   for l in f.limits]
         f = Integral(f.function, *limits)
     return f.subs(x, ab)*dab_dsym

    def _eval_derivative(self, sym):
        """Evaluate the derivative of the current Integral object.
        We follow these steps:

        (1) If sym is not part of the function nor the integration limits,
            return 0

        (2) Check for a possible application of the Fundamental Theorem of
            Calculus [1]

        (3) Derive under the integral sign [2]

        References:
           [1] http://en.wikipedia.org/wiki/Fundamental_theorem_of_calculus
           [2] http://en.wikipedia.org/wiki/Differentiation_under_the_integral_sign
        """

        if not sym in self.atoms(Symbol):
            return S.Zero

        if (sym, None) in self.limits:
            # case undefinite integral
            if len(self.limits) == 1:
                return self.function
            else:
                _limits = list(self.limits)
                _limits.pop(_limits.index((sym, None)))
                return Integral(self.function, *tuple(_limits))

        # diff under the integral sign
        # we do not check for regularity conditions (TODO), see issue 1116
        if len(self.limits) > 1:
            # TODO:implement the multidimensional case
            raise NotImplementedError
        int_var = self.limits[0][0]
        lower_limit, upper_limit = self.limits[0][1]
        if sym == int_var:
            sym = Symbol(str(int_var), dummy=True)
        return (
            self.function.subs(int_var, upper_limit) * diff(upper_limit, sym)
            - self.function.subs(int_var, lower_limit) * diff(lower_limit, sym)
            + integrate(diff(self.function, sym), (int_var, lower_limit, upper_limit))
        )

def test_sho_R_nl():
    omega, r = symbols('omega r')
    l = symbols('l', integer=True)
    u = Function('u')

    # check that it obeys the Schrodinger equation
    for n in range(5):
        schreq =  ( -diff(u(r), r, 2)/2 + ((l*(l+1))/(2*r**2)
                    + omega**2*r**2/2 - E_nl(n, l, omega))*u(r) )
        result = schreq.subs(u(r), r*R_nl(n, l, omega/2, r))
        assert simplify(result.doit()) == 0

def calculate_tangent(f_x_str, f_y_str, t0):
    f_x_sym = sym_parser(f_x_str)
    f_y_sym = sym_parser(f_y_str)
    f_x = parser(f_x_str)
    f_y = parser(f_y_str)

    from sympy.core import diff
    from sympy.abc import t

    df_x_sym = diff(f_x_sym)
    df_y_sym = diff(f_y_sym)
    df_x = lambdify(t, df_x_sym, ['numpy'])
    df_y = lambdify(t, df_y_sym, ['numpy'])

    x = np.array([f_x(t0)],dtype=np.float64)
    y = np.array([f_y(t0)],dtype=np.float64)
    u = np.array([df_x(t0)],dtype=np.float64)
    v = np.array([df_y(t0)],dtype=np.float64)

    return x,y,u,v

    def simplify(self, x):
        """simplify(self, x)

           Compute a simplified representation of the function using
           property number 4.

           x can be:

           - a symbol

           Examples
           ========

           >>> from sympy import DiracDelta
           >>> from sympy.abc import x, y

           >>> DiracDelta(x*y).simplify(x)
           DiracDelta(x)/Abs(y)
           >>> DiracDelta(x*y).simplify(y)
           DiracDelta(y)/Abs(x)

           >>> DiracDelta(x**2 + x - 2).simplify(x)
           DiracDelta(x - 1)/3 + DiracDelta(x + 2)/3

           See Also
           ========

           is_simple, Directdelta

        """
        from sympy.polys.polyroots import roots

        if not self.args[0].has(x) or (len(self.args) > 1 and self.args[1] != 0 ):
            return self
        try:
            argroots = roots(self.args[0], x)
            result = 0
            valid = True
            darg = abs(diff(self.args[0], x))
            for r, m in argroots.items():
                if r.is_real is not False and m == 1:
                    result += self.func(x - r)/darg.subs(x, r)
                else:
                    # don't handle non-real and if m != 1 then
                    # a polynomial will have a zero in the derivative (darg)
                    # at r
                    valid = False
                    break
            if valid:
                return result
        except PolynomialError:
            pass
        return self

    def simplify(self, x):
        """simplify(self, x)

           Compute a simplified representation of the function using
           property number 4.

           x can be:

           - a symbol

           Examples
           ========

           >>> from sympy import DiracDelta
           >>> from sympy.abc import x, y

           >>> DiracDelta(x*y).simplify(x)
           DiracDelta(x)/Abs(y)
           >>> DiracDelta(x*y).simplify(y)
           DiracDelta(y)/Abs(x)

           >>> DiracDelta(x**2 + x - 2).simplify(x)
           DiracDelta(x - 1)/3 + DiracDelta(x + 2)/3

           See Also
           ========

           is_simple, Directdelta

        """
        from sympy.polys.polyroots import roots

        if not self.args[0].has(x) or (len(self.args)>1 and self.args[1] != 0 ):
            return self
        try:
            argroots = roots(self.args[0], x, \
                                                     multiple=True)
            result = 0
            valid = True
            darg = diff(self.args[0], x)
            for r in argroots:
                #should I care about multiplicities of roots?
                if r.is_real and not darg.subs(x,r).is_zero:
                    result = result + DiracDelta(x - r)/abs(darg.subs(x,r))
                else:
                    valid = False
                    break
            if valid:
                return result
        except PolynomialError:
            pass
        return self

    def length(self):
        """The curve length.

        Examples
        ========

        >>> from sympy.geometry.curve import Curve
        >>> from sympy import cos, sin
        >>> from sympy.abc import t
        >>> Curve((t, t), (t, 0, 1)).length
        sqrt(2)
        """
        integrand = sqrt(sum(diff(func, self.limits[0])**2 for func in self.functions))
        return integrate(integrand, self.limits)

    def slope(self):
        """
        Returns a Singularity Function expression which represents
        the slope the elastic curve of the Beam object.

        Examples
        ========
        There is a beam of length 30 meters. A moment of magnitude 120 Nm is
        applied in the clockwise direction at the end of the beam. A pointload
        of magnitude 8 N is applied from the top of the beam at the starting
        point. There are two simple supports below the beam. One at the end
        and another one at a distance of 10 meters from the start. The
        deflection is restricted at both the supports.

        Using the sign convention of upward forces and clockwise moment
        being positive.

        >>> from sympy.physics.continuum_mechanics.beam import Beam
        >>> from sympy import symbols
        >>> E, I = symbols('E, I')
        >>> R1, R2 = symbols('R1, R2')
        >>> b = Beam(30, E, I)
        >>> b.apply_load(-8, 0, -1)
        >>> b.apply_load(R1, 10, -1)
        >>> b.apply_load(R2, 30, -1)
        >>> b.apply_load(120, 30, -2)
        >>> b.bc_deflection = [(10, 0), (30, 0)]
        >>> b.solve_for_reaction_loads(R1, R2)
        >>> b.slope()
        (-4*SingularityFunction(x, 0, 2) + 3*SingularityFunction(x, 10, 2)
            + 120*SingularityFunction(x, 30, 1) + SingularityFunction(x, 30, 2) + 4000/3)/(E*I)
        """
        x = self.variable
        E = self.elastic_modulus
        I = self.second_moment
        if not self._boundary_conditions['slope']:
            return diff(self.deflection(), x)

        C3 = Symbol('C3')
        slope_curve = integrate(self.bending_moment(), x) + C3

        bc_eqs = []
        for position, value in self._boundary_conditions['slope']:
            eqs = slope_curve.subs(x, position) - value
            bc_eqs.append(eqs)

        constants = list(linsolve(bc_eqs, C3))
        slope_curve = slope_curve.subs({C3: constants[0][0]})
        return S(1)/(E*I)*slope_curve

def line_integrate(field, curve, vars):
    """line_integrate(field, Curve, variables)

    Compute the line integral.

    Examples
    ========

    >>> from sympy import Curve, line_integrate, E, ln
    >>> from sympy.abc import x, y, t
    >>> C = Curve([E**t + 1, E**t - 1], (t, 0, ln(2)))
    >>> line_integrate(x + y, C, [x, y])
    3*sqrt(2)

    See Also
    ========

    integrate, Integral
    """
    F = sympify(field)
    if not F:
        raise ValueError(
            "Expecting function specifying field as first argument.")
    if not isinstance(curve, Curve):
        raise ValueError("Expecting Curve entity as second argument.")
    if not is_sequence(vars):
        raise ValueError("Expecting ordered iterable for variables.")
    if len(curve.functions) != len(vars):
        raise ValueError("Field variable size does not match curve dimension.")

    if curve.parameter in vars:
        raise ValueError("Curve parameter clashes with field parameters.")

    # Calculate derivatives for line parameter functions
    # F(r) -> F(r(t)) and finally F(r(t)*r'(t))
    Ft = F
    dldt = 0
    for i, var in enumerate(vars):
        _f = curve.functions[i]
        _dn = diff(_f, curve.parameter)
        # ...arc length
        dldt = dldt + (_dn * _dn)
        Ft = Ft.subs(var, _f)
    Ft = Ft * sqrt(dldt)

    integral = Integral(Ft, curve.limits).doit(deep=False)
    return integral

    def simplify(self, x):
        """simplify(self, x)

           Compute a simplified representation of the function using
           property number 4.

           x can be:

           - a symbol

           Examples
           --------

           >>> from sympy import *
           >>> x, y = symbols('xy')

           >>> DiracDelta(x*y).simplify(x)
           1/abs(y)*DiracDelta(x)
           >>> DiracDelta(x*y).simplify(y)
           1/abs(x)*DiracDelta(y)

           >>> DiracDelta(x**2+x-2).simplify(x)
           1/3*DiracDelta(-1 + x) + 1/3*DiracDelta(2 + x)

        """
        if not self.args[0].has(x) or (len(self.args)>1 and self.args[1] != 0 ):
            return self
        try:
            argroots = sympy.polys.rootfinding.roots(self.args[0],x, \
                                                     multiple=True)
            result = 0
            valid = True
            darg = diff(self.args[0], x)
            for r in argroots:
                #should I care about multiplicities of roots?
                if r.is_real and not darg.subs(x,r).is_zero:
                    result = result + DiracDelta(x - r)/abs(darg.subs(x,r))
                else:
                    valid = False
                    break
            if valid:
                return result
        except Exception,e :
            print e
            raise
            pass

def line_integrate(field, curve, vars):
    """line_integrate(field, Curve, variables)

       Compute the line integral.

       Examples
       --------
       >>> from sympy import *
       >>> x, y, t = symbols('xyt')
       >>> C = Curve([E**t + 1, E**t - 1], (t, 0, ln(2)))
       >>> line_integrate(x + y, C, [x, y])
       3*sqrt(2)

    """
    F = sympify(field)
    if not F:
        raise ValueError("Expecting function specifying field as first argument.")
    if not isinstance(curve, Curve):
        raise ValueError("Expecting Curve entity as second argument.")
    if not isinstance(vars, (list, tuple)):
        raise ValueError("Expecting list/tuple for variables.")
    if len(curve.functions) != len(vars):
        raise ValueError("Field variable size does not match curve dimension.")

    if curve.parameter in vars:
        raise ValueError("Curve parameter clashes with field parameters.")

    # Calculate derivatives for line parameter functions
    # F(r) -> F(r(t)) and finally F(r(t)*r'(t))
    Ft = F
    dldt = 0
    for i, var in enumerate(vars):
        _f = curve.functions[i]
        _dn = diff(_f, curve.parameter)
        # ...arc length
        dldt = dldt + (_dn * _dn)
        Ft = Ft.subs(var, _f)
    Ft = Ft * dldt**(S(1)/2)

    integral = Integral(Ft, curve.limits).doit()
    return integral

 def _eval_derivative(self, x):
     return self.func(*[(diff(e, x), c) for e, c in self.args])

def lie_w(omega, X, args):
    """Return a skew-symmetric tensor of type (0,p).
    Function lie_w calculates all the components of the Lie derivative
    differential forms in a symbolic form.

    Indexes the output tensor will start as well as the input tensor (array)
    "omega". If all the input parameters of the same type, they must be equal
    to the initial indexes.

    Examples:
    =========

    >>> from tensor_analysis.tensor_fields import lie_w
    >>> from sympy import symbols, cos
    >>> from tensor_analysis.arraypy import Arraypy
    >>> x1, x2, x3 = symbols('x1 x2 x3')

    omega - skew-symmetric tensor. Can be a tensor of type (0,p) or an array
    arraypy:

    >>> omega=Arraypy([2,3,1]).to_tensor((-1,-1))
    >>> omega[1,2]=x3
    >>> omega[1,3]=-x2
    >>> omega[2,1]=-x3
    >>> omega[2,3]=x1
    >>> omega[3,1]=x2
    >>> omega[3,2]=-x1

    X - the vector field along which the derivative is calculated:

    >>> X = [x1*x2**3,x2-cos(x3),x3**3-x1]

    args it's a list of symbol arguments.
    It's can be in list, array of arraypy or contravariant tensor:

    >>> arg = [x1, x2, x3]

    Lie derivative of a differential form:

    >>> li = lie_w(omega,X,arg)
    >>> print(li)
    0  x2**3*x3 + x3**3 + x3  -x2**4 - 3*x2*x3**2 - x2 + x3*sin(x3) + cos(x3)  
    -x2**3*x3 - x3**3 - x3  0  -2*x1*x2**3 + 3*x1*x3**2 + x1  
    x2**4 + 3*x2*x3**2 + x2 - x3*sin(x3) - cos(x3)  2*x1*x2**3 - 3*x1*x3**2 - x1  0

    >>> li.type_pq
    (0, 2)

    """
    # Handling of a vector of arguments
    check_vector_of_arguments(args)
    if isinstance(args, list):
        idx_args = 0
    else:
        idx_args = args.start_index[0]

    # Handling of a vector field
    check_the_vector_field(X)
    if isinstance(X, (TensorArray, Arraypy)):
        idx_X = X.start_index[0]
    else:
        idx_X = 0

    # Handling of a differential form
    if not isinstance(omega, (TensorArray, Arraypy)):
        raise ValueError(
            "The type of differential form must be TensorArray or Arraypy")
    if not is_asymmetric(omega):
        raise ValueError("The differential form must be a skew-symmetric")
    idx_omega = omega.start_index[0]


# Define the start index in the output tensor
    if type(omega) == type(X) == type(args):
        if idx_omega != idx_X or idx_omega != idx_args or idx_X != idx_args:
            raise ValueError(
                "The start index of differential form, vector field and \
                vector of argements must be equal")
    if isinstance(omega, type(X)) and idx_omega != idx_X:
        raise ValueError(
            "The start index of differential form and vector field must be \
            equal")
    idx_st = idx_omega

# Creating the output array in accordance with start indexes
    n = omega.shape[0]  # the dimensionality of the input array
    r = len(omega.shape)  # the rank of the input array
    a = Arraypy([r, n, idx_st])
    valence_list = [(-1) for k in range(r)]
    diff_Lie = a.to_tensor(valence_list)

    # Calculation
    idx = diff_Lie.start_index
    if isinstance(args, (TensorArray, Arraypy)):
        args = args.to_list()
    if isinstance(X, (TensorArray, Arraypy)):
        X = X.to_list()
    l_idx = len(idx)
    for p in range(len(diff_Lie)):
        for k in range(l_idx + 1):
#.........这里部分代码省略.........

def dw(omega, args):
    """Return a skew-symmetric tensor of type (0, p+1).
    Indexes the output tensor will start as well as the input tensor (array).
    If the input parameters of the same type, they must be equal to the initial
    indexes.

    Examples:
    =========

    >>> from tensor_analysis.tensor_fields import dw
    >>> from sympy import symbols
    >>> from tensor_analysis.arraypy import Arraypy
    >>> x1, x2, x3 = symbols('x1 x2 x3')

    omega - differential form, differential which is calculated.
    It's can be a skew-symmetric tensor of type (0,p) or an array arraypy:

    >>> omega=Arraypy([2,3,1]).to_tensor((-1,-1))
    >>> omega[1,2]=x3
    >>> omega[1,3]=-x2
    >>> omega[2,1]=-x3
    >>> omega[2,3]=x1
    >>> omega[3,1]=x2
    >>> omega[3,2]=-x1

    args it's a list of symbol arguments of differential form.
    It's can be in list, array of arraypy or contravariant tensor.

    External differential of a differential forms:

    >>> domega=dw(omega, [x1,x2,x3])
    >>> print(domega)
    0 0 0
    0 0 3
    0 -3 0
    0 0 -3
    0 0 0
    3 0 0
    0 3 0
    -3 0 0
    0 0 0
    >>> domega.type_pq
    (0, 3)

    """
    # Handling of a vector of arguments
    check_vector_of_arguments(args)
    # The definition of the start index of args
    if isinstance(args, list):
        idx_args = 0
    else:
        idx_args = args.start_index[0]

    # Handling of a differential form
    if not isinstance(omega, (TensorArray, Arraypy)):
        raise ValueError(
            "The type of differential form must be TensorArray or Arraypy")
    if omega.rank > 1:
        if not is_asymmetric(omega):
            raise ValueError("The differential form must be a skew-symmetric")
    idx_omega = omega.start_index[0]

    # Define the start index in the output tensor
    if isinstance(omega, type(args)) and idx_omega != idx_args:
        raise ValueError("The start index of the differential form and \
        vector of arguments must be equal")
    idx_st = idx_omega

    # Creating the output array in accordance with start indexes
    n = omega.shape[0]  # the dimensionality of the input array
    p = len(omega.shape)  # the rank of the input array
    valence_ind = [(-1) for k in range(p + 1)]
    d_omega = Arraypy([p + 1, n, idx_st]).to_tensor(valence_ind)

    # Calculation
    idx = d_omega.start_index
    if isinstance(args, (TensorArray, Arraypy)):
        args = args.to_list()

    for i in range(len(d_omega)):
        # tuple_list_indx it's list of tuple. Example:[(0, 1), (0, 1), (0, 0)]
        tuple_list_indx = [
            delete_index_from_list(
                idx,
                f) for f in range(
                len(idx))]
        for k in range(p + 1):
            d_omega[idx] += Add(((-1)**k) * diff(omega[tuple_list_indx[k]],
                                                 args[idx[k] - idx_st]))
        idx = d_omega.next_index(idx)

    # Output
    return d_omega