(define ((L-central-polar m V) local)
  (let ((q (coordinate local))
        (qdot (velocity local)))
    (let ((r (ref q 0))
          (phi (ref q 1))
          (rdot (ref qdot 0))
          (phidot (ref qdot 1)))
      (- (* 1/2 m
           (+ (square rdot)
              (square (* r phidot))) )
         (V r)))))

;(show-time
; (lambda ()
;   (series:print
;    (((Lie-transform
;       (Lagrangian->Hamiltonian
;	(L-central-polar 'm (lambda (r) (- (/ 'GM r)))))
;;       'dt)
;      state->q)
;     (->H-state 0
;		 (coordinate-tuple 'r_0 'phi_0)
;		 (momentum-tuple 'p_r_0 'p_phi_0)))
;    4)))
;
;
;(up r_0 phi_0)
;(up (/ (* dt p_r_0) m) (/ (* dt p_phi_0) (* m (expt r_0 2))))
;(up
; (+ (/ (* -1/2 GM (expt dt 2)) (* m (expt r_0 2)))
;    (/ (* 1/2 (expt dt 2) (expt p_phi_0 2)) (* (expt m 2) (expt r_0 3))))
; (/ (* -1 (expt dt 2) p_phi_0 p_r_0) (* (expt m 2) (expt r_0 3))))
;(up
; (+
;  (/ (* 1/3 GM (expt dt 3) p_r_0) (* (expt m 2) (expt r_0 3)))
;  (/ (* -1/2 (expt dt 3) (expt p_phi_0 2) p_r_0) (* (expt m 3) (expt r_0 4))))
; (+ (/ (* (expt dt 3) p_phi_0 (expt p_r_0 2)) (* (expt m 3) (expt r_0 4)))
;    (/ (* 1/3 GM (expt dt 3) p_phi_0) (* (expt m 2) (expt r_0 5)))
;    (/ (* -1/3 (expt dt 3) (expt p_phi_0 3)) (* (expt m 3) (expt r_0 6)))))
;process time: 14740 (13780 RUN + 960 GC); real time: 14827
;;Value: ...
;;
;;;;;; Debugging stuff:
;(define (count-entries table)
;;  (let ((max-depth 0) (num-enders 0) (num-extenders 0)
;	(max-num-branches 0) (max-num-enders 0) (max-num-extenders 0))
;    (let lp ((depth 0) (table table))
;      (set! max-depth (max depth max-depth))
;      (set! num-enders
;	    (fix:+ (weak-length (table-enders table)) num-enders))
;      (set! max-num-enders
;	    (max (weak-length (table-enders table)) max-num-enders))
;      (set! num-extenders
;	    (fix:+ (length (table-extenders table)) num-extenders))
;      (set! max-num-extenders
;	    (max (length (table-extenders table)) max-num-extenders))
;      (set! max-num-branches
;	    (max (length (table-subindex table)) max-num-branches))
;      (for-each (lambda (s)
;		  (lp (fix:+ depth 1) (cdr s)))
;		(table-subindex table)))
;    (write-line
;     `(,max-depth ,num-enders ,num-extenders
;		  ,max-num-enders ,max-num-extenders ,max-num-branches))))
;
;(define (find-big-subindex table num-branches)
;  (if (fix:= num-branches (length (table-subindex table)))
;      (begin (write-line (map car (table-subindex table)))
;	     table)
;      (for-each (lambda (s)
;		  (find-big-subindex (cdr s) num-branches))
;		(table-subindex table))))
;
;;;; Stuff to fix
;
;(load "/usr/local/scmutils/relativity/calculus")
;(load "/usr/local/scmutils/relativity/metric")
;
;|#
 
(define cartesian-plane-basis
  (coordinate-system->basis cartesian-plane))
(define d/dx (coordinate-basis-vector-field cartesian-plane 'd/dx 0))
(define d/dy (coordinate-basis-vector-field cartesian-plane 'd/dy 1))
(define dx (coordinate-basis-1form cartesian-plane 'dx 0))
(define dy (coordinate-basis-1form cartesian-plane 'dy 0))


(define polar-basis (coordinate-system->basis polar))
(define r (compose (component 0) (polar '->coords)))
(define theta (compose (component 1) (polar '->coords)))
(define d/dr (coordinate-basis-vector-field polar 'd/dr 0))
(define d/dtheta (coordinate-basis-vector-field polar 'd/dtheta 1))
(define dr (coordinate-basis-1form polar 'dr 0))
(define dtheta (coordinate-basis-1form polar 'dtheta 1))

(define X
  (components->vector-field 
   (up (literal-function 'X^0 (-> (UP Real Real) Real))
       (literal-function 'X^1 (-> (UP Real Real) Real)))
   cartesian-plane
   'X))

(define V
  (components->vector-field 
   (up (literal-function 'V^0 (-> (UP Real Real) Real))
       (literal-function 'V^1 (-> (UP Real Real) Real)))
   cartesian-plane
   'V))

(define (polar-metric v1 v2)
  (let ((1form-basis (basis->1form-basis polar-basis))
	(r (compose (component 0) (polar '->coords))))
    (+ (* ((ref 1form-basis 0) v1)
	  ((ref 1form-basis 0) v2))
       (* (square r)
	  ((ref 1form-basis 1) v1)
	  ((ref 1form-basis 1) v2)))))

(set! *divide-out-terms* #f)

(pe ((metric->first-Christoffel polar-basis polar-metric)
     ((polar '->point) (up 'r^0 'theta^0))))
(down (down (down 0 0) (down 0 r^0))
      (down (down 0 r^0) (down (* -1 r^0) 0)))

;;; This runs out of memory, without intermediate simp.
(pe ((metric->second-Christoffel polar-basis polar-metric)
     ((polar '->point) (up 'r^0 'theta^0))))
(down (down (up 0 0) (up 0 (/ 1 r^0)))
      (down (up 0 (/ 1 r^0)) (up (* -1 r^0) 0)))

(define (clean-weak-list weak-list)
  (if (weak-pair? weak-list)
      (if (weak-pair/car? weak-list)
	  (weak-cons (weak-car weak-list)
		     (clean-weak-list (weak-cdr weak-list)))
	  (clean-weak-list (weak-cdr weak-list)))
      weak-list))

(define (clean-weak-list weak-list)
  (if (weak-pair? weak-list)
      (let scanlp ((scan weak-list))
	(let ((rest (weak-cdr scan)))
	  (if (weak-pair? rest)
	      (if (weak-pair/car? rest)
		  (scanlp rest)
		  (begin
		    (weak-set-cdr! scan (weak-cdr rest))
		    (scanlp scan)))
	      (let frontlp ((front weak-list))
		(if (weak-pair/car? front)
		    front
		    (let ((rest (weak-cdr front)))
		      (if (weak-pair? rest)
			  (frontlp rest)
			  rest)))))))
      weak-list))

(define (clean-weak-alist weak-alist)
  (if (pair? weak-alist)
      (cond ((not (car weak-alist))
	     (clean-weak-alist (cdr weak-alist)))
	    ((weak-pair/car? (car weak-alist))
	     (cons (car weak-alist)
		   (clean-weak-alist (cdr weak-alist))))
	    (else (clean-weak-alist (cdr weak-alist))))
      weak-alist))


(define (clean-subtable-alist alist)
  (if (pair? alist)
      (if (clean-expression-table (cdar alist))
	  (cons (car alist)
		(clean-subtable-alist (cdr alist)))
	  (clean-subtable-alist (cdr alist)))
      alist))