program tabprnf
      implicit none
      integer*4 np0,ic,npoint,i,i1,iargc,nm,column,length
      PARAMETER (np0=2000000)
      real*8 wsmax,wsmin,s(np0),t(np0),x,dw,dd,hw,xx,const,wo,w(np0),
     1sf,temp,cm,cmo,sfc,wm,excnm,exccm,coefinv(5),refrindex,roainv(5),
     2 dscs_ww,dscs_pw,hck,tc(np0)
      character*80 fn,s80
      character*1 pt
      character*30 sexpmode,sunit 
      logical LG,lnm,lwe
c     lwe - weather to use constant bandwidth in cm-1 scale
c
      write(6,6009)
6009  format(' Usage: tabprnf name type npoints min max G T width ri',/,
     1' Arguments: 1 name - name of the .tab file ',/,
     1'            2 type - type of spectrum (column # in name)',/,
     1'            3 npoints - number of spectral points',/,
     1'            4 min - wmin',/,
     1'            5 max -wmax',/,
     1'            6 G - Gaussian peaks instead of Lorentzians',/,
     1'            7 T - temperature',/,
     1'            8 width - full width at half height',/,
     1'              negative number: constant bandwidth in reciprocal'
     1//' scale (cm-1 or nm)',/,
     1'            9 ri - refractive index',/)
      write(6,6008)
6008  format(
     1' Type options:',/,
     1'    1 abs,2 vcd, 3 ram-old, 8 roa-old, 20 mcd, 30 ord, 40 nmr',/,
     1'    negative number above: cm-1 converted to nm or nm to cm-1',/,
     2'    -RAMROAEXP/UNITTYPE:',/,
     2'        RAMROAEXP: RAM_SCP_180, RAM_DCPI_180, RAM_DCPII_180, ',
     2'ROA_SCP_180',/,
     2'        UNITTYPE: KA4amu, WWcm6, WWcm2, pWcm2J, ppcm2',/)
c     3'    -INV/d1,d2,d3,d4,d5/common_denominator/UNITTYPE',/,
c     3'        d1,d2,d3,d4,d5: tensor invariant coefficients',/,
c     3'        common_denominator: com. den. for ten.inv.coef',/,
c     3'        UNITTYPE: as above',/)
c
      lnm=.false.
      if(iargc().lt.1)then
       write(6,*)'.TAB filename:'
       read(5,500)fn
500    format(a80)
      else
       call getarg(1,fn)
       write(6,*)'generating spectra from'
       write(6,*)fn
      endif
c
      if(iargc().gt.1)then
        call getarg(2,s80)
      else
        write(6,*)'Type spectra (see above for options):'
        read(*,'(a)')s80
      endif
      call upper_case(s80)
      write(*,*)'debug>',s80(1:length(s80))
      if(s80(1:2).eq.'-R')then ! option is string (Raman definition)
        call readramopt(s80,coefinv,sexpmode,sunit)
        ic=100
        write(*,*)'Raman experiment: ',sexpmode
        write(*,*)' ic: ',ic
        write(*,*)' invariant coefficients: ',(coefinv(i1),i1=1,5)
        write(*,*)' output unit: ',sunit
      else ! option is number
        read(s80,*)ic
        if(ic.lt.0)then
          ic=-ic
          lnm=.true.
          write(6,*)'cm-1 converted to nm !!!'
         endif
         write(6,*)ic,' column after frequency'
         if(ic.eq.1 )write(6,*)'absorption'
         if(ic.eq.2 )write(6,*)'circular dichroism'
         if(ic.eq.3 )write(6,*)'Raman scattering'
         if(ic.eq.8 )write(6,*)'Raman optical activity'
         if(ic.eq.20)write(6,*)'magnetic circular dichroism'
         if(ic.eq.30)write(6,*)'ORD from Rs in the second column'
         if(ic.eq.40)write(6,*)'NMR from shifts in the second column'
       endif
c
      if(iargc().gt.2)then
       call getarg(3,s80)
       read(s80,*)npoint
       write(6,*)npoint,' points'
      else
       write(6,*)'number of points:'
       read(5,*)npoint
      endif
      if(npoint.gt.np0)call report('too many spectral points')
      if(iargc().gt.3)then
       call getarg(4,s80)
       read(s80,*)wsmin
       call getarg(5,s80)
       read(s80,*)wsmax
       write(6,6000)wsmin,wsmax
6000   format(' frequency interval: ',f12.2,' - ',f12.2)
      else
       write(6,*)'Wmin:'
       read(5,*)wsmin
       write(6,*)'Wmax:'
       read(5,*)wsmax
      endif
      if(iargc().gt.5)then
       call getarg(6,pt)
       if(pt.eq.'g'.or.pt.eq.'G')then
        write(6,*)'Gaussian peaks'
        LG=.true.
       else
        write(6,*)'Lorentzian peaks'
        LG=.false.
       endif
      else
       write(6,*)'peak type (G/L):'
       read(5,501)pt
501    format(a1)
       if(pt.eq.'G'.or.pt.eq.'g')then
        LG=.true.
       else
        LG=.false.
       endif
      endif
c
      if(iargc().gt.6)then
       call getarg(7,s80)
       read(s80,*)temp
       write(6,6001)temp
6001   format(' Temperature: ',f5.1,' K')
      else
       write(6,*)'temperature (K):'
       read(5,*)temp
      endif
c
      if(iargc().gt.7)then
       call getarg(8,s80)
       read(s80,*)hw
      else
       write(6,*)'peak width (FWHH):'
       read(5,*)hw
      endif
c
      if(hw.lt.0)then
       lwe=.true.
       hw=-hw
       write(6,502)hw
c502    format(' FWHH ',f10.3,' cm-1')
502    format(' FWHH ',f10.3,' in reciprocal unit')
      else 
       lwe=.false.
       write(6,503)hw
c503    format(' FWHH ',f10.3,' nm')
503    format(' FWHH ',f10.3,' constant in same unit')
      endif
c
      if(iargc().gt.8)then
       call getarg(9,s80)
       read(s80,*)refrindex
       write(6,6002)refrindex
6002   format(' Refractive index: ',f5.1)
      else
       write(6,*)'Refractive index:'
       read(5,*)refrindex
      endif
      refrindex=((refrindex**2+2.0d0)/3.0d0)**4
      write(*,*)' Lorentz factor^2 = ',refrindex
cc
      column=ic
      if(ic.ge.10)column=ic/10
      if(ic.eq.30)column=2
      if(ic.eq.40)column=1
      if(ic.eq.13)column=3
      if(ic.eq.18)column=8
c
      open(9,file=fn)
      nm=0
      if(ic.ne.100)then ! simple read from one column
        do 1 i=1,3
1       read(9,*)
2       read(9,*,end=99,err=99)wo,wo,(x,i=1,column)
        nm=nm+1
        if(nm.gt.np0)call report('Too many transitions')
        if(lnm)wo=1.0d7/wo
        w(nm)=wo
        t(nm)=x
        goto 2
      else ! ic.eq.100: multi columnn read (roa invariants)
        read(9,'(a)')s80
        read(s80(28:length(s80)),*)excnm
        read(9,*)
        read(9,*)
        exccm=1.d7/excnm
22      read(9,*,end=99,err=99)wo,wo,(roainv(i),i=1,5)
        nm=nm+1
        if(nm.gt.np0)call report('Too many transitions')
        if(lnm)wo=1.0d7/wo
        w(nm)=wo
        t(nm)=1.d0*(roainv(1)*coefinv(1)+roainv(2)*coefinv(2))+
     1        1.d0*(roainv(3)*coefinv(3)+roainv(4)*coefinv(4)+
     2              roainv(5)*coefinv(5))
c       Important:note that ROA invariants(3-5) are multiplied by 10**4
c       either Raman or ROA invariants can be used in same expression
c       defined in subroutine readramopt()
        goto 22
      endif
99    close(9)
      write(6,*)nm, ' modes'
c
      do 3 i=1,npoint
3     s(i)=0.0d0
      dw=(wsmax-wsmin)/dble(npoint-1)
      if(LG)then
       dd=hw/2.0d0/dsqrt(log(2.0d0))
      else
       dd=hw/2.0d0
      endif
c
      do 4 i=1,npoint
      xx=wsmin+dw*dble(i-1)
      do 4 i1=1,nm
      wo=w(i1)
      const=t(i1)
      if(lwe)then
c      lambda in cm-1 (generally reciprocal unit):
c      'cm' is not used for output S.PRN, only 'sfc'
c!!!   normalization of bands is not done properly (wrong unit)!!!!
       cm=1.0d7/xx
       cmo=1.0d7/wo
       sfc=sf(dd,LG,cm,cmo,ic)
       wm=1.0d7/wo
      else
       sfc=sf(dd,LG,xx,wo,ic)
       wm=wo
      endif
c!!!  Following lines for ic=1,2,20,3,8,13,18,30 are not optimal,
c!!!  since they are recalculated every 'i' step
c!!!  it is better to calculate them only once and save them into 
c!!!  field tc() as for ic=100
c     abs:
      if(ic.eq.1)const=108.7d0*t(i1)*wm
c     vcd:
      if(ic.eq.2)const=435.0d0*t(i1)*wm
c     mcd - B-terms:
      if(ic.eq.20)const=-5.98442d-3*t(i1)*wm
c     Raman:
      if(ic.eq.3)const=t(i1)/wo/(1.0d0-exp(-(1.44d0*wo/temp)))
c     roa:
      if(ic.eq.8)const=t(i1)/wo/(1.0d0-exp(-(1.44d0*wo/temp)))
c     Raman - absolute:
      if(ic.eq.13)const=t(i1)
c     roa - absolute:
      if(ic.eq.18)const=t(i1)
c     ord:
      if(ic.eq.30)const=t(i1)*8.35d18
c     nmr:
      if(ic.eq.40)const=t(i1)
      if(ic.eq.100)then
        if(i.eq.1)then
          dscs_ww = 2.6273379522352304d-44 ! cm5.A-4.amu
          dscs_pw = 1.3226326663297376d-21 ! cm4.J-1.A-4.amu
          hck = 1.4387772526779683d0 ! cm.K
         if(sunit(1:6).eq.'KA4AMU')then ! K*A**4/amu unit - Gaussian16
          tc(i1)=t(i1)
         elseif(sunit(1:5).eq.'WWCM6')then ! 
          tc(i1)=t(i1)*dscs_ww/wo/(1.d0-exp(-(hck*wo/temp)))
         elseif(sunit(1:5).eq.'WWCM2')then
          tc(i1)=t(i1)*dscs_ww*(exccm-wo)**4/wo/
     1          (1.d0-exp(-(hck*wo/temp)))
         elseif(sunit(1:6).eq.'PWCM2J')then
          tc(i1)=t(i1)*dscs_pw*(exccm-wo)**3/wo/
     1          (1.d0-exp(-(hck*wo/temp)))
c         elseif(sunit(1:6).eq.'PWCM5J')then ! same as  'WWCM6/h/c'
c          tc(i1)=t(i1)*dscs_pw/wo/(1.d0-exp(-(hck*wo/temp)))
         elseif(sunit(1:5).eq.'PPCM2')then
          tc(i1)=t(i1)*dscs_ww*exccm*(exccm-wo)**3/wo/
     1          (1.d0-exp(-(hck*wo/temp)))
c         elseif(sunit(1:5).eq.'PPCM6')then ! same as 'WWCM6'
c          tc(i1)=t(i1)*dscs_ww/wo/(1.d0-exp(-(hck*wo/temp)))
         else
          write(*,*)'ERROR: wrong sunit:',sunit(1:length(sunit))
         endif
         tc(i1)=tc(i1)*refrindex
        endif
        const=tc(i1)
      endif
4     s(i)=s(i)+sfc*const 
      if(ic.eq.100)then
        open(10,file='S.TAB')
        do i1=1,nm
808       format(g25.12,' ',g25.12)
          write(10,808)w(i1),tc(i1)
        enddo
        close(10)
        write(6,*)'S.TAB written'
        write(6,*)'  units:'
        if(sunit(1:6).eq.'KA4AMU')then ! K*A**4/amu unit - Gaussian16
          write(*,*)'Wavenumber(cm-1) '//
     1    'Polarizability_invariant_combinationr(K*A**4/amu)'
        elseif(sunit(1:5).eq.'WWCM6')then ! 
          write(*,*)'Wavenumber(cm-1) '//
     1    'R.s.c.s/w**4(cm**6/sr)'
        elseif(sunit(1:5).eq.'WWCM2')then
          write(*,*)'Wavenumber(cm-1) '//
     1    'R.s.c.s(cm**2/sr)'
        elseif(sunit(1:6).eq.'PWCM2J')then
          write(*,*)'Wavenumber(cm-1) '//
     1    'R.s.c.s(cm**2*photons/J/sr)'
        elseif(sunit(1:5).eq.'PPCM2')then
          write(*,*)'Wavenumber(cm-1) '//
     1    'R.s.c.s(cm**2/sr)'
        endif
        write(*,*)'  ROA intensity multiplied by 10**4'
        write(*,*)'S.PRN unit: as above * cm'
      endif
c
      call wrspec2('S.PRN',npoint,wsmin,dw,s)
      end
c     ==================================================================
      subroutine wrspec2(sn,nw,wmin,dw,spec)
      implicit none
      character*(*) sn
      integer*4 nw,iw
      real*8 wmin,dw,spec(*),w
      open(8,file=sn)
c     format according to the interval between the points:
      if(dabs(dw).gt.1.0d-1)then
       do 1 iw=1,nw
       w=wmin+dble(iw-1)*dw
1      write(8,805)w,spec(iw)
805    format(f12.2,' ',g25.12)
      else
       if(dabs(dw).gt.1.0d-3)then
        do 2 iw=1,nw
        w=wmin+dble(iw-1)*dw
2       write(8,806)w,spec(iw)
806     format(f14.4,' ',g25.12)
       else
        do 3 iw=1,nw
        w=wmin+dble(iw-1)*dw
3       write(8,807)w,spec(iw)
807     format(g25.12,' ',g25.12)
       endif
      endif
      close(8)
      write(6,*)sn,' written'
      return
      end
c     ==================================================================
      function sf(d,g,x,x0,ic)
      implicit none
      logical g
      integer*4 ic
      real*8 sf,d,x,x0,dd,ab,normL,normG
c      pi=4.0d0*atan(1.0d0)
c      spi=dsqrt(pi)
      normL=0.63661977236758138d0 ! 2/pi 
      normG=0.93943727869965132d0 ! 2*sqrt(log(2)/pi)
      dd=((x-x0)/d)**2
c     write(6,*)d,x,x0,g
      if(ic.eq.30)then
c      ORD, real parts, only Lorentz possible:
       ab=(x-x0)*(x+x0)*x**2*x0**2
       sf=ab/(ab**2+(d**2)**2)
      else
c      dispersion spectra, imaginary part
       if(g)then
        if(dd.lt.20.0d0)then
         sf=normG*exp(-dd)/d
c         sf=exp(-dd)/d/spi
        else
         sf=0.0d0
        endif
       else
        if(dd.lt.1.0d10)then
         sf=normL/(d*(dd+1.0d0))
c         sf=1.0d0/d/(dd+1.0d0)/pi
        else
         sf=0.0d0
        endif
       endif
      endif
      return
      end
c     ==================================================================
      subroutine report(s)
      character*(*)s
      write(6,*)s
      stop
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c        1         2         3         4         5         6         7  2
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      integer*4 function length(s)
      character s*(*)
      integer*4 i
      do i=len(s),1,-1
        if (s(i:i).ne.' ')goto 99999
      enddo
99999 length=i
      return
      end
c
      subroutine upper_case(s)
      character s*(*)
      integer*4 i,i1,length
      do i=1,length(s)
        i1=iachar(s(i:i))
        if(i1.ge.97.and.i1.le.122) s(i:i)=char(i1-32)
      enddo
      return
      end
c
      subroutine lower_case(s)
      character s*(*)
      integer*4 i,i1,length
      do i=1,length(s)
        i1=iachar(s(i:i))
        if(i1.ge.65.and.i1.le.90) s(i:i)=char(i1+32)
      enddo
      return
      end
c
      subroutine readramopt(s80,coefinv,sexpmode,sunit)
      integer*4 i1,length
      real*8 coefinv(5)
      character*80 s80
      character*30 sexpmode,sunit 
      i1=1
10    i1=i1+1
      if(i1.ge.len(s80)) then
        write(*,*)'Error 1 in readramopt:end of str reached,no / found'
        return
      endif
      if(s80(i1:i1).ne.'/') goto 10 
      sexpmode=s80(2:i1-1)
      sunit=s80(i1+1:length(s80))
      do i1=1,5
        coefinv(i1)=0.0d0
      enddo
      if(sexpmode(1:11).eq.'RAM_SCP_180')then
        coefinv(1)=180.d0
        coefinv(2)=28.d0
      elseif(sexpmode(1:12).eq.'RAM_DCPI_180')then
        coefinv(2)=24.d0
      elseif(sexpmode(1:13).eq.'RAM_DCPII_180')then
        coefinv(1)=180.d0
        coefinv(2)=4.d0
      elseif(sexpmode(1:11).eq.'ROA_SCP_180')then
        coefinv(4)=96.d0
        coefinv(5)=32.d0
      elseif(sexpmode(1:9).eq.'ROA_SCP_0')then
        coefinv(3)=720.d0
        coefinv(4)=16.d0
        coefinv(5)=-16.d0
      else
        write(*,*)'Error: unknown sexpmode: ',
     1   sexpmode(1:length(sexpmode))
      endif
      return
      end