program ffit4
      implicit none
      integer*4 nb,ib,NP,i,NExc,iex,ip,ic,j,nd
      real*8 f,ff,dV,er
      reaL*8,allocatable::c(:),n(:),x0(:),LX(:),dx(:),x(:)
      integer*4,allocatable::state(:),si(:)

c     f(fi,psi)=sum (k) ck gk(fi,psi)
      open(9,file='FFIT4.INP')
      read(9,*)nb,nd
      allocate(x0(nd),dx(nd),LX(nd),c(nb**nd),n(nb**nd),x(nd),
     1state(nd),si((nb-1)*nd))
      read(9,*)(x0(i),i=1,nd)
      read(9,*)(dx(i),i=1,nd)
      read(9,*)(LX(i),i=1,nd)
      read(9,*)NP

      dV=1.0d0
      do 31 i=1,nd
31    dV=dV*dx(i)
      
      do 30 ib=1,nb**nd
      c(ib)=0.0d0
30    n(ib)=0.0d0

      do 2 ip=1,NP
      read(9,*)x(1),(x(j),j=1,nd),f

c     ground state:
      ib=1
      call vzi(state,nd)
      call digest(state,ib,n,nd,c,x,x0,LX,dV,f)

c     states Nexc excited:
      do 3 Nexc=1,(nb-1)*nd
c     distribute Nexc excitations upon centers 1..nd
c     initial indices - all to the first center:
      do 42 iex=1,Nexc
42    si(iex)=1
5     call vzi(state,nd)
      do 4 iex=1,Nexc
4     state(si(iex))=state(si(iex))+1
c     check that the state is allowed:
      do 11 I=1,nd
11    if(state(I).gt.nb-1)goto 12

      ib=ib+1
      if(ib.gt.nb**nd)call report('ib > nb**nd')
      call digest(state,ib,n,nd,c,x,x0,LX,dV,f)
c
c     find index to be changed
12    do 8 ic=Nexc,1,-1
8     if(si(ic).lt.nd)goto 9
      goto 3
9     do 10 i=ic+1,Nexc
10    si(i)=si(ic)+1
      si(ic)=si(ic)+1
c
      goto 5
c
3     continue
C
2     continue

      close(9)

      write(6,*)np,' points'
      write(6,*)nd,' variables'
      write(6,*)nb**nd,' basis functions'

      open(9,file='FFIT4.PAR')
      write(9,608)nd,nb**nd
608   format(i10,' variables',/,
     1i10,' basis functions, coefficients, norms:')
      ib=1
      call vzi(state,nd)
      call wrstate(ib,state,nd,LX,c,n)
      do 103 Nexc=1,(nb-1)*nd
      do 142 iex=1,Nexc
142   si(iex)=1
105   call vzi(state,nd)
      do 104 iex=1,Nexc
104   state(si(iex))=state(si(iex))+1
      do 111 I=1,nd
111   if(state(I).gt.nb-1)goto 112
      ib=ib+1
      call wrstate(ib,state,nd,LX,c,n)
112   do 108 ic=Nexc,1,-1
108   if(si(ic).lt.nd)goto 109
      goto 103
109   do 110 i=ic+1,Nexc
110   si(i)=si(ic)+1
      si(ic)=si(ic)+1
      goto 105
103   continue
      close(9)
      write(6,*)'Coefficients written into FFIT4.PAR'

      open(19,file='FIT.TXT')
      er=0.0d0
      open(9,file='FFIT4.INP')
      read(9,*)
      read(9,*)
      read(9,*)
      read(9,*)
      read(9,*)
      do 40 ip=1,NP
      read(9,*)x(1),(x(j),j=1,nd),f
      ff=0.0d0

      ib=1
      call vzi(state,nd)
      call addstate(state,nd,x,x0,LX,c,ff,ib)
      do 203 Nexc=1,(nb-1)*nd
      do 242 iex=1,Nexc
242   si(iex)=1
205   call vzi(state,nd)
      do 204 iex=1,Nexc
204   state(si(iex))=state(si(iex))+1
      do 211 I=1,nd
211   if(state(I).gt.nb-1)goto 212
      ib=ib+1
      call addstate(state,nd,x,x0,LX,c,ff,ib)
212   do 208 ic=Nexc,1,-1
208   if(si(ic).lt.nd)goto 209
      goto 203
209   do 210 i=ic+1,Nexc
210   si(i)=si(ic)+1
      si(ic)=si(ic)+1
      goto 205
203   continue

      er=er+(f-ff)**2
      do 41 j=1,nd
41    write(19,403)x(j)
403   format(f10.2,$)
40    write(19,400)f,ff
400   format(4f12.4)
      er=dsqrt(er/dble(NP))
      close(9)
      close(19)
      write(6,*)'Old and fitted values written into FIT.TXT'
      write(6,401)er
401   format('RMS Error: ',f12.6)
      end

      function fik(k,x,x0,L)
c       K:  1      2           3           4           5           6
c     fik:  const  sin(2pix/L) cos(2pix/L) sin(4pix/L) cos(4pix/L) sin(6pix/L)
      implicit none
      integer*4 k
      real*8 x,x0,L,y,sf,fik,pi
      pi=4.0d0*datan(1.0d0)
      y=x-x0
      if(k.eq.1)then
       sf=1.0d0/dsqrt(L)
       fik=sf
      else
       sf=dsqrt(2.0d0/L)
       if(mod(k,2).ne.0)then
        fik=sf*cos(dble(2*((k-1)/2))*pi*y/L)
       else
        fik=sf*sin(dble(2*(k/2))*pi*y/L)
       endif
      endif
      return
      end
c     ============================================================
      subroutine setk(k,fs,snk,pixl,a,L)
      implicit none
      integer*4 k,nk
      real*8 a,L
      character*5 fs
      character*3 snk
      character*7 pixl
      if(k.eq.1)then
       fs=' one '
       snk='   '
       pixl='       '
       a=1.0d0/dsqrt(L)
      else
       a=dsqrt(2.0d0/L)
       if(mod(k,2).ne.0)then
        fs=' cos('
        nk=2*((k-1)/2)
       else
        fs=' sin('
        nk=2*(k/2)
       endif
       write(snk,900)nk
900    format(i3)
       if(nk.eq.1)snk='   '
       pixl=' pix/L)'
      endif
      return
      end

      subroutine vzi(s,n)
      integer*4 s(*),n,i
      DO  1 i=1,n
1     s(i)=0
      return
      end

      subroutine digest(state,ib,n,nd,c,x,x0,LX,dV,f)
      implicit none
      integer*4 state(*),ib,i,nd
      real*8 n(*),c(*),p,fik,x(*),x0(*),LX(*),dV,f
      p=1.0d0
      do 131 i=1,nd
131   p=p*fik(state(i)+1,x(i),x0(i),LX(i))
      n(ib)=n(ib)+p*p*dV
      c(ib)=c(ib)+f*p*dV
      return
      end

      subroutine wrstate(ib,state,nd,LX,c,n)
      implicit none
      integer*4 ib,state(*),nd,i
      real*8 LX(*),a,c(*),n(*),ax
      character*5 fs1
      character*3 snk1
      character*7 pixl1
      write(9,301)ib
301   format(i4,1x,$)
      a=1.0d0
      do 112 i=1,nd
      call setk(state(i)+1,fs1,snk1,pixl1,ax,LX(i))
      a=a*ax
112   write(9,609)fs1//snk1//pixl1//' '
609   format((a),$)
      write(9,619)a*c(ib),n(ib)
619   format(2f12.4)
      return
      end

      subroutine addstate(state,nd,x,x0,LX,c,ff,ib)
      implicit none
      integer*4 state(*),nd,i,ib
      real*8 x(*),x0(*),LX(*),c(*),p,fik,ff
      p=1.0d0
      do 302 i=1,nd
302   p=p*fik(state(i)+1,x(i),x0(i),LX(i))
      ff=ff+c(ib)*p
      return
      end

      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end