diff --git a/EBSDAnalysis/@EBSD/mean.m b/EBSDAnalysis/@EBSD/mean.m
deleted file mode 100644
index aad77cefa..000000000
--- a/EBSDAnalysis/@EBSD/mean.m
+++ /dev/null
@@ -1,18 +0,0 @@
-function varargout = mean( ebsd,varargin)
-% returns mean, kappas and eigenvector of ebsd object
-%
-% Syntax
-% [m lambda v kappa] = mean(ebsd) - 
-%
-% Input
-%  ebsd     - @EBSD
-%
-% Output
-%  m        - one equivalent mean @orientation
-%  lambda   - eigenvalues of orientation tensor
-%  v        - eigenvectors of orientation tensor
-%  kappa    - parameters of bingham distribution
-%
-
-[varargout{1:nargout}]  = mean(ebsd.orientations,'weights',ebsd.weights,varargin{:});
-
diff --git a/EBSDAnalysis/@EBSD/private/subsind.m b/EBSDAnalysis/@EBSD/private/subsind.m
index 70f50acb1..a1bde5010 100644
--- a/EBSDAnalysis/@EBSD/private/subsind.m
+++ b/EBSDAnalysis/@EBSD/private/subsind.m
@@ -67,9 +67,8 @@
     
   elseif isnumeric(subs{i})
     
-    iind = false(size(ind));
-    iind(subs{i}) = true;
-    ind = ind & iind;
+    ind = subs{i};
+    return
     
   end
 end
diff --git a/EBSDAnalysis/@grain2d/private/subsind.m b/EBSDAnalysis/@grain2d/private/subsind.m
index 450eeb259..3bec9b8bd 100644
--- a/EBSDAnalysis/@grain2d/private/subsind.m
+++ b/EBSDAnalysis/@grain2d/private/subsind.m
@@ -38,10 +38,9 @@
       error('Out of range; index must be a positive integer or logical.')
     end
     
-    iind = false(size(ind));
-    iind(subs{i}) = true;
-    ind = ind & iind;
-    
+    ind = subs{i};
+    return
+        
   elseif isa(subs{i},'polygon')
     
     ind = ind & inpolygon(grains,subs{i})';
diff --git a/EBSDAnalysis/@grainBoundary/intersect.m b/EBSDAnalysis/@grainBoundary/intersect.m
new file mode 100644
index 000000000..c4d987a09
--- /dev/null
+++ b/EBSDAnalysis/@grainBoundary/intersect.m
@@ -0,0 +1,73 @@
+function  [x,y,segLength] = intersect(gB,xy1,xy2,varargin)
+% length of a boundary segment
+%
+% Syntax
+%   [x,y] = intersect(gB,xy1,xy2)
+%   [x,y,segLength] = intersect(gB,xy1,xy2)
+%
+% Input
+%  gb - @grainBoundary
+%  xy1, xy2 coordinates of the endpoints of the line
+%
+% Output
+%  x,y - list of intersection points
+%
+% Example
+%  mtexdata csl
+%  grains = calcGrains(ebsd)
+%  plot(grains.boundary)
+%  % define some line
+%  xy1 = [0,10];  % staring point
+%  xy2 = [31,41]; % end point
+%  line([xy1(1);xy2(1)],[xy1(2);xy2(2)],'linewidth',1.5,'color','g')
+%  [x,y] = grains.boundary.intersect(xy1,xy2);
+%  hold on
+%  scatter(x,y,'red')
+%  hold off
+%  % find the number of intersection points  
+%  sum(~isnan(x))
+
+n_rows_1 = size(xy1,1);
+n_rows_2 = length(gB);
+
+% end points of the lines
+X1 = repmat(xy1(:,1),1,n_rows_2);
+Y1 = repmat(xy1(:,2),1,n_rows_2);
+X2 = repmat(xy2(:,1),1,n_rows_2);
+Y2 = repmat(xy2(:,2),1,n_rows_2);
+
+% end points boundary segments
+X3 = repmat(gB.V(gB.F(:,1),1).',n_rows_1,1);
+Y3 = repmat(gB.V(gB.F(:,1),2).',n_rows_1,1);
+X4 = repmat(gB.V(gB.F(:,2),1).',n_rows_1,1);
+Y4 = repmat(gB.V(gB.F(:,2),2).',n_rows_1,1);
+
+X4_X3 = X4-X3;
+Y1_Y3 = Y1-Y3;
+Y4_Y3 = Y4-Y3;
+X1_X3 = X1-X3;
+X2_X1 = X2-X1;
+Y2_Y1 = Y2-Y1;
+
+numerator_a = X4_X3 .* Y1_Y3 - Y4_Y3 .* X1_X3;
+numerator_b = X2_X1 .* Y1_Y3 - Y2_Y1 .* X1_X3;
+denominator = Y4_Y3 .* X2_X1 - X4_X3 .* Y2_Y1;
+
+u_a = numerator_a ./ denominator;
+u_b = numerator_b ./ denominator;
+inside = (u_a >= 0) & (u_a <= 1) & (u_b >= 0) & (u_b <= 1);
+
+% Find the adjacency matrix A of intersecting lines.
+x = X1 + X2_X1 .* u_a;
+y = Y1 + Y2_Y1 .* u_a;
+x(~inside) = NaN;
+y(~inside) = NaN;
+
+% sort by distance to starting point
+d = sqrt((x(~isnan(x))-xy1(1)).^2 + (y(~isnan(x))-xy1(1)).^2);
+[~,ind] = sort(d);
+
+segLength = diff(d(ind));
+
+
+
diff --git a/EBSDAnalysis/@grainBoundary/private/subsind.m b/EBSDAnalysis/@grainBoundary/private/subsind.m
index 6a758a81c..3c9313cde 100644
--- a/EBSDAnalysis/@grainBoundary/private/subsind.m
+++ b/EBSDAnalysis/@grainBoundary/private/subsind.m
@@ -37,9 +37,8 @@
       error('Out of range; index must be a positive integer or logical.')
     end
     
-    iind = false(size(ind));
-    iind(subs{i}) = true;
-    ind = ind & iind;
+    ind = subs{i};
+    return
     
   elseif isa(subs{i},'polygon')
     
diff --git a/Makefile b/Makefile
index 4fd5c9321..b3e7bd7a6 100644
--- a/Makefile
+++ b/Makefile
@@ -60,7 +60,7 @@ clean:
 
 
 # rule for making release
-RNAME = mtex-4.0.18
+RNAME = mtex-4.0.20
 RDIR = ../releases
 release:
 	rm -rf $(RDIR)/$(RNAME)*
diff --git a/PoleFigureAnalysis/@PoleFigure/cat.m b/PoleFigureAnalysis/@PoleFigure/cat.m
index 3d1948e7b..141728a60 100644
--- a/PoleFigureAnalysis/@PoleFigure/cat.m
+++ b/PoleFigureAnalysis/@PoleFigure/cat.m
@@ -4,6 +4,8 @@
 % concatenate properties
 pf = cat@dynProp(dim,varargin{:});
 
+varargin(cellfun(@isempty,varargin)) = []; 
+
 warning('off','MATLAB:structOnObject');
 for k=1:numel(varargin)
   s(k) = struct(varargin{k}); %#ok<AGROW>
diff --git a/VERSION b/VERSION
index 2e9571a51..4476a33ee 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-MTEX 4.0.18
+MTEX 4.0.20
diff --git a/geometry/@sphericalRegion/eq.m b/geometry/@sphericalRegion/eq.m
index e56e0c9d5..5fc543a27 100644
--- a/geometry/@sphericalRegion/eq.m
+++ b/geometry/@sphericalRegion/eq.m
@@ -5,7 +5,7 @@
 sR2 = cleanUp(sR2);
 
 % number of restrictions should be the same
-if numel(sR1) ~= numel(sR2), return; end
+if length(sR1.N) ~= length(sR2.N), return; end
 
 % normals should fit
 fit = dot_outer(sR1.N,sR2.N)>1-1e-6;
diff --git a/geometry/@sphericalRegion/sphericalRegion.m b/geometry/@sphericalRegion/sphericalRegion.m
index 7cf1ee306..4af1a6f88 100644
--- a/geometry/@sphericalRegion/sphericalRegion.m
+++ b/geometry/@sphericalRegion/sphericalRegion.m
@@ -1,11 +1,11 @@
 classdef sphericalRegion
-  %sphericalRegion implements a region region on the sphere
-  %   The region is bounded by small circles given by there normal vectors
-  %   and the maximum inner product, i.e., all points v inside a region
-  %   satisfy the conditions dot(v, N) <= alpha
+  %sphericalRegion implements a region on the sphere
+  % The region is bounded by small circles given by there normal vectors
+  % and the maximum inner product, i.e., all points v inside a region
+  % satisfy the conditions dot(v, N) <= alpha
   
   properties
-    N = vector3d     % the nornal vectors of the bounding circles
+    N = vector3d     % the normal vectors of the bounding circles
     alpha = [] % the cosine of the bounding circle
     antipodal = false
   end
diff --git a/geometry/@symmetry/fundamentalSector.m b/geometry/@symmetry/fundamentalSector.m
index da42f9a8c..c87431677 100644
--- a/geometry/@symmetry/fundamentalSector.m
+++ b/geometry/@symmetry/fundamentalSector.m
@@ -7,10 +7,10 @@
 %   sector of rotational axes for a specific rotation angle omega
 %
 % Input
-%  cs - symmetry
+%  cs - @symmetry
 %
 % Ouput
-%  sR - spherical Region
+%  sR - @sphericalRegion
 %
 % Options
 %  antipodal - include [[AxialDirectional.html,antipodal symmetry]]
diff --git a/geometry/@vector3d/scatter.m b/geometry/@vector3d/scatter.m
index 93e7586ef..8d13fdc1d 100644
--- a/geometry/@vector3d/scatter.m
+++ b/geometry/@vector3d/scatter.m
@@ -35,6 +35,16 @@
   % check that there is something left to plot
   if all(isnan(x) | isnan(y)), continue; end
     
+  % add some nans if lines are plotted
+  if check_option(varargin,'edgecolor')
+    d = sqrt(diff(x).^2 + diff(y).^2);
+    ind = find(d > 2.5);
+    for k = 1:numel(ind)
+      x = [x(1:ind(k)+k-1);nan;x(ind(k)+k:end)];
+      y = [y(1:ind(k)+k-1);nan;y(ind(k)+k:end)];
+    end
+  end
+  
   % default arguments
   patchArgs = {'parent',sP(i).ax,...
     'vertices',[x(:) y(:)],...
diff --git a/interfaces/loadPoleFigure_beartex.m b/interfaces/loadPoleFigure_beartex.m
index b61b2a679..b074d6d2b 100644
--- a/interfaces/loadPoleFigure_beartex.m
+++ b/interfaces/loadPoleFigure_beartex.m
@@ -2,17 +2,18 @@
 % import data fom BeaTex file
 %
 % Syntax
-%   pf = loadPoleFigure_beartex(fname,<options>)
+%   pf = loadPoleFigure_beartex(fname)
 %
 % Input
-%  fname    - filename
+%  fname - filename
 %
 % Output
-%  pf - vector of @PoleFigure
+%  pf - @PoleFigure
 %
 % See also
 % ImportPoleFigureData loadPoleFigure
 
+pf = PoleFigure;
 fid = efopen(fname);
 
 ipf = 1;
@@ -27,7 +28,7 @@
       % c = textscan(fid,'%s',7,'delimiter','\n','whitespace','');
       comment = deblank(c{1}(1:50));
     catch
-      if ~exist('pf','var')
+      if ~exist('data','var')
         error('format BearTex does not match file %s',fname);
       else
         break
@@ -43,6 +44,7 @@
     hkl = sscanf(c{7},'%f',3);
     h{ipf} = Miller(hkl(1),hkl(2),hkl(3),cs);
     
+       
     info = str2num(reshape(c{7}(11:40),5,[])');
     
     %  theta = 0:info(3):90-info(3);
@@ -53,26 +55,24 @@
       l = fgetl(fid);
       data{ipf}(:,k) = str2num( reshape(l(2:end),4,[]).' );
     end
+  
+    % restrict to the mesured region
+    ind = r.theta < info(1)*degree-eps | r.theta > info(2)*degree+eps;
+    allR{ipf} = r(~ind);
+    data{ipf}(ind) = [];
     
     fgetl(fid);
-        
-    
-    % mintheta = info(1);  maxtheta = info(2);
-        
-    
+            
+    % mintheta = info(1);  maxtheta = info(2);            
     ipf = ipf+1;
   end
   
-  pf = PoleFigure(h,r,data,cs,ss,'comment',comment,varargin{:});
+  pf = PoleFigure(h,allR,data,cs,ss,'comment',comment,varargin{:});
+
+  fclose(fid);
   
 catch
   if ~exist('pf','var')
     interfaceError(fname,fid);
   end
 end
-
-pf(pf.r.theta < info(1)*degree-eps | pf.r.theta > info(2)*degree+eps) = [];
-
-
-fclose(fid);
-
diff --git a/plotting/sphericalPlot.m b/plotting/sphericalPlot.m
index 308117a18..930ea7688 100644
--- a/plotting/sphericalPlot.m
+++ b/plotting/sphericalPlot.m
@@ -1,5 +1,5 @@
 classdef sphericalPlot < handle
-  %sphericalProjection
+  % sphericalPlot is responsible for visualizing spherical data  
   
   properties
     proj = sphericalProjection
@@ -163,7 +163,6 @@ function doGridInFront(sP)
         isgrid = ismember(childs,[sP.grid(:);sP.boundary(:)]);
         istext = strcmp(get(childs,'type'),'text');
   
-        % TODO: this crahes on MATLAB 2014b
         set(sP.ax,'Children',[childs(istext); sP.boundary(:); sP.grid(:);childs(~isgrid & ~istext)]);
       end
     end    
@@ -174,12 +173,9 @@ function doGridInFront(sP)
     function plotPlainGrid(sP,varargin)
       
       % the ticks
-      %dgrid = get_option(varargin,'grid_res',30*degree);
       polarRange = sP.sphericalRegion.polarRange;
       theta = round(linspace(polarRange(1),polarRange(3),4)/degree);
-      rho = round(linspace(polarRange(2),polarRange(4),4)/degree);      
-      %theta = round((polarRange(1):dgrid:polarRange(3))/degree);
-      %rho = round((polarRange(2):dgrid:polarRange(4))/degree);
+      rho = round(linspace(polarRange(2),polarRange(4),4)/degree);            
       
       set(sP.ax,'XTick',rho);
       set(sP.ax,'YTick',theta);
@@ -205,8 +201,6 @@ function plotPolarGrid(sP,varargin)
       for i = 1:length(theta), circ(sP,theta(i)); end
       
       % draw meridians
-      %rho = 0:dgrid:2*pi-dgrid;
-      %for i = 1:length(rho), plotMeridian(sP,rho(i)); end
       plotMeridians(sP,0:dgrid:pi-dgrid);
 
     end
@@ -247,29 +241,3 @@ function circ(sP,theta,varargin)
   end
 end
 
-% 
-% 
-% % control legend entry
-% try
-%   hAnnotation = get(l,'Annotation');
-%   hLegendEntry = get([hAnnotation{:}],'LegendInformation');
-%   set([hLegendEntry{:}],'IconDisplayStyle','off')
-% catch %#ok<CTCH>
-% end
-% 
-% % labels
-% 
-% 
-% if any(isnan(X)), return;end
-% if check_option(varargin,'ticks'), v = 'on';else v = 'off';end
-% 
-% % set back color index
-% if isappdata(gca,'PlotColorIndex')
-%   if isempty(colorIndex)
-%     setappdata(gca,'PlotColorIndex',1);
-%   else
-%     setappdata(gca,'PlotColorIndex',colorIndex);
-%   end
-% end
-% 
-% end