adding some misc files I had sitting around in my working dir, not

really important but might as well hold on to them

3 files changed, 300 insertions, 0 deletions

diff --git a/other/henon-huge.scm b/other/henon-huge.scm
new file mode 100644
index 0000000..8109e24
--- /dev/null
+++ b/other/henon-huge.scm
@@ -0,0 +1,21 @@
+
+;;; want to plot roughly (-1,1), (-1,1), alpha = 1.32
+(define ((henon-map alpha) x y return failure)
+  (if (or (> x 1) (< x -1) (> y 1) (< y -1))
+      failure)
+  (return (- (* x (cos alpha))
+	     (* (- y (square x)) (sin alpha)))
+	  (+ (* x (sin alpha))
+	     (* (- y (square x)) (cos alpha)))))
+
+(
+
+;gnuplot('set terminal png; set xrange[-4:4]; set yrange[-100:100]; set output "/home/bnewbold/fromsage.png"; plot "/home/bnewbold/output";')
+
+(start-gnuplot "huge-output.data")
+
+(define window (frame -1. 1. -1. 1.))
+
+;(explore-map window (henon-map 1.32) 2000)
+
+
diff --git a/other/ipython_log.py b/other/ipython_log.py
new file mode 100644
index 0000000..fa1fd29
--- /dev/null
+++ b/other/ipython_log.py
@@ -0,0 +1,101 @@
+#log# Automatic Logger file. *** THIS MUST BE THE FIRST LINE ***
+#log# DO NOT CHANGE THIS LINE OR THE TWO BELOW
+#log# opts = Struct({'__allownew': True,
+ 'c': '\nimport sage.misc.misc; print sage.misc.misc.branch_current_hg_notice(sage.misc.misc.branch_current_hg()); from sage.misc.interpreter import preparser; preparser(True);import sage.all_cmdline; sage.all_cmdline._init_cmdline(globals());from sage.all import Integer, RealNumber;import os; os.chdir("/home/bnewbold/thesis/other");import sage.misc.interpreter;from sage.misc.interpreter import attached_files;from sage.all_cmdline import *;_=sage.misc.interpreter.load_startup_file("/home/bnewbold/.sage//init.sage");',
+ 'logfile': 'ipython_log.py'})
+#log# args = []
+#log# It is safe to make manual edits below here.
+#log#-----------------------------------------------------------------------
+import sage.misc.misc; print sage.misc.misc.branch_current_hg_notice(sage.misc.misc.branch_current_hg()); from sage.misc.interpreter import preparser; preparser(True);import sage.all_cmdline; sage.all_cmdline._init_cmdline(globals());from sage.all import Integer, RealNumber;import os; os.chdir("/home/bnewbold/thesis/other");import sage.misc.interpreter;from sage.misc.interpreter import attached_files;from sage.all_cmdline import *;_=sage.misc.interpreter.load_startup_file("/home/bnewbold/.sage//init.sage");
+
+
+from pylab import *
+
+d = loadtxt('./huge-output.data')
+plot(d)
+d
+d = d[Integer(0)]
+d
+d = loadtxt('./huge-output.data')
+transpose(d)
+d = transpose(d)
+d = transpose(d)
+#?plot
+
+d
+d = transpose(d)
+d
+p = plot(d[Integer(0)],d[Integer(1)],'b.')
+p
+show(p)
+p
+p = p[Integer(0)]
+p
+show(p)
+_ip.system("rm huge-output.png")
+savefig('huge-output.png')
+d
+max(d[Integer(0)])
+xrange([-Integer(1),Integer(1)])
+#?xrange
+p
+p.xaxis()
+f = figure()
+gca()
+a = gca()
+a.xaxis
+a.xaxis()
+#?a.set_xlim
+#?a.set_xlim
+a.set_xlim((-Integer(1),Integer(1)))
+a.set_ylim((-Integer(1),Integer(1)))
+savefig('huge-output.png')
+p = plot(d[Integer(0)],d[Integer(1)],'b.')
+savefig('huge-output.png')
+a = gca()
+a.set_ylim((-Integer(1),Integer(1)))
+a.set_xlim((-Integer(1),Integer(1)))
+min(d[Integer(0)])
+max(d[Integer(0)])
+min(d[Integer(1)])
+max(d[Integer(1)])
+size(d)
+d
+mean(d)
+mean(d[Integer(0)])
+median(d[Integer(0)])
+median(d[Integer(1)])
+p = plot(d[Integer(0)][:Integer(10)],d[Integer(1)][:Integer(10)],'k.')
+savefig('huge-output.png')
+a = gca()
+a.set_xlim((-Integer(1),Integer(1)))
+a.set_ylim((-Integer(1),Integer(1)))
+savefig('huge-output.png')
+#?plot
+figure()
+plot(d[Integer(0)],d[Integer(1)],'k',markersize=RealNumber('.0001'))
+def doaxes(): a = gca(); a.set_xlim((-Integer(1),Integer(1))); a.set_ylim((-Integer(1),Integer(1)));
+
+doaxes()
+savefig('huge-output.png')
+#?savefig
+savefig('huge-output.png', dpi=Integer(1000))
+savefig('huge-output.png', dpi=Integer(200))
+savefig('huge-output.png', dpi=Integer(200), format='png')
+#?savefig
+savefig('huge-output.png')
+a = gca()
+a.set_xlim[-RealNumber('.001'),RealNumber('.001')]
+a.set_xlim((-RealNumber('.001'),RealNumber('.001')))
+savefig('huge-output.png')
+figure()
+plot(d[Integer(0)][Integer(5000)],d[Integer(1)][Integer(5000)],'k',markersize=Integer(1))
+plot(d[Integer(0)][Integer(5000)],d[Integer(1)][Integer(5000)],'k',markersize=Integer(1))
+d
+size(d)
+size(d[Integer(0)])
+size(d[Integer(1)])
+plot(d[Integer(0)][Integer(5000)],d[Integer(1)][Integer(5000)],'k',markersize=RealNumber('.1'))
+_ip.magic("logon ")
+_ip.magic("logstart ")
+
diff --git a/src/from-canon.scm b/src/from-canon.scm
new file mode 100644
index 0000000..a2e06e2
--- /dev/null
+++ b/src/from-canon.scm
@@ -0,0 +1,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))