#=============================================================================
#
# Plotting 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 PLOTTING_READ <> 1 then
PLOTTING_READ := 1:

read`d:/optics/AESOP/zseries.p`:
readlib(optimize):
readlib(`optimize/makeproc`):

plothelp := proc()
	print(`plotting: [ `
	.`coordplot2D, coordplot3D, curveplot2D, curveplot3D, `
	.`\nimplicitplot2D, implicitplot3D, `
	.`jpgplot, ODEplot, PlotZernikeWavefront, \nscalarplot, surfaceplot, `
	.`tubecurve, tubecurveplot, `
	.`]`);
end:

#plothelp();


#---------------------------------------------------------------
# jpgplot( plotcommand, filename, [width,height] )
#
# Send a plot to a .jpg file.
#
# 'plotcommand' and 'filename' are required; '[width,height]'
# is optional and defaults to [640,480].
#---------------------------------------------------------------
jpgplot := proc( plotcommand, filename )

	if nargs > 1 then
		dojpgplot(args);
	fi;

	#
	# a bug prevents us from doing this here
	#
#	if op([3,2],[plotsetup()])<>`inline` or nargs>0 and args[1]=`reset` then
#		resetplot();
#	fi;
end:

dojpgplot := proc( plotcommand, filename )
	local pargs, h, w;
	pargs := [ plotdevice=jpeg, plotoutput=filename ];
	if nargs = 3 then
		w := args[3][1];
		h := args[3][2];
		pargs := [ op(pargs),
				 plotoptions=`height=`.(convert(h,string)).
							 `,width=`.(convert(w,string)) ];
	else
		pargs := [ op(pargs), plotoptions=`height=480,width=640` ];
	fi;
	plotsetup( op(pargs) );
	eval(plotcommand);
end:

#---------------------------------------------------------------
# Reset the plotting device to the screen.  A bug in Maple
# prevents including this in jpgplot() or dojpgplot().
#---------------------------------------------------------------
resetplot := proc()
	if op([3,2],[plotsetup()])<>`inline` or nargs>0 and args[1]=`reset` then
		plotsetup( plotdevice=inline, plotoutput=terminal, plotoptions=`` ):
	fi;
end:



#=============================================================================
# ODE Plotting
#=============================================================================
with(DEtools):

#---------------------------------------------------------------
# ODE solution curves
#---------------------------------------------------------------
ODEplot := proc( ODEs::{`=`,list,set}, depvars::{list,set},
				 ivar_range::{name=range}, ICs::list )

	DEplot( ODEs, depvars, ivar_range, [eval(ICs)], axes=normal,
			linecolor=navy, thickness=1, axes=normal,
			titlefont=[TIMES,ITALIC,18], axesfont=[HELVETICA,12],
			labelfont=[TIMES,ROMAN,14], scaling=unconstrained,
			method=classical[rk4],
#			method=gear[bstoer],
			args[5..nargs] );

end:



#=============================================================================
# 2D Plotting
#=============================================================================
#---------------------------------------------------------------
# curveplot2D( curve, t=t1..t2, plot_options )
# where curve = [x(t),y(t)]
#---------------------------------------------------------------
curveplot2D := proc( curve, t_limits )

	plot( [op(curve), t_limits],
		  axes=frame, numpoints=100, args[3..nargs] );

end:



#---------------------------------------------------------------
# implicitplot2D( eqn, x=x1..x2, y=y1..y2, plot_options )
#---------------------------------------------------------------
implicitplot2D := proc( eqn, x_limits, y_limits )

	plots[implicitplot]( eqn, x_limits, y_limits,
						 axes=frame, grid=[50,50],
						 args[4..nargs] );

end:



#---------------------------------------------------------------
# scalarplot( curve, t=t1..t2, plot_options )
#---------------------------------------------------------------
scalarplot := proc( curve, t_limits )

	plot( curve, t_limits, numpoints=100,
		  args[3..nargs] );

end:



#=============================================================================
# 3D Plotting
#=============================================================================
#------------------------------------------------------------
# Plot a 3D space curve
# curveplot3D( curve, t=t1..t2, plot3d_options )
# where curve = [x(t),y(t),z(t)]
#------------------------------------------------------------
curveplot3D := proc( curve, t_limits )

	plots[spacecurve]( curve, t_limits, scaling=constrained,
			numpoints=100, projection=normal, args[3..nargs] );

end:



#---------------------------------------------------------
# tubecurve( curve, t, radius, theta )
#
# Create a "tube" surrounding a space curve, including
# torsion or "twisting" effects.
#---------------------------------------------------------
tubecurve := proc( curve, t, r, theta )
	local normvec, binorm, tube;
	normvec := UnitNormal(curve,t);
	binorm  := Binormal(curve,t);
	if has(curve,sin(t)) or has(curve,cos(t)) then
		normvec := map( simplify, evalm(normvec), trig );
		binorm  := map( simplify, evalm(binorm), trig );
	fi;
	tube := evalm( curve +
			r*( cos(theta)*normvec + sin(theta)*binorm ) );
end:



#---------------------------------------------------------
# Plot a tube curve
#
# tubecurveplot( curve, radius, t=t1..t2,
#                theta=theta1..theta2,
#                plot3d_options )
#
# Use this plotting function, instead of the Maple
# tubeplot(), to explicitely show torsion effects.
#---------------------------------------------------------
tubecurveplot := proc( curve, radius, t_limits, theta_limits )

	local tube, t, theta, r;

	t     := lhs(t_limits);
	theta := lhs(theta_limits);

	#temporarily replace numerical "radius" with name "r"
	tube  := tubecurve( curve, t, r, theta );
	tube  := map( radsimp, tube );
	tube  := map( collect, tube,
					[ r, cos(theta), sin(theta) ] );
	tube  := subs( r=radius, evalm(tube) );

	plot3d( tube, t_limits, theta_limits,
			scaling=constrained, grid=[100,10], args[5..nargs] );
end:



#------------------------------------------------------------
# Plot a 3D surface
#
# surfaceplot( [x(u,v),y(u,v),z(u,v)], u=u1..u2, v=v1..v2,
#              plot3d_options )
#------------------------------------------------------------
surfaceplot := proc( eqns::{list,set}, u_limits, v_limits )

	plot3d( eqns, u_limits, v_limits,
			grid=[20,20], args[4..nargs] );

end:



#------------------------------------------------------------
# Plot a 3D surface implicitly
#
# implicitplot3D( eqn, x=x1..x2, y=y1..y2, z=z1..z2,
#                 plot3d_options )
#------------------------------------------------------------
implicitplot3D := proc( eqn, x_limits, y_limits, z_limits )

	plots[implicitplot3d]( eqn, x_limits, y_limits, z_limits,
		grid=[20,20,20], args[5..nargs] );

end:



#------------------------------------------------------------
# Plot a 3D contour plot
#
# contourplot3D( eqn, x=x1..x2, y=y1..y2, plot3d_options )
#
# additional plot options are filled=true and contours=n
#------------------------------------------------------------
contourplot3D := proc( eqn, x_limits, y_limits )

	plots[contourplot3d]( eqn, x_limits, y_limits,
		filled=true, coloring=[yellow,red],
		contours=10, args[4..nargs] );

end:



#------------------------------------------------------------
# Plot n constant curves in a given coordinate system
#
# xyeqs = [x(u,v),y(u,v)]
#------------------------------------------------------------
coordplot2D := proc( xyeqs::list, urange::{name=range},
					 vrange::{name=range}, n::integer )

	local uplots, vplots, umin, umax, vmin, vmax, i,
		  U, V, plotargs, x, y;

	umin := op([2,1],urange);
	umax := op([2,2],urange);
	vmin := op([2,1],vrange);
	vmax := op([2,2],vrange);
	plotargs := args[5..nargs];

	for i from 1 to n do
		U := umin + (i-1)/(n-1)*(umax-umin);
		uplots[i] := plot( [op(subs(lhs(urange)=U,xyeqs)), vrange], color=blue );
	od;
	for i from 1 to n do
		V := vmin + (i-1)/(n-1)*(vmax-vmin);
		vplots[i] := plot( [op(subs(lhs(vrange)=V,xyeqs)), urange], color=red );
	od;

	plots[display]( [ op(convert(uplots,list)),
					  op(convert(vplots,list)) ],
					  scaling=constrained, labels=[`x`,`y`], plotargs );

end:



#------------------------------------------------------------
# Plot three constant surfaces in a given coordinate system
#
# xyz = [x(u,v,w),y(u,v,w),z(u,v,w)]
#------------------------------------------------------------
coordplot3D := proc( xyz::list, uvaleqn::{name=numeric},
					 vvaleqn::{name=numeric}, wvaleqn::{name=numeric},
					 urange::{name=range}, vrange::{name=range},
					 wrange::{name=range} )
	local p, q, r, g, X, Y, Z, u, v, w, uval, vval, wval;
	uval := rhs(uvaleqn);
	vval := rhs(vvaleqn);
	wval := rhs(wvaleqn);
	g := 30;
	p := surfaceplot( subs(lhs(urange)=uval,xyz), vrange, wrange,
					  grid=[g,g], shading=zhue );
	q := surfaceplot( subs(lhs(vrange)=vval,xyz), urange, wrange,
					  grid=[g,g], shading=zhue ):
	r := surfaceplot( subs(lhs(wrange)=wval,xyz), urange, vrange,
					  grid=[g,g], shading=zhue ):
	plots[display]( [p,q,r], scaling=constrained,
					labels=[`X`,`Y`,`Z`], args[8..nargs] );
end:



#------------------------------------------------------------
# Plot a wavefront with specified Zernike terms subtracted
#
# PlotZernikeWavefront( opd, Zcoeffs, sublist, units, pert,
#                       f, R, plot3d_options )
#
#------------------------------------------------------------
PlotZernikeWavefront := proc( opd::algebraic, epsilon::{name,list}, c::table,
							  sublist::list, units::string, pert::{float,list},
							  focal_length::float, beam_radius::float )

	local G, F, i, p, n, m, k, scale, plotargs, tlist:

	#
	# Build up the terms to subtract from the wavefront.
	#
	p := 0:
	for i from 1 to nops(sublist) do
		n := sublist[i][1]:
		m := sublist[i][2]:
		p := p + ZernikeTerm( n, m, c );
	od:

	#
	# Subtract the desired terms and simplify.
	#
	if type(epsilon,list) then
		tlist := [op(epsilon),rho];
	else
		tlist := [epsilon,rho];
	fi;
	G := collect( opd-p, tlist, factor ):
	G := combine( G, trig ):
	G := collect( G, tlist, factor ):
	debug_print( procname, `wavefront:`, 1, G ):

	#
	# Make an optimized procedure out of what's left...
	#
	if units=`cm` then
		scale := 1:
	elif units=`mm` then
		scale := 10;
	elif units=`microns` or units=`micron` or units=`um` then
		scale := 10000:
	elif units=`nm` then
		scale := 10000000:
	elif units=`pm` then
		scale := 10000000000:
	elif units=`fm` then
		scale := 10000000000000:
	else
		ERROR(` Unknown units!`):
	fi:
	if type(epsilon,list) then
		tlist := [rho,phi,op(epsilon),f,R];
	else
		tlist := [rho,phi,epsilon,f,R];
	fi;
	F := `optimize/makeproc`( [optimize(-scale*G)], parameters=tlist ):

	#
	# ...and create a color plot from it.
	#
	if type(epsilon,list) then
		tlist := [rho,phi,op(pert),focal_length,beam_radius];
	else
		tlist := [rho,phi,pert,focal_length,beam_radius];
	fi;
	plotargs := [[rho*cos(phi), rho*sin(phi), F(op(tlist))],
				 phi=0..2*Pi, rho=0..beam_radius,
				 grid=[35,20], orientation=[110,65],
				 light=[60,90,2,2,2], ambientlight=[0.4,0.4,0.4],
				 labels=[`radius (cm)`,``,`OPD `],
				 title=`Wavefront Error (`.units.`)`];
	if nargs > 8 then
		for i from 9 to nargs do
			plotargs := [op(plotargs),args[i]];
		od;
	fi;
	plot3d( op(plotargs) );

end:




#---------------------------------------------------------------------------
# From Robert Israel
#---------------------------------------------------------------------------
evl := t -> convert(evalm(t), list):

myarrow := proc(base::{list,vector},vect::{list,vector})

## Plots an arrow as a line with a pyramidal head.
## usage: myarrow(base, vect, options)
## "base" is point for base of arrow, "vect" is tip - base
## These can be lists or vectors.
## "options" are the usual plot3d options.
## Colour is red by default.
## More options could be added to control the proportions of the arrow.

local b1,v1,stem,r,hbase,hpoint,v2,v3,head,opts;
    b1 := evalf(evl(base));
    v1 := evalf(evl(vect));
    if nops(b1) <> 3 or nops(v1) <> 3 then
        ERROR(`points need three dimensions`)
    fi;
    if not type(b1,[realcons,realcons,realcons]) and
        type(v1,[realcons,realcons,realcons]) then
        ERROR(`points must evaluate to real constants`)
    fi;
    hbase := evl(b1+.8*v1);
    hpoint := evl(b1+v1);
    opts := args[3 .. nargs];
    if indets({opts}) intersect {color,colour} = {} then
        opts := opts,colour = red
    fi;
    stem := plots[spacecurve]([b1,hbase],thickness = 3,opts);
    if v1 = [0,0,0] then RETURN(stem) fi;
    r := sqrt(linalg[norm](v1,2));
    v2 := evl(linalg[crossprod]([1,0,0],v1));
    if linalg[dotprod](v2,v2) < .01*linalg[dotprod](v1,v1) then
        v2 := evl(linalg[crossprod]([0,1,0],v1))
    fi;
    v2 := evl(.05*r*v2/linalg[norm](v2,2));
    v3 := linalg[crossprod](v1,v2);
    v3 := evl(.05*r*v3/linalg[norm](v3,2));
    head := plots[polygonplot3d]([hpoint,evl(hbase+v2),
        evl(hbase+v3), evl(hbase-v2),
        evl(hbase-v3), evl(hbase+v2)]
        ,style = patch,opts);
    plots[display]({stem,head})
end:

surfarrows := proc(fld,surf,r1,r2)
## Plots a 3d surface (parametric or cartesian) with arrows from
## a vector field at points on the surface
## Usage:
##   surfarrows(F(s,t), [x(s,t),y(s,t),z(s,t)], s=a1 .. b1, t=a2 .. b2, options)
##     for vectors F(s,t) at points [x(s,t), y(s,t), z(s,t)]
##     where F(s,t) is a list expression
##   surfarrows(F(x,y), z(x,y), x=a1 .. b1, y=a2 .. b2, options)
##     for vectors F(x,y) at points [x, y, z(x,y)]
##     where F(x,y) is a list expression
##   surfarrows(F, [x(s,t),y(s,t),z(s,t)], s=a1 .. b1, t=a2 .. b2, options)
##     for vectors F(x(s,t),y(s,t),z(s,t)) at [x(s,t),y(s,t),z(s,t)]
##     where F is a list-valued function of three variables
##   surfarrows(F, z(s,t), x=a1 .. b1, y=a2 .. b2, options)
##     for vectors F(x,y,z(s,t)] at [x,y,z(x,y)]
##     where F is a list-valued function of three variables
##   "options" are the usual 3d plot options, plus
##         arrowgrid= [n1, n2]
##     for arrows in an n1 by n2 grid, equally spaced in terms of the parameters
##     of the surface (default is arrowgrid=[8,8])
## Caution: inputs are not checked for validity

local opts,gspec,p1,p2,s2,f2,v1,v2,a1,a2,b1,b2,n1,n2,h1,h2,i,j;
    opts := args[5 .. nargs];
    v1 := op(1,r1);
    v2 := op(1,r2);
    a1 := op(1,op(2,r1));
    b1 := op(2,op(2,r1));
    a2 := op(1,op(2,r2));
    b2 := op(2,op(2,r2));
    gspec := indets({opts},identical('arrowgrid') = list);
    if gspec = {} then gspec := arrowgrid = [8,8]
    else opts := op({opts} minus gspec); gspec := gspec[1]
    fi;
    p1 := plot3d(surf,r1,r2,opts);
    opts :=
        op({opts} minus indets({opts},identical('grid') = list))
        ;
    n1 := op(2,gspec)[1]-1;
    n2 := op(2,gspec)[2]-1;
    h1 := (b1-a1)/n1;
    h2 := (b2-a2)/n2;
    if type(surf,list) then s2 := surf
    elif type(surf,vector) then s2 := evl(surf)
    else s2 := [v1,v2,surf]
    fi;
    if type(fld,procedure) then f2:= fld(op(s2))
    else f2:= evl(fld)
    fi;
    p2 := seq(seq(
     myarrow(op(subs(v1 = a1+h1*i,v2 = a2+j*h2,[s2,f2])),opts),
     i = 0 .. n1),j = 0 .. n2);
    plots[display]({p2,p1},opts)
end:

## Examples:
#
#> surfarrows((x,y,z) -> [-y, x, 0], x^2 + y^2, x=-1 .. 1, y=-1 .. 1);
#
#> torus:= [(2+cos(u))*cos(v), (2+cos(u))*sin(v), sin(u)]:
#> with(linalg):
#> unit:= V -> evl(V/norm(V,2)):
#> tnormal:= unit(crossprod(diff(torus,v),diff(torus,u))):
#> surfarrows( tnormal/2, torus, u=0 .. 2*Pi, v=0 .. 2*Pi, scaling=constrained);

#---------------------------------------------------------------------------
#---------------------------------------------------------------------------


fi: