[交流] 【求助】帮忙改写展宽程序已有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
      
      wbegin = 1000.0
      wend = 1800.0

      call readfreq(wavenr,IRint,weigenv,wbegin,wend)
      call convolution(wbegin,wend,wavenr,IRint,weigenv)
      stop
!     end program
      end
      
      subroutine readfreq(wavenr,IRint,weigenv,wbegin,wend)
!     print *,'beginning of subroutine...'
!     read in eigenfrequencies and IR intensities
      implicit none
      integer arbitrary
      parameter(arbitrary=3500)
      double precision wavenr(arbitrary),IRint(arbitrary)
      double precision wbegin,wend
      integer i,reason,weigenv
      double precision wavetmp, inttmp
      i=0
!     print *,'entering do statement'
      do
      read(*,*,IOSTAT=reason) wavetmp,inttmp
      if (reason>0) stop
      if (reason<0) return
      if(wavetmp.gt.wbegin .and. wavetmp .lt. wend) then
      i=i+1
      wavenr(i)=wavetmp
      IRint(i)=inttmp
!     print *,'eigenfreq: ',wavenr(i),' IR intensity: ',
!    &        IRint(i),' i = ',i
      weigenv=i
      endif
      end do
      return
      end subroutine
      
      subroutine convolution(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
      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
      
      do k=1,noofscanpoints
      write(*,*) w1(k),conv(k)
      enddo
      
      return
      end subroutine

      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

Originally posted by sci_papers at 2010-05-16 15:42:16:

再麻烦问下您,我看我程序的展宽语句怎么和Swizard 中的高斯或者洛伦茨展宽方程式不一样呀,难道我这个展宽程序不是那两个?之前别人告诉我说这个就是高斯展宽.可高斯展宽是
[img]http://pic.muchong.com/201005/16 ...