[½»Á÷] ¡¾ÇóÖú¡¿°ïæ¸Äдչ¿í³ÌÐò ÒÑÓÐ3È˲ÎÓë

      program CONVOLUTE
      implicit none
      integer arbitrary
      parameter(arbitrary=3500)
      double precision wavenr(arbitrary),IRint(arbitrary)
      integer weigenv
!     wavenumber start and finish
      double precision wbegin,wend
      character*80 input, output
      
      write(*,'(a,$)') 'input file name: '
      read(*, '(a)') input
      write(*,'(a,$)') 'output file name: '
      read(*, '(a)') output
      write(*,'(a,$)') 'start of wavenumber: '
      read(*,*) wbegin
      write(*,'(a,$)') 'end of wavenumber: '
      read(*,*) wend
      
      call readfreq(input,wavenr,IRint,weigenv,wbegin,wend)
      
      call convolution(output,wbegin,wend,wavenr,IRint,weigenv)
      stop
      end
      
      subroutine readfreq(fname,wavenr,IRint,weigenv,wbegin,wend)
!      read in eigenfrequencies and IR intensities
      implicit none
      integer arbitrary
      parameter(arbitrary=3500)
      double precision wavenr(arbitrary),IRint(arbitrary)
      double precision wbegin,wend
      character*80 fname
!      print *,'opening the file'
      integer i,reason,weigenv
      double precision wavetmp, inttmp
      open (unit=3,file=fname,status='OLD')
      i=0
!      print *,'entering do statement'
      do
      read(3,*,IOSTAT=reason) wavetmp,inttmp
      if (reason>0) then
        write(*,*) 'reason>0, i.e. somethings wrong'
        stop
      endif
      if (reason<0) then
        write(*,*) 'reason<0, reached the end of file'
        close(3)
        return
      endif
      if(wavetmp.gt.wbegin .and. wavetmp .lt. wend) then
        i=i+1
        wavenr(i)=wavetmp
        IRint(i)=inttmp
        write(*,*) 'eigenfreq: ',wavenr(i),' IR intensity: ',
     &             IRint(i),' i = ',i
        weigenv=i
      endif
      end do
      return
!      print *,'end of subroutine'
      end subroutine
      
      
      subroutine convolution(fname,wbegin,wend,wavenr,IRint,weigenv)
      implicit none
      integer bignumber,arbitrary
      parameter(bignumber=10000000,arbitrary=3500)
      double precision wbegin,wend,wavenr(arbitrary),IRint(arbitrary)
      integer noofscanpoints,n,m,k,weigenv
      parameter(noofscanpoints=100000)
      double precision deltaw,w1(bignumber),conv(bignumber)
      double precision mygamma,wave2ev
      character*80 fname
      parameter(wave2ev=8065.46)
!      mygamma=4.0d-3*wave2ev
      mygamma=10.0
!      mygamma=4.0d-3 in eV unit, convert to wavenumber by multiple wave2ev
      
      deltaw=(wend-wbegin)/dfloat(noofscanpoints)
      do n=1,noofscanpoints
      w1(n)=wbegin+n*deltaw
      conv(n)=0.0
      enddo
      do m=1,weigenv
      do k=1,noofscanpoints
!      print *,'m=',m,'k=',k
      conv(k)=conv(k)+IRint(m)/((w1(k)-wavenr(m))**2+mygamma**2)
!      print *,'w1=',w1(k),' conv(k)=',conv(k)
      enddo
      enddo
      
      open(unit=1,file=fname,status='unknown')
      do k=1,noofscanpoints
      write(1,*) w1(k),conv(k)
      enddo
      close(1)
      
      return
      end subroutine