;;;  review aspects of SO3

;;; "rotations" is the manifold of rotations in SO3

(define Euler-angles
  (coordinate-system-at 'Euler 'Euler-patch SO3))
(define Euler-angles-chi (Euler-angles '->coords))
(define Euler-angles-chi-inverse (Euler-angles '->point))
(define-coordinates (up theta phi psi) Euler-angles)
(define Euler-angles-basis (coordinate-system->basis Euler-angles))

;;; for the moment define
(define alternate-angles
  (coordinate-system-at 'alternate 'alternate-patch SO3))
(define alternate-angles-chi (alternate-angles '->coords))
(define alternate-angles-chi-inverse (alternate-angles '->point))
(define-coordinates (up vartheta varphi varpsi) alternate-angles)
(define alternate-angles-basis (coordinate-system->basis alternate-angles))

(define time-line the-real-line)
(define time-chi (time-line '->coords))
(define time-chi-inverse (time-line '->point))
(define-coordinates t time-line)


#|
;;;----------------------------------------------------------------
;;; Derive basis for rotations about rectangular coordinate axes.
;;; Of course, these have the generators.

(pe ((D (lambda (xyz)
	  (((D rotate-x) 'theta)
	   ((rotate-x (- 'theta)) xyz))))
     (up 'x 'y 'z)))
(down (up 0 0 0)
      (up 0 0 1)
      (up 0 -1 0))

(pe ((D (lambda (xyz)
	  (((D rotate-y) 'theta)
	   ((rotate-y (- 'theta)) xyz))))
     (up 'x 'y 'z)))
(down (up 0 0 -1)
      (up 0 0 0)
      (up 1 0 0))

(pe ((D (lambda (xyz)
	  (((D rotate-z) 'theta)
	   ((rotate-z (- 'theta)) xyz))))
     (up 'x 'y 'z)))
(down (up 0 1 0)
      (up -1 0 0)
      (up 0 0 0))
|#

;;; tuple multipliers for rotation are obtained with the procedures
;;; rotate-x-tuple, rotate-y-tuple, rotate-z-tuple

#|
want to define vector fields in the group parameter manifold
that correspond to particular spatial rotations

(define (equation-x p q)
  (let ((theta (ref q 0))
	(phi (ref q 1))
	(psi (ref q 2))
	(a (ref p 0))
	(b (ref p 1))
	(c (ref p 2)))
    ((D (lambda (eps)
	  (- (* (rotate-z-tuple (+ phi (* a eps)))
		(rotate-x-tuple (+ theta (* b eps)))
		(rotate-z-tuple (+ psi (* c eps))))
	     (* (rotate-x-tuple eps)
		(rotate-z-tuple phi)
		(rotate-x-tuple theta)
		(rotate-z-tuple psi)))))
     0)))
;;; kind of big

(define (equation2-x p q)
  (let ((theta (ref q 0))
	(phi (ref q 1))
	(psi (ref q 2))
	(a (ref p 0))
	(b (ref p 1))
	(c (ref p 2)))
    ((D (lambda (eps)
	  (- (* (rotate-z-tuple (+ phi (* a eps)))
		(rotate-x-tuple (+ theta (* b eps)))
		(rotate-z-tuple (+ psi (* c eps)))
		(rotate-z-tuple (- psi))
		(rotate-x-tuple (- theta))
		(rotate-z-tuple (- phi)))
	     (rotate-x-tuple eps))))
     0)))

(pe (equation2-x (up 'a 'b 'c) (up 'theta 'phi 'psi)))
(down
 (up 0
     (+ (* c (cos theta)) a)
     (+ (* c (cos phi) (sin theta)) (* -1 b (sin phi))))
 (up (+ (* -1 c (cos theta)) (* -1 a))
     0
     (+ -1 (* c (sin phi) (sin theta)) (* b (cos phi))))
 (up (+ (* -1 c (cos phi) (sin theta)) (* b (sin phi)))
     (+ 1 (* -1 c (sin phi) (sin theta)) (* -1 b (cos phi)))
     0))

;;; mechanical solution

(load "/usr/local/scmutils/src/solve/linreduce")

(define foo-x
  (solve
   (lambda (p)
     (list->vector
      (map simplify
	   (ultra-flatten (equation2-x p (up 'theta 'phi 'psi))))))
   3 9 list))

(pe ((cadr foo-x) #()))
(up (/ (* -1 (sin phi) (cos theta)) (sin theta))
    (cos phi)
    (/ (sin phi) (sin theta)))

;;; e_x is therefore
;;; therefore

(define e_x
  (+ (* (/ (* -1 (sin phi) (cos theta)) (sin theta)) d/dphi)
     (* (cos phi) d/dtheta)
     (* (/ (sin phi) (sin theta)) d/dpsi)))

;;; checks with MTW p.243 (watch notation: phi<->psi )


;;;----------------------------------------------------------------
;;; now for the other two

(define (equation2-z p q)
  (let ((theta (ref q 0))
	(phi (ref q 1))
	(psi (ref q 2))
	(a (ref p 0))
	(b (ref p 1))
	(c (ref p 2)))
    ((D (lambda (eps)
	  (- (* (rotate-z-tuple (+ phi (* a eps)))
		(rotate-x-tuple (+ theta (* b eps)))
		(rotate-z-tuple (+ psi (* c eps)))
		(rotate-z-tuple (- psi))
		(rotate-x-tuple (- theta))
		(rotate-z-tuple (- phi)))
	     (rotate-z-tuple eps))))
     0)))

(pe (equation2-z (up 'a 'b 'c) (up 'theta 'phi 'psi)))
(down
 (up 0
     (+ -1 (* c (cos theta)) a)
     (+ (* c (cos phi) (sin theta)) (* -1 b (sin phi))))
 (up (+ 1 (* -1 c (cos theta)) (* -1 a))
     0
     (+ (* c (sin phi) (sin theta)) (* b (cos phi))))
 (up (+ (* -1 c (cos phi) (sin theta)) (* b (sin phi)))
     (+ (* -1 c (sin phi) (sin theta)) (* -1 b (cos phi)))
     0))

(pe ((cadr foo-z) #()))
(up 1 0 0)

so

(define e_z d/dphi)

;;;----------------------------------------------------------------

(define (equation2-y p q)
  (let ((theta (ref q 0))
	(phi (ref q 1))
	(psi (ref q 2))
	(a (ref p 0))
	(b (ref p 1))
	(c (ref p 2)))
    ((D (lambda (eps)
	  (- (* (rotate-z-tuple (+ phi (* a eps)))
		(rotate-x-tuple (+ theta (* b eps)))
		(rotate-z-tuple (+ psi (* c eps)))
		(rotate-z-tuple (- psi))
		(rotate-x-tuple (- theta))
		(rotate-z-tuple (- phi)))
	     (rotate-y-tuple eps))))
     0)))

(pe (equation2-y (up 'a 'b 'c) (up 'theta 'phi 'psi)))
(down
 (up 0
     (+ (* c (cos theta)) a)
     (+ 1 (* c (cos phi) (sin theta)) (* -1 b (sin phi))))
 (up (+ (* -1 c (cos theta)) (* -1 a))
     0
     (+ (* c (sin phi) (sin theta)) (* b (cos phi))))
 (up (+ -1 (* -1 c (cos phi) (sin theta)) (* b (sin phi)))
     (+ (* -1 c (sin phi) (sin theta)) (* -1 b (cos phi)))
     0))

(pe ((cadr foo-y) #()))
(up (/ (* (cos theta) (cos phi)) (sin theta))
    (sin phi)
    (/ (* -1 (cos phi)) (sin theta)))

(define e_y
  (+ (* (/ (* (cos phi) (cos theta)) (sin theta)) d/dphi)
     (* (sin phi) d/dtheta)
     (* (/ (* -1 (cos phi)) (sin theta)) d/dpsi)))

;;; interchanging psi phi to compare to mtw, checks

;;; summary
|#

(define e_x
  (+ (* (/ (* -1 (sin phi) (cos theta)) (sin theta)) d/dphi)
     (* (cos phi) d/dtheta)
     (* (/ (sin phi) (sin theta)) d/dpsi)))
(define e_y
  (+ (* (/ (* (cos phi) (cos theta)) (sin theta)) d/dphi)
     (* (sin phi) d/dtheta)
     (* (/ (* -1 (cos phi)) (sin theta)) d/dpsi)))
(define e_z d/dphi)

#|

;;;----------------------------------------------------------------

;;; [e_x , e_y] + e_z = 0

(pe (((+ (commutator e_x e_y) e_z)
      (literal-manifold-function 'f Euler-angles))
     (Euler-angles-chi-inverse (up 'theta 'phi 'psi))))
0

;;; [e_y , e_z] + e_x = 0

(pe (((+ (commutator e_y e_z) e_x)
      (literal-manifold-function 'f Euler-angles))
     (Euler-angles-chi-inverse (up 'theta 'phi 'psi))))
0

(pe (((+ (commutator e_z e_x) e_y)
      (literal-manifold-function 'f Euler-angles))
     (Euler-angles-chi-inverse (up 'theta 'phi 'psi))))
0

|#

(define so3-vector-basis
  (down e_x e_y e_z))

(define so3-dual-basis
  (vector-basis->dual so3-vector-basis Euler-angles))

(define so3-basis
  (make-basis so3-vector-basis so3-dual-basis))

#|

(pe ((so3-dual-basis so3-vector-basis)
     (Euler-angles-chi-inverse (up 'theta 'phi 'psi))))
(up (down 1 0 0) (down 0 1 0) (down 0 0 1))

;;; the structure constants
(pe (s:map
     (lambda (e~k)
       (s:map
	(lambda (e_i)
	  (s:map
	   (lambda (e_j)
	     ((e~k (commutator e_i e_j))
	      (Euler-angles-chi-inverse (up 'theta 'phi 'psi))))
	   so3-vector-basis))
	so3-vector-basis))
     so3-dual-basis))
(up (down (down 0 0 0) (down 0 0 -1) (down 0 1 0))
    (down (down 0 0 1) (down 0 0 0) (down -1 0 0))
    (down (down 0 -1 0) (down 1 0 0) (down 0 0 0)))

|#

(define so3-structure-constants
  (up (down (down 0 0 0) (down 0 0 -1) (down 0 1 0))
      (down (down 0 0 1) (down 0 0 0) (down -1 0 0))
      (down (down 0 -1 0) (down 1 0 0) (down 0 0 0))))

;;;----------------------------------------------------------------
;;; Euler velocities as components of differential...

;;; differential of d/dt takes derivatives of functions on rotations
;;; along paths in rotations
#|
(pe
 (let* ((gamma (compose
		Euler-angles-chi-inverse
		(up (literal-function 'f^theta)
		    (literal-function 'f^phi)
		    (literal-function 'f^psi))
		time-chi)))
   ((((differential gamma) d/dt)
     (literal-manifold-function 'f Euler-angles))
    (time-chi-inverse 't))))

(+ (* ((D f^phi) t)
      (((partial 1) f) (up (f^theta t) (f^phi t) (f^psi t))))
   (* ((D f^theta) t)
      (((partial 0) f) (up (f^theta t) (f^phi t) (f^psi t))))
   (* ((D f^psi) t)
      (((partial 2) f) (up (f^theta t) (f^phi t) (f^psi t)))))


;;; components of the differential of d/dt on the Euler basis
(pe (let* ((gamma (compose
		   Euler-angles-chi-inverse
		   (up (literal-function 'f^theta)
		       (literal-function 'f^phi)
		       (literal-function 'f^psi))
		   time-chi))
	   (Euler-basis-over-gamma
	    (basis->basis-over-map gamma Euler-angles-basis))
	   (1form-basis-over-gamma
	    (basis->1form-basis Euler-basis-over-gamma)))
      (s:map (lambda (w) ((w ((differential gamma) d/dt))
			  (time-chi-inverse 't)))
	     1form-basis-over-gamma)))
(up ((D f^theta) t) ((D f^phi) t) ((D f^psi) t))
|#

#|

computing components of quasi-velocity on so3-basis

(pe
 (let* ((gamma (compose
		Euler-angles-chi-inverse
		(up (literal-function 'f^theta)
		    (literal-function 'f^phi)
		    (literal-function 'f^psi))
		time-chi))
	(basis-over-gamma (basis->basis-over-map gamma so3-basis))
	(1form-basis (basis->1form-basis basis-over-gamma))
	(vector-basis (basis->vector-basis basis-over-gamma)))
   (s:map (lambda (w)
	    ((w ((differential gamma) d/dt))
	     (time-chi-inverse 't)))
	  1form-basis)))
(up
 (+ (* (sin (f^phi t)) (sin (f^theta t)) ((D f^psi) t))
    (* (cos (f^phi t)) ((D f^theta) t)))
 (+ (* -1 (cos (f^phi t)) (sin (f^theta t)) ((D f^psi) t))
    (* (sin (f^phi t)) ((D f^theta) t)))
 (+ (* (cos (f^theta t)) ((D f^psi) t)) ((D f^phi) t)))

compare to angular velocity components

(pe ((Euler->omega (up (literal-function 'f^theta)
		       (literal-function 'f^phi)
		       (literal-function 'f^psi)))
     't))
(matrix-by-rows
 (list
  (+ (* (sin (f^phi t)) (sin (f^theta t)) ((D f^psi) t))
     (* (cos (f^phi t)) ((D f^theta) t))))
 (list
  (+ (* -1 (cos (f^phi t)) (sin (f^theta t)) ((D f^psi) t))
     (* (sin (f^phi t)) ((D f^theta) t))))
 (list (+ (* (cos (f^theta t)) ((D f^psi) t)) ((D f^phi) t))))

it checks.

|#

#|

what are the vector fields for which the quasivelocities are the
body components of the angular velocities?

omega_body = M^T omega = M^-1 omega
M = Rz(phi) Rx(theta) Rz(psi)
M^T = M^-1 = Rz(-psi) Rx(-theta) Rz(-phi)

spatial components of vectors rotate by M^-1 to get body components

e_x, e_y, e_z are the basis vector fields that correspond
to rotation about the x, y, z axes
they are vector fields on the space of euler angles

e_xp, e_yp, e_zp are basis vectors that correspond to rotation
about body xp, yp, zp axes

ep(f) = e(f) M

(pe (* (rotate-z-tuple 'phi)
       (* (rotate-x-tuple 'theta)
	  (rotate-z-tuple 'psi))))
(down
 (up (+ (* -1 (sin psi) (cos theta) (sin phi)) (* (cos psi) (cos phi)))
     (+ (* (sin psi) (cos theta) (cos phi)) (* (sin phi) (cos psi)))
     (* (sin theta) (sin psi)))
 (up (+ (* -1 (cos theta) (sin phi) (cos psi)) (* -1 (sin psi) (cos phi)))
     (+ (* (cos theta) (cos psi) (cos phi)) (* -1 (sin psi) (sin phi)))
     (* (sin theta) (cos psi)))
 (up (* (sin theta) (sin phi)) (* -1 (sin theta) (cos phi)) (cos theta)))

(pe (* (down 'e_x 'e_y 'e_z)
	(* (rotate-z-tuple 'phi)
	   (* (rotate-x-tuple 'theta)
	      (rotate-z-tuple 'psi)))))
(down
 (+ (* -1 (sin psi) (cos theta) (sin phi) e_x)
    (* (sin psi) (cos theta) (cos phi) e_y)
    (* (cos psi) (cos phi) e_x)
    (* (sin phi) (cos psi) e_y)
    (* (sin theta) (sin psi) e_z))
 (+ (* -1 (cos theta) (sin phi) (cos psi) e_x)
    (* (cos theta) (cos psi) (cos phi) e_y)
    (* -1 (sin psi) (cos phi) e_x)
    (* -1 (sin psi) (sin phi) e_y)
    (* (sin theta) (cos psi) e_z))
 (+ (* (sin theta) (sin phi) e_x)
    (* -1 (sin theta) (cos phi) e_y)
    (* (cos theta) e_z)))

(pe ((dtheta
      (down
       (+ (* -1 (sin psi) (cos theta) (sin phi) e_x)
	  (* (sin psi) (cos theta) (cos phi) e_y)
	  (* (cos psi) (cos phi) e_x)
	  (* (sin phi) (cos psi) e_y)
	  (* (sin theta) (sin psi) e_z))
       (+ (* -1 (cos theta) (sin phi) (cos psi) e_x)
	  (* (cos theta) (cos psi) (cos phi) e_y)
	  (* -1 (sin psi) (cos phi) e_x)
	  (* -1 (sin psi) (sin phi) e_y)
	  (* (sin theta) (cos psi) e_z))
       (+ (* (sin theta) (sin phi) e_x)
	  (* -1 (sin theta) (cos phi) e_y)
	  (* (cos theta) e_z))))
     (Euler-angles-chi-inverse
      (up 'theta 'phi 'psi))))
(down (cos psi) (* -1 (sin psi)) 0)

(pe ((dphi
      (down
       (+ (* -1 (sin psi) (cos theta) (sin phi) e_x)
	  (* (sin psi) (cos theta) (cos phi) e_y)
	  (* (cos psi) (cos phi) e_x)
	  (* (sin phi) (cos psi) e_y)
	  (* (sin theta) (sin psi) e_z))
       (+ (* -1 (cos theta) (sin phi) (cos psi) e_x)
	  (* (cos theta) (cos psi) (cos phi) e_y)
	  (* -1 (sin psi) (cos phi) e_x)
	  (* -1 (sin psi) (sin phi) e_y)
	  (* (sin theta) (cos psi) e_z))
       (+ (* (sin theta) (sin phi) e_x)
	  (* -1 (sin theta) (cos phi) e_y)
	  (* (cos theta) e_z))))
     (Euler-angles-chi-inverse
      (up 'theta 'phi 'psi))))
(down (/ (sin psi) (sin theta)) (/ (cos psi) (sin theta)) 0)

(pe ((dpsi
      (down
       (+ (* -1 (sin psi) (cos theta) (sin phi) e_x)
	  (* (sin psi) (cos theta) (cos phi) e_y)
	  (* (cos psi) (cos phi) e_x)
	  (* (sin phi) (cos psi) e_y)
	  (* (sin theta) (sin psi) e_z))
       (+ (* -1 (cos theta) (sin phi) (cos psi) e_x)
	  (* (cos theta) (cos psi) (cos phi) e_y)
	  (* -1 (sin psi) (cos phi) e_x)
	  (* -1 (sin psi) (sin phi) e_y)
	  (* (sin theta) (cos psi) e_z))
       (+ (* (sin theta) (sin phi) e_x)
	  (* -1 (sin theta) (cos phi) e_y)
	  (* (cos theta) e_z))))
     (Euler-angles-chi-inverse
      (up 'theta 'phi 'psi))))
(down (/ (* -1 (cos theta) (sin psi)) (sin theta))
      (/ (* -1 (cos psi) (cos theta)) (sin theta))
      1)

#|	  
d/dtheta (down (cos psi)
	       (* -1 (sin psi))
	       0)
d/dphi (down (/ (sin psi) (sin theta))
	     (/ (cos psi) (sin theta))
	     0)
d/dpsi (down (/ (* -1 (sin psi) (cos theta)) (sin theta))
      (/ (* -1 (cos psi) (cos theta)) (sin theta))
      1)

deduce that
|#
|#

(define ep_x
  (+ (* (cos psi) d/dtheta)
     (* (/ (sin psi) (sin theta)) d/dphi)
     (* (/ (* -1 (sin psi) (cos theta)) (sin theta)) d/dpsi)))
(define ep_y
  (+ (* (* -1 (sin psi)) d/dtheta)
     (* (/ (cos psi) (sin theta)) d/dphi)
     (* (/ (* -1 (cos psi) (cos theta)) (sin theta)) d/dpsi)))
(define ep_z d/dpsi)

#|

what are the commutators?

;;; [ep_x , ep_y] - ep_z = 0

(pe (((- (commutator ep_x ep_y) ep_z)
      (literal-manifold-function 'f Euler-angles))
     (Euler-angles-chi-inverse
      (up 'theta 'phi 'psi))))
0

;;; [ep_y , ep_z] - ep_x = 0

(pe (((- (commutator ep_y ep_z) ep_x)
      (literal-manifold-function 'f Euler-angles))
     (Euler-angles-chi-inverse
      (up 'theta 'phi 'psi))))
0

;;; [ep_z , ep_x] - ep_y = 0

(pe (((- (commutator ep_z ep_x) ep_y)
      (literal-manifold-function 'f Euler-angles))
     (Euler-angles-chi-inverse
      (up 'theta 'phi 'psi))))
0

the signs of the structure constants are opposite
looks promising.

|#

(define so3p-vector-basis
  (down ep_x ep_y ep_z))

(define so3p-dual-basis
  (vector-basis->dual so3p-vector-basis Euler-angles))

(define so3p-basis
  (make-basis so3p-vector-basis so3p-dual-basis))

#|

;;; the structure constants
(pe (s:map
     (lambda (e~k)
       (s:map
	(lambda (e_i)
	  (s:map
	   (lambda (e_j)
	     ((e~k (commutator e_i e_j))
	      (Euler-angles-chi-inverse
	       (up 'theta 'phi 'psi))))
	   so3p-vector-basis))
	so3p-vector-basis))
     so3p-dual-basis))

(up (down (down 0 0 0) (down 0 0 1) (down 0 -1 0))
    (down (down 0 0 -1) (down 0 0 0) (down 1 0 0))
    (down (down 0 1 0) (down -1 0 0) (down 0 0 0)))

|#

(define so3p-structure-constants
  (up (down (down 0 0 0) (down 0 0 1) (down 0 -1 0))
      (down (down 0 0 -1) (down 0 0 0) (down 1 0 0))
      (down (down 0 1 0) (down -1 0 0) (down 0 0 0))))

#|

;;quasivelocities on the so3p basis

(pe
 (let* ((gamma (compose
		Euler-angles-chi-inverse
		(up (literal-function 'f^theta)
		    (literal-function 'f^phi)
		    (literal-function 'f^psi))
		time-chi))
	(basis-over-gamma (basis->basis-over-map gamma so3p-basis))
	(1form-basis (basis->1form-basis basis-over-gamma))
	(vector-basis (basis->vector-basis basis-over-gamma)))
   (s:map (lambda (w)
	    ((w ((differential gamma) d/dt))
	     (time-chi-inverse 't)))
	  1form-basis)))

(up
 (+ (* ((D f^phi) t) (sin (f^theta t)) (sin (f^psi t)))
    (* (cos (f^psi t)) ((D f^theta) t)))
 (+ (* ((D f^phi) t) (sin (f^theta t)) (cos (f^psi t)))
    (* -1 (sin (f^psi t)) ((D f^theta) t)))
 (+ (* ((D f^phi) t) (cos (f^theta t))) ((D f^psi) t)))

;;compare to omega body

(pe ((Euler->omega-body
      (up (literal-function 'theta)
	  (literal-function 'phi)
	  (literal-function 'psi))) 't))
(matrix-by-rows
 (list
  (+ (* (sin (theta t)) (sin (psi t)) ((D phi) t))
     (* ((D theta) t) (cos (psi t)))))
 (list
  (+ (* (sin (theta t)) (cos (psi t)) ((D phi) t))
     (* -1 ((D theta) t) (sin (psi t)))))
 (list (+ (* (cos (theta t)) ((D phi) t)) ((D psi) t))))

worked.

check the determining equation

(pe (-
     (*
      ((so3-vector-basis
	(literal-manifold-function 'f Euler-angles))
       (Euler-angles-chi-inverse
	(up 'theta 'phi 'psi)))
      (* (rotate-z-tuple 'phi)
	 (* (rotate-x-tuple 'theta)
	    (rotate-z-tuple 'psi))))
     ((so3p-vector-basis
       (literal-manifold-function 'f Euler-angles))
      (Euler-angles-chi-inverse
	(up 'theta 'phi 'psi)))))
(down 0 0 0)

|#

#|

;;;****************************************************************
;;; Instead of Euler use an alternate set of angles ...
;;;****************************************************************

|#

(define (Angles->M q)
  (let ((vartheta (ref q 0))
	(varphi (ref q 1))
	(varpsi (ref q 2)))
    (* (rotate-z-tuple varphi)
       (rotate-x-tuple vartheta)
       (rotate-y-tuple varpsi))))

(define (Angles->M^-1 q)
  (let ((vartheta (ref q 0))
	(varphi (ref q 1))
	(varpsi (ref q 2)))
    (* (rotate-y-tuple (- varpsi))
       (rotate-x-tuple (- vartheta))
       (rotate-z-tuple (- varphi)))))
#|

;;; (load "tracetranspose")
(define ((equation0 vartheta varphi varpsi) p)
  (let ((a (ref p 0)) (b (ref p 1)) (c (ref p 2)))
    (let ((M ((D (lambda (t) (- (* (Angles->M (up (+ vartheta (* a t))
						  (+ varphi (* b t))
						  (+ varpsi (* c t))))
				   (Angles->M^-1 (up vartheta varphi varpsi)))
				(rotate-x-tuple t))))
		 0)))
	 (up (ref M 0 1)
	     (ref M 0 2)
	     (ref M 1 2)))))
(pe ((equation0 'vartheta 'varphi 'varpsi) (up 'a 'b 'c)))
(up (+ (* c (sin vartheta)) b)
    (+ (* -1 c (cos varphi) (cos vartheta)) (* -1 a (sin varphi)))
    (+ -1 (* -1 c (cos vartheta) (sin varphi)) (* a (cos varphi))))

(pe ((cadr (solve (equation0 'vartheta 'varphi 'varpsi) 3 3 list)) #()))
(up (cos varphi) (* (tan vartheta) (sin varphi)) (/ (* -1 (sin varphi)) (cos vartheta)))

(define ((equation1 vartheta varphi varpsi) p)
  (let ((a (ref p 0)) (b (ref p 1)) (c (ref p 2)))
    (let ((M ((D (lambda (t) (- (* (Angles->M (up (+ vartheta (* a t))
						  (+ varphi (* b t))
						  (+ varpsi (* c t))))
				   (Angles->M^-1 (up vartheta varphi varpsi)))
				(rotate-y-tuple t))))
		 0)))
	 (up (ref M 0 1)
	     (ref M 0 2)
	     (ref M 1 2)))))
(pe ((cadr (solve (equation1 'vartheta 'varphi 'varpsi) 3 3 list)) #()))
(up (sin varphi) (* -1 (tan vartheta) (cos varphi)) (/ (cos varphi) (cos vartheta)))

(define ((equation2 vartheta varphi varpsi) p)
  (let ((a (ref p 0)) (b (ref p 1)) (c (ref p 2)))
    (let ((M ((D (lambda (t) (- (* (Angles->M (up (+ vartheta (* a t))
						  (+ varphi (* b t))
						  (+ varpsi (* c t))))
				   (Angles->M^-1 (up vartheta varphi varpsi)))
				(rotate-z-tuple t))))
		 0)))
	 (up (ref M 0 1)
	     (ref M 0 2)
	     (ref M 1 2)))))
(pe ((cadr (solve (equation2 'vartheta 'varphi 'varpsi) 3 3 list)) #()))
(up 0 1 0)

|#

(define ea_x
  (+ (* (cos varphi) d/dvartheta)
     (* (* (sin varphi) (tan vartheta)) d/dvarphi)
     (* (/ (* -1 (sin varphi)) (cos vartheta)) d/dvarpsi)))

(define ea_y
  (+ (* (sin varphi) d/dvartheta)
     (* (* -1 (cos varphi) (tan vartheta)) d/dvarphi)
     (* (/ (cos varphi) (cos vartheta)) d/dvarpsi)))

(define ea_z d/dvarphi)

(define so3a-vector-basis
  (down ea_x ea_y ea_z))

(define so3a-dual-basis
  (vector-basis->dual so3a-vector-basis alternate-angles))

(define so3a-basis
  (make-basis so3a-vector-basis so3a-dual-basis))

#|

(pe (((+ (commutator ea_x ea_y) ea_z)
      (literal-manifold-function 'f alternate-angles))
     (alternate-angles-chi-inverse
      (up 'vartheta 'varphi 'varpsi))))

0

(pe (((+ (commutator ea_y ea_z) ea_x)
      (literal-manifold-function 'f alternate-angles))
     (alternate-angles-chi-inverse
      (up 'vartheta 'varphi 'varpsi))))
0

(pe (((+ (commutator ea_z ea_x) ea_y)
      (literal-manifold-function 'f alternate-angles))
     (alternate-angles-chi-inverse
      (up 'vartheta 'varphi 'varpsi))))
0

(pe ((so3a-dual-basis so3a-vector-basis)
     (alternate-angles-chi-inverse
      (up 'vartheta 'varphi 'varpsi))))
(up (down 1 0 0) (down 0 1 0) (down 0 0 1))

;;; the structure constants
(pe (s:map
     (lambda (e~k)
       (s:map
	(lambda (e_i)
	  (s:map
	   (lambda (e_j)
	     ((e~k (commutator e_i e_j))
	      (alternate-angles-chi-inverse
	       (up 'vartheta 'varphi 'varpsi))))
	   so3a-vector-basis))
	so3a-vector-basis))
     so3a-dual-basis))
(up (down (down 0 0 0) (down 0 0 -1) (down 0 1 0))
    (down (down 0 0 1) (down 0 0 0) (down -1 0 0))
    (down (down 0 -1 0) (down 1 0 0) (down 0 0 0)))

|#

(define so3a-structure-constants
  (up (down (down 0 0 0) (down 0 0 -1) (down 0 1 0))
      (down (down 0 0 1) (down 0 0 0) (down -1 0 0))
      (down (down 0 -1 0) (down 1 0 0) (down 0 0 0))))

#|

(pe (* (down 'ea_x 'ea_y 'ea_z)
       (Angles->M (up 'vartheta 'varphi 'varpsi))))
(down
 (+ (* -1 ea_x (sin varpsi) (sin vartheta) (sin varphi))
    (* ea_y (sin varpsi) (sin vartheta) (cos varphi))
    (* ea_x (cos varpsi) (cos varphi))
    (* ea_y (sin varphi) (cos varpsi))
    (* -1 ea_z (cos vartheta) (sin varpsi)))
 (+ (* -1 ea_x (cos vartheta) (sin varphi))
    (* ea_y (cos vartheta) (cos varphi))
    (* ea_z (sin vartheta)))
 (+ (* ea_x (sin vartheta) (sin varphi) (cos varpsi))
    (* -1 ea_y (sin vartheta) (cos varpsi) (cos varphi))
    (* ea_x (sin varpsi) (cos varphi))
    (* ea_y (sin varpsi) (sin varphi))
    (* ea_z (cos vartheta) (cos varpsi))))

(pe ((dvartheta (down
	      (+ (* -1 ea_x (sin varpsi) (sin vartheta) (sin varphi))
		 (* ea_y (sin varpsi) (sin vartheta) (cos varphi))
		 (* ea_x (cos varpsi) (cos varphi))
		 (* ea_y (sin varphi) (cos varpsi))
		 (* -1 ea_z (cos vartheta) (sin varpsi)))
	      (+ (* -1 ea_x (cos vartheta) (sin varphi))
		 (* ea_y (cos vartheta) (cos varphi))
		 (* ea_z (sin vartheta)))
	      (+ (* ea_x (sin vartheta) (sin varphi) (cos varpsi))
		 (* -1 ea_y (sin vartheta) (cos varpsi) (cos varphi))
		 (* ea_x (sin varpsi) (cos varphi))
		 (* ea_y (sin varpsi) (sin varphi))
		 (* ea_z (cos vartheta) (cos varpsi)))))
     (alternate-angles-chi-inverse
      (up 'vartheta 'varphi 'varpsi))))
(down (cos varpsi) 0 (sin varpsi))

(pe ((dvarphi (down
	      (+ (* -1 ea_x (sin varpsi) (sin vartheta) (sin varphi))
		 (* ea_y (sin varpsi) (sin vartheta) (cos varphi))
		 (* ea_x (cos varpsi) (cos varphi))
		 (* ea_y (sin varphi) (cos varpsi))
		 (* -1 ea_z (cos vartheta) (sin varpsi)))
	      (+ (* -1 ea_x (cos vartheta) (sin varphi))
		 (* ea_y (cos vartheta) (cos varphi))
		 (* ea_z (sin vartheta)))
	      (+ (* ea_x (sin vartheta) (sin varphi) (cos varpsi))
		 (* -1 ea_y (sin vartheta) (cos varpsi) (cos varphi))
		 (* ea_x (sin varpsi) (cos varphi))
		 (* ea_y (sin varpsi) (sin varphi))
		 (* ea_z (cos vartheta) (cos varpsi)))))
     (alternate-angles-chi-inverse
      (up 'vartheta 'varphi 'varpsi))))
(down (/ (* -1 (sin varpsi)) (cos vartheta)) 0 (/ (cos varpsi) (cos vartheta)))

(pe ((dvarpsi (down
	      (+ (* -1 ea_x (sin varpsi) (sin vartheta) (sin varphi))
		 (* ea_y (sin varpsi) (sin vartheta) (cos varphi))
		 (* ea_x (cos varpsi) (cos varphi))
		 (* ea_y (sin varphi) (cos varpsi))
		 (* -1 ea_z (cos vartheta) (sin varpsi)))
	      (+ (* -1 ea_x (cos vartheta) (sin varphi))
		 (* ea_y (cos vartheta) (cos varphi))
		 (* ea_z (sin vartheta)))
	      (+ (* ea_x (sin vartheta) (sin varphi) (cos varpsi))
		 (* -1 ea_y (sin vartheta) (cos varpsi) (cos varphi))
		 (* ea_x (sin varpsi) (cos varphi))
		 (* ea_y (sin varpsi) (sin varphi))
		 (* ea_z (cos vartheta) (cos varpsi)))))
     (alternate-angles-chi-inverse
      (up 'vartheta 'varphi 'varpsi))))
(down (* (sin varpsi) (tan vartheta)) 1 (* -1 (cos varpsi) (tan vartheta)))

d/dvartheta
(down (cos varpsi) 0 (sin varpsi))
d/dvarphi
(down (/ (* -1 (sin varpsi)) (cos vartheta)) 0 (/ (cos varpsi) (cos vartheta)))
d/dvarpsi
(down (* (sin varpsi) (tan vartheta)) 1 (* -1 (cos varpsi) (tan vartheta)))

|#

(define eap_x
  (+ (* (cos varpsi) d/dvartheta)
     (* (/ (* -1 (sin varpsi)) (cos vartheta)) d/dvarphi)
     (* (* (sin varpsi) (tan vartheta)) d/dvarpsi)))

(define eap_y
  d/dvarpsi)

(define eap_z
  (+ (* (sin varpsi) d/dvartheta)
     (* (/ (cos varpsi) (cos vartheta)) d/dvarphi)
     (* (* -1 (cos varpsi) (tan vartheta)) d/dvarpsi)))

#|

(pe (((- (commutator eap_x eap_y) eap_z)
      (literal-manifold-function 'f alternate-angles))
     (alternate-angles-chi-inverse
      (up 'vartheta 'varphi 'varpsi))))
0

(pe (((- (commutator eap_y eap_z) eap_x)
      (literal-manifold-function 'f alternate-angles))
     (alternate-angles-chi-inverse
      (up 'vartheta 'varphi 'varpsi))))
0

(pe (((- (commutator eap_z eap_x) eap_y)
      (literal-manifold-function 'f alternate-angles))
     (alternate-angles-chi-inverse
      (up 'vartheta 'varphi 'varpsi))))
0

|#

(define so3ap-vector-basis
  (down eap_x eap_y eap_z))

(define so3ap-dual-basis
  (vector-basis->dual so3ap-vector-basis alternate-angles))

(define so3ap-basis
  (make-basis so3ap-vector-basis so3ap-dual-basis))

#|

(pe ((s:map
      so3ap-dual-basis
      so3ap-vector-basis)
     (alternate-angles-chi-inverse
      (up 'vartheta 'varphi 'varpsi))))
(down (up 1 0 0) (up 0 1 0) (up 0 0 1))

;;; structure-constants
(pe (s:map
     (lambda (e~k)
       (s:map
	(lambda (e_i)
	  (s:map
	   (lambda (e_j)
	     ((e~k (commutator e_i e_j))
	      (alternate-angles-chi-inverse
	       (up 'vartheta 'varphi 'varpsi))))
	   so3ap-vector-basis))
	so3ap-vector-basis))
     so3ap-dual-basis))
(up (down (down 0 0 0) (down 0 0 1) (down 0 -1 0))
    (down (down 0 0 -1) (down 0 0 0) (down 1 0 0))
    (down (down 0 1 0) (down -1 0 0) (down 0 0 0)))

|#

(define so3ap-structure-constants
  (up (down (down 0 0 0) (down 0 0 1) (down 0 -1 0))
      (down (down 0 0 -1) (down 0 0 0) (down 1 0 0))
      (down (down 0 1 0) (down -1 0 0) (down 0 0 0))))

#|

;;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

check angular velocities in alternate angles

(pe (let ((q (up (literal-function 'vartheta)
		 (literal-function 'varphi)
		 (literal-function 'varpsi))))
      (wcross->w (* (Angles->M^-1 (q 't))
		    ((D (lambda (t) (Angles->M (q t)))) 't)))))
(up
 (+ (* -1 (cos (vartheta t)) (sin (varpsi t)) ((D varphi) t))
    (* ((D vartheta) t) (cos (varpsi t))))
 (+ (* (sin (vartheta t)) ((D varphi) t)) ((D varpsi) t))
 (+ (* (cos (vartheta t)) (cos (varpsi t)) ((D varphi) t))
    (* ((D vartheta) t) (sin (varpsi t)))))

(pe
 (let* ((gamma (compose
		alternate-angles-chi-inverse
		(up (literal-function 'vartheta)
		    (literal-function 'varphi)
		    (literal-function 'varpsi))
		time-chi))
	(basis-over-gamma (basis->basis-over-map gamma so3ap-basis))
	(1form-basis (basis->1form-basis basis-over-gamma))
	(vector-basis (basis->vector-basis basis-over-gamma)))
   (s:map (lambda (w)
	    ((w ((differential gamma) d/dt))
	     (time-chi-inverse 't)))
	  1form-basis)))
(up
 (+ (* -1 (cos (vartheta t)) (sin (varpsi t)) ((D varphi) t))
    (* ((D vartheta) t) (cos (varpsi t))))
 (+ (* (sin (vartheta t)) ((D varphi) t)) ((D varpsi) t))
 (+ (* (cos (vartheta t)) (cos (varpsi t)) ((D varphi) t))
    (* ((D vartheta) t) (sin (varpsi t)))))

|#