* Cget.F
* the three-point tensor coefficients
* this file is part of LoopTools
* last modified 18 Oct 01 th

#include "defs.h"


	integer function Cget(p1, p2, p1p2, m1, m2, m3)
	implicit none
	double precision p1, p2, p1p2, m1, m2, m3

	integer cachelookup
	external cachelookup, Cijk

	double precision para(6)

	double precision Ccache(2)
	common /Cbase/ Ccache

* The following data statement initializes two *pointers* to NULL
* (see cache.c). This is sneaky but ok because IEEE says 0D0 is
* represented as 8 0-bytes.
	data Ccache /0, 0/

	para(1) = p1
	para(2) = p2
	para(3) = p1p2
	para(4) = m1
	para(5) = m2
	para(6) = m3

	Cget = cachelookup(para, Ccache, Cijk, 6, 13)
	end

************************************************************************

	subroutine Cijk(P, C)
	implicit none
	double precision mtest,mtest2,gtest,nogw
	double complex mwc
	common /resonant/ mwc

	double precision P(6)
	double complex C(13)

	double complex a0i(2), b0p, b1p
	common /bsave/ a0i, b0p, b1p

	double complex B11, C0, B11res, B1res, B0res,B00,B00res
	double complex B1resWA,B11resWA
	external B11, C0, B11res, B1res, B0res,B00,B00res
	external B1resWA,B11resWA

	double complex C0resonant,C0onshIR,prova
	external C0resonant,C0onshIR


	XREAL M11, M12, M22, det2, f1, f2
	double complex b1123, b023, b123
	double complex b1113, b113, b1112, b112
	double complex s1, s2

#ifndef LANGUAGE_FORTRAN_90
	double complex lc1, lc2
	lc1() = dcmplx(M11*dble(s1) + M12*dble(s2),
     +    M11*dimag(s1) + M12*dimag(s2))
	lc2() = dcmplx(M12*dble(s1) + M22*dble(s2),
     +    M12*dimag(s1) + M22*dimag(s2))
#endif

	mtest = dreal(mwc)
	gtest = -dimag(mwc)*2d0
	mtest2 = mtest*mtest

c	mtest = 80.35d0
c	mtest2=80.35d0*80.35d0
c	gtest = 2.08699d0

	M12 = .5D0*(QEXT(P(2)) - QEXT(P(3)) - P(1))
	det2 = 2*(QEXT(P(1))*P(3) - M12*M12)
	M12 = M12/det2
	M11 = P(3)/det2
	M22 = P(1)/det2
	f1 = QEXT(P(4)) - QEXT(P(5)) + P(1)
	f2 = QEXT(P(4)) - QEXT(P(6)) + P(3)

c	print*,P(1),P(2),P(3),P(4),P(5),P(6),mw2

	   nogw = 1d0

	if (P(1).eq.0d0.and.P(3).eq.P(4).and.P(5).eq.mtest2.and.
     -      P(6).eq.0d0.and.P(2).eq.mtest2) then
           C(cc0) = C0onshIR(P(3),mtest,gtest)
	   print*,'in effetti ci passa IR 3',C(cc0)

	else if (P(3).eq.0d0.and.P(2).eq.P(6).and.P(4).eq.mtest2.and.
     -      P(5).eq.0d0.and.P(1).eq.mtest2)  then
           C(cc0) = C0onshIR(P(2),mtest,gtest)
	   print*,'in effetti ci passa IR 4',C(cc0)

ccccccccc  caso 1
	else if (P(2).eq.P(6).and.P(4).eq.mtest2.and.P(5).eq.0d0) then

! 100 GeV^2 is for sure above any threshold!
	   if (P(1).lt.100d0) nogw = 0d0
           C(cc0) = C0resonant(P(2), P(3), P(1), P(5), P(6),mtest,nogw*gtest)
	   nogw = 1d0

*****	   prova = C0(P(1), P(2), P(3), P(4), P(5), P(6))
c	   print*,'in effetti ci passa 1',C(cc0)
c	   print*,'prova 1              ',prova
	b1123 = B11(P(2), P(6), P(5))
	b023 = b0p
	b123 = b1p
	b1113 = B11(P(3), P(4), P(6))
	b113 = b1p
	s1 = b0p - b023 - f1*C(cc0)
	b1112 = B11resWA(P(1), mtest,nogw*gtest)
	b112 = B1resWA(P(1), mtest,nogw*gtest)
        b0p = B0res(P(1), mtest, nogw*gtest)

c	print*,'b1123=',b1123
c	print*,'b023 =',b023
c	print*,'b123 =',b123
c	print*,'b1113=',b1113
c	print*,'b113 =',b113
c	print*,'b1123=',b1112
c	print*,'b123 =',b112
c	print*,'s1   =',s1


ccccccccc  caso 2
	else if (P(2).eq.P(5).and.P(4).eq.mtest2.and.P(6).eq.0d0) then

	   if (P(3).lt.100d0) nogw = 0d0
           C(cc0) = C0resonant(P(1), P(2), P(3), P(6), P(5),mtest,nogw*gtest)
	   nogw = 1d0

******	   prova = C0(P(1), P(2), P(3), P(4), P(5), P(6))
c	   print*,'in effetti ci passa 2',C(cc0)
c	   print*,'prova 2              ',prova
	b1123 = B11(P(2), P(6), P(5))
	b023 = b0p
	b123 = b1p
	b1113 = B11resWA(P(3), mtest, nogw*gtest)
	b113 = B1resWA(P(3), mtest, nogw*gtest)
        b0p = B0res(P(3), mtest, nogw*gtest)
	s1 = b0p - b023 - f1*C(cc0)
	b1112 = B11(P(1), P(4), P(5))
	b112 = b1p


ccccccccc  caso 3
	else if (P(3).eq.P(6).and.P(5).eq.mtest2.and.P(4).eq.0d0) then

	   if (P(1).lt.100d0) nogw = 0d0
           C(cc0) = C0resonant(P(2), P(3), P(1), P(4), P(6),mtest,nogw*gtest)
	   nogw = 1d0
*****	   prova = C0(P(1), P(2), P(3), P(4), P(5), P(6))
c	   print*,'in effetti ci passa 3',C(cc0)
c	   print*,'prova 3              ',prova
	b1123 = B11(P(2), P(6), P(5))
	b023 = b0p
	b123 = b1p
	b1113 = B11(P(3), P(4), P(6))
	b113 = b1p
	s1 = b0p - b023 - f1*C(cc0)
	b1112 = B11res(P(1), mtest, nogw*gtest)
	b112 = B1res(P(1), mtest, nogw*gtest)
        b0p = B0res(P(1), mtest, nogw*gtest)

c	print*,'b1123=',b1123
c	print*,'b023 =',b023
c	print*,'b123 =',b123
c	print*,'b1113=',b1113
c	print*,'b113 =',b113
c	print*,'b1123=',b1112
c	print*,'b123 =',b112
c	print*,'s1   =',s1


ccccccccc  caso 4
	else if (P(3).eq.P(4).and.P(5).eq.mtest2.and.P(6).eq.0d0) then

	   if (P(2).lt.100d0) nogw = 0d0
           C(cc0) = C0resonant(P(1), P(3), P(2), P(6), P(4),mtest,nogw*gtest)
	   nogw = 1d0

****	   prova = C0(P(1), P(2), P(3), P(4), P(5), P(6))
c	   print*,'in effetti ci passa 4',C(cc0)
c	   print*,'prova 4              ',prova
	b1123 = B11res(P(2), mtest, nogw*gtest)
	b023 = B0res(P(2), mtest, nogw*gtest)
	b123 = B1res(P(2), mtest, nogw*gtest)
	b1113 = B11(P(3), P(4), P(6))
	b113 = b1p
	s1 = b0p - b023 - f1*C(cc0)
	b1112 = B11(P(1), P(4), P(5))
	b112 = b1p

c	print*,'b1123=',b1123
c	print*,'b023 =',b023
c	print*,'b123 =',b123
c	print*,'b1113=',b1113
c	print*,'b113 =',b113
c	print*,'b1123=',b1112
c	print*,'b123 =',b112
c	print*,'s1   =',s1


ccccccccc  caso 5
	else if (P(1).eq.P(5).and.P(6).eq.mtest2.and.P(4).eq.0d0) then

	   if (P(3).lt.100d0) nogw = 0d0
           C(cc0) = C0resonant(P(1), P(2), P(3), P(4), P(5),mtest,nogw*gtest)
	   nogw = 1d0

*****	   prova = C0(P(1), P(2), P(3), P(4), P(5), P(6))
c	   print*,'in effetti ci passa 5',C(cc0)
c	   print*,'prova 5              ',prova
	b1123 = B11(P(2), P(6), P(5))
	b023 = b0p
	b123 = b1p
	b1113 = B11res(P(3), mtest,nogw*gtest)
	b113 = B1res(P(3), mtest,nogw*gtest)
	b0p = B0res(P(3), mtest,nogw*gtest)
	s1 = b0p - b023 - f1*C(cc0)
	b1112 = B11(P(1), P(4), P(5))
	b112 = b1p


ccccccc caso 6
         else if (P(1).eq.P(4).and.P(6).eq.mtest2.and.P(5).eq.0d0)  then

	   if (P(2).lt.100d0) nogw = 0d0
           C(cc0) = C0resonant(P(1), P(3), P(2), P(5), P(4),mtest,nogw*gtest)
	   nogw = 1d0

******	   prova = C0(P(1), P(2), P(3), P(4), P(5), P(6))
c	   print*,'in effetti ci passa 6',C(cc0)
c	   print*,'prova 6              ',prova
	b1123 = B11resWA(P(2), mtest, nogw*gtest)
	b023 = B0res(P(2), mtest, nogw*gtest)
	b123 = B1resWA(P(2), mtest, nogw*gtest)
	b1113 = B11(P(3), P(4), P(6))
	b113 = b1p
	s1 = b0p - b023 - f1*C(cc0)
	b1112 = B11(P(1), P(4), P(5))
	b112 = b1p

c	print*,'S=',P(2)
c	print*,'b1123=',b1123
c	print*,'b023 =',b023
c	print*,'b123 =',b123
c	print*,'b1113=',b1113
c	print*,'b113 =',b113
c	print*,'b1123=',b1112
c	print*,'b123 =',b112
c	print*,'s1   =',s1

	   else

	C(cc0) = C0(P(1), P(2), P(3), P(4), P(5), P(6))
c	print*,'valore normale=',C(cc0)
	b1123 = B11(P(2), P(6), P(5))
	b023 = b0p
	b123 = b1p
	b1113 = B11(P(3), P(4), P(6))
	b113 = b1p
	s1 = b0p - b023 - f1*C(cc0)
	b1112 = B11(P(1), P(4), P(5))
	b112 = b1p

	end if

c	b1123 = B11(P(2), P(6), P(5))
c	b023 = b0p
c	b123 = b1p
c	b1113 = B11(P(3), P(4), P(6))
c	b113 = b1p
c	s1 = b0p - b023 - f1*C(cc0)
c	b1112 = B11(P(1), P(4), P(5))
c	b112 = b1p

c	print*,'senza larghezza',P(2),P(6),P(5)
c	print*,'b1123=',b1123
c	print*,'b023 =',b023
c	print*,'b123 =',b123
c	print*,'b1113=',b1113
c	print*,'b113 =',b113
c	print*,'b1123=',b1112
c	print*,'b123 =',b112
c	print*,'s1   =',s1


	s2 = b0p - b023 - f2*C(cc0)
	C(cc1) = lc1()
	C(cc2) = lc2()

	C(cc00) = .5D0*(P(4)*C(cc0) +
     +    .5D0*(f1*C(cc1) + f2*C(cc2) + b023)) + .25D0
	s1 = -(f1*C(cc1) + b123) - 2*C(cc00)
	s2 = -(f2*C(cc1) + b123 - b112)
	C(cc11) = lc1()
	C(cc12) = lc2()
	b023 = b023 + b123
	s1 = b023 + b113 - f1*C(cc2)
	s2 = b023 - f2*C(cc2) - 2*C(cc00)
	C(cc22) = lc2()


!	print*,'Css ',C(cc1),C(cc2),C(cc00)

	C(cc001) = (P(4)*C(cc1) + .5D0*(f1*C(cc11) +
     +    f2*C(cc12) + b123))/3D0 - 1/18D0
	C(cc002) = (P(4)*C(cc2) + .5D0*(f1*C(cc12) +
     +    f2*C(cc22) - b023))/3D0 - 1/18D0
	s1 = -(b1123 + f1*C(cc11)) - 4*C(cc001)
	s2 = -(b1123 + f2*C(cc11) - b1112)
	C(cc111) = lc1()
	C(cc112) = lc2()
	b1123 = b1123 + b023 + b123
	s1 = -(b1123 + f1*C(cc22) - b1113)
	s2 = -(b1123 + f2*C(cc22)) - 4*C(cc002)
	C(cc122) = lc1()
	C(cc222) = lc2()

#ifdef LANGUAGE_FORTRAN_90
	contains

	double complex function lc1
	lc1 = dcmplx(M11*dble(s1) + M12*dble(s2),
     +    M11*dimag(s1) + M12*dimag(s2))
	end function lc1

	double complex function lc2
	lc2 = dcmplx(M12*dble(s1) + M22*dble(s2),
     +    M12*dimag(s1) + M22*dimag(s2))
	end function lc2
#endif
	end

	double complex function C0resonant(p12,p1p2,p22,m02,m12,MW,GW)
	double precision p12,p1p2,p22,m02,m12,pi2
	double precision MW,GW
	double complex cli2,dummy

	integer ier,ipi12
	parameter (pi2=3.1415926**2)

	call ffzxdl(cli2,ipi12,dummy,(mw**2-p22)/mw**2,1,ier)


	C0resonant = 1d0/p22*
     -                    (log(p22/m12)*log((MW**2-dcmplx(0d0,GW*MW)-
     -                                       p22)/MW**2)+
     -                     dble(cli2)+ipi12*pi2/12d0-pi2/6d0)

	end function C0resonant

	double complex function C0onshIR(m,MW,GW)
	double precision m,MW,GW,lambda2
	external getlambda
	lambda2 = getlambda()
	print*,'cutoff=',lambda2,m
	C0onshIR = 1d0/MW**2*log(MW/sqrt(m))*log(lambda2/sqrt(m)/MW)


	end function C0onshIR