function [results]=CalcRmAnova(anovadata,factor,betweenfactor,n,pointwisestatus,...
                                    CellSum,CellGroupSum,factorstructure,betweenfactorstructure,...
                                    betweensubjectcell,nroffactors,nrofbetweenfactors,nrofsubjects,...
                                    SubjectDf,TotQSf,SubjectQSf)
%   '[results] = calc_rm_anova_fast( anovadata, factor, betweenfactor, n , pointwisestatus,...
%   CellSum, CellGroupSum, factorstructure, betweenfactorstructure,...
%   betweensubjectcell, nroffactors, nrofbetweenfactors, nrofsubjects, SubjectQS,...
%   SubjectDf, TotQSf, SubjectQSf)'
%   
%   calculates n-way-interaction and main-effects (if n==1) for the
%   data given in anovadata.
%
%   The number of rows in anovadata corresponds the number of 
%   subjects. 'factor' is a cell containing the name of the factors
%   in the first column and the number of gradations in the second
%   column. 'n' specifies the order of the interactions to be 
%   calculated. Factors must be in hierarchical order, starting 
%   from the highest level. 
%
%   betweenfactor can contain one grouping factor.
%   
%   'results' is a cell array containing the following information in the
%   columns with the denoted indices:
%
%   1:Name 2:ConditionSumSquare 3:ConditionDf 4:ConditionMeanSquare 5:ErrorSumSquare 
%   6:ErrorDf 7:ErrorMeanSquare 8:F-Value 9:p-Value
%
%   'results{3,2}' for instance contains the ConditionQS for the third
%   tested effect.
%   
%   Lower order effects are subtracted afterwards in RMANOVA.
%
%   See also RMANOVA, ANOVATOAVGFILES.


%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION.                                                                      

%print sums

if nargin<5;pointwisestatus = 0;end
PrintSumStatus = 0;

if PrintSumStatus
    fprintf(1,['Calculating ',num2str(n),'-order-effects:',char(10)]);
end

if isempty(nrofbetweenfactors)
    nrofbetweenfactors = 0;
end

totalnroffactors = nroffactors+nrofbetweenfactors;
combinations = nchoosek(linspace(1,totalnroffactors,totalnroffactors),n);
results = cell(size(combinations,1),6);

for combinationzaehler=1:size(combinations,1)
    
    if PrintSumStatus
        fprintf(1,['current effect: ''',num2str(combinations(combinationzaehler,:)),'''',char(10)]);
        fprintf(1,['----------------------------------------------------------------',char(10)]);
    end
    %store current effect label
    %==================================
    results{combinationzaehler,1}=combinations(combinationzaehler,:);     
    
    
    factorstructure_ia = factorstructure;
    betweenfactorstructure_ia = betweenfactorstructure;
    
    %trim factorstructure
    %====================================
    komponenten = cell(n,2);
    nrofinteractioncells = 1;
    nrofwithininteractioncells = 1;
    nrofbetweeninteractioncells = 1;
    sizestring = [];
    conditionindexstring = ''; 
    keeparray = [];
    withincount = 1;
    betweencount = 1;
    for komponentenzaehler=1:n
        if combinations(combinationzaehler,komponentenzaehler)>nroffactors
            keeparray = [keeparray , combinations(combinationzaehler,komponentenzaehler)-nroffactors];
            %Name der Komponente
            komponenten{komponentenzaehler,1}=betweenfactor{combinations(combinationzaehler,komponentenzaehler)-nroffactors,1};
            %nrofgradations der Komponente
            komponenten{komponentenzaehler,2}=betweenfactor{combinations(combinationzaehler,komponentenzaehler)-nroffactors,2};
            if ~isempty(sizestring)
                sizestring = [sizestring,',',num2str(komponenten{komponentenzaehler,2})];
                conditionindexstring = [conditionindexstring,',','betweenfactorstructure_ia(q,',num2str(betweencount),')'];
            else
                sizestring = [num2str(komponenten{komponentenzaehler,2})];
                conditionindexstring = ['betweenfactorstructure_ia(q,',num2str(betweencount),')'];
            end
            nrofinteractioncells = nrofinteractioncells * komponenten{komponentenzaehler,2};
            nrofbetweeninteractioncells = nrofbetweeninteractioncells * komponenten{komponentenzaehler,2};
            betweencount = betweencount+1;
        else
            %Name der Komponente
            komponenten{komponentenzaehler,1}=factor{combinations(combinationzaehler,komponentenzaehler),1};
            %nrofgradations der Komponente
            komponenten{komponentenzaehler,2}=factor{combinations(combinationzaehler,komponentenzaehler),2};
            if ~isempty(sizestring)
                sizestring = [sizestring,',',num2str(komponenten{komponentenzaehler,2})];
                conditionindexstring = [conditionindexstring,',','factorstructure_ia(',num2str(withincount),',r)'];
            else
                sizestring = [num2str(komponenten{komponentenzaehler,2})];
                conditionindexstring = ['factorstructure_ia(',num2str(withincount),',r)'];
            end
            nrofinteractioncells = nrofinteractioncells * komponenten{komponentenzaehler,2};
            nrofwithininteractioncells = nrofwithininteractioncells * komponenten{komponentenzaehler,2};
            withincount = withincount+1;
        end
    end
    
    if n==1
        sizestring = [sizestring,',1'];
    end
   
    evalstring = ['InteractionSum = zeros(',sizestring,');']; %--->InteractionMean = zeros(.....);
    eval(evalstring);
    
    %conditionmat
    %--------------
    conditionmat = zeros(size(InteractionSum));
    conditionmat(1:end)=linspace(1,nrofinteractioncells,nrofinteractioncells);
    
    %conditionvec
    %---------------
    conditionmat = zeros(size(InteractionSum));
    conditionmat(1:end)=linspace(1,nrofinteractioncells,nrofinteractioncells);
    if length(conditionmat)>1
        betweenfactorstructure_ia = betweenfactorstructure_ia(:,keeparray);
        if ~isempty(betweenfactorstructure_ia)
            conditionvec = zeros(size(betweenfactorstructure_ia,1),size(anovadata,2));
        else
            conditionvec = zeros(1,size(anovadata,2));
        end
        %restzeilen löschen
        deletevec=[];
        deletezaehler=1;
        for r=1:nroffactors
           if isempty(find(combinations(combinationzaehler,:)==r))
               deletevec(deletezaehler)=r;
               deletezaehler=deletezaehler+1;
           end
        end
        factorstructure_ia(deletevec,:)=[];
        if isempty(factorstructure_ia)
            conditionindexstring = [conditionindexstring,',1'];
        end
        evalstring = ['conditionvec(q,r)=conditionmat(',conditionindexstring,');'];
        if ~isempty(factorstructure_ia)
            for q=1:size(conditionvec,1)
                for r=1:size(factorstructure_ia,2)
                    eval(evalstring);
                end
            end
        elseif ~isempty(betweenfactorstructure_ia)
            for q=1:size(conditionvec,1)
                for r=1:size(betweenfactorstructure_ia,2)
                    eval(evalstring);
                end
            end
        end
    else
        conditionvec = 1;
    end
    
    
    
    %calculate conditionmeans, QS & df |  calculate test-qs
    %-------------------------------------------------------
    ConditionQS = 0;
    TestQS = 0;
    InteractionSubjectSum=zeros(nrofsubjects,size(InteractionSum,2));
    
    %reines between design
    %========================
    if isempty(factorstructure) & ~isempty(betweenfactorstructure) 
        
        
        ConditionQS = 0;
        TestQS = 0;
        TestDf = 0;
        for r=1:nrofinteractioncells
            relevantcells = find(conditionvec(:,1)==r);
            InteractionSum(r)=sum(CellGroupSum(find(conditionvec==r)));
            TestDf = TestDf + (length([betweensubjectcell{relevantcells}])-1);
            InteractionSubjectSum(betweensubjectcell{r},1) = sum(anovadata(betweensubjectcell{r},:),2);
        end
        results{combinationzaehler,3} = 1;
        for z=1:size(betweenfactorstructure_ia,2)
            results{combinationzaehler,3} = results{combinationzaehler,3} * (length(finddiffel(betweenfactorstructure_ia(:,z)))-1);
        end
        ConditionMS = [];
        TestMS = [];
        if PrintSumStatus
            InteractionSum = InteractionSum
            InteractionSubjectSum = InteractionSubjectSum
            InteractionSumSquare = sum(InteractionSum(:).^2) / ((size(anovadata,2))*(size(anovadata,1)/nrofinteractioncells));
            fprintf(1,['EffectSumSquare = ',num2str(InteractionSumSquare,'%5.5f'),char(10)]);
            ErrorSumSquare = SubjectQSf ;
            fprintf(1,['ErrorSumSquare = ',num2str(ErrorSumSquare,'%5.5f'),char(10)]);
            fprintf(1,'\n\n');
        else
            InteractionSumSquare = sum(InteractionSum(:).^2) / ((size(anovadata,2))*(size(anovadata,1)/nrofinteractioncells));
            ErrorSumSquare = SubjectQSf ; 
        end
        ConditionQSf = InteractionSumSquare - TotQSf;
		if length(ConditionQSf)>1
            asdfasdf;
		end
		ErrorQSf = ErrorSumSquare;
         if length(ErrorQSf)>1
            asdfasdf;
		end
   
     
    %reines within desing
    %========================   
    elseif isempty(betweenfactorstructure) 
        
        for r=1:nrofinteractioncells
            InteractionSum(r)=sum(CellSum(find(conditionvec==r)));
            InteractionSubjectSum(:,r)=sum(anovadata(:,find(conditionvec==r)),2);
        end
         %ConditionDf
        if n~=1
            results{combinationzaehler,3}=prod(size(InteractionSum)-1);%df
        else
            results{combinationzaehler,3}=length(InteractionSum)-1;
        end
        ConditionQSf = sum(InteractionSum(:).^2) / ((size(anovadata,2)/nrofinteractioncells)*size(anovadata,1)) - TotQSf;
		if length(ConditionQSf)>1
            asdfasdf;
		end
		ErrorQSf = sum(InteractionSubjectSum(:).^2)/ ((size(anovadata,2)/nrofinteractioncells)) - TotQSf;
         if length(ErrorQSf)>1
            asdfasdf;
		end
        
        
    %gemischtes design
    %=======================
    else
        if isempty(factorstructure_ia) %reiner between effect
            ConditionQS = 0;
            TestQS = 0;
            TestDf = 0;
            for r=1:nrofinteractioncells
                relevantcells = find(conditionvec(:,1)==r);
                InteractionSum(r)=sum(CellGroupSum(find(conditionvec==r)));
                TestDf = TestDf + (length([betweensubjectcell{relevantcells}])-1);
                InteractionSubjectSum(betweensubjectcell{r},1) = sum(anovadata(betweensubjectcell{r},:),2);
            end
            results{combinationzaehler,3} = 1;
            for z=1:size(betweenfactorstructure_ia,2)
                results{combinationzaehler,3} = results{combinationzaehler,3} * (length(finddiffel(betweenfactorstructure_ia(:,z)))-1);
            end
            ConditionMS = [];
            TestMS = [];
            if PrintSumStatus
                InteractionSum = InteractionSum
                InteractionSubjectSum = InteractionSubjectSum
                InteractionSumSquare = sum(InteractionSum(:).^2) / ((size(anovadata,2))*(size(anovadata,1)/nrofinteractioncells));
                fprintf(1,['EffectSumSquare = ',num2str(InteractionSumSquare,'%5.5f'),char(10)]);
                ErrorSumSquare = SubjectQSf ;
                fprintf(1,['ErrorSumSquare = ',num2str(ErrorSumSquare,'%5.5f'),char(10)]);
                fprintf(1,'\n\n');
            else
                InteractionSumSquare = sum(InteractionSum(:).^2) / ((size(anovadata,2))*(size(anovadata,1)/nrofinteractioncells));
                ErrorSumSquare = SubjectQSf ; 
            end
            ConditionQSf = InteractionSumSquare - TotQSf;
			if length(ConditionQSf)>1
                asdfasdf;
			end
			ErrorQSf = ErrorSumSquare;
             if length(ErrorQSf)>1
                asdfasdf;
			end
            
            
        elseif isempty(betweenfactorstructure_ia) %reiner within effect
            
             
            InteractionSubjectSum=zeros(nrofsubjects,size(InteractionSum,1));  
            InteractionBetweenSubjectSum = zeros(nrofsubjects,size(InteractionSum,1));
            for r=1:nrofinteractioncells
                InteractionSum(r)=sum(CellSum(find(conditionvec==r)));
            end
            for r=1:nrofwithininteractioncells
                InteractionSubjectSum(:,r) = sum(anovadata(:,find(conditionvec(1,:)==r)),2);
            end
            %ConditionDf
            if n~=1
                results{combinationzaehler,3}=prod(size(InteractionSum)-1);%df
            else
                results{combinationzaehler,3}=length(InteractionSum)-1;
            end
            ConditionDf = results{combinationzaehler,3};
            if PrintSumStatus
                InteractionSum = InteractionSum
                InteractionSubjectSum = InteractionSubjectSum
                InteractionSumSquare = sum(InteractionSum(:).^2) / ((size(anovadata,2)/nrofinteractioncells)*size(anovadata,1));
                fprintf(1,['EffectSumSquare = ',num2str(InteractionSumSquare,'%5.5f'),char(10)]);
                ErrorSumSquare = sum(InteractionSubjectSum(:).^2)/ ((size(anovadata,2)/nrofinteractioncells)) ;
                fprintf(1,['ErrorSumSquare = ',num2str(ErrorSumSquare,'%5.5f'),char(10)]);
                fprintf(1,'\n\n');
            else
                InteractionSumSquare = sum(InteractionSum(:).^2) / ((size(anovadata,2)/nrofinteractioncells)*size(anovadata,1));
                ErrorSumSquare = sum(InteractionSubjectSum(:).^2)/ ((size(anovadata,2)/nrofinteractioncells)) ;
            end
            ConditionQSf = InteractionSumSquare - TotQSf;
			if length(ConditionQSf)>1
                asdfasdf;
			end
			ErrorQSf = ErrorSumSquare - TotQSf;
             if length(ErrorQSf)>1
                asdfasdf;
			end
        
            
        else %gemischter effect
            
            for r=1:nrofinteractioncells
                InteractionSum(r)=sum(CellGroupSum(find(conditionvec==r)));
            end
            for q=1:nrofbetweeninteractioncells
                for r=1:nrofwithininteractioncells
                    InteractionSubjectSum(betweensubjectcell{q},r) = sum(anovadata(betweensubjectcell{q},find(conditionvec(1,:)==r)),2);
                end
            end
            %ConditionDf
            if n~=1
                results{combinationzaehler,3}=prod(size(InteractionSum)-1);%df
            else
                results{combinationzaehler,3}=length(InteractionSum)-1;
            end
            ConditionDf = results{combinationzaehler,3};
            if PrintSumStatus
                InteractionSum = InteractionSum
                InteractionSubjectSum = InteractionSubjectSum
                InteractionSumSquare = sum(InteractionSum(:).^2) / ((size(anovadata,2)/nrofwithininteractioncells)*(size(anovadata,1)/nrofbetweeninteractioncells));
                fprintf(1,['EffectSumSquare = ',num2str(InteractionSumSquare,'%5.5f'),char(10)]);
                ErrorSumSquare = sum(InteractionSubjectSum(:).^2)/ ((size(anovadata,2)/nrofwithininteractioncells)) ;
                fprintf(1,['ErrorSumSquare = ',num2str(ErrorSumSquare,'%5.5f'),char(10)]);
                fprintf(1,'\n\n');
            else
                InteractionSumSquare = sum(InteractionSum(:).^2) / ((size(anovadata,2)/nrofwithininteractioncells)*(size(anovadata,1)/nrofbetweeninteractioncells));
                ErrorSumSquare = sum(InteractionSubjectSum(:).^2)/ ((size(anovadata,2)/nrofwithininteractioncells)) ;
            end
            ConditionQSf = InteractionSumSquare - TotQSf;
			if length(ConditionQSf)>1
                asdfasdf;
			end
			ErrorQSf = ErrorSumSquare - TotQSf;
             if length(ErrorQSf)>1
                asdfasdf;
			end  
        end
    end
    
    %ConditionQS
    results{combinationzaehler,2}=ConditionQSf;
    %ErrorQs
    results{combinationzaehler,5}=ErrorQSf;
    %ErrorDf
    if ~isempty(betweensubjectcell)
        tmpdf = length([betweensubjectcell{:}]) - length(betweensubjectcell);
        if n~=1
            results{combinationzaehler,6}=prod([(size(InteractionSum)-1),tmpdf]);%errordf
        else
            results{combinationzaehler,6}=prod([length(InteractionSum)-1,tmpdf]);%errordf
        end
        testdf =  results{combinationzaehler,6}; 
    else
        if n~=1
            results{combinationzaehler,6}=prod([(size(InteractionSum)-1),nrofsubjects-1]);%errordf
        else
            results{combinationzaehler,6}=prod([length(InteractionSum)-1,nrofsubjects-1]);%errordf
        end
        testdf =  results{combinationzaehler,6}; 
    end  
end


return