title 'Polynomial contrasts';
proc iml;
	reset print log fw=6;
	
	*-- vector of treatment levels;
	X = {1, 2, 3, 4 };
	V = j(4,1) || x || x##2 || x##3;

	*-- orthogonalize with Gram-Schmidt;	
	call gsorth(C, T, flag, V);
	
	*-- or, generate them directly;
	C = orpol( X );

	*-- 'unnormalize' to get the usual values';
	C = C * diag( sqrt({4 20 4 20}) );
	
	*-- can use unequally-spaced treatment levels;
	X = {1 2 4 9};
	C = orpol( X ) * sqrt(4);
	
	*-- and/or unequal sample sizes;
	N = {10 10 20 20};
	X = {1 2 4 9};
	C = orpol( X, 3, N ) * sqrt(60);