program checkstate
      implicit none
      integer*4 nstate,nc,i,j,it,iargc,nd,ix,ty,iy,iz
      real*8 w,g,excnm,wi,l,ev,cm,t,s,al(3,3),f2,f1,g2,m2,st,
     1gp(3,3),A(3,3,3),cum,dd,dm,cuma,cumg,p2,ff,gg
      real*8,allocatable:: d(:,:),e(:),c(:),ali(:,:,:),m(:,:),
     1q(:,:),qm(:,:,:),gpi(:,:,:),ai(:,:,:,:)
      integer*4,allocatable::ind(:)
      character*80 fn,s80
      if(iargc().ne.4)then
       write(6,601)
 601   format(' Usage: checkstate <filename>',
     1 ' <excnm> <Gammacm-1> <type>')
       stop
      endif
      call getarg(1,fn)
      call getarg(2,s80)
      read(s80,*)excnm
      call getarg(3,s80)
      read(s80,*)g
      g=g/219474.5d0
      g2=g**2
      call getarg(4,s80)
      read(s80,*)ty
      w=(1.0d7/excnm)/219474.5d0
      nc=nstate(fn)
      write(6,600)nc
600   format(i6,' states')
      allocate(e(nc),d(nc,3),c(nc),ind(nc),ali(nc,3,3),m(nc,3),q(nc,6),
     1qm(nc,3,3),gpi(nc,3,3),ai(nc,3,3,3))
      d=0.0d0
      e=0.0d0
      c=0.0d0
      ali=0.0d0

      open(8,file=fn)
1     read(8,80,end=88,err=88)s80
80    format(a80)
      if(s80(2:65).eq.'Ground to excited state transition electric dipol
     1e moments (Au):')then
       read(8,*)
       do 2 i=1,nc
       read(8,80)s80
2      read(s80(11:len(s80)),*)(d(i,j),j=1,3)
      endif
      if(s80(2:65).eq.'Ground to excited state transition magnetic dipol
     1e moments (Au):')then
       read(8,*)
       do 21 i=1,nc
       read(8,80)s80
21     read(s80(11:len(s80)),*)(m(i,j),j=1,3)
      endif
      if(s80(2:69).eq.'Ground to excited state transition velocity quadr
     1upole moments (Au):')then
c     XX          YY          ZZ          XY          XZ          YZ
       read(8,*)
       do 22 i=1,nc
       read(8,80)s80
22     read(s80(11:len(s80)),*)(q(i,j),j=1,6)
      endif
      if(s80(2:15).eq.'Excited State '.and.s80(16:16).ne.'s')then
       do 3 i=1,79
       if(s80(i:i)  .eq.':') read(s80(15:i-1),*)nd
3      if(s80(i:i+1).eq.'eV')read(s80(i-10:i-2),*)ev
       e(nd)=ev/27.211384205943d0
      endif
      goto 1
 88   close(8)

      m= 0.5d0*m
      do 23 i=1,nc
      do 231 ix=1,6
231   q(i,ix)=-q(i,ix)/e(i)
      qm(i,1,1)=(3.0d0*q(i,1)-q(i,1)-q(i,2)-q(i,3))/2.0d0
      qm(i,2,2)=(3.0d0*q(i,2)-q(i,1)-q(i,2)-q(i,3))/2.0d0
      qm(i,3,3)=(3.0d0*q(i,3)-q(i,1)-q(i,2)-q(i,3))/2.0d0
      qm(i,1,2)=q(i,4)*1.5d0
      qm(i,2,1)=q(i,4)*1.5d0
      qm(i,1,3)=q(i,5)*1.5d0
      qm(i,3,1)=q(i,5)*1.5d0
      qm(i,2,3)=q(i,6)*1.5d0
23    qm(i,3,2)=q(i,6)*1.5d0
 

      cm=0.0d0
      do 4 i=1,nc
      wi=e(i)
      m2=(wi-w)**2
      p2=(wi+w)**2
      dd=d(i,1)*d(i,1)+d(i,2)*d(i,2)+d(i,3)*d(i,3)
      dm=d(i,1)*m(i,1)+d(i,2)*m(i,2)+d(i,3)*m(i,3)
      ff=(wi-w)/(m2+g2)
      gg=(wi+w)/(p2+g2)
      f2=ff+gg
      f1=ff-gg
      do 41 ix=1,3
      do 41 iy=1,3
      ali(i,ix,iy)=d(i,ix)*d(i,iy)*f2
      gpi(i,ix,iy)=d(i,ix)*m(i,iy)*f1
      ai(i,ix,iy,1) =d(i,ix)*qm(i,iy,1)*f2
      ai(i,ix,iy,2) =d(i,ix)*qm(i,iy,2)*f2
41    ai(i,ix,iy,3) =d(i,ix)*qm(i,iy,3)*f2
      if(ty.eq.1)then
       c(i) =dabs(dd*f2)
      else
       if(ty.eq.2)then
        c(i)=dabs(dd*( (m2-g2)/(m2+g2)**2 +(p2-g2)/(p2+g2)**2) )
       else
        if(ty.eq.3)then
         c(i)=dabs(dd*dm/(m2+g2)**2)
        else
         write(6,*)ty,' unknown type'
         stop
        endif
       endif
      endif
 4    if(c( i).gt.cm )cm =c( i)
c     static limit:      
      gpi=gpi/w

      st=0.0d0
      do 51 i=1,nc
51    st=st+c(i)
      al=0.0d0
      gp=0.0d0
      a=0.0d0
      write(6,6022)
      s=0.0d0
      open(9,file='1.txt')
      do 5 i=1,nc
      ind(i)=i
      s=s+c(i)
      l=1.0d7/(e(i)*219474.5d0)
      dd=d(i,1)*d(i,1)+d(i,2)*d(i,2)+d(i,3)*d(i,3)
      dm=d(i,1)*m(i,1)+d(i,2)*m(i,2)+d(i,3)*m(i,3)
      do 82 ix=1,3
      do 82 iy=1,3
      a(ix,iy,1)=a(ix,iy,1)+ai(i,ix,iy,1)
      a(ix,iy,2)=a(ix,iy,2)+ai(i,ix,iy,2)
      a(ix,iy,3)=a(ix,iy,3)+ai(i,ix,iy,3)
      gp(ix,iy)=gp(ix,iy)+gpi(i,ix,iy)
82    al(ix,iy)=al(ix,iy)+ali(i,ix,iy)
      cuma=(al(1,1)+al(2,2)+al(3,3))/3.0d0
      cumg=(gp(1,1)+gp(2,2)+gp(3,3))/3.0d0
      cum=0.0d0
      do 101 ix=1,3
      do 101 iy=1,3
      do 101 iz=1,3
101   cum=cum+2.0d0*a(ix,ix,iy)*(a(iy,iz,iz)+2.0d0*a(iz,iy,iz))
     1             +a(ix,iy,iy)*(a(ix,iz,iz)+2.0d0*a(iz,ix,iz))
     1       +2.0d0*a(ix,iy,iz)*(a(ix,iy,iz)+a(iy,ix,iz)+a(iz,ix,iy))
      cum=cum/105.0d0
      write(9,602)i,l,dd,dm,c(i)/cm*100.0d0,s/st*100.0d0,cuma,cumg,cum
  5   write(6,602)i,l,dd,dm,c(i)/cm*100.0d0,s/st*100.0d0,cuma,cumg,cum
602   format(i6,f9.1,2f12.6,2f6.1,2f12.3,f12.1)
6022  format(' state     l/nm           D           R   C/%  sC/%',
     1       '       alpha           G           A')
      close(9)

      do 6 i=1,nc
      do 6 j=i+1,nc
      if(c(j).gt.c(i))then
       it=ind(i)
       ind(i)=ind(j)
       ind(j)=it
       t=c(i)
       c(i)=c(j)
       c(j)=t
       t=e(i)
       e(i)=e(j)
       e(j)=t
       do 8 ix=1,3
       t=d(i,ix)
       d(i,ix)=d(j,ix)
       d(j,ix)=t
       t=m(i,ix)
       m(i,ix)=m(j,ix)
       m(j,ix)=t
       do 8 iy=1,3
       t=ali(i,ix,iy)
       ali(i,ix,iy)=ali(j,ix,iy)
       ali(j,ix,iy)=t
       t=gpi(i,ix,iy)
       gpi(i,ix,iy)=gpi(j,ix,iy)
8      gpi(j,ix,iy)=t
      endif
6     continue

      write(6,603)

      s=0.0d0
      al=0.0d0
      gp=0.0d0
      open(9,file='2.txt')
      do 7 i=1,nc
c     original index:
      j=ind(i)
      dd=d(i,1)*d(i,1)+d(i,2)*d(i,2)+d(i,3)*d(i,3)
      dm=d(i,1)*m(i,1)+d(i,2)*m(i,2)+d(i,3)*m(i,3)
      l=1.0d7/(e(i)*219474.5d0)
      s=s+c(i)
      do 81 ix=1,3
      do 81 iy=1,3
      gp(ix,iy)=gp(ix,iy)+gpi(i,ix,iy)
81    al(ix,iy)=al(ix,iy)+ali(i,ix,iy)
      cuma=(al(1,1)+al(2,2)+al(3,3))/3.0d0
      cumg=(gp(1,1)+gp(2,2)+gp(3,3))/3.0d0
      write(9,604)i,j,l,dd,dm,cuma,cumg,c(i)/cm*100.0d0 ,s/st*100.0d0
  7   write(6,604)i,j,l,dd,dm,cuma,cumg,c(i)/cm*100.0d0 ,s/st*100.0d0
604   format(i6,' / ',i5,f9.1,2f12.6,2f12.3,2f7.1)
603   format('Most constributing:',/,
     1' State / Orig.     l/nm           d           r',
     1'       alpha           G    p/%   pr/% ')
      close(9)

      end


      function nstate(filename)
      implicit none
      character*(*) filename
      integer*4 nstate,r
      character*(80) s80
      r=0
      open(8,file=filename)
1     read(8,80,end=88,err=88)s80
80    format(a80)
      if(s80(2:14).eq.'Excited State')r=r+1
      goto 1
88    close(8)
      nstate=r
      return
      end