program smooth
      implicit none
c     one-purpusse program to smooth 2D PES in angles
      real*8,allocatable::fi(:),psi(:),g(:),ak(:),sum(:)
      integer*4 i,n,ifi,fimin,fimax,delfi,ipsi,psimin,psimax,delpsi,
     1nbf,kmax,i1,i2,iw,jw,is,js,w,ibf,ii
      real*8 delfio,delpsio,wi,wj,f1,f2,pi,gmin
      logical lex
      character*4 key

      pi=4.0d0*datan(1.0d0)
      n=0
c     read angles
      open(9,file='ANGLE.LST')
1     read(9,*,end=99,err=99)
      n=n+1
      goto 1
99    rewind 9
      allocate(fi(n),psi(n),g(n))
      do 2 i=1,n
2     read(9,*)fi(i),fi(i),psi(i)
      close(9)
      write(6,*)n,' angle values'

      open(9,file='test.txt')
      do 3 i=1,n
3     read(9,*)g(i),g(i),g(i)
      close(9)
      write(6,*)n,' energies'

c     maximum frequency:
      kmax=3
      delfio=18.0d0
      delpsio=18.0d0
      delfi=5
      delpsi=5
      psimin=-185
      fimin=-185
      fimax=185
      psimax=185
      inquire(file='SMOOTH.OPT',exist=lex)
      if(lex)then
       open(9,file='SMOOTH.OPT')
4      read(9,90,end=98,err=98)key
90     format(a4)
       if(key.eq.'KMAX')read(9,*)kmax
       if(key.eq.'DFI ')read(9,*)delfio
       if(key.eq.'DPSI')read(9,*)delpsio
       if(key.eq.'IFI ')read(9,*)delfi
       if(key.eq.'IPSI')read(9,*)delpsi
       if(key.eq.'PSII')read(9,*)psimin
       if(key.eq.'PSIA')read(9,*)psimax
       if(key.eq.'FII ')read(9,*)fimin
       if(key.eq.'FIA ')read(9,*)fimax
       goto 4
98     close(9)
       write(6,*)' SMOOTH.OPT read'
      endif
      nbf=(kmax+kmax+1)**2
      write(6,60)kmax,nbf,delfio,delpsio,delfi,delpsi,fimin,fimax,
     1psimin,psimax
60    format(' kmax ',i4,/,
     1' number of basis functions ',i4,/,
     2' delfi, delpsi ',2f10.2,/,
     3' new pes:',/,
     4' delfi, delpsi ',2i10,/,
     4' fimin fimax ',2i10,/,
     4' psimin psimax ',2i10,/)
      allocate(ak(nbf))
      ibf=0
      do 5 i1=0,2*kmax
c     i1:0    1    2    3    4
c     cos0 sin1 cos1 sin2 cos2 ...
      iw=(i1+1)/2
      wi=dble(iw)
      if(mod(i1,2).eq.0)then
       is=0
      else
       is=1
      endif
      do 5 i2=0,2*kmax
      jw=(i2+1)/2
      wj=dble(jw)
      if(mod(i2,2).eq.0)then
       js=0
      else
       js=1
      endif
      ibf=ibf+1
      ak(ibf)=0.0d0
      do 51 w=1,n
      if(is.eq.1)then
       f1=sin(wi*fi(w)*pi/180.0d0)
      else
       f1=cos(wi*fi(w)*pi/180.0d0)
      endif
      if(js.eq.1)then
       f2=sin(wj*psi(w)*pi/180.0d0)
      else
       f2=cos(wj*psi(w)*pi/180.0d0)
      endif
51    ak(ibf)=ak(ibf)+g(w)*f1*f2
5     ak(ibf)=ak(ibf)/pi**2*delfio*delpsio*(pi/180.0d0)**2

      allocate(sum(((fimax-fimin)/delfi+1)*((psimax-psimin)/delpsi+1)))
      ii=0
      gmin=0.0d0
      do 6 ifi=fimin,fimax,delfi
      do 6 ipsi=psimin,psimax,delpsi
      ii=ii+1
      sum(ii)=0.0d0
      ibf=0
      do 7 i1=0,2*kmax
      iw=(i1+1)/2
      wi=dble(iw)
      if(mod(i1,2).eq.0)then
       is=0
      else
       is=1
      endif
      do 7 i2=0,2*kmax
      jw=(i2+1)/2
      wj=dble(jw)
      if(mod(i2,2).eq.0)then
       js=0
      else
       js=1
      endif
      ibf=ibf+1
      if(is.eq.1)then
       f1=sin(wi*dble(ifi)*pi/180.0d0)
      else
       f1=cos(wi*dble(ifi)*pi/180.0d0)
      endif
      if(js.eq.1)then
       f2=sin(wj*dble(ipsi)*pi/180.0d0)
      else
       f2=cos(wj*dble(ipsi)*pi/180.0d0)
      endif
7     sum(ii)=sum(ii)+f1*f2*ak(ibf)
      if(ii.eq.1)gmin=sum(ii)
      if(sum(ii).lt.gmin)gmin=sum(ii)
6     continue

      open(9,file='pes.txt')
      ii=0
      do 8 ifi=fimin,fimax,delfi
      do 8 ipsi=psimin,psimax,delpsi
      ii=ii+1
8     write(9,900)ifi,ipsi,sum(ii)-gmin
900   format(2i10,f12.5)
      close(9)
      write(6,*)' smoothed function in pes.txt'

      stop

      end