function [efs,F,cdfs,p,eps,dfs,b,y2,sig,eta2]=repanova(d,D,fn,gg,alpha);	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% N-way repeated measures ANOVAs 
% (all factors must be repeated measures)
%
% Input:
%
% d = data	A matrix with rows = replications (eg subjects) and
%		columns = conditions 
% D = factors	A vector with as many entries as factors, each entry being
%		the number of levels for that factor
%
%		Data matrix d must have as many columns (conditions) as
%		the product of the elements of the factor matrix D
%
%		First factor rotates slowest; last factor fastest
% 
% 	Eg, in a D=[2 3] design: factor A with 2 levels; factor B with 3:
%	    data matrix d must be organised:
%
%		A1B1	A1B2	A1B3	A2B1	A2B2	A2B3
% 	rep1
%	rep2
%	...
%
% fn = factor names (string array) [optional]
% gg = Greenhouse-Geiser correction flag [default=1]
% alpha = nominal alpha value [default=.05]
%	
% Output:
%
% efs 	= effect, eg [1 2] = interaction between factor 1 and factor 2
% F   	= F value
% cdfs 	= corrected df's (using Greenhouse-Geisser)
% p     = p-value
% eps   = epsilon
% dfs   = original dfs
% b     = betas
% eta2  = partial eta-squared (measure of effect size)
%
% Rik.Henson@mrc-cbu.cam.ac.uk, 17/3/03
%
% Added eta^2 R.Henson, 17/8/08
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<5
	alpha=0.05;
end

if nargin<4
	gg=1;		% G-G correction
end

if nargin<3		% No naming of factors provided (added 26/3/03)
   for f=1:length(D)
	fn{f}=sprintf('%d',f);
   end
end

Nf = length(D);		% Number of factors
Nd = prod(D);		% Number of conditions
Ne = 2^Nf - 1;		% Number of effects
Nr = size(d,1);		% Number of replications (eg subjects)

sig=cell(2,1);

if size(d,2) ~= Nd
	error(sprintf('data has %d conditions; design only %d',size(d,2),Nd))
end

sc = cell(Nf,2);	% create main effect/interaction component contrasts
for f = 1 : Nf
	sc{f,1} = ones(D(f),1);
	sc{f,2} = detrend(eye(D(f)),0);
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Added
sy = cell(Nf,2);	% create main effect/interaction component contrasts
for f = 1 : Nf
	sy{f,1} = ones(D(f),1)/D(f);
	sy{f,2} = eye(D(f));
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tss = sum(detrend(d(:),0).^2);    % total sum of squares (only needed for eta2 below)

for e = 1 : Ne

	cw = num2binvec(e,Nf)+1;

	c  = sc{1,cw(Nf)};	% create full contrasts
	for f = 2 : Nf
		c = kron(c,sc{f,cw(Nf-f+1)});
	end

	y = d * c;		% project data to contrast sub-space

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Added
	cy  = sy{1,cw(Nf)};	% create full contrasts
	for f = 2 : Nf
		cy = kron(cy,sy{f,cw(Nf-f+1)});
	end
	y2{e} = mean(d * cy);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	nc = size(y,2);
	df1 = rank(c);
	df2 = df1*(Nr-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% long, SPM way (GG doesn't work):
%
%	y = y(:);
%	X = kron(eye(nc),ones(Nr,1));
%	b{e} = pinv(X)*y;
%	Y = X*b{e};
%	R = eye(nc*Nr)- X*pinv(X);
%	r = y - Y;
%%	V = r*r';
%	V = y*y';
%	eps(e) = trace(R*V)^2 / (df1 * trace((R*V)*(R*V)));
%
%%	ss = Y'*y;
%%	mse = (y'*y - ss)/df2;
%%	mss = ss/df1;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% computationally simpler way

	b{e} = mean(y);			% betas

	ss   =  sum(y*b{e}');
	mse  = (sum(diag(y'*y)) - ss)/df2;
	mss  =  ss/df1;

	if gg
		V      = cov(y);			% sample covariance
		eps(e) = trace(V)^2 / (df1*trace(V'*V));% Greenhouse-Geisser 
	else
		eps(e) = 1;
	end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        eta2(e,1) = ss/tss;             % Approx measure of effect size (total eta^2)
        eta2(e,2) = ss/(ss+mse*df2);    % Another approx measure of effect size (partial eta^2)
	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	efs{e} = Nf+1-find(cw==2);			% codes which effect 

	F(e)   = mss/mse;

	dfs(e,:)  = [df1 df2];
	cdfs(e,:) = eps(e)*dfs(e,:);

	p(e) = 1-spm_Fcdf(F(e),cdfs(e,:));

	if p(e) < alpha; sig{1}=[sig{1} e]; end
	if p(e) < alpha/Ne; sig{2}=[sig{2} e]; end

	en=fn{efs{e}(1)};	% Naming of factors modified (added 21/6/06)
	for f = 2:length(efs{e})
		en = [fn{efs{e}(f)} en];
	end

	disp(sprintf('Effect %02d: %-18s F(%3.2f,%3.2f)=%4.3f,\tp=%4.3f (eta2=%4.3f)',...
		e,en,cdfs(e,1),cdfs(e,2),F(e),p(e),eta2(e,1)))

%	disp(sprintf('Effect %02d: %-18s F(%3.2f,%3.2f)=%4.3f,\tp=%4.3f',...
%		e,en,dfs(e,1),dfs(e,2),F(e),1-spm_Fcdf(F(e),dfs(e,:))))
end

disp(sprintf('\n'))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Added
%for e = 1 : Ne
%	en=fn{efs{e}(1)};	
%	for f = 2:length(efs{e})
%		en = [en fn{efs{e}(f)}];
%	end
%	disp(sprintf('Effect %02d: %-18s \n',e,en))
%	disp(y2{e})
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sub-function to code all main effects/interactions

function b = num2binvec(d,p)

if nargin<2
	p = 0;		% p is left-padding with zeros option
end

d=abs(round(d));

if(d==0)
	b = 0;
else
	b=[];
 	while d>0
		b=[rem(d,2) b];
		d=floor(d/2);
 	end
end

b=[zeros(1,p-length(b)) b];

