*---------------------------------------------------------------*
* This program calculates  Z^* Higgs Gamma 3-point function     *
* The notation  and other explanations can be found at the      *
* page:                                                         *
*                                                               *
*       http://porthos.ist.utl.pt/CTQFT                         *  
*      (See also Nuclear Physics B267(1986)509)                 *
*                                                               *
*    Needs: LoopTools, CUBA integration package, gauss.f for    * 
*    gaussian integration.                                      *
*                                                               *
*    Author: Jorge Crispim Romao                                *
*    email:jorge.romao@ist.utl.pt                               *
*    Last Change:20/10/2008                                     *
*---------------------------------------------------------------*

      program ghgtest
      implicit none


**************  Standard Model Physical Constants ******
*
      INCLUDE 'smpc.h'
*
********************************************************

*
*** Variables in the program
*
      integer i,j,iboson
      real*8 lbd2    ! For LoopTools
      real*8 g,G_F,mW,mZ,thetaW,mH,GW,GZ,pi
      real*8 ee,alpha,alphaz
      real*8 mlep(3),mup(3),mdown(3)
      real*8 q2,beta,betaH,mH2,mW2,mf,mf2,betap,betapH,tw2
      complex*16 cf,d1,d3,d1paper,d3paper,cfpaper,S,J1,J2,L
      complex*16 c8c13,c8c13paper,d4,d4paper

      common/data/g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi,iboson
      common/data2/mlep,mup,mdown

*
*** Initialize LoopTools
*
#include "looptools.h"

      call ffini

      lbd2=1.d0
      call setlambda(lbd2)


      pi=acos(-1.d0)
      mW=SM_MW
      GW=SM_GAW
      mZ=SM_MZ
      GZ=SM_GAZ
      G_F=SM_GFERMI
      mlep(1)=SM_EMASS
      mlep(2)=SM_MUMASS
      mlep(3)=SM_TAUMASS
      mup(3)=SM_TOPMASS
      mup(2)=SM_CHARMMASS
      mup(1)=SM_UPMASS
      mdown(1)=SM_DOWNMASS
      mdown(2)=SM_STRANGEMASS
      mdown(3)=SM_BOTTOMMASS
      alpha=SM_ALPHA
      alphaz=SM_ALPHAZ
      thetaW=asin(sqrt(SM_S2TW))
      ee=sqrt(4*pi*alpha)
      g=ee/sin(thetaW)

      mW=mZ*cos(thetaW)
      mW2=mW*mW
      call setmudim(mW**2)

      mH=50d0
      q2=1d5
      mf=170d0
      mf2=mf*mf
      tw2=tan(thetaW)**2

      do i=1,5

         mH=mH+30d0
         mH2=mH*mH
         betap=q2/mf2
         betapH=mH2/mf2
         beta=q2/mW2
         betaH=mH2/mW2

        cf=-(4*mf2*(2*q2*B0(mH2,mf2,mf2) 
     &        - 2*q2*B0(q2,mf2,mf2) - 
     &        (mH2 - q2)*(-2 + (-4*mf2 + mH2 - q2)*
     &        C0(0,q2,mH2,mf2,mf2,mf2))))/(mH2 - q2)**2
        cfpaper=-4*(J1(betap,betapH)-4d0*J2(betap,betapH))

       write(*,98)mH,cf,cfpaper

        c8c13=-((mH2*(-1 + tw2)*(q2*B0(mH2,mW2,mW2) 
     &       - q2*B0(q2,mW2,mW2) + 
     &       (mH2 - q2)*(1 + 2*mW2*C0(0,q2,mH2,mW2,mW2,mW2))))/
     &       (mH2 - q2)**2)

        c8c13paper=-2d0*(1d0-tw2)*mH2/mW2*J2(beta,betaH)

       write(*,98)mH,c8c13,c8c13paper


       d1=-((mH2*(-mH2 + q2) + 10*mW2*(-mH2 + q2) 
     &      - mH2*(-mH2 + q2)*tw2 - 
     &      2*mW2*(-mH2 + q2)*tw2 - 10*mW2*q2*B0(mH2,mW2,mW2) + 
     &      mH2*q2*(-1 + tw2)*B0(mH2,mW2,mW2) + 
     &      2*mW2*q2*tw2*B0(mH2,mW2,mW2) + mH**4*B0(q2,mW2,mW2) + 
     &      mH2*(-mH2 + q2)*B0(q2,mW2,mW2) - 
     &      2*mW2*q2*(-5 + tw2)*B0(q2,mW2,mW2) - 
     &      mH**4*tw2*B0(q2,mW2,mW2) - 
     &      mH2*(-mH2 + q2)*tw2*B0(q2,mW2,mW2) + 
     &      2*mH2*mW2*(-mH2 + q2)*C0(0,q2,mH2,mW2,mW2,mW2) + 
     &      20*mW**4*(-mH2 + q2)*C0(0,q2,mH2,mW2,mW2,mW2) + 
     &      12*mW2*(-mH2 + q2)**2*C0(0,q2,mH2,mW2,mW2,mW2) - 
     &      2*mH2*mW2*(-mH2 + q2)*tw2*C0(0,q2,mH2,mW2,mW2,mW2) - 
     &      4*mW**4*(-mH2 + q2)*tw2*C0(0,q2,mH2,mW2,mW2,mW2) - 
     &      4*mW2*(-mH2 + q2)**2*tw2*C0(0,q2,mH2,mW2,mW2,mW2))/
     &      (-mH2 + q2)**2)


       d1paper=4d0*(3d0-tw2)*J1(beta,betaH)
     &      +4d0*(-5d0+tw2-0.5d0*betaH*(1d0-tw2))*J2(beta,betaH)
         
       write(*,98)mH,d1,d1paper

       d3=(mW2*(-((1 + tw2)*B0(mH2,mW2,mW2)) + 
     &      (1 + tw2)*B0(q2,mW2,mW2) + 
     &      2*(-mH2 + q2)*(-2 + tw2)*C0(0,q2,mH2,mW2,mW2,mW2)))/
     &      (2.*(-mH2 + q2))


       d3paper=-1d0/(beta-betaH)*(0.5d0*(1d0+tw2)*(L(beta)-L(betaH))
     &      +(2d0-tw2)*(S(beta)-S(betaH)))


       write(*,98)mH,d3,d3paper

       d4=-3*mW2*C0(0,q2,mH2,mW2,mW2,mW2)

       d4paper=-3d0/(beta-betaH)*(S(beta)-S(betaH))

       write(*,98)mH,d4,d4paper

      enddo


 98   format(f8.2,4x,2("(",1p,e12.5,",",1p,e12.5,")",2x))


      end

*-------------------------------------------------------------------*
      complex*16 function J1(x,y)
      implicit none
      real*8 x,y
      complex*16 S

      J1=-(S(x)-S(y))/(x-y)
    
      return
      end

*-------------------------------------------------------------------*
      complex*16 function J2(x,y)
      implicit none
      real*8 x,y
      complex*16 S,L

      J2=0.5d0/(x-y)+1d0/(x-y)**2*(S(x)-S(y)-0.5d0*x*(L(x)-L(y)))
    
      return
      end

*-------------------------------------------------------------------*
      complex*16 function S(x)
      implicit none
      real*8 x,pi,ReS,ImS

      pi=acos(-1d0)

      if(x.lt.4d0)then
         S=-2d0*(asin(sqrt(x)/2d0))**2
      else
         reS=-pi**2/2d0+0.5d0*(log((x-2d0+sqrt(x**2-4d0*x))/2d0))**2
         ImS=-2d0*pi*log(sqrt(x)/2d0+sqrt(x/4d0-1d0))
         S=cmplx(reS,ImS)
      endif
      return
      end



*-------------------------------------------------------------------*
      complex*16 function L(x)
      implicit none
      real*8 x,ReL,ImL,pi
      
      pi=acos(-1d0)

      if(x.lt.4d0)then
         L=-2d0+2d0*sqrt((4d0-x)/x)*asin(sqrt(x)/2d0)
      else
         ReL=-2d0+sqrt((x-4d0)/x)*log((x-2d0+sqrt(x*(x-4d0)))/2d0)
         ImL=-pi*sqrt((x-4d0)/x)
         L=cmplx(ReL,ImL)
      endif
      return
      end



*-------------------------------------------------------------------*
      complex*16 function F12(x)
      implicit none
      real*8 x
      complex*16 ftau

      F12=-2d0*x*(1d0+(1d0-x)*ftau(x))

      return
      end

*-------------------------------------------------------------------*
      complex*16 function ftau(x)
      implicit none
      real*8 x,realftau,imgftau,etaplus,etaminus,pi

      pi=acos(-1d0)

      if(x.gt.1d0)then
         ftau=asin(sqrt(1/x))**2
      else
         etaplus=1d0+sqrt(1d0-x)
         etaminus=1d0-sqrt(1d0-x)
         realftau=pi**2/4d0-1d0/4d0*log(etaplus/etaminus)**2
         imgftau=pi/2d0*log(etaplus/etaminus)
         ftau=cmplx(realftau,imgftau)
      endif
      return
      end
*-------------------------------------------------------------------*

*-------------------------------------------------------------------*
      complex*16 function gtau(x)
      implicit none
      real*8 x,realgtau,imggtau,etaplus,etaminus,pi

      pi=acos(-1d0)

      if(x.gt.1d0)then
         gtau=sqrt(x-1d0)*asin(sqrt(1/x))
      else
         etaplus=1d0+sqrt(1d0-x)
         etaminus=1d0-sqrt(1d0-x)
         realgtau=0.5d0*sqrt(1d0-x)*log(etaplus/etaminus)
         imggtau=-pi/2d0*sqrt(1d0-x)
         gtau=cmplx(realgtau,imggtau)
      endif
      return
      end
*-------------------------------------------------------------------*

*-------------------------------------------------------------------*
      complex*16 function I1(x,y)
      implicit none
      real*8 x,y
      complex*16 ftau,gtau

      I1=x*y/2d0/(x-y)+(x*y)**2/2d0/(x-y)**2*(ftau(x)-ftau(y))+
     &     x**2*y/(x-y)**2*(gtau(x)-gtau(y))
    
      return
      end

*-------------------------------------------------------------------*
      complex*16 function I2(x,y)
      implicit none
      real*8 x,y
      complex*16 ftau

      I2=-x*y/2d0/(x-y)*(ftau(x)-ftau(y))
    
      return
      end

*-------------------------------------------------------------------*
