diff --git a/scripts/extract_on_line.m b/scripts/extract_on_line.m
index 55b0d074dd6f5cb7f4b2c6b050fa1aedeb928e75..c7abf946cf4118dc608bc2097d598e3fd1b91e21 100644
--- a/scripts/extract_on_line.m
+++ b/scripts/extract_on_line.m
@@ -1,92 +1,95 @@
-function [ cell_points_on, p_on, cell_points_nw_on, cell_face_on,p_on2 ] ...
- = extract_on_line( t, x, nw,on, p )
-%UNTITLED4 returns cell_points (pointers to x) and p for each cell
-%on the line. cell_on_line is the logical vector s.t. p(c_o_l)=p_on.
-%cell_points_nw_on has a cell for each grid cell with the logical information
-%whether the cell_points_nw_on{c}(i) pointed to by cell_points{c}(i) lies nw
-% NB for generality, t is assumed to be a cell. Note that in the case
-% that a face lies on the line, a matrix is filled, because we know the
-% exact number of points (2) in this case.
-
-count=1;
-count2=1;
-
-cell_points_on = cell(0,1);
-cell_points_nw_on = cell(0,1);
-cell_face_on = zeros(0,4);
-p_on = zeros(0,1);
-p_on2 = zeros(0,1);
-noton = ~on;
-if iscell(t);
-
- for c = 1:length(t)
- cells = t{c};
- points_nw = nw(cells);
- cell_on_line = not(all(points_nw) | all(not(points_nw)));
- vert_on = on(cells);
- % two vertices on the line
- if sum(vert_on) == 2
- vert_not = noton(cells);
- not_on_nw = points_nw(vert_not);
- % check that it is really a face, not e.g. point 1 and 3 of a
- % quadrilateral
- if all(not_on_nw) || all(~not_on_nw)
- ind = find(vert_on);
- cell_face_on(count2,1:2) = x(cells(ind(1)),:);
- cell_face_on(count2,3:4) = x(cells(ind(2)),:);
- p_on2(count2,1) = p(c);
- count2 = count2+1;
- %else it should be included in the main cell:
- else
- cell_points_on{count} = x(cells,:);
- cell_points_nw_on{count} = points_nw;
- p_on(count,1) = p(c);
- count = count + 1;
- end
- % cells intersected by the line, not two vertices on it
- elseif cell_on_line
- cell_points_on{count} = x(cells,:);
- cell_points_nw_on{count} = points_nw;
- p_on(count,1) = p(c);
- count = count + 1;
-
- end
- end
-else
- for c = 1:length(t)
- cells = t(c,:);
- points_nw = nw(cells);
- cell_on_line = not(all(points_nw) | all(not(points_nw)));
- vert_on = on(cells);
- % two vertices on the line
- if sum(vert_on) == 2
- vert_not = noton(cells);
- not_on_nw = points_nw(vert_not);
- % check that it is really a face, not e.g. point 1 and 3 of a
- % quadrilateral
- if all(not_on_nw) || all(~not_on_nw)
- ind = find(vert_on);
- cell_face_on(count2,1:2) = x(cells(ind(1)),:);
- cell_face_on(count2,3:4) = x(cells(ind(2)),:);
- p_on2(count2,1) = p(c);
- count2 = count2+1;
- %else it should be included in the main cell:
- else
- cell_points_on{count} = x(cells,:);
- cell_points_nw_on{count} = points_nw;
- p_on(count,1) = p(c);
- count = count + 1;
- end
- % cells intersected by the line, not two vertices on it
- elseif cell_on_line
- cell_points_on{count} = x(cells,:);
- cell_points_nw_on{count} = points_nw;
- p_on(count,1) = p(c);
- count = count + 1;
-
- end
- end
-end
-
-end
-
+function [ cellPointsOn, pOn, cellPointsNorthWestOn, cellFaceOn, pOn2 ] ...
+ = find_on_line( t, x, northWest,on, p )
+%find_on_line returns cellPoints (pointers to x) and p for each cell
+%on the line.
+%cellOnLine is the logical vector s.t. p(c_o_l)=p_on.
+%cellPointsNorthWestOn has a cell for each grid cell with the logical information
+%whether the cellPointsNorthWestOn{c}(i) pointed to by cell_points{c}(i) lies
+%nw of the line. This (and the two last outputs) are given to avoid performing
+%calculations multiple times.
+% NB for generality, t is assumed to be a cell. Note that in the case
+% that a face lies on the line, a matrix is filled, because we know the
+% exact number of points (2) in this case.
+
+count=1;
+count2=1;
+
+cellPointsOn = cell(0,1);
+cellPointsNorthWestOn = cell(0,1);
+cellFaceOn = zeros(0,4);
+pOn = zeros(0,1);
+pOn2 = zeros(0,1);
+notOn = ~on;
+if iscell(t);
+
+ for c = 1:length(t)
+ cells = t{c};
+ pointsNorthWest = northWest(cells);
+ cellOnLine = not(all(pointsNorthWest) || all(not(pointsNorthWest)));
+ verticesOn = on(cells);
+ % two vertices on the line
+ if sum(verticesOn) == 2
+ verticesNot = notOn(cells);
+ notOnNorthWest = pointsNorthWest(verticesNot);
+ % check that it is really a face, not e.g. point 1 and 3 of a
+ % quadrilateral
+ if all(notOnNorthWest) || all(~notOnNorthWest)
+ ind = find(verticesOn);
+ cellFaceOn(count2,1:2) = x(cells(ind(1)),:);
+ cellFaceOn(count2,3:4) = x(cells(ind(2)),:);
+ pOn2(count2,1) = p(c);
+ count2 = count2+1;
+ %else it should be included in the main cell:
+ else
+ cellPointsOn{count} = x(cells,:);
+ cellPointsNorthWestOn{count} = pointsNorthWest;
+ pOn(count,1) = p(c);
+ count = count + 1;
+ end
+ % cells intersected by the line, not two vertices on it
+ elseif cellOnLine
+ cellPointsOn{count} = x(cells,:);
+ cellPointsNorthWestOn{count} = pointsNorthWest;
+ pOn(count,1) = p(c);
+ count = count + 1;
+
+ end
+ end
+else
+ for c = 1:length(t)
+ cells = t(c,:);
+ pointsNorthWest = northWest(cells);
+ cellOnLine = not(all(pointsNorthWest) | all(not(pointsNorthWest)));
+ verticesOn = on(cells);
+ % two vertices on the line
+ if sum(verticesOn) == 2
+ verticesNot = notOn(cells);
+ notOnNorthWest = pointsNorthWest(verticesNot);
+ % check that it is really a face, not e.g. point 1 and 3 of a
+ % quadrilateral
+ if all(notOnNorthWest) || all(~notOnNorthWest)
+ ind = find(verticesOn);
+ cellFaceOn(count2,1:2) = x(cells(ind(1)),:);
+ cellFaceOn(count2,3:4) = x(cells(ind(2)),:);
+ pOn2(count2,1) = p(c);
+ count2 = count2+1;
+ %else it should be included in the main cell:
+ else
+ cellPointsOn{count} = x(cells,:);
+ cellPointsNorthWestOn{count} = pointsNorthWest;
+ pOn(count,1) = p(c);
+ count = count + 1;
+ end
+ % cells intersected by the line, not two vertices on it
+ elseif cellOnLine
+ cellPointsOn{count} = x(cells,:);
+ cellPointsNorthWestOn{count} = pointsNorthWest;
+ pOn(count,1) = p(c);
+ count = count + 1;
+
+ end
+ end
+end
+
+end
+