Huynh-Feldt epsilon general procedure

Calc. the Huynh-Feldt epsilon for any arbitrary combination of between and within subject factors

%function [EpsHF EpsList EpsGG]=GenCalcHFEps(Y,BTFacs,WInFacs,S)
%
% This will calculate the Geisser-Greenhouse and Huynh-Feldt epsilon value
% for the general case when given a univariate dataset with any amount of
% between or within subject factors.
%
% After calling this function, the intention is to multiply the degrees of
% freedom for tests in a repeated measures ANOVA by the corresponding
% epsilon value given by this code to obtain a corrected p-value. The
% F-statistic does not change as a result of this.
%
% Inputs:
% Y- A column vector of the dependent variable (the value
% measured) for each data point.
% BTFacs- A matrix of between subject factors. Must have the same
% number of rows as Y. Each factor is down a different column
% of BTFacs, and the value in each row denotes the level of
% that factor for each corresponding datapoint in Y. If there
% are no between subject factors, input an empty matrix, [].
% WInFacs- A matrix of within subject factors. Must have the same
% number of rows as Y. Each factor is down a different column
% of WInFacs, and the value in each row denotes the level of
% that factor for each corresponding datapoint in Y.
% S- A column vector of subject numbers corresponding to each
% datapoint in Y. Must have the same number of rows as Y. If
% left empty, the program assumes all subjects are entered in
% the same order for all combinations of factors.
%
% Outputs:
% EpsHF- A row vector of all the possible values for the Huynh-Feldt
% epsilon, correspondding to all main effects for ecah within
% subject factor, and all interactions of within subject
% factors.
% EpsList- A cell array of a text list of the effects corresponding to
% the positions in EpsHF. For example 'A' means a main effect
% for the first (within subject) factor, 'AB' means a 2way
% interaction for the first two (within subject) factors,
% etc.
% EpsGG- The Geisser-Greenhouse epsilon values, in teh same form as
% EpsHF.
%
% The Huynh-Feldt epsilon value is less conservative while maintaining the
% proper Type I error rate, which is why the program title focuses on that.
% The Geisser-Greenhouse epsilon is calculated first as part of the
% process of calculating the Huynh-Feldt value, so I've allowed it to be
% optionally returned as well.
%
% This code follows a procedure described in Huynh(1978), but includes a
% modification of the last step, as described in Chen & Dunlap(1994).
%
% I've tested this with results from SPSS for a couple of sample datasets
% I've found in various websites, and it gives results that match
% According to what I've seen and read though, SPSS and SAS give results
% for EpsHF that are wrong when both between and with-in factors are
% present, and thus they slightly differ from the results this program,
% although the results for EpsGG are identical. To reproduce the SPSS
% results, uncomment line 153. See other comments near there for references
% on why I think SPSS/SAS is wrong.
%
%
% References:
%
% Huynh H. "Some approximate tests for repeated measurement designs",
% Psychometrika (1978)
%
% Chen, RS and Dunlap, WP "A MonteCarlo STudy on the Performance of a
% Corrected Formula for eps(tilda) suggested by Lecoutre", Journal of
% Educational Statistics (1994)