--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lqg/lqg_julien1_loop.py	Fri May 27 14:24:59 2011 +0530
@@ -0,0 +1,36 @@
+#!/usr/bin/python
+
+import scipy as sp
+from lqg import lqg
+from cl import cl
+
+A, dA = sp.convolve([1, -0.7], [1, 0.9]), 2
+B, dB = sp.convolve([0.079, 0.221], [1, 0.9]), 2
+C, dC = sp.array([1, -0.7]), 1
+k = 2
+int1 = 1
+F, dF = sp.array([1, -1]), 1
+V, dV = 1, 0
+W, dW = 1, 0
+u_lqg = []
+y_lqg =[]
+uy_lqg = []
+no_points = 101
+rhovector = sp.logspace(-1.63, 1.2, no_points)
+
+for i, rho in enumerate(rhovector):
+    print i, rho
+    R1, dR1, Sc, dSc = lqg (A, dA, B, dB, C, dC, k, rho, V, dV, W, dW, F, dF)
+    print R1, dR1, Sc, dSc
+    Nu, dNu, Du, dDu, Ny, dNy, Dy, dDy, yvar, uvar = \
+        cl(A, dA, B, dB, C, dC, k, Sc, dSc, R1, dR1, int1)
+    uvar = sp.atleast_1d(uvar.squeeze())
+    yvar = sp.atleast_1d(yvar.squeeze())
+    u_lqg = sp.concatenate((u_lqg, uvar))
+    y_lqg = sp.concatenate((y_lqg, yvar))
+    if sp.size(uy_lqg) == 0:
+        uy_lqg = sp.array([rho, uvar[0], yvar[0]])
+    else:
+        uy_lqg = sp.vstack((uy_lqg, sp.array([rho, uvar[0], yvar[0]])))
+
+plot(u_lqg, y_lqg)