(* ::Package:: *)

(************************************************************************)
(* 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.                                                         *)
(************************************************************************)



Needs["Quaternions`"]
HJtensorProduct[x_,y_]:= Module[
{QZero=Quaternion[0,0,0,0],Q,QTmp,X=Apply[Quaternion,x],Y=Apply[Quaternion,y]},
Table[
Q=QZero; Q[[i]]=1.0;
QTmp=QZero; QTmp[[i]]=1.0;
QTmp = X**QTmp;
QTmp=QTmp**Conjugate[Y];
Apply[List,QTmp],
{i,1,4}
]
]
HamiltonJacobi4Sym[S_]:=Module[
{T=ConstantArray[0,{3,3}],u,\[Sigma],v},
T[[1,1]]=0.25*(S[[1,1]]+S[[2,2]]-S[[3,3]]-S[[4,4]]);
T[[2,1]]=0.5*(S[[2,3]]+S[[1,4]]);
T[[3,1]]=0.5*(S[[2,4]]-S[[1,3]]);
T[[1,2]]=0.5*(S[[2,3]]-S[[1,4]]);
T[[2,2]]=0.25*(S[[1,1]]-S[[2,2]]+S[[3,3]]-S[[4,4]]);
T[[3,2]]=0.5*(S[[3,4]]+S[[1,2]]);
T[[1,3]]=0.5*(S[[2,4]]+S[[1,3]]);
T[[2,3]]=0.5*(S[[3,4]]-S[[1,2]]);
T[[3,3]]=0.25*(S[[1,1]]-S[[2,2]]-S[[3,3]]+S[[4,4]]);
{u,\[Sigma],v}=SingularValueDecomposition[T,1];{u,\[Sigma],v}=Map[Flatten,{u,\[Sigma],v}];
u=Prepend[u,0];
v=Prepend[v,0];
u[[2]]+=1;
v[[2]]+=1;
u=u/Norm[u];
v=v/Norm[v];
HJtensorProduct[u,v]
]

Clear[HamiltonJacobi4Skw]
HamiltonJacobi4Skw[K_]:=Module[
{p,q,i=Quaternion[0,1,0,0],pQ,qQ,Np,Nq},
(* The following are Niloufer M 6.6 and 6.7 *)
p=0.5 {-(K[[1,2]]+K[[3,4]]),(K[[2,4]]-K[[1,3]]),-(K[[1,4]]+K[[2,3]])};
q= 0.5{   K[[1,2]]-K[[3,4]],     K[[2,4]]+K[[1,3]],      K[[1,4]]-K[[2,3]]};
{pQ,qQ}={Apply[Quaternion,Join[{0},p]],Apply[Quaternion,Join[{0},q]]};
(* The stuff below does not work *)
(* Not sure why but it almost has to be in the Tensor Product*)
(*{p,q}={Apply[Quaternion,p],Apply[Quaternion,q]};
{x,y}={(Abs[p]- i**p)/Abs[Abs[p]-i**p],(Abs[q]- i **q)/Abs[Abs[q]-i **q]};
{x,y}={Apply[List,x],Apply[List,y]};
HJtensorProduct[x,y]*)
(* direct implementation *)
{Np,Nq}=Map[Norm,{p,q}];
Qp=1/Abs[Abs[pQ]-i**pQ] ({
 {Np+p[[1]], 0, -p[[3]], p[[2]]},
 {0, Np+p[[1]], p[[2]], p[[3]]},
 {p[[3]], -p[[2]], Np+p[[1]], 0},
 {-p[[2]], -p[[3]], 0, Np+p[[1]]}
});
Qq=1/Abs[Abs[qQ]-i**qQ] ({
 {Nq+q[[1]], 0, q[[3]], -q[[2]]},
 {0, Nq+q[[1]], q[[2]], q[[3]]},
 {-q[[3]], -q[[2]], Nq+q[[1]], 0},
 {q[[2]], -q[[3]], 0, Nq+q[[1]]}
});
Qp.Qq
]

NonUpperTriangularity[A_List]:=Module[
{AUpper},
AUpper=UpperTriangularize[ A];
1-Norm[AUpper,"Frobenius"]/Norm[A,"Frobenius"]
]
NonUpperTriangularity[A_SparseArray]:=Module[
{AUpper},
AUpper=UpperTriangularize[ A];
1-Norm[AUpper,"Frobenius"]/Norm[A,"Frobenius"]
]


Needs["Quaternions`"]
HJtensorProduct[x_,y_]:= Module[
{QZero=Quaternion[0,0,0,0],Q,QTmp,X=Apply[Quaternion,x],Y=Apply[Quaternion,y]},
Table[
Q=QZero; Q[[i]]=1.0;
QTmp=QZero; QTmp[[i]]=1.0;
QTmp = X**QTmp;
QTmp=QTmp**Conjugate[Y];
Apply[List,QTmp],
{i,1,4}
]
]


HamiltonJacobi4Sym[S_]:=Module[
{T=ConstantArray[0,{3,3}],u,\[Sigma],v},
T[[1,1]]=0.25*(S[[1,1]]+S[[2,2]]-S[[3,3]]-S[[4,4]]);
T[[2,1]]=0.5*(S[[2,3]]+S[[1,4]]);
T[[3,1]]=0.5*(S[[2,4]]-S[[1,3]]);
T[[1,2]]=0.5*(S[[2,3]]-S[[1,4]]);
T[[2,2]]=0.25*(S[[1,1]]-S[[2,2]]+S[[3,3]]-S[[4,4]]);
T[[3,2]]=0.5*(S[[3,4]]+S[[1,2]]);
T[[1,3]]=0.5*(S[[2,4]]+S[[1,3]]);
T[[2,3]]=0.5*(S[[3,4]]-S[[1,2]]);
T[[3,3]]=0.25*(S[[1,1]]-S[[2,2]]-S[[3,3]]+S[[4,4]]);
TOut=T;
{u,\[Sigma],v}=SingularValueDecomposition[T,1];
{u,\[Sigma],v}=Map[Flatten,{u,\[Sigma],v}];
u=Prepend[u,0];
v=Prepend[v,0];
u[[2]]+=1;
v[[2]]+=1;
u=u/Norm[u];
v=v/Norm[v];
HJtensorProduct[u,v]
]



