(* Kinematics Generator *)
Print["by K.Kajda ver:1.2"]
Print["last modification 01 Feb 2007"]


BeginPackage["KinematicsGen`"]


Kinematics::usage="Usage: Kinematics[kinematics,choice], kinematics is a list of squared external momenta and choice sets independent Mandelstam variables e.g.\n kinematics={p1^2->m^2,p2^2->m^2,p3^2->m^2,p4^2->m^2}\n choice={s,t}\n\n NOTE 1: All external momenta are set as IN and clocwise, use ?s, ?u, ?t to see Mandelstam variables.\n\n NOTE 2: For 5pt function use:\n {s12, s23, s34, s45, s15}.\n\n NOTE 3: For 6pt function use:\n {s12, s23, s34, s45, s56, s16, s123, s234, s345}.\n\n NOTE 4: For 5 and 6 point functions independent kinematical variables must be chosen in a cyclic way. See examples."
s::usage="Mandelstam variable for 4-point function:\n s=(p1+p2)^2=(p3+p4)^2"
u::usage="Mandelstam variable for 4-point function:\n t=(p2+p3)^2=(p1+p4)^2"
t::usage="Mandelstam variable for 4-point function:\n u=(p1+p3)^2=(p2+p4)^2"
p1::usage="First external momenta"
p2::usage="Second external momenta"
p3::usage="Third external momenta"
p4::usage="Fourth external momenta"
p5::usage="Fifth external momenta"
p6::usage="Sixth external momenta"


Begin["`Private`"]


Kinematics[kin_,choice_]:=
  Module[{momlst, Lm, prodlst, grul, urul, 
	  EST, EST2, SLCheck, FLS, 
	  inv, g, scalarNEW, b2,
	  expr, Aux, equat, chlst, Sleft,
	  variab, chrules, leftprods, momsum, 
          sd, M2, M1, result},


	 (* functions *)

    EST[x_]:= 
      ToExpression[
      "p"<>StringTake[ToString[x],{2}]<>"*"<>
      "p"<>StringTake[ToString[x],{3}]];

    FLS[m_]:=
      Flatten[
       Solve[Plus @@ momlst==0 /.m->f[m],f[m]] /.f[x_]->x];

    EST2[x_]:= 
      ToExpression[
      "p"<>StringTake[ToString[x],{2}]<>"*"<> 
      "p"<>StringTake[ToString[x],{3}]<>"*"<> 
      "p"<>StringTake[ToString[x],{4}]];

    SLCheck[x_]:=StringLength[ToString[x]];

	 (* end functions *)

    momlst=kin/.(a_^_->_)->a;
    Lm=Length[momlst];

    prodlst=
      DeleteCases[
	Union[Flatten[
	 Outer[List,momlst,momlst]/.{a_,b_}->a*b
      ]],_^2];

    grul=#->g[#]& /@ prodlst;


(******************************************************)
(* 4-point function *)
If[Lm==4,
    {

    M2=prodlst 
	 /.{(p1*p2|p3*p4)->s,(p2*p3|p1*p4)->t,(p1*p3|p2*p4)->u};
    M1=Expand[#^2& /@ (prodlst/.Times->Plus)]/.kin/.grul;

  
    If[Length[choice]==2,
    {
      Aux=Plus @@ choice;
      equat=s+t+u-Aux==Plus @@ (#^2& /@ momlst)-Aux;
      urul=Flatten[Solve[equat,equat /.(x_== y_)->x] /.kin]
    },{urul={}}];

    result=Solve[M1==M2,g[#]& /@ prodlst]/.urul/. g[x_]->x;
    result=Join[kin,result]//Flatten;
    Goto[KinEnd]

    }];

(******************************************************)
(* 5-point function *)
If[Lm==5,
    {
 
    inv=Flatten[
      Solve[(Expand[(Plus @@ EST[#])^2] /.grul /.kin)==#,
        g[Times @@ EST[#]]]& /@ choice];
    scalarNEW=
      DeleteCases[prodlst,#|##2& @@ (Times @@ EST[#]& /@ choice)];
	       
 (* Prepare Equations *)

    M1=(scalarNEW /.a_*b_->Expand[(a+b)^2]) /.kin /.grul;
    b2=scalarNEW /. px_*py_->px+py;
    M2=(Expand[(#)^2/.FLS[#]]/.kin& /@ b2) /.grul /.inv;
    expr= #==0& /@ Inner[Plus,M2,-M1,List];

 (* End Prepare Equations *)

    result=
      Join[kin,inv,
       Flatten[Solve[expr,g /@ scalarNEW]]/.g[x_]->x]/.g[x_]->x;


    }];


(******************************************************)
(* 6-point function *)

(******************************************************)
If[Lm==6,
    {
   
    chlst=If[SLCheck[#]==3,EST[#],EST2[#]]& /@ choice;
    Sleft=(Expand[(# /.Times->Plus)^2]& /@ chlst) /.kin /.grul;
    variab=Cases[Flatten[Sleft /.Plus->List],a_.*g[x_]->g[x]];
    chrules=Solve[#==0& /@ (Sleft-choice),variab]// Flatten;
    leftprods=Complement[prodlst /.grul,variab,{}]/.g[x_]->x;
    momsum=#/.Times->Plus& /@ leftprods;


 (* Prepare Equations *)

    M1=(Expand[(p1+p2+p3+p4+p5+p6/.#->0)^2]& /@ momsum)
    /.kin/.grul/.chrules;

    M2=(Expand[(#/.Times->Plus)^2]& /@ leftprods)/.kin/.grul;

 (* End Prepare Equations*)

    sd=Solve[#==0& /@ (M1-M2),g[#]& /@ leftprods];
    result=Join[kin,chrules,sd]/. g[a_] -> a //Flatten;

    }
];


Label[KinEnd];

 Return[result]]

End[]

EndPackage[]