(* Mathematica program for Problem 8.8

  This program is used to calculate the SO (4) group in 3D Coulomb
 potential \ V(r) = -1/r 

  date : 11/11/2009
  author: Jorge C. Romão
  email : jorge.romao@ist.utl.pt
*)

(*  Definitions *)
(*
Comumutator3D is a function to calculate commutators in 3
dimensions, where f and g are variables to be computed 
pi and qi (i=1,2,3)  are canonical coordinates and their
conjugate momenta respectively 
H is Hamiltonian of system, h is Planck constant, l is
angular momentum and A is Rungz-Lenz  vector in units where
m=1, \hbar=1, c=1 and \alpha=1
*)

Commutator3D[f_, g_]:= I h (D[f, q1] D[g, p1]- D[f,p1] D[g, q1]
			    +D[f, q2] D[g, p2] -D[f, p2] D[g,q2]
			    + D[f, q3] D[g, p3]- D[f, p3] D[g, q3] );

H=(p1^2+p2^2+p3^2)/2 - 1/Sqrt[q1^2+q2^2+q3^2];

h1=(p1^2+p2^2+p3^2)/2;
h2=- 1/Sqrt[q1^2+q2^2+q3^2];

(* Angular Momenta *)
L1=q2 p3-q3 p2 ;
L2=q3 p1-q1 p3;
L3=q1 p2-q2 p1;

(* Runge-Lenz Vector *)
A1=p2 L3 - p3 L2 - q1/Sqrt[q1^2+q2^2+q3^2];
A2=p3 L1 - p1 L3 - q2/Sqrt[q1^2+q2^2+q3^2];
A3=p1 L2 - p2 L1 - q3/Sqrt[q1^2+q2^2+q3^2];

(* Check usual Commutation [q_i,p_j] relations *)
Print["[q1,p1] = ",Commutator3D[q1,p1]//FullSimplify]
Print["[q1,p2] = ",Commutator3D[q1,p2]//FullSimplify]
Print["[q1,p3] = ",Commutator3D[q1,p3]//FullSimplify]
Print["[q2,p1] = ",Commutator3D[q2,p1]//FullSimplify]
Print["[q2,p2] = ",Commutator3D[q2,p2]//FullSimplify]
Print["[q2,p3] = ",Commutator3D[q2,p3]//FullSimplify]
Print["[q3,p1] = ",Commutator3D[q3,p1]//FullSimplify]
Print["[q3,p2] = ",Commutator3D[q3,p2]//FullSimplify]
Print["[q3,p3] = ",Commutator3D[q3,p3]//FullSimplify]

(* Check Angular Momentum Commutation relations *)
Print["[L1,L2] = ",Commutator3D[L1,L2]//FullSimplify]
Print["[L2,L3] = ",Commutator3D[L2,L3]//FullSimplify]
Print["[L3,L1] = ",Commutator3D[L3,L1]//FullSimplify]

(* Check  [L_i,q_j]  *)
Print["[L1,q1] = ",Commutator3D[L1,q1]//FullSimplify]
Print["[L1,q2] = ",Commutator3D[L1,q2]//FullSimplify]
Print["[L1,q3] = ",Commutator3D[L1,q3]//FullSimplify]
Print["[L2,q1] = ",Commutator3D[L2,q1]//FullSimplify]
Print["[L2,q2] = ",Commutator3D[L2,q2]//FullSimplify]
Print["[L2,q3] = ",Commutator3D[L2,q3]//FullSimplify]
Print["[L3,q1] = ",Commutator3D[L3,q1]//FullSimplify]
Print["[L3,q2] = ",Commutator3D[L3,q2]//FullSimplify]
Print["[L3,q3] = ",Commutator3D[L3,q3]//FullSimplify]

(* Check  [L_i,p_j]  *)
Print["[L1,p1] = ",Commutator3D[L1,p1]//FullSimplify]
Print["[L1,p2] = ",Commutator3D[L1,p2]//FullSimplify]
Print["[L1,p3] = ",Commutator3D[L1,p3]//FullSimplify]
Print["[L2,p1] = ",Commutator3D[L2,p1]//FullSimplify]
Print["[L2,p2] = ",Commutator3D[L2,p2]//FullSimplify]
Print["[L2,p3] = ",Commutator3D[L2,p3]//FullSimplify]
Print["[L3,p1] = ",Commutator3D[L3,p1]//FullSimplify]
Print["[L3,p2] = ",Commutator3D[L3,p2]//FullSimplify]
Print["[L3,p3] = ",Commutator3D[L3,p3]//FullSimplify]

(* Commutation with the Hamiltonian *)
Print["[H,A1] = ",Commutator3D[H,A1]//FullSimplify]
Print["[H,A2] = ",Commutator3D[H,A1]//FullSimplify]
Print["[H,A3] = ",Commutator3D[H,A1]//FullSimplify]

(* Commutation [L_i, A_j] *)
Print["[L1,A1] = ",Commutator3D[L1,A1]//FullSimplify]
Print["[L1,A2] = ",Commutator3D[L1,A2]//FullSimplify]
Print["[L1,A3] = ",Commutator3D[L1,A3]//FullSimplify]
Print["[L2,A1] = ",Commutator3D[L2,A1]//FullSimplify]
Print["[L2,A2] = ",Commutator3D[L2,A2]//FullSimplify]
Print["[L2,A3] = ",Commutator3D[L2,A3]//FullSimplify]
Print["[L3,A1] = ",Commutator3D[L3,A1]//FullSimplify]
Print["[L3,A2] = ",Commutator3D[L3,A2]//FullSimplify]
Print["[L3,A3] = ",Commutator3D[L3,A3]//FullSimplify]

(* Check [A_i,A_j] Commutation relations *)
Print["[A1,A1] = ",Commutator3D[A1,A1]//FullSimplify]
Print["[A1,A2] = ",Commutator3D[A1,A2]//FullSimplify]
Print["[A1,A3] = ",Commutator3D[A1,A3]//FullSimplify]
Print["[A2,A1] = ",Commutator3D[A2,A1]//FullSimplify]
Print["[A2,A2] = ",Commutator3D[A2,A2]//FullSimplify]
Print["[A2,A3] = ",Commutator3D[A2,A3]//FullSimplify]
Print["[A3,A1] = ",Commutator3D[A3,A1]//FullSimplify]
Print["[A3,A2] = ",Commutator3D[A3,A2]//FullSimplify]
Print["[A3,A3] = ",Commutator3D[A3,A3]//FullSimplify]

(* Check that [L^2,A_i] != 0 *)
Print["[L^2,A1] = ",Commutator3D[L1^2+L2^2+L3^2,A1]//FullSimplify]
Print["[L^2,A2] = ",Commutator3D[L1^2+L2^2+L3^2,A2]//FullSimplify]
Print["[L^2,A3] = ",Commutator3D[L1^2+L2^2+L3^2,A3]//FullSimplify]