home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / GRAPHICS.PAK / PLOTFIEL.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  8.0 KB  |  269 lines
  1.  
  2. (*:Version: Mathematica 2.0 *)
  3.  
  4. (*:Name: Graphics`PlotField` *)
  5.  
  6. (*:Title: Vector Field Plots of 2D Vector Functions *)
  7.  
  8. (*:Author:
  9.     Kevin McIsaac, Wolfram Research, Inc.
  10.     Updated by Mike Chan and Bruce Sawhill, Wolfram Research, Inc.,
  11.     September 1990
  12.     Modified April, 1991, by John M. Novak
  13. *)
  14.  
  15. (*:Keywords:
  16.     vector field plot, 2D Vector Functions, Polya representation
  17. *)
  18.  
  19. (*:Requirements: none. *)
  20.  
  21. (*:Warnings: none. *)
  22.  
  23. (*:Sources:
  24.  
  25. *)
  26.  
  27. (*:Summary: This package does plots of vector fields in the plane.
  28. You give the vector field as a pair of functions in the function
  29. PlotVectorField. If you give
  30. a scalar function, you can use PlotGradientField or
  31. PlotHamiltonianField to plot its gradient vector field or
  32. Hamiltonian vector field respectively.
  33. *)
  34.  
  35. BeginPackage["Graphics`PlotField`",
  36.     "Geometry`Rotations`", 
  37.     "Utilities`FilterOptions`"]
  38.  
  39. PlotVectorField::usage =
  40. "PlotVectorField[f, {x, x0, x1, (xu)}, {y, y0, y1, (yu)}, (options)]
  41. produces a vector field plot of the two-dimensional vector function f.";
  42.  
  43. PlotGradientField::usage = "PlotGradientField[f, {x, x0, x1, (xu)},
  44. {y, y0, y1, (yu)}, (options)] produces a vector field plot of the
  45. gradient vector field of the scalar function f by calculating its
  46. derivatives analytically.";
  47.  
  48. PlotHamiltonianField::usage = "PlotHamiltonianField[f, {x, x0, x1,
  49. (xu)}, {y, y0, y1, (yu)}, (options)] produces a vector field plot
  50. of the Hamiltonian vector field of the scalar-valued function f by
  51. calculating its derivatives analytically.";
  52.  
  53. PlotPolyaField::usage= "PlotPolyaField[f[x + I y], {x, x0, x1,
  54. (dx)}, {y, y0, y1, (dy)}, options] plots the function f in the
  55. complex plane using Polya representation.";
  56.  
  57. ListPlotVectorField::usage = "ListPlotVectorField[{{vec11,vec12,..},...}]
  58. accepts a rectangular array of vectors (larger than 2x2) and displays
  59. them, with each vector positioned in the same location graphically
  60. as the matrix would be (ie, vector 1,1 in the upper right corner).
  61. ListPlotVectorField[{{pt,vec},{pt,vec},...] displays a list of
  62. vectors, each based at the corresponding point."
  63.  
  64. ScaleFactor::usage=
  65. "ScaleFactor is an option for the PlotField functions that
  66. scales the vectors to a specified length.  Default is Automatic;
  67. at this setting, those functions that use a coordinate grid
  68. (PlotVectorField, etc.) have the vectors scaled to this grid.
  69. This scaling is applied after ScaleFunction and MaxArrowLength.";
  70.  
  71. ScaleFunction::usage=
  72. "ScaleFunction rescales each vector to a length determined by applying
  73. a pure function to the current length of that vector.  It will ignore
  74. vectors of 0 magnitude. Note that because this is applied before the
  75. ScaleFactor, this is most useful for resizing the relative lengths of the
  76. vectors.  This is also applied before MaxArrowLength."
  77.  
  78. MaxArrowLength::usage=
  79. "MaxArrowLength is an option for the PlotField functions
  80. that determines the largest vector to be drawn. The 
  81. value is compared to the magnitudes of all 
  82. the vectors and causes all longer vectors to not be
  83. drawn. This is applied after the ScaleFunction but before the 
  84. ScaleFactor. The initial setting is MaxArrowLength->None (no maximum).";
  85.  
  86. ColorFunction::usage=
  87. "ColorFunction is an option for the PlotField functions that
  88. determines the color and style used to display the vectors. It
  89. is a pure function that accepts a value of 0 to 1; 0 corresponds
  90. to the shortest vector, 1 the longest.";
  91.  
  92. Begin["`Private`"]
  93.  
  94. rot45 = N[RotationMatrix2D[7Pi/8]];
  95. rotm45 = N[RotationMatrix2D[-7Pi/8]];
  96. vector[{x_, y_},{dx_,dy_}, 
  97.             {{x0_,x1_},{y0_,y1_}}, ratio_] :=
  98.     Block[{xr, yr, p1, p2, end, start, 
  99.             direction},
  100.         xr = x1 - x0;
  101.         yr = y1 - y0;
  102.         start = {(x-x0)/xr,(y-y0)/yr}{1,ratio};
  103.         end = {(x+dx-x0)/xr,(y+dy-y0)/yr}{1,ratio};
  104.         direction = (end-start)/3;
  105.         p1 = end + (rot45 . direction);
  106.         p2 = end + (rotm45 . direction);
  107.         If[NumberQ[dx] && NumberQ[dy],
  108.             {Line[{Scaled[start], Scaled[end]}],
  109.              Line[{Scaled[p1], Scaled[end],
  110.                        Scaled[p2]}]},
  111.          (* else *)
  112.          {PointSize[.015],Point[{x,y}]}
  113.     ]
  114. ]
  115. automatic[x_, value_] :=
  116.     If[x === Automatic, value, x]
  117.  
  118. magnitude[v_List] := Sqrt[Plus @@ (v^2)]
  119.  
  120. ListPlotVectorField::ragged =
  121.     "ListPlotVectorField requires a rectangular array of vectors."
  122.  
  123. Options[ListPlotVectorField] = 
  124.     {ScaleFactor->Automatic, 
  125.      ScaleFunction->None,
  126.      MaxArrowLength->None,
  127.      ColorFunction->None,
  128.      AspectRatio->1};
  129.  
  130. ListPlotVectorField[ vects:{{_?VectorQ,_?VectorQ}..}, opts___] :=
  131.     Module[{maxsize,scale,scalefunct,colorfunct,points,
  132.             vectors,colors,mags,scaledmag,allvecs,
  133.             vecs = N[vects]},
  134.         {maxsize,scale,scalefunct,colorfunct} =
  135.             {MaxArrowLength,ScaleFactor,ScaleFunction,
  136.             ColorFunction}/.{opts}/.Options[ListPlotVectorField];
  137.  
  138.         If[Not[NumberQ[N[maxsize]] && maxsize =!= None],
  139.             maxsize = None,
  140.             maxsize = N[maxsize]];
  141.         If[Not[NumberQ[N[scale]]] &&
  142.                 scale =!= Automatic,
  143.             scale = Automatic,
  144.             scale = N[scale]];
  145.  
  146.         vecs = Cases[vecs,{_,_?(VectorQ[#,NumberQ]&)}];
  147.         {points, vectors} = Transpose[vecs];
  148.         mags = Map[magnitude,vectors];
  149.         If[colorfunct == None, colorfunct = {}&];
  150.         If[Max[mags - Min[mags]] == 0,
  151.             colors = Map[colorfunct,Table[0,{Length[mags]}]],
  152.             colors = Map[colorfunct,
  153.                 (mags - Min[mags])/Max[mags - Min[mags]]]
  154.         ];
  155.  
  156.         If[scalefunct =!= None,
  157.              scaledmag = (If[# == 0, 0, scalefunct[#]]&) /@ mags;
  158.              vectors = MapThread[If[#2 == 0, {0,0}, #1 #2/#3]&,
  159.                 {vectors,scaledmag,mags}];
  160.              mags = scaledmag
  161.            ];
  162.  
  163.         allvecs = Transpose[{colors,points,vectors,mags}];  
  164.  
  165.         If[maxsize =!= None,
  166.              allvecs = Select[allvecs, (#[[4]]<=maxsize)&]
  167.            ];
  168.         
  169.         scale = automatic[scale,Max[mags]]/Max[mags];
  170.  
  171.          allvecs = Map[{#[[1]],#[[2]],scale #[[3]]}&,
  172.              allvecs];
  173.          
  174.         pr = PlotRange[ Graphics[
  175.                 Flatten[Apply[Line[{#2,#2+#3}]&,allvecs,{1}]]]];
  176.         ratio = AspectRatio /.{opts}/.Options[PlotVectorField];
  177.         
  178.         Show[Graphics[
  179.                  Flatten[Apply[{#1,vector[#2,#3,pr,ratio]}&,
  180.                      allvecs,{1}]],
  181.                  PlotRange->pr,
  182.                  FilterOptions[Graphics, opts], 
  183.                 AspectRatio -> ratio]]
  184.              
  185.     ]/; Last[Dimensions[vects]] == 2
  186.  
  187.     
  188. ListPlotVectorField[ vects:{{(_?VectorQ)..}..},opts___] :=
  189.     Module[{dim = Dimensions[vects],flag = True,vecs = Reverse[vects]},
  190.         flag = And @@ Map[SameQ[Dimensions[#],
  191.             Dimensions[vecs[[1]]]]&, vecs];
  192.         If[flag,
  193.             ar = Array[{#2,#1}&,Evaluate[Take[dim,2]]];
  194.                 ar = Flatten[MapThread[List,{ar,vecs},2],1],
  195.             Message[ListPlotVectorField::ragged]
  196.         ];
  197.     ListPlotVectorField[ar,opts]/;flag
  198.     ]/; Last[Dimensions[vects]] == 2
  199.  
  200. Options[PlotVectorField] =
  201.     Join[Options[ListPlotVectorField],{PlotPoints->15}]
  202.  
  203. SetAttributes[PlotVectorField, HoldFirst]
  204.  
  205. PlotVectorField[f_, {u_, u0_, u1_, du_:Automatic},
  206.              {v_, v0_, v1_, dv_:Automatic}, opts___] :=
  207.     Module[{plotpoints,dua,dva,vecs},
  208.         {plotpoints} = {PlotPoints}/.{opts}/.Options[PlotVectorField];
  209.         dua = automatic[du,(u1 - u0)/(plotpoints-1)];
  210.         dva = automatic[dv,(v1 - v0)/(plotpoints-1)];
  211.         vecs = Flatten[Table[{N[{u,v}],N[f]},
  212.             Evaluate[{u,u0,u1,dua}],Evaluate[{v,v0,v1,dva}]],1];
  213.         ListPlotVectorField[vecs,
  214.             FilterOptions[ListPlotVectorField,opts],
  215.             FilterOptions[Graphics,opts],
  216.             ScaleFactor->N[Min[dua,dva]]]
  217.     ]
  218.  
  219. PlotGradientField[function_, 
  220.         {u_, u0__}, 
  221.         {v_, v0__},
  222.         options___] :=
  223.     PlotVectorField[Evaluate[{D[function, u], D[function, v]}],
  224.                 {u, u0},
  225.                 {v, v0},
  226.                 options]
  227.  
  228. PlotHamiltonianField[function_, 
  229.         {u_, u0__}, 
  230.         {v_, v0__},
  231.         options___] :=
  232.     PlotVectorField[Evaluate[{D[function, v], -D[function, u]}],
  233.                 {u, u0},
  234.                 {v, v0},
  235.                 options]
  236.  
  237.                 
  238. SetAttributes[PlotPolyaField, HoldFirst]
  239. PlotPolyaField[f_, x_List, y_List, opts___] :=
  240.     PlotVectorField[{Re[#], -Im[#]} & @ f, x, y, opts,
  241.       ScaleFunction->(Log[#+1]&)]
  242.                 
  243. End[]   (* Graphics`PlotField`Private` *)
  244.  
  245. EndPackage[]    (* Graphics`PlotField` *)
  246.  
  247. (*:Limitations: none known. *)
  248.  
  249.  
  250. (*:Examples:
  251.  
  252. PlotVectorField[ {Sin[x],Cos[y]},{x,0,Pi},{y,0,Pi}]
  253.  
  254. PlotVectorField[ { Sin[x y], Cos[x y] },{x,0,Pi},{y,0,Pi}]
  255.  
  256. PlotGradientField[ x^3 + y^4,{x,0,10},{y,0,10}]
  257.  
  258. PlotPolyaField[ (x+I y)^4,{x,5,10},{y,5,10}]
  259.  
  260.  
  261. *)
  262.  
  263.  
  264.  
  265.  
  266.  
  267.  
  268.  
  269.