      program gmul
      implicit none
      integer          MAXKPTS
      parameter(MAXKPTS=20)
      integer          level,spin,atom,nkpts,i,k
      double precision w(MAXKPTS),kpt(MAXKPTS,3),pop(10),avepop(10)
      double precision stot,ptot,dtot,eta,wtot

      character*10     dump_str
      character        dump_char
      integer          dump_int
      integer          iargc,nargs

      character*80     resfile,levelstr,spinstr,atomstr

c     check arguments
      nargs=iargc()
      if (nargs.eq.3) then
         write(*,*) "I am expecting a spin averaged res.file..."
         call getarg(1,resfile)
         call getarg(2,levelstr)
         call getarg(3,atomstr)
         spin=1
      else if (nargs.eq.4) then
         write(*,*) "I am expecting a spin polarised res.file..."
         call getarg(1,resfile)
         call getarg(2,levelstr)
         call getarg(3,spinstr)
         call getarg(4,atomstr)
         read(spinstr,*)  spin
      else
         write(*,*) "syntax:"
         write(*,*) "gmul res.file level [spin] atom"
         write(*,*) " "
         write(*,*) "* level values must be present in the ",
     &              "analysis res.file;"
         write(*,*) "* spin argument is optional, and may take ",
     &              "the following values:"
         write(*,*) "   1 - spin up"
         write(*,*) "   2 - spin down"
         stop
      endif
      read(levelstr,*) level
      read(atomstr,*)  atom
      open(unit=10,file=resfile)
      open(unit=99,file="dat.bandst")
      
      do k=1,10
        avepop(k)=0.0d0
      end do
      stot=0.0d0
      ptot=0.0d0
      dtot=0.0d0
      wtot=0.0d0
      read(99,*) nkpts
      read(99,*) (w(i),(kpt(i,k),k=1,3),i=1,nkpts)
      do i=1,nkpts
        wtot=wtot+w(i)
      end do
      do i=1,nkpts
        w(i)=w(i)/wtot
      end do
      write(*,*) "n kpoints: ",nkpts
      write(*,*) "weights  : ",(w(i),i=1,nkpts)
      do i=1,nkpts
        call getpos(level,spin,atom,pop)
        write(*,100) i,(pop(k),k=1,10)
        do k=1,10
          avepop(k)=avepop(k)+w(i)*pop(k)
        end do
      end do
      stot=avepop(1)
      do i=1,3
        ptot=ptot+avepop(i+1)
      end do
      do i=1,6
        dtot=dtot+avepop(i+4)
      end do
      eta=stot+ptot+dtot
      write(*,50) atom,(avepop(i),i=1,10)
      write(*,*) "l       pop(l)      eta(l)"
      write(*,110)
      write(*,70) stot/eta,stot
      write(*,80) ptot/eta,ptot
      write(*,90) dtot/eta,dtot
      write(*,110)
      write(*,60) stot/eta+ptot/eta+dtot/eta,stot+ptot+dtot
50    format("atom :",i3,"   ",10f9.5)
60    format("tot  :",f9.5,"    ",f9.5)
70    format("s    :",f9.5,"    ",f9.5)
80    format("p    :",f9.5,"    ",f9.5)
90    format("d    :",f9.5,"    ",f9.5)
100   format("kpt  :",i3,"   ",10f9.5)
110   format("----------------------------")
      close(unit=10)
      close(unit=99)
      stop
      end


      subroutine getpos(level,spin,atom,pop)
      double precision pop(10)
      integer          level,spin,atom,i,tmp1,tmp2
      character*80     s1,s2,s3,ss,ls,ks
      logical          gotit
      gotit=.false.
      read(10,30) s1,ks,s2,ss,s3,ls
      do while (.not.gotit)
        read(10,30) s1,ks,s2,ss,s3,ls
        if (s1.eq."WAVEFUNCTION ANALYSIS: kpt") then
          read(ls,*) tmp1
          read(ss,*) tmp2
          if ((level.eq.tmp1).and.(spin.eq.tmp2)) then
            gotit=.true.
          end if
        end if
      end do    
      do i=1,3+atom
        read(10,*) s1
      end do
      read(10,40) s1,(pop(i),i=1,10)
30    format(a26,a4,a5,a2,a6,a4)
40    format(a12,10f9.5)
      return
      end


