program nvtau
c     supplementary model of periodic twisted systems
      implicit none
      character*80 s80
      real*8 t,V,w0,AN,rq,pi,f1,ds,rs,sr1p,sr1m,maxd,maxr,
     1sr2p,sr2m,l,ecm,wsmin,wsmax,fwhm,debye,clight,FD,FC,AMU,CM,BOHR,
     1WR,si1p,si1m,si2m,si2p,mxi,mxr,uxi,uxr,myi,myr,uyi,uyr
      integer*4 N,q,np,i
      real*8,allocatable::sd(:,:)
      logical lg
      lg=.false.
      fwhm=15.0d0
      w0=1650.0d0
      wsmin=1550.0d0
      wsmax=1750.0d0
      np=201
      BOHR=0.52917705993d0
      l=5.08d0/BOHR/4.0d0
      pi=4.0d0*datan(1.0d0)
      debye=2.541732d0
      AMU=1822.887D0
      CM=219474.63d0
      clight=137.03604d0
      FD =debye*dsqrt(CM/AMU/2.0d0)
      FC=debye*dsqrt(2.0d0/CM/AMU)/clight
c      spectra: 1 absorbtion
c               2 CD
      allocate(sd(2,np))
      do 2 i=1,np
      sd(1,i)=0.0d0
2     sd(2,i)=0.0d0
      call getarg(1,s80)
      read(s80,*)N
      call getarg(2,s80)
      read(s80,*)V
      call getarg(3,s80)
      read(s80,*)t
      t=pi*t/180.0d0
      AN=dble(N)

      OPEN(18,FILE='DOG.TAB')
      WRITE (18,501)
  501 FORMAT (/,'          FREQ.     ',
     1'  D [D^2]     ','  R [D^-2]DO         R - CO      I[km/mol]',
     2/,'-----------------------------------------------------------')

      do 1 q=0,N-1
      rq=dble(q)
      f1=2.0d0*pi*rq/AN
      ecm=w0+2.0d0*V*cos(f1)
      WR=dsqrt(ecm)
      call s1s(f1+t,sr1p,si1p,AN)
      call s1s(f1-t,sr1m,si1m,AN)
      call s2s(f1+t,sr2p,si2p,AN)
      call s2s(f1-t,sr2m,si2m,AN)
      uxr=sr1p+sr1m
      uxi=si1p+si1m
      uyr=  si1p-si1m
      uyi=-(sr1p-sr1m)
      mxr= (si2m-si2p)*l
      mxi=-(sr2m-sr2p)*l
      myr= (sr2m+sr2p)*l
      myi= (si2m+si2p)*l
      write(6,*)q
      write(6,*)'sr1p si1p',sr1p,si1p
      write(6,*)'sr1m si1m',sr1m,si1m
      write(6,*)'ux:',uxr,uxi,'uy:',uyr,uyi
      write(6,*)'sr2p si2p',sr2p,si2p
      write(6,*)'sr2m si2m',sr2m,si2m
      write(6,*)'mx:',mxr,mxi,'my:',myr,myi
      ds=uxr*uxr+uxi*uxi+uyr*uyr+uyi*uyi
      rs=uxr*mxr+uxi*mxi+uyr*myr+uyi*myi
      ds=ds*FD*FD/WR/WR/dble(N)
      rs=rs*FD*FC/WR*WR/dble(N)
      call sspread(ecm,ds,rs,sd,np,wsmin,wsmax,fwhm,lg)
1     WRITE(18,523)q,ecm,ds,rs
523   FORMAT(I5,f16.4,4G12.4)
      WRITE(18,1819)
1819  FORMAT('-------------------------------------------------')
      close(18)
      call wspectra('d.prn',wsmin,wsmax,sd,np,1)
      call wspectra('r.prn',wsmin,wsmax,sd,np,2)
      maxr=0.0d0
      maxd=0.0d0
      do 3 i=1,np
      if(sd(1,i).gt.maxd)maxd=sd(1,i)
3     if(dabs(sd(2,i)).gt.dabs(maxr))maxr=sd(2,i)
      write(6,600)maxd,maxr,maxr/maxd
600   format(3e13.4)
      end

      subroutine wspectra(f,wsmin,wsmax,s,np,j)
      implicit none
      real*8 s(2,np),wsmin,wsmax,w,dw
      integer*4 i,np,j
      character*(*) f
      OPEN(30,FILE=f)
      dw=(wsmax-wsmin)/dble(np-1)
      w=wsmin-dw
      do 1 i=1,np
      w=w+dw
 1    write(30,300)w,s(j,i)
300   format(f12.2,' ',g12.4)
      close(30)
      return
      end

      subroutine period(k,kk)
      implicit none
      real*8 k,kk,tpi
      integer*4 n
      tpi=8.0d0*datan(1.0d0)
      n=int(kk/tpi)
      k=kk-dble(n)*tpi
      if(dabs(k-tpi).lt.dabs(k))k=k-tpi
      if(dabs(k+tpi).lt.dabs(k))k=k+tpi
      return
      end


      subroutine s1s(kk,r,i,N)
      implicit none
      real*8 k,kk,r,i,N,down,l
      l=1.0d0
      call period(k,kk)
      if(dabs(k).gt.1.0d-6)then
c      sum(I=0..N-1) 1/sqrt(N) exp[i(2piq/N+kt)I]
       down=2.0d0*dsqrt(N)*(l-cos(k))
       r=(l+cos(k*(N-l))-cos(k)-cos(k*N))/down
       i=(  sin(k*(N-l))+sin(k)-sin(k*N))/down
      else
       r=dsqrt(N)
       i=0.0d0
      endif
      return
      end

      subroutine s2s(kk,r,i,N)
      implicit none
      real*8 k,kk,r,i,N,down,l
      call period(k,kk)
      l=1.0d0
      if(dabs(k).gt.1.0d-6)then
c      sum(I=0..N-1) 1/sqrt(N) I exp[i(2piq/N+kt)I]
       down=4.0d0*dsqrt(N)*(l-cos(k))**2
       r=(    (N-l)*cos(k*(N+l))-(3.0d0*N-2.0d0)*cos(k*N)
     1 +(3.0d0*N-l)*cos(k*(N-l))-N*cos(k*(N-2.0d0))
     1 +2.0d0*cos(k)-2.0d0)/down
       i=(    (N-l)*sin(k*(N+l))-(3.0d0*N-2.0d0)*sin(k*N)
     1 +(3.0d0*N-l)*sin(k*(N-l))-N*sin(k*(N-2.0d0)))/down
      else
       r=dsqrt(N)*(N-l)/2.0d0
       i=0.0d0
      endif
      return
      end

      subroutine sspread(ecm,d1,d2,s,np,wsmin,wsmax,fwhm,lg)
      implicit none
      real*8 s(2,np),wsmin,wsmax,w,dw,ecm,pd,fwhm,sf,d1,
     1d2,dd(2),prod
      integer*4 i,np
      logical lg
      dd(1)=d1*108.7d0
      dd(2)=d2*435.0d0
      dw=(wsmax-wsmin)/dble(np-1)
      w=wsmin-dw
      if(lg)then
       pd=fwhm/2.0d0/dsqrt(log(2.0d0))
      else
       pd=fwhm/2.0d0
      endif
      do 1 i=1,np
      w=w+dw
      prod=w*sf(pd,lg,w,ecm,1)
      s(1,i)=s(1,i)+prod*dd(1)
 1    s(2,i)=s(2,i)+prod*dd(2)
      return
      end
      
      function sf(d,g,x,x0,ic)
      implicit none
      logical g
      integer*4 ic
      real*8 pi,spi,sf,d,x,x0,dd,ab
      pi=4.0d0*atan(1.0d0)
      spi=dsqrt(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=exp(-dd)/d/spi
        else
         sf=0.0d0
        endif
       else
        if(dd.lt.1.0d10)then
         sf=1.0d0/d/(dd+1.0d0)/pi
        else
         sf=0.0d0
        endif
       endif
      endif
      return
      end