(*********************** Program VacPol.m **************************)
(*
 Version compatible with FeynCalc 10.1.0. There are big changes with
 respect to version 9.2 and older versions.

 Date: 09/12/2024
 Author: Jorge C. Romao
 email: jorge.romao@tecnico.ulisboa.pt
*)

(* First input FeynCalc *)
(* Uncomment below if you want to call from this program. If open a new
   mathematica notebook and load FeynCalc from there you should not load
   it again
*)
(*
<< FeynCalc`
*)

(* Now write the numerator of the Feynman diagram. We define the 
   constant

       C=alpha/(4 pi)
 
   I also use the FCE notation available since FeynCalc 6. See manual
   for explanations.
 
*)

(* Everything has to be given in Dimension D *)

num:= - C TR[GAD[mu] . (GSD[q] + m) . GAD[nu] .  (GSD[q]+GSD[k]+m)]
  
(* Set some Options *)

SetOptions[PaVeReduce,A0ToB0->True];
   
(* Define the amplitude *)

amp:=num * FAD[{q+k,m},{q,m}]

(* Calculate the result *)

(* The function OneLoop is deprecated. Now one should use TID, see
   the manual for explanations *)

res:=(-I / Pi^2) TID[amp,q,UsePaVeBasis -> True, ToPaVe -> True] 
Ans=PaVeReduce[res] //Simplify ;
Print["The result in dimension D:\n"]
Print["Ans=",Ans,"\n"]

(* Now we have to take D->4 but one should not forget the finite part
 that comes from eps Div[B0] where D=4-eps and Div[B0]=2/eps *)

(* First the D->4 part *)  

AnsD4 = Ans /. D -> 4;

Print["The D-> 4 part:\n"]
Print["AnsD4=",AnsD4,"\n"]

(* Now the Finite Part *)

AnsFinitePart =  Normal[Series[Ans /. {D -> 4 - eps, PaVe[0,{x_},{y_, y_}] -> 2/eps}, {eps, 0, 0}]] /. eps -> Infinity ;

Print["The Finite part:\n"]
Print["AnsFinitePart=",AnsFinitePart,"\n"]

(* Finally the full Result *)

Print["Final Result:\n"]
FinalAns=AnsD4 + AnsFinitePart //Simplify;

Print["VacPol=",FinalAns]  
  
(******************** End of Program VacPol.m **********************)
