program match

      integer m,fit,cvm,list(15),i,n,fit

      real x(100),y(100),sg(100),e(15),cvr(100,100)
      real chi
      external fun

      n=50
      m=15
      fit=7 
      cvm=100
      
      list(1)=1
      list(2)=3
      list(3)=5
      list(4)=7
      list(5)=9
      list(6)=11
      list(7)=13
      list(8)=1
      list(9)=3
      list(10)=5
      list(11)=7
      list(12)=9            
      list(13)=11
      list(14)=13
      list(15)=15

      e(1)=0.0
      e(2)=0.0
      e(3)=0.0
      e(4)=0.0
      e(5)=0.0
      e(6)=0.0
      e(7)=0.0
      e(8)=0.0
      e(9)=0.0
      e(10)=0.0
      e(11)=0.0
      e(12)=0.0
      e(13)=0.0
      e(14)=0.0
      e(15)=1.0
      

      open(unit=20,file='match.in')

      do 10 i=1,n

        read(20,*) x(i),y(i)

        sg(i)=1.0

c        print*, x(i),y(i)

   10 continue

      call lfit(x,y,sg,n,e,m,list,fit,cvr,cvm,chi,fun)

      print*, 'chi =', chi

      do 20 i=1,m

        print*, i,'th coefficient =',e(i)

   20 continue



      

      close(20)

      stop
      end




      subroutine fun(xx,basis,mm)

      integer mm
      real xx,basis(mm),kl,ka,b

      kl=0.5
      ka=0.01

        b=(kl+xx)/sqrt(ka**2 + (kl+xx)**2)

        
        basis(1)= 1.0-b**(0.5)
        basis(2)= 1.0-b
        basis(3)= 1.0-b**(1.5)
        basis(4)= 1.0-b**(2.0)
        basis(5)= 1.0-b**(2.5)
        basis(6)= 1.0-b**(3.0)
        basis(7)= 1.0-b**(3.5)
        basis(8)= 1.0-b**(4.0)
        basis(9)= 1.0-b**(4.5)
        basis(10)= 1.0-b**(5.0)
        basis(11)= 1.0-b**(5.5)
        basis(12)= 1.0-b**(6.0)
        basis(13)= 1.0-b**(6.5)
        basis(14)= 1.0-b**(7.0)
        basis(15)= 1.0

        return
        end

      
      SUBROUTINE LFIT(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,COVAR,NCVM,CHISQ,
     *FUNCS)
      PARAMETER (MMAX=50)
      DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),A(MA),LISTA(MA),
     *    COVAR(NCVM,NCVM),BETA(MMAX),AFUNC(MMAX)
      KK=MFIT+1
      DO 12 J=1,MA
        IHIT=0
        DO 11 K=1,MFIT
          IF (LISTA(K).EQ.J) IHIT=IHIT+1
11      CONTINUE
        IF (IHIT.EQ.0) THEN
          LISTA(KK)=J
          KK=KK+1
        ELSE IF (IHIT.GT.1) THEN
          PAUSE 'Improper set in LISTA'
        ENDIF
12    CONTINUE
      IF (KK.NE.(MA+1)) PAUSE 'Improper set in LISTA'
      DO 14 J=1,MFIT
        DO 13 K=1,MFIT
          COVAR(J,K)=0.
13      CONTINUE
        BETA(J)=0.
14    CONTINUE
      DO 18 I=1,NDATA
        CALL FUNCS(X(I),AFUNC,MA)
        YM=Y(I)
        IF(MFIT.LT.MA) THEN
          DO 15 J=MFIT+1,MA
            YM=YM-A(LISTA(J))*AFUNC(LISTA(J))
15        CONTINUE
        ENDIF
        SIG2I=1./SIG(I)**2
        DO 17 J=1,MFIT
          WT=AFUNC(LISTA(J))*SIG2I
          DO 16 K=1,J
            COVAR(J,K)=COVAR(J,K)+WT*AFUNC(LISTA(K))
16        CONTINUE
          BETA(J)=BETA(J)+YM*WT
17      CONTINUE
18    CONTINUE
      IF (MFIT.GT.1) THEN
        DO 21 J=2,MFIT
          DO 19 K=1,J-1
            COVAR(K,J)=COVAR(J,K)
19        CONTINUE
21      CONTINUE
      ENDIF
      CALL GAUSSJ(COVAR,MFIT,NCVM,BETA,1,1)
      DO 22 J=1,MFIT
        A(LISTA(J))=BETA(J)
22    CONTINUE
      CHISQ=0.
      DO 24 I=1,NDATA
        CALL FUNCS(X(I),AFUNC,MA)
        SUM=0.
        DO 23 J=1,MA
          SUM=SUM+A(J)*AFUNC(J)
23      CONTINUE
        CHISQ=CHISQ+((Y(I)-SUM)/SIG(I))**2
24    CONTINUE
      CALL COVSRT(COVAR,NCVM,MA,LISTA,MFIT)
      RETURN
      END