*23456789012345678901234567890123456789012345678901234567890123456789012
      function funmin(p,n)
      implicit double precision (a-h,o-z)
      DIMENSION P(N)
      save
      include 'para'
*      parameter (llw=5,llm=1,lls=1,nnd=1,nn=llw+lls+2+nnd*(llm+1)) 
* mm is the number of the inpout data points
*      parameter (mm=50)
      dimension polm(nnd),x(nnd),gfw(nn)
      dimension xset(mm,nnd),y(mm)
      data itric/0/
*
      if (itric.eq.0) then
         itric=1
      if (n.ne.nn) then
         write(*,*)' SUBSUM: # of Variables N not equal to parameter nn'
         write(*,*)' N = ',N,' nn = ', nn
         stop
      endif
*
* Fetch Data: 
*******      USER PROVIDED SUBROUTINE:   FETDAT 
        call fetdat(xset,mm,nnd,y)
*******      
      endif
*
      funmin=0.
      do i=1,mm
        do j=1,nnd
           x(j)=xset(i,j)
        enddo
      call fwann(x,nnd,p,nn,llw,llm,lls,polw,polm,pols,fwnn,gfw)
        RES= fwnn-y(i)
        funmin=funmin+res**2
      enddo
      call getfla(1,xf1)
      if (xf1.eq.1.) then
       call setfla(1,0.)
       itap = nunit()
       open(itap,file='res')
       rewind itap
      do i=1,mm
        do j=1,nnd
           x(j)=xset(i,j)
        enddo
      call fwann(x,nnd,p,nn,llw,llm,lls,polw,polm,pols,fwnn,gfw)
      write(itap,*) (x(j),j=1,nnd), fwnn, y(i)
      enddo
      close(itap)
      endif

      end
*=======================================================================
      subroutine fetdat(xset,mm,nnd,y)
      implicit double precision (a-h,o-z)
      dimension xset(mm,nnd),y(mm)
      nt=nunit()
      open(nt,file='datal')
      rewind(nt)
      do m=1,mm
        read(nt,*)(xset(m,j),j=1,nnd)
      enddo
      close(nt)
      ntt=nunit()
      open(ntt,file='datat')
      rewind(ntt)
      do m=1,mm
        read(ntt,*)y(m)
      enddo
      close(ntt)
      end
*=======================================================================
      subroutine granal(n,p,grad)
      implicit double precision (a-h,o-z)
      DIMENSION P(N), GRAD(N)
      save
      include 'para'
*     parameter (llw=5,llm=1,lls=1,nnd=1,nn=llw+lls+2+nnd*(llm+1)) 
* mm is the number of the inpout data points
*     parameter (mm=50)
      dimension polm(nnd),x(nnd),gfw(nn)
      dimension xset(mm,nnd),y(mm)
      data itrc/0/
*
      if (itrc.eq.0) then
          itrc=1
          if(n.ne.nn) STOP 'GRANAL:'
*
* Fetch Data:     USER PROVIDED SUBROUTINE:   FETDAT 
*******      
        call fetdat(xset,mm,nnd,y)
*******      
      endif
*
        do ix = 1,n
           grad(ix)= 0.d0
        enddo
      do i=1,mm
        do j=1,nnd
           x(j)=xset(i,j)
        enddo
      call fwann(x,nnd,p,n,llw,llm,lls,polw,polm,pols,fwnn,gfw)

        do ix = 1,n
           grad(ix)= grad(ix) + 2*(fwnn-y(i))*gfw(ix)
        enddo
      enddo
      end

*=======================================================================
      subroutine polynom(p,n,lw,lm,ls,nd,s,polw,polm,pols)
      implicit double precision (a-h,o-z)
      save 
      dimension p(n),polm(nd)
* W(0:LW) coef: p(1:LW+1)
* MU(0:LM,ND) coef: P(LW+2:LW+1+ND*(LM+1))
* SIG(i) coef: P(LW+2+ND*(LM+1):LW+1+ND*(LM+1)+LS+1)
      data itr/0/
      if(itr.eq.0) then
         if(N.ne.LW+LS+2+ND*(LM+1)) then
            write(*,*)'Parameter count inconsistent'
            stop
         endif
         itr=1
      endif
*
      polw = p(lw+1)
      do i=1,lw
         polw=polw*s+p(lw+1-i)
      enddo
*
      pols=p(n)
      do i=1,ls
         pols = pols*s+p(n-i)
      enddo
*
      do j=1,nd
         polm(j) = p(lw+1+j*(lm+1))
         do i=1,lm
            polm(j)=polm(j)*s+p(lw+1+j*(lm+1)-i)
         enddo
      enddo
      end
*=======================================================================
*23456789012345678901234567890123456789012345678901234567890123456789012
*=======================================================================
      subroutine fwann(x,nd,p,n,lw,lm,ls,polw,polm,pols,fwnn,gfw)
      implicit double precision (a-h,o-z)
      parameter (mi=100)
      save
      dimension x(nd),p(n),polm(nd),gfw(n),tse(mi)
      data itrr/0/
      if (itrr.eq.0) then
      itrr=1
      pi=acos(-1.d0)
      do i=1,mi
         tse(i) = cos((2.d0*i-1.d0)/(2.d0*mi)*pi)
      enddo
      endif
*     
      fwnn = 0.
      do i=1,n
        gfw(i)=0.
      enddo
      do mx=1,mi
        s=tse(mx) 
        call polynom(p,n,lw,lm,ls,nd,s,polw,polm,pols)
        anor=0.d0
        do i=1,nd
          anor=anor+(x(i)-polm(i))**2
        enddo
*
        expn=-anor/pols**2/2.d0
        exfac = exp(expn)
        fs=polw*exfac
        fwnn = fwnn + fs
        do i=0,lw
          k=i+1
          gfw(k)=gfw(k)+s**i*exfac
        enddo
        do j=1,nd
          do ia=0,lm
            k=lw+1+ia+(lm+1)*j-lm
            gfw(k)=gfw(k)+fs*s**ia*(x(j)-polm(j))/pols**2
          enddo
        enddo
        do i = 0,ls
           k=lw+1+nd*(lm+1)+1+i
           gfw(k) = gfw(k)+fs*s**i*anor/pols**3
        enddo
      enddo
      fwnn = fwnn*pi/mi
      do i=1,n
        gfw(i)=gfw(i)*pi/mi
      enddo
      end
*=======================================================================
*=======================================================================
      subroutine janal(m,n,p,fj,ld)
      implicit double precision (a-h,o-z)
      DIMENSION P(N),FJ(LD,n) 
      save
      include 'para'
*     parameter (llw=5,llm=1,lls=1,nnd=1,nn=llw+lls+2+nnd*(llm+1)) 
* mm is the number of the inpout data points
*     parameter (mm=50)
      dimension polm(nnd),x(nnd),gfw(nn)
      dimension xset(mm,nnd),y(mm)
      data itrc/0/
*
      if (itrc.eq.0) then
          itrc=1
*
* Fetch Data:     USER PROVIDED SUBROUTINE:   FETDAT 
*******      
        call fetdat(xset,mm,nnd,y)
*******      
      endif
*
      do i=1,mm
        do j=1,nnd
           x(j)=xset(i,j)
        enddo
      call fwann(x,nnd,p,n,llw,llm,lls,polw,polm,pols,fwnn,gfw)
        do ka = 1,n
         fj(i,ka)=gfw(ka)
        enddo
      enddo
      end

