Index: trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py
===================================================================
--- trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py	(revision 8359)
+++ trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py	(revision 8375)
@@ -144,4 +144,6 @@
         # self.check_integrity()
 
+        #print 'Into Evolve'
+
         msg = 'Attribute self.beta_w must be in the interval [0, 2]'
         assert 0 <= self.beta_w <= 2.0, msg
@@ -155,9 +157,10 @@
         if self.store is True and self.get_time() == self.get_starttime():
             self.initialise_storage()
-
+        #print 'Into Generic_Domain Evolve'
         # Call basic machinery from parent class
         for t in Generic_Domain.evolve(self, yieldstep=yieldstep,
                                        finaltime=finaltime, duration=duration,
                                        skip_initial_step=skip_initial_step):
+            #print 'Out of Generic_Domain Evolve'
             # Store model data, e.g. for subsequent visualisation
             if self.store is True:
@@ -208,4 +211,6 @@
         from swb2_domain_ext import compute_fluxes_ext_central \
                                       as compute_fluxes_ext
+
+        #print "."
 
         # Shortcuts
@@ -364,5 +369,7 @@
 
         # Remove very thin layers of water
+        #print 'Before protect'
         domain.protect_against_infinitesimal_and_negative_heights()
+        #print 'After protect'
 
         #for name in domain.conserved_quantities:
@@ -378,5 +385,7 @@
         #        raise Exception('Unknown order')
         
+        #print 'Before extrapolate'
         domain.extrapolate_second_order_edge_sw()
+        #print 'After extrapolate'
 
         #balance_deep_and_shallow(domain)
Index: trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c
===================================================================
--- trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c	(revision 8359)
+++ trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c	(revision 8375)
@@ -787,7 +787,16 @@
   double dqv[3], qmin, qmax, hmin, hmax, bedmax, stagemin;
   double hc, h0, h1, h2, beta_tmp, hfactor;
-  double xmom_centroid_store[number_of_elements], ymom_centroid_store[number_of_elements], dk, dv0, dv1, dv2, de[3];
-  double stage_centroid_store[number_of_elements];
-    
+  double dk, dv0, dv1, dv2, de[3];
+  //double xmom_centroid_store[number_of_elements], ymom_centroid_store[number_of_elements]
+  //double stage_centroid_store[number_of_elements];
+  
+  double *xmom_centroid_store, *ymom_centroid_store, *stage_centroid_store;
+
+  // Use malloc to avoid putting these variables on the stack, which can cause
+  // segfaults in large model runs
+  xmom_centroid_store = malloc(number_of_elements*sizeof(double));
+  ymom_centroid_store = malloc(number_of_elements*sizeof(double));
+  stage_centroid_store = malloc(number_of_elements*sizeof(double));
+ 
   if(extrapolate_velocity_second_order==1){
       // Replace momentum centroid with velocity centroid to allow velocity
@@ -819,4 +828,5 @@
 
     if (number_of_boundaries[k]==3)
+    //if (0==0)
     {
       // No neighbours, set gradient on the triangle to zero
@@ -1366,4 +1376,8 @@
       } 
   }
+
+  free(xmom_centroid_store);
+  free(ymom_centroid_store);
+  free(stage_centroid_store);
 
   return 0;
@@ -1859,4 +1873,6 @@
   number_of_elements = stage_centroid_values -> dimensions[0];  
 
+  //printf("In C before Extrapolate");
+  //e=1;
   // Call underlying computational routine
   e = _extrapolate_second_order_edge_sw(number_of_elements,
@@ -1884,4 +1900,5 @@
                    extrapolate_velocity_second_order);
 
+  //printf("In C before edges-to-vertices");
 
   //Extrapolate from edges to vertices
@@ -1896,4 +1913,5 @@
                                       (double *) elevation_vertex_values -> data, 
                                       extrapolate_velocity_second_order);
+  //printf("In C after edges-to-vertices");
 
   if (e == -1) {
Index: trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py	(revision 8359)
+++ trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py	(revision 8375)
@@ -4,5 +4,5 @@
 
 # Time-index to plot outputs from
-index=75
+index=900
 
 #p2 = util.get_output('channel_floodplain1_bal_dev.sww', minimum_allowed_height=0.01)
Index: trunk/anuga_work/development/gareth/tests/util.py
===================================================================
--- trunk/anuga_work/development/gareth/tests/util.py	(revision 8359)
+++ trunk/anuga_work/development/gareth/tests/util.py	(revision 8375)
@@ -73,32 +73,4 @@
 ##############
 
-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,2 +170,91 @@
 #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
