#=============================================================================
#
# Some useful Maple procedures.
#
# Marc A. Murison
# Astronomical Applications Dept.
# U.S. Naval Observatory
# 3450 Massachusetts Ave., NW
# Washignton, DC 20392
# murison@riemann.usno.navy.mil
# http://riemann.usno.navy.mil/AESOP/
#
#=============================================================================

if UTILS_READ <> 1 then
UTILS_READ := 1;

utils_help := proc()
	print(`utils: [ annular_average, applyfunc, CombineTrig, CombTrigOpt, `
		 .`common_factor, \ncosfix, coshfix, csquare, debug_print, expansion, `
		 .`funcops, FunctionCalls, getsqrts, \ninvalgsubs, invsubs, is_real_parts, polytest, `
		 .`pullout, rootfunc, RotAxis, signsqrt, sinfix, sinhfix, \nsmall_divisors, `
		 .`termfunc, tmsg, topsqrt, topsqrts ]`);
end:

#utils_help();


#=============================================================================
#
# Print messages of the form
#
# <proc_name> [<elapsed time>]: <msg>
#
# to the terminal screen, where elapsed time is the current
# time minus the time stored in the optional third argument.
#
#=============================================================================
tmsg := proc( proc_name, msg )
    local dt;
	if nargs=3 then
    	dt := time() - args[3];
	else
		dt := time();
	fi;
	if proc_name <> `` then
		if dt < 1 then
			printf(proc_name.`[0]: %s\n`,msg);
		else
			printf(proc_name.`[%.0f]: %s\n`,dt,msg);
		fi;
	else
		if dt < 1 then
			printf(`[0] %s\n`,msg);
		else
			printf(`[%.0f] %s\n`,dt,msg);
		fi;
	fi;
end:



#=============================================================================
#
# pname = procname in calling procedure
#
#=============================================================================
debug_print := proc( pname, msg, verbosity_level )

	global verbosity, time0;

	if verbosity_level=0 then
		_debug_print( pname, msg, args[3..nargs] );
	elif assigned( verbosity ) then
		if verbosity_level <= verbosity then
			_debug_print( pname, msg, args[3..nargs] );
		fi;
	fi;
end:

_debug_print := proc( pname, msg );

	if assigned( time0 ) then
		tmsg( pname, msg, time0 );
	else
		tmsg( pname, msg );
	fi;
	if nargs > 3 then
		print( eval(args[4..nargs]) );
	else
		print();					#force a flush
	fi;

end:



#=============================================================================
#
# Remove all top-level subscripts, e.g. v[x] -> vx.
#
#=============================================================================
remove_subscripts := proc( expr )

	local	p, p1, p2, loc, Q;
	global	_subscripts_, clear_subscripts_, rlevel;

	if not assigned(clear_subscripts_) or clear_subscripts_=true then
		_subscripts_ := {};
	fi;

	Q := copy(expr);
	for p in Q do
		if type(p,integer) then
			next;
		fi;
		loc := location(Q,p);
		if loc=[] then
			debug_print(procname,`Unable to locate `,0,p);
			next;
		fi;
		if type(p,indexed) then
			p1 := op(0,p);
			p2 := op(1,p);
			if nops(p)=1 and nops(p2)=1 and not type(p2,indexed) then
				_subscripts_ := _subscripts_ union {``.p1.p2=p};
				p := ``.p1.p2;
			else
				next;
			fi;
		elif nops(p) > 1 then
			clear_subscripts_ := false;
			p := procname(p);
			clear_subscripts_ := true;
		else
			next;
		fi;
		Q := subsop( loc=p, Q );
	od;

	eval(Q);

end:



#=============================================================================
#
# Restore top-level subscripts, e.g. vx -> v[x] [see remove_subscripts()].
#
#=============================================================================
restore_subscripts := proc( expr )

	global	_subscripts_;

	subs( _subscripts_, expr );

end:



#=============================================================================
#
# A better subs.  If optional arguments are passed, then they are interpreted
# as a function (and further args for that function) to be applied to the
# expression after the substitution has been made.  If that is the case,
# the substitution is undone after the function has been applied.  If no
# optional function args are given, then a normal subs() is performed.
#
# Example: Subs( foo, x[0]=x0, expansion, x0, 4 );
#
#=============================================================================
Subs := proc( expr::anything, subslist::{`=`,list(`=`),set(`=`)} )

	local fargs, Q;

	Q := subs( subslist, expr );

	fargs := [args[3..nargs]];
	if nops(fargs) > 0 then
		Q := fargs[1](Q,op(fargs[2..nops(fargs)]));
		Q := invsubs( subslist, Q );
	fi;

	Q;
end:



#=============================================================================
#
# Do an "inverse" subs().  That is, for each x=y in subslist, do a
# substitution in expr of y=x.
#
#=============================================================================
invsubs := proc( subslist::{`=`,list(`=`),set(`=`)}, expr::anything )

	local oldsubs, newsubs, p;

	oldsubs := copy(subslist);
	if type(oldsubs,`=`) then
		oldsubs := [oldsubs];
	fi;

	newsubs := [];
	for p in oldsubs do
		newsubs := [op(newsubs),rhs(p)=lhs(p)];
	od;

	if type(oldsubs,set) then
		newsubs := convert(newsubs,set);
	fi;

	subs( newsubs, expr );
end:



#=============================================================================
#
# Do an "inverse" algsubs().  That is, for each x=y in subslist, do an
# "algebraic" substitution in expr of y=x.  An extension of the
# algsubs() syntax is that for invalgsubs() you can specify more than
# one substitution, as in subs() and invsubs().
#
#=============================================================================
invalgsubs := proc( subslist::{`=`,list(`=`),set(`=`)}, expr::anything )

	local oldsubs, newsubs, p, newexpr;

	oldsubs := copy(subslist);
	if type(oldsubs,`=`) then
		oldsubs := [oldsubs];
	fi;

	newsubs := [];
	for p in oldsubs do
		newsubs := [op(newsubs),rhs(p)=lhs(p)];
	od;

	newexpr := copy(expr);
	for p in newsubs do
		newexpr := algsubs( p, newexpr );
	od;

	newexpr;
end:



#=============================================================================
#
# Make column vector substitutions.  Example:
#
# > conicoid := z-(x^2+y^2)/(2*f+sqrt(4*f^2 - epsilon*(x^2+y^2)));
#
#                                         2    2
#                                        x  + y
#     conicoid := z - ---------------------------------------------
#                               2            2            2 1/2
#                     2 f + (4 f  - epsilon x  - epsilon y )
#
# > mat(x,y,z) = mat(x0,y0,z0) + t*mat(vx,vy,vz);
#
#                               [x]   [x0]     [vx]
#                               [y] = [y0] + t [vy]
#                               [z]   [z0]     [vz]
#
# > rootfunc( matsubs( conicoid, " ), collect, epsilon );
#
#                                       2              2
#                            (x0 + t vx)  + (y0 + t vy)
#   z0 + t vz - --------------------------------------------------------
#                                   2              2               2 1/2
#               2 f + ((-(x0 + t vx)  - (y0 + t vy) ) epsilon + 4 f )
#
#=============================================================================
matsubs := proc( expr, subsvec::`=` )

	local i, k, tmp, s, n;

	s := evalm(subsvec);

	if not type(lhs(s),'array'(1..integer,1..1))
	or not type(rhs(s),'array'(1..integer,1..1)) then
		ERROR(`: 2nd arg must be of the form array(1..n,1..1)=array(1..n,1..1)`);
	fi;

	n := rows(lhs(s));
	for i from 1 to n do
		tmp[i] := lhs(s)[i,1] = rhs(s)[i,1];
	od;

	subs( seq(tmp[k],k=1..n), expr );
end:



#=============================================================================
#
# Find the op() location of term in expr.  location returns a Maple list
# suitable for use in op().  Hence, if it is possible to return
# term via the op() function, then
#
#    op( location(expr,term), expr ) = term
#
#=============================================================================
location := proc( expr::anything, term::anything )

	local L, M, k, parts, i;

	if has(expr,term) then

		if type(expr,table) then
			L := [3];
			parts := [op(op(3,expr))];
		else
			L := [];
			parts := [op(expr)];
		fi;

		#check each whole term of this level of op(expr)
		if member( term, parts, k ) then
			RETURN( [k] );
		else
			#descend to each term of this level from left to right
			for i from 1 to nops(parts) do
				M := procname( parts[i], term );
				if nops(M) > 0 then
					RETURN( [op(L),i,op(M)] );
				fi;
			od;
		fi;
	fi;

	RETURN( [] );
end:



#=============================================================================
#
# A better collect().  Optional argument is 'loc' or 'alg'.  'loc' forces
# use of location() method, while 'alg' forces use of algsubs() method.  The
# default is 'alg'.
#
# Side effects:  global _T_ cleared.
#
#=============================================================================
Collect := proc( expr::{algebraic,list,set,relation,array},
				 subexpr::{algebraic,list,set} )

	local	Z, s, k, j, p, L, form, fargs, skip;
	global	_T_;

	if not type( eval(expr), algebraic ) then
		RETURN(map(procname,args));
	fi;

	if not type(subexpr,{list,set}) then
		s := [subexpr];
	else
		s := subexpr;
	fi;

	#
	# extract collection method and adjust function args
	#
	form := 'alg';
	for p in args[3..nargs] do
		if p='loc' then
			form := 'loc';
		fi;
	od;
	fargs := op( remove( has, {args[3..nargs]}, {'alg','loc'} ) );

	#
	# check to see if we can use Maple's collect()
	#
	select( hastype, s, {`+`,`*`} );
	if " = [] then
		skip := false;
		for p in s do
			if type(p,fractional_power) then
				skip := true;
				break;
			fi;
		od;
		if not skip then
			RETURN( collect(expr,subexpr,fargs) );
		fi;
	fi;

	Z   := expr;
	_T_ := '_T_';

	#
	# dummy substitutions
	#
	for k from 1 to nops(s) do
		if form='loc' then
			Z := `Collect/locsubs`( Z, s, k );
		else
			Z := `Collect/algsubs`( Z, s, k );
		fi;
	od;

	#
	# collect on the dummies
	#
	Z := collect( Z, [seq(_T_[j],j=1..nops(s))], fargs );

	#
	# subs back
	#
	for k from 1 to nops(s) do
		Z := subs( _T_[k]=s[k], Z );
	od;

	Z;

end:

`Collect/locsubs` := proc( Z, s, k )
	local	L, tmp;
	global	_T_;
	#
	# repeat subsop until no more subexpressions in Z
	#
	tmp := Z;
	L   := location( tmp, s[k] );
	while nops(L) > 0 do
		tmp := subsop( L=_T_[k], tmp );
		L   := location( tmp, s[k] );
	od;
	tmp;
end:

`Collect/algsubs` := proc( Z, s, k )
	local	Q, p, H, M, S, j;
	global	_T_;

	Q := Z;

	#
	# if any of the subexpressions is the argument of a diff,
	# then we have to substitute a dummy or else a bug
	# in algsubs will delete the derivative from the expression!
	#
	if hastype(Q,diff) then
		M := MakeCoeffs(Q);
		S := [];
		j := 0;
		for p in M[derivs] do
			if isdiff(p) and has(p,s[k]) then
				j := j + 1;
				S := [op(S),p];		#save the substituted diffs
				Q := subs(p=H[j],Q);
			fi;
		od;
	fi;

	#
	# now do the algebraic substitutions
	#
	Q := traperror( algsubs( s[k]=_T_[k], Q, exact ) );

	#
	# check if algsubs failed
	#
	if member(`cannot compute degree of pattern in`,{Q})  #algsubs failed
	or member(`no variables in pattern`,{Q}) then
		debug_print(procname,
					`subexpression too complicated for algsubs:`,
					1, lasterror, s[k] );
		debug_print(procname,`trying location() approach...`,1);
		RETURN( `Collect/locsubs`( Z, s, k ) );
	fi;

	#
	# subs back for the diffs
	#
	if nops(S) <> 0 then
		for j from 1 to nops(S) do
			Q := subs( H[j]=S[j], Q );
		od;
	fi;

	Q;
end:



#=============================================================================
#
# Side relations simplification with
#   - temporary replacement of subscripts
#   - collection on the side relations
#   - optional function and function args to apply during collection
#
#=============================================================================
srsimp := proc( expr::{`=`,algebraic,array,list,set},
				siderelslist::{`=`,algebraic,list,set} )

	local srlist, sublist, A, k;

	if not type(expr,algebraic) then
		RETURN( map(procname,args) );
	fi;

	if not type(siderelslist,{list,set}) then
		srlist := [siderelslist];
	elif type(siderelslist,set) then
		srlist := convert(siderelslist,list);
	else
		srlist := copy(siderelslist);
	fi;

	sublist := [seq(srlist[k]=A.k,k=1..nops(srlist))];
	simplify(expr,convert(sublist,set));
	remove_subscripts(");
	Collect( ", [seq(A.k,k=1..nops(srlist))], args[3..nargs] );
	restore_subscripts(");
	invsubs( sublist, " );
end:



#=============================================================================
#
# Define a "fractional_power" type.  An expression is of this type if
# it contains a subexpression raised to a fractional power.
#
#=============================================================================
`type/fractional_power` := proc( expr::{algebraic,equation,list,set,indexed,array} )
	local p;
	if type(expr,`^`) then
		if type(op(2,expr),fraction) then
			RETURN(true);
		fi;
	else
		if nops(p) > 1 then
			for p in expr do
				if procname(p)=true then
					RETURN(true);
				fi;
			od;
		fi;
	fi;
	false;
end:



#=============================================================================
#
# Extracts the top-level terms containing fractional powers and
# returns them in a list.  For example,
#
# > f := solve( a + b*t + c*t^2, t );
# > g := (sqrt(f[2]*A+B)+C)^2;
#
#          /    /          2         1/2         \1/2    \2
#          |    |  (-b - (b  - 4 c a)   ) A      |       |
#      g = |1/2 |2 ------------------------ + 4 B|    + C|
#          \    \            c                   /       /
#
# > getsqrts(g);
#             2         1/2  /               2         1/2      \1/2
#         A (b  - 4 c a)     |   A b     A (b  - 4 c a)         |
#  [- 1/2 -----------------, |-2 --- - 2 ----------------- + 4 B|    C]
#                 c          \    c              c              /
#
#=============================================================================
topsqrts := proc( expr::{algebraic,equation} )
	local x, p, rootlist;
	if type(expr,equation) then
		RETURN( [op(procname(lhs(expr))),op(procname(rhs(expr)))] );
	fi;
	if not hastype(expr,fractional_power) then
		debug_print(procname,`No sqrt terms in expression!`,3);
		RETURN( [] );
	fi;
	rootlist := [];
	x := expand(expr);
	for p in x do
		if hastype(p,fractional_power) then
			rootlist := [op(rootlist),p];
		fi;
	od;
	RETURN(rootlist);
end:



#=============================================================================
#
# Extract the first fractional power term encountered in an expression.
#
#=============================================================================
topsqrt := proc( expr::{algebraic,`=`} )
	local p;
	if not hastype(expr,fractional_power) then
		debug_print(procname,`No sqrt term in expression!`,3,expr);
		RETURN(FAIL);
	fi;
	if type( expr, fractional_power ) then
		expr;
	else
		# check top level terms
		for p in expr do
			if type( p, fractional_power ) then
				RETURN(p);
			fi;
		od;
		# descend subexpressions
		for p in expr do
			if hastype(p,fractional_power) then
				RETURN( procname(p) );
			fi;
		od;
	fi;
	RETURN(FAIL);
end:



#=============================================================================
#
# Return the sign of a term involving a fractional power.
# Use getsqrts() to isolate fractional power terms.  For example,
#
# > f := solve( a + b*t + c*t^2, t );
# > g := (sqrt(f[2]*A+B)+C)^2;
# > h := topsqrts(g);
#
#                   2         1/2
#               A (b  - 4 c a)
#   h := [- 1/2 -----------------,
#                      c
#
#       /               2         1/2      \1/2
#       |   A b     A (b  - 4 c a)         |
#       |-2 --- - 2 ----------------- + 4 B|    C]
#       \    c             c               /
#
# > signsqrt(h[1]);
#
#                                   -1
#
#=============================================================================
signsqrt := proc( expr::{`*`,`^`} )
	local p, cursign;
	if not hastype(expr,fractional_power) then
		debug_print(procname,`No sqrt terms in expression!`,3);
		RETURN(FAIL);
	fi;
	cursign := 1;
	for p in expr do
		cursign := cursign*sign(p);
	od;
	RETURN(cursign);
end:



#=============================================================================
#
# Extracts all fractional power terms in expr and returns them in a
# list structure.  The list structure is of the form [[s1,loc1],[s2,loc2],...]
# where s1, etc. are the terms and loc1, etc. are the expression parse tree
# locations as defined by the op function.  Hence, s1=op(loc1,expr) and so
# on.  For example,
#
# > f := solve( a + b*t + c*t^2, t );
# > g := (sqrt(f[2]*A+B)+C)^2;
#
#          /    /          2         1/2         \1/2    \2
#          |    |  (-b - (b  - 4 c a)   ) A      |       |
#      g = |1/2 |2 ------------------------ + 4 B|    + C|
#          \    \            c                   /       /
#
# > getsqrts(g);
#
#=============================================================================
getsqrts := proc( expr::{algebraic,equation,list,set,indexed,array} )

	local k, p, oplist1, oplist2;

	if nargs > 1 then
		oplist1 := args[2];
		oplist2 := args[3];
	else
		oplist1 := [];
		oplist2 := [];
	fi;

	if type(expr,fractional_power) then
		if oplist1 <> [] then
			oplist2 := [ op(oplist2), oplist1 ];
		fi;
		if hastype( op(1,expr), fractional_power ) then
			oplist2 := [ op(oplist2), procname(op(1,expr),oplist2,[]) ];
		fi;
		RETURN( oplist2 );
	elif hastype(expr,fractional_power) then
		for k from 1 to nops(expr) do
			p := op(k,expr);
			if type(p,fractional_power) then
				oplist2 := [ op(oplist2), [op(oplist1),k] ];
				if hastype( op(1,p), fractional_power ) then
					oplist2 := [ op(oplist2), procname(op(1,p),oplist2,[op(oplist1),k]) ];
				fi;
			elif hastype(p,fractional_power) then
				if oplist2 <> [] then
					oplist2 := [ op(oplist2), procname(p,oplist2,[op(oplist1),k]) ];
				else
					oplist2 := procname(p,oplist2,[op(oplist1),k]);
				fi;
			fi;
		od;
	fi;

	RETURN(oplist2);
end:



#=============================================================================
#
# Recursively apply a procedure, func(), to subexpressions that are raised
# to a fractional power.  This is handy for collecting terms under square
# roots.  For example,
#
#   rootfunc( expr, collect, [var1, var2], factor )
#
# Any arguments beyond the second argument, expr, are treated as additional
# arguments of func().  In the above example, subexpressions are passed
# to collect as collect( expr, [var1, var2], factor ).  rootfunc() is able
# to handle multiply-nested fractional power expressions.
#
#=============================================================================
rootfunc := proc( expr::{algebraic,equation,list,set,indexed,array},
				  func::procedure )
	local p, z, loc;

	if type(expr,array) then
	  RETURN( map(procname,args) );
	fi;

	if type(expr,fractional_power) then
		RETURN( applyop(func,1,expr,args[3..nargs]) );
	fi;

	z := expr;
	for p in expr do
		if nops(p)>1 and hastype(p,fractional_power) then
			loc := location(z,p);
			if type(p,fractional_power) then
				p := applyop(func,1,p,args[3..nargs]);
			else
				p := procname(p,func,args[3..nargs]);
			fi;
			z := subsop(loc=p,z);
		fi;
	od;
	z;

end:



#=============================================================================
#
# Force 1 - cos^2 ==> sin^2.
#
#=============================================================================
cosfix := proc( expr::{algebraic,array,`=`}, var::{algebraic,list} )

	local p, S, tmp, _H_, k, q, Q, i;

	if not type(expr,algebraic) then
		RETURN( map( procname, args ) );
	fi;

	tmp := expr;
	if type(var,list) then
		S := {};
		k := 0;
		Q := 0;
		for p in var do
			if type(p,`+`) then
				k    := k + 1;
				tmp  := subs( p=_H_[k], tmp );
				S    := { op(S), 1-cos(_H_[k])^2=sin(_H_[k])^2 };
				q[k] := p;
				Q    := 1;
			else
				S := { op(S), 1-cos(p)^2=sin(p)^2 };
			fi;
		od;
		tmp := simplify( tmp, S );
		if Q=1 then
			for i from 1 to k do
				tmp := subs( _H_[i]=q[i], tmp );
			od;
		fi;
		tmp;
	elif type(var,`+`) then
		tmp := subs( var=_H_, expr );
		S   := {1-cos(_H_)^2=sin(_H_)^2};
		simplify( tmp, S );
		subs( _H_=var, " );
	else
		S := {1-cos(var)^2=sin(var)^2};
		simplify( tmp, S );
	fi;

end:



#=============================================================================
#
# Force 1 - sin^2 ==> cos^2.
#
#=============================================================================
sinfix := proc( expr::{algebraic,array,`=`}, var::{algebraic,list} )

	local p, S, tmp, _H_, k, q, Q, i;

	if not type(expr,algebraic) then
		RETURN( map( procname, args ) );
	fi;

	tmp := expr;
	if type(var,list) then
		S := {};
		k := 0;
		Q := 0;
		for p in var do
			if type(p,`+`) then
				k    := k + 1;
				tmp  := subs( p=_H_[k], tmp );
				S    := { op(S), 1-sin(_H_[k])^2=cos(_H_[k])^2 };
				q[k] := p;
				Q    := 1;
			else
				S := { op(S), 1-sin(p)^2=cos(p)^2 };
			fi;
		od;
		tmp := simplify( tmp, S );
		if Q=1 then
			for i from 1 to k do
				tmp := subs( _H_[i]=q[i], tmp );
			od;
		fi;
		tmp;
	elif type(var,`+`) then
		tmp := subs( var=_H_, expr );
		S   := {1-sin(_H_)^2=cos(_H_)^2};
		simplify( tmp, S );
		subs( _H_=var, " );
	else
		S := {1-sin(var)^2=cos(var)^2};
		simplify( tmp, S );
	fi;

end:



#=============================================================================
#
# Force cosh^2 - 1 ==> sinh^2.
#
#=============================================================================
coshfix := proc( expr::{algebraic,array,`=`}, var::{algebraic,list} )

	local p, S, tmp, _H_, k, q, Q, i;

	if not type(expr,algebraic) then
		RETURN( map( procname, args ) );
	fi;

	tmp := expr;
	if type(var,list) then
		S := {};
		k := 0;
		Q := 0;
		for p in var do
			if type(p,`+`) then
				k    := k + 1;
				tmp  := subs( p=_H_[k], tmp );
				S    := { op(S), cosh(_H_[k])^2-1=sinh(_H_[k])^2 };
				q[k] := p;
				Q    := 1;
			else
				S := { op(S), cosh(p)^2-1=sinh(p)^2 };
			fi;
		od;
		tmp := simplify( tmp, S );
		if Q=1 then
			for i from 1 to k do
				tmp := subs( _H_[i]=q[i], tmp );
			od;
		fi;
		tmp;
	elif type(var,`+`) then
		tmp := subs( var=_H_, expr );
		S   := {cosh(_H_)^2-1=sinh(_H_)^2};
		simplify( tmp, S );
		subs( _H_=var, " );
	else
		S := {cosh(var)^2-1=sinh(var)^2};
		simplify( tmp, S );
	fi;

end:



#=============================================================================
#
# Force 1 + sinh^2 ==> cosh^2.
#
#=============================================================================
sinhfix := proc( expr::{algebraic,array,`=`}, var::{algebraic,list} )

	local p, S, tmp, _H_, k, q, Q, i;

	if not type(expr,algebraic) then
		RETURN( map( procname, args ) );
	fi;

	tmp := expr;
	if type(var,list) then
		S := {};
		k := 0;
		Q := 0;
		for p in var do
			if type(p,`+`) then
				k    := k + 1;
				tmp  := subs( p=_H_[k], tmp );
				S    := { op(S), 1+sinh(_H_[k])^2=cosh(_H_[k])^2 };
				q[k] := p;
				Q    := 1;
			else
				S := { op(S), 1+sinh(p)^2=cosh(p)^2 };
			fi;
		od;
		tmp := simplify( tmp, S );
		if Q=1 then
			for i from 1 to k do
				tmp := subs( _H_[i]=q[i], tmp );
			od;
		fi;
		tmp;
	elif type(var,`+`) then
		tmp := subs( var=_H_, expr );
		S   := {1+sinh(_H_)^2=cosh(_H_)^2};
		simplify( tmp, S );
		subs( _H_=var, " );
	else
		S := {sinh(var)^2-1=cosh(var)^2};
		simplify( tmp, S );
	fi;

end:



#=============================================================================
#
# Given an expression composed of various trigonometric terms involving
# at least two angles theta[1] and theta[2], combine the terms containing
# sines and cosines of theta[1] and theta[2] into single sine and cosine
# terms.  Other sines and cosines involving other angles are protected
# from the Maple combine() procedure.  The input expression may also be a
# one or two dimensional array.
#
#=============================================================================
CombineTrig := proc( expr::{algebraic,`=`,array}, var1::algebraic, var2::algebraic )

	local p, k, j, sublist, A, kmax, f, f1, f2, g, isfun;

	f := copy(expr);
	if type(f,{'array'(1),'array'(2),`=`}) then
		RETURN( map( procname, f, var1, var2 ) );
	fi;
	if not ( has(var1,indets(f)) and has(var2,indets(f)) ) then
		RETURN(eval(f));
	fi;

	# f is a function, such as arcsin(stuff)
	if type(f,function) then
		g := f;
		f := op(f);			#do the simplification on "stuff"
		isfun := YES;
	else
		isfun := NO;
	fi;

	f := simplify( expand(f) );
	if not hasfun(f,sin) and not hasfun(f,cos) then
		if isfun=YES then
			RETURN(eval(g));
		else
			RETURN(eval(f));
		fi;
	fi;

	#
	# Substitute for sines and cosines of other indeterminates
	# in order to hide them from combine().  Must be careful, since
	# if var1 or var2 are themselves functions of other variables
	# [e.g., lambda(t)] then Maple will classify them as functions;
	# hence the "if hasfun..." condition.
	#
	sublist := array;
	k := 1;
	for p in indets(f,function) do
		if not p=sin(var1) and not p=cos(var1) and
		   not p=sin(var2) and not p=cos(var2) and
		   not p=var1 and not p=var2 then
			if hasfun(p,cos) or hasfun(p,sin) then
				sublist[k] := p;
				k := k + 1;
			fi;
		fi;
	od;
	kmax := k-1;
	f := subs( seq( sublist[k]=A[k], k=1..kmax ), f );

	#
	# Simplify the remaining trig terms with combine()
	#
	f := combine(f);

	#
	# Restore previously substituted sines and cosines
	#
	f := subs( seq(A[k]=sublist[k], k=1..kmax ), f );

	if isfun=YES then
		f := op(0,g)(f);
	fi;

	f;

end:



#=============================================================================
#
# Returns the number of function calls required by expr.
#
#=============================================================================
readlib(cost);
FunctionCalls := proc( expr )
	local t, p;
	t := cost(expr);
	if not has(t,`functions`) then
		RETURN(0);
	fi;
	if type( t, `+` ) then
		for p in t do
			if has( p, `functions` ) then
				if is( op(1,p), integer ) then	#e.g. p = 5 `functions`
					RETURN( op(1,p) );
				else							#p = `functions`
					RETURN(1);
				fi;
			fi;
		od;
	elif type( t, `*` ) then
		RETURN( op(1,t) );
	elif t=`functions` then
		RETURN(1);
	else
		RETURN(0);
	fi;
end:



#=============================================================================
#
# aliased to CombTrigOpt
#
# Use CombineTrig() and FunctionCalls() to find the combination of
# angles, two at a time, that produces the minimum number of function
# calls and return that most-efficient form of expr.
#
#=============================================================================
`CombineTrig/optimize` := proc( expr )

	local p, j, k, kmax, vars, L1, L2, c0, c1, c2, f, f1, f2,
		  n, verbose;

	if nargs > 1 and (args[2]=ON or args[3]=ON) then
		verbose := ON;
	else
		verbose := OFF;
	fi;

	f := copy(expr);
	if type(f,'array'(1)) then
		RETURN( map( procname, f, args[2..nargs] ) );
	fi;
	if type(f,'array'(2)) then
		RETURN( map( procname, f, args[2..nargs] ) );
	fi;
	if type(f,`=`) then
		RETURN( map( procname, f, args[2..nargs] ) );
	fi;

	#
	# fiducial cost
	#
	c0 := FunctionCalls(f);
	if verbose=ON then
		printf(`fiducial cost = %d functions`, c0 );
		print();	#force the printf to the screen
	fi;
	if c0=1 then
		RETURN(eval(f));
	fi;

	#
	# pick out the angles (or at least the indeterminates)
	#
	if nargs=1 or (nargs=2 and not type(args[2],list)) or
	   (not type(args[2],list) and not type(args[3],list)) then
		k := 1;
		for p in indets(f) do
			if not hasfun(p,cos) and not hasfun(p,sin) and not isdiff(p) then
				vars[k] := p;
				k := k + 1;
			else		# include already-existing multiple angles
				if type( op(1,p), `+` ) then
					vars[k] := op(1,p);
					k := k + 1;
				fi;
			fi;
		od;
	elif type(args[2],list) then
		vars := copy(args[2]);
	elif nargs=3 and type(args[3],list) then
		vars := copy(args[3]);
	else
		ERROR(` wrong arg type(s)`);
	fi;
	vars := convert(vars,set);	# weed out any repeats
	kmax := nops(vars);
	if verbose=ON then
		printf(`Trying %d combinations...`, kmax*(kmax-1)/2 );
		print();
	fi;

	#
	# loop over angle combinations
	#
	L1 := CombineTrig( f, vars[1], vars[2] );
	c1 := FunctionCalls( L1 );
	if c1=0 then
		RETURN( eval(L1) );
	fi;
	if verbose=ON then
		printf( `  1: cost[%a][%a] = %d functions`, vars[1], vars[2], c1 );
		lprint(` `);
	fi;
	n := 2;
	for j from 1 to kmax do
		for k from j+1 to kmax do
			if k <> j and not (j=1 and k=2) then
				L2 := CombineTrig( f, vars[j], vars[k] );
				c2 := FunctionCalls( L2 );
				if verbose=ON then
					printf( `%3d: cost[%a][%a] = %d functions`,
							 n, vars[j], vars[k], c2 );
					print();
					n := n + 1;
				fi;
				if c2 < c1 then
					L1 := L2;
					c1 := c2;
				fi;
			fi;
		od;
	od;

	if c0 <= c1 then
		if verbose=ON then
			printf( `Best form is original: cost = %d functions\n`, c0 );
		fi;
		RETURN(eval(f));  #no more efficient form than original
	else
		if verbose=ON then
			printf( `Best cost = %d functions\n`, c1 );
		fi;
		RETURN(eval(L1));
	fi;

end:
alias(CombTrigOpt=`CombineTrig/optimize`);



#==============================================================================
#
# Calculate the average of expr over a circular annulus of
# inner radius a and outer radius b.  The coordinates used
# must be contained in the list coords, radius first then
# azimuth -- e.g., coords := [rho,phi].  The list paramlist
# is a list of ordering variables for collect.
#
#                        b  2*Pi
#                        /\ /\
#                1       |  |
# <expr> := ------------ |  | r * expr(r,phi) dphi dr
#           Pi*(b^2-a^2) |  |
#                       \/ \/
#                       r=a phi=0
#
#==============================================================================
annular_average := proc( expr, coords, a, b, paramlist )

    local I_phi, I_rho, I1, I2, expr_avg;

    debug_print( procname, `Calculating average over annulus...`, 1 );
    I_phi := int( expr, coords[2] );
    I_phi := subs( coords[2]=2*Pi, I_phi ) - subs( coords[2]=0, I_phi );
    I_rho := int( I_phi*coords[1], coords[1] );
    I1    := subs( coords[1]=b, I_rho );
    I2    := subs( coords[1]=a, I_rho );
    expr_avg := 1/Pi/(b^2-a^2)*(I1-I2);
    expr_avg := collect( expr_avg, paramlist, 'recursive', simplify );

    RETURN( expr_avg );
end:



#=============================================================================
#
# Determine if any of the subexpressions of expr are explicitly
# multiplied by I=sqrt(-1).
#
# Returns true if I does not explicitly appear anywhere in expr,
# false if it does.  WARNING: indexed quantities are treated as real,
# since Maple does not allow indexed names to be assumed.
#
#=============================================================================
is_real_parts := proc( expr )

    local p, result;
    option remember;

    for p in expr do
        if p = (-1)^(1/2) then
            RETURN(false)
        elif nops(p) > 1 then
            result := procname(p)
        else
            result := is(p,real)
        fi;
        if result = false then
            RETURN(false)
        elif result = FAIL then
			if is(p,indexed) then
				RETURN(true);
			fi;
            debug_print( procname,
						 `Unable to determine if `.(convert(p,string)).
						 ` is real or not!`, 0 );
            ERROR();
        fi
    od;

    RETURN(true);
end:



#=============================================================================
#
# Return a 3x3 matrix which will rotate a vector ccw around an axis.
# Positive rotation is ccw when viewing from the positive axis direction
# toward the negative (right-hand rule).
#
#=============================================================================
rotaxis := proc( theta::algebraic, axis::string )
	local R;
	if axis=`x` or axis=`X` then
		R := array( 1..3, 1..3, [
					[ 1,          0,           0 ],
					[ 0, cos(theta), -sin(theta) ],
					[ 0, sin(theta),  cos(theta) ]
				   ] );
	 elif axis=`y` or axis=`Y` then
		R := array( 1..3, 1..3, [
					[  cos(theta), 0, sin(theta) ],
					[           0, 1,          0 ],
					[ -sin(theta), 0, cos(theta) ]
				   ] );
	 elif axis=`z` or axis=`Z` then
		R := array( 1..3, 1..3, [
					[ cos(theta), -sin(theta), 0 ],
					[ sin(theta),  cos(theta), 0 ],
					[          0,           0, 1 ]
				   ] );
	 else
		 ERROR(`:  axis must be "x", "y", or "z"`);
	 fi;
	 RETURN( op(R) );
end:



#=============================================================================
#
# Rotate a 3D vector by angle phi around the axis defined by w.  Positive
# rotation is defined by the right-hand rule (clockwise when viewing
# along w).
#
# Neither w nor vec needs to be a unit vector.
#
#=============================================================================
rotvec := proc( v::{vector,matrix}, w::{vector,matrix}, phi::algebraic )

	local   newvec, w_hat, V, W;

	if type(v,'array'(1..3,1..1)) then
		V := convert(eval(v),vector);
	else
		V := v;
	fi;
	if type(w,'array'(1..3,1..1)) then
		W := convert(eval(w),vector);
	else
		W := w;
	fi;
	if nops(eval(V))<>3 or nops(eval(W))<>3 then
		ERROR(`: input vectors must contain 3 elements only`);
	fi;

	w_hat  := evalm(W/mag(W));
	newvec := eval(V)*cos(phi) + (1-cos(phi))*(dot(V,w_hat))*eval(w_hat)
			  - cross(V,w_hat)*sin(phi);

	if type(v,matrix) then
		eval(mat(evalm(newvec)));
	else
		eval(newvec);
	fi;
end:



#=============================================================================
#
# Apply func to the specified operands (as a group) of the expression.
#
#=============================================================================
funcops := proc( expr::`+`, func::procedure, oplist::{integer,list} )
	local part1, part2, i, k;
	part1 := 0;
	for i from 1 to nops(oplist) do
		if type(oplist,integer) then
			k := oplist;
		else
			k := oplist[i];
		fi;
		part1 := part1 + op(k,expr);
	od;
	part2 := expr - part1;
	part1 := func( part1, args[4..nargs] );
	RETURN( part2 + part1 );
end:



#=============================================================================
#
# Complete the square in the (optional) variable var in the (2nd degree in
# var) polynomial expr.
#
# var can be a name, a function, a `+`, `*`, or `^` expression, or a
# set or list of functions, variable names, and `+`, `*`, `^` expressions.
#
# An additional argument -- a Maple procedure name -- may be included and,
# if present, will be applied to the remainder terms.  For example,
#
#    csquare( expr, [x,y], factor );
#
#=============================================================================
csquare := proc( expr::{algebraic,array,indexed,equation,numeric} )

    local	A, B, C, newexpr, var, t, i, L, R,
			subsflag, savevar, srvars, F, func;
	global	_H_;

    if type(expr,'numeric') then
        RETURN( expr );
    fi;

    #
    # grab the indeterminate if no second arg; check for function to apply
    #
	func := NULL;
    if nargs=1 then
        if nops(indets(expr))=0 then
            RETURN(expr);
        elif nops(indets(expr))=1 then
            var := op(indets(expr));
        else
            ERROR( `: Can't determine the indeterminate!` );
        fi;
    elif nargs=2 then
		if type(args[2],procedure) then
			func := args[2];
		else
			var := args[2];
		fi;
    elif nargs=3 then
		if type(args[2],procedure) then
			func := args[2];
			var  := args[3];
		elif type(args[3],procedure) then
			func := args[3];
			var  := args[2];
		else
			ERROR(`: wrong argument type(s)`);
		fi;
    else
        ERROR(`: wrong number arguments!`);
    fi;
	if not type(var,{name,set,list,`+`,`*`,`^`,function}) then
		ERROR( `: var must be one of `
			   .`{name,set,list,``+``,``*``,``^``,function}.` );
	fi;

	newexpr := copy(expr);
	if type(var,{set,list}) then
		for t in var do
			newexpr := procname( newexpr, t, func );
		od;
		RETURN( eval(newexpr) );
	fi;

    #
    # complete the square
    #
	_H_      := '_H_';
	subsflag := OFF;					#used to flag algebraic var
	if type( newexpr, {array,indexed,equation} ) then
		if type(var,{`+`,`*`,`^`,function}) then
			subsflag := ON;
			savevar  := var;
			var      := _H_;
			if type(savevar,function) then
				srvars := [savevar];
			else
				srvars := [op(indets(savevar))];
			fi;
			#
			# Define a mapping function to apply to each element of newexpr.
			# F does side relation simplification so as to substitute the
			# dummy variable _H_ for whatever expression is contained in var.
			#
			F := (_arg,_v,_srv)->
				if has(_arg,_v) then
					simplify( _arg, {_v=_H_}, _srv );
				else
					_arg;
				fi;
			newexpr := map( F, newexpr, savevar, srvars );
		fi;
		newexpr := map( procname, newexpr, var, func );
		if subsflag=ON then
			newexpr := subs( _H_=savevar, eval(newexpr) );
		fi;
	elif type( newexpr, polynom(anything,var) ) then
		if type(var,{`+`,`*`,`^`,function}) then
			subsflag := ON;
			savevar  := var;
			var      := _H_;
			newexpr  := simplify( newexpr, {savevar=_H_}, [op(indets(savevar))] );
		fi;
		newexpr := collect( newexpr, var );
		A := coeff( newexpr, var, 2 );	#form: A*var^2 + B*var + C
		B := coeff( newexpr, var, 1 );
		if A <> 0 and B <> 0 then
			C := coeff( newexpr, var, 0 );
			if func <> NULL and type(func,procedure) then
				newexpr := A*(var + B/(2*A))^2 + func(C-(B^2/(4*A)) );
			else
				newexpr := A*(var + B/(2*A))^2 + (C - (B^2/(4*A)) );
			fi;
		fi;
		if subsflag=ON then
			newexpr := subs( _H_=savevar, newexpr );
		fi;
	else
        RETURN( eval(expr) );
	fi;

	RETURN( eval(newexpr) );
end:



#=============================================================================
#
# Apply a function to an expression n times.
#
#=============================================================================
applyfunc := proc( expr, func::{procedure,function}, n::nonnegint )
	if n=0 then
		expr;
	else
		if nargs > 3 then
			func( procname(expr,func,n-1,args[4..nargs]), args[4..nargs] );
		else
			func( procname(expr,func,n-1) );
		fi;
	fi;
end:



#=============================================================================
#
# Apply func to individual (top-level) terms in expression expr.
#
#=============================================================================
termfunc := proc( expr::{algebraic,array(1),`=`}, func::procedure )

    local   tlist, N, i, newexpr, tmp;

    newexpr := copy( expr );

    if type(newexpr,{'array'(1),`=`}) then
		RETURN( map( procname, newexpr, func, args[3..nargs] ) );
    fi;

    if type(newexpr,{`+`,`*`}) then
		tlist := map( func, [op(newexpr)], args[3..nargs] );
        if type(newexpr,`*`) then
            newexpr := convert( tlist, `*` );
        else
            newexpr := convert( tlist, `+` );
        fi;
    elif type(newexpr,`^`) then
        N := nops(newexpr);
        for i from 1 to N do
            newexpr := subsop( [i]=procname(op(i,newexpr),func,args[3..nargs]),
							   newexpr );
            if newexpr=1 then break fi;
        od;
    else
        newexpr := func( newexpr, args[3..nargs] );
    fi;

    gc();
    RETURN( newexpr );
end:



#=============================================================================
#
# Selectively factor out a subexpression subexpr from an expression expr.
#
#=============================================================================
pullout := proc( expr::{name,`+`,`^`,`*`,`=`,function},
				 subexpr::{name,function,`+`,`*`,`^`,list} )

	local newexpr1, newexpr2, p, q, flag, r, allpunts;

	if type(expr,name) then
		RETURN(expr);
	fi;

	#
	# special case: expr is an equation
	#
	if type(expr,`=`) then
		RETURN( map(procname,expr,subexpr) );
	fi;

	#
	# special case: we have a list of subexpressions -- do each in order
	#
	if type(subexpr,list) then
		newexpr1 := expr;
		for p in subexpr do
			if not type(p,{name,function,`+`,`^`,`*`}) then
				ERROR(`: subexpression `.(convert(p,string)).` not of type name, function, +, ^, or *`);
			fi;
			newexpr1 := procname( newexpr1, p );
		od;
		RETURN(newexpr1);
	fi;

	#
	# special case: expr is type function
	#
	if type(expr,function) then
		RETURN(expr);
	fi;

	#
	# special case: expr is type `*`
	#
	if type(expr,`*`) then
		newexpr1 := expr;
		#
		# make sure expr contains the entire subexpr
		#
		flag := ON;
		if type(subexpr,{`+`,`^`}) then
			if not has(expr,subexpr) then
				flag := OFF;
			fi;
		else
			for q in subexpr do
				if not has(expr,q) then
					flag := OFF;
					break;
				fi;
			od;
		fi;
		if flag=ON then
			#
			# factor the first term of expr that contains subexpr
			#
			for p in expr do
				if type(p,{`+`,`^`}) then
					newexpr1 := expr/p * procname(p,subexpr);
					break;
				fi;
			od;
		fi;
		RETURN(newexpr1);
	fi;

	#
	# special case: expr is type `^`
	#
	if type(expr,`^`) then
		newexpr1 := expr;
		if type(op(1,expr),`+`) then
			newexpr1 := procname(op(1,expr),subexpr)^op(2,expr);
		fi;
		RETURN(newexpr1);
	fi;
	if not type(expr,`+`) then
		RETURN(expr);
	fi;

	#
	# general case: expr is type `+`
	#
	# first, some loop setup
	#
	newexpr1 := 0;
	newexpr2 := 0;
	allpunts := YES;
	#
	# loop over the parts of expr
	#
	for p in expr do
		#
		# first make sure each part contains the entire subexpr
		#
		flag := ON;
		if type(subexpr,{name,function,`+`,`^`}) then
			if not has(p,subexpr) then
				flag := OFF;
			fi;
		else	# is a `*`
			for q in subexpr do
				if not has(p,q) then
					flag := OFF;
					break;
				fi;
			od;
		fi;
		#
		# now pull out the subexpr from those terms of
		# expr that contain it
		#
		if flag=ON then
			if type(p,{`*`,`^`}) then
				r := procname(p,subexpr);
				if factor(r-p)=0 then
					newexpr1 := newexpr1 + r/subexpr;
				else
					newexpr1 := newexpr1 + r;
				fi;
				allpunts := NO;
			#
			# no dice -- punt on this p
			#
			else
				newexpr1 := newexpr1 + (p/subexpr);
			fi;
		#
		# punt on this p since it doesn't contain subexpr
		#
		else
			newexpr2 := newexpr2 + p;
		fi;
	od;

	if allpunts=YES then
		expr;
	else
		newexpr2 + subexpr*newexpr1;
	fi;

end:



#=============================================================================
#
# Test expr for the presence of a restricted form of polynomial in
# any of the variables in varlist.  Returns 'true' or 'false'.
#
# The restriction in the definition of "polynomial" is that the
# expression or any of the parts of the expression must contain 2
# or more terms separated by an addition operator ('+' or '-'), at
# least one of which contains a power of at least one of the
# variables in varlist.  For example, suppose varlist:=X.  Then
#
#      2                  3
#   a*X    and   (a + b) X
#
# are NOT "polynomials" as far as polytest is concerned, but
#
#      2                        3
#   a*X  + b   and   (a + b X) X
#
# are.  polytest() is designed to be used by subxfunc(), which
# makes use of this restriction.
#
#=============================================================================
polytest := proc( expr::algebraic, varlist::{list,string} )

    local p;

    if type( expr, `+` ) then
        for p in varlist do
            if has(expr,p) and degree(expr,p) <> 0 then
                RETURN( true );
            fi;
        od;
    elif type( expr, {`*`,`^`} ) then
        for p in expr do
            if polytest(p,varlist) = true then
                RETURN( true );
            fi;
        od;
    fi;

    false;

end:



#=============================================================================
#
# Check an expression or a list of expressions for small divisors in
# the variables given in varlist.
#
# For example, if varlist:=[x,b] and A:=(a+x)/(x*c+b^2), then
#
#   small_divisors( A, varlist ) = true
#
# But if varlist:=[x,c], then
#
#   small_divisors( A, varlist ) = false
#
# The argument varlist is a list (e.g., [theta,f,C]) or a vector.
#
# small_divisors() makes use of common_factor().
#
#=============================================================================
small_divisors := proc( expr::{algebraic,list(algebraic),'array'(1),`=`},
						varlist::{list,string} )

    local den, num, count, _H_, p, equ;

    assume( _H_, real );

    #if expr is a list of equations, do each equation
    if not type(expr,algebraic) then
        for equ in expr do
            if procname(equ,varlist) then
                RETURN(true);
            fi;
        od;
        RETURN(false);
    fi;

    count := 0;
    for p in varlist do
        if not has( expr, p ) then
            count := count + 1;
        fi;
    od;
    if count = nops(varlist) then       #no instances of varlist variables
        RETURN(false);
    fi;

    equ := normal( expr );

    #multiply all varlist variables by a dummy variable _H_ to catch cross terms
    for p in varlist do
        equ := subs( p=_H_*p, equ );
    od;

    for p in equ do
        den := denom( p );
        num := numer( p );
        if den <> 1 then
            if common_factor( den, _H_ ) = true then
                RETURN(true);
            fi;
        fi;
        if denom(num) <> 1 then
            if procname( num, _H_ ) = true then
                RETURN(true);
            fi;
        fi;
    od;

    false;

end:



#=============================================================================
#
# Utility function for small_divisors().
#
# Check an expression for a common multiplicative factor in the variables
# contained in varlist.  It is assumed that negative powers in the list
# variables have been removed (e.g., by feeding common_factor the numerator
# or denominator of an equation in reduced normal form).
#
# For example, if varlist:=x, then
#
#   common_factor( a*(1+x^2+x^3)+b*x, varlist ) = false
#
# But if varlist:=[x,a], then
#
#   common_factor( a*(1+x^2+x^3)+b*x, varlist ) = true
#
# The argument varlist is a Maple list (e.g., [theta,f,C]).
#
#=============================================================================
common_factor := proc( expr::algebraic, varlist::{list,string} )

    local equ, _H_, p;
    assume( _H_, real );

    equ := copy( expr );

    #multiply all variables in varlist by dummy variable _H_ to catch
	#cross terms
    for p in varlist do
        equ := subs( p=_H_*p, equ );
    od;

    if ldegree( collect(equ,_H_), _H_ ) > 0 then
        RETURN(true);
    fi;

    RETURN(false);
end:



#=============================================================================
#
# Expand 'expression' into a series up to and including terms of order
# 'expansion_order' in the small parameters 'params'.
#
# If _order < 1, then no expansions are performed.
#
# params is a Maple list (e.g., [theta,f,C]).
#
# Multiple independent expansions can be done in one step.  For example,
#
#   expansion( expr, theta, 3, [epsilon,psi,alpha], 2, [phi,beta], 5 );
#
# NOTE: there is a Maple bug such that sometimes s is not characterized
# as a "series" type if there are noninteger powers of params
# involved in expression. I haven't been able to figure out a
# workaround for this.  I have tried convert(s,series,params),
# which fails.  When the bug appears, conversion of s to type
# polynom does not remove the O() term.
#
#=============================================================================
# Return Value
# ------------
# The series expansion of expression, converted to the
# Maple type polynom.
#=============================================================================
expansion := proc( expression::{list,set,'array'(1),'array'(2),matrix,`=`,algebraic},
				   params::{set,name,list,function}, expansion_order::integer )

    local   s, s1, ord1, _H_, p, plist, degree_s1, i;

	if type(nargs,even) then
		ERROR(`: wrong number of arguments`);
	fi;

	if not type( eval(expression), algebraic ) then
		RETURN( map( procname, args ) );
	fi;

	if not has(expression,params) then
		RETURN(expression);
	fi;

	#
	# Handle multiple independent expansions
	#
	if nargs > 3 then
		s := procname( expression, params, expansion_order );
		for i from 1 to (nargs-3)/2 do
			s := procname( s, args[2*i+2], args[2*i+3] );
		od;
		RETURN(eval(s));
	fi;


	if expansion_order < 1 then
		debug_print( procname, `WARNING!  Expansion order is less than 1`, 4 );
		debug_print( procname, `Expression NOT expanded:`, 5, expression );
		RETURN(expression);
	fi;

    assume( _H_, real );
    s    := copy(expression);
	ord1 := expansion_order + 1;

	#
	# Expand in a series in the small parameters, taking into account
	# cross terms by expanding on a dummy variable _H_.
	#
	if not type(params,list) then
		plist := [params];
	else
		plist := params;
	fi;
	for p in plist do
		s := subs( p=_H_*p, s );
	od;

	convert( series(s,_H_,ord1), polynom );
	s1 := collect( ", [_H_,op(params)] );

	#
	# Version V release 4 of Maple does not expand the expression
	# in a series with the intended *output* order.  (It only
	# guarantees *internal* calculations to the correct order.)
	# Hence sometimes we must keep increasing the order until the
	# degree of the output polynomial is correct.
	#
	degree_s1 := degree(s1,_H_);
	if degree_s1 = FAIL then
		debug_print( procname, `WARNING!  Unable to determine degree `
					 .`of expanded expression.`, 2 );
		debug_print( procname, `Problem expression:`, 2, s1 );
	else
		#
		# The FAIL test is so that expressions whose degree in {params} is
		# already less than _order won't trigger the while loop.
		#
		if degree(s,convert(params,set))=FAIL and nops(s1) > 1
		and degree_s1 < expansion_order then
			while degree_s1 < expansion_order do
				ord1 := ord1 + 1;
				debug_print( procname, `WARNING!  Increasing expansion order to `
							 .(convert(ord1-1,string)).`...`, 4 );
				convert( series(s,_H_,ord1), polynom );
				s1 := collect( ", [_H_,op(params)] );
				degree_s1 := degree(s1,_H_);
				if degree_s1 = FAIL then
					debug_print( procname, `WARNING!  Unable to determine degree `
								 .`of expanded expression.`, 2 );
					debug_print( procname, `Problem expression:`, 2, s1 );
					break;
				fi;
				if ord1 > expansion_order+7 then
					debug_print(procname,`WARNING! Too many order increases!`,4,ord1);
					debug_print( procname, `Problem expression:`, 5, s1 );
					break;
				fi;
			od;
		fi;

		#
		# For some expressions (like stuff under square roots), Maple's
		# series() expands beyond the power asked for by 2 or more orders.
		# This seems to be a bug.  Redoing the expansion a second time
		# produces the correct expansion.
		#
		if degree_s1 > expansion_order then
			debug_print( procname, `WARNING!  Encountered series() bug...`, 4 );
			debug_print( procname, `Buggered expression:`, 5, s1 );
			s1 := convert( series(s,_H_,ord1), polynom );
		fi;
	fi;

	#
	# Maple gives a bogus O(1) for the series expansion
	# under certain obscure conditions.  A workaround is to force
	# it to find the leading term; then something about the
	# execution history that Maple remembers will trigger the
	# correct series expansion.  It's a time waster, but it's
	# the only way I can think of for getting around the problem.
	#
	if s1=0 then
		debug_print( procname, `WARNING!  Fixing O(1) bug in series()...`, 4 );
		series( leadterm(s), _H_, 999 );
		s1 := convert( series(s,_H_,ord1), polynom );
	fi;

	subs( _H_=1, s1 );
	s := collect( ", params );

    gc();
    RETURN( s );
end:


fi;