Changeset 8375 for trunk/anuga_work/development/gareth/tests/util.py

Timestamp:

Mar 31, 2012, 2:47:12 PM (13 years ago)

Author:

davies

Message:

balanced_dev: fixed stack overflow problem, & new 'transect plot'
routine in util.py

File:

• trunk/anuga_work/development/gareth/tests/util.py (modified) (2 diffs)

trunk/anuga_work/development/gareth/tests/util.py

@@ -73 +73 @@
 ##############
 
-def animate_1D(time, var, x, ylab=' '): #, x=range(var.shape[1]), vmin=var.min(), vmax=var.max()):
-    # Input: time = one-dimensional time vector;
-    #        var =  array with first dimension = len(time) ;
-    #        x = (optional) vector width dimension equal to var.shape[1];
-    
-    import pylab
-    import numpy
-   
-    
-
-    pylab.close()
-    pylab.ion()
-
-    # Initial plot
-    vmin=var.min()
-    vmax=var.max()
-    line, = pylab.plot( (x.min(), x.max()), (vmin, vmax), 'o')
-
-    # Lots of plots
-    for i in range(len(time)):
-        line.set_xdata(x)
-        line.set_ydata(var[i,:])
-        pylab.draw()
-        pylab.xlabel('x')
-        pylab.ylabel(ylab)
-        pylab.title('time = ' + str(time[i]))
-    
-    return
 
 
@@ -198 +170 @@
 #pylab.quiver(x,y,xvel[n,:],yvel[n,:])
 #pylab.savefig('Plot1.png')
+
+def animate_1D(time, var, x, ylab=' '): #, x=range(var.shape[1]), vmin=var.min(), vmax=var.max()):
+    # Input: time = one-dimensional time vector;
+    #        var =  array with first dimension = len(time) ;
+    #        x = (optional) vector width dimension equal to var.shape[1];
+    
+    import pylab
+    import numpy
+   
+    
+
+    pylab.close()
+    pylab.ion()
+
+    # Initial plot
+    vmin=var.min()
+    vmax=var.max()
+    line, = pylab.plot( (x.min(), x.max()), (vmin, vmax), 'o')
+
+    # Lots of plots
+    for i in range(len(time)):
+        line.set_xdata(x)
+        line.set_ydata(var[i,:])
+        pylab.draw()
+        pylab.xlabel('x')
+        pylab.ylabel(ylab)
+        pylab.title('time = ' + str(time[i]))
+    
+    return
+
+def near_transect(p, point1, point2, tol=1.):
+    # Function to get the indices of points in p less than 'tol' from the line
+    # joining (x1,y1), and (x2,y2)
+    # p comes from util.get_output('mysww.sww')
+    #
+    # e.g.
+    # import util
+    # #import transect_interpolate
+    # from matplotlib import pyplot
+    # p=util.get_output('merewether_1m.sww',0.01)
+    # p2=util.get_centroids(p,velocity_extrapolation=True)
+    # #xxx=transect_interpolate.near_transect(p,[95., 85.], [120.,68.],tol=2.)
+    # xxx=util.near_transect(p,[95., 85.], [120.,68.],tol=2.)
+    # pyplot.scatter(xxx[1],p.vel[140,xxx[0]],color='red')
+
+    x1=point1[0]
+    y1=point1[1]
+
+    x2=point2[0]
+    y2=point2[1]
+
+    # Find line equation a*x + b*y + c = 0
+    # based on y=gradient*x +intercept
+    if x1!=x2:
+        gradient= (y2-y1)/(x2-x1)
+        intercept = y1 - gradient*x1
+
+        a = -gradient
+        b = 1.
+        c = -intercept
+    else:
+        #print 'FIXME: Still need to treat 0 and infinite gradients'
+        #assert 0==1
+        a=1.
+        b=0.
+        c=-x2 
+
+    # Distance formula
+    inv_denom = 1./(a**2 + b**2)**0.5
+    distp = abs(p.x*a + p.y*b + c)*inv_denom
+
+    near_points = (distp<tol).nonzero()[0]
+    
+    # Now find a 'local' coordinate for the point, projected onto the line
+    # g1 = unit vector parallel to the line
+    # g2 = vector joining (x1,y1) and (p.x,p.y)
+    g1x = x2-x1 
+    g1y = y2-y1
+    g1_norm = (g1x**2 + g1y**2)**0.5
+    g1x=g1x/g1_norm
+    g1y=g1x/g1_norm
+
+    g2x = p.x[near_points] - x1
+    g2y = p.y[near_points] - y1
+    
+    # Dot product = projected distance == a local coordinate
+    local_coord = g1x*g2x + g1y*g2y
+
+    return near_points, local_coord