program ffit
      implicit none
      integer*4 nb,ib,k1,k2,NP,i,jb,l1,l2
      character*5 fs1,fs2
      character*3 snk1,snk2
      character*7 pixl1,pixl2
      real*8 x0,y0,LX,LY,x,y,f,fik,ff,dx,dy,dV,er,a,ax,ay
      reaL*8,allocatable::c(:),n(:),b(:),m(:,:),ai(:,:)
      logical ls
c     f(fi,psi)=sum (k) ck gk(fi,psi)
      open(9,file='FFIT.INP')
      read(9,*)nb
      ls=nb.lt.0
      if(nb.gt.0)then
       write(6,*)'integration'
      else
       nb=-nb
       write(6,*)'least square fit'
      endif
      read(9,*)x0,y0
      read(9,*)dx,dy
      read(9,*)LX,LY
      read(9,*)NP

      if(ls)then
       allocate(m(nb**2,nb**2),b(nb**2),ai(nb**2,nb**2))
       do 6 ib=1,nb**2
       b(ib)=0.0d0
       do 6 jb=1,nb**2
6      m(ib,jb)=0.0d0
      endif
      allocate(c(nb**2),n(nb**2))
      do 3 ib=1,nb**2
      c(ib)=0.0d0
3     n(ib)=0.0d0
      dV=dx*dy
      do 2 i=1,NP
      read(9,*)x,y,f
      if(ls)then
       ib=0
       do 7 k1=1,nb
       do 7 k2=1,nb
       ib=ib+1
       b(ib)=b(ib)+fik(k1,x,x0,LX)*fik(k2,y,y0,LY)*f
       jb=0
       do 7 l1=1,nb
       do 7 l2=1,nb
       jb=jb+1
7      m(ib,jb)=m(ib,jb)+fik(k1,x,x0,LX)*fik(k2,y,y0,LY)
     1 *                 fik(l1,x,x0,LX)*fik(l2,y,y0,LY)
      endif
      ib=0
      do 2 k1=1,nb
      do 2 k2=1,nb
      ib=ib+1
      n(ib)=n(ib)+  fik(k1,x,x0,LX)**2*fik(k2,y,y0,LY)**2*dV
2     c(ib)=c(ib)+f*fik(k1,x,x0,LX)*fik(k2,y,y0,LY)*dV
      close(9)

      if(ls)then
       call INV(m,ai,nb**2,1.0d-25)
       do 8 ib=1,nb**2
       c(ib)=0.0d0
       do 8 jb=1,nb**2
8      c(ib)=c(ib)+ai(ib,jb)*b(jb)
      endif

      write(6,*)NP,' points'
      write(6,*)nb**2,' 1D basis functions, coefficients, norms:'
      ib=0
      do 1 k1=1,nb
      call setk(k1,fs1,snk1,pixl1,ax,LX)
      do 1 k2=1,nb
      ib=ib+1
      call setk(k2,fs2,snk2,pixl2,ay,LY)
      a=ax*ay
      if(k1.eq.1.and.k2.eq.1)fs2='     '
      if(k1.eq.1.and.k2.gt.1)fs1='     '
      if(k1.gt.1.and.k2.eq.1)fs2='     '
1     write(6,609)ib,fs1//snk1//pixl1//' '//fs2//snk2//pixl2,
     1a*c(ib),n(ib)
609   format(i4,1x,(a),2f14.6)

      write(6,*)'fit:'
      er=0.0d0
      open(9,file='FFIT.INP')
      read(9,*)
      read(9,*)
      read(9,*)
      read(9,*)
      read(9,*)
      do 4 i=1,NP
      read(9,*)x,y,f
      ff=0.0d0
      ib=0
      do 5 k1=1,nb
      do 5 k2=1,nb
      ib=ib+1
5     ff=ff+c(ib)*fik(k1,x,x0,LX)*fik(k2,y,y0,LY)
      er=er+(f-ff)**2
4     write(6,400)x,y,f,ff
400   format(4f12.4)
      er=dsqrt(er/dble(NP))
      close(9)
      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 INV(A,AI,N,TOL)                                     
      IMPLICIT none
      integer*4 N,oo,ii,jj,iw,io,kk,i2,j2
      REAL*8 TOL, A(N,N),AI(N,N),w                           
      real*8, allocatable ::E(:,:)

      allocate(E(N,2*N))

      DO 1  ii=1,N                                                 
      DO 1  jj=1,N                                                 
      e(ii,jj)=a(ii,jj)
      E(II,JJ+N)=0.0D0                                             
1     if (ii.EQ.jj)e(ii,jj+N)=1.0D0                                

      DO 2 ii=1,N-1                                                
      iw=ii
      if (DABS(e(iw,iw)).LT.TOL) then                                 
       DO 3 io=iw+1,N                                              
       oo=IO
       if (io.GT.N)oo=io-N                                    
3      if (DABS(e(oo,iw)).GE.TOL) goto 11                             
       write(6,*)'Inverse cannot be done'
       stop
11     CONTINUE                                                       
       DO 4 kk=1, 2*N                                              
       w=e(iw,kk)                                                     
       e(iw,kk)=e(oo,kk)                                              
4      e(oo,kk)=w                                                     
      ENDIF                                                           
      DO 2 jj=ii+1,N                                               
      DO 6 kk=ii+1, 2*N
6     e(jj,kk)=e(jj,kk)-e(ii,kk)*e(jj,ii)/e(ii,ii)                    
2     e(jj,ii)=0.0D0                                                  

      DO 7  ii=1, N-1
      i2=N-ii+1                                                    
      DO 7 jj=1,i2-1                                                  
      j2=i2-jj                                                        
      DO 9 kk=1, N                                                 
9     e(j2,kk+N)=e(j2,kk+N)-e(i2,kk+N)*e(j2,i2)/e(i2,i2)
7     e(j2,i2)=0.0D0
 
      DO 10 ii=1,N                                                 
      DO 10 jj=1,N                                                
10    AI(ii,jj)=e(ii,jj+N)/e(ii,ii)                                

      deallocate(E)
      RETURN                                                          
      END

      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='const'
       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