(*******************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which
typically will have the same name as this file except ending in
".nb" instead of ".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent
to using the Evaluate Initialization Cells menu command in the front
end.

DO NOT EDIT THIS FILE.  This entire file is regenerated
automatically each time the parent Notebook file is saved in the
Mathematica front end.  Any changes you make to this file will be
overwritten.
***********************************************************************)









BeginPackage["BasicGeo`"]



Diff::usage="Compute the differential/gradient of a function \[Lambda] w.r.t. a set of variables\n
Diff[\[Lambda],variables]"
LieDerivative::usage="Compute the Lie-derivate of either a scalar \[Lambda] along a vector field f, or a co-vector field \[Omega] along a vector field f\n
  LieDerivative[\[Lambda],f,vars]\n
  LieDerivative[\[Omega],f,vars]"
JacobianMatrix::usage="Compute the Jacobian of a vector field\n
  JacobianMatrix[f,vars]"
LieProduct::usage="Compute the Lie-product between two vector fields\n
  LieProduct[f,g,vars]"
ChangeCoordinates::usage="Change coordinates in a control-affine non-linear system.\n
  x'=f(x)+g(x)u\n
  y=h(x)\n
  z=\[CapitalPhi](x)\n
  ChangeCoordinates[f,g,h,\[CapitalPhi],vars]\n
  \n
  The vectorfield g and function h are optional."
Flow::usage="Compute the flow of a vectorfield f, with independent variable t, starting at point x0\n
  Flow[f,t,x0,vars]\n
  Note: Highly experimental"
FrobeniusInverse::usage="Compute the inverse transformation obtained when following the flows in F\n
  FrobeniusInverse[F,x0,origvars,newvars]\n
  Note: Highly experimental"
FrobSolve::usage="Compute all solutions to the Frobenius PDE\n
  FrobSolve[F,x0,vars,d]\n
  Note: Highly experimental"
DistMemberQ::usage=
  "Determine if a vectorfield f is a member of a non-singular distribution \[CapitalDelta]=span{f1,...,fd}. DistMemberQ[f,\[CapitalDelta]]"



















Begin["`Private`"]









Diff[\[Lambda]_,vars_]:=Map[D[\[Lambda],#]&,vars];
Diff[f_List,vars_]:=Map[Diff[#,vars]&,f];

JacobianMatrix[f_, vars_]:=Diff[f,vars];



LieDerivative[\[Lambda]_, f_,vars_] := Diff[\[Lambda], vars].f;
LieDerivative[w_List, f_,vars_] := 
    Simplify[{f}.Transpose[JacobianMatrix[w,vars]]+
          {w}.JacobianMatrix[f,vars]][[1]];



LieProduct[f_,g_, vars_]:=
    JacobianMatrix[g, vars].f - JacobianMatrix[f, vars].g;



ChangeCoordinates[f_,\[CapitalPhi]_,oldvar_]:=Module[{newsys},
      newsys=ChangeCoordinates[f,0*f,0,\[CapitalPhi],oldvar];
      Return[newsys[[1]]]];

ChangeCoordinates[f_,g_,\[CapitalPhi]_,oldvar_]:=Module[{newsys},
      newsys=ChangeCoordinates[f,g,0,\[CapitalPhi],oldvar];
      Return[newsys[[{1,2}]]]];

ChangeCoordinates[f_,g_,h_,\[CapitalPhi]_,oldvar_]:=
    Module[{J,\[CapitalPhi]inv,f2,g2,h2,newvar},
      newvar=\[CapitalPhi]/.{Rule[x_,y_]\[Rule]x};
      
      (* Compute Jacobian *)
      
      J=JacobianMatrix[newvar/.\[CapitalPhi],oldvar];
      
      (* Invert transformation. 
            Assume unique inverse exists *)
      \[CapitalPhi]inv =
        Solve[(#\[Equal](#/.\[CapitalPhi]))&/@newvar,
            oldvar][[1]];
      
      (* Convert f *)
      f2=Simplify[J.f/.\[CapitalPhi]inv];
      
      (* Convert g. Map over all g:s if more than one *)
      
      g2=If[Head[g[[1]]]===List,
          Map[Simplify[J.#/.\[CapitalPhi]inv]&,g],
          Simplify[J.g/.\[CapitalPhi]inv]];
      
      (* Convert h. Map over all h:s if more than one *)
      
      h2=If[Head[h]===List,
          Map[Simplify[#/.\[CapitalPhi]inv]&,h],
          Simplify[h/.\[CapitalPhi]inv]];
      
      Return[{f2,g2,h2}];
      ];



DistMemberQ[f_,\[CapitalDelta]_]:=Module[{n,d,z,rhs,eq},
      d=Length[\[CapitalDelta]];
      n=Length[f];
      consts=(C[#]&/@Range[1,d]);
      rhs=consts.\[CapitalDelta];
      eq=(f[[#]]\[Equal]rhs[[#]]&/@Range[1,n]);
      Return[Length[Solve[eq,consts]]>0]];





Flow[f_,z_,x_,vars_]:=
    DSolveWithX0[CreateFlowDiff[f,vars],
        x,vars]/.{t\[Rule]z};

FindConstants[sol_,x0_]:=Module[{solconstants},
      (* Hitta vilka konstantern fr\[ARing]n ODE-
          l\[ODoubleDot]saren som finns med *)
      
      solconstants=Select[Variables[sol],MatchQ[#,C[n_Integer]]&];
      
      (* Invertera nu systemet och skicka tillbaka resultatet *)
      
      Return[Solve[x0[[#]]\[Equal]sol[[#]]&/@Range[1,Length[x0]],
            solconstants][[1]]];
      ];

CreateFlowDiff[f_,vars_]:=Module[{totvars,primevars},
      totvars=#\[Rule]#[t]&/@vars;
      primevars=#'[t]&/@vars;
      Return[
        primevars[[#]]\[Equal](f[[#]]/.totvars)&/@Range[1,Length[vars]]];
      ];

DSolveWithX0[diffeq_,x0_,vars_]:=Module[{totvars,diffsol,consts},
      totvars=#\[Rule]#[t]&/@vars;
      (*Assume a solution exists, grab first*)
      
      diffsol=(vars/.totvars)/.(DSolve[diffeq,vars/.totvars,t][[1]]);
      (* Compute constants such that x(0)=x0 *)
      
      consts=FindConstants[diffsol/.{t\[Rule]0},x0];
      Return[diffsol/.consts];
      ];



FrobeniusInverse[F_,x0_,origvars_,newvars_]:=
    Module[{x,n},
      n=Length[origvars];
      x=x0;
      For[k=n,k>0,k--,
        x=
            Flow[F[[k]],newvars[[k]],origvars,
                origvars]/.(origvars[[#]]\[Rule]x[[#]]&/@Range[1,n]);
        ];
      Return[Simplify[x]];
      ];



Global`FrobSilent=True;

FrobSolve[F_,x0_,vars_]:=FrobSolve[F,x0,vars,0];
FrobSolve[F_,x0_,vars_,d_]:=
    Module[{n,foovars,\[CapitalPsi],\[CapitalPsi]inv},
      If[Global`FrobSilent,
        Off[Solve::incnst,Solve::ifun,InverseFunction::ifun]];
      n=Length[F[[1]]];
      foovars=z[#]&/@Range[1,n];
      \[CapitalPsi]=FrobeniusInverse[F,x0,vars,foovars];
      \[CapitalPsi]inv=
        Solve[vars[[#]]\[Equal]\[CapitalPsi][[#]]&/@Range[1,n],foovars][[1]];
      If[Global`FrobSilent,
        On[Solve::incnst,Solve::ifun,InverseFunction::ifun]];
      Return[Take[foovars/.\[CapitalPsi]inv,-(n-d)]];
      ];



End[];

EndPackage[];