program call_get_err_band
      real*8 w,q2,a1fit,a1fiterr
      integer first
c
      
      q2 = 1.3
      first = 1
      do i=1,92
         w = 1.080+(i-1)*.010
      call get_FIT2(w,q2,a1fit,a1fiterr,first)
      first = 2
      write(*,*) w,a1fit,a1fiterr
      enddo
      end
c
      subroutine  get_FIT2(w,q2,a1fit,a1fiterr,first)
      IMPLICIT  NONE

      INTEGER   NX,NUM_TERMS,IA1A2,nterms_fit
      INTEGER   J,NP,NPAR,NRC,i,k
      LOGICAL ERR_fit
      character*40 datfile,matfile
      character*40 fitoutfile,fit2outfile,paroutfile,parinfile

      PARAMETER (NX         = 100  )
c      PARAMETER (NUM_TERMS  = 17   )
      PARAMETER (NUM_TERMS  = 17   )
      PARAMETER (NRC        = 2    )
      PARAMETER (NPAR       = 17   )
c      PARAMETER (Nterms_fit      = 5   )
      PARAMETER (Nterms_fit      = 6   )

      integer index(nterms_fit),first
      real*8      temp_large
      integer    temp_i,temp_k
      REAL*8      YFIT(NX)
      REAL*8       PARS(NUM_TERMS)  
      REAL*8      EPARS1(NUM_TERMS),EPARS2(NUM_TERMS)
      REAL*8      ERR_MAT1(Nterms_fit,Nterms_fit),ERR_MAT2(Nterms_fit,Nterms_fit)

      REAL*8      A1A2FIT,A1_OR,A2_OR,CHISQR
      REAL*8      delta,r2,r3,r4,dis,err
      REAL*8      AR(4), WR(4), GR(4), BB(4), AA

      REAL*8      w,Q2,A1(NX,0:NRC),EA1(NX,0:NRC),A2(NX,0:NRC),EA2(NX,0:NRC)
      real*8 functn,a1fit,a1fiterr,a2fit
      real*8 wserr,wsum,wave,werr
      real*8 coverr(nx),a1d_lo,a1d_hi,a2d_lo,a2d_hi,deriv(Nx,Nterms_fit)
      real*8 derr,calc_deriv
      REAL*8 PARS1(NUM_TERMS),PARS1_INIT(NUM_TERMS)
      REAL*8 PARS2(NUM_TERMS),PARS2_INIT(NUM_TERMS)
      COMMON/FITPARS/AR,WR,GR,BB,AA,PARS1,PARS2,PARS1_INIT,PARS2_INIT

      !------------------------------------------------------------------------     
      if (first .eq. 1) then
      write(*,*) 'A1  Error matrix filename:'
      read(*,*) matfile
      OPEN(80,FILE=matfile,STATUS='UNKNOWN')
      do j=1,Nterms_fit
         read(80,'(6f12.7)') (err_mat1(i,j),i=1,Nterms_fit)
      enddo
      close(80)

      write(*,*) ' Starting parameter filename:'
      read(*,*) parinfile
      OPEN(2,FILE=parinfile,STATUS='old')
      DO J=1,NPAR
        READ(2,*) PARS1(J),PARS2(J)
      ENDDO
      CLOSE(2)

      endif
      
      index(1)=1
      index(2)=3
      index(3)=4
      index(4)=5
      index(5)=13
      index(6)=14
      IA1A2=1

    

      ! CUT
      DO J=1,4
        AR(J) = PARS1(J)
        WR(J) = PARS1(J+4)
        GR(J) = PARS1(J+8)
        BB(J) = PARS1(J+12)
      ENDDO
        
      AA      = PARS1(17)


c
c calculate derivatives
      j=1
         do i=1,Nterms_fit
            pars1(index(i))=pars1(index(i))*0.95 
            a1d_lo=functn(w,Q2,pars1,17,1)
            pars1(index(i))=pars1(index(i))/0.95*1.05 
            a1d_hi=functn(w,Q2,pars1,17,1)
            pars1(index(i))=pars1(index(i))/1.05 
            deriv(j,i) = (a1d_hi-a1d_lo)/(pars1(index(i))*.1)
         enddo
c
        j=1
         coverr(j) = 0.0
         do i=1,Nterms_fit
         do k=1,Nterms_fit
            coverr(j) = coverr(j)+deriv(j,i)*err_mat1(i,k)*deriv(j,k) 
         enddo
         enddo
         a1fiterr = sqrt(coverr(j))
c
         a1fit=A1A2FIT(W,Q2,1,4,1,delta,R2,R3,R4,DIS,ERR_fit)


      RETURN
      END   
!======================================================================
C
C        FUNCTION TO OUTPUT A1(W) AND A2(W) FITS TO RSS DATA
C        O. A. RONDON 10/2005
c
c  modify to do just delta resonance + DIS
c
C        INPUTS:
C        W [GEV], Q**2 [GEV**2], ENERGY LOSS NU [GEV]

        REAL*8 FUNCTION FUNCTN(X1,X2,A,NPARS,IA1A2)
!       REAL*8 FUNCTION A1A2FIT(W,ENU,QSQ,N,IMAX,JMAX)
        IMPLICIT NONE
        INTEGER NPARS,IA1A2
        REAL*8 X1,X2
        REAL*8 A(NPARS)
        REAL*8 AR(4), WR(4), GR(4), B(4),ALPHA
        REAL*8 BW(4), SUMBW, DIS
        INTEGER I,IMAX,J,JMAX,K,N
        REAL*8 D,R2,R3,R4
        REAL*8 W,XB,ENU,QSQ,M,MPI
        REAL*8 QCM,QCR(4),G(4),GG(4)
        REAL*8 KCM,KCR(4)
        INTEGER LR(4),JR(4)
        REAL*8 XR(4)
        LOGICAL ERR
        REAL*8 A1A2FIT

        DATA LR/ 1, 2, 3, 3/
        DATA JR/ 1, 1, 2, 2/
        DATA XR/ 0.160D0, 3*0.350D0/

        PARAMETER (M=0.93828D0)
        PARAMETER (MPI=0.134D0)

C       CHECK PI0 THRESHOLD
       
        W   = X1
        QSQ = X2
        
        ERR =.FALSE.
        IF(W.LT.(M+MPI)) THEN
          ERR=.TRUE.
          RETURN
        END IF

C       COPY PARAMETERS TO LOCAL VARIABLES
C       B-W'S

        DO J=1,4
          AR(J)=A(J)
          WR(J)=A(J+4)
          GR(J)=A(J+8)
        END DO
C DIS
        DO J=1,4
          B(J) = A(j+12)
        END DO
        ALPHA = A(17)        

C SUM OF B-W

        QCM = SQRT(((W**2+M**2-MPI**2)/2.D0/W)**2 - M**2)
!        Q1 = SQRT(((WR(I,K)**2+M**2-MPI**2)/2/WR(I,K))**2 - M**2)
        KCM = SQRT(((W**2+M**2+QSQ)/2.D0/W)**2 - M**2)
        
        SUMBW = 0.D0

        DO I=1,4
        QCR(I) = SQRT(((WR(I)**2+M**2-MPI**2)/2.D0/WR(I))**2 - M**2)
        KCR(I) = SQRT(((WR(I)**2+M**2+QSQ)/2.D0/WR(I))**2 - M**2)
        G(I) = GR(I)*(QCM/QCR(I))**(2*LR(I)+1)*
     &         ((QCR(I)**2+XR(I)**2)/(QCM**2+XR(I)**2))**LR(I)
        GG(I) = GR(I)*(KCM/KCR(I))**(2*JR(I))*
     &     ((KCR(I)**2+XR(I)**2)/(KCM**2+XR(I)**2))**JR(I)
        BW(I) = AR(I)*(KCR(I)/KCM)**2*WR(I)**2*G(I)*GG(I)
        BW(I) = BW(I)/((WR(I)**2-W**2)**2+(WR(I)*G(I))**2)
        SUMBW = SUMBW + BW(I)
        END DO

C SAVE BW RESULTS

        D = BW(1)
        R2 = BW(2)
        R3 = BW(3)
        R4 = BW(4)
        
C DIS POLY

        XB = QSQ/(W**2 - M**2 + QSQ)
        DIS = XB**ALPHA*(B(1) + XB*(B(2) + B(3)*XB + B(4)*XB**2))
        IF(IA1A2.EQ.2) THEN        ! ONLY FOR A2 = 1/Q DEPENDENCE
          IF(QSQ.NE.0.D0) THEN
            DIS = DIS/SQRT(QSQ)
          ELSE
            DIS = 0.D0
          END IF
        ENDIF

C A1 OR A2

        A1A2FIT = SUMBW + DIS
        FUNCTN = SUMBW + DIS
        
        RETURN
        END
C
C     FUNCTION TO OUTPUT A1(W) AND A2(W) FITS TO RSS DATA
C     O. A. RONDON 10/2005
C     INPUTS:
C     W [GEV], Q**2 [GEV**2], ENERGY LOSS NU [GEV]
C     N = WHICH ASYMMETRY: 1 = A1, 2 = A2
C     I = NUMBER OF B-W RESONANCES: 3 OR 4
C     J = ORDER OF DIS BACKGROUND POLYNOMIAL: 2 OR 3
C     ONLY VALID N,I,J COMBINATIONS ARE:
C     A1 = 1,4,4
C     A2 = 2,4,4.
C               
C     LOGICAL ERR CAN BE USED TO SKIP W VALUES BELOW PI0 THRESHOLD

      REAL*8 FUNCTION A1A2FIT(W,QSQ,N,IMAX,JMAX,D,R2,R3,R4,DIS,ERR)
      IMPLICIT NONE
      INTEGER NUM_TERMS
      PARAMETER(NUM_TERMS=17)
      REAL*8 W, AR(4), WR(4), GR(4), BB(4), AA
      REAL*8 BW(4), SUMBW, DIS
      INTEGER I,IMAX,J,JMAX,K,N
      REAL*8 D,R2,R3,R4
      REAL*8 XB,ENU,QSQ,M,MPI
      REAL*8 QCM,QCR(4),G(4),GG(4)
      REAL*8 KCM,KCR(4)
      INTEGER LR(4),JR(4)
      REAL*8 XR(4)
      LOGICAL ERR

      !B-W PARAMETERS: AR = AMPLITUDE, WR = CENTROID, GR = WIDTH
      !B : DIS BACKGROUND = X**A*(B0+(B1+B2*X+B3*X**2)*X)
      !A : EXPONENT
      REAL*8 PARS1(NUM_TERMS),PARS1_INIT(NUM_TERMS)
      REAL*8 PARS2(NUM_TERMS),PARS2_INIT(NUM_TERMS)
      COMMON/FITPARS/AR,WR,GR,BB,AA,PARS1,PARS2,PARS1_INIT,PARS2_INIT


C      DATA AR/     -0.5166,      0.2384,      0.4132,      0.1818,
C     &             -0.1172,      0.1880,      0.1630,      0.1498,
C     &              0,           0,           0,           0,
C     &              0,           0,           0,           0        /
C
C
C      DATA WR/      1.2262,      1.3498,      1.5536,      1.7393,
C     &              1.2265,      1.3447,      1.4881,      1.6753,
C     &              0,           0,           0,           0,
C     &              0,           0,           0,           0        /
C
C      DATA GR/      0.2412,      0.0586,      0.2901,      0.1277,
C     &              0.1589,      0.0985,      0.1112,      0.1875,
C     &              0,           0,           0,           0,
C     &              0,           0,           0,           0        /
C
C
C      DATA B /     0.13455,     0.22084,     0.10051,    -0.09316,
C     &             0.08879,     0.07794,    -0.13864,    -0.08338,
C     &             0,           0,           0,           0,
C     &             0,           0,           0,           0        /
C
C      DATA A/-0.18062,0.43254,   0, 0/

      !-------------------------------------------------------
     
     
      DATA LR/ 1, 2, 3, 3/
      DATA JR/ 1, 1, 2, 2/

      DATA XR/ 0.160D0, 3*0.350D0/

      PARAMETER ( M   = 0.93828D0 )
      PARAMETER ( MPI = 0.134D0   )


      !CHECK PI0 THRESHOLD
      ERR = .FALSE.

      
      IF(W.LT.(M+MPI)) THEN
        ERR = .TRUE.
        RETURN
      END IF
      !WHICH FIT?
 
      IF(N.EQ.1) THEN
        K = 1                      !      ONLY FIT FOR A1
      ELSE                         !      A2 FITS
        IF(IMAX.EQ.3) THEN
          K = 3                    !      3 RESONANCES + QUAD DIS
        ELSE
          IF(JMAX.EQ.3) THEN
            K = 4                  !      4 RESONANCES + QUAD DIS
          ELSE
            K = 2                  !      4 RESONANCES + CUBIC DIS
          ENDIF                    !      END POLY ORDER
        ENDIF                      !      END 3 OR 4 B-W
      ENDIF                        !      END A1 OR A2                      

      !SUM OF B-W

      QCM    = SQRT(((W**2+M**2-MPI**2)/2.D0/W)**2 - M**2)
      KCM    = SQRT(((W**2+M**2+QSQ)/2.D0/W)**2 - M**2)

      SUMBW  = 0.D0

      DO I=1,4
         bw(I) = 0.
         ENDDO
C
      DO I=1,IMAX
        QCR(I) = SQRT(((WR(I)**2+M**2-MPI**2)/2.D0/WR(I))**2 - M**2)
        KCR(I) = SQRT(((WR(I)**2+M**2+QSQ)/2.D0/WR(I))**2 - M**2)
        G(I)   = GR(I)*(QCM/QCR(I))**(2*LR(I)+1)*
     &           ((QCR(I)**2+XR(I)**2)/(QCM**2+XR(I)**2))**LR(I)
        GG(I)  = GR(I)*(KCM/KCR(I))**(2*JR(I))*
     &          ((KCR(I)**2+XR(I)**2)/(KCM**2+XR(I)**2))**JR(I)
        BW(I)  = AR(I)*(KCR(I)/KCM)**2*WR(I)**2*G(I)*GG(I)
        BW(I)  = BW(I)/((WR(I)**2-W**2)**2+(WR(I)*G(I))**2)
        SUMBW  = SUMBW + BW(I)
      ENDDO


      !SAVE BW RESULTS

      D  = BW(1)
      R2 = BW(2)
      R3 = BW(3)
      R4 = BW(4)
      
      !DIS POLY

      !XB  = QSQ/(2.D0*M*ENU)
      !WRITE(66,*) XB
      XB = ((W**2-M**2)/QSQ + 1.0)**(-1)
      !WRITE(67,*) XB
      DIS = XB**AA*(BB(1) + XB*(BB(2) + BB(3)*XB + BB(4)*XB**2))


      IF(N.EQ.2) THEN              ! ONLY FOR A2 = 1/Q DEPENDENCE
        IF(QSQ.NE.0.D0) THEN
          DIS = DIS/SQRT(QSQ)
        ELSE
          DIS = 0.D0
        END IF
      ENDIF

      !A1 OR A2

      A1A2FIT = SUMBW + DIS
      RETURN
      END

c
      real*8 function calc_deriv(w,q2,h,pars1,index,err)
c
      implicit none
c
      integer ntab,i,j,index
      real*8 err,h,w,q2,pars1(17),functn
      real*8 con,con2,big,safe,a1d_lo,a1d_hi
      parameter (ntab=10,con=1.4,con2=con*con,big=1.d30,safe=2.)
      real*8 a(ntab,ntab),hh,fac,errt
c
      if (h.eq.0) return
      hh=h
      pars1(index)=pars1(index)-hh
      a1d_lo = functn(w,Q2,pars1,17,1)
      pars1(index)=pars1(index)+2*hh
      a1d_hi = functn(w,q2,pars1,17,1)
      pars1(index)=pars1(index)-hh
      a(1,1)=(a1d_hi-a1d_lo)/2./hh
      err=big
      do i=2,ntab
         hh=hh/con
         pars1(index)=pars1(index)-hh
         a1d_lo = functn(w,Q2,pars1,17,1)
         pars1(index)=pars1(index)+2*hh
         a1d_hi = functn(W,Q2,pars1,17,1)
         pars1(index)=pars1(index)-hh
         a(1,i)=(a1d_hi-a1d_lo)/2./hh
         fac=con2
         do j=2,i
            a(j,i)=(a(j-1,i)*fac-a(j-1,i-1))/(fac-1.)
            errt=max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
            if (errt .le. err) then
               err=errt
               calc_deriv=a(j,i)
            endif
          enddo
          if (abs(a(i,i)-a(i-1,i-1)) .ge. safe*err) return
      enddo
      return
      end