minv/recursion.py
changeset 0 0efde00f9229
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/minv/recursion.py	Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,48 @@
+#!/usr/bin/env python
+# 11.1
+import os, sys
+sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"]
+
+import pylab as pl
+from polyfuncs import poladd
+
+def recursion(A, dA, C, dC, j):
+    Fo, dFo = C, dC
+    Eo, dEo = 1, 0
+    A_z = A[1:dA+1]
+    dA_z = dA-1
+    zi, dzi = 1, 0
+    for i in range(j-1):
+        if dFo == 0:
+            Fn1 = 0
+        else:
+            Fn1 = Fo[1:dFo+1]
+        dFn1 = max(dFo-1, 0)
+        Fn2 = -Fo[0]*A_z
+        dFn2 = dA-1
+        Fn, dFn = poladd(Fn1, dFn1, Fn2, dFn2)
+        Fn = Fn.squeeze()
+        zi = pl.convolve(zi,[0,1])
+        dzi += 1
+        En2 = Fn[0]*zi
+        dEn2 = dzi
+        En, dEn = poladd(Eo, dEo, En2, dEn2)
+        En = En.squeeze()
+        Eo, dEo = En, dEn
+        Fo, dFo = Fn, dFn
+
+    if dFo == 0:
+        Fn1 = 0
+    else:
+        Fn1 = Fo[1:dFo+1]
+
+    dFn1 = max(dFo-1, 0)
+    Fn2 = -Fo[0]*A_z
+    dFn2 = dA-1
+    Fn, dFn = poladd(Fn1, dFn1, Fn2, dFn2)
+    Fn = Fn.squeeze()
+    Ej, dEj = Eo, dEo
+    Fj, dFj = Fn, dFn
+
+
+    return Fj, dFj, Ej, dEj