/* one-way ANOVA */

aov1(x,m,lr,nr,  _ssa,_sse,_sst,_msa,_mse,_f0)

double x[];
int m,lr,nr[];
double *_ssa,*_sse,*_sst,*_msa,*_mse,*_f0;
{
	int i,k,ip,nn;
	double ssa,sse,sst,msa,mse,f0;
	double w,ww,sa[20],st;
	int nua,nue,nut;
	nn=0;
	ww=0.0;
	for(i=0;i<m;i++)
	{
		w=0.0;
		for(k=0;k<nr[i];k++)
		{
			w+=x[i*lr+k];
			nn++;
		}
		sa[i]=w/nr[i];
		ww+=sa[i];
	}
	st=ww/m;
	w=0.0;
	sse=sst=0.0;
	for(i=0;i<m;i++)
	{
		for(k=0;k<nr[i];k++)
		{
			ip=i*lr+k;
			sse+=(x[ip]-sa[i])*(x[ip]-sa[i]);
			sst+=(x[ip]-st)*(x[ip]-st);
		}
	}
	ssa=sst-sse;
	nua=m-1;
	nut=nn-1;
	nue=nut-nua;
	msa=ssa/nua;
	mse=sse/nue;
	f0=msa/mse;
	*_ssa=ssa;
	*_sse=sse;
	*_sst=sst;
	*_msa=msa;
	*_mse=mse;
	*_f0=f0;
	return(0);
}




/***************************************************/
/**  2-way ANOVA for single   experiment          **/
/**                                               **/
/**    function    : aov2()                       **/
/**    input  x    : 2-dimensional array          **/
/**                  x[factor1][factor2]          **/
/**           m,n,r: #level1,#level2,#repeat      **/
/**           ln   : array size of n   direction  **/
/**    output _ssa,_ssb,_sse,_sst :               **/
/**                  squared sum of factor1,2,    **/
/**                            residual,total     **/
/**           _msa,_msb,_mse :                    **/
/**                  mean squared sum             **/
/**           _f0a,_f0b      : F ratios           **/
/**    return    0 : normal return                **/
/**            999 : parameter error              **/
/**                                               **/
/**    programmed by M. Fujita,'93/Nov/08         **/
/***************************************************/

aov2(x,m,n,ln,
	_ssa,_ssb,_sse,_sst, _msa,_msb,_mse,
	_f0a,_f0b)

double x[];
int m,n,ln;
double *_ssa,*_ssb,*_sse,*_sst,*_msa,*_msb,*_mse,
	*_f0a,*_f0b;
{
	int i,j;
	double ssa,ssb,sse,sst,msa,msb,mse,mst,f0a,f0b;
	double tttr[20][20],ttta[20][20];
	double w,ww,sa[20],sb[20],se[20],st;
	int nua,nub,nue,nut;

	if(ln<1||m<1||m>500||n<1||n>500) return(999);

	for(i=0;i<m;i++)                             /* sa[i] */
	{
		w=0.0;
		for(j=0;j<n;j++) w+=x[i*ln+j];
		sa[i]=w/n;
	}
	for(j=0;j<n;j++)                             /* sb[j] */
	{
		w=0.0;
		for(i=0;i<m;i++) w+=x[i*ln+j];
		sb[j]=w/m;
	}
	w=0.0;                                       /* st */
	for(i=0;i<m;i++) w+=sa[i];
	st=w/m;
						 /* squared sum */
	w=0.0;                                       /* ssa */
	for(i=0;i<m;i++) w+=(sa[i]-st)*(sa[i]-st);
	ssa=w*n;
	w=0.0;                                       /* ssb */
	for(j=0;j<n;j++) w+=(sb[j]-st)*(sb[j]-st);
	ssb=w*m;
	w=0.0;                                       /* sst */
	for(i=0;i<m*n;i++) w+=(x[i]-st)*(x[i]-st);
	sst=w;
	sse=sst-ssa-ssb;                        /* sse */

						 /* degree of freedom */
	nua=m-1;
	nub=n-1;
	nut=m*n-1;
	nue=nut-nua-nub;
						 /* mean squared sum */
	msa=ssa/nua;                                 /* msa */
	msb=ssb/nub;                                 /* msb */
	mse=sse/nue;                                 /* mse */
						 /* F ratios */
	f0a=msa/mse;                                 /* f0a */
	f0b=msb/mse;                                 /* f0b */
						 /* return result */
	*_ssa=ssa;
	*_ssb=ssb;
	*_sse=sse;
	*_sst=sst;
	*_msa=msa;
	*_msb=msb;
	*_mse=mse;
	*_f0a=f0a;
	*_f0b=f0b;

	return(0);
}




/***************************************************/
/**  2-way ANOVA for repeated experiment          **/
/**                                               **/
/**    function    : aov2r()                      **/
/**    input  x    : 3-dimensional array          **/
/**                  x[factor1][factor2][repeat]  **/
/**           m,n,r: #level1,#level2,#repeat      **/
/**           ln,lr: array size of n,r direction  **/
/**    output _ssa,_ssb,_ssab,_sse,_sst :         **/
/**                  squared sum of factor1,2,    **/
/**                  factor1*2,residual,total     **/
/**           _msa,_msb,_msab,_mse :              **/
/**                  mean squared sum             **/
/**           _f0a,_f0b,_f0ab : F ratios          **/
/**    return    0 : normal return                **/
/**            999 : parameter error              **/
/**                                               **/
/**    programmed by M. Fujita,'93/Nov/04         **/
/***************************************************/

aov2r(x,m,n,r,ln,lr,
	_ssa,_ssb,_ssab,_sse,_sst, _msa,_msb,_msab,_mse,
	_f0a,_f0b,_f0ab)

double x[];
int m,n,r,ln,lr;
double *_ssa,*_ssb,*_ssab,*_sse,*_sst,*_msa,*_msb,*_msab,*_mse,
	*_f0a,*_f0b,*_f0ab;
{
	int i,j,k;
	double ssa,ssb,ssab,sse,sst,msa,msb,msab,mse,mst,f0a,f0b,f0ab;
	double tttr[20][20],ttta[20][20];
	double w,ww,sa[20],sb[20],se[20],st;
	int nua,nub,nuab,nue,nut;

	if(ln<1||lr<1||m<1||m>500||n<1||n>500||r<1||r>500) return(999);

						 /* sum on each axis */
	for(i=0;i<m;i++)                             /* proj. AB plane */
	for(j=0;j<n;j++)
	{
		w=0.0;
		for(k=0;k<r;k++) w+=x[ln*lr*i+lr*j+k];
		tttr[i][j]=w;
	}
	for(j=0;j<n;j++)                             /* proj. Be plane */
	for(k=0;k<r;k++)
	{
		w=0.0;
		for(i=0;i<m;i++) w+=x[ln*lr*i+lr*j+k];
		ttta[j][k]=w;
	}
	for(i=0;i<m;i++)                             /* sa[i] */
	{
		w=0.0;
		for(j=0;j<n;j++) w+=tttr[i][j];
		sa[i]=w/(n*r);
	}
	for(j=0;j<n;j++)                             /* sb[j] */
	{
		w=0.0;
		for(i=0;i<m;i++) w+=tttr[i][j];
		sb[j]=w/(m*r);
	}
	for(k=0;k<r;k++)                             /* se[k] */
	{
		w=0.0;
		for(j=0;j<r;j++) w+=ttta[j][k];
		se[k]=w/(m*n);
	}
	w=0.0;                                       /* st */
	for(i=0;i<m;i++) w+=sa[i];
	st=w/m;
						 /* squared sum */
	w=0.0;                                       /* ssa */
	for(i=0;i<m;i++) w+=(sa[i]-st)*(sa[i]-st);
	ssa=w*n*r;
	w=0.0;                                       /* ssb */
	for(j=0;j<n;j++) w+=(sb[j]-st)*(sb[j]-st);
	ssb=w*m*r;
	w=0.0;                                       /* ssab */
	for(i=0;i<m;i++)
	for(j=0;j<n;j++)
	{
		ww=tttr[i][j]/r-sa[i]-sb[j]+st;
		w+=ww*ww;
	}
	ssab=w*r;
	w=0.0;                                       /* sst */
	for(k=0;k<m*n*r;k++) w+=(x[k]-st)*(x[k]-st);
	sst=w;
	sse=sst-ssa-ssb-ssab;                        /* sse */

						 /* degree of freedom */
	nua=m-1;
	nub=n-1;
	nuab=nua*nub;
	nut=m*n*r-1;
	nue=nut-nua-nub-nuab;
						 /* mean squared sum */
	msa=ssa/nua;                                 /* msa */
	msb=ssb/nub;                                 /* msb */
	msab=ssab/nuab;                              /* msab */
	mse=sse/nue;                                 /* mse */
						 /* F ratios */
	f0a=msa/mse;                                 /* f0a */
	f0b=msb/mse;                                 /* f0b */
	f0ab=msab/mse;                               /* f0ab */
						 /* return result */
	*_ssa=ssa;
	*_ssb=ssb;
	*_ssab=ssab;
	*_sse=sse;
	*_sst=sst;
	*_msa=msa;
	*_msb=msb;
	*_msab=msab;
	*_mse=mse;
	*_f0a=f0a;
	*_f0b=f0b;
	*_f0ab=f0ab;

	return(0);
}