diff --git a/Contents.m b/Contents.m
index b08fdef..7bfcb86 100644
--- a/Contents.m
+++ b/Contents.m
@@ -60,4 +60,3 @@
 %
 % Author:
 %   Philipp Berens & Marc J. Velasco, 2009
-
diff --git a/circ_ang2rad.m b/circ_ang2rad.m
index 215fcbb..0bb1677 100644
--- a/circ_ang2rad.m
+++ b/circ_ang2rad.m
@@ -8,4 +8,5 @@
 % By Philipp Berens, 2009
 % berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
 
-alpha = alpha * pi /180;
\ No newline at end of file
+alpha = alpha * pi /180;
+end
diff --git a/circ_axial.m b/circ_axial.m
index 57378b3..1cccc51 100644
--- a/circ_axial.m
+++ b/circ_axial.m
@@ -27,3 +27,4 @@
 end
 
 alpha = mod(alpha*p,2*pi);
+end
diff --git a/circ_axialmean.m b/circ_axialmean.m
index ce1f60f..dce22a2 100644
--- a/circ_axialmean.m
+++ b/circ_axialmean.m
@@ -1,4 +1,4 @@
-function [r mu] = circ_axialmean(alphas, m, dim)
+function [r, mu] = circ_axialmean(alphas, m, dim)
 %
 % mu = circ_axialmean(alpha, w)
 %   Computes the mean direction for circular data with axial 
@@ -7,7 +7,7 @@
 %   Input:
 %     alpha	sample of angles in radians
 %     [m		axial correction (2,3,4,...)]
-%     [dim      statistic computed along this dimension, 1]
+%     [dim      statistic computed along this dimension, default: 1st non-singular dimension]
 %
 %   Output:
 %     r		mean resultant length
@@ -27,7 +27,10 @@
 % Distributed under Open Source BSD License
 
 if nargin < 3
-  dim = 1;
+  dim = find(size(alphas) > 1, 1, 'first');
+  if isempty(dim)
+    dim = 1;
+  end
 end
 
 if nargin < 2 || isempty(m)
@@ -38,4 +41,4 @@
 
 r = abs(zbarm);
 mu = angle(zbarm)/m;
-
+end
diff --git a/circ_clust.m b/circ_clust.m
index 37d58b4..b95450f 100644
--- a/circ_clust.m
+++ b/circ_clust.m
@@ -138,7 +138,7 @@ function plotColor(x, y, c, varargin)
 colors={'y', 'b', 'r', 'g', 'c', 'k', 'm'};
 
 hold on;
-for j=1:length(csmall);
+for j=1:length(csmall)
     
    ci = (c == csmall(j));
    plot(x(ci), y(ci), strcat(pstring, colors{mod(j, length(colors))+1}), 'MarkerSize', ms);
@@ -146,6 +146,5 @@ function plotColor(x, y, c, varargin)
 end
 if ~overlay, hold off; end
 figure(figurenum)
-
-
-
+end
+end
diff --git a/circ_cmtest.m b/circ_cmtest.m
index 3e19265..9418b7c 100644
--- a/circ_cmtest.m
+++ b/circ_cmtest.m
@@ -1,4 +1,4 @@
-function [pval med P] = circ_cmtest(varargin)
+function [pval, med, P] = circ_cmtest(varargin)
 %
 % [pval, med, P] = circ_cmtest(alpha, idx)
 % [pval, med, P] = circ_cmtest(alpha1, alpha2)
@@ -68,7 +68,7 @@
 if pval < 0.05
   med = NaN;
 end
-
+end
 
 
 
@@ -88,4 +88,4 @@
 else
   error('Invalid use of circ_wwtest. Type help circ_wwtest.')
 end
-
+end
diff --git a/circ_confmean.m b/circ_confmean.m
index f1d0b48..5799c6e 100644
--- a/circ_confmean.m
+++ b/circ_confmean.m
@@ -10,7 +10,7 @@
 %     [d    spacing of bin centers for binned data, if supplied 
 %           correction factor is used to correct for bias in 
 %           estimation of r, in radians (!)]
-%     [dim  compute along this dimension, default is 1]
+%     [dim  compute along this dimension, default: 1st non-singular dimension]
 %
 %   Output:
 %     t     mean +- d yields upper/lower (1-xi)% confidence limit
@@ -28,7 +28,10 @@
 % berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
 
 if nargin < 5
-  dim = 1;
+  dim = find(size(alpha) > 1, 1, 'first');
+  if isempty(dim)
+    dim = 1;
+  end
 end
 
 if nargin < 4 || isempty(d)
@@ -74,7 +77,4 @@
 
 % apply final transform
 t = acos(t./R);
-  
-
-
-
+end
diff --git a/circ_corrcc.m b/circ_corrcc.m
index 9238849..5fc19bc 100644
--- a/circ_corrcc.m
+++ b/circ_corrcc.m
@@ -1,4 +1,4 @@
-function [rho pval] = circ_corrcc(alpha1, alpha2)
+function [rho, pval] = circ_corrcc(alpha1, alpha2)
 %
 % [rho pval] = circ_corrcc(alpha1, alpha2)
 %   Circular correlation coefficient for two circular random variables.
@@ -50,4 +50,4 @@
 
 ts = sqrt((n * l20 * l02)/l22) * rho;
 pval = 2 * (1 - normcdf(abs(ts)));
-
+end
diff --git a/circ_corrcl.m b/circ_corrcl.m
index da5f01a..8bbdcab 100644
--- a/circ_corrcl.m
+++ b/circ_corrcl.m
@@ -1,4 +1,4 @@
-function [rho pval] = circ_corrcl(alpha, x)
+function [rho, pval] = circ_corrcl(alpha, x)
 %
 % [rho pval] = circ_corrcc(alpha, x)
 %   Correlation coefficient between one circular and one linear random
@@ -47,4 +47,4 @@
 
 % compute pvalue
 pval = 1 - chi2cdf(n*rho^2,2);
-
+end
diff --git a/circ_dist.m b/circ_dist.m
index 1ba0eaf..b9a90ae 100644
--- a/circ_dist.m
+++ b/circ_dist.m
@@ -25,4 +25,5 @@
   error('Input dimensions do not match.')
 end
 
-r = angle(exp(1i*x)./exp(1i*y));
\ No newline at end of file
+r = angle(exp(1i*x)./exp(1i*y));
+end
diff --git a/circ_dist2.m b/circ_dist2.m
index 1de72ef..f260475 100644
--- a/circ_dist2.m
+++ b/circ_dist2.m
@@ -33,4 +33,5 @@
 end
 
 r = angle(repmat(exp(1i*x),1,length(y)) ...
-       ./ repmat(exp(1i*y'),length(x),1));
\ No newline at end of file
+       ./ repmat(exp(1i*y'),length(x),1));
+end
diff --git a/circ_hktest.m b/circ_hktest.m
index fe976dc..17c6745 100644
--- a/circ_hktest.m
+++ b/circ_hktest.m
@@ -1,4 +1,4 @@
-function [pval table] = circ_hktest(alpha, idp, idq, inter, fn)
+function [pval, table] = circ_hktest(alpha, idp, idq, inter, fn)
 
 %
 % [pval, stats] = circ_hktest(alpha, idp, idq, inter, fn)
@@ -238,16 +238,4 @@
     end
     
   end
-
-
 end
-
-
-
-
-
-
-
-
-
-
diff --git a/circ_kappa.m b/circ_kappa.m
index 78d3ffc..76852db 100644
--- a/circ_kappa.m
+++ b/circ_kappa.m
@@ -55,3 +55,4 @@
     kappa = (N-1)^3*kappa/(N^3+N);
   end
 end
+end
diff --git a/circ_ktest.m b/circ_ktest.m
index 408e4ab..5568f00 100644
--- a/circ_ktest.m
+++ b/circ_ktest.m
@@ -51,10 +51,4 @@
   f = 1/f; 
   pval = 2*(1-fcdf(f, n2, n1));
 end
-
-
-
-
-
-
-
+end
diff --git a/circ_kuipertest.m b/circ_kuipertest.m
index 9866dca..4669997 100644
--- a/circ_kuipertest.m
+++ b/circ_kuipertest.m
@@ -45,8 +45,8 @@
 m = length(alpha2(:));
 
 % create cdfs of both samples
-[phis1 cdf1 phiplot1 cdfplot1] = circ_samplecdf(alpha1, res);
-[foo, cdf2 phiplot2 cdfplot2] = circ_samplecdf(alpha2, res); %#ok<ASGLU>
+[phis1, cdf1, phiplot1, cdfplot1] = circ_samplecdf(alpha1, res);
+[foo, cdf2, phiplot2, cdfplot2] = circ_samplecdf(alpha2, res); %#ok<ASGLU>
 
 % maximal difference between sample cdfs
 [dplus, gdpi] = max([0 cdf1-cdf2]);
@@ -56,7 +56,7 @@
 k = n * m * (dplus + dminus);
 
 % find p-value
-[pval K] = kuiperlookup(min(n,m),k/sqrt(n*m*(n+m)));
+[pval, K] = kuiperlookup(min(n,m),k/sqrt(n*m*(n+m)));
 K = K * sqrt(n*m*(n+m));
 
 
@@ -83,14 +83,14 @@
 
 end
 
-function [p K] = kuiperlookup(n, k)
+function [p, K] = kuiperlookup(n, k)
 
 load kuipertable.mat;
 alpha = [.10, .05, .02, .01, .005, .002, .001];
 nn = ktable(:,1);  %#ok<NODEF>
 
 % find correct row of the table
-[easy row] = ismember(n, nn);
+[easy, row] = ismember(n, nn);
 if ~easy
    % find closest value if no entry is present)
    row = length(nn) - sum(n<nn); 
@@ -110,5 +110,4 @@
   p = 1;
 end
 K = ktable(row,idx+1);
-
-end
\ No newline at end of file
+end
diff --git a/circ_kurtosis.m b/circ_kurtosis.m
index 4092c02..23d8e97 100644
--- a/circ_kurtosis.m
+++ b/circ_kurtosis.m
@@ -1,4 +1,4 @@
-function [k k0] = circ_kurtosis(alpha, w, dim)
+function [k, k0] = circ_kurtosis(alpha, w, dim)
 
 % [k k0] = circ_kurtosis(alpha,w,dim)
 %   Calculates a measure of angular kurtosis.
@@ -6,7 +6,7 @@
 %   Input:
 %     alpha     sample of angles
 %     [w        weightings in case of binned angle data]
-%     [dim      statistic computed along this dimension, 1]
+%     [dim      statistic computed along this dimension, default: 1st non-singular dimension]
 %
 %     If dim argument is specified, all other optional arguments can be
 %     left empty: circ_kurtosis(alpha, [], dim)
@@ -25,7 +25,10 @@
 % berens@tuebingen.mpg.de
 
 if nargin < 3
-  dim = 1;
+  dim = find(size(alpha) > 1, 1, 'first');
+  if isempty(dim)
+    dim = 1;
+  end  
 end
 
 if nargin < 2 || isempty(w)
@@ -48,4 +51,5 @@
 theta2 = repmat(theta, size(alpha)./size(theta));
 k = sum(w.*(cos(2*(circ_dist(alpha,theta2)))),dim)./sum(w,dim);
 k0 = (rho2.*cos(circ_dist(mu2,2*theta))-R.^4)./(1-R).^2;    % (formula 2.30)
+end
 
diff --git a/circ_mean.m b/circ_mean.m
index cfce55e..c97916f 100644
--- a/circ_mean.m
+++ b/circ_mean.m
@@ -1,4 +1,4 @@
-function [mu ul ll] = circ_mean(alpha, w, dim)
+function [mu, ul, ll] = circ_mean(alpha, w, dim)
 %
 % [mu ul ll] = circ_mean(alpha, w, dim)
 %   Computes the mean direction for circular data.
@@ -6,7 +6,7 @@
 %   Input:
 %     alpha	sample of angles in radians
 %     [w		weightings in case of binned angle data]
-%     [dim  compute along this dimension, default is 1]
+%     [dim  compute along this dimension, default: 1st non-singular dimension]
 %
 %     If dim argument is specified, all other optional arguments can be
 %     left empty: circ_mean(alpha, [], dim)
@@ -29,7 +29,10 @@
 % berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
 
 if nargin < 3
-  dim = 1;
+  dim = find(size(alpha) > 1, 1, 'first');
+  if isempty(dim)
+    dim = 1;
+  end
 end
 
 if nargin < 2 || isempty(w)
@@ -53,4 +56,5 @@
   t = circ_confmean(alpha,0.05,w,[],dim);
   ul = mu + t;
   ll = mu - t;
-end
\ No newline at end of file
+end
+end
diff --git a/circ_median.m b/circ_median.m
index 72b03f5..aaca36f 100644
--- a/circ_median.m
+++ b/circ_median.m
@@ -5,7 +5,7 @@
 %
 %   Input:
 %     alpha	sample of angles in radians
-%     [dim  compute along this dimension, default is 1, must 
+%     [dim  compute along this dimension, default: 1st non-singular dimension, must 
 %           be either 1 or 2 for circ_median]
 %
 %   Output:
@@ -25,7 +25,10 @@
 % berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
 
 if nargin < 2
-  dim = 1;
+  dim = find(size(alpha) > 1, 1, 'first');
+  if isempty(dim)
+    dim = 1;
+  end  
 end
 
 M = size(alpha);
@@ -48,7 +51,7 @@
 
   dm = abs(m1-m2);
   if mod(n,2)==1
-    [m idx] = min(dm);
+    [m, idx] = min(dm);
   else
     m = min(dm);
     idx = find(dm==m,2);
@@ -70,4 +73,5 @@
 
 if dim == 2
   med = med';
-end
\ No newline at end of file
+end
+end
diff --git a/circ_medtest.m b/circ_medtest.m
index 2b74450..a846886 100644
--- a/circ_medtest.m
+++ b/circ_medtest.m
@@ -40,7 +40,4 @@
 
 % compute p-value with binomial test
 pval = sum(binopdf([0:min(n1,n2) max(n1,n2):n],n,0.5));
-
-
-
-
+end
diff --git a/circ_moment.m b/circ_moment.m
index f930426..beb98c9 100644
--- a/circ_moment.m
+++ b/circ_moment.m
@@ -1,4 +1,4 @@
-function [mp  rho_p mu_p] = circ_moment(alpha, w, p, cent, dim)
+function [mp,  rho_p, mu_p] = circ_moment(alpha, w, p, cent, dim)
 
 % [mp rho_p mu_p] = circ_moment(alpha, w, p, cent, dim)
 %   Calculates the complex p-th centred or non-centred moment 
@@ -9,7 +9,7 @@
 %     [w        weightings in case of binned angle data]
 %     [p        p-th moment to be computed, default is p=1]
 %     [cent     if true, central moments are computed, default = false]
-%     [dim      compute along this dimension, default is 1]
+%     [dim      compute along this dimension, default is 1st non-singular dimension]
 %
 %     If dim argument is specified, all other optional arguments can be
 %     left empty: circ_moment(alpha, [], [], [], dim)
@@ -29,7 +29,10 @@
 % berens@tuebingen.mpg.de
 
 if nargin < 5
-  dim = 1;
+  dim = find(size(alpha) > 1, 1, 'first');
+  if isempty(dim)
+    dim = 1;
+  end
 end
 
 if nargin < 4
@@ -65,5 +68,4 @@
 
 rho_p = abs(mp);
 mu_p = angle(mp);
-
-
+end
diff --git a/circ_mtest.m b/circ_mtest.m
index ef8ae0f..36e7397 100644
--- a/circ_mtest.m
+++ b/circ_mtest.m
@@ -1,4 +1,4 @@
-function [h mu ul ll] = circ_mtest(alpha, dir, xi, w, d)
+function [h, mu, ul, ll] = circ_mtest(alpha, dir, xi, w, d)
 %
 % [h mu ul ll] = circ_mtest(alpha, dir, xi, w, d)
 %   One-Sample test for the mean angle.
@@ -67,3 +67,4 @@
 
 % compute test via confidence limits (example 27.3)
 h = abs(circ_dist2(dir,mu)) > t;
+end
diff --git a/circ_otest.m b/circ_otest.m
index 5c3125f..5c0a0d0 100644
--- a/circ_otest.m
+++ b/circ_otest.m
@@ -1,4 +1,4 @@
-function [pval m] = circ_otest(alpha, sz, w)
+function [pval, m] = circ_otest(alpha, sz, w)
 %
 % [pval, m] = circ_otest(alpha,sz,w)
 %   Computes Omnibus or Hodges-Ajne test for non-uniformity of circular data.
@@ -66,4 +66,5 @@
 else
   % exact formula by Hodges (1955)
   pval = 2^(1-n) * (n-2*m) * nchoosek(n,m);  
-end
\ No newline at end of file
+end
+end
diff --git a/circ_plot.m b/circ_plot.m
index 3a7e418..ca0a7ee 100644
--- a/circ_plot.m
+++ b/circ_plot.m
@@ -89,7 +89,7 @@
     set(gca,'box','off')
     set(gca,'xtick',[])
     set(gca,'ytick',[])
-    text(1.2, 0, '0'); text(-.05, 1.2, '\pi/2');  text(-1.35, 0, '�\pi');  text(-.075, -1.2, '-\pi/2');
+    text(1.2, 0, '0'); text(-.05, 1.2, '\pi/2');  text(-1.35, 0, '�\pi');  text(-.075, -1.2, '-\pi/2');
 
     
   case 'hist'
@@ -141,3 +141,4 @@
 end
 
 a = gca;
+end
diff --git a/circ_r.m b/circ_r.m
index 4a9d447..77fd947 100644
--- a/circ_r.m
+++ b/circ_r.m
@@ -8,7 +8,7 @@
 %     [d    spacing of bin centers for binned data, if supplied 
 %           correction factor is used to correct for bias in 
 %           estimation of r, in radians (!)]
-%     [dim  compute along this dimension, default is 1]
+%     [dim  compute along this dimension, default: 1st non-singular dimension]
 %
 %     If dim argument is specified, all other optional arguments can be
 %     left empty: circ_r(alpha, [], [], dim)
@@ -29,7 +29,10 @@
 % berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
 
 if nargin < 4
-  dim = 1;
+  dim = find(size(alpha) > 1, 1, 'first');
+  if isempty(dim)
+    dim = 1;
+  end
 end
 
 if nargin < 2 || isempty(w) 
@@ -59,4 +62,4 @@
   c = d/2/sin(d/2);
   r = c*r;
 end
-
+end
diff --git a/circ_rad2ang.m b/circ_rad2ang.m
index 76f8750..7f84e3a 100644
--- a/circ_rad2ang.m
+++ b/circ_rad2ang.m
@@ -8,4 +8,5 @@
 % By Philipp Berens, 2009
 % berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
 
-alpha = alpha / pi *180;
\ No newline at end of file
+alpha = alpha / pi *180;
+end
diff --git a/circ_raotest.m b/circ_raotest.m
index 5b60f9f..3303106 100644
--- a/circ_raotest.m
+++ b/circ_raotest.m
@@ -1,4 +1,4 @@
-function [p U UC] = circ_raotest(alpha)
+function [p, U, UC] = circ_raotest(alpha)
 %
 % [p U UC] = circ_raotest(alpha)
 %   Calculates Rao's spacing test by comparing distances between points on
@@ -59,11 +59,11 @@
 U = (1/2)*U;
 
 % get critical value from table
-[p UC] = getVal(n,U);
+[p, UC] = getVal(n,U);
 
 
 
-function [p UC] = getVal(N, U)
+function [p, UC] = getVal(N, U)
 
 % Table II from Russel and Levitin, 1995
 
@@ -122,9 +122,5 @@
   UC = table(ridx,end-1);
   p = .5;
 end
-
-
-
-
-
-
+end
+end
diff --git a/circ_rtest.m b/circ_rtest.m
index 3eba73a..a96b41f 100644
--- a/circ_rtest.m
+++ b/circ_rtest.m
@@ -1,4 +1,4 @@
-function [pval z] = circ_rtest(alpha, w, d)
+function [pval, z] = circ_rtest(alpha, w, d)
 %
 % [pval, z] = circ_rtest(alpha,w,d)
 %   Computes Rayleigh test for non-uniformity of circular data.
@@ -65,11 +65,4 @@
 %    (24*z - 132*z^2 + 76*z^3 - 9*z^4) / (288*n^2));
 % end
 
-
-
-
-
-
-
-
-
+end
diff --git a/circ_samplecdf.m b/circ_samplecdf.m
index d14c2ef..d3f9394 100644
--- a/circ_samplecdf.m
+++ b/circ_samplecdf.m
@@ -60,10 +60,11 @@
 cdfplottable = [];
 phisplottable = [];
 
-for j=1:length(phis);
+for j=1:length(phis)
    phisplottable = [phisplottable phis(j) phis2(j)]; %#ok<AGROW>
    cdfplottable = [cdfplottable cdf2(j) cdf(j)]; %#ok<AGROW>
 end
 
 phiplot = [phisplottable 2*pi];
-cdfplot = [cdfplottable 1];
\ No newline at end of file
+cdfplot = [cdfplottable 1];
+end
diff --git a/circ_skewness.m b/circ_skewness.m
index 1a8d661..2684e7f 100644
--- a/circ_skewness.m
+++ b/circ_skewness.m
@@ -1,4 +1,4 @@
-function [b b0] = circ_skewness(alpha, w, dim)
+function [b, b0] = circ_skewness(alpha, w, dim)
 %
 % [b b0] = circ_skewness(alpha,w,dim)
 %   Calculates a measure of angular skewness.
@@ -6,7 +6,7 @@
 %   Input:
 %     alpha     sample of angles
 %     [w        weightings in case of binned angle data]
-%     [dim      statistic computed along this dimension, 1]
+%     [dim      statistic computed along this dimension, default: 1st non-singular dimension]
 %
 %     If dim argument is specified, all other optional arguments can be
 %     left empty: circ_skewness(alpha, [], dim)
@@ -25,7 +25,10 @@
 % berens@tuebingen.mpg.de
 
 if nargin < 3
-  dim = 1;
+  dim = find(size(alpha) > 1, 1, 'first');
+  if isempty(dim)
+    dim = 1;
+  end  
 end
 
 if nargin < 2 || isempty(w)
@@ -47,4 +50,5 @@
 % compute skewness 
 theta2 = repmat(theta, size(alpha)./size(theta));
 b = sum(w.*(sin(2*(circ_dist(alpha,theta2)))),dim)./sum(w,dim);
-b0 = rho2.*sin(circ_dist(mu2,2*theta))./(1-R).^(3/2);    % (formula 2.29)
\ No newline at end of file
+b0 = rho2.*sin(circ_dist(mu2,2*theta))./(1-R).^(3/2);    % (formula 2.29)
+end
diff --git a/circ_stats.m b/circ_stats.m
index 2a336e0..44a2250 100644
--- a/circ_stats.m
+++ b/circ_stats.m
@@ -52,11 +52,12 @@
 stats.var = circ_var(alpha,w,d);
 
 % standard deviation
-[stats.std stats.std0] = circ_std(alpha,w,d);
+[stats.std, stats.std0] = circ_std(alpha,w,d);
 
 
 % skewness
-[stats.skewness stats.skewness0] = circ_skewness(alpha,w);
+[stats.skewness, stats.skewness0] = circ_skewness(alpha,w);
 
 % kurtosis
-[stats.kurtosis stats.kurtosis0] = circ_kurtosis(alpha,w);
\ No newline at end of file
+[stats.kurtosis, stats.kurtosis0] = circ_kurtosis(alpha,w);
+end
diff --git a/circ_std.m b/circ_std.m
index c79821d..86225bd 100644
--- a/circ_std.m
+++ b/circ_std.m
@@ -1,4 +1,4 @@
-function [s s0] = circ_std(alpha, w, d, dim)
+function [s, s0] = circ_std(alpha, w, d, dim)
 % [s, s0] = circ_std(alpha, w, d, dim)
 %   Computes circular standard deviation for circular data 
 %   (equ. 26.20, Zar).   
@@ -9,7 +9,7 @@
 %     [d    spacing of bin centers for binned data, if supplied 
 %           correction factor is used to correct for bias in 
 %           estimation of r]
-%     [dim  compute along this dimension, default is 1]
+%     [dim  compute along this dimension, default: 1st non-singular dimension]
 %
 %     If dim argument is specified, all other optional arguments can be
 %     left empty: circ_std(alpha, [], [], dim)
@@ -29,7 +29,10 @@
 % berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
 
 if nargin < 4
-  dim = 1;
+  dim = find(size(alpha) > 1, 1, 'first');
+  if isempty(dim)
+    dim = 1;
+  end  
 end
 
 if nargin < 3 || isempty(d)
@@ -51,4 +54,5 @@
 r = circ_r(alpha,w,d,dim);
 
 s = sqrt(2*(1-r));      % 26.20
-s0 = sqrt(-2*log(r));    % 26.21
\ No newline at end of file
+s0 = sqrt(-2*log(r));   % 26.21
+end
diff --git a/circ_symtest.m b/circ_symtest.m
index 01fa68e..d154021 100644
--- a/circ_symtest.m
+++ b/circ_symtest.m
@@ -33,4 +33,5 @@
 d = circ_dist(alpha,md);
 
 % compute wilcoxon sign rank test
-pval = signrank(d);
\ No newline at end of file
+pval = signrank(d);
+end
diff --git a/circ_var.m b/circ_var.m
index 8c469da..6a0f33f 100644
--- a/circ_var.m
+++ b/circ_var.m
@@ -1,4 +1,4 @@
-function [S s] = circ_var(alpha, w, d, dim)
+function [S, s] = circ_var(alpha, w, d, dim)
 % [S, s] = circ_var(alpha, w, d, dim)
 %   Computes circular variance for circular data 
 %   (equ. 26.17/18, Zar).   
@@ -9,7 +9,7 @@
 %     [d    spacing of bin centers for binned data, if supplied 
 %           correction factor is used to correct for bias in 
 %           estimation of r]
-%     [dim  compute along this dimension, default is 1]
+%     [dim  compute along this dimension, default: 1st non-singular dimension]
 %
 %     If dim argument is specified, all other optional arguments can be
 %     left empty: circ_var(alpha, [], [], dim)
@@ -31,7 +31,10 @@
 % berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
 
 if nargin < 4
-  dim = 1;
+  dim = find(size(alpha) > 1, 1, 'first');
+  if isempty(dim)
+    dim = 1;
+  end    
 end
 
 if nargin < 3 || isempty(d)
@@ -54,4 +57,5 @@
 
 % apply transformation to var
 S = 1 - r;
-s = 2 * S;
\ No newline at end of file
+s = 2 * S;
+end
diff --git a/circ_vmpar.m b/circ_vmpar.m
index 2a81c4b..f247db4 100644
--- a/circ_vmpar.m
+++ b/circ_vmpar.m
@@ -1,4 +1,4 @@
-function [thetahat kappa] = circ_vmpar(alpha,w,d)
+function [thetahat, kappa] = circ_vmpar(alpha,w,d)
 %
 % [thetahat, kappa] = circ_vmpar(alpha, w, d)
 %   Estimate the parameters of a von Mises distribution.
@@ -36,3 +36,4 @@
 kappa = circ_kappa(r);
 
 thetahat = circ_mean(alpha,w);
+end
diff --git a/circ_vmpdf.m b/circ_vmpdf.m
index 524d6a1..7d158c3 100644
--- a/circ_vmpdf.m
+++ b/circ_vmpdf.m
@@ -1,4 +1,4 @@
-function [p alpha] = circ_vmpdf(alpha, thetahat, kappa)
+function [p, alpha] = circ_vmpdf(alpha, thetahat, kappa)
 %
 % [p alpha] = circ_vmpdf(alpha, thetahat, kappa)
 %   Computes the circular von Mises pdf with preferred direction thetahat 
@@ -44,3 +44,4 @@
 % evaluate pdf
 C = 1/(2*pi*besseli(0,kappa));
 p = C * exp(kappa*cos(alpha-thetahat));
+end
diff --git a/circ_vmrnd.m b/circ_vmrnd.m
index 913241e..a948657 100644
--- a/circ_vmrnd.m
+++ b/circ_vmrnd.m
@@ -79,9 +79,4 @@
 if exist('m','var')
   alpha = reshape(alpha,m(1),m(2));
 end
-
-
-
-
-
-
+end
diff --git a/circ_vtest.m b/circ_vtest.m
index 5abbb64..1410ccd 100644
--- a/circ_vtest.m
+++ b/circ_vtest.m
@@ -1,4 +1,4 @@
-function [pval v] = circ_vtest(alpha, dir, w, d)
+function [pval, v] = circ_vtest(alpha, dir, w, d)
 %
 % [pval, v] = circ_vtest(alpha, dir, w, d)
 %   Computes V test for non-uniformity of circular data with a specified 
@@ -75,3 +75,4 @@
 
 % compute p-value from one tailed normal approximation
 pval = 1 - normcdf(u);
+end
diff --git a/circ_wwtest.m b/circ_wwtest.m
index ba833d0..4d7793b 100644
--- a/circ_wwtest.m
+++ b/circ_wwtest.m
@@ -1,4 +1,4 @@
-function [pval table] = circ_wwtest(varargin)
+function [pval, table] = circ_wwtest(varargin)
 %
 % [pval, table] = circ_wwtest(alpha, idx, [w])
 % [pval, table] = circ_wwtest(alpha1, alpha2, [w1, w2])
@@ -159,4 +159,4 @@ function checkAssumption(rw,n)
   else
     error('Invalid use of circ_wwtest. Type help circ_wwtest.')
   end
-end
\ No newline at end of file
+end