      program home2000
      implicit none
      integer i,ieh,iteta,iteta1,iphi1,npoints,npointseh
      integer intmh,centenas,dezenas,unidades,iwrite
      real*8 ga,gv
      real*8 eh,teta,teta1,phi1
      real*8 eh1,eh2,ehpole
      real*8 sig_s1,sig_s2,sig_s3,sig_t,sig_t2,sig_st1,sig_st2
      real*8 dsig_s,dsig_sa,dsig_sb,dsig_t,dsig_st1,dsig_st2
      real*8 dsigdeh_t,dsigdeh_st,dsig_t2
      real*8 rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2
      real*8 ehmax,tetamin,tetamax
      real*8 teta1min,teta1max,phi1min,phi1max
      real*8 xeh1(64),weh1(64),xteta(64),wteta(64)
      real*8 xeh2(64),weh2(64),xeh3(64),weh3(64)
      real*8 xteta1(64),wteta1(64),xphi1(64),wphi1(64)
      character*3 ext
      character*9 difout
      common/data/rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2
      
      open(10,file='res.out',status='unknown')

      difout(1:2)='mh'
      difout(6:9)='.out'
      
      npoints=6
      npointseh=16
      pi=acos(-1.d0)
      rs=205.d0
      gf=1.16639d-5
      mz=91.187d0
      sw2=0.23d0
      gz=2.5d0
      factor=0.389379d12*gf**3*mz**4/rs**2
      ga=-0.25d0
      gv=-0.25d0+sw2
      c1=ga**4+gv**4+6.d0*ga**2*gv**2
      c2=(ga**2-gv**2)**2


      mh=69.0d0

      do i=1,61
      mh=mh+1.d0
      
      intmh=int(mh)
      centenas=int(intmh/100)
      dezenas=int((intmh-100*centenas)/10)
      unidades=intmh-centenas*100-dezenas*10
      ext(1:1)=char(centenas+48)
      ext(2:2)=char(dezenas+48)
      ext(3:3)=char(unidades+48)
      difout(3:5)=ext(1:3)

      iwrite=0
      if(abs(int(mh/10.d0)-mh/10.d0).lt.0.001d0)then
      iwrite=1
      open(40,file=difout,status='unknown')
      endif



***   Nota: A integracao em eh e partida em tres intervalos por causa da
***         ressonancia dum dos propagadores do Z. As integracoes angulares
***         nao tem problema.

      
      ehmax=(rs**2+mh**2)/(2.d0*rs)
      ehpole=((rs**2+mh**2-mz**2)/(2.d0*rs))
      eh1=max(mh,ehpole-gz)
      eh2=min(max(mh,ehpole+gz),ehmax)


      tetamin=0.d0
      tetamax=pi
      teta1min=0.d0
      teta1max=pi
      phi1min=0.d0
      phi1max=2.d0*pi


      call dsetgauss(mh,eh1,npointseh,xeh1,weh1)
      call dsetgauss(eh1,eh2,npointseh,xeh2,weh2)
      call dsetgauss(eh2,ehmax,npointseh,xeh3,weh3)
      call dsetgauss(tetamin,tetamax,npoints,xteta,wteta)
      call dsetgauss(teta1min,teta1max,npoints,xteta1,wteta1)
      call dsetgauss(phi1min,phi1max,npoints,xphi1,wphi1)

      sig_s1=0.d0
      sig_s2=0.d0
      sig_s3=0.d0
      sig_t=0.d0
      sig_t2=0.d0
      sig_st1=0.d0
      sig_st2=0.d0
      if(eh1.eq.mh)goto 20


      do ieh=1,npointseh
         eh=xeh1(ieh)
         sig_s3=sig_s3+dsig_sb(eh)*weh1(ieh)
         dsigdeh_t=0.d0
         dsigdeh_st=0.d0
         do iteta=1,npoints
            teta=xteta(iteta)
            sig_s2=sig_s2+dsig_sa(eh,teta)*weh1(ieh)*wteta(iteta)
            do iteta1=1,npoints
               teta1=xteta1(iteta1)
               do iphi1=1,npoints
                  phi1=xphi1(iphi1)
                  sig_s1=sig_s1+dsig_s(eh,teta,teta1,phi1)*weh1(ieh)*
     #                wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  sig_t=sig_t+dsig_t(eh,teta,teta1,phi1)*weh1(ieh)*
     #                wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  sig_t2=sig_t2+dsig_t2(eh,teta,teta1,phi1)*weh1(ieh)*
     #                wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  sig_st1=sig_st1+dsig_st1(eh,teta,teta1,phi1)*
     #            weh1(ieh)*wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  sig_st2=sig_st2+dsig_st2(eh,teta,teta1,phi1)*
     #            weh1(ieh)*wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  if(iwrite.eq.1)then
                  dsigdeh_t=dsigdeh_t+dsig_t(eh,teta,teta1,phi1)*
     #                wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  dsigdeh_st=dsigdeh_st+dsig_st1(eh,teta,teta1,phi1)*
     #                wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  endif
               end do
            end do
         end do
         if(iwrite.eq.1)then
         write(40,*)mh,eh,dsig_sb(eh),dsigdeh_t,dsigdeh_st
         endif
      end do

 20   if(eh2.eq.mh)goto 30

      do ieh=1,npointseh
         eh=xeh2(ieh)
         dsigdeh_t=0.d0
         dsigdeh_st=0.d0
         sig_s3=sig_s3+dsig_sb(eh)*weh2(ieh)
         do iteta=1,npoints
            teta=xteta(iteta)
            sig_s2=sig_s2+dsig_sa(eh,teta)*weh2(ieh)*wteta(iteta)
            do iteta1=1,npoints
               teta1=xteta1(iteta1)
               do iphi1=1,npoints
                  phi1=xphi1(iphi1)
                  sig_s1=sig_s1+dsig_s(eh,teta,teta1,phi1)*weh2(ieh)*
     #                wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  sig_t=sig_t+dsig_t(eh,teta,teta1,phi1)*weh2(ieh)*
     #                wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  sig_t2=sig_t2+dsig_t2(eh,teta,teta1,phi1)*weh2(ieh)*
     #                wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  sig_st1=sig_st1+dsig_st1(eh,teta,teta1,phi1)*
     #            weh2(ieh)*wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  sig_st2=sig_st2+dsig_st2(eh,teta,teta1,phi1)*
     #            weh2(ieh)*wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  if(iwrite.eq.1)then
                  dsigdeh_t=dsigdeh_t+dsig_t(eh,teta,teta1,phi1)*
     #                wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  dsigdeh_st=dsigdeh_st+dsig_st1(eh,teta,teta1,phi1)*
     #                wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  endif
               end do
            end do
         end do
         if(iwrite.eq.1)then
         write(40,*)mh,eh,dsig_sb(eh),dsigdeh_t,dsigdeh_st
         endif
      end do

 30   continue

      do ieh=1,npointseh
         eh=xeh3(ieh)
         dsigdeh_t=0.d0
         dsigdeh_st=0.d0
         sig_s3=sig_s3+dsig_sb(eh)*weh3(ieh)
         do iteta=1,npoints
            teta=xteta(iteta)
            sig_s2=sig_s2+dsig_sa(eh,teta)*weh3(ieh)*wteta(iteta)
            do iteta1=1,npoints
               teta1=xteta1(iteta1)
               do iphi1=1,npoints
                  phi1=xphi1(iphi1)
                  sig_s1=sig_s1+dsig_s(eh,teta,teta1,phi1)*weh3(ieh)*
     #                wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  sig_t=sig_t+dsig_t(eh,teta,teta1,phi1)*weh3(ieh)*
     #                wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  sig_t2=sig_t2+dsig_t2(eh,teta,teta1,phi1)*weh3(ieh)*
     #                wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  sig_st1=sig_st1+dsig_st1(eh,teta,teta1,phi1)*
     #            weh3(ieh)*wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  sig_st2=sig_st2+dsig_st2(eh,teta,teta1,phi1)*
     #            weh3(ieh)*wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)

                  if(iwrite.eq.1)then
                  dsigdeh_t=dsigdeh_t+dsig_t(eh,teta,teta1,phi1)*
     #                wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  dsigdeh_st=dsigdeh_st+dsig_st1(eh,teta,teta1,phi1)*
     #                wteta(iteta)*wteta1(iteta1)*wphi1(iphi1)
                  endif
               end do
            end do
         end do
         if(iwrite.eq.1)then
         write(40,*)mh,eh,dsig_sb(eh),dsigdeh_t,dsigdeh_st
         endif
      end do


***   sigz*: Diagrama do canal s
***   sig_s1: Integracao completamente numerica nas 4 variaveis
***   sig_s2: Integracao analitica em duas e numerica em 2 variaveis(eh e teta)
***   sig_s3: Integracao analitica em tres e numerica em eh.
***   sig_t:  Diagrama do canal t
***   sig_t2:  Diagrama do canal t ERRADO para comparacao
***   sig_st1:  Interferencia so com larguras no segundo Z
***   sig_st2:  Interferencia com larguras em todos os propagadores

      write(10,99)mh,sig_s1,sig_s3,sig_t,sig_t2,sig_st1,sig_st2



      close(40)

      end do
 99   format(f4.0,6(e11.5,1x))

      end



      real*8 function dsig_s(eh,teta,teta1,phi1)
      implicit none
      real*8 eh,teta,teta1,phi1
      real*8 p1(4),p2(4),q1(4),q2(4)
      real*8 c_s,delta0,m12,gama,g,beta,b,ct,st,ct1,st1,cp1,sp1
      real*8 sp,dzs2,dzsp2,dot,ph
      real*8 rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2
      common/data/rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2


      
      c_s=factor*32.d0*mz**4/sqrt(2.d0)/pi**4/8.d0
      dzs2=(rs**2-mz**2)**2
      sp=rs**2+mh**2-2.d0*rs*eh
      dzsp2=(sp-mz**2)**2+mz**2*gz**2
      
           
      p1(4)=rs/2.d0
      p1(1)=0.d0
      p1(2)=0.d0
      p1(3)=rs/2.d0

      p2(4)=rs/2.d0
      p2(1)=0.d0
      p2(2)=0.d0
      p2(3)=-rs/2.d0
      
      delta0=m12(eh)
      g=gama(eh)
      b=beta(eh)
      ct=cos(teta)
      st=sin(teta)
      ct1=cos(teta1)
      st1=sin(teta1)
      sp1=sin(phi1)
      cp1=cos(phi1)

      ph=sqrt(eh**2-mh**2)

      q1(4)=delta0/2.0d0*(g+b*g*ct1)
      q1(1)=delta0/2.d0*(-ct*st1*cp1-st*(g*ct1+g*b))
      q1(2)=delta0/2.d0*st1*sp1
      q1(3)=delta0/2.d0*(st*st1*cp1-ct*(g*ct1+g*b))
      
      q2(4)=delta0/2.0d0*(g-b*g*ct1)
      q2(1)=delta0/2.d0*(ct*st1*cp1-st*(-g*ct1+g*b))
      q2(2)=delta0/2.d0*(-st1*sp1)
      q2(3)=delta0/2.d0*(-st*st1*cp1-ct*(-g*ct1+g*b))
      

      dsig_s=c_s/dzs2/dzsp2*ph*sin(teta)*sin(teta1)*
     #        (c1*dot(p1,q2)*dot(p2,q1)+c2*dot(p1,q1)*dot(p2,q2))


      return
      end


      real*8 function dsig_sa(eh,teta)
      implicit none
      real*8 eh,teta
      real*8 c_s
      real*8 sp,dzs2,dzsp2,ph
      real*8 rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2
      common/data/rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2

      c_s=factor*32.d0*mz**4/sqrt(2.d0)/pi**4*pi/24.d0
      dzs2=(rs**2-mz**2)**2
      sp=rs**2+mh**2-2.d0*rs*eh
      dzsp2=(sp-mz**2)**2+mz**2*gz**2
      
      ph=sqrt(eh**2-mh**2)


      dsig_sa=c_s/dzs2/dzsp2*ph*sin(teta)*(c1+c2)*
     #         rs**2*(rs**2+mh**2/2.d0+eh**2/2.d0-2.d0*rs*eh
     #         -ph**2/2.d0*cos(teta)**2)

      return
      end

      real*8 function dsig_sb(eh)
      implicit none
      real*8 eh
      real*8 c_s
      real*8 sp,dzs2,dzsp2,ph
      real*8 rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2
      common/data/rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2

      c_s=factor*32.d0*mz**4/sqrt(2.d0)/pi**4*pi/24.d0
      dzs2=(rs**2-mz**2)**2
      sp=rs**2+mh**2-2.d0*rs*eh
      dzsp2=(sp-mz**2)**2+mz**2*gz**2
      
      ph=sqrt(eh**2-mh**2)


      dsig_sb=c_s/dzs2/dzsp2*ph*rs**2*(c1+c2)*
     #          (2.d0*rs**2+mh**2+eh**2-4.d0*rs*eh-ph**2/3.d0)


      return
      end


      real*8 function dsig_t(eh,teta,teta1,phi1)
      implicit none
      real*8 eh,teta,teta1,phi1
      real*8 p1(4),p2(4),q1(4),q2(4)
      real*8 c_t,delta0,m12,gama,g,beta,b,ct,st,ct1,st1,cp1,sp1
      real*8 d11,d22,dot,ph
      real*8 rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2
      common/data/rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2


      
      c_t=factor*32.d0*mz**4/sqrt(2.d0)/pi**4/8.d0
      
           
      p1(4)=rs/2.d0
      p1(1)=0.d0
      p1(2)=0.d0
      p1(3)=rs/2.d0

      p2(4)=rs/2.d0
      p2(1)=0.d0
      p2(2)=0.d0
      p2(3)=-rs/2.d0
      
      delta0=m12(eh)
      g=gama(eh)
      b=beta(eh)
      ct=cos(teta)
      st=sin(teta)
      ct1=cos(teta1)
      st1=sin(teta1)
      sp1=sin(phi1)
      cp1=cos(phi1)

      ph=sqrt(eh**2-mh**2)

      q1(4)=delta0/2.0d0*(g+b*g*ct1)
      q1(1)=delta0/2.d0*(-ct*st1*cp1-st*(g*ct1+g*b))
      q1(2)=delta0/2.d0*st1*sp1
      q1(3)=delta0/2.d0*(st*st1*cp1-ct*(g*ct1+g*b))
      
      q2(4)=delta0/2.0d0*(g-b*g*ct1)
      q2(1)=delta0/2.d0*(ct*st1*cp1-st*(-g*ct1+g*b))
      q2(2)=delta0/2.d0*(-st1*sp1)
      q2(3)=delta0/2.d0*(-st*st1*cp1-ct*(-g*ct1+g*b))
      
      d11=mz**2+2.d0*dot(p1,q1)
      d22=mz**2+2.d0*dot(p2,q2)

      dsig_t=c_t/d11**2/d22**2*ph*sin(teta)*sin(teta1)*
     #        (c1*dot(p1,q2)*dot(p2,q1)+c2*dot(p1,p1)*dot(q1,q2))

      return
      end

      real*8 function dsig_t2(eh,teta,teta1,phi1)
      implicit none
      real*8 eh,teta,teta1,phi1
      real*8 p1(4),p2(4),q1(4),q2(4)
      real*8 c_t,delta0,m12,gama,g,beta,b,ct,st,ct1,st1,cp1,sp1
      real*8 d11,d22,dot,ph,sp
      real*8 rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2
      common/data/rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2


      
      c_t=factor*32.d0*mz**4/sqrt(2.d0)/pi**4/8.d0
      sp=rs**2+mh**2-2.d0*rs*eh      
           
      p1(4)=rs/2.d0
      p1(1)=0.d0
      p1(2)=0.d0
      p1(3)=rs/2.d0

      p2(4)=rs/2.d0
      p2(1)=0.d0
      p2(2)=0.d0
      p2(3)=-rs/2.d0
      
      delta0=m12(eh)
      g=gama(eh)
      b=beta(eh)
      ct=cos(teta)
      st=sin(teta)
      ct1=cos(teta1)
      st1=sin(teta1)
      sp1=sin(phi1)
      cp1=cos(phi1)

      ph=sqrt(eh**2-mh**2)

      q1(4)=delta0/2.0d0*(g+b*g*ct1)
      q1(1)=delta0/2.d0*(-ct*st1*cp1-st*(g*ct1+g*b))
      q1(2)=delta0/2.d0*st1*sp1
      q1(3)=delta0/2.d0*(st*st1*cp1-ct*(g*ct1+g*b))
      
      q2(4)=delta0/2.0d0*(g-b*g*ct1)
      q2(1)=delta0/2.d0*(ct*st1*cp1-st*(-g*ct1+g*b))
      q2(2)=delta0/2.d0*(-st1*sp1)
      q2(3)=delta0/2.d0*(-st*st1*cp1-ct*(-g*ct1+g*b))
      
      d11=mz**2+2.d0*dot(p1,q1)
      d22=mz**2+2.d0*dot(p2,q2)

      dsig_t2=c_t/d11**2/d22**2*ph*sin(teta)*sin(teta1)*
     #       rs**2*(c1/24.d0*(2.d0*rs**2+mh**2+eh**2-4.d0*rs*eh
     #       -ph**2*cos(teta)**2)+c2/4.d0*sp)

      return
      end


      real*8 function dsig_st1(eh,teta,teta1,phi1)
      implicit none
      real*8 eh,teta,teta1,phi1
      real*8 p1(4),p2(4),q1(4),q2(4)
      real*8 c_st,delta0,m12,gama,g,beta,b,ct,st,ct1,st1,cp1,sp1
      real*8 sp,dzs,dot,ph,d11,d22,realpropzsp
      real*8 rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2
      common/data/rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2


      
      c_st=factor*64.d0*mz**4/sqrt(2.d0)/pi**4/8.d0
      dzs=(rs**2-mz**2)
      sp=rs**2+mh**2-2.d0*rs*eh
      realpropzsp=(sp-mz**2)/((sp-mz**2)**2+mz**2*gz**2)
      
           
      p1(4)=rs/2.d0
      p1(1)=0.d0
      p1(2)=0.d0
      p1(3)=rs/2.d0

      p2(4)=rs/2.d0
      p2(1)=0.d0
      p2(2)=0.d0
      p2(3)=-rs/2.d0
      
      delta0=m12(eh)
      g=gama(eh)
      b=beta(eh)
      ct=cos(teta)
      st=sin(teta)
      ct1=cos(teta1)
      st1=sin(teta1)
      sp1=sin(phi1)
      cp1=cos(phi1)

      ph=sqrt(eh**2-mh**2)

      q1(4)=delta0/2.0d0*(g+b*g*ct1)
      q1(1)=delta0/2.d0*(-ct*st1*cp1-st*(g*ct1+g*b))
      q1(2)=delta0/2.d0*st1*sp1
      q1(3)=delta0/2.d0*(st*st1*cp1-ct*(g*ct1+g*b))
      
      q2(4)=delta0/2.0d0*(g-b*g*ct1)
      q2(1)=delta0/2.d0*(ct*st1*cp1-st*(-g*ct1+g*b))
      q2(2)=delta0/2.d0*(-st1*sp1)
      q2(3)=delta0/2.d0*(-st*st1*cp1-ct*(-g*ct1+g*b))

      
      d11=mz**2+2.d0*dot(p1,q1)
      d22=mz**2+2.d0*dot(p2,q2)
      


      dsig_st1=c_st/dzs/d11/d22*realpropzsp*ph*sin(teta)*sin(teta1)*
     #        c1*dot(p1,q2)*dot(p2,q1)


      return
      end


      real*8 function dsig_st2(eh,teta,teta1,phi1)
      implicit none
      real*8 eh,teta,teta1,phi1
      real*8 p1(4),p2(4),q1(4),q2(4)
      real*8 c_st,delta0,m12,gama,g,beta,b,ct,st,ct1,st1,cp1,sp1
      real*8 sp,dot,ph
      real*8 rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2
      real*8 a1,a2,a3,a4,b,num,den,realpart
      common/data/rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2


      
      c_st=factor*64.d0*mz**4/sqrt(2.d0)/pi**4/8.d0
      sp=rs**2+mh**2-2.d0*rs*eh

           
      p1(4)=rs/2.d0
      p1(1)=0.d0
      p1(2)=0.d0
      p1(3)=rs/2.d0

      p2(4)=rs/2.d0
      p2(1)=0.d0
      p2(2)=0.d0
      p2(3)=-rs/2.d0
      
      delta0=m12(eh)
      g=gama(eh)
      b=beta(eh)
      ct=cos(teta)
      st=sin(teta)
      ct1=cos(teta1)
      st1=sin(teta1)
      sp1=sin(phi1)
      cp1=cos(phi1)

      ph=sqrt(eh**2-mh**2)

      q1(4)=delta0/2.0d0*(g+b*g*ct1)
      q1(1)=delta0/2.d0*(-ct*st1*cp1-st*(g*ct1+g*b))
      q1(2)=delta0/2.d0*st1*sp1
      q1(3)=delta0/2.d0*(st*st1*cp1-ct*(g*ct1+g*b))
      
      q2(4)=delta0/2.0d0*(g-b*g*ct1)
      q2(1)=delta0/2.d0*(ct*st1*cp1-st*(-g*ct1+g*b))
      q2(2)=delta0/2.d0*(-st1*sp1)
      q2(3)=delta0/2.d0*(-st*st1*cp1-ct*(-g*ct1+g*b))

      
      a1=rs**2-mz**2
      a2=sp-mz**2
      a3=-mz**2-2.d0*dot(p1,q1)
      a4=-mz**2-2.d0*dot(p2,q2)
      b=mz+gz

      num=a1*a2*a3*a4-b**2*(a1*a2+a1*a3+a1*a4+a2*a3+a2*a4+a3*a4)
     #    +b**4
      den=(a1**2+b**2)*(a2**2+b**2)*(a3**2+b**2)*(a4**2+b**2)

      realpart=num/den

      dsig_st2=c_st*realpart*ph*sin(teta)*sin(teta1)*
     #        c1*dot(p1,q2)*dot(p2,q1)


      return
      end




      real*8 function m12(eh)
      implicit none
      real*8 eh
      real*8 rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2
      common/data/rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2
      m12=sqrt(rs**2+mh**2-2.d0*rs*eh)
      return
      end

      real*8 function e12(eh)
      implicit none
      real*8 eh,m12
      real*8 rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2
      common/data/rs,mh,pi,gf,mz,sw2,factor,gz,c1,c2
      e12=(rs**2+m12(eh)**2-mh**2)/2.d0/rs
      return
      end

      real*8 function gama(eh)
      implicit none
      real*8 e12,m12,eh
      gama=e12(eh)/m12(eh)
      return
      end

      real*8 function beta(eh)
      implicit none
      real*8 gama,eh
      beta=sqrt(1.d0-1.d0/gama(eh)**2)
      return
      end

      real*8 function dot(p,q)
      implicit none
      real*8 p(4),q(4)
      dot=p(4)*q(4)-p(1)*q(1)-p(2)*q(2)-p(3)*q(3)
      return
      end



c
c       This subroutine stores the values for
c       the gaussian integration
c       Only works with n=4,6,8,10,16,20,24,32
c
      subroutine dsetgauss(a,b,n,x,w)
      implicit none
      integer i,n,nm,kiplus,kiminus
      real*8 a,b,gauss4(4),gauss6(6),gauss8(8),work(128)
      real*8 gauss10(10),gauss16(16),gauss20(20)
      real*8 gauss24(24),gauss32(32)
      real*8 x(64),w(64),xuni(64),wuni(64)
      data gauss4/0.339981043584856d0,0.652145154862546d0,  
     #0.861136311594053d0,0.347854845137454d0/
      data gauss6/0.238619186083197d0,0.467913934572691,
     #0.661209386466265d0,0.360761573048139,
     #0.932469514203152d0,0.171324492379170d0/
      data gauss8/0.183434642495650d0,0.362683783378362d0,
     #0.525532409916329d0,0.313706645877887d0,
     #0.796666477413627d0,0.222381034453374d0,
     #0.960289856497536d0,0.101228536290376d0/
      data gauss10/0.148874338981631d0,0.295524224714753d0,
     #0.433395394129247d0,0.269266719309996d0,
     #0.679409568299024d0,0.219086362515982d0,
     #0.865063366688985d0,0.149451349150581d0,
     #0.973906528517172d0,0.066671344308688d0/
      data gauss16/0.095012509837637d0,0.189450610455068d0,
     #0.281603550779258d0,0.182603415044923d0,
     #0.458016777657227d0,0.169156519395002d0,
     #0.617876244402643d0,0.149595988816576d0,
     #0.755404408355003d0,0.124628971255533d0,
     #0.865631202387831d0,0.095158511682492d0,
     #0.944575023073232d0,0.062253523938647d0,
     #0.989400934991649d0,0.027152459411754d0/
      data gauss20/0.076526521133497d0,0.152753387130725d0,
     #0.227785851141645d0,0.149172986472603d0,
     #0.373706088715419d0,0.142096109318382d0,
     #0.510867001950827d0,0.131688638449176d0,
     #0.636053680726515d0,0.118194531961518d0,
     #0.746331906460150d0,0.101930119817240d0,
     #0.839116971822218d0,0.083276741576704d0,
     #0.912234428251325d0,0.062672048334109d0,
     #0.963971927277913d0,0.040601429800386d0,
     #0.993128599185094d0,0.017614007139152d0/
      data gauss24/0.064056892862605d0,0.127938195346752d0,
     #0.191118867473616d0,0.125837456346828d0,
     #0.315042679696163d0,0.121670472927803d0,
     #0.433793507626045d0,0.115505668053725d0,
     #0.545421471388839d0,0.107444270115965d0,
     #0.648093651936975d0,0.097618652104113d0,
     #0.740124191578554d0,0.086190161531953d0,
     #0.820001985973902d0,0.073346481411080d0,
     #0.886415527004401d0,0.059298584915436d0,
     #0.938274552002732d0,0.044277438817419d0,
     #0.974728555971309d0,0.028531388628933d0,
     #0.995187219997021d0,0.012341229799987d0/
      data gauss32/0.048307665687738d0,0.096540088514727d0,
     #0.144471961582796d0,0.095638720079274d0,
     #0.239287362252137d0,0.093844399080804d0,
     #0.331868602282127d0,0.091173878695763d0,
     #0.421351276130635d0,0.087652093004403d0,
     #0.506899908932229d0,0.083311924226946d0,
     #0.587715757240762d0,0.078193895787070d0,
     #0.663044266930215d0,0.072345794108848d0,
     #0.732182118740289d0,0.065822222776361d0,
     #0.794483795967942d0,0.058684093478535d0,
     #0.849367613732569d0,0.050998059262376d0,
     #0.896321155766052d0,0.042835898022226d0,
     #0.934906075937739d0,0.034273862913021d0,
     #0.964762255587506d0,0.025392065309262d0,
     #0.985611511545268d0,0.016274394730905d0,
     #0.997263861849481d0,0.007018610009470d0/


      nm=n/2

      if(n.eq.4)then
      do i=1,n
      work(i)=gauss4(i)
      enddo
      endif
      if(n.eq.6)then
      do i=1,n
      work(i)=gauss6(i)
      enddo
      endif
      if(n.eq.8)then
      do i=1,n
      work(i)=gauss8(i)
      enddo
      endif
      if(n.eq.10)then
      do i=1,n
      work(i)=gauss10(i)
      enddo
      endif
      if(n.eq.16)then
      do i=1,n
      work(i)=gauss16(i)
      enddo
      endif
      if(n.eq.20)then
      do i=1,n
      work(i)=gauss20(i)
      enddo
      endif
      if(n.eq.24)then
      do i=1,n
      work(i)=gauss24(i)
      enddo
      endif
      if(n.eq.32)then
      do i=1,n
      work(i)=gauss32(i)
      enddo
      endif



      

      do 10 i=1,nm
      kiplus=nm+i
      kiminus=nm+1-i
      xuni(kiplus)=work(2*i-1)
      xuni(kiminus)=-work(2*i-1)
      wuni(kiplus)=work(2*i)
      wuni(kiminus)=work(2*i)
 10   continue

      do i=1,n
      x(i)=(b-a)/2.d0*xuni(i)+(b+a)/2.
      w(i)=(b-a)/2.d0*wuni(i)
      end do

      return
      end



