c **** x4fti used in FEVAL  - lx1 set to LT/2 ****

c KLUDGE in "ne funct1" to set x7 x8 ==0

c cxxx has kludge for x5 > 0   - not used
c care ifail=2 in refits - change use 17-02-03
** version that uses fitex.h to set parameter sizes
c****************************************************************
c Start of General correlated fit code here
c****************************************************************
C  FIT T CORRELATED DATA WITH PARAMETERS X
c     written Chris Michael 6-7-1994  (uses NAG - library 18)
C
C  INPUT: DAT(it,ip,nl) storage size itst npst nlst
c      Fits itmin to itmax (Convention 0,..,itmax) for 
c      operators ip=1,..npth. There are nlat data samples (assumed
c      statistically independent)
c
c      Correlated chi**2 fits to DAT using Inverse Correlation  
c      CINV as evaluated by CORFIT. The Fit function (to DAT) is 
c      obtained from FEVAL - given NPAR parameters X 
c      Initial values for X are needed as input
c
c      FIT details - see FEVAL -  
c          This version uses
c          nexp exponentials, npth operators
c          if npar=nexp*(npth+1) then m(nexp) free (idv=1 in TFIT)
c          else if npar=nexp*(npth)-1 then m(nexp) fixed to xdv input
c                                          ( if idv=0 in TFIT)
c          LT needed in FEVAL (for cosh fit - ifexp gives cosh/exp choice)
c          see /factf/ for factorising information in fit

c
c      iprin (print out amount) (none if iprin=0)
c
c      CORRELATION CHOICES: icorrt, icorrp (see CORFIT prologue)
c                           idelt ievals  xlmin
c   
c  OUTPUT:  Final X() parameters, printout, chi2
c
c  CALLS 
c       CORFIT deals with correlations in data 
c          (calls CDISPLAY, INVERC, EIGEN which call F01ABF and F02ABF)
c          (also calls if icorrt=2 COR5, FMIN2, E04FDF, LSFUN1)       
c       
c       NAG minimisation routines (E04JAF)
c          (calls FUNCT1, FEVAL) 
c          Fit function is defined in FEVAL
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE FITEX(dat,itmin,itmax,npth,x,npars,chi2)
      IMPLICIT double precision (A-H,O-Z) 
      Include 'fitex.h'
      PARAMETER (npst=npsth,itst=itsth)

c ** fit workspace - see dispx, feval also for x size

      PARAMETER ( LIW=NPARSv+4, LW=NPARSv*(NPARSv+13) )
      double precision   W(LW),xl(nparsv),xu(nparsv),xx(nparsv)
     &,xorig(nparsv)
      INTEGER IW(LIW)

c***input from calling program*****
      DIMENSION dat(0:itst,npst,nlst),x(npars)
      COMMON/PARMS/ xdv,xlmin,npar,nexp,npthv,lt,
     & icorrt,icorrp,idelt,ievals,itminv,itmaxv,nlat,iprin

      common/factf/nfac1,nfac2,ip1v(npst),ip2v(npst),ifs(npst)  

c***local storage
      COMMON/correl/FD(0:itst,npst),CINV(0:itst,0:itst,npst,npst)
      DIMENSION cdum(0:itst,0:itst,npst,npst)
     & ,c(0:itst,0:itst,npst,npst),fita(0:itst,npst)


c Data input verify
       if(nparsv.lt.np*(nexp+1) ) stop ' nparsv too small'
       if(icorrt.eq.2.and.idelt.le.1) stop 'idelt bad choice'
       if(npth.gt.npst) stop 'fitex storage error: npst'
       if(itmin.lt.0.or.itmax.lt.0) stop 'fitex error: it < 0'
       if(itmin.gt.itst.or.itmax.gt.itst) 
     &      stop 'fitex storage error: itst'
       if(itmin.gt.itmax) stop 'fitex  error: itmax'

      ITminv=ITmin
      ITmaxv=ITmax
      npthv=npth
 
      DO 111 It1=itmin,itmax
      do 111 ip1=1,npth 
        FD(It1,ip1)=0.0
      DO 111 It2=itmin,itmax
      do 111 ip2=1,npth 
       CDUM(It1,It2,ip1,ip2)=0.0
 111  CONTINUE
      
      XNL=1.0/FLOAT(NLAT)
 
      do 100 nl=1,nlat
       DO 70 Ip1=1,npth
        do 71 it1=itmin,itmax
         FD(It1,ip1)=FD(It1,ip1)+DAT(it1,ip1,nl)*XNL
          DO 72 ip2=1,npth
           do 73 it2=itmin,itmax
            CDUM(It1,It2,ip1,ip2)=CDUM(It1,It2,ip1,ip2)
     &      +dat(it1,ip1,nl)*dat(it2,ip2,nl)*XNL
 73        continue     
 72       CONTINUE
 71     continue     
 70    CONTINUE
 
 100  CONTINUE

c assume data set independent in NL (no autocorrelation)
      xf=1.0/float(nlat-1)
      DO 80 Ip1=1,npth
      do 81 it1=itmin,itmax
      DO 82 ip2=1,npth
      do 83 it2=itmin,itmax
      C(It1,It2,ip1,ip2)=(CDUM(It1,It2,ip1,ip2)-
     &    fd(it1,ip1)*fd(it2,ip2))*xf

 83   continue     
 82   CONTINUE
 81   continue     
 80   CONTINUE

C c is data correl matrix 
      call corfit(c,itmin,itmax,itst,
     &     npth,npst,cinv,icorrt,icorrp,idelt,ievals,xlmin,iprin)


C
C       MINIMISATION  FAIRLY CONTINUOUS DERIVATIVES
C
C CARE X(N) GETS OVERWRITTEN IF NPARS>N>NPAR   !
        do ix=1,npar
         xorig(ix)=x(ix)
        enddo
c ifail=0 is OK but ifail=1 suppresses error massages
        IFAIL = 1
C JAF: LIW.GT.N+2  LW.GT.MAX(13,12*N+N(N-1)/2) (KCF less)
C       IBOUND = 0   NEEDS MAX AND MIN FOR ALL X
        IBOUND = 1
cXXX
c        Ibound=0
c        Do ipdx=1,npar
c         xl(ipdx)=-1000.0
c         xu(ipdx)=1000.0
c        Enddo
c        xl(5)=0.0
c        xu(5)=2.0

        ndx=0
        CALL E04JAF( NPAR-ndx,IBOUND, XL,XU,X(1+ndx),
     &   FSUMSQ, IW,LIW,W,LW, IFAIL)

c ifail=1 wrong inputs
c       2 400*n iterations and not converged
c       3 no minimum found - lowest value given
c       4 overflow
c       5-8 doubt if true minimum (5 likely,...8 unlikely)

        IF (IFAIL.GT.0) THEN
         WRITE(6,992) IFAIL,fsumsq,(x(i),i=1,npar)
 992      format(' first attempt,  ifail=',i2,' chi2=',f15.5,
     &    /,15f6.3/,15f6.3/,15f6.3)
         if(iprin.gt.0) call dispx(x,xx,npars)
         ifail=1

        CALL E04JAF( NPAR-ndx,IBOUND, XL,XU,X(1+ndx),
     &   FSUMSQ, IW,LIW,W,LW, IFAIL)


c 17-2-03 swap here to keep orig if ifail=2  was .eq. 2
        if(ifail.eq.4)then
        do ix=1,npar
         x(ix)=xorig(ix)
        enddo
        endif

         if(iprin.gt.0) write(6,933)ifail, fsumsq
 933      format(' redone fit with new ifail=',i2,' chi2=',f15.5,
     &   ' using orig x')
        ENDIF

        chi2=fsumsq
        if(iprin.gt.0)then
C PRINTOUT HERE (care in calling prgm too)

         call dispx(x,xx,npars)

         ndata=(itmax-itmin+1)*npth
         WRITE(6,99) FSUMSQ,ITmin,itmax,npar,ndata,nlat

 99   FORMAT(' FULL CHI2=',F13.3,' IT1=',I2,' it2=',i2,' npars=',i2,
     &' ndata=',i3,' nconfigs=',i4/
     &'  IP1 IP2  IT          FIT        DATA       ERROR        CHI  ')
 
        call feval(fita,itst,Npst,x,npar)

       DO 31 IP=1,npth
        DO 32 IT=ITmin,ITmax
         fit=fita(it,ip)
         ERR=SQRT(C(It,it,ip,ip))
         CHI=(FIT-FD(It,ip))/ERR
         WRITE(6,98) IP1v(ip),IP2v(ip),IT,FIT,FD(It,ip),ERR,CHI
 32     CONTINUE
 31    CONTINUE

 98    FORMAT(2I3,I8,4F13.6)

      endif
      
      RETURN
      END
C******************************************************
c
c written Chris Michael 16-3-94 to 6-7-94
c for general description see C. Michael and A. McKerrell July 1994
c
c input  corr(it1,it2,ip1,ip2)  it=itmin,..itmax ip=1,..npth
c                      
c output cinv(it1,it2,ip1,ip2) according to models:
c NOTE not all elements set to zero (or symmetry cases filled)
c
c flags icorrt icorrp
c         0      0      diagonal inversion
c         1      0      NN model inverse in t, diag in p 
c         1      2      NN model inverse in t, correl in p 
c         0      2      general correl in p, diag in t
c         2      0      5-diag inverse in t, diag  in p
c         2      2      NN model inverse in t, correl in p 
c         3      0      inversion (ievals eigenmodes) in t, diag in p
c         3      2      inversion (ievals eigenmodes) in t, correl in p
c
c       iprin     level of printout 0-8 (use 0 for bootstrap calls)
c         (eigen has iprin 2: evals, 6: evects, 4: display trunc M )
c         (cdisplay has iprin>4 : off diagonal averages)
c   
c       idelt     abs(it1-it2) used for fitting fall-off of corrt=1
c                   for icorrt=2 (5-diag) fit 1 and idelt
c       ievals    number of exact eigenvalues for approximate inverse
c                 remainder of eigenvalues modified, see below:
c                 Use ievals=full size   for exact inverse
c                 If ievals=99  fix minimum eigenvalue to xlmin
c                 If ievals<0 use SVD inverse with |ievals| exact terms
c                 If ievals>0 use ievals  exact eigenmodes and remainder 
c                   same average eigenvalue:  returned as xlmin
c       xlmin     minimum eigenvalue of Correlation matrix (used if 
c                    ievals =99)
c
c CALLS: cdisplay  inverc  eigen
c        which in turn call NAG routines to INVERT and 
c        find eigenvalues of a real symmetric matrix
c
c      (also calls if icorrt=2 COR5, FMIN2, E04FDF, LSFUN1)      
c
c
c SENSIBLE CHOICES
c
c icorrt=0 icorrp=0    is uncorrelated fit
c
c        1        0    fits 1 exponential to correlation in t
c                      for each ip. Use idelt >0 and ~(itmax-itmin)/4
c
c        2        0    fits 2 exponential to corelation in t
c                      for each ip. Use idelt >1 and ~(itmax-itmin)/4
c
c        3        0    uses ievals eigenmodes of correlation in t
c                      exactly for each ip. Use 2< ievals <= itmax-itmin+1
c                      Larger ievals are OK if more data points - 
c                      typically use ievals < sqrt (nlat)                  
c                      Full exact inversion of correlation in t: 
c                        ievals=itmax-itmin+1 
c
c        3        2    uses ievals eigenmodes of correlation between 
c                      data at all t and ip. Use ievals >>2 
c                      <= np*(itmax-itmin+1)   
c                      Larger ievals are OK if more data points - 
c                      typically use ievals < sqrt (nlat)  
c                      Full exact inversion of correlation in t and p: 
c                        ievals=np*(itmax-itmin+1) 
c
C******************************************************
      subroutine corfit(corr,itmin,itmax,itst,
     &     npth,npst,cinv,icorrt,icorrp,idelt,ievals,xlmin,iprin)
        IMPLICIT double precision (A-H,O-Z)

        Include 'fitex.h'
        parameter (npthl=npsth,itstl=itsth*npthl)
       
       dimension corr(0:itst,0:itst,npst,npst),
     &           cinv(0:itst,0:itst,npst,npst),
     &           cnorm(0:itstl),cdum(0:itstl,0:itstl),
     &           xas(npthl),xrr(npthl,npthl),xl(0:itstl)

c Data input verify
       if(itstl.lt.itmax.or.npthl.lt.npth)
     &                     stop 'corfit storage error'

	itval=itmax-itmin  
        itvalp=itval+1
         itstp=itstl+1

       if(iprin.ge.4)then
C DISPLAY input correlation blocks
           do 50 ip1=1,npth
           do 51 ip2=ip1,npth
            do 56 it1=itmin,itmax
            do 57 it2=itmin,itmax
              cdum(it1,it2)=corr(it1,it2,ip1,ip2)
     &        /sqrt( corr(it1,it1,ip1,ip1)* corr(it2,it2,ip2,ip2) )
 57         continue
 56         continue
            write(6,910)ip1,ip2
 910     format(/' input normalised correlation t*t for paths=',2i2)

            call cdisplay(cdum,itmin,itmax,itstl,5)

c evaluate estimate of smallest eigenvalue from alternating eigenmode

           v=0.0
           s1=-1.0
           do 65 it1=itmin,itmax
            s1=-s1
            s2=-1.0
            do 66 it2=itmin,itmax
             s2=-s2
             v=v+s1*cdum(it1,it2)*s2
 66         continue
 65        continue
c             normalise vector +-+-+-,..
           v=v/itvalp
           write(6,971)v
 971   format(/' Estimates of smallest eigenvalues by various methods'
     & /' projection of norm correl on vector +-+-,.. is',f14.7)

c evaluate estimates of smallest eigenvector from submatrices
           if(ip1.eq.ip2.and.icorrt.ge.2)then
c maximum submatrix given by statistics: ievals?
           idtmax=min(itval,ievals)
           do 68 itd=1,idtmax
           itdiff=itd+1
           itb=itmax-itd
           do 67 ita=itmin,itb

c use cinv as dummy to receive output not needed
c use iprin=-1 specially to return least eigenvalue 
c        (no output to cinv)
           call eigen(cdum(ita,ita),itstp,itdiff,1,xl(ita),cinv,-1)

 67        continue
           write(6,939)itdiff,itdiff
 939       format(i2,' by ',i2,' submatrix has min evals:')

            write(6,9)(xl(it),it=itmin,itb)
 9          format(5f13.6)
 68        continue
       endif

 51    continue
 50    continue
       endif

       if(icorrp.eq.0.or.npth.eq.1)then
c**************BLOCK DIAGONAL IN PATHS ****************************
        do 1 ip=1,npth
        if(icorrt.eq.0) then
c---diagonal in t ------------------------------------
         do 20 it=itmin,itmax
          cinv(it,it,ip,ip)=1.0/corr(it,it,ip,ip)
 20      continue

        else

c       measure corr normalised by diagonal: model it
c       icorrt.ne.0 ------------------------------------

        do 30 it=itmin,itmax
	cnorm(it)=1.0/sqrt(corr(it,it,ip,ip))
 30     continue

calculate average normalised correlation at idelt  

        if(idelt+itmin.gt.itmax)stop ' idelt too big'
        nrav=0
	rav=0.0 
        nrav1=0
	rav1=0.0 
        do 301 it1=itmin,itmax
         it2=it1+idelt
           if(it2.le.itmax)then
	     nrav=nrav+1
             rav=rav+corr(it1,it2,ip,ip)*cnorm(it1)*cnorm(it2)
           endif
         it2=it1+1
           if(it2.le.itmax)then
	     nrav1=nrav1+1
             rav1=rav1+corr(it1,it2,ip,ip)*cnorm(it1)*cnorm(it2)
           endif
 301    continue
        rav=rav/float(nrav)
        rav1=rav1/float(nrav1)
   
        if(iprin.ge.2)then
        write(6,930)ip,1,rav1,idelt,rav
 930    format(' norm correl ip=',i2,' average  C(t1-t2=',i2,')=',f8.5
     &  /'                   average  C(t1-t2=',i2,')=',f8.5)
        endif
         

        if(icorrt.eq.1)then
c---nearest neighbour in cinv in t-----------
c... option to guess NN model to try and fit C


        ravv=(rav)**(1.0/float(idelt))
       xdiag0=1.0/(1.0-ravv*ravv)
       xdiag=(1.0+ravv*ravv)*xdiag0
       xnn=-ravv*xdiag0

c  speed up by just diagonal and off diagonal initialised (see FUNCT1)

       do 430 it1=itmin,itmax-1
        it2=it1+1
      cinv(it1,it1,ip,ip)=xdiag*(cnorm(it1)*cnorm(it1))
      cinv(it1,it2,ip,ip)=xnn*(cnorm(it1)*cnorm(it2))
      cinv(it2,it1,ip,ip)=xnn*(cnorm(it1)*cnorm(it2))

 430   continue
      cinv(itmin,itmin,ip,ip)=xdiag0*(cnorm(itmin)*cnorm(itmin))
      cinv(itmax,itmax,ip,ip)=xdiag0*(cnorm(itmax)*cnorm(itmax))


         if(iprin.ge.4) then

         do 410 it1=itmin,itmax
         do 411 it2=itmin,itmax
          cdum(it1,it2)=0.0
          if(abs(it1-it2).eq.1) cdum(it1,it2)=xnn
          if(it1.eq.it2) cdum(it1,it2)=xdiag      
          if(it1.eq.itmax.and.it2.eq.itmax)cdum(it1,it2)=xdiag0
          if(it1.eq.itmin.and.it2.eq.itmin)cdum(it1,it2)=xdiag0
          cinv(it1,it2,ip,ip)=cdum(it1,it2)*(cnorm(it1)*cnorm(it2))
 411     continue
 410     continue

          if(iprin.ge.8) then 
             call cdisplay(cdum,itmin,itmax,itstl,iprin)
c to get more info on eigenvectors:
             call eigen(cdum(itmin,itmin),itstp,itvalp,
     &           itvalp,xlmin,cdum(itmin,itmin),2)
          else
             call inverc(cdum(itmin,itmin),cdum(itmin,itmin),
     &            itstp,itvalp)
          endif
            write(6,911)ip
 911     format(/'  model normalised correlation t*t for path=',i2)
              call cdisplay(cdum,itmin,itmax,itstl,5)
         endif

       endif
	
        if(icorrt.eq.2)then

      call cor5(rav1,1,rav,idelt,dd0,xn0,dd1,dd,xn,xnn,iprin)


c  speed up by just diagonal and 4 off diagonals initialised (see FUNCT1)
c  could speed up more by just upper diagonals used (symmetric)

        if(iprin.ge.4)then 
         do 510 it1=itmin,itmax
         do 511 it2=itmin,itmax
          cinv(it1,it2,ip,ip)=0.0
 511     continue
 510     continue
        endif

       if(itmax-itmin.lt.4) stop ' 5-diagonal wont fit'

      cinv(itmin,itmin,ip,ip)=dd0*(cnorm(itmin)*cnorm(itmin))
      cinv(itmin+1,itmin+1,ip,ip)=dd1*(cnorm(itmin+1)*cnorm(itmin+1))
      do 500 it1=itmin+2,itmax-2
      cinv(it1,it1,ip,ip)=dd*(cnorm(it1)*cnorm(it1))
 500   continue
      cinv(itmax-1,itmax-1,ip,ip)=dd1*(cnorm(itmax-1)*cnorm(itmax-1))
      cinv(itmax,itmax,ip,ip)=dd0*(cnorm(itmax)*cnorm(itmax))

      cinv(itmin,itmin+1,ip,ip)=xn0*(cnorm(itmin)*cnorm(itmin+1))
      cinv(itmin+1,itmin,ip,ip)=xn0*(cnorm(itmin)*cnorm(itmin+1))
       do 530 it1=itmin+1,itmax-2
        it2=it1+1
      cinv(it1,it2,ip,ip)=xn*(cnorm(it1)*cnorm(it2))
      cinv(it2,it1,ip,ip)=xn*(cnorm(it1)*cnorm(it2))
 530   continue
      cinv(itmax,itmax-1,ip,ip)=xn0*(cnorm(itmax)*cnorm(itmax-1))
      cinv(itmax-1,itmax,ip,ip)=xn0*(cnorm(itmax)*cnorm(itmax-1))

       do 501 it1=itmin,itmax-2
        it2=it1+2
      cinv(it1,it2,ip,ip)=xnn*(cnorm(it1)*cnorm(it2))
      cinv(it2,it1,ip,ip)=xnn*(cnorm(it1)*cnorm(it2))
 501   continue


      
        if(iprin.ge.4)then 

         do 520 it1=itmin,itmax
         do 521 it2=itmin,itmax
          cdum(it1,it2)=cinv(it1,it2,ip,ip)/(cnorm(it1)*cnorm(it2))
 521     continue
 520     continue


          if(iprin.ge.8) then 
           write(6,913)ip
 913      format(/'  model norm inverse correlation t*t for path=',i2)
          call cdisplay(cdum,itmin,itmax,itstl,iprin)  
c to get more info on eigenvectors:
             call eigen(cdum(itmin,itmin),itstp,itvalp,
     &           itvalp,xlmin,cdum(itmin,itmin),2)
          else
             call inverc(cdum(itmin,itmin),cdum(itmin,itmin),
     &            itstp,itvalp)
          endif
           write(6,912)ip
 912     format(/'  model normalised correlation t*t for path=',i2)
           call cdisplay(cdum,itmin,itmax,itstl,5)
         
        endif

	endif
c-----expand Corr in eigenvalues, and modify if needed---------------
        if(icorrt.eq.3)then

       do 610 it1=itmin,itmax
       do 611 it2=itmin,itmax
        cdum(it1,it2)=corr(it1,it2,ip,ip)*cnorm(it1)*cnorm(it2)
 611   continue
 610   continue

c inverse (truncated if ievals relevant)
             call eigen(cdum(itmin,itmin),itstp,itvalp,
     &           ievals,xlmin,cdum(itmin,itmin),min(iprin,2))


       do 620 it1=itmin,itmax
       do 621 it2=itmin,itmax
        cinv(it1,it2,ip,ip)=cdum(it1,it2)*(cnorm(it1)*cnorm(it2))
 621   continue
 620   continue

        if(iprin.ge.8)then
          write(6,915)ip
 915     format(/' inverse normalised correlation t*t for path=',i2)
          call cdisplay(cdum,itmin,itmax,itstl,5)  
        endif

        if(iprin.ge.4) then
         write(6,919)ip
 919     format(/' model normalised correlation t*t for path=',i2)
       call inverc(cdum(itmin,itmin),cdum(itmin,itmin),itstp,itvalp)
          call cdisplay(cdum,itmin,itmax,itstl,5)  
        endif
	 

        endif	
        endif
1       continue
        RETURN

        else if(icorrp.ge.2.and.npth.gt.1)then
c**************OFF DIAGONAL SYMMETRIC IN PATHS ********************
        if(icorrt.eq.0)then
         if(iprin.ge.8)    write(6,977)
 977        format(/'   it    inverse norm corr(ip1,ip2)  ,....')
         do 700 it=itmin,itmax
          do 702 ip1=1,npth
          do 704 ip2=1,npth
           cdum(ip1,ip2)=corr(it,it,ip1,ip2)
 704      continue
 702      continue
         call inverc(cdum(1,1),cdum(1,1),itstp,npth)
          do 712 ip1=1,npth
          do 714 ip2=1,npth
           cinv(it,it,ip1,ip2)=cdum(ip1,ip2)
           if(iprin.ge.8) cdum(ip1,ip2)=cdum(ip1,ip2)*
     &     sqrt(corr(it,it,ip1,ip1)*corr(it,it,ip2,ip2))
 714      continue
 712      continue

          if(iprin.ge.8)then
            write(6,978)it, ((cdum(i,j),j=1,npth),i=1,npth)
 978        format(i4,3(6f12.5/))
          endif

 700     continue

         return
         endif

c hereon ip and it correlated:

         if(itstl.lt.npth*itvalp) stop 'corfit: storage 2,2 problem'

          do 732 ip1=1,npth
          do 730 it1=itmin,itmax
           iv1=(ip1-1)*itvalp+it1-itmin
           cnorm(iv1)=1.0/sqrt(corr(it1,it1,ip1,ip1))
 730      continue
 732      continue

         if(icorrt.eq.1.or.icorrt.eq.2)then
c model diagonal and off diagonal correlations
calculate average normalised correlation at idelt  (and at 0 )
c positivity of model needed: one way is to fix slope as shown

        do 735 ip1=1,npth
        do 734 ip2=ip1,npth
        nrav=0
	rav=0.0 
        rrr=0.0
        do 736 it1=itmin,itmax
           iv1=(ip1-1)*itvalp+it1-itmin
           iw2=(ip2-1)*itvalp+it1-itmin
           if(ip1.ne.ip2)
     &            rrr=rrr+corr(it1,it1,ip1,ip2)*cnorm(iv1)*cnorm(iw2)
         it2=it1+idelt
           iv2=(ip2-1)*itvalp+it2-itmin
           if(it2.le.itmax)then
	     nrav=nrav+1
             rav=rav+corr(it1,it2,ip1,ip2)*cnorm(iv1)*cnorm(iv2)
           endif
         it2=it1-idelt
           iv2=(ip2-1)*itvalp+it2-itmin
           if(it2.ge.itmin)then
	     nrav=nrav+1
             rav=rav+corr(it1,it2,ip1,ip2)*cnorm(iv1)*cnorm(iv2)
           endif
 736    continue
        rav=rav/float(nrav)
        if(ip1.ne.ip2) then
         xrr(ip1,ip2)=rrr/itvalp
        else
c****icorrt=1 ****here*****
         xas(ip1)=-log(rav)/float(idelt)
         xrr(ip1,ip2)=1.0
        endif
        rav=rav/xrr(ip1,ip2)
        if(iprin.ge.2)then
          write(6,931)ip1,ip2,idelt,rav,xrr(ip1,ip2)
 931      format(' norm correl ip1,ip2=',2i2,' C(',i2,
     &    ')/C(0)=',f8.5,'  C(0)=',f8.5)
        endif

 734    continue
 735    continue

c model correlation here
         do 748 ip1=1,npth
         do 749 ip2=ip1,npth
         xa=0.5*(xas(ip1)+xas(ip2))

   
        do 740 it1=itmin,itmax
           iv1=(ip1-1)*itvalp+it1-itmin
           ivx1=(ip2-1)*itvalp+it1-itmin
        do 741 it2=itmin,itmax
           iv2=(ip2-1)*itvalp+it2-itmin
           ivx2=(ip1-1)*itvalp+it2-itmin

         cdum(iv1,iv2)=xrr(ip1,ip2)*exp(-xa*abs(it2-it1))
c assume symmetric in ip1, ip2
                        cdum(ivx1,ivx2)=cdum(iv1,iv2)
 741    continue
 740    continue       
 749    continue
 748    continue

c cdum set up icort = 1 or 2
        endif

        if (icorrt.eq.3)then
c**************FULL MATRIX  IN TIME AND PATHS **************************

         do 800 it1=itmin,itmax
         do 801 it2=itmin,itmax
          do 802 ip1=1,npth
           iv1=(ip1-1)*itvalp+it1-itmin
          do 804 ip2=1,npth
           iv2=(ip2-1)*itvalp+it2-itmin
           cdum(iv1,iv2)=corr(it1,it2,ip1,ip2)*cnorm(iv1)*cnorm(iv2)
 804      continue
 802      continue
 801     continue
 800     continue

         endif

c invert cdum for icorrt = 1,2 or 3
c need full inversion even for icorrt=1 because of off diag in ip
care negative eigenvalues possible if model is dud
c       call inverc(cdum,cdum,itstp,npth*itvalp)
        iev=ievals
        if(icorrt.eq.1) iev=npth*itvalp
        call eigen(cdum,itstp,npth*itvalp,iev,
     &                xlmin,cdum,min(2,iprin))
        do 810 it1=itmin,itmax
         do 811 it2=itmin,itmax
          do 812 ip1=1,npth
           iv1=(ip1-1)*itvalp+it1-itmin
          do 814 ip2=1,npth
           iv2=(ip2-1)*itvalp+it2-itmin
           cinv(it1,it2,ip1,ip2)=cdum(iv1,iv2)*cnorm(iv1)*cnorm(iv2)
 814      continue
 812      continue
 811     continue
 810    continue


c display inverse model correlation (normalised)
          if(iprin.ge.8)then
         do 870 ip1=1,npth
         do 871 ip2=ip1,npth

           write(6,985)ip1,ip2
 985       format(/' inverse norm correl for ip=',i2,',',i2)
           call cdisplay(cdum((ip1-1)*itvalp,(ip2-1)*itvalp)
     &     ,0,itval,itstl,iprin) 

 871     continue
 870     continue
          endif
   
c display truncated input - reinvert cdum (not needed icorrt=1)
         if(iprin.ge.4.and.iev.ne.npth*itvalp)then

          call inverc(cdum,cdum,itstp,npth*itvalp)

         do 873 ip1=1,npth
         do 874 ip2=ip1,npth

          write(6,988)ip1,ip2
 988       format(/' model norm correl for ip=',i2,',',i2)
          call cdisplay(cdum((ip1-1)*itvalp,(ip2-1)*itvalp)
     &    ,0,itval,itstl,iprin)

 874     continue
 873     continue

         endif
         RETURN


        else
        stop ' corfit: wrong icorrt icorrp choices'

        endif
        return
        end
c******************************************************************
c  CALLED by CORFIT
c Code to find sensible 5-diagonal matrix inverse 
c    to given (reduced) correlations
c    C(idt=it1)=p1, C(idt=it2)=p2  (assumes C(0)=1.0 )
c
c  CALLS FMIN2
c
c 5-diagonal convention:
c  [x y c 0 0 0 0
c   y z b c 0 0 0
c   c b a b c 0 0
c   0 c b a b c 0
c   0 0 c b a b c
c   0 0 0 c b z y
c   0 0 0 0 c y x]
c
c******************************************************************
      subroutine cor5(p1in,it1in,p2in,it2in,x,y,z,a,b,c,iprin)
      implicit double precision (a-h,o-z)
      common/abfit/p1,p2,it1,it2

      if(p1.gt.1.0.or.p2.gt.1.0) stop 'cor5: corr > 1'
      p1=p1in
      p2=p2in
      it1=it1in
      it2=it2in
c assume 1 exponential to start
      xl1=(p2)**(1.0/float(it2))
      call fmin2(xl1,0.0,xf,yf,iprin)


      xl1=xf
      xl2=yf

      if(abs(xl1-xl2).lt.1.0D-10) xl2=xl1-1.0D-6
      do 10 ir=1,19
      ctw=( xl1**(ir+1)*(1.0-xl2**2) -xl2**(ir+1)*(1.0-xl1**2) )
     &  /( (xl1-xl2)*(1.0+xl1*xl2) )
      if(ir.eq.it1)c1v=ctw
      if(ir.eq.it2)c2v=ctw
      if(ir.eq.1)p1v=ctw
      if(ir.eq.2)p2v=ctw
      if(iprin.ge.8) write(6,900)ir,ctw
 900  format(' idt=',i2,'  corr~ =',f12.6)
 10   continue

      xd=1.0/( (xl1-xl2)*(1.0+xl1*xl2)  )
      xc1=xd*(1.0-xl2**2)
      xc2=xd*(1.0-xl1**2)

      FFF=(xc1*xl1**(it1+1)-xc2*xl2**(it1+1)-p1)**2

      FFF=FFF+(xc1*xl1**(it2+1)-xc2*xl2**(it2+1)-p2)**2

      if(iprin.ge.4) write(6,92)xl1,xl2,FFF,it1,c1v,it2,c2v
 92   format(' Eigenvalues=',2g14.6,' chi^2 fit=',g15.7 /
     &  ' with c(',i2,')=',f12.6,' and c(',i2,')=',f12.6)

      if(abs(min(xl1,xl2)).le.1.0D-10) then
       xl=max(xl1,xl2)
       xx1=xl+1.0/xl
       bbya=-1.0/xx1
       a=1.0/(1.0+p1v*2.0*bbya)
       b=bbya*a
       c=0.0
       x=a/(xx1*xl)
       z=a
       y=b
      else

       xx1=xl1+1.0/xl1
       xx2=xl2+1.0/xl2
       bbyc=-xx1-xx2
       cbya=1.0/(2.0+xx1*xx2)
       bbya=bbyc*cbya
       a=1.0/(1.0+p1v*2.0*bbya+p2v*2.0*cbya)
       b=bbya*a
       c=cbya*a
       x=c/(xl1*xl2)
       z=a-c*xl1*xl2
       y=-c*(xl1+xl2)/(xl1*xl2)
      endif
      if(iprin.ge.8) write(6,1)p1v,p2v,x,y,z,a,b,c
    1 format(2f12.5,3x,3f15.5/3f15.5)
      return
      end
c****************************************************************
c  CALLED by COR5
c NAG minimising code F(x,y) start x0,y0  final xf,yf
c LEAST SQUARES ( tested 30-6-94 )
c
c THIS VERSION FIXES 0.0 <= x,y <= 1.0
c by mapping xl= tanh(x)*.5+.5
c  CALLS LSFUN1, E04FDF
c****************************************************************
      subroutine fmin2(x0,y0,xf,yf,iprin)
      implicit double precision (a-h,o-z)
      external LSFUN1
      dimension xv(2),W(50),iw(50)

c      xv(1)=x0
c      xv(2)=y0
       xv(1)=0.5 * log ( x0 /(1.0-x0) )
       xv(2)=-1.0
      ifail=1
      call e04fdf(2,2,xv,fsumsq,iw,50,w,50,ifail)
      if(ifail.gt.1) write(6,99) ifail, fsumsq
 99   format('cor5: icorrt=2 fit  ifail=',i8,' chi2=',g15.6)
      xf=tanh(xv(1))*0.5+0.5
      yf=tanh(xv(2))*0.5+0.5
      return
      end
c******************************************************************
c called by e04fdf
c i=1,m  Sum_of_f(i)**2 fitted  parameters xm(n)
c******************************************************************
      subroutine lsfun1(m,n,xm,f)
      implicit double precision (a-h,o-z)
      dimension xm(n),f(m),x(4)
      common/abfit/p1,p2,it1,it2

      x(1)=tanh(xm(1))*0.5+0.5
      x(2)=tanh(xm(2))*0.5+0.5

      xd=1.0/( (x(1)-x(2))*(1.0+x(1)*x(2))  )
      xc1=xd*(1.0-x(2)**2)
      xc2=xd*(1.0-x(1)**2)

      F(1)=xc1*x(1)**(it1+1)-xc2*x(2)**(it1+1)-p1

      F(2)=xc1*x(1)**(it2+1)-xc2*x(2)**(it2+1)-p2
      return
      end
c******************************************************************
c  CALLED by CORFIT
c  display c(it1,it2) it=itmin,..itmax and calculate average  
c   off axis diagonal values (if iprin > 4)
c    c normalised for display by scale XSC
c
c
c******************************************************************
        subroutine cdisplay(c,itmin,itmax,itst,iprin)
        IMPLICIT double precision (A-H,O-Z)
        Include 'fitex.h'
        parameter (itstv=npsth*itsth)
        dimension c(0:itst,0:itst),irv(itstv),inum(0:itstv),rav(0:itstv)

        if(itmax.gt.itstv)stop 'cdisplay: storage error'

        nit=itmax-itmin
        do 1 it=0,itmax
	 inum(it)=it
 1      continue
c fix SCALE to display as I3 format
        xsc=100.0
        if(abs(c(itmin,itmin)).gt.1.01) xsc=10.0
        if(abs(c(itmin,itmin)).gt.10.01) xsc=1.0
        if(abs(c(itmin,itmin)).gt.100.1) xsc=0.1
        if(abs(c(itmin,itmin)).gt.1000.1) xsc=0.01
        write(6,991)itmin,itmax,itmin,itmax,xsc
991      format(' cdisplay: c(',i2,':',i2,',',i2,':',i2,') with scale'
     &   ,' xsc=',f8.3) 
        write(6,993) itst,(inum(it),it=itmin,itmax)
 993      format(i4,':',25i3)
        write(6,994)
 994      format(' ')
	do 2 idt=0,nit
         rav(idt)=0.0
 2      continue

 
        do 31 it1=itmin,itmax
         do 32 it2=itmin,itmax

	  idt=abs(it2-it1)
          rav(idt)=rav(idt)+c(it1,it2)
          irv(it2)=nint(xsc*c(it1,it2))
 32      continue

         write(6,99) it1,(irv(it2),it2=itmin,itmax)
 99      format(i4,':',25i3)
 31     continue
 
        if(iprin.gt.4)then

         write(6,907)
 907      format(' Average off-diagonal elements')

        do 3 idt=1,nit
         rav(idt)=0.5*rav(idt)/(nit-idt+1)
 3      continue
         rav(0)=rav(0)/(nit+1)

        write(6,902)(inum(it),rav(it),it=0,nit)
 902     format(5(i3,':',f11.5))
        endif
        return
        end 
C======================================================================C
C                                                                      C
      SUBROUTINE FUNCT1( NPAR , X, FC )
C                                                                      C
C     LEAST SQUARES FUNCTION CALLED BY E04JAF
C     Calls Feval to get F(X); 
C     needs data FD and correlation CINV 
c     needs flags icorrt icorrp
C     CALLS FEVAL
C======================================================================C
C
      IMPLICIT double precision (A-H,O-Z)
      Include 'fitex.h'
      parameter (itst=itsth,npst=npsth)

c***input from calling program*****
      DIMENSION  X(NPAR) 
      COMMON/correl/FD(0:itst,npst),CINV(0:itst,0:itst,npst,npst)
      COMMON/PARMS/xdv,xlmin,nparv,nexp,npth,lt,
     & icorrt,icorrp,idelt,ievals,itmin,itmax,nlat,iprin

c**local***
      Logical ifp
      Dimension  FIT(0:itst,npst)

      ifp=icorrp.ne.0.and.npth.gt.1

      call Feval(fit,itst,npst,X,npar)

      FC = 0.

c     g4 gi   gi  g5 gi  - want a1 zero to g4 gi   - say 7,8
c KLUDGE here to force a parameter to fixed value (assumed zero here)
      if(ikludge .eq. 123) FC=100000.0*( x(7)**2+x(8)**2 )
      if(ikludge .eq. 246) FC=100000.0*( x(5)**2+x(6)**2 )
      if(ikludge .eq. 987) FC=100000.0*( x(9)**2+x(10)**2 )

       npthm=npth-1

c       itminl=9       itminll=10  from fitex.h

      do 21 ip1=1,npth  
       DO 20 it1 = itmin,itmax
         ifdo=1
        if(ikludge .eq. 666 .and. ip1 .gt. 10 .and. it1 .lt. itminl)
     &   ifdo=0  
        if(ikludge .eq. 666 .and. (ip1 .gt. 19 .or. ip1 .eq. 15) 
     &   .and. it1 .lt. itminll)      ifdo=0  
         if (ifdo .eq. 1) 
     &  fc=fc+cinv(it1,it1,ip1,ip1)*(fit(it1,ip1)-fd(it1,ip1))**2
 20    continue
 21   continue      
 
      if(.not.ifp)then  

      if(icorrt.eq.0) return

c---Nearest neighbour in t (assume symmetric)------

      if(icorrt.ge.1) then
      do 31 ip1=1,npth  
       DO 30 it1 = itmin,itmax-1
         it2=it1+1
          fc=fc+cinv(it1,it2,ip1,ip1)*(fit(it1,ip1)-fd(it1,ip1))*
     &                            2.0*(fit(it2,ip1)-fd(it2,ip1))
 30    continue
 31   continue      
      endif

      if(icorrt.ge.2)then
c also add further diagonal in t (assume symmetric)

      do 300 ip1=1,npth
       DO 303 it1 = itmin,itmax-2
        it2=it1+2
          fc=fc+cinv(it1,it2,ip1,ip1)*(fit(it1,ip1)-fd(it1,ip1))*
     &                        2.0*    (fit(it2,ip1)-fd(it2,ip1))
 303    continue
 300   continue   
      endif
      
c---General correlation in t (assume symmetric)------

      if(icorrt.gt.2.and.(itmax-itmin).gt.2) then
      do 41 ip1=1,npth  
       DO 40 it1 = itmin,itmax-3
       DO 42 it2 = it1+3,itmax
          fc=fc+cinv(it1,it2,ip1,ip1)*(fit(it1,ip1)-fd(it1,ip1))*
     &                        2.0*    (fit(it2,ip1)-fd(it2,ip1))
 42    continue
 40    continue
 41   continue      
      endif

      else
c+++++++several paths, correlated here ++++++
c add to diagonal correlation term first

      if(icorrt.eq.0)then
c offdiagonal in p
       do 50 ip1=1,npthm
        do 51 ip2=ip1+1,npth
         do 55 it=itmin,itmax
          fc=fc+cinv(it,it,ip1,ip2)*(fit(it,ip1)-fd(it,ip1))*
     &                          2.0*(fit(it,ip2)-fd(it,ip2))
 55     continue
 51     continue
 50    continue

      else

c offdiagonal in t 
      do 471 ip1=1,npth  
       DO 470 it1 = itmin,itmax-1
       DO 472 it2 = it1+1,itmax
          fc=fc+cinv(it1,it2,ip1,ip1)*(fit(it1,ip1)-fd(it1,ip1))*
     &                        2.0*    (fit(it2,ip1)-fd(it2,ip1))
 472    continue
 470    continue
 471   continue   
   
c offdiagonal in p
c assume symmetric under ip1 ip2 swap not it1 it2
      do 400 ip1=1,npthm
      DO 401 ip2 =ip1+1,npth
       DO 402 it1 = itmin,itmax
       do 403 it2 = itmin,itmax
          fc=fc+cinv(it1,it2,ip1,ip2)*(fit(it1,ip1)-fd(it1,ip1))*
     &                        2.0*    (fit(it2,ip2)-fd(it2,ip2))
 403    continue
 402    continue
 401   continue      
 400   continue   

       endif
      endif


      if(FC.lt.-1.0D-8) STOP ' Funct1: chi**2 is negative'      

      RETURN
      END

c******************************************************************
c  CALLED by FUNCT1, FITEX, TFIT, GRAPH
c  Function for fits: F(it,ip)  it convention 0,..
c   given NPAR parameters X  
c   as written    NEXP exponential fit to NPTH operators ip
c   options: cosh (needs lt ) ifexp=0
c            exp  ifexp=-1 (Could speed up this case in code)
c            sinh (needs lt) ifexp=1
c          m(nexp) fixed at xdv if npar < nexp*(npth+1)
c
c Conventions   if factorising:
c    F=x1p*x2p exp(-xm1*t)+ x1p*x2p exp(-(xm1+xm2)*t)
c else (nfac1=1)
c    F=x2p exp(-xm1*t)+ x2p exp(-(xm1+xm2)*t)
c
c   ORDER x path1 m0;  path2 m0 ,  etc
c   normalise coefficients to 1.0 at t=tmin for ip=1 (by XD in TFIT)
c
c  changes here need initial values set appropriately (TFIT FITEX)
c     and display modified (DISPX)
c npar set by funct1 if called E04
c******************************************************************
        subroutine feval(f,itst,npst,xv,npar)
        IMPLICIT double precision (A-H,O-Z)
        Include 'fitex.h'
        PARAMETER (npstv=npsth)
        dimension xm(nparsv),f(0:itst,npst),xv(npar),x(nparsv)

c***input from calling program*****
      COMMON/PARMS/ xdv,xlmin,nparv,nexp,npth,lt,
     & icorrt,icorrp,idelt,ievals,itmin,itmax,nlat,iprin


       common/factf/nfac1,nfac2,ip1v(npstv),ip2v(npstv),ifs(npstv)


        np=max(nfac1,nfac2)
        nm=nexp*np
        nmm=nm+nexp

         do ix=1,npar
          x(ix)=xv(ix)
         enddo


        xm(1)=x(nm+1)

       if( ifs(1).eq. 88) then
       do ip = 1,np
	do  it=itmin,itmax
         f(it,ip)=xm(1)+x(ip)*(float(it-lt/2))**2
        enddo
       enddo
       return
       endif
	
	do 10 it=itmin,itmax
	 ca=exp(-xm(1)*(it-itmin))
         cb= exp(-xm(1)*(lt-it-itmin))
	 cc=ca+cb
         cs=ca-cb
         do 11 ip=1,npth

         if(ifs(ip).eq.0)then
            ce=cc
          elseif(ifs(ip).eq.1)then
            ce=cs
          elseif(ifs(ip).eq.-1)then
            ce=ca
          elseif(ifs(ip).eq.9)then
c point-point  (code in tmfit24 - so care if not linked)
          irv=2
c          lx1=24
           lx1=lt/2
c use interpolation here now (need lx1=lt/2) -- code in tmfit48 only
c BUT STUB HERE ALSO FOR COMPILER
          ce=x4fti(xm(1),it,irv,lx1,lt)/x4fti(xm(1),itmin,irv,lx1,lt)
          endif
           if(nfac1.eq.1)then
          f(it,ip)=x(ip2v(ip))*ce
           else
          f(it,ip)=x(ip1v(ip))*x(ip2v(ip))*ce
           endif
 11      continue
 10     continue
 
        if(nexp.gt.1)then
	
        do 20 im=2,nexp
	  xm(im)=x(nm+im)
 20     continue
        if(npar.eq.nmm-1) xm(nexp)=xdv
	do 30 it=itmin,itmax
	 do 33 im=2,nexp
	 ca=exp(-xm(im)*(it-itmin))
         cb= exp(-xm(im)*(lt-it-itmin))
	 cc=ca+cb
         cs=ca-cb
         do 31 ip=1,npth

         if(ifs(ip).eq.0)then
            ce=cc
          elseif(ifs(ip).eq.1)then
            ce=cs
          elseif(ifs(ip).eq.-1)then
            ce=ca
          elseif(ifs(ip).eq.9)then
c point-point  (code in tmfit24 - so care if not linked)
          irv=2
c          lx1=24
           lx1=lt/2
          ce=x4fti(xm(im),it,irv,lx1,lt)/x4fti(xm(im),itmin,irv,lx1,lt)
          endif
           if(nfac1.eq.1)then
          f(it,ip)=f(it,ip)+x(np*(im-1)+ip2v(ip))*ce
           else
          f(it,ip)=f(it,ip)+x(np*(im-1)+ip1v(ip))*
     &                      x(np*(im-1)+ip2v(ip))*ce
           endif
 31      continue
 33      continue
 30     continue
       
        endif 

        Return
        end
c--------------------------------------------------------------------
c  CALLED by CORFIT 
c  invert C (in place allowed)  dimension ndstv size ndata
c  assumes C is positive definite and symmetric
c  cdum has only i>=j elements with ci(i,j) filled
c  CALLS F01ABF
c--------------------------------------------------------------------
      subroutine inverc(c,ci,ndstv,ndata)
      IMPLICIT double precision (A-H,O-Z)
      Include 'fitex.h'
      PARAMETER (npst=npsth,itst=itsth,NDST=NPST*ITST)
      dimension c(ndstv,ndata),ci(ndstv,ndata)
      double precision WKSPCE(NDST),cdum(ndst,ndst)
      if(ndst.lt.ndata.or.ndstv.lt.ndata) stop ' inverc: storage error' 
      IFAIL=0
      CALL F01ABF(c,NDSTv,NDATA,cdum,NDST,WKSPCE,IFAIL)
      if(ifail.eq.1) write(6,98)ifail
 98   format(' inversion in inverc failed ifail=',i2,' not positive 
     & definite')
      if(ifail.eq.2) write(6,99)ifail
 99   format(' inversion in inverc failed ifail=',i2,' ill conditioned')
      do 1 i=1,ndata
       do 2 j=1,i
        ci(i,j)=cdum(i,j)
        ci(j,i)=cdum(i,j)
 2     continue
 1    continue
      RETURN
      END
C***********************************************************************
c   CALLED by CORFIT
C   CALCULATE EIGENVECTORS  and hence INVERSE
c If ieval=99  fix minimum eigenvalue to xlmin
c If ieval<0 use SVD inverse keping |ieval| exact terms
c If ieval>0 use ieval  exact eigenmodes and remainder 
c  same average eigenvalue:  returned as xlmin
c
c output cinv (with modified eigenmodes) can be in place
c
c   NOTE iprin=-1 returns least eigenvalue as xlmin
c        iprin>0  negative eigenvalue warning and inversion summary
c        iprin>1  eigenvalues and truncation info
c        iprin>3  display (truncated) input matrix
c        iprin>7  eigenvectors
c
c MODS 2-3-94 to "smooth" xl more
c  CALLS cdisplay, f02abf
C***********************************************************************
      SUBROUTINE EIGEN(C,NPVST,NPV,ieval,xlmin,cinv,iprin)
      IMPLICIT double precision (A-H,O-Z)
      Include 'fitex.h'
      PARAMETER (NPSTS=itsth*npsth)

      DIMENSION C(NPVST,NPVst),cinv(npvst,npvst)
      Dimension EVECTS(NPSTS,NPSTS),xev(NPSTS),
     &        TEMP(NPSTS,NPSTS),EVALS(NPSTS),work(NPSTS)

      ievals=iabs(ieval)
      if(NPSTS.lt.npv.or.npvst.lt.npv) stop 'Eigen: storage error'
      if(ievals.gt.npv.and.ieval.ne.99) stop 'Eigen: ievals too big'

      DO 10 I=1,NPV
      DO 10 J=1,NPV
  10  TEMP(I,J)=C(I,J)

      IFAIL=0
c Real symmetric (uses lower triangle): eigenvalues ascending order
c  and normalised eigenvectors

      Call f02abf(temp,NPSTS,npv,evals,evects,NPSTS,work,ifail)
       if(ifail.ne.0.and.iprin.gt.0)write(6,919)ifail
 919   format(' error in Eigen: Ifail=',i2)

      if(iprin.eq.-1) then
       xlmin=evals(1)
       return
      endif
      if(evals(1).le.0.0.and.iprin.gt.0)then
       write(6,95)npv,evals(1)
 95    format(' CARE:  ',i2,'th eigenvalue =',f12.5,' NEGATIVE')
      endif
      
      if(iprin.gt.1)write(6,981)
 981   format('  eigenvalues ')
      IF(IPRIN.GT.1)WRITE(6,98)NPV,(EVALS(IP),IP=1,NPV)
      if(iprin.gt.7)write(6,982)
 982   format('  eigenvectors ')
      do 483 i=1,npv
      IF(IPRIN.GT.7)WRITE(6,98)I,(EVECTS(IP,I),IP=1,NPV)
 483  continue   
 98   FORMAT(I6,6F12.6,7(/,6X,6f12.6))

      do 200 ip1=1,npv
      do 201 ip2=1,npv
       temp(ip1,ip2)=0.0
       cinv(ip1,ip2)=0.0
 201  continue
 200  continue 

      if(ieval.lt.npv.and.ieval.ge.0)then
       xevs=0.0
        do 150 ie=1,npv-ievals
         xevs=xevs+evals(ie)
 150    continue
       xlmin=xevs/(npv-ievals)
       if(iprin.ge.1) write(6,97)npv-ievals,evals(1),
     &   evals(npv-ievals),xlmin
 97    format(' Inversion:',i3,' smallest eigenvalues ',f10.6,
     &  ' to', f10.6,/' replaced by max(value, average): ',  f10.7)
       if(xevs.le.0.0.and.iprin.ge.1)write(6,902)
 902   format(' warning: average eigenvalue zero !')
      else if(ieval.eq.npv) then
c full inversion
       if(iprin.ge.1) write(6,92) evals(1)
 92     format('  Full inversion: smallest eigenvalue=',f10.7)
       xlmin=0.0
      else if(ieval.eq.99)then
c  ievals < 0 so just use xlmin here 
       if(iprin.ge.1) write(6,93) xlmin
 93     format(' Inversion: eigenvalues <',f10.6,' replaced by  ', 
     &  'that value')
      else if(ieval.lt.0)then
       if(iprin.ge.1) write(6,94)npv-ievals,evals(1),evals(npv-ievals)  
 94     format(' SVD inverse: ',i2,' smallest eigenvalues ',f10.6,
     &' to ',f10.6,/' replaced by zero')
      else
       stop ' eigen: wrong ievals selection'  
      endif

c  modify eigenvales as needed
        xnsum=0.0
        xosum=0.0
      do 155 ie=1,npv
         xl=evals(ie)
c replace xl < xlmin by xlmin
         if(ieval.eq.99.and.xl.lt.xlmin) xl=xlmin
c replace npv-ievals smallest evals by average if less
         if(ieval.ge.0.and.ie.le.(npv-ievals).and.xl.lt.xlmin)
     &                     xl=xlmin
c replace npv-ievals smallest eigenvalues by SVD (zero in inverse)
         if(ieval.lt.0.and.ie.le.(npv-ievals)) xl=0.0
         xev(ie)=xl
         xnsum=xnsum+xl
         xosum=xosum+evals(ie)
 155  continue

c sum of eigenvalues of 
         xnorm=xosum/xnsum

       if(xosum.ne.xnsum.and.iprin.ge.1) write(6,909) 
     &     (xev(ie)*xnorm,ie=1,npv)
 909    format(' New eigenvalues: ',7(/,6X,6f12.6)  )

      do 100 ie=1,npv
       do 110 ip1=1,npv
        do 120 ip2=1,npv
         xl=xev(ie)*xnorm
         if(xl.le.0.0) then 
c replace npv-ievals smallest eigenvalues by SVD (zero in inverse)
          xlinv=0.0
         else
          xlinv=1.0/xl
         endif
         temp(ip1,ip2)=temp(ip1,ip2)+
     &       evects(ip1,ie)*evects(ip2,ie)*xl

         cinv(ip1,ip2)=cinv(ip1,ip2)+
     &       evects(ip1,ie)*evects(ip2,ie)*xlinv
 120    continue
 110   continue
 100  continue
      npstm=NPSTS-1
      npvm=npv-1
      if(iprin.ge.4) call cdisplay(temp,0,npvm,npstm,5)
      RETURN
      END
c********************************************************
c
c  FT (on lattice) of 1/(p^2+m^2)
c
c also include form-factor exp(-a p^2) - but should not be necessary
c
c*********************************************************
      function x4ft(xmpi,it,ir,l,lx4)
      implicit double precision (a-h,o-z)
       twopi=8.0d0*atan(1.0d0)
       anorm=1.0d0/dfloat(l*l*l*lx4)
       c=0.0
        lx4h=lx4/2
        lh=l/2
         xpf=twopi/dfloat(l)
       do itv=0,lx4h
       xtct=2.0d0
       if(itv .eq. 0 .or. itv .eq. lx4h) xtct=1.0d0
        xpt=itv*twopi/dfloat(lx4)
        xden=xmpi*xmpi+(2.0d0*sin(xpt*0.5d0) )**2
c cos is even - so fold 
       do ixv=0,lh
        xxct=2.0d0
       if(ixv .eq. 0 .or. ixv .eq. lh) xxct=1.0d0
        xpx=ixv*xpf
        xden1=xden+(2.0d0*sin(xpx*0.5d0) )**2
         ccv=xtct*xxct*cos(xpt*dfloat(it))*cos(xpx*dfloat(ir)) 

c fold to speed up sin vanishes at 0,lh - so double [was -lh,lh-1 ]
       do iyv=0,lh
       xyct=2.0d0
       if(iyv .eq. 0 .or. iyv .eq. lh) xyct=1.0d0
        xpy=iyv*xpf
        xden2=xden1+(2.0d0*sin(xpy*0.5d0) )**2
       do izv=0,lh
       xzct=2.0d0
       if(izv .eq. 0 .or. izv .eq. lh) xzct=1.0d0
        xpz=izv*xpf
        xden3=xden2+(2.0d0*sin(xpz*0.5d0) )**2

        c=c+ xyct*xzct*ccv /xden3

       Enddo
       Enddo
       Enddo

       Enddo

       
       x4ft=c*anorm

       return
       End

cccccccccccccccccccccccccccccccccc
c
c  interpolated x4ft  -- for speed
c
cccccccccccccccccccccccccccccccccc
      function x4fti(xm,itv,ir,lx1,lx4)
      implicit double precision (a-h,o-z)
      integer it,im,ir,itv,init,nmass,im0,im1,itmax,lx1,lx4
 
      parameter (nmass=200, xm0=0.05, xmax=0.55, itmax=24)
      dimension xval(nmass,itmax)
      save xval, init
      data init /0/

       dm=(xmax-xm0)/float(nmass-1)
      if (init .eq. 0 ) then
      init=1
      do im=1,nmass 
       xmpi=xm0+(im-1)*dm

       do it=1,itmax
        val= x4ft(xmpi,it,ir,lx1,lx4)
         xtr=sqrt(float(it*it+ir*ir) ) 
c ad hoc interpolation formula
         xnaive=xtr**(-1.5)*exp(-xmpi*xtr)*sqrt(xmpi)
         valrat=val/xnaive
c       write(6,90)it,xmpi,val,valrat
90      format(i3,f12.5,e15.5,f12.5)
        xval(im,it)=valrat
       enddo

       enddo
      endif



c find grid points and linear interpolate if inside box
       xtr=sqrt(float(itv*itv+ir*ir) ) 
       xmint = (xm-xm0)/dm

       if( xmint .lt. nmass .and. xmint  .gt. 0 .and. 
     &   itv .le. itmax .and. itv .gt. 0) then

       im0=xmint+1
         xmp0=xm0+(im0-1)*dm
         xml=xval(im0,itv)*xtr**(-1.5)*exp(-xmp0*xtr)*sqrt(xmp0)

       im1=im0+1
         xmp1=xm0+(im1-1)*dm
         xmu=xval(im1,itv)*xtr**(-1.5)*exp(-xmp1*xtr)*sqrt(xmp1)
      
c         write(*,*) xmp0,xm,xmp1
c         write(*,*) xml,(dxm/dm),xmu

         dxm=(xmint-float(im0-1))*dm
         val=(xml*(dm-dxm)+xmu*(dxm) ) /dm

c         valx4= x4ft(xm,itv,ir,lx1,lx4)
c         write(6,91)itv,xm,xval,valx4
91       format(i3,3g15.6)

         else

         val=x4ft(xm,itv,ir,lx1,lx4)

         endif

         x4fti=val

         return
         End

