function [results]=calc_rm_anova(anovadata,factor,betweenfactor,n,pointwisestatus,
                                    CellMean,CellGroupMean,SubjectMean,GrandMean,factorstructure,betweenfactorstructure,
                                    betweensubjectcell,nroffactors,nrofbetweenfactors,nrofsubjects,TotQS,TotDf,SubjectQS,SubjectDf,GroupMean)
%   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.                                                                      


if nargin<5;pointwisestatus = 0;end
global hPlot3dMenuFig hPlot2dMenuFig
hPlot3dList=get(hPlot3dMenuFig,'Userdata');

%   'calc_nway_interaction(anovadata,nrofsubjects,validfactor,betweenfactor,n)'
%   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 grouping factors
%   
%   'results' is a cell array containing the following information in the
%   columns with the denoted indices:
%
%   1:Name 2:ConditionQS(SS) 3:ConditionDf 4:ConditionMeanSquare(MS) 5:ErrorQS** 
%   6:ErrorDf 7:ErrorMeanSquare(MS)** 8:F-Value* 9:p-Value*
%   10:epsilon(huynh-feldt) 11: p-Value(huynh-feldt)
%
%   *Values are added in 'rmanova' after calculation of all SSs.
%   **ErrorQS is the SS for the 'tested effect' x 'subject'. From this,
%   group-factor-main-effect-SS and subject-SS will be subtractet
%   (see Winer (1971)).
%
%   'results{3,2}' for instance contains the ConditionQS for the third
%   tested effect.

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)
    
    results{combinationzaehler,1}=combinations(combinationzaehler,:);
    betweenstatus = 0;
    for e=1:nrofbetweenfactors
        if any(combinations(combinationzaehler,:)==(nroffactors+e))
            betweenstatus = 1;
        end
    end
    
    
    %effect without a between component
    if ~betweenstatus
        
        factorstructure_ia = factorstructure;
        komponenten = cell(n,2);
        nrofinteractioncells = 1;
        sizestring = '';
        conditionindexstring = '';
        for komponentenzaehler=1:n
            %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(komponentenzaehler),',r)'];
            else
                sizestring = [num2str(komponenten{komponentenzaehler,2})];
                conditionindexstring = ['factorstructure_ia(',num2str(komponentenzaehler),',r)'];
            end
            nrofinteractioncells = nrofinteractioncells * komponenten{komponentenzaehler,2};
        end
        
        if n==1
            sizestring = ['1,',sizestring];
        end
        evalstring = ['InteractionMean = zeros(',sizestring,');']; %--->InteractionMean = zeros(.....);
        eval(evalstring);
        
        %conditionmat
        %--------------
        conditionmat = zeros(size(InteractionMean));
        conditionmat(1:end)=linspace(1,nrofinteractioncells,nrofinteractioncells);
        
        %conditionvec
        %---------------
        conditionvec = zeros(1,size(anovadata,2));
        %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,:)=[];
        evalstring = ['conditionvec(1,r)=conditionmat(',conditionindexstring,');'];
        for r=1:size(factorstructure_ia,2)
            eval(evalstring);
        end
        
        
        
        %calculate conditionmeans, QS & df |  calculate test-qs
        %-------------------------------------------------------
        ConditionQS = 0;
        TestQS = 0;
        InteractionSubjectMean=zeros(nrofsubjects,size(InteractionMean,2));
        for r=1:nrofinteractioncells
            %InteractionMean
            %-------------------
            InteractionMean(r)=mean(CellMean(find(conditionvec==r)));
            %InteractionSubjectMean
            %----------------------------
            InteractionSubjectMean(:,r)=mean(anovadata(:,find(conditionvec==r)),2);
            %Calculate SS
            %---------------
            ConditionQS = ConditionQS + nrofsubjects*length(find(conditionvec==r))*((InteractionMean(r)-GrandMean)^2);
            for s=1:nrofsubjects
                TestQS = TestQS + length(find(conditionvec==r))*(InteractionSubjectMean(s,r)-GrandMean)^2;
            end
        end
        
          
        %ConditionQS
        results{combinationzaehler,2}=ConditionQS;
        
        %ConditionDf
        if n~=1
            results{combinationzaehler,3}=prod(size(InteractionMean)-1);%df
        else
            results{combinationzaehler,3}=length(InteractionMean)-1;
        end
        conditiondf = results{combinationzaehler,3};
        
        %ConditionMS
        results{combinationzaehler,4}=results{combinationzaehler,2}/results{combinationzaehler,3};
        %ErrorQs
        results{combinationzaehler,5}=TestQS;
        
        %ErrorDf
        if ~isempty(betweensubjectcell)
            tmpdf = length([betweensubjectcell{:}]) - length(betweensubjectcell);
            if n~=1
                results{combinationzaehler,6}=prod([(size(InteractionMean)-1),tmpdf]);%errordf
            else
                results{combinationzaehler,6}=prod([length(InteractionMean)-1,tmpdf]);%errordf
            end
            testdf =  results{combinationzaehler,6}; 
        else
            if n~=1
                results{combinationzaehler,6}=prod([(size(InteractionMean)-1),nrofsubjects-1]);%errordf
            else
                results{combinationzaehler,6}=prod([length(InteractionMean)-1,nrofsubjects-1]);%errordf
            end
            testdf =  results{combinationzaehler,6}; 
        end
       
        %ErrorMS
        %results{combinationzaehler,7}=results{combinationzaehler,5}/results{combinationzaehler,6};
        %F-Value
        %results{combinationzaehler,8}=(ConditionQS/conditiondf)/(TestQS/testdf);% F-Value
        %f_value = results{combinationzaehler,8};
        %p-Value
        %results{combinationzaehler,9}=fcdf(f_value,conditiondf,testdf); %p
        
        
        
     %effect with a between component
     else
        
        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 = ['InteractionMean = zeros(',sizestring,');']; %--->InteractionMean = zeros(.....);
        eval(evalstring);
        
        
        %conditionvec
        %---------------
        conditionmat = zeros(size(InteractionMean));
        conditionmat(1:end)=linspace(1,nrofinteractioncells,nrofinteractioncells);
        if length(conditionmat)>1
            betweenfactorstructure_ia = betweenfactorstructure_ia(:,keeparray);
            conditionvec = zeros(size(betweenfactorstructure_ia,1),size(anovadata,2));
            %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,');'];
            for q=1:size(conditionvec,1)
                for r=1:size(factorstructure_ia,2)
                    eval(evalstring);
                end
            end
        else
            conditionvec = 1;
        end
        
        

        %calculate conditionmeans, QS & df |  calculate test-qs
        %-------------------------------------------------------
        % testqs for between interactions is (effect without between component X subject)
        % testqs for between maineffects is (SubjectQS - MaineffectQS)
        
        if isempty(factorstructure_ia) %reiner between
            ConditionQS = 0;
            TestQS = 0;
            TestDf = 0;
            for r=1:nrofinteractioncells
                relevantcells = find(conditionvec(:,1)==r);
                InteractionMean(r)=mean(CellGroupMean(find(conditionvec==r)));
                ConditionQS = ConditionQS + length([betweensubjectcell{relevantcells}])*length(find(conditionvec==r))*((InteractionMean(r)-GrandMean)^2);
                TestDf = TestDf + (length([betweensubjectcell{relevantcells}])-1)
            end
            ConditionDf = 1;
            for z=1:size(betweenfactorstructure_ia,2)
                ConditionDf = ConditionDf * (length(finddiffel(betweenfactorstructure_ia(:,z)))-1);
            end
            TestQS = SubjectQS;
            ConditionMS = [];
            TestMS = [];
        else
            ConditionQS = 0;
            ConditionDf = 0;
            TestQS = 0;
            TestDf = 0;
            for r=1:nrofinteractioncells
                [relevantcellsrows,relevantcellscols] = find(conditionvec==r);
                messwerte = 0;
                for t=1:length(relevantcellsrows)
                    messwerte = messwerte + length(betweensubjectcell{relevantcellsrows(t)});
                end
                %InteractionMean
                %-------------------
                InteractionMean(r)=mean(CellGroupMean(find(conditionvec==r)));
                ConditionQS = ConditionQS + messwerte*((InteractionMean(r)-GrandMean)^2);
                
            end
            ConditionDf = prod(size(InteractionMean)-1);
            ConditionMS = [];
            TestQS = []; %fill later with non-between analogon
            TestDf = 0; %different degrees of freedom than pure within effect!!!!!
            for z=1:size(betweenfactorstructure_ia,1)
                TestDf = TestDf + (length(betweensubjectcell{z})-1);
            end
            TestDf = TestDf * ConditionDf;
            
            TestMS = [];
        end
        
        %ConditionQS
        results{combinationzaehler,2}=ConditionQS;
        %ConditionDf
        results{combinationzaehler,3}=ConditionDf;
        %ConditionMS
        results{combinationzaehler,4}=ConditionMS;
        %ErrorQs
        results{combinationzaehler,5}=TestQS;
        %ErrorDf
        results{combinationzaehler,6}=TestDf;
        %ErrorMS
        results{combinationzaehler,7}=TestMS;
        %F-Value
        results{combinationzaehler,8}=[];% F-Value
        %p-Value
        results{combinationzaehler,9}=[]; %p
        
       
         
     end
         
end
    
    
    
        



return