      PROGRAM decay
*
      IMPLICIT NONE
*
      REAL*8 scalar
      REAL*8, dimension(3) :: vector
*
      INTERFACE
       FUNCTION SHvec(amh,amt,amb,als,sct,scb,scw) RESULT(value)
       IMPLICIT NONE
       REAL*8 amh,amt,amb,sct,scb,scw,als
       REAL*8, dimension(3) :: value
       END FUNCTION SHvec
      END INTERFACE
*       
      INTERFACE
       FUNCTION SHgamgam(amh,amt,amb,sct,scb,scw)  
       IMPLICIT NONE
       REAL*8 amh,amt,amb,sct,scb,scw,SHgamgam
       END FUNCTION SHgamgam
      END INTERFACE
* 
      INTERFACE
       FUNCTION SHgluglu(amh,amt,amb,als,sct,scb)  
       IMPLICIT NONE
       REAL*8 amh,amt,amb,sct,scb,als,SHgluglu
       END FUNCTION SHgluglu
      END INTERFACE
* 
      INTERFACE
       FUNCTION SHZgam(amh,amt,amb,sct,scb,scw)  
       IMPLICIT NONE
       REAL*8 amh,amt,amb,sct,scb,scw,SHZgam
       END FUNCTION SHZgam
      END INTERFACE
* 
       vector= SHvec(125.d0,172.5d0,4.69d0,0.11d0,1.d0,1.d0,1.d0)
       print 1099,vector(1),vector(2),vector(3)
       scalar= SHgamgam(125.d0,172.5d0,4.69d0,1.d0,1.d0,1.d0)
       print 1098,scalar
       scalar= SHgluglu(125.d0,172.5d0,4.69d0,0.11d0,1.d0,1.d0)
       print 1098,scalar
       scalar= SHZgam(125.d0,172.5d0,4.69d0,1.d0,1.d0,1.d0)
       print 1098,scalar
*
 1099  format(1x,3e20.5)
 1098  format(1x,e20.5)
*
       STOP
*
       END PROGRAM 
*
*--------------------------------------------------------------------------
*
      FUNCTION SHvec(amh,amt,amb,als,sct,scb,scw) RESULT(value) 
*
*     amh       Higgs mass [GeV]      
*     amt       top   mass [GeV]      
*     amb       b     mass [GeV]      
*     als       alpha_s
*
* OUTPUT value(1) gAmma(H --> gamma gamma)
* OUTPUT value(2) gAmma(H --> glue glue)
* OUTPUT value(3) gAmma(H --> Z gamma)
*
      IMPLICIT NONE
*
      REAL*8, parameter :: alpha=137.0359895d0
      REAL*8, parameter :: pi= 3.141592653589793238462643d0
      REAL*8, parameter :: mw= 80.398d0
      REAL*8, parameter :: mz= 91.17667d0
      REAL*8, parameter :: GF= 1.16637d-5
*
      REAL*8 amh,amt,amhs,mhc,amts,amb,ambs,mws,tt,tb,btt,btb,lt,lb,
     #       ftr,fti,fbr,fbi,arg,atr,ati,abr,abi,tw,arw,aiw,sct,scb,scw,
     #       btw,lw,fwr,fwi,ar,ai,arq,aiq,als,rt,rb,rw,gtr,gti,gbr,gbi,
     #       gwr,gwi,fltr,flti,flbr,flbi,flwr,flwi,gltr,glti,glbr,glbi,
     #       glwr,glwi,mzs,cts,sts,tts,cbtr,cbti,cbbr,cbbi,cbwr,cbwi,
     #       i1tr,i1ti,i1br,i1bi,i1wr,i1wi,i2tr,i2ti,i2br,i2bi,i2wr,
     #       cth,sth,tth,i2wi
      REAL*8, dimension(3) :: value
*
      amhs= amh*amh
      mhc= amhs*amh
      amts= amt*amt
      ambs= amb*amb
*
      tt= 0.25d0*amhs/amts
*
      IF(tt > 1.d0) THEN
       btt= sqrt(1.d0-1.d0/tt)
       lt= log((1.d0+btt)/(1.d0-btt))
       ftr= -0.25d0*(lt*lt-pi*pi)
       fti= 0.5d0*pi*lt
      ELSE
       arg= sqrt(tt)
       ftr= asin(arg)*asin(arg)
       fti= 0.d0
      ENDIF
*
      atr= 2.d0*(tt+(tt-1)*ftr)/(tt*tt)
      ati= 2.d0*(tt-1)*fti/(tt*tt)
*
      tb= 0.25d0*amhs/ambs
*
      IF(tb > 1.d0) THEN
       btb= sqrt(1.d0-1.d0/tb)
       lb= log((1.d0+btb)/(1.d0-btb))
       fbr= -0.25d0*(lb*lb-pi*pi)
       fbi= 0.5d0*pi*lb
      ELSE
       arg= sqrt(tb)
       fbr= asin(arg)*asin(arg)
       fbi= 0.d0
      ENDIF
*
      abr= 2.d0*(tb+(tb-1)*fbr)/(tb*tb)
      abi= 2.d0*(tb-1)*fbi/(tb*tb)
*
      mws= mw*mw
      tw= 0.25d0*amhs/mws
*
      IF(tw > 1.d0) THEN
       btw= sqrt(1.d0-1.d0/tw)
       lw= log((1.d0+btw)/(1.d0-btw))
       fwr= -0.25d0*(lw*lw-pi*pi)
       fwi= 0.5d0*pi*lw
      ELSE
       arg= sqrt(tw)
       fwr= asin(arg)*asin(arg)
       fwi= 0.d0
      ENDIF
*
      arw= -((2.d0*tw+3.d0)*tw+3.d0*(2.d0*tw-1.d0)*fwr)/(tw*tw)
      aiw= -3.d0*(2.d0*tw-1.d0)*fwi/(tw*tw)
*
      arq= sct*atr+scb*abr
      aiq= sct*ati+scb*abi
*
      ar= 4.d0/3.d0*sct*atr+1.d0/3.d0*scb*abr+scw*arw
      ai= 4.d0/3.d0*sct*ati+1.d0/3.d0*scb*abi+scw*aiw
*
      value(1)= GF/(alpha*alpha)/128.d0/sqrt(2.d0)/pi**3*
     #     mhc*(ar*ar+ai*ai)
*
      value(2)= GF*als*als/64.d0/sqrt(2.d0)/pi**3*
     #     mhc*(arq*arq+aiq*aiq)*(1.d0+(95.d0/4.d0-35.d0/6.d0)*als)
*
*---
*
      IF(tt < 1.d0) THEN
       btt= sqrt(1.d0/tt-1.d0)
       arg= sqrt(tt)
       gtr= btt*asin(arg)
       gti= 0.d0
      ELSE
       btt= sqrt(1.d0-1.d0/tt)
       lt= log((1.d0+btt)/(1.d0-btt))
       gtr= 0.5d0*btt*lt
       gti= -0.5d0*pi*btt
      ENDIF
*
      IF(tb < 1.d0) THEN
       btb= sqrt(1.d0/tb-1.d0)
       arg= sqrt(tb)
       gbr= btb*asin(arg)
       gbi= 0.d0
      ELSE
       btb= sqrt(1.d0-1.d0/tb)
       lb= log((1.d0+btb)/(1.d0-btb))
       gbr= 0.5d0*btb*lb
       gbi= -0.5d0*pi*btb
      ENDIF
*
      IF(tw < 1.d0) THEN
       btw= sqrt(1.d0/tw-1.d0)
       arg= sqrt(tw)
       gwr= btw*asin(arg)
       gwi= 0.d0
      ELSE
       btw= sqrt(1.d0-1.d0/tw)
       lw= log((1.d0+btw)/(1.d0-btw))
       gwr= 0.5d0*btw*lw
       gwi= -0.5d0*pi*btw
      ENDIF
*
      mzs= mz*mz
*
      rt= 0.25d0*mzs/amts
      rb= 0.25d0*mzs/ambs
      rw= 0.25d0*mzs/mws
*
      IF(rt > 1.d0) THEN
       btt= sqrt(1.d0-1.d0/rt)
       lt= log((1.d0+btt)/(1.d0-btt))
       fltr= -0.25d0*(lt*lt-pi*pi)
       flti= 0.5d0*pi*lt
      ELSE
       arg= sqrt(rt)
       fltr= asin(arg)*asin(arg)
       flti= 0.d0
      ENDIF
*      
      IF(rb > 1.d0) THEN
       btb= sqrt(1.d0-1.d0/rb)
       lb= log((1.d0+btb)/(1.d0-btb))
       flbr= -0.25d0*(lb*lb-pi*pi)
       flbi= 0.5d0*pi*lb
      ELSE
       arg= sqrt(rb)
       flbr= asin(arg)*asin(arg)
       flbi= 0.d0
      ENDIF
*
      IF(rw > 1.d0) THEN
       btw= sqrt(1.d0-1.d0/rw)
       lw= log((1.d0+btw)/(1.d0-btw))
       flwr= -0.25d0*(lw*lw-pi*pi)
       flwi= 0.5d0*pi*lw
      ELSE
       arg= sqrt(rw)
       flwr= asin(arg)*asin(arg)
       flwi= 0.d0
      ENDIF
*
      IF(rt < 1.d0) THEN
       btt= sqrt(1.d0/rt-1.d0)
       arg= sqrt(rt)
       gltr= btt*asin(arg)
       glti= 0.d0
      ELSE
       btt= sqrt(1.d0-1.d0/rt)
       lt= log((1.d0+btt)/(1.d0-btt))
       gltr= 0.5d0*btt*lt
       glti= -0.5d0*pi*btt
      ENDIF
*
      IF(rb < 1.d0) THEN
       btb= sqrt(1.d0/rb-1.d0)
       arg= sqrt(rb)
       glbr= btb*asin(arg)
       glbi= 0.d0
      ELSE
       btb= sqrt(1.d0-1.d0/rb)
       lb= log((1.d0+btb)/(1.d0-btb))
       glbr= 0.5d0*btb*lb
       glbi= -0.5d0*pi*btb
      ENDIF
*
      IF(rw < 1.d0) THEN
       btw= sqrt(1.d0/rw-1.d0)
       arg= sqrt(rw)
       glwr= btw*asin(arg)
       glwi= 0.d0
      ELSE
       btw= sqrt(1.d0-1.d0/rw)
       lw= log((1.d0+btw)/(1.d0-btw))
       glwr= 0.5d0*btw*lw
       glwi= -0.5d0*pi*btw
      ENDIF
*
      i1tr= -0.5d0/(tt-rt)+0.5d0/(tt-rt)**2*(ftr-fltr)+
     #      rt/(tt-rt)**2*(gtr-gltr)
      i1ti= 0.5d0/(tt-rt)**2*(fti-flti)+rt/(tt-rt)**2*(gti-glti)
*
      i1br= -0.5d0/(tb-rb)+0.5d0/(tb-rb)**2*(fbr-flbr)+
     #      rb/(tb-rb)**2*(gbr-glbr)
      i1bi= 0.5d0/(tb-rb)**2*(fbi-flbi)+rb/(tb-rb)**2*(gbi-glbi)
*
      i1wr= -0.5d0/(tw-rw)+0.5d0/(tw-rw)**2*(fwr-flwr)+
     #      rw/(tw-rw)**2*(gwr-glwr)
      i1wi= 0.5d0/(tw-rw)**2*(fwi-flwi)+rw/(tw-rw)**2*(gwi-glwi)
*
      i2tr= 0.5d0/(tt-rt)*(ftr-fltr)
      i2ti= 0.5d0/(tt-rt)*(fti-flti)
*
      i2br= 0.5d0/(tb-rb)*(fbr-flbr)
      i2bi= 0.5d0/(tb-rb)*(fbi-flbi)
*
      i2wr= 0.5d0/(tw-rw)*(fwr-flwr)
      i2wi= 0.5d0/(tw-rw)*(fwi-flwi)
*
      cbtr= i1tr-i2tr
      cbti= i1ti-i2ti
*
      cbbr= i1br-i2br
      cbbi= i1bi-i2bi
*
      cts= mws/mzs
      sts= 1.d0-cts
      cth= sqrt(cts) 
      sth= sqrt(sts)
      tts= sts/cts
      tth= sqrt(tts)
*  
      cbwr= 1.d0/tth*(4.d0*(3.d0-tts))*i2wr+((1.d0+2.d0*tw)*tts-
     #      (5.d0+2.d0*tw))*i1wr
      cbwi= 1.d0/tth*(4.d0*(3.d0-tts))*i2wi+((1.d0+2.d0*tw)*tts-
     #      (5.d0+2.d0*tw))*i1wi
*
      ar= sct*(2.d0*(1.d0-8.d0/3.d0*sts)/(sth*cth)*cbtr)
      ar= ar+scb*(-(-1.d0+4.d0/3.d0*sts)/(sth*cth)*cbbr)
      ar= ar+scw*cbwr
*
      ai= sct*(2.d0*(1.d0-8.d0/3.d0*sts)/(sth*cth)*cbti)
      ai= ai+scb*(-(-1.d0+4.d0/3.d0*sts)/(sth*cth)*cbbi)
      ai= ai+scw*cbwi
*
      value(3)= GF/(alpha*alpha)/(64.d0*sqrt(2.d0)*pi**3)*mhc*
     #          (1.d0-mzs/amhs)**3*(ar*ar+ai*ai)
*
      RETURN
*
      END FUNCTION SHvec
*
*------------------------------------------------------------------------------------
*
      FUNCTION SHgamgam(amh,amt,amb,sct,scb,scw)  
*
*     amh       Higgs mass [GeV]      
*     amt       top   mass [GeV]      
*     amb       b     mass [GeV]      
*
* OUTPUT gAmma(H --> gamma gamma)
*
      IMPLICIT NONE
*
      REAL*8, parameter :: alpha=137.0359895d0
      REAL*8, parameter :: pi= 3.141592653589793238462643d0
      REAL*8, parameter :: mw= 80.398d0
      REAL*8, parameter :: GF= 1.16637d-5
*
      REAL*8 amh,amt,amhs,mhc,amts,amb,ambs,mws,tt,tb,btt,btb,lt,lb,
     #       ftr,fti,fbr,fbi,arg,atr,ati,abr,abi,tw,arw,aiw,sct,scb,scw,
     #       btw,lw,fwr,fwi,ar,ai,arq,aiq,als,rt,rb,rw,gtr,gti,gbr,gbi,
     #       gwr,gwi,fltr,flti,flbr,flbi,flwr,flwi,gltr,glti,glbr,glbi,
     #       glwr,glwi,mzs,cts,sts,tts,cbtr,cbti,cbbr,cbbi,cbwr,cbwi,
     #       i1tr,i1ti,i1br,i1bi,i1wr,i1wi,i2tr,i2ti,i2br,i2bi,i2wr,
     #       cth,sth,tth,i2wi,SHgamgam
*
      amhs= amh*amh
      mhc= amhs*amh
      amts= amt*amt
      ambs= amb*amb
*
      tt= 0.25d0*amhs/amts
*
      IF(tt > 1.d0) THEN
       btt= sqrt(1.d0-1.d0/tt)
       lt= log((1.d0+btt)/(1.d0-btt))
       ftr= -0.25d0*(lt*lt-pi*pi)
       fti= 0.5d0*pi*lt
      ELSE
       arg= sqrt(tt)
       ftr= asin(arg)*asin(arg)
       fti= 0.d0
      ENDIF
*
      atr= 2.d0*(tt+(tt-1)*ftr)/(tt*tt)
      ati= 2.d0*(tt-1)*fti/(tt*tt)
*
      tb= 0.25d0*amhs/ambs
*
      IF(tb > 1.d0) THEN
       btb= sqrt(1.d0-1.d0/tb)
       lb= log((1.d0+btb)/(1.d0-btb))
       fbr= -0.25d0*(lb*lb-pi*pi)
       fbi= 0.5d0*pi*lb
      ELSE
       arg= sqrt(tb)
       fbr= asin(arg)*asin(arg)
       fbi= 0.d0
      ENDIF
*
      abr= 2.d0*(tb+(tb-1)*fbr)/(tb*tb)
      abi= 2.d0*(tb-1)*fbi/(tb*tb)
*
      mws= mw*mw
      tw= 0.25d0*amhs/mws
*
      IF(tw > 1.d0) THEN
       btw= sqrt(1.d0-1.d0/tw)
       lw= log((1.d0+btw)/(1.d0-btw))
       fwr= -0.25d0*(lw*lw-pi*pi)
       fwi= 0.5d0*pi*lw
      ELSE
       arg= sqrt(tw)
       fwr= asin(arg)*asin(arg)
       fwi= 0.d0
      ENDIF
*
      arw= -((2.d0*tw+3.d0)*tw+3.d0*(2.d0*tw-1.d0)*fwr)/(tw*tw)
      aiw= -3.d0*(2.d0*tw-1.d0)*fwi/(tw*tw)
*
      arq= sct*atr+scb*abr
      aiq= sct*ati+scb*abi
*
      ar= 4.d0/3.d0*sct*atr+1.d0/3.d0*scb*abr+scw*arw
      ai= 4.d0/3.d0*sct*ati+1.d0/3.d0*scb*abi+scw*aiw
*
      SHgamgam= GF/(alpha*alpha)/128.d0/sqrt(2.d0)/pi**3*
     #          mhc*(ar*ar+ai*ai)
*
*
      RETURN
*
      END FUNCTION SHgamgam
*
*------------------------------------------------------------------------------
*
      FUNCTION SHgluglu(amh,amt,amb,als,sct,scb) 
*
*     amh       Higgs mass [GeV]      
*     amt       top   mass [GeV]      
*     amb       b     mass [GeV]      
*     als       alpha_s
*
* OUTPUT gAmma(H --> glue glue)
*
      IMPLICIT NONE
*
      REAL*8, parameter :: alpha=137.0359895d0
      REAL*8, parameter :: pi= 3.141592653589793238462643d0
      REAL*8, parameter :: GF= 1.16637d-5
*
      REAL*8 amh,amt,amhs,mhc,amts,amb,ambs,mws,tt,tb,btt,btb,lt,lb,
     #       ftr,fti,fbr,fbi,arg,atr,ati,abr,abi,tw,arw,aiw,sct,scb,scw,
     #       btw,lw,fwr,fwi,ar,ai,arq,aiq,als,rt,rb,rw,gtr,gti,gbr,gbi,
     #       gwr,gwi,fltr,flti,flbr,flbi,flwr,flwi,gltr,glti,glbr,glbi,
     #       glwr,glwi,mzs,cts,sts,tts,cbtr,cbti,cbbr,cbbi,cbwr,cbwi,
     #       i1tr,i1ti,i1br,i1bi,i1wr,i1wi,i2tr,i2ti,i2br,i2bi,i2wr,
     #       cth,sth,tth,i2wi,SHgluglu
*
      amhs= amh*amh
      mhc= amhs*amh
      amts= amt*amt
      ambs= amb*amb
*
      tt= 0.25d0*amhs/amts
*
      IF(tt > 1.d0) THEN
       btt= sqrt(1.d0-1.d0/tt)
       lt= log((1.d0+btt)/(1.d0-btt))
       ftr= -0.25d0*(lt*lt-pi*pi)
       fti= 0.5d0*pi*lt
      ELSE
       arg= sqrt(tt)
       ftr= asin(arg)*asin(arg)
       fti= 0.d0
      ENDIF
*
      atr= 2.d0*(tt+(tt-1)*ftr)/(tt*tt)
      ati= 2.d0*(tt-1)*fti/(tt*tt)
*
      tb= 0.25d0*amhs/ambs
*
      IF(tb > 1.d0) THEN
       btb= sqrt(1.d0-1.d0/tb)
       lb= log((1.d0+btb)/(1.d0-btb))
       fbr= -0.25d0*(lb*lb-pi*pi)
       fbi= 0.5d0*pi*lb
      ELSE
       arg= sqrt(tb)
       fbr= asin(arg)*asin(arg)
       fbi= 0.d0
      ENDIF
*
      abr= 2.d0*(tb+(tb-1)*fbr)/(tb*tb)
      abi= 2.d0*(tb-1)*fbi/(tb*tb)
*
      arq= sct*atr+scb*abr
      aiq= sct*ati+scb*abi
*
      SHgluglu= GF*als*als/64.d0/sqrt(2.d0)/pi**3*
     #     mhc*(arq*arq+aiq*aiq)*(1.d0+(95.d0/4.d0-35.d0/6.d0)*als)
*
*
      RETURN
*
      END FUNCTION SHgluglu
*
*-----------------------------------------------------------------------------
*
      FUNCTION SHZgam(amh,amt,amb,sct,scb,scw) 
*
*     amh       Higgs mass [GeV]      
*     amt       top   mass [GeV]      
*     amb       b     mass [GeV]      
*
* OUTPUT gAmma(H --> Z gamma)
*
      IMPLICIT NONE
*
      REAL*8, parameter :: alpha=137.0359895d0
      REAL*8, parameter :: pi= 3.141592653589793238462643d0
      REAL*8, parameter :: mw= 80.398d0
      REAL*8, parameter :: mz= 91.17667d0
      REAL*8, parameter :: GF= 1.16637d-5
*
      REAL*8 amh,amt,amhs,mhc,amts,amb,ambs,mws,tt,tb,btt,btb,lt,lb,
     #       ftr,fti,fbr,fbi,arg,atr,ati,abr,abi,tw,arw,aiw,sct,scb,scw,
     #       btw,lw,fwr,fwi,ar,ai,arq,aiq,als,rt,rb,rw,gtr,gti,gbr,gbi,
     #       gwr,gwi,fltr,flti,flbr,flbi,flwr,flwi,gltr,glti,glbr,glbi,
     #       glwr,glwi,mzs,cts,sts,tts,cbtr,cbti,cbbr,cbbi,cbwr,cbwi,
     #       i1tr,i1ti,i1br,i1bi,i1wr,i1wi,i2tr,i2ti,i2br,i2bi,i2wr,
     #       cth,sth,tth,i2wi,alz,SHZgam
*
      amhs= amh*amh
      mhc= amhs*amh
      amts= amt*amt
      ambs= amb*amb
*
      tt= 0.25d0*amhs/amts
*
      IF(tt > 1.d0) THEN
       btt= sqrt(1.d0-1.d0/tt)
       lt= log((1.d0+btt)/(1.d0-btt))
       ftr= -0.25d0*(lt*lt-pi*pi)
       fti= 0.5d0*pi*lt
      ELSE
       arg= sqrt(tt)
       ftr= asin(arg)*asin(arg)
       fti= 0.d0
      ENDIF
*
      tb= 0.25d0*amhs/ambs
*
      IF(tb > 1.d0) THEN
       btb= sqrt(1.d0-1.d0/tb)
       lb= log((1.d0+btb)/(1.d0-btb))
       fbr= -0.25d0*(lb*lb-pi*pi)
       fbi= 0.5d0*pi*lb
      ELSE
       arg= sqrt(tb)
       fbr= asin(arg)*asin(arg)
       fbi= 0.d0
      ENDIF
*
      mws= mw*mw
      tw= 0.25d0*amhs/mws
*
      IF(tw > 1.d0) THEN
       btw= sqrt(1.d0-1.d0/tw)
       lw= log((1.d0+btw)/(1.d0-btw))
       fwr= -0.25d0*(lw*lw-pi*pi)
       fwi= 0.5d0*pi*lw
      ELSE
       arg= sqrt(tw)
       fwr= asin(arg)*asin(arg)
       fwi= 0.d0
      ENDIF
*
*---
*
      IF(tt < 1.d0) THEN
       btt= sqrt(1.d0/tt-1.d0)
       arg= sqrt(tt)
       gtr= btt*asin(arg)
       gti= 0.d0
      ELSE
       btt= sqrt(1.d0-1.d0/tt)
       lt= log((1.d0+btt)/(1.d0-btt))
       gtr= 0.5d0*btt*lt
       gti= -0.5d0*pi*btt
      ENDIF
*
      IF(tb < 1.d0) THEN
       btb= sqrt(1.d0/tb-1.d0)
       arg= sqrt(tb)
       gbr= btb*asin(arg)
       gbi= 0.d0
      ELSE
       btb= sqrt(1.d0-1.d0/tb)
       lb= log((1.d0+btb)/(1.d0-btb))
       gbr= 0.5d0*btb*lb
       gbi= -0.5d0*pi*btb
      ENDIF
*
      IF(tw < 1.d0) THEN
       btw= sqrt(1.d0/tw-1.d0)
       arg= sqrt(tw)
       gwr= btw*asin(arg)
       gwi= 0.d0
      ELSE
       btw= sqrt(1.d0-1.d0/tw)
       lw= log((1.d0+btw)/(1.d0-btw))
       gwr= 0.5d0*btw*lw
       gwi= -0.5d0*pi*btw
      ENDIF
*
      mzs= mz*mz
*
      rt= 0.25d0*mzs/amts
      rb= 0.25d0*mzs/ambs
      rw= 0.25d0*mzs/mws
*
      IF(rt > 1.d0) THEN
       btt= sqrt(1.d0-1.d0/rt)
       lt= log((1.d0+btt)/(1.d0-btt))
       fltr= -0.25d0*(lt*lt-pi*pi)
       flti= 0.5d0*pi*lt
      ELSE
       arg= sqrt(rt)
       fltr= asin(arg)*asin(arg)
       flti= 0.d0
      ENDIF
*      
      IF(rb > 1.d0) THEN
       btb= sqrt(1.d0-1.d0/rb)
       lb= log((1.d0+btb)/(1.d0-btb))
       flbr= -0.25d0*(lb*lb-pi*pi)
       flbi= 0.5d0*pi*lb
      ELSE
       arg= sqrt(rb)
       flbr= asin(arg)*asin(arg)
       flbi= 0.d0
      ENDIF
*
      IF(rw > 1.d0) THEN
       btw= sqrt(1.d0-1.d0/rw)
       lw= log((1.d0+btw)/(1.d0-btw))
       flwr= -0.25d0*(lw*lw-pi*pi)
       flwi= 0.5d0*pi*lw
      ELSE
       arg= sqrt(rw)
       flwr= asin(arg)*asin(arg)
       flwi= 0.d0
      ENDIF
*
      IF(rt < 1.d0) THEN
       btt= sqrt(1.d0/rt-1.d0)
       arg= sqrt(rt)
       gltr= btt*asin(arg)
       glti= 0.d0
      ELSE
       btt= sqrt(1.d0-1.d0/rt)
       lt= log((1.d0+btt)/(1.d0-btt))
       gltr= 0.5d0*btt*lt
       glti= -0.5d0*pi*btt
      ENDIF
*
      IF(rb < 1.d0) THEN
       btb= sqrt(1.d0/rb-1.d0)
       arg= sqrt(rb)
       glbr= btb*asin(arg)
       glbi= 0.d0
      ELSE
       btb= sqrt(1.d0-1.d0/rb)
       lb= log((1.d0+btb)/(1.d0-btb))
       glbr= 0.5d0*btb*lb
       glbi= -0.5d0*pi*btb
      ENDIF
*
      IF(rw < 1.d0) THEN
       btw= sqrt(1.d0/rw-1.d0)
       arg= sqrt(rw)
       glwr= btw*asin(arg)
       glwi= 0.d0
      ELSE
       btw= sqrt(1.d0-1.d0/rw)
       lw= log((1.d0+btw)/(1.d0-btw))
       glwr= 0.5d0*btw*lw
       glwi= -0.5d0*pi*btw
      ENDIF
*
      i1tr= -0.5d0/(tt-rt)+0.5d0/(tt-rt)**2*(ftr-fltr)+
     #      rt/(tt-rt)**2*(gtr-gltr)
      i1ti= 0.5d0/(tt-rt)**2*(fti-flti)+rt/(tt-rt)**2*(gti-glti)
*
      i1br= -0.5d0/(tb-rb)+0.5d0/(tb-rb)**2*(fbr-flbr)+
     #      rb/(tb-rb)**2*(gbr-glbr)
      i1bi= 0.5d0/(tb-rb)**2*(fbi-flbi)+rb/(tb-rb)**2*(gbi-glbi)
*
      i1wr= -0.5d0/(tw-rw)+0.5d0/(tw-rw)**2*(fwr-flwr)+
     #      rw/(tw-rw)**2*(gwr-glwr)
      i1wi= 0.5d0/(tw-rw)**2*(fwi-flwi)+rw/(tw-rw)**2*(gwi-glwi)
*
      i2tr= 0.5d0/(tt-rt)*(ftr-fltr)
      i2ti= 0.5d0/(tt-rt)*(fti-flti)
*
      i2br= 0.5d0/(tb-rb)*(fbr-flbr)
      i2bi= 0.5d0/(tb-rb)*(fbi-flbi)
*
      i2wr= 0.5d0/(tw-rw)*(fwr-flwr)
      i2wi= 0.5d0/(tw-rw)*(fwi-flwi)
*
      cbtr= i1tr-i2tr
      cbti= i1ti-i2ti
*
      cbbr= i1br-i2br
      cbbi= i1bi-i2bi
*
      cts= mws/mzs
      sts= 1.d0-cts
      cth= sqrt(cts) 
      sth= sqrt(sts)
      tts= sts/cts
      tth= sqrt(tts)
*  
      cbwr= 1.d0/tth*(4.d0*(3.d0-tts))*i2wr+((1.d0+2.d0*tw)*tts-
     #      (5.d0+2.d0*tw))*i1wr
      cbwi= 1.d0/tth*(4.d0*(3.d0-tts))*i2wi+((1.d0+2.d0*tw)*tts-
     #      (5.d0+2.d0*tw))*i1wi
*
      ar= sct*(2.d0*(1.d0-8.d0/3.d0*sts)/(sth*cth)*cbtr)
      ar= ar+scb*((1.d0-4.d0/3.d0*sts)/(sth*cth)*cbbr)
      ar= ar+scw*cbwr
*
      ai= sct*(2.d0*(1.d0-8.d0/3.d0*sts)/(sth*cth)*cbti)
      ai= ai+scb*((1.d0-4.d0/3.d0*sts)/(sth*cth)*cbbi)
      ai= ai+scw*cbwi
*
      SHZgam= GF/(alpha*alpha)/(64.d0*sqrt(2.d0)*pi**3)*mhc*
     #        (1.d0-mzs/amhs)**3*(ar*ar+ai*ai)
*
      RETURN
*
      END FUNCTION SHZgam

