[Maxima-commits] CVS: maxima/interfaces/xmaxima/Tkmaxima rk4.tcl, 1.1, 1.2

[Jaime E. V. <vi...@us...>]

Update of /cvsroot/maxima/maxima/interfaces/xmaxima/Tkmaxima
In directory 23jxhf1.ch3.sourceforge.com:/tmp/cvs-serv8166

Modified Files:
	rk4.tcl 
Log Message:
Adds a new function, "orthogonal", to calulate numerically the
orthogonal trajectories of the field.


Index: rk4.tcl
===================================================================
RCS file: /cvsroot/maxima/maxima/interfaces/xmaxima/Tkmaxima/rk4.tcl,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -d -r1.1 -r1.2
--- rk4.tcl	31 Oct 2008 14:50:48 -0000	1.1
+++ rk4.tcl	7 Nov 2008 14:50:47 -0000	1.2
@@ -60,3 +60,44 @@
     }
     return $ans
 }
+
+proc orthogonal { f g t0 x0 y0 sx sy nsteps dir} {
+    set n $nsteps
+    set ans "$x0 $y0"
+    set xn $x0
+    set yn $y0
+    set tn $t0
+    catch {
+	while { [incr nsteps -1] >= 0 } {
+	    set kn1 [expr -1*[$g $tn $xn $yn] ]
+	    set ln1 [$f $tn $xn $yn]
+
+	    set h [expr {$dir/[vectorlength [expr {$kn1/$sx}] [expr {$ln1/$sy}]]}]
+	    set h2 [expr {$h / 2.0 }]
+	    set h6 [expr {$h / 6.0 }]
+
+	    set arg [list [expr {$tn + $h2}] [expr {$xn + $h2 * $kn1}] [expr {$yn + $h2*$ln1}]]
+	    set kn2 [expr -1*[eval $g $arg] ]
+	    set ln2 [eval $f $arg]
+
+	    set arg [list [expr {$tn + $h2}] [expr {$xn + $h2 * $kn2}] [expr {$yn +$h2*$ln2}]]
+	    set kn3 [expr -1*[eval $g $arg] ]
+	    set ln3 [eval $f $arg]
+
+	    set arg [list [expr {$tn + $h}] [expr {$xn + $h * $kn3}] [expr {$yn + $h*$ln3}]]
+	    set kn4 [expr -1*[eval $g $arg] ]
+	    set ln4 [eval $f $arg]
+
+	    set dx [expr {$h6 * ($kn1+2*$kn2+2*$kn3+$kn4)}]
+	    set dy [expr {$h6 * ($ln1+2*$ln2+2*$ln3+$ln4)}]
+
+	    if { [vectorlength $dx $dy] > 5*[vectorlength $sx $sy] } { return $ans }
+	    set xn [expr {$xn + $dx}]
+	    set yn [expr {$yn + $dy}]
+	    set tn [expr {$tn + $h}]
+
+	    lappend ans  $xn $yn
+	}
+    }
+    return $ans
+}