# Full Weber Hamiltonian, N=4, 2D, m_i=1, q=(+1,+1,-1,-1), c=1, kappa=1
# Variables: q = [x1,y1,x2,y2,x3,y3,x4,y4], p = [px1,py1,...,px4,py4]

H = (-1 + (1//2)*(((px2*x2 - px2*x4 - px4*x2 + px4*x4 + py2*y2 - py2*y4 - py4*y2 + py4*y4) / sqrt(x2^2 - 2x2*x4 + x4^2 + y2^2 - 2y2*y4 + y4^2))^2)) / sqrt(x2^2 - 2x2*x4 + x4^2 + y2^2 - 2y2*y4 + y4^2) + (-1 + (1//2)*(((px2*x2 - px2*x3 - px3*x2 + px3*x3 + py2*y2 - py2*y3 - py3*y2 + py3*y3) / sqrt(x2^2 - 2x2*x3 + x3^2 + y2^2 - 2y2*y3 + y3^2))^2)) / sqrt(x2^2 - 2x2*x3 + x3^2 + y2^2 - 2y2*y3 + y3^2) + (-1 + (1//2)*(((px1*x1 - px1*x3 - px3*x1 + px3*x3 + py1*y1 - py1*y3 - py3*y1 + py3*y3) / sqrt(x1^2 - 2x1*x3 + x3^2 + y1^2 - 2y1*y3 + y3^2))^2)) / sqrt(x1^2 - 2x1*x3 + x3^2 + y1^2 - 2y1*y3 + y3^2) + (-1 + (1//2)*(((px1*x1 - px1*x4 - px4*x1 + px4*x4 + py1*y1 - py1*y4 - py4*y1 + py4*y4) / sqrt(x1^2 - 2x1*x4 + x4^2 + y1^2 - 2y1*y4 + y4^2))^2)) / sqrt(x1^2 - 2x1*x4 + x4^2 + y1^2 - 2y1*y4 + y4^2) + (1 - (1//2)*(((px1*x1 - px1*x2 - px2*x1 + px2*x2 + py1*y1 - py1*y2 - py2*y1 + py2*y2) / sqrt(x1^2 - 2x1*x2 + x2^2 + y1^2 - 2y1*y2 + y2^2))^2)) / sqrt(x1^2 - 2x1*x2 + x2^2 + y1^2 - 2y1*y2 + y2^2) + (1 - (1//2)*(((px3*x3 - px3*x4 - px4*x3 + px4*x4 + py3*y3 - py3*y4 - py4*y3 + py4*y4) / sqrt(x3^2 - 2x3*x4 + x4^2 + y3^2 - 2y3*y4 + y4^2))^2)) / sqrt(x3^2 - 2x3*x4 + x4^2 + y3^2 - 2y3*y4 + y4^2) + (1//2)*(px1^2) + (1//2)*(px2^2) + (1//2)*(px3^2) + (1//2)*(px4^2) + (1//2)*(py1^2) + (1//2)*(py2^2) + (1//2)*(py3^2) + (1//2)*(py4^2)
