(* :Name: NumericalMath`OrderStar` *)

(*
   :Title: Order Stars For Approximants To Functions.
*)

(* :Author: Mark Sofroniou *)

(* :Summary:
This package plots the order star of an approximating function,
to an essentially analytic function.  It is common to consider
rational approximants to functions such as Pade approximants.
Various information about a numerical scheme (such as order and stability)
may be ascertained from its order star. For example, Runge-Kutta methods
may be considered as rational approximants to the exponential,
where relative and absolute stability regions are considered in
terms of the linear scalar test problem of Dahlquist.
The zeros, poles and interpolation points convey important additional
information and may also be displayed.
*)

(* :Context: NumericalMath`OrderStar` *)

(* :Package Version: 1.0 *)

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

(* :History:
 Original Version by Mark Sofroniou, January, 1993.
 Updated with suggestions from Jerry Keiper, December 1993.
 Revised for release, July 1995.
*)

(* :Keywords:
 Numerical Integration, Runge-Kutta, ODE, Order Star, Stability.
*)

(* :Source:
 Mark Sofroniou, Ph.D. Thesis (1994), Loughborough University,
 Loughborough, Leicestershire LE11 3TU, England.

 For a comprehensive study see Order Stars, A. Iserles &
 S. P. Norsett, Chapman & Hall, 1991.
*)

(* :Mathematica Version: 3.0 *)

(* :Limitations:
 The package relies on the code for ContourPlot to draw
 the order star. ContourPlot evaluates the function on a fixed
 mesh and, in some cases, may not produce a smooth contour
 unless the option PlotPoints is set to a fairly high number
 (several hundred). This results in small sub-divisions of
 the plot region and computation may be slow. Precision may
 also need to be increased in certain cases. The option 
 OrderStarSubPlots has been added to address these issues.
 The values for the poles and zeros of the function and
 approximant and the interpolation points are found using
 NSolve and may not always give full solutions (or indeed any
 at all). Therefore, values may be specified using options.
*)

BeginPackage["NumericalMath`OrderStar`","Calculus`Pade`",
	"Utilities`FilterOptions`", "Graphics`FilledPlot`"]

OrderStar::usage = "OrderStar[r, f] plots the order star
of the approximating function r to the function f.
OrderStar[r, f, var] may be used to specify the variable explicitly."

OrderStarInterpolation::usage = "OrderStarInterpolation is an
option to OrderStar specifying whether interpolation points of
an approximant to a function should be displayed.
OrderStarInterpolation may evaluate to True, False or a list
of {x,y} coordinates (useful if NSolve fails to detect solutions)."

OrderStarKind::usage = "OrderStarKind is an option to
OrderStar specifying the type of order star to be displayed.
OrderStarKind may be set to First or to Second.
Order stars of the first kind trace out the level curve Abs[R/F]==1.
Order stars of the second kind trace out the level curve Re[R-F]==0."

OrderStarLegend::usage = "OrderStarLegend is an option to OrderStar
specifying whether (or where) to display the legend of symbols
used to represent zeros, poles and interpolation points.
OrderStarLegend may evaluate to True, False or {{xmin,ymin},{xmax,ymax}}
where the x,y values are scaled coordinates ranging from 0 to 1."

OrderStarPoles::usage = "OrderStarPoles is an option to OrderStar
specifying whether poles of an approximant and a function should be
displayed. OrderStarPoles may evaluate to any pair consisting of
True, False or a list of {x,y} coordinates (useful if NSolve fails
to detect solutions)."

OrderStarZeros::usage = "OrderStarZeros is an option to OrderStar
specifying whether zeros of an approximant and a function should be
displayed. OrderStarZeros may evaluate to any pair consisting of
True, False or a list of {x,y} coordinates (useful if NSolve fails
to detect solutions)."

OrderStarSubPlots::usage = "OrderStarSubPlots is an option to OrderStar
used to resolve fine features in specified plot regions. Since ContourPlot
works on a fixed mesh some features may be inaccurately represented -
a rational approximant to the exponential near the origin is one
such example. OrderStarSubPlots must evaluate to False or a list
containing ContourGraphics objects or lists of PlotRange and
PlotPoints option rules."

OrderStarSymbolSize::usage = "OrderStarSymbolSize is an option
of OrderStar specifying the size of the symbols used to represent
poles, zeros and interpolation points."

OrderStarSymbolThickness::usage = "OrderStarSymbolThickness is
an option of OrderStar specifying the thickness of the outline
of the symbols used to represent poles and zeros."

RungeKuttaLinearStabilityFunction::usage =
"RungeKuttaLinearStabilityFunction[amat, bvec, z] gives the linear
stability function of the Runge Kutta method with coefficient matrix
amat and weight vector bvec in the variable z."


Unprotect[OrderStar, OrderStarInterpolation, OrderStarKind,
OrderStarLegend, OrderStarPoles, OrderStarSymbolSize,
OrderStarSymbolThickness, OrderStarZeros, RungeKuttaLinearStabilityFunction];


(* Set default options. *)

Options[OrderStar] =
  {Axes->True, AxesOrigin->{0,0}, ColorFunction->Function[GrayLevel[1-#/2]],
   FrameTicks->None, PlotPoints->80, PlotRange->Automatic, Ticks->None,
   AspectRatio -> Automatic, OrderStarInterpolation->False,
   OrderStarKind->First, OrderStarLegend->False,
   OrderStarPoles->{True,True}, OrderStarZeros->{True,True},
   OrderStarSubPlots -> Automatic, OrderStarSymbolSize->.01,
   OrderStarSymbolThickness->.003
  };


Begin["`Private`"]

(**** Linear stability function for Runge-Kutta methods ****)

ScalarQ[_List]:= False;
ScalarQ[_]:= True;

RungeKuttaLinearStabilityFunction[a_?MatrixQ, b_?VectorQ, z_?ScalarQ]:=
  Module[{id, s},
    s = Length[b];
    id = IdentityMatrix[s];
    Det[ id - z*(a-Table[b,{s}]) ]/Det[ id - z*a ]
  ] /; Length[a]===Length[b];


(* Generic Error message. *)

OrderStar::opts = "The option `1` in OrderStar did not evaluate to `2`";

TFLQ[opt_] := (opt === True || opt === False || ListQ[opt]);
TFLString = "a pair consisting of True, False, or a List of {x,y} coordinates.";

opttest[bool_, mess_] := 
	If[bool, True, Message[OrderStar::opts, Apply[Sequence, mess]]; False];

optmessages = {
	{"OrderStarInterpolation", TFLString},
	{"OrderStarKind", "First or Second."},
	{"OrderStarLegend", "True, False, or a List of scaled \
		coordinates between 0 and 1."},
	{"PlotRange", "Automatic, or a list {{xmin,xmax},{ymin,ymax}}."},
	{"PlotPoints", "an integer>=2 or a list of two such integers."},
	{"OrderStarSymbolSize", "a positive number."},
	{"OrderStarSymbolThickness", "a positive number."},
	{"OrderStarZeros", TFLString},
	{"OrderStarPoles", TFLString},
	{"OrderStarSubPlots","a list of ContourGraphics and/or pairs of \
		PlotRange and PlotPoint option rules."}
  };

OptionTest[opts___]:= 
  Module[{datatypes, optlist, subplots},
    subplots = OrderStarSubPlots /. {opts};
    If[subplots === OrderStarSubPlots,
      subplots = OrderStarSubPlots /. Options[OrderStar]];
    If[MatchQ[subplots, {__Rule}], subplots = {subplots}];
    If[subplots === {}, subplots = False];
    optlist = {
      OrderStarInterpolation, OrderStarKind, OrderStarLegend,
      PlotRange, PlotPoints, OrderStarSymbolSize,
      OrderStarSymbolThickness, OrderStarZeros, OrderStarPoles} /.
        {opts} /. Options[OrderStar];
    AppendTo[optlist, subplots];
    datatypes = {
      TFLQ[optlist[[1]]], (* interpolation *)
      MemberQ[{First, Second}, optlist[[2]]], (* kind *)
      TFLQ[optlist[[3]]] && (* legendcoords *)
        If[ListQ[optlist[[3]]],
          (Apply[And,Map[(0<=#<=1)&,Flatten[optlist[[3]]]]]),
	      True],
      MatchQ[ Map[Union,N[optlist[[4]]]], (* plotrange *)
        Automatic|{{_?NumberQ,_?NumberQ},{_?NumberQ,_?NumberQ}}],
      MatchQ[optlist[[5]],_?(IntegerQ[#]&&#>=2&) | (* plotpoints *)
        {_?(IntegerQ[#]&&#>=2&),_?(IntegerQ[#]&&#>=2&)}],
      TrueQ[Positive[optlist[[6]]]], (* symbolsize *)
      TrueQ[Positive[optlist[[7]]]], (* symbolthickness *)
      MatchQ[optlist[[8]],{_?TFLQ,_?TFLQ}], (* zeros *)
      MatchQ[optlist[[9]],{_?TFLQ,_?TFLQ}], (* poles *)
      optlist[[10]]===Automatic || (* subplots *)
        optlist[[10]]===False ||
          MatchQ[
            optlist[[10]],
            {(_ContourGraphics |
              {PlotRange -> {{_,_},{_,_}}} |
              {PlotRange -> {{_,_},{_,_}},
                 PlotPoints -> (_Integer | {_Integer, _Integer})})...}]
    };
    If[And @@ MapThread[opttest,{datatypes,optmessages}], optlist, $Failed]
  ]; (* End of OptionTest. *)


(* Valid variables are non-numeric symbols and integer indexed functions (arrays). *)

OrderStar::var = "The expressions `1` and `2` are not univariate functions \
  of the same variable.";

varQ[_?NumericQ] = False;
varQ[_Symbol] = True;
varQ[_[(_Integer)..]] = True;
varQ[_] = False;

findz[f_,g_] :=
  Module[{nf},
    nf =
      Select[
        Union[ Join[ Level[{f,g},{-1}], Level[{f,g},{-2}] ] ],
        varQ
      ];
    If[Length[nf]==1,
      nf[[1]],
      Message[OrderStar::var, f, g]; $Failed
    ]
  ];


(* makeCP[ ] draws the main plot and sub plots *)

makeCP[f_, {plotrange_, plotpoints_}, plotoptions___] :=
  Module[{x, y},
    ContourPlot[f[x + I y],
      Evaluate[Prepend[plotrange[[1]], x]],
      Evaluate[Prepend[plotrange[[2]], y]],
      DisplayFunction :> Identity, 
      Evaluate[PlotPoints -> plotpoints],
      Evaluate[plotoptions]
    ]
  ];


(* automaticPR[ ] extracts a plot range from the list of zeros and poles
 using a scaling factor and an offset. The range of the plot aspect ratio
 is also restricted. *)

automaticPR[points_List] :=
  Module[{armax, scale, offset, xdiff, xmid, xmin, xmax, xrange,
      ydiff, ymid, ymin, ymax, yrange},
    {armax,scale,offset} = {2.5, 1.3, 1.};
    {xrange,yrange} = Thread[points];
    xmin = Min[xrange]; xmax = Max[xrange];
    ymin = Min[yrange]; ymax = Max[yrange];
    xmid = (xmin+xmax)/2.; ymid = (ymin+ymax)/2.;
    xdiff = Abs[xmax-xmin]/2.; ydiff = Abs[ymax-ymin]/2.;
    xdiff = Max[offset + scale xdiff, ydiff/armax];
    ydiff = Max[offset + scale ydiff, xdiff/armax];
    {xmid + {-xdiff, xdiff}, ymid + {-ydiff, ydiff}}
  ];

OrderStar::cvar = "The variables `1` and `2` in the expressions \
	`3` and `4` are not the same.";

OrderStar[R_, F_, opts___?OptionQ] :=
  Module[{ans, var = findz[R,F]},
    ans /; (var =!= $Failed &&
            (ans = OrderStar[R, F, var, opts];
             True))
  ];

OrderStar[R_, F_, var:(_Symbol | _Symbol[(_?IntegerQ)..]), opts___?OptionQ] :=
  Module[{ans},
    ans /; ((ans = OptionTest[opts]) =!= $Failed &&
            (ans = mainOrderStar[R, F, var, ans, opts];
             True))
  ];

(* Start of main routine. *)

mainOrderStar[R_, F_, var_, optslist_List, opts___?OptionQ] :=
  Module[{aspect, comp, contour, funct, glfontinfo, groptions, orderstar,
      plotinterp, plotlegend, plotpolesf, plotpolesr, plotzerosf,plotzerosr,
      plotoptions, plotrange, substars, subcgstars, interpolation,
      kind, legendcoords, range, ppoints, symbdata, zerosf, zerosr,
      polesf, polesr, subplots, symb},

    {interpolation, kind, legendcoords, range, ppoints, symbdata,
     symbdata, {zerosr, zerosf}, {polesr, polesf}, subplots} = optslist;


    plotoptions = Flatten[{opts,Contours->{contour},Options[OrderStar]}];
    plotoptions = DeleteCases[plotoptions,(PlotRange | PlotPoints) -> _, 1];
    plotoptions = FilterOptions[ContourPlot, Sequence @@ plotoptions];
    {aspect,comp} = {AspectRatio,Compiled} /. {plotoptions} /. Options[ContourPlot];

    groptions = FilterOptions[Graphics, plotoptions];

(* Get any font information for the legend window. The head of the pattern used
 * could be either Rule or RuleDelayed *)

	glfontinfo =
		Apply[
			Sequence,
			Cases[{groptions}, _[DefaultFont,_] | _[FormatType,_] | _[TextStyle,_]]
		];

(* Order star of first or second kind as a function of a symbol. *)

    funct =
      If[kind===First,
        contour = 1; Abs[R/F],
        contour = 0; Re[R-F]
      ] /. var->symb;

(* Compiled function definition for efficient graphics rendering. *)

    funct =
      If[TrueQ[comp],
        Compile[Evaluate[{{symb,_Complex}}, funct, {{_, _Complex}}] ],
        Function[Evaluate[{symb}, funct] ]
      ];

(* Used to calculate automatic plot range. *)

    plotinterp = findsolution[Numerator[R] - F Denominator[R], var,
			"interpolation points", "approximant", interpolation];
    plotpolesf = findsolution[1/F, var, "poles", "function", polesf];
    plotpolesr = findsolution[Denominator[R], var, "poles", "approximant", polesr];
    plotzerosf = findsolution[F, var, "zeros", "function", zerosf];
    plotzerosr = findsolution[Numerator[R], var, "zeros", "approximant", zerosr];


(* Calculate plot range. *)

    plotrange = Union[plotinterp, plotpolesf, plotpolesr, plotzerosf,
			plotzerosr, {{0, 0}}, SameTest -> Equal];

    plotrange = 
      If[range===Automatic,
        If[plotrange == {{0, 0}},
          {{-10,10},{-10,10}},
        (* else *)
          automaticPR[plotrange]
        ],
      (* else *)
        range
      ];

    aspect =
      Abs[N[#]]& @
        If[aspect === Automatic,
          (Subtract @@ plotrange[[2]])/(Subtract @@ plotrange[[1]]),
        (* else *)
          aspect
        ];


(* symbdata is {size, thickness, aspect} *)

    symbdata = {optslist[[6]], optslist[[7]], aspect};

    plotinterp = makeshape[plotinterp,symbdata,interpsymbol];
    plotpolesf = makeshape[plotpolesf,symbdata,polefsymbol];
    plotpolesr = makeshape[plotpolesr,symbdata,polersymbol];
    plotzerosf = makeshape[plotzerosf,symbdata,zerofsymbol];
    plotzerosr = makeshape[plotzerosr,symbdata,zerorsymbol];


(* Information window for the symbols used. *)

    plotlegend =
      If[legendcoords,
        legendwindow[{{0.01,1-0.25/aspect}, {0.35,1-0.01/aspect}},
          symbdata, glfontinfo ],
        {},
        legendwindow[legendcoords, symbdata, glfontinfo ]
      ];


(* Make the main plot. *)

    orderstar = makeCP[funct, {plotrange, ppoints}, plotoptions];


(* Make the subplots. *)

    ppoints = Min[ Max[ 15, Round[.5 ppoints] ], 50 ];

    If[subplots===Automatic,
      subplots = {{PlotRange -> .1 plotrange,
                   PlotPoints -> ppoints}}
    ];

(* Select ContourGraphics or PlotRange and PlotPoints rules. *)

	subcgstars = Cases[subplots,_ContourGraphics];
	subplots = DeleteCases[subplots,_ContourGraphics];

    substars =
      If[MatchQ[subplots,{__}], (* Non-empty list *)
        Map[makeCP[funct, #, Axes->False, plotoptions]&,
          {PlotRange, PlotPoints} /. subplots /. PlotPoints -> ppoints
        ],
      (* else *)
        {}
      ];


(* graphics for the symbols *)

    symbdata = {plotinterp,plotpolesf,plotpolesr,plotzerosf,plotzerosr};
    symbdata = Map[Graphics, Flatten[Thread[symbdata]]];


(* Combine the Graphics rendering symbol outline after the
 symbol background shape.  *)

    Show[Flatten[{orderstar, substars, subcgstars, symbdata, plotlegend}],
      AxesFront -> True, PlotRange -> plotrange, groptions,
      DisplayFunction:>$DisplayFunction ]

  ]; (* End of mainOrderStar. *)


(* Define functions used in OrderStar. *)

(* Generate a list of solution points. *)

(* Avoid possible division by zero when poles not required. *)

SetAttributes[findsolution,HoldFirst];

findsolution[__,False] = {};

findsolution[eqn_, var_, info_, func_, True]:=
  extractsolutions[eqn,var,info,func];

findsolution[eqn_, var_, info_, func_, points_List]:=
  Union[
    Join[
      extractsolutions[eqn,var,info,func],
      SetPrecision[points,6] /. 0 -> 0.0
    ],
    SameTest -> Equal
  ];

(* Remove infinite solutions and generate solution messages. *)

OrderStar::sols = "Warning: No `1` of `2` found using NSolve. Either inverse
functions or transcendental dependencies were involved. Try specifying omitted
points using options.";

(* No solutions *)

finitesolutions[{},_,_,False]:= {};

(* No solutions, but Solve used inverse fuctions etc *)

finitesolutions[{},info_,func_,True]:=
  (Message[OrderStar::sols,info,func]; {});

(* Finite solutions *)

finitesolrules = {Complex[x_,y_]->{x,y}, x_?NumberQ->{x,0}};

finitesolutions[solutions_,info_,func_,False]:=
  Module[{fsols},
    fsols = DeleteCases[solutions,_DirectedInfinity];
    fsols /. finitesolrules
  ];

(* Generate a message if there were no finite solutions, but Solve
 used inverse fuctions etc *)

finitesolutions[solutions_,info_,func_,True]:=
  Module[{fsols},
    fsols = DeleteCases[solutions,_DirectedInfinity];
    If[fsols === {},
      Message[OrderStar::sols,"finite "<>info,func]
    ];
    fsols /. finitesolrules
  ];

(* NSolve eqn in terms of var. Suppress Solve messages, but
 set a flag if messages were generated. *)

extractsolutions[eqn_,var_,info_,func_]:=
  Module[{sol, msgs=False},
    Block[{$Messages},
		Check[ sol = NSolve[eqn==0, var], msgs = True; sol]
    ];
    sol = If[MatchQ[sol,_NSolve|{}|{{}}], {}, var /. sol];
    SetPrecision[finitesolutions[sol,info,func,msgs], 6] /. 0 -> 0.0
  ];


(* General graphics primitive for the symbfuncs
	interpsymbol, polefsymbol, polersymbol, zerofsymbol,
	and zerorsymbol. *)

makeshape[{},__]:= {{},{}};

makeshape[coords_List, symbdata_, symbfunc_]:= 
  Apply[Sequence,Map[symbfunc[#,symbdata]&,coords]];

makeshape[Scaled[coords_], symbdata_, symbfunc_]:=
  symbfunc[coords, symbdata];


(* Graphics primitives for symbols with graphicsprims Line and Polygon. *)

interpsymbol[coords_, {size_, thick_, ar_}]:=
  {{GrayLevel[1], Disk[coords,Scaled[{size, size/ar}]]},
   {GrayLevel[0], Thickness[thick], Circle[coords,Scaled[{size, size/ar}]]}};

square[coords_, {size_, _, ar_}, graphicsprim_]:=
  graphicsprim[Scaled[#, coords]& /@
    (size {{-1,-1/ar},{1,-1/ar},{1,1/ar},{-1,1/ar},{-1,-1/ar}})];

diamond[coords_, {size_, _, ar_}, graphicsprim_]:=
  graphicsprim[Scaled[#, coords]& /@
    (size {{-1.4,0},{0,-1.4/ar},{1.4,0},{0,1.4/ar},{-1.4,0}})];

polefsymbol[coords_, {size_, thick_, ar_}]:=
  {{}, {GrayLevel[0], Thickness[thick],
        Line[{Scaled[{-size,size/ar},coords],
          Scaled[{size,-size/ar},coords]}],
        Line[{Scaled[{-size,-size/ar},coords],
          Scaled[{size,size/ar},coords]}]}};

polersymbol[coords_, sd:{_,thick_,_}]:=
  {{GrayLevel[1], diamond[coords,sd,Polygon]},
   {Thickness[thick], GrayLevel[0], diamond[coords,sd,Line]}};

zerofsymbol[coords_, {size_, thick_, ar_}]:=
  {{}, {GrayLevel[0], Thickness[thick],
        Line[{Scaled[{0,-size/ar},coords],
          Scaled[{0, size/ar},coords]}],
        Line[{Scaled[{-size,0},coords],
          Scaled[{ size,0},coords]}]}};

zerorsymbol[coords_, sd:{_,thick_,_}]:=
  {{GrayLevel[1], square[coords,sd,Polygon]},
   {Thickness[thick], GrayLevel[0], square[coords,sd,Line]}};



(* Primitives for symbol information window. *)

showsymbols[sd:{size_, _, ar_}, opts___]:=
  With[{symbpos = 2. size, textposn = 6. size},
    Graphics[
      Flatten[{
        makeshape[Scaled[{symbpos, 0.9}],sd,polersymbol ],
       Text[" Poles of approximant",Scaled[{textposn, 0.9}],{-1,0}],
       makeshape[Scaled[{symbpos, 0.7}],sd,zerorsymbol ],
       Text[" Zeros of approximant",Scaled[{textposn, 0.7}],{-1,0}],
       makeshape[Scaled[{symbpos, 0.5}],sd,polefsymbol],
       Text[" Poles of function",Scaled[{textposn, 0.5}],{-1,0}],
       makeshape[Scaled[{symbpos, 0.3}],sd,zerofsymbol],
       Text[" Zeros of function",Scaled[{textposn, 0.3}],{-1,0}],
       makeshape[Scaled[{symbpos, 0.1}],sd,interpsymbol],
       Text[" Interpolation points",Scaled[{textposn, 0.1}],{-1,0}]
       }], AspectRatio -> ar, opts
     ]
   ];

legendwindow[{pt1_List, pt2_List}, {size_, thick_, ar_}, opts___]:=
  Module[{pt1s = Scaled[pt1], pt2s = Scaled[pt2], sd, xr, yr},
      {xr, yr} = 1/Abs[pt2-pt1];
    sd = {xr size, xr thick, ar xr/yr};
    Graphics[
      {GrayLevel[1], Rectangle[pt1s,pt2s],
       Thickness[thick], GrayLevel[0],
        Line[{pt1s, Scaled[{pt1[[1]], pt2[[2]]}],
            pt2s, Scaled[{pt2[[1]], pt1[[2]]}], pt1s}],
        Rectangle[pt1s, pt2s, showsymbols[sd, opts]]}
    ]
  ];

End[];    (* End `Private` Context. *)

(* Protect exported symbols. *)

SetAttributes[
{OrderStar, RungeKuttaLinearStabilityFunction},
ReadProtected];

Protect[OrderStar, OrderStarInterpolation, OrderStarKind,
OrderStarLegend, OrderStarPoles, OrderStarSymbolSize,
OrderStarSymbolThickness, OrderStarZeros,
RungeKuttaLinearStabilityFunction];

EndPackage[];    (* End package Context. *)
