program SIGMA_ZZ

       implicit     NONE

c  constants

       real           alpha
       parameter      (alpha = 1./137.)
       real           G_mu
       parameter      (G_mu =  1.16639e-5)   ! GeV-2
       real           hbarc
       parameter      (hbarc = 0.1973)       ! GeV-fm
       real           M_Z
       parameter      (M_Z = 91.2)
       real           pi
       parameter      (pi = 3.14159265)
       real           sw2
       parameter      (sw2 = 0.23)

c  local variables

       integer        i
       integer        num_points

       real           Ecm, Ecm_min, Ecm_max
       real           s
       real           sigma
       real           x_Z

       real           G4_tpu, L4, s1, s2, ae, ve, lambda

c  begin

       ae = 1.
       ve = 1. - 4.*sw2

       num_points = 101
       Ecm_min = 150.
       Ecm_max = 250.

       do i = 1, num_points

         Ecm = (Ecm_max - Ecm_min)*FLOAT(i-1)/FLOAT(num_points-1) 
     &         + Ecm_min

c  cross section from Bardin et al hep-ph/9406340

           s = Ecm**2
           s1 = M_Z**2
           s2 = M_Z**2
           lambda = s**2 + s1**2 + s2**2 - 2.*s*s1 - 2.*s1*s2 - 2.*s2*s

           if ( lambda .le. 0. ) then

             sigma = 0.

           else

             L4 = (1./sqrt(lambda))  * LOG ( 
     &       ( s - s1 - s2 + sqrt(lambda) ) / 
     &       ( s - s1 - s2 - sqrt(lambda) ) )
             G4_tpu = ( SQRT(lambda) / s ) * ( 
     &       ( ( s**2 + (s1 + s2)**2 ) / (s - s1 - s2 ) ) * L4
     &       - 2. )
             sigma = ( (G_mu * M_Z**2)**2 / (64.*pi*s) ) * 
     &       (ve**4 + 6.*(ve*ae)**2 + ae**4) * G4_tpu

             sigma = sigma * hbarc**2 * 1.e10

           endif

         write (*, *) Ecm, sigma

       end do

       stop
       END