#! /bin/sh -x

me=`basename $0`
name=m/$me-$1

[ -d m ] || mkdir m

echo "$name" | grep -E '^m/(deltaMs|deltaMd)-(Hp|Cha|Neu|Glu|SM|All)$' || {
  echo "invoke me as deltaMs|deltaMd Hp|Cha|Neu|Glu|SM|All" ;
  exit 1
}

rm -f $name.*

math -run me=\"$me\" -run tag=\"$1\" -run name=\"$name\" << \_EOF_
AppendTo[$Echo, "stdout"];

<< FeynArts`

<< FormCalc`

MU = MU2 = 0;
MC = MC2 = 0;
MD = MD2 = 0;
MS = MS2 = 0;
MB = MB2 = 0;

fs = 2;
fd["deltaMd"] = 1;
fd["deltaMs"] = 3;

proc = {-F[4, {fs}], F[4, {fd[me]}]} -> {-F[4, {fd[me]}], F[4, {fs}]};

opt["All"] = Sequence[];
opt["Hp"]  = LastSelections -> S[5];
opt["Cha"] = LastSelections -> F[12];
opt["Neu"] = Sequence[ExcludeParticles -> F[15], LastSelections -> F[11]];
opt["Glu"] = LastSelections -> F[15];
opt["SM"]  = ExcludeParticles -> {F[11|12|15], S[5]};

Attributes[C0Z] = Attributes[C00Z] =
Attributes[D0Z] = Attributes[D00Z] = {Orderless};
C0i[LoopTools`cc0, 0, 0, 0, m__] := C0Z[m];
C0i[LoopTools`cc00, 0, 0, 0, m__] := C00Z[m];
D0i[LoopTools`dd0, 0, 0, 0, 0, 0, 0, m__] := D0Z[m];
D0i[LoopTools`dd00, 0, 0, 0, 0, 0, 0, m__] := D00Z[m];


tops = CreateTopologies[1, 2 -> 2, BoxesOnly];

ins = InsertFields[tops,
  proc, opt[tag],
  InsertionLevel -> {Particles},
  Model -> MSSMQCD];

Paint[ins,
  ColumnsXRows -> {4, 5},
  DisplayFunction -> (Export[name <> ".ps", #, ImageSize -> 72 2 {4, 5}]&)];

amp = CalcFeynAmp[CreateFeynAmp[ins],
  FermionChains -> Chiral,
  FermionOrder -> Colour[1, 4, 2, 3]];

amp = amp /.
  D0i[id_, __] :> 0 /; !MemberQ[{LoopTools`dd0, LoopTools`dd00}, id] /.
  (#1 -> 0 &)@@@ Select[Abbr[], !FreeQ[#, DiracChain[__, _k, __]] &] /.
  S|T|U -> 0;

amp = ApplyUnitarity[amp, CKM, 3];

Put[amp, name <> ".amp"];
Put[Abbr[], name <> ".abbr"];


Alfa2 = (Alfa = Sqrt[2]/Pi GF MW2 SW2)^2;

Attributes[ColorDelta] = Attributes[Op] = {Orderless};

dc = Plus@@ amp //. Abbr[] /.
  { Col1 -> 1, Col3 -> 3,
    Col2 -> 2, Col4 -> 4,
    Spinor[_[i:1 | 4], __] -> s[i],
    Spinor[_[i:2 | 3], __] -> b[i] } /.
  SUNT[i_, j_] -> ColorDelta[i, j] /.
  Mat[ DiracChain[s[i_], g1__, b[j_]] ColorDelta[i_, j_] *
       DiracChain[s[k_], g2__, b[l_]] ColorDelta[k_, l_] ] :>
    -Signature[{i, j, k, l}] Op[{g1}, {g2}] /.
  {Sqrt[2] -> sqrt2, 1/Sqrt[2] -> 1/sqrt2};

Put[dc, name <> ".dc"];


OpName[{6, mu_Lor}, {6, mu_}] = "cVLL";
OpName[{7, mu_Lor}, {7, mu_}] = "cVRR";

OpName[{6, mu_Lor}, {7, mu_}] = "cLR1";
OpName[{6}, {7}] = "cLR2";

OpName[{6}, {6}] = "cSLL1";
OpName[{7}, {7}] = "cSRR1";
OpName[{-6, mu_Lor, nu_Lor}, {-6, mu_, nu_}] = "cSLL2";
OpName[{-7, mu_Lor, nu_Lor}, {-7, mu_, nu_}] = "cSRR2";

ops = Union[Cases[dc, _Op, Infinity]];

coeff = ((OpName@@ #) <> tag -> Coefficient[dc, #])&/@ ops;

chk = dc /. Thread[ops -> 0];
If[ chk =!= 0,
  Print["WARNING: REMAINDER NOT ZERO"];
  Put[chk, name <> ".ZERO"] ];


fin[var_ -> expr_] := (Print[var]; var -> fin[expr]);

fin[expr_] := Simplify @ Collect[expr, {_SumOver, _Log}, FullSimplify];

coeff = fin/@ coeff;

Put[coeff, name <> ".coeff"];


hh = OpenFortran[name <> ".F"];

SetOptions[PrepareExpr,
  Optimize -> True,
  Expensive -> {_C0Z, _C00Z, _D0Z, _D00Z}];

SetOptions[WriteExpr,
  TmpType -> False,
  RealArgs -> {C0Z, C00Z, D0Z, D00Z}];

WriteExpr[hh, coeff];

Close[hh];
_EOF_

gzip -f $name*.ps


