diff options
| -rw-r--r-- | other/henon-huge.scm | 21 | ||||
| -rw-r--r-- | other/ipython_log.py | 101 | ||||
| -rw-r--r-- | src/from-canon.scm | 178 |
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)) |