#=============================================================================
#
# Maple procedure beam_intersect( inputbeam, surface, n, coordframe, small_parms,
#								 order, param_ordering, nodenom_list )
#
# Solve for beam-surface intersection point.
#
# Find the intersection point of a beam and a quadratic surface 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
# ---------
# 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'
#
# The argument surface must be indexed as follows:
#
#	   surface[eqn]   := f(xname,yname,zname);
#	   surface[coord] := [xname,yname,zname];
#	   surface[pos]
#	   surface[dir]
#
# where f(x,y,z) is the equation for the surface,
#
#	   f(x,y,z) = 0
#
# in the local coordinates, and xname, yname, and zname are the labels
# used for the 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 of the surface axis of symmetry in the global coordinate
# system.
#
# The argument coordframe signals whether the coordinate frame in which the
# input beam is expressed is the surface LOCAL frame or the GLOBAL
# reference frame.  If coordframe=LOCAL, no frame conversion is done.  If
# coordframe=GLOBAL, the input beam is transformed to the surface LOCAL
# frame, the intersection is found, then the intersection point is
# transformed to the GLOBAL frame before being returned.
#
# small_parms is a Maple list (e.g., [theta,f,C]) of small parameters to
# expand the expressions in at the intersection point.
#
# param_ordering is a list of "important" variables by which the expressions
# may be organized and then recursively subexpression-simplified via
# collect().
#
# nodenom_list
# ------------
# a list of variables which are not allowed in an
# expression denominator when finding beam intersection
# points.
#=============================================================================
#
# Return Value
# ------------
#
# intersct() returns a table containing the first solution for the
# intersection of the beam with the surface, along with the distance
# from the beam anchor to the intersection point, indexed as follows:
#
#	   return_value[pos]	   an AESOP 'point'
#	   return_value[path]	  the path length
#
# Coordinates are wrt the surface local reference frame if intersct()
# is called with coordframe=LOCAL, or wrt the global frame if
# coordframe=GLOBAL.
#
#=============================================================================
#
# Additional Requirements
# -----------------------
#
# The global variables
#	   GLOBAL := 1;
#	   LOCAL  := 0;
# must be defined before calling this procedure.
#
# The floating point variable time0 must be initialized before calling
# this procedure:  time0:=time();
#
# The file utils.p must have been read by Maple before this routine
# can be used.
#
#=============================================================================
#
# Method
# ------
#
# The input beam must be in the local coordinate frame of the surface.
#
# The parameterized equation for a line coincident with the input beam is
#
#	   r(t) = r0 + t*n
#
# where r0=(x0,y0,z0) is a point on the line (the beam anchor point), t is
# the parameterization of the position on the line (ie, the path length from
# r0), and n is the unit direction vector
#
#		   [ sin(psi)*cos(lambda) ]   [ nx ]
#	   n = [ sin(psi)*sin(lambda) ] = [ ny ]
#		   [	   cos(psi)	   ]   [ nz ]
#
# where (psi,lambda) are the polar and azimuthal angls wrt the +z axis
# and counterclockwise from the +x axis in the local coordinate frame.
#
# The equation of the surface (in the local frame) is of the form
#
#	   f(x,y,z) = 0
#
# We can find the value of t which corresponds to the intersection point
# by substituting
#
#	   x = x0 + t*sin(psi)*cos(lambda) = x0 + t*nx
#
#	   y = y0 + t*sin(psi)*sin(lambda) = y0 + t*ny
#
#	   z = z0 + t*cos(psi)			 = z0 + t*nz
#
# into f(x,y,z) and solving for t.  For a quadratic surface there will
# in general be two solutions.  Then we substitute the solution(s)
# for t back into the equations for the line,
#
#	   r(t) = r0 + t*n
#
# to determine the x, y, and z values of the intersection point.
#
# All calulations and result are in the local reference frame.
#
#=============================================================================

if INTERSCT_READ <> 1 then
INTERSCT_READ := 1:
read`c:/maple/startup.p`:
read`d:/optics/AESOP/opsys.p`:

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

	local  i, _t_, tmp, tmpt, eq, nsol, m, tval,
		   localsurf, localbeam, retval, bmag, p, loc, ord, cc_order;

	global Tval, Solns;

	assume( _t_, real );
	assume( tmpt, real );
	assume( tval, real );
	assume( bmag >= 0 );
	Solns := 'Solns';
	Tval  := 'Tval';

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

	#
	# transform to local frame if necessary
	#
	if coordframe=GLOBAL then
		localbeam[pos] := transform_pos( localbeam[pos], localsurf[pos],
										 localsurf[dir], LOCAL );
		localbeam[dir] := transform_dir( localbeam[dir], localsurf[dir],
										 LOCAL );
	fi:
	debug_print( procname, `input beam in local frame`, 4, eval(localbeam) );

	#
	# substitute x=x0+t*nx and y=y0+t*ny and z=z0+t*nz into the surface
	# equation f(x,y,z)=0
	#
	debug_print( procname, `Equation substitutions...`, 1 );
	tmp := evalm( localbeam[pos] + _t_*localbeam[dir] );
	tmp := map( normal, tmp );
	eq  := subs( localsurf[coord][1]=tmp[1],
				 localsurf[coord][2]=tmp[2],
				 localsurf[coord][3]=tmp[3],
				 localsurf[eqn] );
	debug_print( procname, `raw equation for t:`, 4, subs(_t_='t',eq) );

	#
	# do some simplifying (the expansion seems to be necessary for quadratic
	# surfaces; if not expanded, the sqrt in the soln for t doesn't expand
	# right, and bogus answers show up for small order expansions)
	#
	debug_print( procname, `Expanding substitution result...`, 1 );
	if type(eq,`+`) then
		for p in eq do
			loc := location(eq,p);
			if loc <> [] then
				p  := expansion( p, op(expansion_list) );
				eq := subsop( loc=p, eq );
			fi;
		od;
	fi;
	eq := expansion( eq, op(expansion_list) );
	debug_print( procname, `Simplifying substitution result...`, 1 );
	eq := simpside( eq, [_t_,op(param_ordering)], sublist, siderels, siderels_list );
	debug_print( procname, `simplified equation for t: `, 3, subs(_t_='t',eq) );

	#
	# solve for t
	#
	debug_print( procname, `Solving for t...`, 1 );
	if eq=0 then
		debug_print( procname, `WARNING: equation in t is zero!`, 1 );
		tval := 0;
	else
		#optical flat
		if localsurf[flen]=infinity then
			debug_print( procname, `Entering linear solver...`, 1 );
			tval := `beam_intersect/solve/linear`( eq, _t_, expansion_list,
												   param_ordering, nodenom_list );
		#paraboloid
		elif localsurf[cc]=0 then
			debug_print( procname, `Entering quadratic solver...`, 1 );
			tval := `beam_intersect/solve/quadratic`( eq, _t_, expansion_list,
													  param_ordering, nodenom_list );
		#spheroid
		elif localsurf[cc]=1 then
			debug_print( procname, `Entering spheroid solver...`, 1 );
			tval := `beam_intersect/solve/spheroid`( eq, _t_, expansion_list,
													 param_ordering, nodenom_list );
		#conicoid
		elif type(localsurf[cc],name) then
tval := `beam_intersect/solve/quadratic`( eq, _t_, expansion_list,
										  param_ordering, nodenom_list );
#			debug_print( procname, `Entering nonlinear solver...`, 1 );
#			loc := location(expansion_list,localsurf[cc]);
#			cc_order := expansion_list[loc[1]+1];
#			tval := `beam_intersect/solve/nonlinear`( eq, _t_, localsurf[cc],
#													  cc_order, expansion_list,
#													  param_ordering, nodenom_list );
		else
			ERROR(`: unknown surface type!`);
		fi;
		Tval := tval;	#assign to global for debugging
		debug_print( procname, `Simplified solution for t:`, 3, t=Tval );
	fi:

	#
	# assemble the intersection point
	#
	debug_print( procname,
				 `Adding vector to previous position and simplifying...`, 1 );
	tmp := evalm( localbeam[pos] + tval*localbeam[dir] );
	if expansion_list[2] > 0 then
		debug_print( procname, `Expanding in series...`, 1 );
		tmp := expansion( tmp, op(expansion_list) );
		tmp := Collect( eval(tmp), param_ordering, simplify );
	fi:
	retval[pos] := tmp;
	debug_print( procname, `Intersection point:`, 3, eval(tmp) );

	#
	# transform back to global frame if necessary
	#
	if coordframe=GLOBAL then
		debug_print( procname, `Transforming to GLOBAL frame...`, 1 );
		retval[pos] := transform_pos( retval[pos], localsurf[pos],
									  localsurf[dir], GLOBAL );
		debug_print( procname, `Expanding in series...`, 1 );
		retval[pos] := expansion( retval[pos], op(expansion_list) );
		retval[pos] := Collect( eval(retval[pos]), param_ordering, simplify );
	fi:

	#
	# path length from input beam anchor point to intersection point
	#
	debug_print( procname, `Calculating path length...`, 1 );
	bmag := mag( localbeam[dir] );
	bmag := expansion( bmag, op(expansion_list) );
	bmag := simpside( bmag, param_ordering, sublist, siderels, siderels_list );
	debug_print( procname, `beam dir magnitude:`, 3, bmag );
	retval[path] := collect( expansion( tval*bmag, op(expansion_list) ),
							 param_ordering,
							 simplify );
	debug_print( procname, `path length:`, 3, retval[path] );

	gc();
	RETURN( eval(retval) );
end:



`beam_intersect/solve/linear` := proc( eq, t, expansion_list, param_ordering,
									   nodenom_list )
	debug_print( procname, `Solving for t...`, 1 );
	solve( eq, t );
	debug_print( procname, `Expanding and checking for small divisors...`, 1 );
	`beam_intersect/solve/expansion`( ", expansion_list, param_ordering,
									  nodenom_list );
end:



`beam_intersect/solve/quadratic` := proc( eq, t, expansion_list, param_ordering,
										  nodenom_list )
	local tsols, tsol;

	debug_print( procname, `Solving for t...`, 1 );
	tsols := [solve( eq, t )];

	if nops(tsols)=1 then
		tsol := tsols[1];
	elif nops(tsols)=2 then
		if is_real_parts(tsols[1]) then
			tsol := tsols[1];
		elif is_real_parts(tsols[2]) then
			tsol := tsols[2];
		else
			debug_print(procname,`FATAL: both solutions are imaginary!`,0,eq);
			ERROR(`: quitting now`);
		fi:
	else
		debug_print(procname,`FATAL: not a quadratic equation in `.t.`!`,0,eq);
		ERROR(`: quitting now`);
	fi;

	debug_print( procname, `Expanding and checking for small divisors...`, 1 );
	`beam_intersect/solve/expansion`( tsol, expansion_list, param_ordering,
									  nodenom_list );
end:



`beam_intersect/solve/spheroid` := proc( eq, t, expansion_list, param_ordering,
										 nodenom_list )
	ERROR(`: not implemented yet`);
	debug_print( procname, `Solving for t...`, 1 );
	debug_print( procname, `Expanding and checking for small divisors...`, 1 );
end:



`beam_intersect/solve/nonlinear` := proc( eq, t, cc, cc_order,
										  expansion_list, param_ordering,
										  nodenom_list )
	local tparab, tsol, q, k, A, Asol, j;

	# start with solution for nearest paraboloid
	debug_print( procname, `Solving paraboloidal approximation for t...`, 1 );
	tparab := `beam_intersect/solve/quadratic`( subs(cc=0,eq), t, expansion_list,
												param_ordering, nodenom_list );

	# add power series in conic constant to paraboloid solution
	debug_print( procname, `Expanding nonlinear model...`, 1 );
	q := subs( t=tparab+sum(A[k]*cc^k,k=1..cc_order), eq );
	q := expansion( q, cc, cc_order );

	# solve for the coefficients
	debug_print( procname, `Solving for the coefficients...`, 1 );
	Asol[0] := A[0]=0;
	for k from 1 to cc_order do
		debug_print( procname, `...k=`.k.`...`, 1 );
		debug_print( procname, `......collecting & isolating...`, 2 );
		collect( coeff(q,cc,k), A[k] );
		isolate( ", A[k] );
		debug_print( procname, `......substituting & expanding...`, 2 );
		subs( seq(Asol[j],j=0..k-1), " );
		expansion( ", op(expansion_list) );
		Asol[k] := collect( ", param_ordering, factor );
		debug_print( procname, `...A[`.k.`]:`, 3, Asol[k] );
		gc();
	od;

	# assemble the solution
	debug_print( procname, `Assembling nonlinear solution...`, 1 );
	tsol := tparab;
	for k from 1 to cc_order do
		tsol := tsol + rhs(Asol[k])*cc^k;
	od;

	# collect and test for small divisors
	debug_print( procname, `Collecting and checking for small divisors...`, 1 );
	`beam_intersect/solve/expansion`( tsol, expansion_list, param_ordering,
									  nodenom_list );
end:



`beam_intersect/solve/expansion` := proc( tsol, expansion_list, param_ordering,
										  nodenom_list )
	local t;
	t := tsol;
	if expansion_list[2] > 0 then
		debug_print( procname, `Expanding in series...`, 1 );
		t := expansion( t, op(expansion_list) );
		t := collect( t, param_ordering, factor );
		if small_divisors( t, nodenom_list ) then
			debug_print(procname,`FATAL: solution has small divisors`,0,tsol);
			ERROR(`: quitting now`);
		fi;
	fi;
	t;
end:


fi;