#=============================================================================
#
# reflect( inputbeam, surface, n, simplepath, params, order, param_ordering,
#           nodenom_list )
#
# Find the reflection of a beam from a quadratic surface, as well as the
# cumulative optical path, up to and including order order in the small
# parameters params.
#
# If order < 0, then no expansions are performed.
#
#=============================================================================
#
# 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/
#
#=============================================================================
#
# Arguments
# ---------
#
# inputbeam
# ---------
# The argument inputbeam must be a table of data type 'beam':
#
#       beam[pos]     anchor position of beam (type 'point')
#       beam[dir]     direction of propogation (type 'vector')
#       beam[path]    accumulated optical path
#       beam[type] = 'beam'
#
# surface
# -------
# The argument surface is indexed as follows:
#
#       surface[eqn] := f(x,y,z);
#       surface[coord]
#       surface[pos]
#       surface[dir]    direction of surface symmetry axis
#       surface[type]   'mirror', 'beamsplitter', etc.
#
# where f(x,y,z) is the equation describing the surface
# in the local coordinate system.
#
# surface[coord] are the labels used for the (local) x, y, and z coordinates
# in the equation describing the surface, surface[pos] is the position of
# the surface vertex in the global coordinate system, and surface[dir]
# is the direction (unit vector) of the positive z axis of symmetry in the
# global coordinate system.
#
# n
# -
# The argument n is the desired intersection number.  In most cases
# where the surfaces are conicoids, there will be two possible intersections
# of a beam with the surface.  n is which intersection to use (1 or 2).
#
# simplepath
# ----------
# If simplepath=ON, simplification of the optical path length is
# performed (beware -- this can take an inordinate amount of time!).
#
# params
# ------
# params is a Maple list (e.g., [theta,f,C]) of small parameters to
# expand the expressions in at the intersection point.
#
# param_ordering
# -------
# param_ordering is a list of "important" variables by which the expressions may be
# organized and then recursively subexpression simplified via collect()
# in beam_intersect().
#
# nodenom_list
# ------------
# a list of variables which are not allowed in an
# expression denominator when finding beam intersection
# points.
#=============================================================================
#
# Return Value
# ------------
#
# reflect() returns a table containing the reflection of the incident
# beam with the surface, indexed in the same manner as the incident beam,
# and with respect to the global coordinate frame.
#
#=============================================================================
#
# Additional Requirements
# -----------------------
#
# The global variables
#       GLOBAL := 1;
#       LOCAL  := 0;
# must be defined before calling this procedure (used for transforming
# coordinate systems -- see transfor_dir() and transfor_pos() in utils.p).
#
# The floating point variable time0 must be initialized before calling
# this procedure:  time0:=time();
#
# The files intersct.p and utils.p must have been read by Maple before
# this routine can be used.
#
#=============================================================================
#
# Method
# ------
#
# See intersct.p for the method of finding the intersection point between
# the incident beam and the surface.
#
# In the local coordinate frame, the reflected beam is in the same plane
# as the incident beam and the normal to the surface at the intersection
# point.  Thus, we can write the equation of a line coincident with the
# reflected beam,
#
#       r'(t) = r + t*n
#
# where r is the incident beam, r' the reflected beam, n the normal
# to the surface at the intersection point, and t the parameterization
# of the line, which must be solved for.  Since the incident and reflected
# angles must be equal, we have
#
#       dot( r', n )     dot( r, n )
#       ------------ = - ------------
#         |r'| |n|         |r| |n|
#
# Combine this with
#                                         2
#       dot( r', n ) = dot( r, n ) + t*|n|
#
# and the further condition that the magnitude of the reflected beam
# is equal to the magnitude of the incident beam, and we can solve for t:
#
#              dot(r,n)
#       t = -2 --------
#                   2
#                |n|
#
# After r' is found, the result is transformed back to the global reference
# frame.
#
#=============================================================================

if REFLECT_READ <> 1 then
REFLECT_READ := 1:
read`c:/maple/startup.p`:
read`d:/optics/AESOP/intersct.p`:
read`d:/optics/AESOP/surfnorm.p`:

reflect := proc( inputbeam::beam, surface::algebraic, expansion_list,
				 param_ordering::{list,name}, nodenom_list::{list,name},
				 sublist::{list,name,set}, siderels::set,
				 siderels_list::{list,name} )

	local   reflbeam, reflvec, normvec, tmp, localsurf, ipt,
			x0, y0, z0, r0, v, vx, vy, vp, vz, st, ct:
	global	NEW_T;

	assume( tmp, real );

	reflbeam  := copy( inputbeam );
	localsurf := copy( surface );


	#
	# transform to local coordinate frame
	#
	# DO NOT perform series expansions here -- can lead to slight errors
	#
	debug_print( procname, `Transforming input beam to LOCAL frame...`, 1 );
	reflbeam[pos] := transform_pos( reflbeam[pos], localsurf[pos],
									localsurf[dir], LOCAL );
	reflbeam[dir] := transform_dir( reflbeam[dir], localsurf[dir], LOCAL );
	debug_print( procname, `input beam in local frame`, 4, eval(reflbeam) );


	#
	# find the intersection point (local coord's are returned)
	#
	debug_print( procname, `Finding intersection point...`, 1 );
	if not assigned(NEW_T) or NEW_T=OFF then
		ipt := beam_intersect( reflbeam, localsurf, LOCAL, expansion_list,
							   param_ordering, nodenom_list,
							   sublist, siderels, siderels_list );
	else
		x0 := reflbeam[pos][1];
		y0 := reflbeam[pos][2];
		z0 := reflbeam[pos][3];
		r0 := sqrt(x0^2+y0^2);
		v  := mag(reflbeam[dir]);
		vx := reflbeam[dir][1]/v;
		vy := reflbeam[dir][1]/v;
		vp := sqrt(vx^2+vy^2);
		vz := sqrt(1-vp^2);
#		if localsurf[side] = CONCAVE then
			vz := -vz;
#		fi;
		st := ( x0*vy - y0*vx ) / ( r0*vp );
		ct := ( x0*vx + y0*vy ) / ( r0*vp );
		r0 := expansion( r0, op(expansion_list) );
		vp := expansion( vp, op(expansion_list) );
		vz := expansion( vz, op(expansion_list) );
		st := expansion( st, op(expansion_list) );
		ct := expansion( ct, op(expansion_list) );
		ipt[path] := ConicoidSolve( localsurf[flen], localsurf[cc], r0, z0,
									st, ct, vp, vz, CONCAVE );
		ipt[path] := expansion( ipt[path], op(expansion_list) );
		ipt[path] := simpside( ipt[path], param_ordering, sublist,
								   siderels, siderels_list );
		ipt[pos]  := evalm( reflbeam[pos] +
							ipt[path]*reflbeam[dir]/mag(reflbeam[dir]) );
	fi;
	debug_print( procname, `Expanding intersection point...`, 1 );
	reflbeam[pos] := expansion( ipt[pos], op(expansion_list) );
	debug_print( procname, `Simplifying intersection point...`, 1 );
	reflbeam[pos] := simpside( reflbeam[pos], param_ordering, sublist,
							   siderels, siderels_list );
	debug_print( procname, `intersection point`, 4, reflbeam[pos] );

	#
	# find the normal vector to the concave side of the surface at
	# the intersection point
	#
	debug_print( procname, `Finding surface normal...`, 1 );
	normvec := surfnorm( localsurf, reflbeam[pos] );
	debug_print( procname, `normvec before expansion`, 5, normvec );

	debug_print( procname, `Expanding normal vector...`, 1 );
	normvec := expansion( normvec, op(expansion_list) );
	debug_print( procname, `Simplifying normal vector...`, 1 );
	normvec := simpside( normvec, param_ordering, sublist, siderels, siderels_list );
	debug_print( procname, `simplified normal vector`, 3, normvec );


	#
	# calculate the reflected beam direction
	#
	debug_print( procname, `Calculating reflected beam...`, 1 );
	tmp := 2*dot(reflbeam[dir],normvec) / dot(normvec,normvec);
	reflvec := evalm( reflbeam[dir] - tmp*normvec );

	debug_print( procname, `Expanding reflection vector...`, 1 );
	reflvec := expansion( reflvec, op(expansion_list) );
	debug_print( procname, `Simplifying reflection vector...`, 1 );
	reflvec := simpside( reflvec, param_ordering, sublist, siderels, siderels_list );
	debug_print( procname, `reflected beam in local frame`, 4, eval(reflvec) );


	#
	# transform to the global reference coordinate frame
	#
	debug_print( procname, `Transforming output beam to GLOBAL frame...`, 1 );
	debug_print( procname, `...direction...`, 1 );
	reflbeam[dir] := transform_dir( reflvec, localsurf[dir], GLOBAL );
	debug_print( procname, `...position...`, 1 );
	reflbeam[pos] := transform_pos( reflbeam[pos], localsurf[pos],
									localsurf[dir], GLOBAL );

	debug_print( procname, `Expanding position...`, 1 );
	reflbeam[pos] := expansion( reflbeam[pos], op(expansion_list) );
	debug_print( procname, `Simplifying position...`, 1 );
	reflbeam[pos] := simpside( reflbeam[pos], param_ordering, sublist, siderels, siderels_list );
	debug_print( procname, `Expanding direction...`, 1 );
	reflbeam[dir] := expansion( reflbeam[dir], op(expansion_list) );
	debug_print( procname, `Simplifying direction...`, 1 );
	reflbeam[dir] := simpside( reflbeam[dir], param_ordering, sublist, siderels, siderels_list );


	#
	# calculate the cumulative optical path
	#
	debug_print( procname, `Calculating optical path...`, 1 );
	reflbeam[path] := inputbeam[path] + ipt[path]:
	debug_print( procname, `Expanding optical path...`, 1 );
	reflbeam[path] := expansion( reflbeam[path], op(expansion_list) );
	debug_print( procname, `Simplifying optical path...`, 1 );
	reflbeam[path] := simpside( reflbeam[path], param_ordering, sublist,
								siderels, siderels_list );

	gc();
	RETURN( op(reflbeam) );
end:


fi;