(* :Title: 3D Vector Fields *)

(* :Context: Graphics`PlotField3D` *)

(* :Author: John M. Novak *)

(* :Summary:
This package does plots of vector fields in three dimensions.
PlotVectorField3D allows one to specify the functions describing the
three components of the field.  PlotGradientField3D plots the 
gradient vector field associated with a scalar function. 
ListPlotVectorField3D plots a three-dimensional array of vectors.
*)

(* :Mathematica Version: 2.0 *)
(* :Copyright: Copyright 1990-1999, Wolfram Research, Inc.*)

(* :Package Version: 1.0 *)

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

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

(* :Limitations: *)

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

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 PlotField.m and PlotField3D.m functions
that scales the vectors so that the longest vector displayed is of the
length specified by this option.  The default is Automatic;
at this setting, those functions that use a coordinate grid
(PlotVectorField, etc.) have the vectors scaled so that each fits within
the grid without overlapping the base of another vector. If set to None,
no rescaling will occur. In ListPlotVectorField, Automatic has the same
effect as None. This scaling is applied after ScaleFunction and
MaxArrowLength.";

ScaleFunction::usage=
"ScaleFunction is an option for the PlotField.m and PlotField3D.m functions
that 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, and not for linear scaling.  This is applied before MaxArrowLength.";

MaxArrowLength::usage=
"MaxArrowLength is an option for the PlotField.m and PlotField3D.m
functions that determines the longest 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. The initial setting is MaxArrowLength->None (no maximum).";

If[ValueQ[ColorFunction::usage] &&
  (!StringMatchQ[ColorFunction::usage, "*PlotField*"]),
  ColorFunction::usage = ColorFunction::usage <>
      " In the PlotField.m and PlotField3D.m functions, it sets the style for \
      each vector by its magnitude, with the magnitudes scaled between 0 and 1."
]

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[Graphics3D, 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, sf},
		{plotpoints, sf} = {PlotPoints, ScaleFactor}/.{opts}/.
			Options[PlotVectorField3D];
		dua = automatic[du,(u1 - u0)/(plotpoints-1)];
		dva = automatic[dv,(v1 - v0)/(plotpoints-1)];
		dwa = automatic[dw,(w1 - w0)/(plotpoints-1)];
		If[sf =!= None && !NumberQ[N[sf]],
		   sf = N[Min[dua, dva, dwa]]
        ];
		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[Graphics3D,opts],
			ScaleFactor->sf ]
	]

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: *)
