
(* :Title: 3D Vector Fields *)

(* :Context: Graphics`PlotField3D` *)

(* :Author: John M. Novak *)

(* :Summary:
	Plots vector fields in 3D
*)

(* :Mathematica Version: 2.0 *)

(* :Package Version: 1.0 *)

(* :History:
	V 1.0 April 91 by John M. Novak - based extensively on
		PlotField.m by Kevin McIsaac, Mike Chan, ECM, and John Novak
		VectorField3D.m by Wolfram Research `90 and ECM
*)

(* :Keywords:
	vector fields, gradient field, 3D graphics
*)

(* :Limitations: *)

BeginPackage["Graphics`PlotField3D`","Utilities`FilterOptions`"]

ListPlotVectorField3D::usage =
	"ListPlotVectorField3D[{{pt,vec},...},(options)] plots a
	list of vectors in three dimensions, each vector based at
	a corresponding point pt."

PlotVectorField3D::usage =
	"PlotVectorField3D[{xfunc,yfunc,zfunc},xrange,yrange,zrange]
	plots a vector field designated by the given functions,
	over the given ranges, where a range is described as
	{variable,min,max,(increment)}.  Also accepts options like
	those of ListPlotVectorField3D."

PlotGradientField3D::usage =
	"PlotGradientField3D[function,xrange,yrange,zrange,(options)]
	plots the gradient of the given scalar function, over the
	designated ranges, where a range is given as {variable,
	min,max,(increment)}."

ScaleFactor::usage=
	"ScaleFactor is an option for the PlotField3D functions that
	scales the vectors to a specified length.  Default is Automatic;
	at this setting, those functions that use a coordinate grid
	(PlotVectorField3D, etc.) have the vectors scaled to this grid.
	This scaling is applied after ScaleFunction and MaxArrowLength.";

ScaleFunction::usage=
	"ScaleFunction rescales each vector to a length determined by applying
	a pure function to the current length of that vector.  It will ignore
	vectors of 0 magnitude. Note that because this is applied before the
	ScaleFactor, this is most useful for resizing the relative lengths of the
	vectors.  This is also applied before MaxArrowLength."

MaxArrowLength::usage=
	"MaxArrowLength is an option for the PlotField3D functions
	that determines the largest vector to be drawn. The 
	value is compared to the magnitudes of all 
	the vectors and causes all longer vectors to not be
	drawn. This is applied after the ScaleFunction but before the 
	ScaleFactor. Default is None (no maximum.)";

ColorFunction::usage=
	"ColorFunction is an option for the PlotField3D functions that
	determines the color and style used to display the vectors. It
	is a pure function that accepts a value of 0 to 1; 0 corresponds
	to the shortest vector, 1 the longest.";

VectorHeads::usage =
	"VectorHeads is an option for the PlotField3D functions that
	determines whether the vectors will be displayed with heads.
	Default is VectorHeads->False."

Begin["`Private`"]

cross3[{a1_, a2_, a3_}, {b1_, b2_, b3_}] := 
	{-(a3 b2) + a2 b3, a3 b1 - a1 b3,  -(a2 b1) + a1 b2}

mag[a_] := Sqrt[Apply[Plus, a^2]]

automatic[x_, value_] :=
	If[x === Automatic, value, x]

vector3D[point:{x_, y_, z_}, grad:{dx_, dy_, dz_},False] :=
	Line[{point, point + grad}]

vector3D[point:{x_,y_,z_}, grad:{dx_,dy_,dz_},True] :=
	Point[{x,y,z}]/;grad == {0,0,0}

vector3D[point:{x_, y_, z_}, grad:{dx_, dy_, dz_},True] :=
	Module[{endpoint, perp, perpm, offsetPoint,
		   arrowA, arrowB, arrowC, arrowD},
	  endpoint = point + grad;

	  perp = cross3[grad, {0,0,1}];
	  perpm = mag[perp];
	  If[perpm == 0,
		perp = cross3[grad, {0,1,0}];
		perpm = mag[perp]
	  ];
	  perp = perp mag[grad]/(7 perpm);
	  
	  offsetPoint = point + 4/5 grad;
	  arrowA = offsetPoint + perp;
	  
	  perp = cross3[grad, perp];
	  perp = perp mag[grad]/(7 mag[perp]);
	  arrowB = offsetPoint + perp;
	  
	  perp = cross3[grad, perp];
	  perp = perp mag[grad]/(7 mag[perp]);
	  arrowC = offsetPoint + perp;
	  
	  perp = cross3[grad, perp];
	  perp = perp mag[grad]/(7 mag[perp]);
	  arrowD = offsetPoint + perp;
	  
	  {Line[{point, endpoint}], 			(* 3D arrow shaft *)
	   Line[{arrowA, endpoint, arrowC}], 		(* point of arrow *)
	   Line[{arrowB, endpoint, arrowD}], 		(* point of arrow *)
	   Line[{arrowA, arrowB, arrowC, arrowD, arrowA}] (* base of point *)
	  }
	]

Options[ListPlotVectorField3D] = 
	{ScaleFactor->Automatic, 
	 ScaleFunction->None,
	 MaxArrowLength->None,
	 ColorFunction->None,
	 VectorHeads->False};


ListPlotVectorField3D[vects:{{_?VectorQ,_?VectorQ}..},opts___] :=
	Module[{maxsize,scale,scalefunct,colorfunct,heads,points,
			vectors,mags,colors,scaledmag,allvecs,vecs=N[vects]},
		{maxsize,scale,scalefunct,colorfunct,heads} =
			{MaxArrowLength,ScaleFactor,ScaleFunction,
			ColorFunction,VectorHeads}/.{opts}/.
			Options[ListPlotVectorField3D];
		
		(* option checking *)
		If[Not[NumberQ[maxsize]] && maxsize != None,
			maxsize = None,
			maxsize = N[maxsize]];
		If[Not[NumberQ[scale]] && scale != Automatic,
			scale = Automatic,
			scale = N[scale]];
		heads = TrueQ[heads];
		
		vecs = Cases[vecs,{_,_?(VectorQ[#,NumberQ]&)}];
		{points, vectors} = Transpose[vecs];
		mags = Map[mag,vectors];
		If[colorfunct == None, colorfunct = {}&];
		If[Max[mags - Min[mags]] == 0,
			colors = Map[colorfunct,Table[0,{Length[mags]}]],
			colors = Map[colorfunct,
				(mags - Min[mags])/Max[mags - Min[mags]]]
		];

		If[scalefunct =!= None,
		 	scaledmag = (If[# == 0, 0, scalefunct[#]]&) /@ mags;
		 	vectors = MapThread[If[#2 == 0, {0,0,0}, #1 #2/#3]&,
				{vectors,scaledmag,mags}];
		 	mags = scaledmag
		   ];

		allvecs = Transpose[{colors,points,vectors,mags}];  

		If[maxsize =!= None,
		 	allvecs = Select[allvecs, (#[[4]]<=maxsize)&]
		   ];
		
		If[Max[mags] != 0,
			scale = automatic[scale,Max[mags]]/Max[mags];
 			allvecs = Map[{#[[1]],#[[2]],scale #[[3]]}&,
 				allvecs]
 		];

		(* alternate method of vector generation requires pr.
		pr = PlotRange[ Graphics3D[
				Flatten[Apply[Line[{#2,#2+#3}]&,allvecs,{1}]]]];
		*)
		
		Show[Graphics3D[
		 		Flatten[Apply[{#1,vector3D[#2,#3,heads]}&,
		 			allvecs,{1}]],
		 		FilterOptions[Graphics, opts]]]
	]/; Last[Dimensions[vects]] === 3

Options[PlotVectorField3D] =
	Join[Options[ListPlotVectorField3D],{PlotPoints->7}]

SetAttributes[PlotVectorField3D, HoldFirst]

PlotVectorField3D[f_, {u_, u0_, u1_, du_:Automatic},
			 {v_, v0_, v1_, dv_:Automatic},
			 {w_,w0_,w1_,dw_:Automatic},opts___] :=
	Module[{plotpoints,dua,dva,dwa,vecs},
		{plotpoints} = {PlotPoints}/.{opts}/.
			Options[PlotVectorField3D];
		dua = automatic[du,(u1 - u0)/(plotpoints-1)];
		dva = automatic[dv,(v1 - v0)/(plotpoints-1)];
		dwa = automatic[dw,(w1 - w0)/(plotpoints-1)];
		vecs = Flatten[Table[{N[{u,v,w}],N[f]},
			Evaluate[{u,u0,u1,dua}],Evaluate[{v,v0,v1,dva}],
			Evaluate[{w,w0,w1,dwa}]],2];
		ListPlotVectorField3D[vecs,
			FilterOptions[ListPlotVectorField3D,opts],
			FilterOptions[Graphics,opts],
			ScaleFactor->N[Min[dua,dva,dwa]]]
	]

PlotGradientField3D[function_, 
		{u_, u0__}, 
		{v_, v0__},
		{w_, w0__},
		options___] :=
	PlotVectorField3D[Evaluate[{D[function, u],
					D[function, v],D[function,w]}],
			{u, u0},
			{v, v0},
			{w, w0},
			options]

End[]

EndPackage[]

(* :Tests: *)

(* :Examples: *)
