c      module m_cgwdos 
c      contains
      subroutine cgwdos(volwgt,ec ,emin,emax,dos,nr)
#define WRITEDOS
#undef WRITEDOS
c
c from slinz
c
C- Adds to PDOS for one tetrahedron
C ----------------------------------------------------------------
Ci Inputs
Ci   volwgt, weight on tetrahedron; ec energies at corners of tethdn.;
Ci   emin, emax, energy window; nr, number of bins + 1
ci   ec: energy dependent E_nk(E) including imaginary part from emin to emax
Co Outputs
Co   dos(k),  density of states in kth bin from tethdn.
Cr Remarks
c
c Lambin and Vigneron , PRB 29 , 3430. 
c eq.7 
c  r_i(z) = (z-Ei)**2/ /Pi_{k!=i} (E_k-E_i ) + 
c           Sum_{j!=i} (z-E_j)**3 / Pi_{k!=j} (E_k-E_j) * Ln ( (z-E_j)(z-E_i) )/ (E_i-E_j)
c  
c
Cr
C ----------------------------------------------------------------
C     implicit none
C Passed parameters
      integer nr
      double precision dos(nr),volwgt,emin,emax
      complex(8) :: ec(4,nr)
C Local parameters

      integer ie, i,k,j
      complex(8) :: r(4)
      real(8) :: de,pi
      complex(8) ::z , etmp(4)
      parameter(pi=3.1415926535897932384626433832795d0)

      integer:: iec(4),ithesame(2), indexk(4)
      complex(8),allocatable :: ectmp(:,:)

c initialize 
c      dos = 0.0d0

c reorder
      iec=(/1,2,3,4/)
      do i=1,3
        do j=i+1,4
          if ( real(ec(iec(i),1)) > real(ec(iec(j),1)) ) then
            k=iec(j)
            iec(j)=iec(i)
            iec(i)=k
          endif
        enddo
      enddo
c check really reordered or not
      do i=1,3
         if ( real(ec(iec(i),1)) > real(ec(iec(i+1),1)) ) then
            write(*,*) 'order is wrong ', iec, ec(iec(:),1) 
            stop ' in cgwdos'
         endif
      enddo
c which case?
      icase=0  ! all different
      do i=1,3
        if ( ec(iec(i),1) .eq. ec(iec(i+1),1) ) then
          icase=icase+1
          ithesame(1)=i
          ithesame(2)=i+1
        endif
      enddo
      if (icase>0) then
         allocate( ectmp(4,nr) )
         do i=1,4
           ectmp(i,:) = ec(iec(i),:) 
         enddo
         ec = ectmp 
         deallocate(ectmp) 
      endif


c calc r_j(E) and sum them
      de = (emax-emin)/(nr-1)


      if (icase.eq.0) goto 100 
      if (icase.eq.1) goto 200
      if (icase.eq.2)  then
         if ( ec(1,1).eq.ec(2,1) .and. ec(2,1).eq.ec(3,1) .and.  ec(3,1).ne.ec(4,1) ) goto 300
         if ( ec(1,1).eq.ec(2,1) .and. ec(2,1).ne.ec(3,1) .and.  ec(3,1).eq.ec(4,1) ) goto 400
         if ( ec(1,1).ne.ec(2,1) .and. ec(2,1).eq.ec(3,1) .and.  ec(3,1).eq.ec(4,1) ) goto 500
      endif
      if (icase.eq.3) goto 600


c all different 
 100  continue
      do ie=1,nr
         z= emin+ de*(ie-1)

         call eq7(z,ec(1,ie),r)

         do i=1,4
            dos(ie)= dos(ie)+  abs( imag( r(i) ) )  * volwgt /pi 
         enddo
#ifdef WRITEDOS
        write(*,*)ie, r(1:4)
#endif

      enddo 
      return

c calc r_j(E) and sum them
c     use eq. (A1)
 200  continue
      l = ithesame(1) ; m =ithesame(2) ;ie = 0
      do i=1,4
        if (i.ne.l .and. i.ne.m) then
        do j=1,4
          if ( j.ne.i .and. j.ne.l .and. j.ne.m ) then
            ie=ie+1
            indexk(ie)=j
          endif
        enddo
        endif
      enddo

c     for l,m (E_l=E_m) and  indexk(1), indexk(2)


      do ie=1,nr
         z= emin+ de*(ie-1)

         etmp=(/ec(indexk(1),ie),ec(indexk(2),ie), ec(l,ie), ec(m,ie) /)
         call eqA1and2(z,etmp,r)

         do i=1,4
            dos(ie)= dos(ie)+  abs( imag( r(i) ) )  * volwgt /pi
         enddo
#ifdef WRITEDOS
        write(*,*)ie, r(1:4)
#endif


      enddo
      return

 300  continue
c                      A3
      do ie=1,nr
         z= emin+ de*(ie-1)

         call eqA3(z,ec(1,ie),r)
     
         do i=1,4
            dos(ie)= dos(ie)+  abs( imag( r(i) ) )  * volwgt /pi
         enddo
#ifdef WRITEDOS
        write(*,*) ie,r(1:4)
#endif


      enddo
      return

 400  continue
C                       A4
      do ie=1,nr
         z= emin+ de*(ie-1)

         call eqA4(z,ec(1,ie),r)

         do i=1,4
            dos(ie)= dos(ie)+  abs( imag( r(i) ) )  * volwgt /pi
         enddo
#ifdef WRITEDOS
        write(*,*) ie,r(1:4)
#endif


      enddo
      return

 500  continue
C                      A5
      do ie=1,nr
         z= emin+ de*(ie-1)

         call eqA5(z,ec(1,ie),r)

         do i=1,4
            dos(ie)= dos(ie)+  abs( imag( r(i) ) )  * volwgt /pi
         enddo
#ifdef WRITEDOS
        write(*,*) ie,r(1:4)
#endif


      enddo
      return

 600  continue
C                    A6
      do ie=1,nr
         z= emin+ de*(ie-1)

         call eqA6(z,ec(1,ie),r)

         do i=1,4
            dos(ie)= dos(ie)+  abs( imag( r(i) ) )  * volwgt /pi
         enddo
#ifdef WRITEDOS
        write(*,*) ie,r(1:4)
#endif


      enddo
      return

      end subroutine cgwdos


c-----------------------------------------
      subroutine eq7(z,e,r)
      implicit none
      complex(8) z,e(4),r(4)

      integer:: i, j,k
      complex(8) :: eme,t
      do  i=1,4
        r(i)= (z-e(i))**2 
        eme=1.0d0
        do k=1,4
          if (k.ne.i) then
             eme=eme*(e(k)-e(i))
          endif
        enddo
        r(i)=r(i)/eme

        do j=1,4
          if (j.ne.i) then
           t= (z-e(j))**3 
           t = t* log( (z-e(j))/(z-e(i)) )
           t = t/(e(i)-e(j))
           eme=1.0d0
           do k=1,4
             if (k.ne.j) then
                 eme=eme*(e(k)-e(j))
             endif
           enddo
           r(i) = r(i) + t/eme
          endif
        enddo
      enddo


      end subroutine eq7

c---------------------------------------------------------------------
      subroutine eqA1and2(z,e,r)
      implicit none
      complex(8) z,e(4),r(4)

      complex(8) :: ei,ek, em, ej

      ei=e(1); ek=e(2);! el=e(3); 
      em=e(4)

c         A1
      r(1) =  (z-ei)**2/( (ei-em)**2*(ek-ei) ) * 
     .  (  1.0d0 + ( 2.0d0*( z-em )/(ei-em) + (z-ek)/(ei-ek) )*log(z-ei) )
     .  + (z-em)**2/( (em-ei)**2*(ek-em) ) *
     .  ( 1.0d0 + ( 2.0d0*(z-ei)/(em-ei)+ (z-ek)/( em-ek ) )*log(z-em) ) +
     .   (z-ek)**3/( (ek-ei)**2*(ek-em)**2 )* log(z-ek) 
      r(2) = r(1)

c        A2
      ej=e(1); ek=e(2); !el=e(3);
      em=e(4)

      r(3) =   (z-ej)**3/( (ej-ek)*(ej-em)**3 ) * log(z-ej) + 
     .  (z-ek)**3/( (ek-ej)*(ek-em)**3 ) * log(z-ek)  +
     .   (z-em)/( (em-ej)*(em-ek) )* 
     .    (
     .        0.5 + (z-ej)/(em-ej) + (z-ek)/(em-ek) +
     .       ( 
     .         ((z-ej)/(em-ej))**2+ ((z-ek)/(em-ek))**2 + (z-ej)*(z-ek)/(em-ej)/(em-ek)  
     .       ) * log(z-em)
     .    )
      r(4)=r(3)

      end subroutine eqA1and2

c---------------------------------------------------------------------
      subroutine eqA3(z,e,r)
      implicit none
      complex(8) z,e(4),r(4)

      ! complex(8) :: e1,e2,
      complex(8) ::e3,e4

      ! e1=e(1); e2=e(2); 
      e3=e(3); e4=e(4)

      
         r(1)= (z-e4)**3/(e4-e3)**4*log( (z-e4)/(z-e3) )+
     .     ( 6.0d0*(z-e4)**2-3.0*(e4-e3)*(z-e4)+ 2*(e4-e3)**2 )
     .    / ( 6.0d0*(e4-e3)**3 )
         r(2)=r(1); r(3)=r(1); 

         r(4)= 3.0d0*(z-e4)**2*(z-e3)/(e4-e3)**4 * log( (z-e3)/(z-e4) )
     .     - 1.5d0*(z-e3)*
     .       (
     .          2.0d0*(z-e4)-(e4-e3)
     .       ) / (e4-e3)**3
     .     - 1.0d0/(e4-e3)        


      end subroutine eqA3
c---------------------------------------------------------------------

      subroutine eqA4(z,e,r)
      implicit none
      complex(8) z,e(4),r(4)

      complex(8) :: e2,e3! ,e4,e1
      
      !e1=e(1); 
      e2=e(2); e3=e(3);
      !e4=e(4)
      r(1) = 3.0d0*(z-e3)**2*(z-e2)/(e3-e2)**4*
     .   log( (z-e2)/(z-e3) ) -
     .   1.5d0*(z-e2)*( 2.0d0*(z-e3)- (e3-e2) )/(e3-e2)**3 -
     .   1.0d0/(e3-e2) 
      r(2)=r(1)

      r(3) = 3.0d0*(z-e2)**2*(z-e3)/(e3-e2)**4*
     .   log( (z-e3)/(z-e2) ) +
     .   1.5d0*(z-e3)*( 2.0d0*(z-e2)+ (e3-e2) )/(e3-e2)**3 +
     .   1.0d0/(e3-e2)
      r(4)=r(3)

      end subroutine eqA4
c---------------------------------------------------------------------
      subroutine eqA5(z,e,r)
      implicit none
      complex(8) z,e(4),r(4)

      complex(8) :: e1,e2!,e3,e4

      e1=e(1); e2=e(2)!; e3=e(3); e4=e(4)

      r(1) = 3.0d0*(z-e1)**2*(z-e2)/ (e2-e1)**4* log ( (z-e2)/(z-e1) )
     .     + 1.5d0*(z-e2)*( 2.0d0*(z-e1)+ e2-e1 )/(e2-e1)**3 +
     .     1.0d0/(e2-e1)

      r(2) = (z-e1)**3/(e2-e1)**4 * log( (z-e1)/(z-e2) ) -
     .     ( 6.0d0*(z-e1)**2 + 3.0d0*(z-e1)*(e2-e1) + 2.0d0*(e2-e1)**2 )
     .     / ( 6.0d0*(e2-e1)**3 ) 
      r(3)=r(2)
      r(4)=r(2)

      end subroutine eqA5

c---------------------------------------------------------------------
      subroutine eqA6(z,e,r)
      implicit none
      complex(8) z,e(4),r(4)

      r(1) = 0.25d0/( z- e(1) )
      r(2)=r(1)
      r(3)=r(1)
      r(4)=r(1)

      end  subroutine eqA6
c      end module 
 
