subroutine integral(gamma1,gamma2,gamma3,gamma4,ity1,ity2,a1,a2,
     &  si1,si2,si3,si4,transp,q)
      implicit none
      integer ity1,ity2
      integer ilam1,ilam2,ilam3,ilam4,iincr
      double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
     &  lambda2,lambda3,lambda4
      logical transp
      double precision elocal,ele
      double precision delta,delta2,sum,ene,sumene,boltz
      double precision q,a1(2,2),a2(2,2),si1,si2,si3,si4
      double precision conv /.01745329252d0/,pi /3.141592654d0/
      
      iincr=20
      delta=iincr*conv
      delta2=0.5d0*delta
cd      print *,'iincr',iincr,' delta',delta 
cd      write(2,*) gamma1,gamma2,ity1,ity2,a1,a2,si1,si2,si3,si4,transp

cd      do ilam1=-180,180,5
cd        do ilam2=-180,180,5
cd          lambda1=ilam1*conv+delta2
cd          lambda2=ilam2*conv+delta2
cd          write(2,'(2i5,2f10.5)') ilam1,ilam2,elocal(2,lambda1,lambda2),
cd     &    ele(lambda1,lambda2,a1,1.0d0,1.d00)
cd        enddo
cd      enddo
cd      stop

      sum=0.0d0
      sumene=0.0d0
      do ilam1=-180,179,iincr
        do ilam2=-180,179,iincr
          do ilam3=-180,179,iincr
            do ilam4=-180,179,iincr
              lambda1=ilam1*conv+delta2
              lambda2=ilam2*conv+delta2
              lambda3=ilam3*conv+delta2
              lambda4=ilam4*conv+delta2
cd              write (2,*) ilam1,ilam2,ilam3,ilam4
cd              write (2,*) lambda1,lambda2,lambda3,lambda4
              ene=
     &         -elocal(ity1,lambda1,lambda2,.false.)*
     &          elocal(ity2,lambda3,lambda4,transp)*
     &          ele(si1*lambda1+gamma1,si3*lambda3+gamma3,a1)*
     &          ele(si2*lambda2+gamma2,si4*lambda4+gamma4,a2)
cd              write (2,*) elocal(ity1,lambda1,gamma1-pi-lambda2),
cd     &        elocal(ity2,lambda3,gamma2-pi-lambda4),
cd     &        ele(lambda1,lambda2,a1,si1,si3),
cd     &        ele(lambda3,lambda4,a2,si2,si4) 
              sum=sum+ene
            enddo
          enddo
        enddo
      enddo
      q=sum/(2*pi)**4*delta**4
      write (2,* )'sum',sum,' q',q
      return
      end
c---------------------------------------------------------------------------
      subroutine integral3(gamma1,gamma2,ity1,ity2,ity3,ity4,
     &  a1,koniec,q1,q2,q3,q4)
      implicit none
      integer ity1,ity2,ity3,ity4
      integer ilam1,ilam2,ilam3,ilam4,iincr
      double precision gamma1,gamma2,gamma3,gamma4,beta,lambda1,
     &  lambda2,lambda3,lambda4
      logical koniec
      double precision elocal,ele
      double precision delta,delta2,sum1,sum2,sum3,sum4,
     &  ene1,ene2,ene3,ene4,boltz
      double precision q1,q2,q3,q4,a1(2,2),a2(2,2)
      double precision conv /.01745329252d0/,pi /3.141592654d0/
      
      iincr=60
      delta=iincr*conv
      delta2=0.5d0*delta
cd      print *,'iincr',iincr,' delta',delta 
      write(2,*) gamma1,gamma2,ity1,ity2,ity3,ity4,a1,koniec

cd      do ilam1=-180,180,5
cd        do ilam2=-180,180,5
cd          lambda1=ilam1*conv+delta2
cd          lambda2=ilam2*conv+delta2
cd          write(2,'(2i5,2f10.5)') ilam1,ilam2,elocal(2,lambda1,lambda2),
cd     &    ele(lambda1,lambda2,a1,1.0d0,1.d00)
cd        enddo
cd      enddo
cd      stop

      sum1=0.0d0
      sum2=0.0d0
      sum3=0.0d0
      sum4=0.0d0
      do ilam1=-180,179,iincr
        do ilam2=-180,179,iincr
          do ilam3=-180,179,iincr
            do ilam4=-180,179,iincr
              lambda1=ilam1*conv+delta2
              lambda2=ilam2*conv+delta2
              lambda3=ilam3*conv+delta2
              lambda4=ilam4*conv+delta2
cd              write (2,*) ilam1,ilam2,ilam3,ilam4
cd              write (2,*) lambda1,lambda2,lambda3,lambda4
              if (.not.koniec) then
              ene1=
     &          elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)*
     &          elocal(ity3,lambda3,gamma2-pi-lambda4,.false.)*
     &          ele(lambda2,lambda4,a1)
              else
              ene1=
     &          elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)*
     &          elocal(ity3,lambda3,lambda4,.false.)*
     &          ele(lambda2,-lambda4,a1)
              endif
              ene2=
     &          elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)*
     &          elocal(ity4,lambda3,lambda4,.false.)*
     &          ele(lambda2,lambda3,a1)
              if (.not.koniec) then
              ene3=
     &          elocal(ity2,lambda1,lambda2,.false.)*
     &          elocal(ity3,lambda3,gamma2-pi-lambda4,.false.)*
     &          ele(lambda1,lambda4,a1)
              else
              ene3=
     &          elocal(ity2,lambda1,lambda2,.false.)*
     &          elocal(ity3,lambda3,lambda4,.false.)*
     &          ele(lambda1,-lambda4,a1)
              endif
              ene4=
     &          elocal(ity2,lambda1,lambda2,.false.)*
     &          elocal(ity4,lambda3,lambda4,.false.)*
     &          ele(lambda1,lambda3,a1)
              sum1=sum1+ene1
              sum2=sum2+ene2
              sum3=sum3+ene3
              sum4=sum4+ene4
            enddo
          enddo
        enddo
      enddo
      q1=sum1/(2*pi)**4*delta**4
      q2=sum2/(2*pi)**4*delta**4
      q3=sum3/(2*pi)**4*delta**4
      q4=sum4/(2*pi)**4*delta**4
      write (2,* )'sum',sum1,sum2,sum3,sum4,' q',q1,q2,q3,q4
      return
      end
c-------------------------------------------------------------------------
      subroutine integral5(gamma1,gamma2,gamma3,gamma4,ity1,ity2,ity3,
     &  ity4,ity5,ity6,a1,a2,si1,si2,si3,si4,transp,ene1,ene2,ene3,ene4)
      implicit none
      integer ity1,ity2,ity3,ity4,ity5,ity6
      integer ilam1,ilam2,ilam3,ilam4,ilam5,iincr
      double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
     &  lambda2,lambda3,lambda4,lambda5
      logical transp
      double precision elocal,ele
      double precision eloc1,eloc2,eloc3,eloc4,eloc5,eloc6,ele1,ele2
      double precision delta,delta2,sum,ene,sumene,pom
      double precision ene1,ene2,ene3,ene4,sum1,sum2,sum3,sum4,
     &  a1(2,2),a2(2,2)
      integer si1,si2,si3,si4
      double precision conv /.01745329252d0/,pi /3.141592654d0/
      
      iincr=60
      delta=iincr*conv
      delta2=0.5d0*delta
cd      print *,'iincr',iincr,' delta',delta 
cd      write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2,
cd     &  ' gamma3=',gamma3,' gamma4=',gamma4
cd      write(2,*) ity1,ity2,ity3,ity4,ity5,ity6
cd      write(2,*) 'a1=',a1
cd      write(2,*) 'a2=',a2
cd      write(2,*) si1,si2,si3,si4,transp

      sum1=0.0d0
      sum2=0.0d0
      sum3=0.0d0
      sum4=0.0d0
      do ilam1=-180,179,iincr
        do ilam2=-180,179,iincr
          do ilam3=-180,179,iincr
            do ilam4=-180,179,iincr
              do ilam5=-180,179,iincr
                lambda1=ilam1*conv+delta2
                lambda2=ilam2*conv+delta2
                lambda3=ilam3*conv+delta2
                lambda4=ilam4*conv+delta2
                lambda5=ilam5*conv+delta2
                if (transp) then
                  ele1=ele(lambda1,si4*lambda4,a1)
                  ele2=ele(lambda2,lambda3,a2)
                else
                  ele1=ele(lambda1,lambda3,a1)
                  ele2=ele(lambda2,si4*lambda4,a2)
                endif
                eloc2=elocal(ity2,lambda1,gamma2-pi-lambda2,.false.)
                eloc5=elocal(ity5,lambda3,gamma4-pi-si4*lambda4,.false.)
                pom=ele1*ele2*eloc2*eloc5
                if (si1.gt.0) then
                  eloc1=elocal(ity1,lambda5,gamma1-pi-lambda1,.false.)
                  sum1=sum1+pom*eloc1
                endif
                eloc3=elocal(ity3,lambda2,lambda5,.false.)
                sum2=sum2+pom*eloc3
                eloc4=elocal(ity4,lambda5,gamma3-pi-lambda3,.false.)
                sum3=sum3+pom*eloc4
                if (si4.gt.0) then
                  eloc6=elocal(ity6,lambda4,lambda5,.false.)
                  sum4=sum4+pom*eloc6
                endif
              enddo
            enddo
          enddo
        enddo
      enddo
      pom=1.0d0/(2*pi)**5*delta**5
      ene1=sum1*pom
      ene2=sum2*pom
      ene3=sum3*pom
      ene4=sum4*pom 
c      write (2,* )'sum',sum1,sum2,sum3,sum4,' q',ene1,ene2,ene3,ene4
      return
      end
c-------------------------------------------------------------------------
      subroutine integral_turn6(gamma1,gamma2,gamma3,gamma4,ity1,ity2,
     &  ity3,ity4,ity5,ity6,a1,a2,ene_turn6)
      implicit none
      integer ity1,ity2,ity3,ity4,ity5,ity6
      integer ilam1,ilam2,ilam3,ilam4,ilam5,ilam6,iincr
      double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
     &  lambda2,lambda3,lambda4,lambda5,lambda6
      logical transp
      double precision elocal,ele
      double precision eloc1,eloc2,eloc3,eloc4,eloc41,eloc5,eloc6,
     &  eloc61,ele1,ele2
      double precision delta,delta2,sum,ene,sumene,pom,ene5
      double precision ene_turn6,sum5,a1(2,2),a2(2,2)
      double precision conv /.01745329252d0/,pi /3.141592654d0/
      
      iincr=60
      delta=iincr*conv
      delta2=0.5d0*delta
cd      print *,'iincr',iincr,' delta',delta 
      write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2,
     &  ' gamma3=',gamma3,' gamma4=',gamma4
      write(2,*) ity1,ity2,ity3,ity4,ity5,ity6
      write(2,*) 'a1=',a1
      write(2,*) 'a2=',a2

      sum5=0.0d0
      do ilam1=-180,179,iincr
        do ilam2=-180,179,iincr
          do ilam3=-180,179,iincr
            do ilam4=-180,179,iincr
              do ilam5=-180,179,iincr
                lambda1=ilam1*conv+delta2
                lambda2=ilam2*conv+delta2
                lambda3=ilam3*conv+delta2
                lambda4=ilam4*conv+delta2
                lambda5=ilam5*conv+delta2
                ele1=ele(lambda1,-lambda4,a1)
                ele2=ele(lambda2,lambda3,a2)
                eloc2=elocal(ity2,lambda1,gamma2-pi-lambda2,.false.)
                eloc5=elocal(ity5,lambda3,lambda4,.false.)
                pom=ele1*ele2*eloc2*eloc5
                eloc3=elocal(ity3,lambda2,gamma3-pi-lambda5,.false.)
                eloc4=elocal(ity4,lambda5,gamma4-pi-lambda3,.false.)
                sum5=sum5+pom*eloc3*eloc4
              enddo
            enddo
          enddo
        enddo
      enddo
      pom=-1.0d0/(2*pi)**5*delta**5
      ene_turn6=sum5*pom 
c      print *,'sum6',sum6,' ene6',ene6
      return
      end
c-------------------------------------------------------------------------
      subroutine integral6(gamma1,gamma2,gamma3,gamma4,ity1,ity2,ity3,
     &  ity4,ity5,ity6,a1,a2,si1,si2,si3,si4,transp,ene1,ene2,ene3,ene4,
     &  ene5,ene6)
      implicit none
      integer ity1,ity2,ity3,ity4,ity5,ity6
      integer ilam1,ilam2,ilam3,ilam4,ilam5,ilam6,iincr
      double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
     &  lambda2,lambda3,lambda4,lambda5,lambda6
      logical transp
      double precision elocal,ele
      double precision eloc1,eloc2,eloc3,eloc4,eloc41,eloc5,eloc6,
     &  eloc61,ele1,ele2
      double precision delta,delta2,sum,ene,sumene,pom
      double precision ene1,ene2,ene3,ene4,ene5,ene6,sum1,sum2,sum3,
     &  sum4,sum5,sum6,a1(2,2),a2(2,2)
      integer si1,si2,si3,si4
      double precision conv /.01745329252d0/,pi /3.141592654d0/
      
      iincr=60
      delta=iincr*conv
      delta2=0.5d0*delta
cd      print *,'iincr',iincr,' delta',delta 
cd      write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2,
cd     &  ' gamma3=',gamma3,' gamma4=',gamma4
cd      write(2,*) ity1,ity2,ity3,ity4,ity5,ity6
cd      write(2,*) 'a1=',a1
cd      write(2,*) 'a2=',a2
cd      write(2,*) si1,si2,si3,si4,transp

      sum1=0.0d0
      sum2=0.0d0
      sum3=0.0d0
      sum4=0.0d0
      sum5=0.0d0
      sum6=0.0d0
      eloc1=0.0d0
      eloc6=0.0d0
      eloc61=0.0d0
      do ilam1=-180,179,iincr
        do ilam2=-180,179,iincr
          do ilam3=-180,179,iincr
            do ilam4=-180,179,iincr
              do ilam5=-180,179,iincr
                do ilam6=-180,179,iincr
                lambda1=ilam1*conv+delta2
                lambda2=ilam2*conv+delta2
                lambda3=ilam3*conv+delta2
                lambda4=ilam4*conv+delta2
                lambda5=ilam5*conv+delta2
                lambda6=ilam6*conv+delta2
                if (transp) then
                  ele1=ele(lambda1,si4*lambda4,a1)
                  ele2=ele(lambda2,lambda3,a2)
                else
                  ele1=ele(lambda1,lambda3,a1)
                  ele2=ele(lambda2,si4*lambda4,a2)
                endif
                eloc2=elocal(ity2,lambda1,gamma2-pi-lambda2,.false.)
                eloc5=elocal(ity5,lambda3,gamma4-pi-si4*lambda4,.false.)
                pom=ele1*ele2*eloc2*eloc5
                if (si1.gt.0) then
                  eloc1=elocal(ity1,lambda5,gamma1-pi-lambda1,.false.)
                endif
                eloc3=elocal(ity3,lambda2,lambda6,.false.)
                sum1=sum1+pom*eloc1*eloc3
                eloc4=elocal(ity4,lambda5,gamma3-pi-lambda3,.false.)
                if (si4.gt.0) then
                  eloc6=elocal(ity6,lambda4,lambda6,.false.)
                  eloc61=elocal(ity6,lambda4,lambda5,.false.)
                endif
                sum2=sum2+pom*eloc4*eloc6
                eloc41=elocal(ity4,lambda6,gamma3-pi-lambda3,.false.)
                sum3=sum3+pom*eloc1*eloc41
                sum4=sum4+pom*eloc1*eloc6
                sum5=sum5+pom*eloc3*eloc4
                sum6=sum6+pom*eloc3*eloc61
                enddo
              enddo
            enddo
          enddo
        enddo
      enddo
      pom=-1.0d0/(2*pi)**6*delta**6
      ene1=sum1*pom
      ene2=sum2*pom
      ene3=sum3*pom
      ene4=sum4*pom 
      ene5=sum5*pom 
      ene6=sum6*pom 
c      print *,'sum6',sum6,' ene6',ene6
      return
      end
c-------------------------------------------------------------------------
      subroutine integral3a(gamma1,gamma2,ity1,ity2,a1,si1,ene1)
      implicit none
      integer ity1,ity2,ity3,ity4,ity5,ity6
      integer ilam1,ilam2,ilam3,ilam4,ilam5,ilam6,iincr
      double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
     &  lambda2,lambda3,lambda4,lambda5,lambda6
      logical transp
      double precision elocal,ele
      double precision eloc1,eloc2,eloc3,eloc4,eloc41,eloc5,eloc6,
     &  eloc61,ele1,ele2
      double precision delta,delta2,sum,ene,sumene,pom
      double precision ene1,ene2,ene3,ene4,ene5,ene6,sum1,sum2,sum3,
     &  sum4,sum5,sum6,a1(2,2),a2(2,2)
      integer si1,si2,si3,si4
      double precision conv /.01745329252d0/,pi /3.141592654d0/
      
      iincr=60
      delta=iincr*conv
      delta2=0.5d0*delta
cd      print *,'iincr',iincr,' delta',delta 
cd      write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2
cd      write(2,*) ity1,ity2
cd      write(2,*) 'a1=',a1
cd      write(2,*) si1,

      sum1=0.0d0
      eloc1=0.0d0
      do ilam1=-180,179,iincr
        do ilam2=-180,179,iincr
          do ilam3=-180,179,iincr
            lambda1=ilam1*conv+delta2
            lambda2=ilam2*conv+delta2
            lambda3=ilam3*conv+delta2
            ele1=ele(lambda1,si1*lambda3,a1)
            eloc1=elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)
            if (si1.gt.0) then
              eloc2=elocal(ity2,lambda2,gamma2-pi-lambda3,.false.)
            else
              eloc2=elocal(ity2,lambda2,lambda3,.false.)
            endif
            sum1=sum1+ele1*eloc1*eloc2
          enddo
        enddo
      enddo
      pom=1.0d0/(2*pi)**3*delta**3
      ene1=sum1*pom
      return
      end
c-------------------------------------------------------------------------
      subroutine integral4a(gamma1,gamma2,gamma3,ity1,ity2,ity3,a1,si1,
     &  ene1)
      implicit none
      integer ity1,ity2,ity3,ity4,ity5,ity6
      integer ilam1,ilam2,ilam3,ilam4,ilam5,ilam6,iincr
      double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
     &  lambda2,lambda3,lambda4,lambda5,lambda6
      logical transp
      double precision elocal,ele
      double precision eloc1,eloc2,eloc3,eloc4,eloc41,eloc5,eloc6,
     &  eloc61,ele1,ele2
      double precision delta,delta2,sum,ene,sumene,pom
      double precision ene1,ene2,ene3,ene4,ene5,ene6,sum1,sum2,sum3,
     &  sum4,sum5,sum6,a1(2,2),a2(2,2)
      integer si1,si2,si3,si4
      double precision conv /.01745329252d0/,pi /3.141592654d0/
      
      iincr=60
      delta=iincr*conv
      delta2=0.5d0*delta
cd      print *,'iincr',iincr,' delta',delta 
cd      write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2,
cd     &  ' gamma3=',gamma3
cd      write(2,*) ity1,ity2,ity3
cd      write(2,*) 'a1=',a1
cd      write(2,*) 'si1=',si1
      sum1=0.0d0
      do ilam1=-180,179,iincr
        do ilam2=-180,179,iincr
          do ilam3=-180,179,iincr
            do ilam4=-180,179,iincr
              lambda1=ilam1*conv+delta2
              lambda2=ilam2*conv+delta2
              lambda3=ilam3*conv+delta2
              lambda4=ilam4*conv+delta2
              ele1=ele(lambda1,si1*lambda4,a1)
              eloc1=elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)
              eloc2=elocal(ity2,lambda2,gamma2-pi-lambda3,.false.)
              if (si1.gt.0) then
                eloc3=elocal(ity3,lambda3,gamma3-pi-lambda4,.false.)
              else
                eloc3=elocal(ity3,lambda3,lambda4,.false.)
              endif
              sum1=sum1+ele1*eloc1*eloc2*eloc3
            enddo
          enddo
        enddo
      enddo
      pom=-1.0d0/(2*pi)**4*delta**4
      ene1=sum1*pom
      return
      end
c-------------------------------------------------------------------------
      double precision function elocal(i,x,y,transp)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.TORSION'
      integer i
      double precision x,y,u(2),v(2),cu(2),dv(2),ev(2) 
      double precision scalar2
      logical transp
      u(1)=dcos(x)
      u(2)=dsin(x)
      v(1)=dcos(y)
      v(2)=dsin(y) 
      if (transp) then
        call matvec2(cc(1,1,i),v,cu)
        call matvec2(dd(1,1,i),u,dv)
        call matvec2(ee(1,1,i),u,ev)
        elocal=scalar2(b1(1,i),v)+scalar2(b2(1,i),u)+scalar2(cu,v)+
     &   scalar2(dv,u)+scalar2(ev,v)
      else 
        call matvec2(cc(1,1,i),u,cu)
        call matvec2(dd(1,1,i),v,dv)
        call matvec2(ee(1,1,i),v,ev)
        elocal=scalar2(b1(1,i),u)+scalar2(b2(1,i),v)+scalar2(cu,u)+
     &   scalar2(dv,v)+scalar2(ev,u)
      endif
      return
      end
c-------------------------------------------------------------------------
      double precision function ele(x,y,a)
      implicit none
      double precision x,y,a(2,2),si1,si2,u(2),v(2),av(2)
      double precision scalar2
      u(1)=-cos(x)
      u(2)= sin(x)
      v(1)=-cos(y)
      v(2)= sin(y)
      call matvec2(a,v,av)
      ele=scalar2(u,av) 
      return
      end