program scub
c     FT smoothing cub file
      implicit none
      character*40 fn
      character*80 s80
      integer*4 inat,nat,n1,n2,n3,ia,ix,iy,iz,i,kmax,
     1nbf,k,l,m,nb1,i1,iargc,i2
      real*8 ax,ay,az,px,py,pz,pi,dx,dy,dz,xp,yp,zp,
     1dv,g1,gk,gl,xl,yl,zl,one,two
      real*8,allocatable::x(:,:),orb(:),ckl(:),vx(:),vy(:),vz(:)
      if(iargc().ne.2)then
       write(6,*)' Usage: scub file.cub kmax'
       stop
      endif
      one=1.0d0
      two=2.0d0
      pi=4.0d0*datan(1.0d0)
      call getarg(1,fn)
      open(8,file=fn)
      open(9,file='s'//fn)
      call getarg(2,fn)
      read(fn,*)kmax
      nb1=kmax+kmax+1
      nbf=nb1**3
      allocate(ckl(nbf))
      read(8,2000)s80
2000  format(a80)
      write(9,2000)s80
      read(8,2000)s80
      write(9,2000)s80

      read(8,*)inat,ax,ay,az
      nat=iabs(inat)
      allocate(x(3,nat))
      read(8,*)n1,px
      read(8,*)n2,py,py
      read(8,*)n3,pz,pz,pz
      xl=dble(n1)*px
      yl=dble(n2)*py
      zl=dble(n3)*pz
      dx=2.0d0*pi/xl
      dy=2.0d0*pi/yl
      dz=2.0d0*pi/zl
      dv=px*py*pz

      write(9,2001)nat,ax,ay,az
2001  format(i5,3f12.6)
      write(9,2001)n1,px,0.0d0,0.0d0
      write(9,2001)n2,0.0d0,py,0.0d0
      write(9,2001)n3,0.0d0,0.0d0,pz

      do 902 ia=1,nat
      read(8,2000)s80
      read(s80(18:80),*)(x(i,ia),i=1,3)
902   write(9,2200)s80(1:17),(x(i,ia),i=1,3)
2200  format(a17,3f12.6)
c     for reading orbitals:
c     read(8,905)j,(idum,i=1,j)
c05   format(10i5)

      allocate(orb(n1*n2*n3))
      do 906 ix=1,n1
      do 906 iy=1,n2
906   read( 8,4000)
     1(orb(iz+n3*n2*(ix-1)+n3*(iy-1)),iz=1,n3)

      do 917 k=1,nbf
917   ckl(k)=0.0d0

      allocate(vx(n1*nb1),vy(n2*nb1),vz(n2*nb1))
      do 920 k=0,2*kmax
      xp=px
      do 921 ix=1,n1
      xp=xp+px
921   vx(ix+k*n1)=g1(dsqrt(one/xl),dsqrt(two/xl),k,xp*dx)
      yp=py
      do 922 iy=1,n2
      yp=yp+py
922   vy(iy+k*n2)=g1(dsqrt(one/yl),dsqrt(two/yl),k,yp*dy)
      zp=-pz
      do 920 iz=1,n3
      zp=zp+pz
920   vz(iz+k*n3)=g1(dsqrt(one/zl),dsqrt(two/zl),k,zp*dz)

      do 916 k=0,2*kmax
      do 916 ix=1,n1
      gk=vx(ix+k*n1)
      do 916 l=0,2*kmax
      do 916 iy=1,n2
      gl=gk*vy(iy+l*n2)
      i1=n3*(n2*(ix-1)+iy-1)
      do 916 m=0,2*kmax
      i2=m+1+nb1*(l+nb1*k)
      do 916 iz=1,n3
916   ckl(i2)=ckl(i2)+orb(iz+i1)*gl*vz(iz+m*n3)

      do 918 k=1,nbf
918   ckl(k)=ckl(k)*dv

      do 930 iz=1,n1*n2*n3
930   orb(iz)=0.0d0

      do 927 ix=1,n1
      do 927 k=0,2*kmax
      gk=vx(ix+k*n1)
      do 927 iy=1,n2
      i1=n3*(n2*(ix-1)+iy-1)
      do 927 l=0,2*kmax
      gl=gk*vy(iy+l*n2)
      i2=nb1*(l+nb1*k)
      do 927 iz=1,n3
      do 927 m=0,2*kmax
927   orb(iz+i1)=orb(iz+i1)+ckl(m+1+i2)*gl*vz(iz+m*n3)

      do 928 ix=1,n1
      do 928 iy=1,n2
      i1=n3*(n2*(ix-1)+iy-1)
928   write(9,4000)(orb(iz+i1),iz=1,n3)
4000  format(6E13.5)

      close(8)
      close(9)
      end

      function  g1(o,t,k,fi)
      implicit none
      real*8 o,t,gk,wj,g1,fi
      integer*4 k,jw,js
      jw=(k+1)/2
      wj=dble(jw)
      if(mod(k,2).eq.0)then
       js=0
      else
       js=1
      endif
      if(js.eq.1)then
       gk=sin(wj*fi)
      else
       gk=cos(wj*fi)
      endif
      if(k.eq.0)then
       g1=gk*o
      else
       g1=gk*t
      endif
      return
      end