diff git a/raincloud_plots/rm_raincloud.m b/raincloud_plots/rm_raincloud.m
index dd1fe83c59d6c04c830367675167908f2e0404cd..5e0caf8a38da45cfae105a4dd7845b5275444d8b 100644
 a/raincloud_plots/rm_raincloud.m
+++ b/raincloud_plots/rm_raincloud.m
@@ 1,17 +1,17 @@
%% Raincloud lineplot for repeatedmeasures data
% it might be that 'patch' is a much better plotting function to plot the
% densities than 'area'; it seems to handle overplotting of different
% areas with different zero lines very well. This would avoid the horrible
% gymnastics that the repeated measures raincloudplot currently has to deal
% with.

% This is a workinprogress example inspired by this post on Matlab
% answers https://uk.mathworks.com/matlabcentral/answers/1387areaplotwithgradient

% Use like: h = rm_raincloud(data, colours)
+% Use like: h = rm_raincloud(data, colours, plot_top_to_bottom, density_type, bandwidth)
% Where 'data' is an M x N cell array of N data series and M measurements
% And 'colours' is an N x 3 array defining the colour to plot each data series
% h is a cell array containing handles of the various figure parts
+% plot_top_to_bottom: Default plots lefttoright, set to 1 to rotate.
+% density_type: 'ks' (default) or 'RASH'. 'ks' uses matlab's inbuilt 'ksdensity' to
+% determine the shape of the rainclouds. 'RASH' will use the 'rst_RASH'
+% method from Cyril Pernet's Robust Stats toolbox, if that function is on
+% your matlab path.
+% bandwidth: If density_type == 'ks', determines bandwidth of density estimate
+% h is a cell array containing handles of the various figure parts:
+% h.p{i,j} is the handle to the density plot from data{i,j}
+% h.s{i,j} is the handle to the 'raindrops' (individual datapoints) from data{i,j}
+% h.m(i,j) is the handle to the single, large dot that represents mean(data{i,j})
+% h.l(i,j) is the handle for the line connecting h.m(i,j) and h.m(i+1,j)
%% TODO:
% Patch can create colour gradients using the 'interp' option to 'FaceColor'. Allow this?
@@ 51,6 +51,9 @@ for i = 1:nper
case 'ks'
[ks{i,j}, x{i,j}] = ksdensity(data{i,j}, 'NumPoints', nbins(i,j), 'bandwidth', bandwidth);
case 'rash'
+ % check for rst_RASH function (from Robust stats toolbox) in path, fail if not found
+ assert(exist('rst_RASH','file') == 2, 'Could not compute density using RASH method. Do you have the Robust Stats toolbox on your path?');
+
[x{i,j}, ks{i,j}] = rst_RASH(data{i,j});
% override default 'nbins' as rst_RASH determines number of bins
nbins(i,j) = size(ks{i,j},2);