      PROGRAM GFITJ ! MDISPERSION, Bernd Berg May 15, 2010.
C
C Jackknife Version of the LEVENBERG-MARQUARDT FIT PROGRAM. 
C
C TUNED FOR FITTING MASSES FROM CORRELATION FUNCTIONS.
C Explanation of output:
C   CJ="c0j": XMM  mass estimate (momentum k=0).
C   CJ="c1j": EML  energy estimate from correlation function.
C             EL   energy estimate from m and lattice momentum 2*sin(k1/2).
C             E1   energy estimate from m and k1=2pi/N momentum.
C             XM   mass used on input.
C             DEL  deviation (EL-EML)/EL in %.
C             XK1M momentum extrapolation k1=sqrt(EML**2-m**2).
C             XKL  lattice momentum l2=2*sin(k1/2).
C             XK1  k1=2pi/N momentum.
C             DEL  deviation (XKL-XK1m)/XKL in %.
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'
      include 'mc.par'
      CHARACTER CJ*3,CB*3
C Set input data file name: CJ.
      PARAMETER(NFIT=2, CJ="c1j",MDAT=100)
      PARAMETER(LFINAL=.false.,XFINAL=mth*ONE)
C     PARAMETER(LFINAL=.true., XFINAL=mth*ONE)
      DIMENSION X(MDAT),Y(MDAT),YE(MDAT),A(NFIT),AE(NFIT),
     &  COVA(NFIT,NFIT),AMAT(NFIT,NFIT),EJ(0:NRPT),XMJ(0:NRPT),
     &  XK1J(0:NRPT)
      EXTERNAL SUBG
      cfile="unused"
C
      DO IRPT=0,NRPT ! IRPT=0: Estimate from all data (for bias).
        WRITE(CB,'(I3.3)') IRPT
        OPEN(IUD1,FILE=CJ//CB//".d",FORM='FORMATTED',STATUS='OLD')
        READ(IUD1,*) K ! READ DATA:
        READ(IUD1,*) NDAT,(A(I),I=1,NFIT)
        IF(IRPT.EQ.0) WRITE(IUO,'(/," NUMBER OF DATA:",I12)') NDAT
        IF(NDAT.GT.MDAT) STOP "GFITJ: MDAT too small."
        LPRI=.FALSE.
        IF(IRPT.EQ.0) THEN
          IF(K.EQ.1) THEN
            INQUIRE(FILE="xmj.d",EXIST=LEXIST)
            IF(.NOT.LEXIST) THEN
              PRINT*,"xmj.d file does not exist. Run CJ=c0j first."
            STOP "xmj.d file does not exist. Run CJ=c0j first."
            ENDIF
            OPEN(IUD2,FILE="xmj.d",FORM='FORMATTED',STATUS='OLD')
            READ(IUD2,*) XMJ
            CLOSE(IUD2)
          ELSE
            PRINT*,"This program runs only for K=1."
          STOP "This program runs only for K=1."
          ENDIF
          LPRI=.TRUE.
          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
        ENDIF
        DO I=1,NDAT
          READ(IUD1,*) X(I),Y(I),YE(I)
          IF(IRPT.EQ.0) WRITE(IUO,'(3G16.7)') X(I),Y(I),YE(I)
        END DO
        CLOSE(IUD1)
C PERFORM FIT:
        CALL FIT_G(NDAT,X,Y,YE,NFIT,A,COVA,AMAT,CHI2,Q,SUBG,LPRI)
        IF(IRPT.EQ.0) WRITE(IUO,'(/," RESULTS FOR FIT PARAMETERS:")')
        DO I=1,NFIT
          AE(I)=SQRT(COVA(I,I))
          IF(IRPT.EQ.0) WRITE(IUO,
     &      '(I6,".  A(I) =",G16.7,"  +/-",1G16.7)') I,A(I),AE(I)
        END DO
        IF(IRPT.EQ.0) WRITE(IUO,'(/," From all data:",/)')
        IF(IRPT.EQ.1) WRITE(IUO,'(/," Jacknife estimates:",/)')
C
        IF(K.EQ.1) THEN ! Momentum k1 state analysis,
          EJ(IRPT)=A(2)
          XK1J2=A(2)**2-XMJ(IRPT)**2
          IF(XK1J2.LT.ZERO) STOP "K=1: XK1J2.LT.ZERO"
          XK1J(IRPT)=SQRT(XK1J2)
          WRITE(IUO,'(4X,"IRPT,EJ,XK1J:",I6,2G16.7)') 
     &                    IRPT,EJ(IRPT),XK1J(IRPT)
        ELSE
          XMJ(IRPT)=A(2)
          WRITE(IUO,'(4X,"IRPT,XMJ:",I6,G16.7)') IRPT,XMJ(IRPT)
        ENDIF
      ENDDO
C
      WRITE(IUO,'(/,"  LAST CORRELATED FIT CHI2 =",G16.7)') CHI2
      IF(NDAT.GT.NFIT) WRITE(IUO,'("  CHI2 p.d.f.:",14X,G16.7)') 
     &                                CHI2/(NDAT-NFIT)
      WRITE(IUO,'("  GOODNESS Q OF THE FIT:",4X,G16.7,/)') Q
C
      CALL STEBJ0(NRPT,XMJ(1),XJM,XMV,XMME)
      CALL BIAS(NRPT,XMJ(1),XMJ(0),XMM,BIASM)
      CALL STEBJ0(NRPT,EJ(1),EJM,EV,EMLE)
      CALL BIAS(NRPT,EJ(1),EJ(0),EML,BIASE)
      CALL STEBJ0(NRPT,XK1J(1),XK1JM,XV,XK1ME)
      CALL BIAS(NRPT,XK1J(1),XK1J(0),XK1M,BIASK)
      XK1=TPI/ML         ! Lowest discrete momentum.
      XKL=TWO*SIN(XK1/2) ! Momentum on lattice.
      EL=SQRT(XMM**2+XKL**2)
      E1=SQRT(XMM**2+XK1**2)
      DEL=100*(EL-EML)/EL
      DEL=100*(XKL-XK1M)/XKL
      WRITE(IUO,'(2X,"Mass parameter on input XM = ",F11.6)') XM
      WRITE(IUO,'(2X,"MC mass     XMM, XMME, Bias =",2F11.6,F12.6)')
     &                            XMM, XMME, BIASM
      WRITE(IUO,'(2X,"MC energy 1 EML, EMLE, Bias =",2F11.6,F12.6)')
     &                            EML, EMLE, BIASE
      WRITE(IUO,'(2X,"Momentum K1, K1L: XK1, XKL = ",2F11.6)') XK1,XKL
      WRITE(IUO,'(2X,"Mass parameter squared XM2 = ",F11.6)') XM2
      DisXM2=EML**2-XKL**2
      DisU=(EML+EMLE)**2-XKL**2
      DisD=(EML-EMLE)**2-XKL**2
      DisXM2E=MAX(DisU-DisXM2,DisXM2-DisD)
      WRITE(IUO,'(2X,"Dispersion DisXM2=E**2-KL**2:",2F11.6)') 
     &                           DisXM2,DisXM2E
C
      STOP
      END

      include '../../Libs/Fortran/bias.f'
      include '../../Libs/Fortran/fit_g.f'
      include '../../Libs/Fortran/error_f.f'
      include '../../Libs/Fortran/fit_ggnu.f'
      include '../../Libs/Fortran/fit_gfinal.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/stebj0.f'
      include '../../Libs/Fortran/stmean.f'
      include '../../ForProg/Fit_g/subg_mass.f'
