def _eval_as_leading_term(self, x):
        from sympy import expand_mul, factor_terms

        old = self

        self = expand_mul(self)
        if not self.is_Add:
            return self.as_leading_term(x)

        unbounded = [t for t in self.args if t.is_unbounded]

        self = self.func(*[t.as_leading_term(x) for t in self.args]).removeO()
        if not self:
            # simple leading term analysis gave us 0 but we have to send
            # back a term, so compute the leading term (via series)
            return old.compute_leading_term(x)
        elif self is S.NaN:
            return old.func._from_args(unbounded)
        elif not self.is_Add:
            return self
        else:
            plain = self.func(*[s for s, _ in self.extract_leading_order(x)])
            rv = factor_terms(plain, fraction=False)
            rv_fraction = factor_terms(rv, fraction=True)
            # if it simplifies to an x-free expression, return that;
            # tests don't fail if we don't but it seems nicer to do this
            if x not in rv_fraction.free_symbols:
                if rv_fraction.is_zero and plain.is_zero is not True:
                    return (self - plain)._eval_as_leading_term(x)
                return rv_fraction
            return rv

    def __new__(cls, b, e, evaluate=True):
        from sympy.functions.elementary.exponential import exp_polar
        from sympy.functions import log
        # don't optimize "if e==0; return 1" here; it's better to handle that
        # in the calling routine so this doesn't get called
        b = _sympify(b)
        e = _sympify(e)
        if evaluate:
            if e is S.Zero:
                return S.One
            elif e is S.One:
                return b
            elif S.NaN in (b, e):
                if b is S.One:  # already handled e == 0 above
                    return S.One
                return S.NaN
            else:
                if e.func == log:
                    if len(e.args) == 2:
                        lbase = e.args[1]
                    else:
                        lbase = S.Exp1
                    if lbase == b:
                        return e.args[0]

                if e is Mul and e.args[1].func == log:
                    if len(e.args[1].args) == 2:
                        lbase = e.args[1].args[1]
                    else:
                        lbase = S.Exp1
                    if lbase == b:
                        return e.args[1].args[0] ** e.args[0]

                # recognize base as E
                if not e.is_Atom and b is not S.Exp1 and b.func is not exp_polar:
                    from sympy import numer, denom, log, sign, im, factor_terms
                    c, ex = factor_terms(e, sign=False).as_coeff_Mul()
                    den = denom(ex)
                    if den.func is log and den.args[0] == b:
                        return S.Exp1**(c*numer(ex))
                    elif den.is_Add:
                        s = sign(im(b))
                        if s.is_Number and s and den == \
                                log(-factor_terms(b, sign=False)) + s*S.ImaginaryUnit*S.Pi:
                            return S.Exp1**(c*numer(ex))

                obj = b._eval_power(e)
                if obj is not None:
                    return obj
        obj = Expr.__new__(cls, b, e)
        obj.is_commutative = (b.is_commutative and e.is_commutative)
        return obj

    def _eval_as_leading_term(self, x):
        from sympy import expand_mul, factor_terms

        old = self

        expr = expand_mul(self)
        if not expr.is_Add:
            return expr.as_leading_term(x)

        infinite = [t for t in expr.args if t.is_infinite]

        expr = expr.func(*[t.as_leading_term(x) for t in expr.args]).removeO()
        if not expr:
            # simple leading term analysis gave us 0 but we have to send
            # back a term, so compute the leading term (via series)
            return old.compute_leading_term(x)
        elif expr is S.NaN:
            return old.func._from_args(infinite)
        elif not expr.is_Add:
            return expr
        else:
            plain = expr.func(*[s for s, _ in expr.extract_leading_order(x)])
            rv = factor_terms(plain, fraction=False)
            rv_simplify = rv.simplify()
            # if it simplifies to an x-free expression, return that;
            # tests don't fail if we don't but it seems nicer to do this
            if x not in rv_simplify.free_symbols:
                if rv_simplify.is_zero and plain.is_zero is not True:
                    return (expr - plain)._eval_as_leading_term(x)
                return rv_simplify
            return rv

    def contact_forces(self):
        """Returns contact forces on each wheel."""

        self.conForceNoncontri = self._kane.auxiliary_eqs.applyfunc(
                            lambda w: factor_terms(signsimp(w))).subs(self._kdd)

        return self.conForceNoncontri

    def __new__(cls, b, e, evaluate=None):
        if evaluate is None:
            evaluate = global_evaluate[0]
        from sympy.functions.elementary.exponential import exp_polar

        b = _sympify(b)
        e = _sympify(e)
        if evaluate:
            if e is S.Zero:
                return S.One
            elif e is S.One:
                return b
            elif e.is_integer and _coeff_isneg(b):
                if e.is_even:
                    b = -b
                elif e.is_odd:
                    return -Pow(-b, e)
            if b is S.One:
                if e in (S.NaN, S.Infinity, -S.Infinity):
                    return S.NaN
                return S.One
            elif S.NaN in (b, e):  # XXX S.NaN**x -> S.NaN under assumption that x != 0
                return S.NaN
            else:
                # recognize base as E
                if not e.is_Atom and b is not S.Exp1 and b.func is not exp_polar:
                    from sympy import numer, denom, log, sign, im, factor_terms
                    c, ex = factor_terms(e, sign=False).as_coeff_Mul()
                    den = denom(ex)
                    if den.func is log and den.args[0] == b:
                        return S.Exp1**(c*numer(ex))
                    elif den.is_Add:
                        s = sign(im(b))
                        if s.is_Number and s and den == \
                                log(-factor_terms(b, sign=False)) + s*S.ImaginaryUnit*S.Pi:
                            return S.Exp1**(c*numer(ex))

                obj = b._eval_power(e)
                if obj is not None:
                    return obj
        obj = Expr.__new__(cls, b, e)
        obj.is_commutative = (b.is_commutative and e.is_commutative)
        return obj

def test_issue_12557():
    '''
    # 3200 seconds to compute the fourier part of issue
    import sympy as sym
    x,y,z,t = sym.symbols('x y z t')
    k = sym.symbols("k", integer=True)
    fourier = sym.fourier_series(sym.cos(k*x)*sym.sqrt(x**2),
                                 (x, -sym.pi, sym.pi))
    assert fourier == FourierSeries(
    sqrt(x**2)*cos(k*x), (x, -pi, pi), (Piecewise((pi**2,
    Eq(k, 0)), (2*(-1)**k/k**2 - 2/k**2, True))/(2*pi),
    SeqFormula(Piecewise((pi**2, (Eq(_n, 0) & Eq(k, 0)) | (Eq(_n, 0) &
    Eq(_n, k) & Eq(k, 0)) | (Eq(_n, 0) & Eq(k, 0) & Eq(_n, -k)) | (Eq(_n,
    0) & Eq(_n, k) & Eq(k, 0) & Eq(_n, -k))), (pi**2/2, Eq(_n, k) | Eq(_n,
    -k) | (Eq(_n, 0) & Eq(_n, k)) | (Eq(_n, k) & Eq(k, 0)) | (Eq(_n, 0) &
    Eq(_n, -k)) | (Eq(_n, k) & Eq(_n, -k)) | (Eq(k, 0) & Eq(_n, -k)) |
    (Eq(_n, 0) & Eq(_n, k) & Eq(_n, -k)) | (Eq(_n, k) & Eq(k, 0) & Eq(_n,
    -k))), ((-1)**k*pi**2*_n**3*sin(pi*_n)/(pi*_n**4 - 2*pi*_n**2*k**2 +
    pi*k**4) - (-1)**k*pi**2*_n**3*sin(pi*_n)/(-pi*_n**4 + 2*pi*_n**2*k**2
    - pi*k**4) + (-1)**k*pi*_n**2*cos(pi*_n)/(pi*_n**4 - 2*pi*_n**2*k**2 +
    pi*k**4) - (-1)**k*pi*_n**2*cos(pi*_n)/(-pi*_n**4 + 2*pi*_n**2*k**2 -
    pi*k**4) - (-1)**k*pi**2*_n*k**2*sin(pi*_n)/(pi*_n**4 -
    2*pi*_n**2*k**2 + pi*k**4) +
    (-1)**k*pi**2*_n*k**2*sin(pi*_n)/(-pi*_n**4 + 2*pi*_n**2*k**2 -
    pi*k**4) + (-1)**k*pi*k**2*cos(pi*_n)/(pi*_n**4 - 2*pi*_n**2*k**2 +
    pi*k**4) - (-1)**k*pi*k**2*cos(pi*_n)/(-pi*_n**4 + 2*pi*_n**2*k**2 -
    pi*k**4) - (2*_n**2 + 2*k**2)/(_n**4 - 2*_n**2*k**2 + k**4),
    True))*cos(_n*x)/pi, (_n, 1, oo)), SeqFormula(0, (_k, 1, oo))))
    '''
    x = symbols("x", real=True)
    k = symbols('k',  integer=True)
    abs2 = lambda x: Piecewise((-x, x <= 0), (x, x > 0))
    assert integrate(abs2(x), (x, -pi, pi)) == pi**2
    # attempting a Piecewise rewrite is not automatic so this needs
    # help
    func = cos(k*x)*sqrt(x**2)
    assert integrate(func, (x, -pi, pi)
        ) == Integral(cos(k*x)*Abs(x), (x, -pi, pi))
    assert factor_terms(integrate(func.rewrite(Piecewise), (x, -pi, pi)
        )) == Piecewise(
        (pi**2, Eq(k, 0)), (((-1)**k - 1)/k**2*2, True))

 def simp_hyp(expr):
     return factor_terms(expand_mul(expr)).rewrite(sin)

    def __call__(self, equations, variables=None, method_options=None):
        method_options = extract_method_options(method_options,
                                                {'simplify': True})

        if equations.is_stochastic:
            raise UnsupportedEquationsException('Cannot solve stochastic '
                                                'equations with this state '
                                                'updater.')
        if variables is None:
            variables = {}

        # Get a representation of the ODE system in the form of
        # dX/dt = M*X + B
        varnames, matrix, constants = get_linear_system(equations, variables)

        # No differential equations, nothing to do (this occurs sometimes in the
        # test suite where the whole model is nothing more than something like
        # 'v : 1')
        if matrix.shape == (0, 0):
            return ''

        # Make sure that the matrix M is constant, i.e. it only contains
        # external variables or constant variables
        t = Symbol('t', real=True, positive=True)

        # Check for time dependence
        dt_value = variables['dt'].get_value()[0] if 'dt' in variables else None

        # This will raise an error if we meet the symbol "t" anywhere
        # except as an argument of a locally constant function
        for entry in itertools.chain(matrix, constants):
            if not is_constant_over_dt(entry, variables, dt_value):
                raise UnsupportedEquationsException(
                    ('Expression "{}" is not guaranteed to be constant over a '
                     'time step').format(sympy_to_str(entry)))

        symbols = [Symbol(variable, real=True) for variable in varnames]
        solution = sp.solve_linear_system(matrix.row_join(constants), *symbols)
        if solution is None or set(symbols) != set(solution.keys()):
            raise UnsupportedEquationsException('Cannot solve the given '
                                                'equations with this '
                                                'stateupdater.')
        b = sp.ImmutableMatrix([solution[symbol] for symbol in symbols])

        # Solve the system
        dt = Symbol('dt', real=True, positive=True)
        try:
            A = (matrix * dt).exp()
        except NotImplementedError:
            raise UnsupportedEquationsException('Cannot solve the given '
                                                'equations with this '
                                                'stateupdater.')
        if method_options['simplify']:
            A = A.applyfunc(lambda x:
                            sp.factor_terms(sp.cancel(sp.signsimp(x))))
        C = sp.ImmutableMatrix(A * b) - b
        _S = sp.MatrixSymbol('_S', len(varnames), 1)
        updates = A * _S + C
        updates = updates.as_explicit()

        # The solution contains _S[0, 0], _S[1, 0] etc. for the state variables,
        # replace them with the state variable names 
        abstract_code = []
        for idx, (variable, update) in enumerate(zip(varnames, updates)):
            rhs = update
            if rhs.has(I, re, im):
                raise UnsupportedEquationsException('The solution to the linear system '
                                                    'contains complex values '
                                                    'which is currently not implemented.')
            for row_idx, varname in enumerate(varnames):
                rhs = rhs.subs(_S[row_idx, 0], varname)

            # Do not overwrite the real state variables yet, the update step
            # of other state variables might still need the original values
            abstract_code.append('_' + variable + ' = ' + sympy_to_str(rhs))

        # Update the state variables
        for variable in varnames:
            abstract_code.append('{variable} = _{variable}'.format(variable=variable))
        return '\n'.join(abstract_code)

 def Symplify(A):
     return(factor_terms(simplify(A)))

#u5--rear wheel ang. vel., u6--front wheel ang. vel.

(fr, frstar)= kane.kanes_equations(forceList, bodyList)
kdd = kane.kindiffdict()

#path
path1 = '/home/stefenyin'
path2 = '/home/stefenstudy'

if os.path.exists(path1):
    path = path1
else:
    path = path2

con_forces_noncontri = kane.auxiliary_eqs.applyfunc(
        lambda w: factor_terms(signsimp(w))).subs(kdd)

CF = con_forces_noncontri
"""
#==============================================================================
print('building a file and writing the equations into it')

#-----------------
print('building')
path = path + '/bicycle/bi_equations_writing/con_force_nonslip_non_para_input.txt'

try:
    f = open(path,'w')
    f.write('')
    f.close()
    del f