      SUBROUTINE LINCOR(X,Y,N,TITLE)     
C     
C     ROUTINE TO FIT A STRAIGHT LINE OR A PARABOLA OR A CUBIC
C
C     NUM IS THE NUMBER OF TERMS.
C     
      COMMON /PLOTFL/PLNAME
      CHARACTER PLNAME*12
      REAL*8 CM(6,6),SIG
      dimension iwork2(6)
      DOUBLE PRECISION CONST,SLOPE,SECOND,THIRD,YVAL(1000),XVAL(1000)
      DOUBLE PRECISION COV(6,6),COVINV(6,6),UNIT(6,6),SIGY(1000)
      DOUBLE PRECISION DERIV(6),YVEC(6),SOLN(6),RESID,E,C
      DOUBLE PRECISION XTEMP,YTEMP,SGTEMP
      DIMENSION JUNK(6),Y(1),X(1),XBIN(1000),YBIN(1000),NBIN(1000),
     *                            SIGBIN(1000)
      CHARACTER FNAME*25,TITLE*80,FORMAT*80
      CHARACTER YTITLE*40,XTITLE*40 
      FORMAT=' '
      IDEBUG=0
      CONST=0.0D0 
      SLOPE=0.0D0  
      SECOND=0.0D0  
      THIRD=0.0D0
      INTER=0
    4 READ 5,FNAME
    5 FORMAT(A)
      IF(FNAME.EQ.'FORMAT'.OR.FNAME.EQ.'format')THEN
C
C       READ IN THE DATA FILE INPUT FORMAT
C
        READ 5,FORMAT
C
C       NOW READ THE DATA FILE NAME (AGAIN!)
C
        GO TO 4
      END IF
      IF(FNAME.EQ.'REVERS'.OR.FNAME.EQ.'revers')THEN
C
C       INTERCHANGE ORDER OF X AND Y IN DATA FILE
C
        INTER=1
C
C       NOW READ THE DATA FILE NAME (AGAIN!)
C
        GO TO 4
      END IF
      IBIN=0
      READ *,NUM  
      IF(NUM.LT.0)THEN
        IBIN=1
        NUM=-NUM
        READ *,BIN1,BIN2,NB
        IF(NB.GT.1000)STOP 'TOO MANY BINS REQUESTED'
      END IF
      IF(NUM.LE.0.OR.NUM.GE.4)NUM=2
      IPOLY=NUM   
      PRINT 10,FNAME,IPOLY-1    
   10 FORMAT(12X,'DATA TO BE READ FROM FILE ',A/
     *       12X,'POLYNOMIAL OF ORDER ',I2,' FITTED') 
      IF(IBIN.EQ.1)PRINT 11,BIN1,BIN2,NB
   11 FORMAT(12X,'DATA TO BE BINNED INTO THE ABSCISSA INTERVAL'/
     *       12X,F9.2,' TO ',F9.2,' IN ',I4,' BINS')
      PRINT 13    
   13 FORMAT(//)  
      CLOSE(11)   
      OPEN(11,FILE=FNAME,STATUS='OLD')     
      READ(11,*)N 
      IF(INTER.EQ.0)THEN
        IF(FORMAT.EQ.' ')THEN
          READ(11,*)(YVAL(J),SIGY(J),XVAL(J),J=1,N)
        ELSE
          READ(11,FORMAT)(YVAL(J),SIGY(J),XVAL(J),J=1,N)
        END IF
      ELSE
        IF(FORMAT.EQ.' ')THEN
          READ(11,*)(XVAL(J),YVAL(J),J=1,N)
        ELSE
          READ(11,FORMAT)(XVAL(J),YVAL(J),J=1,N)
        END IF
      END IF
      REWIND 11   
      CLOSE(11)   
C
C     BUBBLE SORT THE DATA
C
   15 IFLAG=0
      DO J=2,N
        IF(XVAL(J).LT.XVAL(J-1))THEN
          IFLAG=1
          XTEMP=XVAL(J)
          YTEMP=YVAL(J)
          SGTEMP=SIGY(J)
          XVAL(J)=XVAL(J-1)
          YVAL(J)=YVAL(J-1)
          SIGY(J)=SIGY(J-1)
          XVAL(J-1)=XTEMP
          YVAL(J-1)=YTEMP
          SIGY(J-1)=SGTEMP
        END IF
      END DO
      IF(IFLAG.EQ.1)GO TO 15

      IF(IBIN.EQ.1)THEN
C
C       BIN THE DATA
C
        DO J=1,1000
          XBIN(J)=0.0
          YBIN(J)=0.0
          SIGBIN(J)=0.0
          NBIN(J)=0.0
        END DO
C
C       START BINNING
C
        DO J=1,N
          INDEX=1+NINT((XVAL(J)-BIN1)/(BIN2-BIN1)*(NB-1))
          IF(INDEX.GT.0.AND.INDEX.LE.NB)THEN
            XBIN(INDEX)=BIN1+(INDEX-1.0)/(NB-1.0)*(BIN2-BIN1)
            YBIN(INDEX)=YBIN(INDEX)+YVAL(J)
            NBIN(INDEX)=NBIN(INDEX)+1
            SIGBIN(INDEX)=SIGBIN(INDEX)+YVAL(J)*YVAL(J)
          END IF
        END DO
        NEWN=0
        DO J=1,NB
C
C         CALCULATE THE MEAN OF EACH BIN AND ITS ACTUAL ST DEV
C         AND COPY TO THE XVAL AND YVAL DATA ARRAYS
C
          IF(NBIN(J).GT.0)THEN
            NEWN=NEWN+1
            YBIN(J)=YBIN(J)/NBIN(J)
            IF(NBIN(J).GT.1)
     *        SIGY(NEWN)=SQRT((SIGBIN(J)-NBIN(J)*YBIN(J)*YBIN(J))/
     *                                  (NBIN(J)-1)              )
            XVAL(NEWN)=XBIN(J)
            YVAL(NEWN)=YBIN(J)
            IF(SIGY(NEWN).GT.50.0)PRINT 77,J,SIGBIN(J),NBIN(J),
     *                            YBIN(J),XBIN(J)
   77       FORMAT(1X,'J=',I4,'  SIGBIN(J)=',F12.2,'  NBIN(J)= ',I5,
     *             '   YBIN(J)=',F12.2,'  XBIN(J)=',F7.1)
          END IF
        END DO
        N=NEWN
        PRINT 20,N
   20   FORMAT(12X,'BINNING INTO ',I4,' BINS COMPLETED')
      END IF

C     
C     CALCULATE THE NORMAL EQUATION MATRIX
C     
      DO 30 J=1,NUM     
        YVEC(J)=0.0  
        DO 30 K=1,NUM   
          COV(J,K)=0.0  
          CM(J,K)=0.0   
   30 CONTINUE    
      DO 50 I=1,N 
        E=XVAL(I)
        DERIV(1)=1.0    
        DO 32 J=2,NUM   
          DERIV(J)=DERIV(J-1)*E   
   32   CONTINUE  
        DO 40 J=1,NUM   
          DO 35 K=1,NUM 
            COV(J,K)=COV(J,K)+DERIV(J)*DERIV(K) 
   35     CONTINUE
          YVEC(J)=YVEC(J)+DERIV(J)*YVAL(I)    
   40   CONTINUE  
   50 CONTINUE    
C     
C     INVERT THE NORMAL MATRIX
C     
      DO 60 J=1,NUM     
        DO 60 K=1,NUM   
          COVINV(J,K)=COV(J,K)
   60 CONTINUE    
C     
C     INVERT THE NORMAL MATRIX USING DINV.FOR
C     
      call dinv(num,covinv,6,iwork2,ifail)
      if (ifail.ne.0) then
        print *,'MATRIX INVERSION FAILURE'
        print *,'ifail = ',ifail
        ierror=1
        idebug=1
      endif
C     
C     FROM HERE DOWN TO ABOUT LABEL 135 ARE DEBUG CALCULATIONS & PRINTS 
C     
   65 IF(IDEBUG.EQ.1)THEN     
        PRINT 70  
   70   FORMAT(//' COVARIANCE MATRIX IS : '/)   
        DO 90 I=1,NUM   
          PRINT 80,(COV(I,J),J=1,NUM)     
   80     FORMAT(1X,1P10E13.4)
          PRINT 80
   90   CONTINUE  
        PRINT 95  
   95   FORMAT(///' INVERSE IS : '/)
        DO 100 I=1,NUM  
          PRINT 80,(COVINV(I,J),J=1,NUM)  
          PRINT 80
  100   CONTINUE  
        PRINT 105 
  105   FORMAT(///' THE Y VECTOR IS : '/) 
        PRINT 80,(YVEC(J),J=1,NUM)     
        PRINT 80  
C     
C       TEST FOR GOOD MATRIX INVERSION BY PRODUCING A UNIT MATRIX 
C       AND PRINTING IT AND THEN CALCULATE THE CONTRIBUTION TO THE
C       SOLUTIONS BY CALCULATING A MATRIX GIVEN BY COVINV(I,J)*YVEC(J)     
C     
        PRINT 120 
        DO 125 I=1,NUM  
          DO 115 K=1,NUM
            UNIT(I,K)=0.0     
            DO 110 J=1,NUM    
              UNIT(I,K)=UNIT(I,K)+COVINV(I,J)*COV(J,K)
  110       CONTINUE    
  115     CONTINUE
  120     FORMAT(///' PRODUCT OF COVARIANCE AND ITS INVERSE IS : '/)    
          PRINT 80,(UNIT(I,K),K=1,NUM)    
          PRINT 80
  125   CONTINUE  
        PRINT 128 
  128   FORMAT(///' PRODUCT OF COVARIANCE INVERSE AND Y VECTOR :'/)     
        DO 135 I=1,NUM  
          DO 130 J=1,NUM
            UNIT(I,J)=COVINV(I,J)*YVEC(J)    
  130     CONTINUE
          PRINT 80,(UNIT(I,J),J=1,NUM)    
          PRINT 80
  135   CONTINUE  
        IF(IERROR.EQ.1)STOP   
      END IF
C     
C     CALCULATE THE SOLUTION  
C     
      DO 140 I=1,NUM    
        SOLN(I)=0.0     
        DO 140 J=1,NUM  
          SOLN(I)=SOLN(I)+COVINV(I,J)*YVEC(J)
  140 CONTINUE    
      SIG=0.0     
C     
C     CALCULATE THE NEW SIGMA OF THE RESIDUALS  
C     
      AVEX=0.0
      AVEY=0.0
      DO 150 I=1,N
        X(I)=XVAL(I)
        Y(I)=YVAL(I)
        AVEX=AVEX+XVAL(I)
        AVEY=AVEY+YVAL(I)
        RESID=YVAL(I) 
        C=0.0     
        DO 145 J=1,NUM  
            IF(XVAL(I).EQ.0.0.AND.J.EQ.1)RESID=RESID-SOLN(1)     
            IF(XVAL(I).NE.0.0)RESID=RESID-SOLN(J)*(XVAL(I)**(J-1))    
            IF(XVAL(I).EQ.0.0)C=SOLN(1)  
            IF(XVAL(I).NE.0.0)C=C+SOLN(J)*(XVAL(I)**(J-1))
  145   CONTINUE  
        IF(I.EQ.1)PRINT 13    
        IF(I.EQ.1)PRINT 146
  146   FORMAT(
     *       12X,'         X                   Y           PREDICTED',  
     *       5X,'     RESID'/ 
     *       12X,'----------        ------------        ------------',  
     *       5X,'----------') 
        PRINT 148,XVAL(I),YVAL(I),C,RESID    
  148   FORMAT(12X,F10.0,2(5X,F15.6),5X,F10.6)  
        SIG=SIG+RESID**2
  150 CONTINUE    
      AVEX=AVEX/N
      AVEY=AVEY/N
      SSX=0.0
      SSY=0.0
      SXY=0.0
      DO 160 I=1,N
        SSX=SSX+(XVAL(I)-AVEX)**2
        SSY=SSY+(YVAL(I)-AVEY)**2
        SXY=SXY+(XVAL(I)-AVEX)*(YVAL(I)-AVEY)
  160 CONTINUE
      RCORR=SXY/SQRT(SSX*SSY)

      SIG=SIG/(N-NUM)   
      SIGMA=SQRT(SIG)   
C     
C     CALCULATE THE CORRELATION MATRIX FOR THE ERRORS 
C     
      DO 200 I=1,NUM    
        DO 200 J=1,NUM  
  200     CM(I,J)=COVINV(I,J)*SIG   
      PRINT 250,SIGMA,RCORR,(SOLN(I),SQRT(CM(I,I)),I=1,IPOLY)     
  250 FORMAT(/12X,'SIGMA OF RESIDUALS      ',1PE9.1/
     *        12X,'CORRELATION COEFFICIENT ',0PF9.2//
     *        12X,33X,'SIGMA'/12X,38('-')/
     *        12X,'CONST',4X,0PF15.6,5X,F9.6/   
     *        12X,'SLOPE ',1X,1PE17.10,2X,E12.1/ 
     *        12X,'SECOND',3X,1PE15.3,5X,E9.1/  
     *        12X,'THIRD',4X,1PE14.3,5X,E9.1)   
      PRINT 260   
  260 FORMAT(12X,53('-'))     
      CONST=SOLN(1)     
      SLOPE=SOLN(2)
      SECOND=SOLN(3)
      THIRD=SOLN(4)    
      IEB=1
      READ *,XL,XMIN,XMAX,XTITLE
      READ *,YL,YMIN,YMAX,YTITLE
      NX=0
      NY=0
      DO J=40,1,-1
        IF(NX.EQ.0.AND.XTITLE(J:J).NE.' ')NX=J
        IF(NY.EQ.0.AND.YTITLE(J:J).NE.' ')NY=J
      END DO
C
C     INITIALISE THE PLOTTING
C
c     CALL SETCC
c     CALL PLFILE(PLNAME)
c     CALL PLOTS(IERR,IDUMMY,28)
c     IF(IERR.NE.0) STOP 'ERROR ON OPENING SPOOLTEK CALL'
c     CALL FACTOR(2.0) 
c     CALL FRAME(XL+4.0,YL+4.0)          
c     CALL PLOT(1.5,1.5,-3)   
c     CALL SYMBOL(0.,YL+1.5,0.2,TITLE,0.,60)  
      X(N+1)=XMIN
      Y(N+1)=YMIN
      X(N+2)=(XMAX-XMIN)/XL
      Y(N+2)=(YMAX-YMIN)/YL
      TT=X(N+1)    
      TTT=X(N+2)   
      YY=Y(N+1)
      YYY=Y(N+2)
c     CALL AXIS(0.0,0.0,XTITLE,-NX,XL,0.,TT,TTT) 
c     CALL AXIS(0.0,YL,' ',1,XL,0.,TT,TTT)   
c     CALL AXIS(0.,0.,YTITLE,NY,YL,90.0,YY,YYY) 
c     CALL AXIS(XL,0.,' ',-1,YL,90.0,YY,YYY)
      DO 300 J=1,N  
        XPOS=(X(J)-TT)/TTT 
        YPOS=(Y(J)-YY)/YYY
        ERROR=SIGY(J)/YYY     
c       CALL SYMBOL(XPOS,YPOS,0.05,1,0.0,-1)    
        IF (ERROR.EQ.0.0) GOTO 300
        TOP=YPOS+ERROR     
        BOTTOM=YPOS-ERROR  
        IF(TOP.GT.YL)TOP=YL   
        IF(BOTTOM.LT.0.0)BOTTOM=0.0 
c       CALL PLOT(XPOS+0.05,TOP,3)     
c       CALL PLOT(XPOS-0.05,TOP,2)     
c       CALL PLOT(XPOS,TOP,3)    
c       CALL PLOT(XPOS,BOTTOM,2) 
c       CALL PLOT(XPOS+0.05,BOTTOM,3)  
c       CALL PLOT(XPOS-0.05,BOTTOM,2)  
  300 CONTINUE    
c     CALL PLOT(0.0,0.0,999) 
      RETURN
      END   
