SUBROUTINE CHISQ ( X, Y, N, SLOPE1, YSEPT1, CORCOE )
C
C       Performs a least squares fit of a line to some data
C       where the assumed uncertainties in the input data are
C       the same on both axes.
C
      IMPLICIT UNDEFINED (A-Z)
      SAVE
      INTEGER *4 N, I
      REAL *4 MDUST, MGAS, LFIR, LMID, X(*), Y(*)
      REAL *4 XSUM, YSUM, XXSUM, XYSUM, YYSUM, SIGMA, TDUST
      REAL *4 XXVAR, YYVAR, XYVAR, SLOPE1, SLOPE2, SLOPE3, YSEPT1, 
     $        YSEPT2, YSEPT3, XSEPT, 
     $        DELTA1, DELTA2, DELTA3, B, DELT1X, DELT2X, DELT3X
      REAL *4 SLOPE, YSEPT
      REAL *4 XMEAN, YMEAN, XSTDV, YSTDV, CORCOE,
     $        SIGM, DELTA, SIGFIT, YCALC, DELY
C.............................................................
      XSUM  = 0.
      YSUM  = 0.
      XXSUM = 0.
      XYSUM = 0.
      YYSUM = 0.
      DO 100 I = 1, N
          XSUM  = XSUM  + X(I)
          YSUM  = YSUM  + Y(I)
          XXSUM = XXSUM + X(I)*X(I)
          XYSUM = XYSUM + X(I)*Y(I)
          YYSUM = YYSUM + Y(I)*Y(I)
  100 CONTINUE
  150 CONTINUE
      XMEAN = XSUM  / FLOAT(N)
      YMEAN = YSUM  / FLOAT(N)
      XXVAR = XXSUM / FLOAT(N) - XMEAN*XMEAN
      YYVAR = YYSUM / FLOAT(N) - YMEAN*YMEAN
      XYVAR = XYSUM / FLOAT(N) - YMEAN*XMEAN
      XSTDV = SQRT ( XXVAR )
      YSTDV = SQRT ( YYVAR )
      SIGMA = ( XXVAR - YYVAR ) / ( 2.* XYVAR )
C      WRITE(6,*) '        ', N, ' POINTS '
C      WRITE(6,*) ' MEAN X       = ', XMEAN
C      WRITE(6,*) ' STD DEV OF X = ', XSTDV
C      WRITE(6,*) ' MEAN Y       = ', YMEAN
C      WRITE(6,*) ' STD DEV OF Y = ', YSTDV
C      WRITE(6,*) '  S I G M A   = ', SIGMA
      CORCOE =  XYVAR / SQRT ( XXVAR*YYVAR )
C      WRITE(6,*) ' CORRELATION COEFFIEICNT = ',CORCOE
C
C    A linear fit, with equal noise in x and y.
C
      SLOPE1 = - SIGMA + SQRT ( 1. + SIGMA*SIGMA )
      IF ( SLOPE1*XYVAR .LT. 0 ) THEN
          SLOPE1 = -SIGMA - SQRT ( 1. + SIGMA*SIGMA )
          IF ( SLOPE1*XYVAR .LT. 0 ) THEN
              WRITE(6,*) ' no solution is possible for ',
     $              ' equal weights in X and Y '
              slope1 = 0.
              ysept1 = 0.
              corcoe = 0.
              return
          END IF
      END IF
      YSEPT1 = YMEAN - SLOPE1 * XMEAN
      B = ( XXVAR - YYVAR ) / XYVAR
      DELTA1 = ( YYVAR + SLOPE1*XXVAR - 2.*SLOPE1*XYVAR ) / 
     $          SQRT( 1. + SLOPE1*SLOPE1 )
C      WRITE(6,*) ' The slope ( equal weights is ) ',SLOPE1
C      WRITE(6,*) ' The  y - intersept is ',YSEPT1
C      WRITE(6,*) ' The calculated rms is ', DELTA1
C      WRITE(6,*) '  '
C                      Fit line with standard least squares
C                      and fit y(x)
      YSEPT2 = ( XXSUM*YMEAN - XMEAN*XYSUM ) / (FLOAT(N)*XXVAR)
      SLOPE2 = ( XYSUM/FLOAT(N) - XMEAN*YMEAN ) / XXVAR
      DELTA2 = ( YYVAR + SLOPE2*XXVAR - 2.*SLOPE2*XYVAR )
C      WRITE(6,*) ' Std lst sqrs fit for y(x) gives '
C      WRITE(6,*) ' SLOPE = ', SLOPE2
C      WRITE(6,*) ' YSEPT = ', YSEPT2
C      WRITE(6,*) ' DELTA = ', DELTA2
C                    Fit line with standard least squares 
C                    but fit x(y)
      SLOPE3 = ( XYSUM/FLOAT(N) - XMEAN*YMEAN )/ YYVAR
      XSEPT  = ( YYSUM*XMEAN - YMEAN*XYSUM ) / (FLOAT(N)*YYVAR)
      YSEPT3 = - XSEPT / SLOPE3
      SLOPE3 = 1./ SLOPE3
      DELTA3 = ( YYVAR + SLOPE3*XXVAR - 2.*SLOPE3*XYVAR )
C      WRITE(6,*) ' Std lst sqrs fit forx(y) gives '
C      WRITE(6,*) ' slope = ', SLOPE3
C      WRITE(6,*) ' ysept = ', YSEPT3
C      WRITE(6,*) ' delta = ', DELTA3
      SIGFIT = 0.
      DELTA = 0.
      DELTA1 = 0.
      DO 200 I = 1, N
          DELTA1 = DELTA1 + ( Y(I) - SLOPE1*X(I) - YSEPT1) **2
          DELTA2 = DELTA2 + ( Y(I) - SLOPE2*X(I) - YSEPT2) **2
          DELTA3 = DELTA3 + ( Y(I) - SLOPE3*X(I) - YSEPT3) **2
          DELT1X = DELT1X + ( X(I) - (Y(I)-YSEPT1)/SLOPE1) **2
          DELT2X = DELT2X + ( X(I) - (Y(I)-YSEPT2)/SLOPE2) **2
          DELT3X = DELT3X + ( X(I) - (Y(I)-YSEPT3)/SLOPE3) **2
  200 CONTINUE
      DELTA1 = SQRT ( DELTA1 /FLOAT(N) )
      DELTA2 = SQRT ( DELTA2 /FLOAT(N) )
      DELTA3 = SQRT ( DELTA3 /FLOAT(N) )
      DELT1X = SQRT ( DELT1X /FLOAT(N) )
      DELT2X = SQRT ( DELT2X /FLOAT(N) )
      DELT3X = SQRT ( DELT3X /FLOAT(N) )
C      WRITE(6,*) ' observed deviations from the fits are '
C      WRITE(6,*) '                rms in x         rms in y '
C      WRITE(6,*) '  fit y(x)  ', DELT2X, DELTA2
C      WRITE(6,*) '  fit x(y)  ', DELT3X, DELTA3
C      WRITE(6,*) '  equal wts ', DELT1X, DELTA1
      RETURN
 1010 FORMAT ( A10, 5X, 3F15.0, F10.0, F7.0 )
 1020 FORMAT ( 5X, A10, 4F6.2 )
      END