import sympy
sympy.init_printing()
from sympy import sin, pi
q, p = sympy.symbols("q p")
K = sympy.symbols("K")
mapping = sympy.Matrix([q + p,
p + K/(2*pi) * sin(q+p)])
diff_vars = sympy.Matrix([q, p])
mapping
# Compute linearized map
DP = mapping.jacobian(diff_vars)
DP
# Check that it is area preserving
sympy.det(DP)