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/README.md b/README.md
index dfcfbf9..d1008fe 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,86 @@
-circstat-matlab
-===============
+CircStat for Matlab
+=======================
 
-Matlab Circular Statistics Toolbox
+Toolbox for circular statistics with Matlab. 
+
+## Authors:
+Philipp Berens & Marc J. Velasco
+
+*Email*: berens@tuebingen.mpg.de
+
+*Homepage*: http://www.kyb.tuebingen.mpg.de/~berens/circStat.html
+
+## Contributors:
+Tal Krasovsky
+
+## Reference:
+P. Berens, CircStat: A Matlab Toolbox for Circular Statistics, Journal of Statistical Software, Volume 31, Issue 10, 2009
+http://www.jstatsoft.org/v31/i10
+
+Please cite this paper when the provided code is used. See licensing terms for details.
+
+## Contents:
+- `circ_r` 				Resultant vector length
+- `circ_mean` 			Mean direction of a sample of circular data
+- `circ_axial`			Mean direction for axial data
+- `circ_median`			Median direction of a sample of circular data
+- `circ_std` 			Dispersion around the mean direction (std, mardia)
+- `circ_var` 			Circular variance
+- `circ_skewness`		Circular skewness
+- `circ_kurtosis`		Circular kurtosis
+- `circ_moment`			Circular p-th moment
+- `circ_dist`			Distances around a circle
+- `circ_dist2`			Pairwise distances around a circle
+- `circ_confmean` 		Confidence intervals for mean direction
+- `circ_stats`			Summary statistics
+
+ 
+- `circ_rtest`			Rayleigh's test for nonuniformity
+- `circ_otest`			Hodges-Ajne test (omnibus test) for nonuniformity
+- `circ_raotest`		Rao's spacing test for nonuniformity
+- `circ_vtest`			V-Test for nonuniformity with known mean direction
+- `circ_medtest`		Test for median angle
+- `circ_mtest`			One-sample test for specified mean direction
+- `circ_wwtest`			Multi-sample test for equal means, one-factor ANOVA
+- `circ_hktest` 		Two-factor ANOVA
+- `circ_ktest`      	Test for equal concentration parameter
+- `circ_symtest`		Test for symmetry around median angle
+- `circ_kuipertest`	 	Test whether two distributions are identical (like KS test)
+
+ 
+- `circ_corrcc`			Circular-circular correlation coefficient
+- `circ_corrcl`			Circular-linear correlation coefficient
+
+ 
+- `circ_kappa` 			Compute concentration parameter of a von Mises distribution
+- `circ_plot`			Visualization for circular data
+- `circ_clust`    		Simple clustering for circular data
+- `circ_samplecdf`	 	Evaluate CDF of a sample of angles
+
+ 
+- `rad2ang`				Convert radian to angular values
+- `ang2rad`				Convert angular to radian values
+
+All functions take arguments in radians (expect for `ang2rad`). For a detailed description of arguments and outputs consult the help text in the files.
+
+Since 2010, most functions for descriptive statistics can be used in Matlab style matrix computations. As a last argument, add the dimension along which you want to average. This changes the behavior slightly from previous relaeses, in that input is not reshaped anymore into vector format. Per default, all computations are performed columnwise (along dimension 1). If you prefer to use the old functions, for now they are contained in the subdirectory 'old'.
+
+## References:
+- E. Batschelet, Circular Statistics in Biology, Academic Press, 1981
+- N.I. Fisher, Statistical analysis of circular data, Cambridge University Press, 1996
+- S.R. Jammalamadaka et al., Topics in circular statistics, World Scientific, 2001
+- J.H. Zar, Biostatistical Analysis, Prentice Hall, 1999
+
+
+The implementation follows in most cases 'Biostatistical Analysis' and all referenced equations and tables are taken from this book, if not otherwise noted. In some cases, the other books were preferred for implementation was more straightforward for solutions presented there.
+
+If you have suggestions, bugs or feature requests or want to contribute code, please email us.
+
+## Disclaimer:
+All functions in this toolbox were implemented with care and tested on the examples presented in 'Biostatistical Analysis' were possible. Nevertheless, they may contain errors or bugs, which may affect the outcome of your analysis. We do not take responsibility for any harm coming from using this toolbox, neither if it is caused by errors in the software nor if it is caused by its improper application. Please email us any bugs you find.
+
+By Philipp Berens and Marc J. Velasco, 2009
+
+berens@tuebingen.mpg.de , velasco@ccs.fau.edu - www.kyb.mpg.de/~berens/circStat.html
+
+Distributed under Open Source BSD License
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 c0bb94c..1cccc51 100644
--- a/circ_axial.m
+++ b/circ_axial.m
@@ -20,10 +20,11 @@
 % By Philipp Berens, 2009
 % berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
 
-% TESTING GIT
+
 
 if nargin < 2
     p = 1;
 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 c2d4cda..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)
@@ -55,7 +55,8 @@
 end
 
 if any(n<10)
-  warning('Test not applicable. Sample size in at least one group to small.') %#ok<WNTAG>
+  warning('CIRCSTAT:circ_cmtest:sampleSizeTooSmall', ...
+      'Test not applicable. Sample size in at least one group too small.') %#ok<WNTAG>
 end
   
 
@@ -67,7 +68,7 @@
 if pval < 0.05
   med = NaN;
 end
-
+end
 
 
 
@@ -87,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 44a9b9b..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)
@@ -67,13 +70,11 @@
     t(i) = sqrt(n(i)^2-(n(i)^2-R(i)^2)*exp(c2/n(i)));      % equ. 26.25
   else 
     t(i) = NaN;
-    warning('Requirements for confidence levels not met.');
+    warning('CIRCSTAT:circ_confmean:requirementsNotMet', ...
+        'Requirements for confidence levels not met.');
   end
 end
 
 % apply final transform
 t = acos(t./R);
-  
-
-
-
+end
diff --git a/circ_corrcc.m b/circ_corrcc.m
index 2071e65..5fc19bc 100644
--- a/circ_corrcc.m
+++ b/circ_corrcc.m
@@ -1,6 +1,6 @@
-function [rho pval] = circ_corrcc(alpha1, alpha2)
+function [rho, pval] = circ_corrcc(alpha1, alpha2)
 %
-% [rho pval ts] = circ_corrcc(alpha1, alpha2)
+% [rho pval] = circ_corrcc(alpha1, alpha2)
 %   Circular correlation coefficient for two circular random variables.
 %
 %   Input:
@@ -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 e4dd7fd..8bbdcab 100644
--- a/circ_corrcl.m
+++ b/circ_corrcl.m
@@ -1,6 +1,6 @@
-function [rho pval] = circ_corrcl(alpha, x)
+function [rho, pval] = circ_corrcl(alpha, x)
 %
-% [rho pval ts] = circ_corrcc(alpha, x)
+% [rho pval] = circ_corrcc(alpha, x)
 %   Correlation coefficient between one circular and one linear random
 %   variable.
 %
@@ -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 d541110..f260475 100644
--- a/circ_dist2.m
+++ b/circ_dist2.m
@@ -1,6 +1,6 @@
 function r =  circ_dist2(x,y)
 %
-% r = circ_dist(alpha, beta)
+% r = circ_dist2(alpha, beta)
 %   All pairwise difference x_i-y_j around the circle computed efficiently.
 %
 %   Input:
@@ -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 25ede86..5568f00 100644
--- a/circ_ktest.m
+++ b/circ_ktest.m
@@ -39,7 +39,8 @@
 rbar = (R1+R2)/(n1+n2);
 
 if rbar < .7
-    warning('resultant vector length should be > 0.7') %#ok<WNTAG>
+    warning('CIRCSTAT:circ_ktest:vectorTooShort', ...
+        'Resultant vector length should be > 0.7') %#ok<WNTAG>
 end
 
 % calculate test statistic
@@ -50,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 37e6b1f..4669997 100644
--- a/circ_kuipertest.m
+++ b/circ_kuipertest.m
@@ -1,6 +1,6 @@
 function [pval, k, K] = circ_kuipertest(alpha1, alpha2, res, vis_on)
 
-% [pval, k, K] = circ_kuipertest(sample1, sample2, res, vis_on)
+% [pval, k, K] = circ_kuipertest(alpha1, alpha2, res, vis_on)
 %
 %   The Kuiper two-sample test tests whether the two samples differ 
 %   significantly.The difference can be in any property, such as mean 
@@ -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,21 +83,22 @@
 
 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); 
    if row == 0
        error('N too small.');
    else
-      warning('N=%d not found in table, using closest N=%d present.',n,nn(row)) %#ok<WNTAG>
+      warning('CIRCSTAT:circ_kuipertest:nNotFound', ...
+          'N=%d not found in table, using closest N=%d present.',n,nn(row)) %#ok<WNTAG>
    end
 end
 
@@ -109,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 c9108e0..c97916f 100644
--- a/circ_mean.m
+++ b/circ_mean.m
@@ -1,12 +1,12 @@
-function [mu ul ll] = circ_mean(alpha, w, dim)
+function [mu, ul, ll] = circ_mean(alpha, w, dim)
 %
-% mu = circ_mean(alpha, w)
+% [mu ul ll] = circ_mean(alpha, w, dim)
 %   Computes the mean direction for circular data.
 %
 %   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 4288442..aaca36f 100644
--- a/circ_median.m
+++ b/circ_median.m
@@ -1,15 +1,15 @@
 function med = circ_median(alpha,dim)
 %
-% med = circ_median(alpha)
+% med = circ_median(alpha, dim)
 %   Computes the median direction for circular data.
 %
 %   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:
-%     mu		median direction
+%     med		median direction
 %
 %   circ_median can be slow for large datasets
 %
@@ -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,14 +51,15 @@
 
   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);
   end
 
   if m > 1
-    warning('Ties detected.') %#ok<WNTAG>
+    warning('CIRCSTAT:circ_median:tiesDetected', ...
+        'Ties detected.') %#ok<WNTAG>
   end
 
   md = circ_mean(beta(idx));
@@ -69,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 0950ed9..a846886 100644
--- a/circ_medtest.m
+++ b/circ_medtest.m
@@ -1,6 +1,6 @@
 function pval = circ_medtest(alpha,md)
 %
-% [pval, z] = circ_medtest(alpha,w)
+% pval = circ_medtest(alpha, md)
 %   Tests for significance of the median.
 %   H0: the population has median angle md
 %   HA: the population has not median angle md
@@ -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 241750f..beb98c9 100644
--- a/circ_moment.m
+++ b/circ_moment.m
@@ -1,6 +1,6 @@
-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 cbar sbar] = 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 
 %   of the angular data in angle.
 %
@@ -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 d1e82f7..36e7397 100644
--- a/circ_mtest.m
+++ b/circ_mtest.m
@@ -1,6 +1,6 @@
-function [h mu ul ll] = circ_mtest(alpha, dir, xi, w, d)
+function [h, mu, ul, ll] = circ_mtest(alpha, dir, xi, w, d)
 %
-% [pval, z] = circ_mtest(alpha, dir, w, d)
+% [h mu ul ll] = circ_mtest(alpha, dir, xi, w, d)
 %   One-Sample test for the mean angle.
 %   H0: the population has mean dir.
 %   HA: the population has not mean dir.
@@ -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 b5db47e..c072482 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.
@@ -13,7 +13,7 @@
 %     alpha	sample of angles in radians
 %     [sz   step size for evaluating distribution, default 1 degree
 %     [w		number of incidences in case of binned angle data]
-
+%
 %   Output:
 %     pval  p-value 
 %     m     minimum number of samples falling in one half of the circle
@@ -65,17 +65,7 @@
   pval = sqrt(2*pi) / A * exp(-pi^2/8/A^2);
 else
   % exact formula by Hodges (1955)
-  pval = 2^(1-n) * (n-2*m) * nchoosek(n,m);  
+  % pval = 2^(1-n) * (n-2*m) * nchoosek(n,m);  % revised below for numerical stability
+  pval = exp((1-n)*log(2) + log(n-2*m) + gammaln(n+1) - gammaln(m+1) - gammaln(n-m+1));
+end
 end
-
-  
-  
-
-
-
-
-
-
-
-
-
diff --git a/circ_plot.m b/circ_plot.m
index c07a0fe..ca0a7ee 100644
--- a/circ_plot.m
+++ b/circ_plot.m
@@ -1,6 +1,6 @@
 function a = circ_plot(alpha, format, formats, varargin)
 %
-% r = circ_plot(alpha, ...)
+% a = circ_plot(alpha, ...)
 %   Plotting routines for circular data.
 %
 %   Input:
@@ -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 d2768dc..77fd947 100644
--- a/circ_r.m
+++ b/circ_r.m
@@ -1,5 +1,5 @@
 function r = circ_r(alpha, w, d, dim)
-% r = circ_r(alpha, w, d)
+% r = circ_r(alpha, w, d, dim)
 %   Computes mean resultant vector length for circular data.
 %
 %   Input:
@@ -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 9a0ae4d..7f84e3a 100644
--- a/circ_rad2ang.m
+++ b/circ_rad2ang.m
@@ -1,6 +1,6 @@
 function alpha = circ_rad2ang(alpha)
 
-% alpha = circ-rad2ang(alpha)
+% alpha = circ_rad2ang(alpha)
 %   converts values in radians to values in degree
 %
 % Circular Statistics Toolbox for Matlab
@@ -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 ec31576..3303106 100644
--- a/circ_raotest.m
+++ b/circ_raotest.m
@@ -1,5 +1,5 @@
-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
 %   a circle to those expected from a uniform distribution.
@@ -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 eaaafd8..a96b41f 100644
--- a/circ_rtest.m
+++ b/circ_rtest.m
@@ -1,6 +1,6 @@
-function [pval z] = circ_rtest(alpha, w, d)
+function [pval, z] = circ_rtest(alpha, w, d)
 %
-% [pval, z] = circ_rtest(alpha,w)
+% [pval, z] = circ_rtest(alpha,w,d)
 %   Computes Rayleigh test for non-uniformity of circular data.
 %   H0: the population is uniformly distributed around the circle
 %   HA: the populatoin is not distributed uniformly around the circle
@@ -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 39df433..d3f9394 100644
--- a/circ_samplecdf.m
+++ b/circ_samplecdf.m
@@ -1,5 +1,5 @@
 function [phis, cdf, phiplot, cdfplot] = circ_samplecdf(thetas, resolution)
-
+%
 % [phis, cdf, phiplot, cdfplot] = circ_samplecdf(thetas, resolution)
 %
 %   Helper function for circ_kuipertest.
@@ -60,15 +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];
-
-
-
-
-
+end
diff --git a/circ_skewness.m b/circ_skewness.m
index 8db086a..2684e7f 100644
--- a/circ_skewness.m
+++ b/circ_skewness.m
@@ -1,12 +1,12 @@
-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.
 %
 %   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)
@@ -48,5 +51,4 @@
 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)
-
-
+end
diff --git a/circ_stats.m b/circ_stats.m
index 6194907..44a2250 100644
--- a/circ_stats.m
+++ b/circ_stats.m
@@ -1,6 +1,6 @@
 function stats = circ_stats(alpha, w, d)
 %
-% stats = circ_stats(alpha, w)
+% stats = circ_stats(alpha, w, d)
 %   Computes descriptive statistics for circular data.
 %
 %   Input:
@@ -52,15 +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);
-
-
-
-
+[stats.kurtosis, stats.kurtosis0] = circ_kurtosis(alpha,w);
+end
diff --git a/circ_std.m b/circ_std.m
index 115ebd8..86225bd 100644
--- a/circ_std.m
+++ b/circ_std.m
@@ -1,5 +1,5 @@
-function [s s0] = circ_std(alpha, w, d, dim)
-% s = 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,7 +54,5 @@
 r = circ_r(alpha,w,d,dim);
 
 s = sqrt(2*(1-r));      % 26.20
-s0 = sqrt(-2*log(r));    % 26.21
-
-
-
+s0 = sqrt(-2*log(r));   % 26.21
+end
diff --git a/circ_symtest.m b/circ_symtest.m
index 3c81372..d154021 100644
--- a/circ_symtest.m
+++ b/circ_symtest.m
@@ -1,10 +1,10 @@
 function pval = circ_symtest(alpha)
 %
-% [pval, z] = circ_symtest(alpha,w)
+% pval = circ_symtest(alpha)
 %   Tests for symmetry about the median.
 %   H0: the population is symmetrical around the median
 %   HA: the population is not symmetrical around the median
-
+%
 %
 %   Input:
 %     alpha	sample of angles in radians
@@ -34,7 +34,4 @@
 
 % compute wilcoxon sign rank test
 pval = signrank(d);
-
-
-
-
+end
diff --git a/circ_var.m b/circ_var.m
index 21013f7..6a0f33f 100644
--- a/circ_var.m
+++ b/circ_var.m
@@ -1,5 +1,5 @@
-function [S s] = circ_var(alpha, w, d, dim)
-% 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 75c40c8..f247db4 100644
--- a/circ_vmpar.m
+++ b/circ_vmpar.m
@@ -1,6 +1,6 @@
-function [thetahat kappa] = circ_vmpar(alpha,w,d)
-
-% r = 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.
 %
 %   Input:
@@ -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 4636a66..11cc2ff 100644
--- a/circ_vmpdf.m
+++ b/circ_vmpdf.m
@@ -1,6 +1,6 @@
-function [p alpha] = circ_vmpdf(alpha, thetahat, kappa)
-
-% [p alpha] = circ_vmpdf(alpha, w, p)
+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 
 %   and concentration kappa at each of the angles in alpha
 %
@@ -42,5 +42,5 @@
 alpha = alpha(:);
 
 % evaluate pdf
-C = 1/(2*pi*besseli(0,kappa));
-p = C * exp(kappa*cos(alpha-thetahat));
+p = exp( kappa*(cos(alpha-thetahat)-1) ) / (2*pi*besseli(0,kappa,1));
+end
diff --git a/circ_vmrnd.m b/circ_vmrnd.m
index c0e3e55..00285bc 100644
--- a/circ_vmrnd.m
+++ b/circ_vmrnd.m
@@ -1,8 +1,8 @@
 function alpha = circ_vmrnd(theta, kappa, n)
-
+%
 % alpha = circ_vmrnd(theta, kappa, n)
 %   Simulates n random angles from a von Mises distribution, with preferred 
-%   direction thetahat and concentration parameter kappa.
+%   direction theta and concentration parameter kappa.
 %
 %   Input:
 %     [theta    preferred direction, default is 0]
@@ -47,7 +47,7 @@
 
 % if kappa is small, treat as uniform distribution
 if kappa < 1e-6
-    alpha = 2*pi*rand(n,1);
+    alpha = 2*pi*rand(n,1)-pi;
     return
 end
 
@@ -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 a92ff30..4d7793b 100644
--- a/circ_wwtest.m
+++ b/circ_wwtest.m
@@ -1,4 +1,5 @@
-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])
 %   Parametric Watson-Williams multi-sample test for equal means. Can be
@@ -106,13 +107,17 @@
 function checkAssumption(rw,n)
 
   if n >= 11 && rw<.45
-    warning('Test not applicable. Average resultant vector length < 0.45.') %#ok<WNTAG>
+    warning('CIRCSTAT:circ_wwtest:vectorTooShort', ...
+        'Test not applicable. Average resultant vector length < 0.45.') %#ok<WNTAG>
   elseif n<11 && n>=7 && rw<.5
-    warning('Test not applicable. Average number of samples per population 6 < x < 11 and average resultant vector length < 0.5.') %#ok<WNTAG>
+    warning('CIRCSTAT:circ_wwtest:sampleSize6x11AndVectorTooShort', ...
+        'Test not applicable. Average number of samples per population 6 < x < 11 and average resultant vector length < 0.5.') %#ok<WNTAG>
   elseif n>=5 && n<7 && rw<.55
-    warning('Test not applicable. Average number of samples per population 4 < x < 7 and average resultant vector length < 0.55.') %#ok<WNTAG>
+    warning('CIRCSTAT:circ_wwtest:sampleSize4x7AndVectorTooShort', ...
+        'Test not applicable. Average number of samples per population 4 < x < 7 and average resultant vector length < 0.55.') %#ok<WNTAG>
   elseif n < 5
-    warning('Test not applicable. Average number of samples per population < 5.') %#ok<WNTAG>
+    warning('CIRCSTAT:circ_wwtest:sampleSizeTooSmall', ...
+        'Test not applicable. Average number of samples per population < 5.') %#ok<WNTAG>
   end
 
 end
@@ -155,4 +160,3 @@ function checkAssumption(rw,n)
     error('Invalid use of circ_wwtest. Type help circ_wwtest.')
   end
 end
-