      PROGRAM GFIT ! MASSFIT VERSION, Bernd Berg October 16 2009.
C
C LEVENBERG-MARQUARDT FIT PROGRAM. 
C
C TUNED FOR FITTING MASSES FROM CORRELATION FUNCTIONS.
C
C  NFIT:    NUMBER OF FIT PARAMETERS.
c  MDAT:    MAXIMUM NUMBER OF DATA.
c  NDAT:    ACTUAL NUMBER OF DATA (READ FROM DATA FILE).
C  A:       ARRAY OF FIT PARAMETERS (INITIAL GUESS READ FROM DATA FILE).
C  X,Y,YE:  DATA WITH THEIR STATISTICAL ERRORS.
C  LPRI:    PRINT CONTROL PARAMETER.
C  LFINAL,XFINAL: To extend the plot range of the fit.
C
      include '../../Libs/Fortran/implicit.sta'
      include '../../Libs/Fortran/constants.par'
      include 'lat.par'
      PARAMETER(NFIT=2,LPRI=.TRUE., IUO=6,IUD=10,MDAT=1000)
      PARAMETER(LFINAL=.false.,XFINAL=mth*ONE)
C     PARAMETER(LFINAL=.true., XFINAL=mth*ONE)
      DIMENSION X(MDAT),Y(MDAT),YE(MDAT),A(NFIT),AE(NFIT)
      DIMENSION COVA(NFIT,NFIT),AMAT(NFIT,NFIT)
      EXTERNAL SUBG
      cfile="unused"
C
C READ DATA:
C ==========
      OPEN(IUD,FORM='FORMATTED',STATUS='OLD')
      READ(IUD,*) K
      READ(IUD,*) NDAT,(A(I),I=1,NFIT)
      WRITE(IUO,'(/," NUMBER OF DATA:",I12)') NDAT
      IF(NDAT.GT.MDAT) STOP "GFIT: MDAT too small."
      WRITE(IUO,'(/," INITIAL GUESS FOR THE FIT PARAMETERS:")')
      DO I=1,NFIT
        WRITE(IUO,'(" I,A(I):",I6,G18.7)') I,A(I)
      END DO
      DO I=1,NDAT
        READ(IUD,*)           X(I),Y(I),YE(I)
        WRITE(IUO,'(3G16.7)') X(I),Y(I),YE(I)
      END DO
      CLOSE(IUD)
C
C PERFORM FIT:
C ============
      CALL FIT_G(NDAT,X,Y,YE,NFIT,A,COVA,AMAT,CHI2,Q,SUBG,LPRI)
C
      WRITE(IUO,'(/," RESULTS FOR THE FIT PARAMETERS:")')
      DO I=1,NFIT
        AE(I)=SQRT(COVA(I,I))
        WRITE(IUO,'(I6,".  A(I) =",G16.7,"  +/-",1G16.7)') I,A(I),AE(I)
      END DO
C
      IF(K.EQ.1) THEN
        XK1=TPI/ML         ! Lowest discrete momentum.
        XKL=TWO*SIN(XK1/2) ! Momentum on lattice.
        WRITE(IUO,'(/,"  XK1, XKL, EK =",3G15.6)') XK1,XKL,A(2)
        IF(A(2).LT.XKL) THEN
          WRITE(IUO,'("  Rest mass^2 negativ.",/)')
        ELSE
          XM=SQRT(A(2)**2-XKL**2)
          WRITE(IUO,'("  Rest mass positiv: XM =",3X,G15.6,/)') XM
        ENDIF
      ENDIF
C
      WRITE(IUO,'("  CHI2    OF THE FIT:",G16.7)') CHI2
      WRITE(IUO,'("  PROB. Q OF THE FIT:",G16.7)') Q
      IF(NDAT.GT.NFIT) WRITE(IUO,'("  CHI2 p.d.f.:",7X,G16.7)') 
     &                                CHI2/(NDAT-NFIT)
      WRITE(IUO,*) '     '
C
C PREPARE GNUPLOT PLOT:
C =====================
      CALL FIT_GGNU(IUO,IUD,NDAT,NFIT,X,Y,YE,A)
      IF(LFINAL) CALL FIT_GFINAL(IUO,IUD,NFIT,X(1),XFINAL,A)
C
      STOP
      END

      include '../../Libs/Fortran/fit_g.f'
      include '../../Libs/Fortran/error_f.f'
      include '../../Libs/Fortran/gamma_ln.f' 
      include '../../Libs/Fortran/gamma_p.f' 
      include '../../Libs/Fortran/gaudif.f' 
      include '../../Libs/Fortran/mat_gau.f'
      include '../../Libs/Fortran/fit_ggnu.f'
      include '../../Libs/Fortran/fit_gfinal.f'
      include '../../ForProg/Fit_g/subg_mass.f'
