program guvcd3
      character*80 fo,s80,output,ext*10
      character*1 ok
      parameter (nm0=1000)
      real*8 d(nm0,3),r(nm0,3),e(nm0),ev(nm0),f(nm0),r0(nm0),ecm(nm0)
c
c     30. 1. 2006 zmenen zpusob vypoctu dipolove a rotacni sily
c     a zavedeno vytvareni kontrolniho souboru 
      write(6,*)' Reads UVCD from G03 output'
      write(6,*)' and writes it in a TAB file.'
      write(6,*)
      if(iargc().eq.0)then
       write(6,*)' Name of the output:'
       read(5,'(a)')fo
       write(6,*)'Length or velocity formalism or both (L/V/B) ?'
       read(5,'(a)')ok
      else
       call getarg(1,fo)
       call getarg(2,ok)
      endif
      if(ok.eq.'l')ok='L'
      if(ok.eq.'b')ok='B'

14    open(2,file=fo)
      n=0
      line=0
c
1     read(2,2000,end=2,err=2)s80
c     write(6,2000)s80
      line=line+1
2000  format(a80)

      if((s80(2:44).eq.'Ground to excited state transition electric'
     1.and.ok.eq.'L').or.
     2   (s80(2:51).eq.
     3'Ground to excited state transition velocity dipole'
     3.and.ok.ne.'L'))then
       write(6,2000)s80
       read(2,*)
       line=line+1
3      read(2,*,end=21,err=21)nd,ax,ay,az,osc
       line=line+1
       n=n+1
       if(n.gt.nm0)then
        write(6,*)' Too many transitions ',nm0
        close(2)
        stop
       endif
c      14. 12. 2005 sem dodana korekce znamenek v delkovem formalismu
       if(ok.eq.'L')then
         d(n,1)=-ax
         d(n,2)=-ay
         d(n,3)=-az
        else
         d(n,1)=ax
         d(n,2)=ay
         d(n,3)=az
       endif
       f(n)=osc
       goto 3
21     backspace 2
       write(6,*)n,' electric transitions'
      endif
c
      if(s80(2:44).eq.'Ground to excited state transition magnetic')
     1then
       write(6,2000)s80
       read(2,*)
       line=line+1
       do 5 i=1,n
       line=line+1
       read(2,*)nd,(r(i,j),j=1,3)
c      korekce Gaussianu - ten nedeli magneticke momenty cislem 2
       do 12 j=1,3
12     r(i,j)=r(i,j)/2
5      continue
       write(6,*)n,' magnetic transitions'
       ie=0
      endif
c      
      if(s80(53:63).eq.'R(velocity)'.and.ok.ne.'L') then
c      Gaussian 98 does not write R(length) nor R(velocity)
       do 8 i=1,n
       line=line+1
       read(2,2000)s80
       read(s80,*,err=99)nd,b1,b2,b3,r0(i)
       goto 991
99     write(6,2000)s80
       write(6,*)i
       stop
991    line=line+1
       read(2,2000)s80
       if(s80(2:6).eq.'Total')then
        do 81 k=1,4
        line=line+1
81      read(2,2000)s80
       else
        backspace 2
       endif
       read(2,2000)s80
       if(s80(2:3).eq.'R(')then
        do 82 k=1,4
        line=line+1
82      read(2,2000)s80
       else
        backspace 2
       endif
8      continue
      endif

      if(s80(54:62).eq.'R(length)'.and.ok.eq.'L') then
       write(6,2000)s80
       line=line+1
       do 9 i=1,n
       line=line+1
9      read(2,*,end=2,err=2)nd,b1,b2,b3,r0(i)
      endif
c
      if(s80(2:15).eq.'Excited State '.and.s80(16:16).ne.'s')then
       write(6,2000)s80
       read(s80(15:18),*)nd
       do i=1,79
        if(s80(i:i+1).eq.'eV')read(s80(i-10:i-2),*)ev(nd)
        if(s80(i:i+1).eq.'nm')read(s80(i- 9:i-2),*)e(nd)
       enddo
       close(8)
       ie=ie+1
       if(nd.eq.n)then
        write(6,*)ie,' energies'
c
        if(ok.eq.'L')then
         ext='.l.TAB'
        else
         ext='.v.TAB'
        endif
c
        call titles(fo,ext,output)
        write(6,*)output
	open(3,file=output)
c
	if(ok.eq.'L')then
         write(3,3000)
3000    format('          UVCD spectrum from G03 in length formalism',
     1   /,'   n wavelength (nm) dipole strength (D^2)',
     1   ' rotatory strength (cgs/10**-36)',/,80(1h-))
3200  format('     check of UVCD spectrum from G03 in length formalism',
     1   /,'   n wavelength (nm) dipole strength (D^2)',
     1   ' rotatory strength (cgs/10**-36)',/,80(1h-))
	else
         write(3,3100)
3100    format('          UVCD spectrum from G03 in velocity formalism',
     1   /,'   n wavelength (nm) dipole strength (D^2)',
     1   ' rotatory strength (cgs/10**-36)',/,80(1h-))
3300  format('    check of UVCD spectrum from G03 in velocity formalism',
     1   /,'   n wavelength (nm) dipole strength (D^2)',
     1   ' rotatory strength (cgs/10**-36)',/,80(1h-))
	endif
c
        do 6 i=1,n
	ds=2.12689677295232d-30*f(i)/ev(i)/8065.54476345045d0/1e-36
c       korekce momentu ve velocity formalismu - deleni frekvenci
        if(ok.ne.'L') then
	 do 11 j=1,3
	  d(i,j)=d(i,j)/ev(i)*27.211384205943d0
11	 continue
        endif
        ds1=(d(i,1)*d(i,1)+d(i,2)*d(i,2)+d(i,3)*d(i,3))*2.541765*2.541765
	rs=r0(i)*1e-4
        rs1=(d(i,1)*r(i,1)+d(i,2)*r(i,2)+d(i,3)*r(i,3))
     1 	*2*2.541765*9.2740154d-21*1e-18/1e-36
c
c       angle between electric and magnetic moment:
        de=dsqrt(d(i,1)*d(i,1)+d(i,2)*d(i,2)+d(i,3)*d(i,3))
        dm=dsqrt(r(i,1)*r(i,1)+r(i,2)*r(i,2)+r(i,3)*r(i,3))
        ds=     (d(i,1)*r(i,1)+d(i,2)*r(i,2)+d(i,3)*r(i,3))
        angle=0.0d0
        if(de*dm.gt.0.0d0)then
         angle=acos(ds/de/dm)*180.0d0/4.0d0/atan(1.0d0)
        endif

6       write(3,3001)i,e(i),ds1,rs,angle
c       6. 2. 2006 zmeneno psani vypoctene dipolove sily - do souboru TAB psan
c       vypocet z momentu, do kontrolniho souboru vypocet z oscilatorove sily
3001    format(i4,f12.2,2f20.10,f8.1)
        write(3,3002)
3002    format(80(1h-))
3003    format
     1  ('   n    energy (eV) wavenumber (cm**-1)  oscillator strength')
3004    format(i4,f12.4,f20.6,f17.4)
        do 13 i=1,n
13      ecm(i)=ev(i)*8065.54476345045d0
        close(3)
        write(6,*)'Output TAB file written.'
        goto 2
       endif
      endif
      goto 1
c
2     close(2)
      write(6,*)
      write(6,*)line,' lines read'
      if(ok.eq.'B')then
       ok='L'
       goto 14
      endif
      end
c      
      subroutine titles(ti,ext,out)
c     this subroutine changes the names of output files      
      character*(*) ti,ext
      character*80 out
c      
      do 1 i=1,80
       if (ti(i:i).eq.'.'.or.ti(i:i).eq.' ') then
         k=i-1
	 goto 2
       endif
1     continue
c
2     out=ti(1:k)//ext
      return
      end