#=============================================================================
#
# System-level AESOP routines.
#
#=============================================================================

#=============================================================================
#
# 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 OPSYS_READ <> 1 then
OPSYS_READ := 1:
read`c:/maple/startup.p`:
read`d:/optics/AESOP/objects.p`:

#==============================================================================
#
# raytrace( beam0, surfaces, small_params1, _order1,
#           small_params2, _order2, param_ordering, nodenom_list,
#           sublist, siderels, siderels_list )
#
# Trace an input ray beam0 through the optical system whose elements
# are contained in the table surfaces.  For example,
#    surfaces[1] := primary_mirror;
#    surfaces[2] := secondary_mirror;
#    surfaces[3] := detector;
# The last element of surfaces must be the detector end plane.
#
# Returns a beam with the intersection point on the detector in [pos],
# the beam direction (as reflected from the second to last surface) in
# [dir], and the total optical path length in [path].
#
#==============================================================================
raytrace := proc( beam0, surfaces, expansion_list, param_ordering,
				  nodenom_list, sublist, siderels, siderels_list )

	local i, beams, outbeam, N:

	#
	# Loop over all surfaces but the last.
	#
	beams[0] := beam0:
	N := nops(op(op(surfaces))) - 1;    #last surface is detector
	for i from 1 to N do

		print();
		debug_print( procname,
					 `Calculating reflection `.(convert(i,'string')).`...`, 1 );

		beams[i] := beamtrace( beams[i-1], surfaces[i],
							   expansion_list, param_ordering,
							   nodenom_list, sublist, siderels, siderels_list );
		debug_print( procname, `Simplifying beam...`, 1 );
		beams[i] := simpside( beams[i], param_ordering,
							  sublist, siderels, siderels_list );

		debug_print( procname,
					 `Beam after reflection `.(convert(i,'string')).`:`, 2,
					 beams[i] );

	od:

	#
	# Find the detector intersection point.
	#
	print();
	debug_print( procname, `Calculating detector intersection point...`, 1 );
	outbeam := beam_intersect( beams[N], surfaces[N+1], GLOBAL,
							   expansion_list, param_ordering, nodenom_list,
							   sublist, siderels, siderels_list );

	debug_print( procname, `Simplifying...`, 1 );
	outbeam['type'] := 'beam':
	outbeam[pos]    := simpside( outbeam[pos], param_ordering, sublist,
								 siderels, siderels_list );

	#
	# Copy output beam direction from second-to-last surface.
	#
	outbeam[dir] := copy( beams[N][dir] );

	#
	# Total path length.
	#
	debug_print( procname, `Calculating and simplifying total path length...`, 1 );
	outbeam[path] := simpside( outbeam[path]+beams[N][path], param_ordering,
							   sublist, siderels, siderels_list );

	debug_print( procname, `Detector intersection point:`, 3, outbeam[pos]  );
	debug_print( procname, `Output beam direction:`,       3, outbeam[dir]  );
	debug_print( procname, `Total optical path length:`,   3, outbeam[path] );

	RETURN( op(outbeam) );

end:



#=============================================================================
#
# beamtrace( inputbeam, opelmnt, simplepath, params,
#            order, param_ordering, nodenom_list )
#
# Pass an input beam through an optical element.
#
# Arguments:
#   inputbeam     input table of data type 'beam'
#   opelmnt       optical element (mirror, flat, etc.)
#   simplepath    if on, simplify the answer for the optical path length
#   params        small parameter list for series expansions (a Maple list)
#   order         expansion order (if order < 0, no expansions performed)
#   param_ordering       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.
#
#=============================================================================
beamtrace := proc( inputbeam::beam,
				   opelmnt::{mirror,beamcompressor,lens,beamsplitter,HOE},
				   expansion_list,
				   param_ordering::{name,list}, nodenom_list::{name,list},
				   sublist::{name,list,set}, siderels::set,
				   siderels_list::{name,list} )

	local newbeam:

	newbeam := copy( inputbeam );

	if type(opelmnt,{'mirror','beamsplitter'}) then
		newbeam := reflect( newbeam, opelmnt, expansion_list,
							param_ordering, nodenom_list,
							sublist, siderels, siderels_list );
	elif type(opelmnt,'HOE') then
		newbeam := traverse_pHOE( newbeam, opelmnt, expansion_list,
								  param_ordering, nodenom_list );
	elif type(opelmnt,'lens') then
		ERROR(`Lenses currently not implemented in AESOP`);
	else
		ERROR(`Unknown optical element type.`);
	fi:

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



#=============================================================================
#
# traverse_pHOE
#
# Diffract an input beam off a pHOE
#
#=============================================================================
traverse_pHOE := proc( inputbeam, surface, n, expansion_list,
					   param_ordering, nodenom_list )

	local   newbeam, diffvec, normvec, tmp, localsurf, ipt,
			U, sq, t1, t2, t3, c1, c2, Dc, Dc0, c10, c20, Q:

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

	#transform to local coordinate frame
	debug_print( procname, `Transforming input beam to LOCAL frame...`, 1 );
	newbeam[pos] := transform_pos( newbeam[pos], localsurf[pos],
								   localsurf[dir], LOCAL );
	newbeam[dir] := transform_dir( newbeam[dir], localsurf[dir], LOCAL );

	#find the intersection point (local coord's are returned)
	debug_print( procname, `Finding intersection point...`, 1 );
	ipt := beam_intersect( newbeam, localsurf, n, LOCAL, expansion_list,
						   param_ordering, nodenom_list );
	newbeam[pos] := copy(ipt[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, newbeam[pos] );
	normvec := evalm( normvec/mag(normvec) );

	#calculate the diffracted beam direction
	debug_print( procname, `Calculating & simplifying diffracted beam...`, 1 );
	t1 := evalm( newbeam[pos] - localsurf[cp1] );
	t1 := evalm( t1/mag(t1) );
	t1 := expansion( t1, params1, order1, params2, order2, simplepath );
	t2 := evalm( newbeam[pos] - localsurf[cp2] );
	t2 := evalm( t2/mag(t2) );
	t2 := expansion( t2, op(expansion_list) );
	t3 := expansion( evalm(newbeam[dir]/mag(newbeam[dir])),
					  op(expansion_list) );
	#*add* t2 to t1 because construction point cp2 is virtual
	U  := evalm( cross(normvec,t3) - cross(normvec,t1+t2) );
	U  := expansion( U, params1, op(expansion_list) );
	sq := sqrt( 1 - dot(U,U) );
	sq := expansion( sq, op(expansion_list) );
	diffvec := evalm( sq*normvec - cross(normvec,U) );
	diffvec := expansion( diffvec, op(expansion_list) );

	#optical path correction
	debug_print( procname, `Calculating OPL correction...`, 1 );
	c1  := mag( evalm( newbeam[pos] - localsurf[cp1] ) );
	c2  := mag( evalm( newbeam[pos] - localsurf[cp2] ) );
	Dc  := c1 + c2:
	Dc  := expansion( Dc, op(expansion_list) );
	c10 := mag( evalm( point(0,0,0) - localsurf[cp1] ) );
	c20 := mag( evalm( point(0,0,0) - localsurf[cp2] ) );
	Dc0 := c10 + c20:
	Q   := Dc - Dc0:
	Q   := expansion( Q, op(expansion_list) );

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

	newbeam[dir] := expansion( newbeam[dir], op(expansion_list) );
	newbeam[pos] := expansion( newbeam[pos], op(expansion_list) );

	#calculate the cumulative optical path
	newbeam[path] := inputbeam[path] + ipt[path] - Q:
	newbeam[path] := expansion( newbeam[path], op(expansion_list) );

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



#=============================================================================
#
# rotate_object( element, pivot, diraxis, theta )
#
# Rotate an optical element, beam, or other object.  The axis of rotation
# is the vector 'diraxis', anchored at pivot point 'pivot'.  Rotate by
# an angle 'theta'.  Direction of positive rotation is cw when viewing
# along 'diraxis'.
#
#=============================================================================
rotate_object := proc( element::{mirror,beamcompressor,lens,beamsplitter,beam,HOE},
					   pivot::point, diraxis::vector, theta::algebraic,
					   simplflag::integer )

	local new_element, tmp:

	new_element := copy( element );
	tmp := RotRay( element[pos], element[dir], pivot, diraxis, theta, simplflag );
	new_element[pos] := copy( tmp[pos] );
	new_element[dir] := copy( tmp[dir] );

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



#=============================================================================
#
# translate_object( element, dir, delta )
#
# Translate an optical element, beam, or other object.
#
#=============================================================================
translate_object := proc( element::{mirror,beamcompressor,lens,beamsplitter,beam,HOE},
						  v::vector, delta::algebraic )

	local new_element, unit_v:
	if not type( [args],
				 [ {'mirror','beamcompressor','lens','beamsplitter','beam','HOE'},
				   'vector','algebraic' ] ) then
		argserr( procname, [args] );
		ERROR();
	fi:

	new_element := copy( element );
	unit_v := evalm( v / mag(v) );
	new_element[pos] := evalm( new_element[pos] + delta*unit_v );

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



fi;