blob: a2e06e2f1889bcd1c7ce6e83fef6bfc61ef38992 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
|
(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))
|
