(* :Title:FazPort *)
(* :Author: Miroslava Dubcova, VSCHT *)

(* :Summary: 
      Fazove portrety pro soustavu dvou obycejnych diferencialnich rovnic 

*)
(* :Context: Matematika`FazPort` *)
(* :Package Version: 1.0 *)
(* :Mathematica Version: 2.2 3.0 *)
(* :Copyright: Copyright Miroslava Dubcova *)
(* :History:
*)
(* :Keywords:
   LFazPort, sol, NFazPort, nsol
*)
(* :Warning:  *)

BeginPackage["Matematika`FazPort`","Graphics`Arrow`"]
LFazPort::usage=
      "LFazPort[t1,t2,a11,a12,a21,a22,param] Fazovy portrety soustavy obycejnych 
diferencialnich rovnic x'=Ax. t1 je pocatecni cas integrace, t2 je koncovy 
cas integrace, aij jsou prvky matice A. Param je nepovinny parametr pro vykresleni grafu. "

NFazPort::usage=
      "NFazPort[f1,f2,t1,a1,a2,b1,b2,data,param] Fazovy portrety soustavy obycejnych 
diferencialnich rovnic x'=f(x). f1 a f2 jsou slozky prave strany, t1 je koncovy 
cas integrace, a1, a2, b1, b2 jsou meze grafu a data jsou pocatecni podminky ve tvaru
{{x0,y0},{x1,y1}...}.param je nepovinny parametr pro funkci ParametricPlot."


Begin["`Private`"]
sol[t_,a_,b_,a11_,a12_,a21_,
    a22_]:={x[t],
      y[t]}/.DSolve[{x'[t]==a11* x[t]+a12* y[t],y'[t]==a21* x[t]+a22 *y[t],
        x[0]==a,y[0]==b},{x[t],y[t]},t];
LFazPort[t1_,t2_,a11_,a12_,a21_,a22_,param___]:=Module[{a1,a2,range,asp,alfa},
    alfa=Eigenvalues[{{a11,a12},{a21,a22}}];
    vec=Eigenvectors[{{a11,a12},{a21,a22}}];
    Options[LFazPort]={PlotRange->{{-2,2},{-2,2}},AspectRatio->Automatic};
    range=PlotRange /. {param} /. Options[LFazPort];
    asp=AspectRatio /. {param} /. Options[LFazPort];
   a1=N[alfa[[1]]];a2=N[alfa[[2]]];
   Print["Vlastni cisla:   ",alfa[[1]]," , ",alfa[[2]]];
   Print["vlastni vektory:   ",vec[[1]]," , ",vec[[2]]];
  If[a1==0 || a2==0,Print["Vlastni cislo je nulove"],
   If[Re[alfa[[1]]]<0 && Re[alfa[[2]]]<0,
  ParametricPlot[
    Evaluate[Join[Flatten[Table[sol[z,a,b,a11,a12,a21,a22],{a,-2,2,4},{b,-2,2,1}],2],
      Flatten[Table[sol[z,a,b,a11,a12,a21,a22],{a,-2,2,1},{b,-2,2,4}],2]],
      {z,t1,t2},AspectRatio->asp,PlotRange->range,
Epilog->{Table[Arrow[{a,b},{a+0.5*(a11*a+a12*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2],
b+0.5*(a21*a+a22*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2]},HeadCenter->0],{a,-2,-1},{b,-2,2}],
Table[Arrow[{a,b},{a+0.5*(a11*a+a12*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2],
b+0.5*(a21*a+a22*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2]},HeadCenter->0],{a,1,2},{b,-2,2}]},param]],
    If[Re[a1]>0 && Re[a2]>0,
  ParametricPlot[
    Evaluate[Join[Flatten[Table[sol[z,a,b,a11,a12,a21,a22],{a,-2,2,4},{b,-2,2,1}],2],
      Flatten[Table[sol[z,a,b,a11,a12,a21,a22],{a,-2,2,1},{b,-2,2,4}],2]],
      {z,t1,-t2},AspectRatio->asp,PlotRange->range,
Epilog->{Table[Arrow[{a,b},{a+0.5*(a11*a+a12*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2],
b+0.5*(a21*a+a22*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2]},HeadCenter->0],{a,-2,-1},{b,-2,2}],
Table[Arrow[{a,b},{a+0.5*(a11*a+a12*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2],
b+0.5*(a21*a+a22*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2]},HeadCenter->0],{a,1,2},{b,-2,2}]},param]],
(*  ParametricPlot[
    Evaluate[Flatten[Table[sol[z,0.1*Cos[2Pi*a],0.1*Sin[2Pi*a],a11,a12,a21,a22],
            {a,0,1,0.05}],1],
      {z,t1,t2},AspectRatio->asp,PlotRange->range,
Epilog->Table[Arrow[{a,b},{a+0.5*(a11*a+a12*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2],
b+0.5*(a21*a+a22*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2]},HeadCenter->0],{a,-1,1,2},{b,-1,1,2}],param]],*)
      If[Re[a1]==0,ParametricPlot[
        Evaluate[Join[Flatten[Table[sol[z,a,a,a11,a12,a21,a22],{a,-2,2,0.5}],1],
         Flatten[Table[sol[z,a,-a,a11,a12,a21,a22],{a,-2,2,0.5}],1]
          ],
          {z,t1,t2},AspectRatio->asp,PlotRange->range,
Epilog->{Table[Arrow[{a,b},{a+0.5*(a11*a+a12*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2],
b+0.5*(a21*a+a22*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2]},HeadCenter->0],{a,-2,-1},{b,-2,2}],
Table[Arrow[{a,b},{a+0.5*(a11*a+a12*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2],
b+0.5*(a21*a+a22*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2]},HeadCenter->0],{a,1,2},{b,-2,2}]},param]],
    ParametricPlot[
      Evaluate[Join[Flatten[Table[sol[z,a,b,a11,a12,a21,a22],{a,-2,2,4},{b,-2,2,0.4}],2],
        Flatten[Table[sol[z,a,b,a11,a12,a21,a22],{a,-2,2,0.4},{b,-2,2,4}],2],
          sol[z,vec[[2,1]]*0.01,vec[[2,2]]*0.01,a11,a12,a21,a22],
           sol[z,-vec[[2,1]]*0.01,-vec[[2,2]]*0.01,a11,a12,a21,a22],
            sol[z,2*vec[[1,1]]/vec[[1,2]],2,a11,a12,a21,a22],
             sol[z,-2*vec[[1,1]]/vec[[1,2]],-2,a11,a12,a21,a22]],
        {z,t1,t2},AspectRatio->asp,PlotRange->range,
Epilog->{Table[Arrow[{a,b},{a+0.5*(a11*a+a12*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2],
b+0.5*(a21*a+a22*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2]},HeadCenter->0],{a,-2,-1},{b,-2,2}],
Table[Arrow[{a,b},{a+0.5*(a11*a+a12*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2],
b+0.5*(a21*a+a22*b)/Sqrt[(a11*a+a12*b)^2+(a21*a+a22*b)^2]},HeadCenter->0],{a,1,2},{b,-2,2}]},param]]]]
]]]
nsol[t_,f_,g_,a_,b_,
    t0_]:={x[t],
      y[t]}/.NDSolve[{x'[t]==f[x[t],y[t]],y'[t]==g[x[t],y[t]],x[0]==a,
        y[0]==b},{x[t],y[t]},{t,0,t0}]
NFazPort[f1_,f2_,t1_,a1_,a2_,b1_,b2_,data_,param___]:=Module[{},{
  ParametricPlot[
    Evaluate[Flatten[Apply[nsol,Apply[{z,f1,f2,#1,#2,t1}&,data,-1],{1}],1]],{z,0,
      t1},AspectRatio->(b2-b1)/(a2-a1),PlotRange->{{a1,a2},{b1,b2}},param]}]
End[]
EndPackage[]