#=============================================================================
#
# Definitions of AESOP objects and their corresponding type test functions.
#
#=============================================================================

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


#=============================================================================
#
# point( vx, vy, vz )  or  point( vector )
#
#=============================================================================

# The Maple object 'point' is singularly useless yet protected.  Hence
# we must first undo the protection.
unprotect(point);

# Now we can define our own 'point' type:
point := proc( w )
	if nargs=3 then
		array([args]);
	elif type(w,list) and size(w)=3 then
		array(w);
	elif type(w,'array'(1)) and size(w)=3 then
		array([w[1],w[2],w[3]]);
	else
		ERROR(`An AESOP 'point' must be a 3-element list or array`);
	fi;
end:



#=============================================================================
# type test for a 'point'
#
# NOTE: This will override Maple's 'point' data type.  In Maple Vr4, the routine
# limit() checks its second argument for a Maple 'point' data type, which
# after this routine is read will in general cause an error.  The lazy Maple
# programmer who wrote limit() tested for 'point' as a shortcut to test for
# an expression of the form name=expression.  This is a Maple bug.  The
# most serious effect is that when trying to do definite integration with
# Maple's int() the limit() bug will be triggered in the routine `int/DefInt`
# (you have to readlib(discont) first in order to get at `int/DefInt`).
#
# The workaround is to do definite integrations by hand: do the indefinite
# integration with int(), then manually substitute the integration limits.
#=============================================================================
`type/point` := proc( p )
	if type(p,'array'(1)) and size(p)=3 then
		true
	else
		false
	fi:
end:

read`c:/maple/limit.p`:	  #This is the fixed limit().



#=============================================================================
#
# beam( Vpos, Vdir )
#
# Create an optical "beam" originating at position Vpos and traveling
# in direction Vdir.
#
# A beam is indexed as follows:
#
#   beam[pos]	anchor position of beam (type 'point')
#   beam[dir]	direction of propogation (type 'vector')
#   beam[path]   accumulated optical path
#   beam[type] = 'beam'
#
#=============================================================================
beam := proc( Vpos::point, Vdir::vector )

	local b:

	b[pos]    := copy( Vpos );
	b[dir]    := copy( Vdir );
	b[path]   := 0:
	b['type'] := 'beam':

	RETURN( op(b) );
end:



#=============================================================================
# type test for a 'beam'
#=============================================================================
`type/beam` := proc( p )
	if type(p,table) and
	( {indices(p)} = {['pos'],['dir'],['path'],['type']} ) and
	(p['type'] = 'beam') then
		true
	else
		false
	fi:
end:



#=============================================================================
#
# optical_flat( Vpos, Vdir, coords )
#
# Create an optical flat (mirror) at position Vpos
# and pointing towards direction Vdir.
#
# A flat is indexed like this:
#
#   flat[eqn]	local coordinates equation for a plane (i.e., z=0)
#   flat[pos]	position of center of flat
#   flat[dir]	look direction vector
#   flat[coord]  labels used for the x, y, and z coordinates
#   flat[type] = 'mirror'
#
#=============================================================================
optical_flat := proc( Vpos::point, Vdir::vector, coords::list )

	local flat:

	flat[eqn]    := plane( [0,0,1], coords );
	flat[dir]    := copy( Vdir );
	flat[pos]    := copy( Vpos );
	flat[coord]  := point( coords );
	flat['type'] := 'mirror':
	flat[flen]   := infinity;
	flat[cc]	 := 0;

	RETURN( op(flat) );
end:



#=============================================================================
#
# spheroid( f, Vpos, Vdir, coords )
#
# Create a spheroid of focal length f at position Vpos
# and pointing towards direction Vdir.
#
#=============================================================================
spheroid := proc( f::algebraic, Vpos::point, Vdir::vector, coords::list )

	local p:

	p[eqn]    := quadric( f, 1, coords );
	p[pos]    := copy( Vpos );
	p[dir]    := copy( Vdir );
	p[coord]  := point( coords );
	p['type'] := 'mirror':
	p[flen]   := f;
	p[cc]     := 1;

	RETURN( op(p) );
end:



#=============================================================================
#
# paraboloid( f, Vpos, Vdir )
#
# Create a paraboloid of focal length f at position Vpos and pointing
# towards direction Vdir.
#
#=============================================================================
paraboloid := proc( f::algebraic, Vpos::point, Vdir::vector, coords::list )

	local p:

	p[eqn]    := quadric( f, 0, coords );
	p[pos]    := copy( Vpos );
	p[dir]    := copy( Vdir );
	p[coord]  := point( coords );
	p['type'] := 'mirror':
	p[flen]   := f;
	p[cc]     := 0;

	RETURN( op(p) );
end:



#=============================================================================
#
# conicoid( f, conconst, Vpos, Vdir )
#
# Create a general conicoid of focal length f at position Vpos and pointing
# towards direction Vdir.  The conic constant is conconst.
#
#=============================================================================
conicoid := proc( f::algebraic, conconst::algebraic,
				  Vpos::point, Vdir::vector, coords::list )

	local p:

	p[eqn]    := rootfunc( quadric( f, conconst, coords ), collect, conconst );
	p[pos]    := copy( Vpos );
	p[dir]    := copy( Vdir );
	p[coord]  := point( coords );
	p['type'] := 'mirror':
	p[flen]   := f;
	p[cc]     := conconst;

	RETURN( op(p) );
end:



#=============================================================================
#
# asphere( f, conconst, A, Vpos, Vdir )
#
#=============================================================================
asphere := proc( f::algebraic, conconst::algebraic, A::list,
				 Vpos::point, Vdir::vector, coords::list )

	local p, r, r2, r4, x, y, z;

	x := coords[1];
	y := coords[2];
	z := coords[3];

	r2        := x^2 + y^2;
	r4        := r2^2;
	r         := sqrt(r2);
	p[eqn]    := quadric( f, conconst, coords )
				 + sum( A[n]*r^n, n=1..nops(A) );
	p[pos]    := copy( Vpos );
	p[dir]    := copy( Vdir );
	p[coord]  := point( coords );
	p['type'] := 'mirror';
	p[flen]   := f;
	p[cc]     := conconst;

	RETURN( op(p) );
end:



#=============================================================================
# type test for a 'mirror'
#=============================================================================
`type/mirror` := proc( p )
	if type(p,table) and
	   ( {indices(p)} = {['pos'],['dir'],['eqn'],['coord'],['type'],['flen'],['cc']} ) and
	   (p['type'] = 'mirror') then
		true
	else
		false
	fi:
end:



#=============================================================================
#
# pHOE( f, Vpos, Vdir, c1, c2 )
#
# Create a simplified, focussing HOE on a paraboloidal mirror of focal length
# f at position Vpos and pointing towards direction Vdir.  It will be assumed
# that c2 is virtual, so that a beam starting from c1 will diffract to a
# focus at c2 (see TM93-04).  Further simplifications/restrictions are that
# the diffraction order m=1, and the readout wavelength is equal to the
# construction wavelength.  c1 and c2 must be specified in the LOCAL frame.
#
#=============================================================================
pHOE := proc( f::algebraic, Vpos::point, Vdir::vector, c1::point, c2::point,
			  coords::list )

	local p:

	p[eqn]    := quadric( f, 0, coords );
	p[pos]    := copy( Vpos );
	p[dir]    := copy( Vdir );
	p[coord]  := point( coords );
	p['type'] := 'HOE':
	p[cp1]    := copy( c1 );
	p[cp2]    := copy( c2 );
	p[flen]   := f;
	p[cc]     := 0;

	RETURN( op(p) );
end:



#=============================================================================
# type test for a 'HOE'
#=============================================================================
`type/HOE` := proc( p )
	if type(p,table) and
	   ( {indices(p)} = {['pos'],['dir'],['eqn'],['coord'],
						 ['type'],['cp1'],['cp2']} ) and
	   (p['type'] = 'HOE') then
		true
	else
		false
	fi:
end:



#=============================================================================
#
# beam_splitter( Vpos, Vdir )
#
# Create a beam splitter at position Vpos and pointing towards direction Vdir.
#
#=============================================================================
beam_splitter := proc( Vpos::point, Vdir::vector, coords::list )

	local bs:

	bs[eqn]    := plane( [0, 0, 1], coords );
	bs[dir]    := copy( Vdir );
	bs[pos]    := copy( Vpos );
	bs[coord]  := point( coords );
	bs['type'] := 'beamsplitter':

	RETURN( op(bs) );
end:



#=============================================================================
# type test for a 'beamsplitter'
#=============================================================================
`type/beamsplitter` := proc( p )
	if type(p,table) and
	   ( {indices(p)} = {['pos'],['dir'],['eqn'],['coord'],['type']} ) and
	   (p['type'] = 'beamsplitter') then
		true
	else
		false
	fi:
end:



#=============================================================================
#
# beam_compressor( M, Vpos, Vdir )
#
# Create a beam compressor with magnification factor M which collects
# light from (i.e., it "looks" towards) direction Vdir.
#
# A beam compressor is indexed like this:
#
#   bc[pos]	position of center of beam compressor
#   bc[dir]	look direction vector
#   bc[mag]	magnification factor
#   bc[type] = 'beamcompressor'
#
#=============================================================================
beam_compressor := proc( M::algebraic, Vpos::point, Vdir::vector )

	local bc:

	bc[pos]    := copy( Vpos );
	bc[dir]    := copy( Vdir );
	bc['mag']  := copy( M );
	bc['type'] := 'beamcompressor':

	RETURN( op(bc) );
end:



#=============================================================================
# type test for a 'beamcompressor'
#=============================================================================
`type/beamcompressor` := proc( p )
	if type(p,table) and
	( {indices(p)} = {['pos'],['dir'],['mag'],['type']} ) and
	(p['type'] = 'beamcompressor') then
		true
	else
		false
	fi:
end:



#=============================================================================
# type test for a 'lens'
#=============================================================================
`type/lens` := proc( p )
	if type(p,table) and ( {indices(p)} = {['x'],['y'],['z'],['type']} ) and
	(p['type'] = 'lens') then
		true
	else
		false
	fi:
end:



fi;