(* :Package Version: 2.4 *)

(* :Mathematica Version: 2.2 *)

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

(* :Name: Graphics`PlotField` *)

(* :Title: Vector Field Plots of 2D Vector Functions *)

(* :Author:
    Kevin McIsaac, Wolfram Research, Inc.
    Updated by Mike Chan and Bruce Sawhill, Wolfram Research, Inc.,
    September 1990.
    Modified April 1991, by John M. Novak.
    V2.3, July/October 1992, by John M. Novak--use Arrow.m package.
    V2.4 January 1994 by John M. Novak. Various revisions and improvements.
*)

(* :Keywords:
    vector field plot, 2D Vector Functions, Polya representation
*)

(* :Requirements: None. *)

(* :Warnings: None. *)

(* :Sources: *)

(*:Summary:
This package does plots of vector fields in the plane.
PlotVectorField allows one to specify the functions describing the
two components of the field.  PlotGradientField and PlotHamiltonianField
plot the respective vector fields associated with a scalar function.
PlotPolyaField plots the field associated with a complex-valued
function.  ListPlotVectorField plots a rectangular array of vectors.
*)

BeginPackage["Graphics`PlotField`", 
    "Utilities`FilterOptions`",
    "Graphics`Common`GraphicsCommon`",
    "Graphics`Arrow`"]

PlotVectorField::usage =
"PlotVectorField[f, {x, x0, x1, (xu)}, {y, y0, y1, (yu)}, (options)]
produces a vector field plot of the two-dimensional vector function f.";

PlotGradientField::usage = "PlotGradientField[f, {x, x0, x1, (xu)},
{y, y0, y1, (yu)}, (options)] produces a vector field plot of the
gradient vector field of the scalar function f by calculating its
derivatives analytically.";

PlotHamiltonianField::usage = "PlotHamiltonianField[f, {x, x0, x1,
(xu)}, {y, y0, y1, (yu)}, (options)] produces a vector field plot
of the Hamiltonian vector field of the scalar-valued function f by
calculating its derivatives analytically.";

PlotPolyaField::usage= "PlotPolyaField[f[x + I y], {x, x0, x1,
(dx)}, {y, y0, y1, (dy)}, options] plots the function f in the
complex plane using Polya representation.";

ListPlotVectorField::usage =
"ListPlotVectorField[{{vec11,vec12,..},...}] accepts a rectangular array
of two-dimensional vectors (larger than 2x2) and displays
them, with each vector positioned in the same location graphically
as the matrix would be (ie, vector 1,1 in the upper left corner).
ListPlotVectorField[{{pt,vec},{pt,vec},...}] displays a list of
vectors, each based at the corresponding point.";

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."
]

Begin["`Private`"]

(* --- utilities --- *)
(* local numberQ - supercede by NumericQ after V3.0. *)
numberQ[n_] := NumberQ[N[n]]

(* little utility that is really a holdover from an early version of
   the package; returns the second arg if the first arg is Automatic,
   otherwise return the first arg. *)
automatic[x_, value_] :=
    If[x === Automatic, value, x]

(* utility to compute magnitude of vector *)
magnitude[v_List] := Sqrt[v . v]

(* --- plotting functions --- *)
(* Plot a list of {base, vector} pairs or a matrix of vectors (assumed
   to be placed on an integer base grid) *)
   
ListPlotVectorField::lpvf =
"ListPlotVectorField requires a rectangular array of vectors or a list
of {base, vector} pairs.";

Options[ListPlotVectorField] = 
    Join[{ScaleFactor->Automatic, 
     ScaleFunction->None,
     MaxArrowLength->None,
     ColorFunction->None,
     AspectRatio->1}, Options[Arrow], Options[Graphics]];

SetOptions[ListPlotVectorField,
            PlotRange -> All,
            AspectRatio -> Automatic,
            HeadScaling -> Automatic,
            HeadLength -> 0.02];

ListPlotVectorField[ vects:{{{_,_},{_,_}}..}, opts___?OptionQ] :=
    Module[{maxsize,scale,scalefunct,colorfunct,points,
            vectors,colors,mags,scaledmag,allvecs,
            vecs = N[vects], arropts},
      (* -- get option values -- *)
        {maxsize,scale,scalefunct,colorfunct} =
            {MaxArrowLength,ScaleFactor,ScaleFunction,
            ColorFunction}/.Flatten[{opts}]/.Options[ListPlotVectorField];
      (* select things that can only be vectors from the input *)
        vecs = Cases[vecs,
               {{_?numberQ, _?numberQ}, {_?numberQ, _?numberQ}},
               Infinity];
        {points, vectors} = Transpose[vecs];
        mags = Map[magnitude,vectors];
      (* -- determine the colors -- *)
      (* if the colorfunction is None, cause it to generate empty lists *)
        If[colorfunct == None, colorfunct = {}&];
      (* if all vectors are the same size, make list of colorfunct[0],
          else map the color function across the magnitudes *)
        If[Equal @@ mags,
            colors = Table[Evaluate[colorfunct[0]],{Length[mags]}],
            colors = Map[colorfunct,
                (mags - Min[mags])/Max[mags - Min[mags]]]
        ];
      (* -- scale vectors by scale function -- *)
        If[scalefunct =!= None,
             scaledmag = Map[If[# == 0, 0, scalefunct[#]]&, mags];
             {vectors, mags} = Transpose[MapThread[
                  If[#3 == 0 || !numberQ[#2], {{0,0}, 0}, {#1 #2/#3, #2}]&,
                  {vectors,scaledmag,mags}
              ]]
        ];
      (* regroup colors, points, and mags with the associated vectors *)
        allvecs = Transpose[{colors,points,vectors,mags}];  
      (* pull all vectors with magnitude greater than MaxArrowLength *)
        If[numberQ[maxsize],
             allvecs = Select[allvecs, (#[[4]] <= N[maxsize])&]
        ];
      (* calculate scale factor *)
        If[numberQ[scale],
            scale = scale/Max[mags],
            scale = 1
        ];
      (* filter out arrow options *)
         arropts = FilterOptions[Arrow, ##]& @@
            Flatten[{opts, Options[ListPlotVectorField]}];
      (* turn the vectors into Arrow objects *)
        allvecs = Apply[
            Flatten[{#1, Arrow[#2, #2 + scale #3, arropts]}]&,
            allvecs, {1}];
      (* -- show the vector field plot -- *)
      (* note that line thickness is forced to 0.0001 (thin lines);
         this can be overridden by use of ColorFunction option *)
        Show[Graphics[
             {Thickness[0.0001], allvecs},
             FilterOptions[Graphics, ##]& @@
                Flatten[{opts, Options[ListPlotVectorField]}] 
        ]]
    ]

(* given a matrix of vectors with no specified bases, generate base points
   on an integer grid, and pass back to ListPlotVectorField. *)    
ListPlotVectorField[ vects_List?(TensorRank[#] === 3 &),opts___] :=
    ListPlotVectorField[
        Flatten[MapIndexed[{Reverse[#2], #1}&, Reverse[vects], {2}], 1],
        opts
    ]

ListPlotVectorField[v_, ___] := Null/;(
    Message[ListPlotVectorField::lpvf]; False
    )
    
(* PlotVectorField takes a function that generates vectors *)
Options[PlotVectorField] =
    Join[Options[ListPlotVectorField], {PlotPoints -> Automatic}];

SetAttributes[PlotVectorField, HoldFirst];

(* Note: the slightly odd specification of range increment is for backwards
  compatibility. Users should preferentially use PlotPoints, and not the
  range increment arguments; however, so that it will usually still work
  the way it used to, we allow PlotPoints to take the value Automatic;
  when set, the range argument is used. (If the range argument is also
  Automatic, the increment is treated as if PlotPoints were set to 15.)
  Otherwise, the value of PlotPoints will override whatever is specified
  in the range argument. Eventually, we should phase out the range
  argument completely. --JMN Jan. 94 *)
PlotVectorField[f_, {u_, u0_?numberQ, u1_?numberQ, du_:Automatic},
             {v_, v0_?numberQ, v1_?numberQ, dv_:Automatic}, opts___?OptionQ] :=
    Module[{plotpoints, dua, dva, vecs, xpp, ypp, sf},
      (* -- grab options -- *)
        {plotpoints, sf} = {PlotPoints, ScaleFactor}/.Flatten[{opts}]/.
             Options[PlotVectorField];
        If[Head[plotpoints] === List,
            xpp = First[plotpoints];
            ypp = Last[plotpoints],
          (* else *)
            xpp = ypp = plotpoints
        ];
      (* determine interval between bases of vectors *)
        If[!IntegerQ[xpp],
            dua = automatic[du,(u1 - u0)/14],
            dua = (u1 - u0)/(xpp - 1)
        ];
        If[!IntegerQ[ypp],
            dva = automatic[dv,(v1 - v0)/14],
            dva = (v1 - v0)/(ypp - 1)
        ];
      (* set the scaling factor based on the intervals if it is not
            explicitly None or a number *)
        If[ sf =!= None && !numberQ[sf],
            sf = N[Min[dua, dva]]
        ];
      (* -- determine the vectors -- *)
        vecs = Flatten[Table[{N[{u,v}],N[f]},
            Evaluate[{u,u0,u1,dua}],Evaluate[{v,v0,v1,dva}]],1];
      (* call ListPlotVectorField *)
        ListPlotVectorField[vecs,
            FilterOptions[ListPlotVectorField,##]& @@
                Flatten[{ScaleFactor -> sf, opts, Options[PlotVectorField]}]
        ]/;MatchQ[vecs, {{_?VectorQ, _?VectorQ}..}]
    ]

(* PlotGradientField - computes the gradient of a scalar function and
   calls PlotVectorField on the result *)
Options[PlotGradientField] = Options[PlotVectorField];

PlotGradientField[function_, 
        {u_, u0__}, 
        {v_, v0__},
        options___] :=
    PlotVectorField[Evaluate[{D[function, u], D[function, v]}],
                {u, u0},
                {v, v0},
                options, Options[PlotGradientField]]

(* PlotHamiltonianField - computes the hamiltonian field from a scalar
   function and passes it on to PlotVectorField *)
Options[PlotHamiltonianField] = Options[PlotVectorField];

PlotHamiltonianField[function_, 
        {u_, u0__}, 
        {v_, v0__},
        options___] :=
    PlotVectorField[Evaluate[{D[function, v], -D[function, u]}],
                {u, u0},
                {v, v0},
                options, Options[PlotHamiltonianField]]


(* PlotPolyaField takes a complex scalar function and produces a vector
   field in the complex plane, with vector magnitudes scaled as indicated. *) 
Options[PlotPolyaField] = Options[PlotVectorField];

SetOptions[PlotPolyaField, ScaleFunction -> (Log[# + 1]&) ]

SetAttributes[PlotPolyaField, HoldFirst]

PlotPolyaField[f_, x_List, y_List, opts___] :=
    PlotVectorField[{Re[#], -Im[#]} & @ f, x, y, opts,
      Options[PlotPolyaField] ]
                
End[]   (* Graphics`PlotField`Private` *)

EndPackage[]    (* Graphics`PlotField` *)

(*:Limitations: None known. *)


(*:Examples:

PlotVectorField[ {Sin[x],Cos[y]},{x,0,Pi},{y,0,Pi}]

PlotVectorField[ { Sin[x y], Cos[x y] },{x,0,Pi},{y,0,Pi}]

PlotGradientField[ x^3 + y^4,{x,0,10},{y,0,10}]

PlotPolyaField[ (x+I y)^4,{x,5,10},{y,5,10}]


*)







