#=============================================================================
#
# Some Maple procedures useful, among other things, for AESOP.
#
# 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 AESOPUTL_READ <> 1 then
AESOPUTL_READ := 1:
read`c:/maple/utils.p`:
read`d:/optics/AESOP/objects.p`:

#=============================================================================
#
# transform_dir( invec, global_dir, xform_direction )
#
# Transform the input direction vector to global or local coordinates.
#
#=============================================================================
#
# Arguments
# ---------
#
# The input direction angles are transformed according to xform_direction:
#
#	   xform_direction := 1;	   from local to global coordinates
#	   xform_direction := 0;	   from global to local coordinates
#
# The input direction is of type 'vector'.
#
# global_dir (type 'vector') contains the components of the unit direction
# vector of the z axis of the local frame expressed in the global frame.
#
#=============================================================================
#
# Return Value
# ------------
#
# The transformed direction vector as an AESOP 'vector' data type.
#
#=============================================================================
transform_dir := proc( invec::vector, global_dir::vector,
					   xform_direction::integer )

	local   g, out, h, rot, vglob, vloc:

	assume( h >= 0 ):

	g   := copy( global_dir ):
	out := copy( invec ):
	h   := sqrt( g[1]^2 + g[2]^2 ):

	if xform_direction = LOCAL then	 #global to local

		if h = 0 then				   #axial -- no rotation
			rot := array( [ [1, 0, 0],
							[0, 1, 0],
							[0, 0, 1] ] ):
		else
			rot := array( [ [g[1]*g[3]/h, g[2]*g[3]/h,   -h],
								[-g[2]/h,	  g[1]/h,	0],
								   [g[1],		g[2], g[3]] ] ):
		fi;

		vloc := evalm( rot &* invec ):
#	   out  := map( normal, vloc ):
		out  := vloc:

	elif xform_direction = GLOBAL then  #local to global

		if h = 0 then				   #axial -- no rotation
			rot := array( [ [1, 0, 0],
							[0, 1, 0],
							[0, 0, 1] ] ):
		else
			rot := array( [ [g[1]*g[3]/h, -g[2]/h, g[1]],
							[g[2]*g[3]/h,  g[1]/h, g[2]],
							[-h, 0, g[3]] ] ):
		fi;

		vglob := evalm( rot &* invec ):
#	   out   := map( normal, vglob ):
		out   := vglob:

	else									#direction specification error

		ERROR(`Invalid transform direction.`):

	fi:

	gc():
	RETURN( op(out) ):
end:



#=============================================================================
#
# transform_pos( inpos, global_pos, global_dir, xform_direction )
#
# Transform the input position to global or local coordinates.
#
#=============================================================================
#
# Arguments
# ---------
#
# The input position is transformed according to xform_direction:
#
#	   xform_direction := 1;	   from local to global coordinates
#	   xform_direction := 0;	   from global to local coordinates
#
# The input position is of type 'point'.
#
# global_dir (type 'vector') contains the components of the unit
# direction vector of the z axis of the local frame expressed in
# the global frame.
#
# global_pos (type 'point') contains the position in the global frame of the
# origin of the local frame.
#
#=============================================================================
#
# Return Value
# ------------
#
# The transformed position, of type 'point'.
#
#=============================================================================
transform_pos := proc( inpos::point, global_pos::point, global_dir::vector, xform_direction::integer )

	local   out, v:

	if xform_direction = LOCAL then		 #global to local

		v   := evalm( inpos - global_pos ):
		v   := transform_dir( v, global_dir, xform_direction ):
		out := point( v ):

	elif xform_direction = GLOBAL then	  #local to global

		v   := convert(inpos,vector);
		v   := transform_dir( v, global_dir, xform_direction ):
		out := evalm( global_pos + v ):

	else									#direction specification error

		ERROR(`Invalid transform direction in transfor_pos`):

	fi:

	RETURN( op(out) ):
end:



#=============================================================================
#
# simpx( expr, siderels, varlist )
#
# Replacement for Maple's simplify/siderel().
# Knows about AESOP 'point', 'vector', and 'beam' data types.
#
#=============================================================================
simpx := proc( expr::{algebraic,point,vector,array(1),beam},
			   siderels::set, varlist::{list,string} )

	local   newexpr:

	#AESOP point or vector
	newexpr := copy( expr ):
	if type(newexpr,'beam') then
		newexpr[pos]  := procname( newexpr[pos],  siderels, varlist ):
		newexpr[dir]  := procname( newexpr[dir],  siderels, varlist ):
		newexpr[path] := procname( newexpr[path], siderels, varlist ):
		RETURN( op(newexpr) ):
	fi:
	if type(newexpr,{'point','vector','array'(1)}) then
		newexpr := map( procname, newexpr, siderels, varlist );
		RETURN( op(newexpr) ):
	fi:

	newexpr := simplify( newexpr, siderels, varlist ):

	gc():
	RETURN( newexpr ):
end:



#==============================================================================
#
# simpside( expr, param_list, sublist, siderels, siderels_list )
#
# Handles AESOP beam, point, or vector objects.
#
# Perform "extended" side relation simplification on expr.  This routine
# provides you with a way to "undo" the side relations, then simplify via
# subxfunc(simplify...), then simplify with side relations.
#
# param_list = parameter list for use in collect( expr, param_list, simplify )
# sublist	 = list of substitutions to make before side relation simplification
# siderels   = list of side relations
# siderels_list = side relation variables ordering list:
#				  simplify( expr, siderels, siderels_list )
#
#==============================================================================
simpside := proc( expr, param_list, sublist, siderels, siderels_list )

	local newexpr:

	newexpr := copy( expr );

	if type(newexpr,'beam') then
		newexpr[pos]  := procname( newexpr[pos], param_list, sublist,
								   siderels, siderels_list );
		newexpr[dir]  := procname( newexpr[dir], param_list, sublist,
								   siderels, siderels_list );
		newexpr[path] := procname( newexpr[path], param_list, sublist,
								   siderels, siderels_list );
		RETURN( op(newexpr) ):
	elif type(newexpr,{'point','vector'}) then
		newexpr := map( procname, newexpr, param_list, sublist,
						siderels, siderels_list );
		RETURN( op(newexpr) ):
	else
		_simpside( newexpr, param_list, sublist, siderels, siderels_list );
	fi:

end:

_simpside := proc( expr, param_list, sublist, siderels, siderels_list )
	local p, q, loc;
	subx( expr, sublist );
	p := collect( ", param_list, simplify );
	if type(p,`+`) then
		for q in p do
			loc := location(p,q);
			if loc <> [] then
				q := simpx( q, siderels, siderels_list );
				p := subsop( loc=q, p );
			fi;
		od;
	else
		p := simpx( p, siderels, siderels_list );
	fi;
	collect( p, param_list, simplify );
end:


#=============================================================================
#
# subx( expr, exprlist )
#
# Replacement for Maple's subs().
# Knows about AESOP 'point', 'vector', and 'beam' data types.
#
#=============================================================================
subx := proc( expr::{algebraic,point,vector,beam},
			  exprlist::{list,set,string,`=`} )

	local   newexpr:

	newexpr := copy( expr ):
	if type(newexpr,'beam') then
		newexpr[pos]  := procname( newexpr[pos], exprlist ):
		newexpr[dir]  := procname( newexpr[dir], exprlist ):
		newexpr[path] := procname( newexpr[path], exprlist ):
		RETURN( op(newexpr) ):
	fi:
	if type(newexpr,{'point','vector'}) then
		newexpr := map( procname, newexpr, exprlist );
		RETURN( op(newexpr) ):
	fi:

	newexpr := subs( exprlist, newexpr ):

	gc():
	RETURN( newexpr ):
end:



#=============================================================================
#
# RotRay( position, direction, p, axis, phi, simplflag )
#
# Rotate a position vector "position" and a direction vector "direction"
# by "phi" around the pivot point "p" and direction vector "axis".  Positive
# rotation is defined by the right-hand rule (clockwise when viewing
# along direction).  If simplflag=ON, partial simplification of
# the answer is performed.
#
#=============================================================================
RotRay := proc( position::point, direction::vector, p::point,
				axis::vector, phi::algebraic, simplflag::integer )

	local   newvec, pivot_vec:

	#
	# rotate a vector from the pivot point to the anchor of the ray
	#
	newvec[pos] := copy( position );
	if position[1] <> p[1]  and  position[2] <> p[2]  and
	   position[3] <> p[3]  then
		pivot_vec   := evalm( position - p );
		newvec[pos] := evalm( rotvec( pivot_vec, axis, phi ) );
		newvec[pos] := evalm( p + newvec[pos] );
	fi:

	#
	# rotate the direction vector
	#
	newvec[dir] := evalm( rotvec( direction, axis, phi ) );

	gc():
	RETURN( op(newvec) ):
end:


fi;