SUBROUTINE SYSVIS ( DEV )
C
C  This subroutine averages the observed visibility from several stars,
C  as a function of either time of night, or hour angle to determine
C  the variation of system visibility.
C
C    Modified  1 Nov 1988 to provide better smoothing
C             17 JUL 1989 to do rms about new fit
C             10 Mar 1991 fixed weighting in ftype 4 (was garbage)
C
      IMPLICIT UNDEFINED (A-Z)
      SAVE
      REAL *4 PI, RAD
      PARAMETER ( PI = 3.1415926535 )
      PARAMETER ( RAD= 180. / PI )
      INCLUDE 'VPLOT.INC'
      INTEGER *4 NS, IPLOT(50), I, J, K, L, ICOLOR, SYMBOL(5), LINE
      INTEGER *4 I1, I2, ISYS, N, IMARK(NFIT), MASK(4), IBASE, ISTAR
      INTEGER *4 ICHAN, NDATA, INEAR, IERR
      INTEGER *4 ISMOOTH, ORDER, LTYPE, NPTS, YTYPE, FTYPE, IFIRST
      INTEGER *4 ILAST, J1, RTYPE, STYPE, CALTYPE
      INTEGER *4 WHITE, RED, GREEN, YELLOW, BLUE, LIGHT_BLUE
      REAL *4 CMIN, CMAX, XMIN, XMAX, YMIN, YMAX, HMAX, FMAX, NMAX
      REAL *4 XDATA(NFIT), YDATA(NFIT), MINERR(4)
      REAL *4 ERROR(NFIT), FITRMS, FITMEAN, CHISQ, YFIT(NFIT)
      REAL *4 XLO, XHI, YLO, YHI, RX, RY, DISTANCE, D2, R1, RS
      REAL *4 X(NFIT), Y(NFIT), DISP, SCALE, CAL0, CAL1, TEMP
      REAL *4 XPLOT(NFIT), YPLOT(NFIT), ZPLOT(NFIT), ERROR(NFIT)
      REAL *4 VSYS(NFIT), TSYS(NFIT), VECX(2), VECY(2), COEF(NFIT)
      REAL *4 CRMS, CMEAN, CRANGE, DELTA, EXPO
      CHARACTER *(*)  DEV
      CHARACTER *1 CH
      CHARACTER *5 STARNUMB
      CHARACTER *75 STARLIST, TITLE(6)
      LOGICAL HAPLOT, CAL, HCOPY
      INCLUDE 'HARDCOPY.INC'
C-----------------------------------------------------------------------
      DATA SYMBOL  / 16, 3, 5, 6, 13 /
      DATA LINE, EXPO    / 29, 2.0 /
      DATA WHITE, RED, GREEN, YELLOW, BLUE, LIGHT_BLUE
     $       / 1,   2,     3,      6,    4,         11  /
      DATA TCAL0, TCAL1, HCAL0, HCAL1, ZCAL0, ZCAL1, MCAL0, MCAL1
     $    /   1.,   14.,   -8.,   16.,    0.,   70.,   0.,  180. /
      DATA MASK / 1, 2, 4, 8 /
      DATA ISMOOTH, ORDER, CALTYPE, YTYPE, LTYPE, FTYPE, RTYPE, STYPE
     $         / 1,     1,       1,     1,     2,     0,     1,   2  /
      HCOPY = .FALSE.
C-----------------------------------------------------------------------
C CALTYPE: 1 => Time,         2 => Seeing,  3 => Zenith Distance,
C          4 => Mirror angle, 5 => F.A.S.N, 6 => Photon rate.
C   LTYPE: 1 => POINTS,  2 => LINES,   3 => ERROR BARS
C   YTYPE: 1 => OBS VIS, 2 => CAL VIS, 3 => OBS/EST VIS 4=> CAL/EST VIS
C   FTYPE: 1 => BOXCAR,  2 => SMOOFT   3 => POLYNOMIAL  4=> SYMM POLY
C   STYPE: 1 => LASER,   2 => PHASE    RMS USED IN SEEING CALIBRATION.
C-----------------------------------------------------------------------
  100 CONTINUE
C CHANNEL SELECT
          WRITE(6,*) ' enter channel number '
          READ (5,*,END=900) ICHAN
          ICHAN = MIN ( MAX(1,ICHAN), 4 )
C BASELINE SELECT
          IF ( NBASE .EQ. 1 ) THEN
              IBASE = BLIST(1)
              WRITE(6,*) ' BASELINE NUMBER = ', IBASE
          ELSE
              WRITE(6,'(A,(5I5))') ' BASELINES ARE ',
     $                              (BLIST(I),I=1,NBASE)
              WRITE(6,*) ' ENTER BASELINE: '
              READ (5,*,END=900) IBASE
              DO J = 1, NBASE
                  IF ( BLIST(J) .EQ. IBASE ) GO TO 115
              END DO
              WRITE(6,*) ' BAD BASELINE '
              GO TO 100
  115         CONTINUE
          END IF
C STAR LIST SELECT
          NS = 0
          WRITE(6,'(5I5)')  ( SLIST(L),L=1,NSTAR )
  130     CONTINUE
          WRITE(6,*) ' Enter STAR NUMBER, 0 => end list  '
          WRITE(6,*) '                   -1 => use all stars '
          WRITE(6,*) '                   -2 => use last list '
          READ (5,*) ISTAR
          IF ( ISTAR .EQ. -2 ) THEN
              NS = NCALSTAR
              DO J = 1, NCALSTAR
                  IPLOT(J) = CALSTAR(J)
              END DO
          ELSE IF ( ISTAR .EQ. -1 ) THEN
              NS = NSTAR
              NCALSTAR = NS
              DO J = 1, NSTAR
                  IPLOT(J)   = SLIST(J)
                  CALSTAR(J) = SLIST(J)
              END DO
          ELSE IF ( ISTAR .GT. 0 ) THEN
              DO 140 J = 1, NSTAR
                  IF ( ISTAR .NE. SLIST(J) ) GO TO 140
                  NS = NS + 1
                  NCALSTAR    = NS
                  IPLOT(NS)   = ISTAR
                  CALSTAR(NS) = ISTAR
                  GO TO 130
  140         CONTINUE
              WRITE(6,*) ' STAR NOT IN STAR LIST '
              GO TO 130
          END IF
          IF ( ISTAR .NE. -2 ) THEN
              DO J = 1, ISCAN
                  BAD(J) = 0
              END DO
          END IF
C-----------------------------------------------------------------------
C    Return here after every plot
C
  145 CONTINUE
      IF ( HCOPY ) THEN
          NPLOTS = NPLOTS + 1
          WRITE( PLOTFILE(6:8),'(I3.3)') NPLOTS
          CALL PGBEGIN(0, PLOTFILE//PLOT_DEV, 1, 1)
      ELSE
          CALL PGBEGIN(0, '/EGA', 1, 1)
      END IF
C
  150 CONTINUE
C-----------------------------------------------------------------------
C  Make a list of scans to include
C
      NDATA = 0
      CMEAN = 0.
      CRMS  = 0.
      XMIN =  1.E6
      XMAX = -1.E6
      HMAX = 0.
      FMAX = 0.
      NMAX = 0.
      DO 450 J = 1, ISCAN
          IF ( IAND(OKSCAN(J),MASK(ICHAN)) .NE. 0 ) GO TO 450
          HMAX = MAX ( PHI(J,STYPE), HMAX )
          FMAX = MAX ( 0.5*BDEN(J,1)*VIS(J,1), FMAX )
          NMAX = MAX ( BDEN(J,ICHAN), NMAX )
          IF ( IAND(   BAD(J),MASK(ICHAN)) .NE. 0 ) GO TO 450
          IF (   BASE(J) .NE. IBASE   ) GO TO 450
          DO 440 K = 1, NS
              IF ( IPLOT(K) .EQ. STAR(J) ) THEN
                  NDATA = NDATA + 1
                  IF ( CALTYPE .EQ. 1 ) THEN
                      XDATA(NDATA) = HOURS(J)
                  ELSE IF ( CALTYPE .EQ. 2 ) THEN
                      XDATA(NDATA) = PHI(J,STYPE)
                  ELSE IF ( CALTYPE .EQ. 3 ) THEN
                      XDATA(NDATA) = ZD(J)
                  ELSE IF ( CALTYPE .EQ. 4 ) THEN
                      XDATA(NDATA) = MA(J)
                  ELSE IF ( CALTYPE .EQ. 5 ) THEN
                      XDATA(NDATA) = 0.5*BDEN(J,1)*VIS(J,1)
                  ELSE IF ( CALTYPE .EQ. 6 ) THEN
                      XDATA(NDATA) = BDEN(J,ICHAN)
                  ELSE
                      WRITE(6,*) ' BOZO CITY, YOU CANNOT GET HERE '
                  END IF
C                  XMIN = MIN ( XMIN, XDATA(NDATA) )
C                  XMAX = MAX ( XMAX, XDATA(NDATA) )
C
C   NOTE:  VIS is really visibility squared !
C
                  IF ( YTYPE .EQ. 1 ) THEN
                      YDATA(NDATA) = VIS (J,ICHAN)
                      ERROR(NDATA) = VSIG(J,ICHAN)
                  ELSE IF ( YTYPE .EQ. 2 ) THEN
                      YDATA(NDATA) = VCAL(J,ICHAN)
                      ERROR(NDATA) = VSIG(J,ICHAN) / SCAL(J,ICHAN)
                  ELSE IF ( YTYPE .EQ. 3 ) THEN
                      YDATA(NDATA) = VIS (J,ICHAN) / VEST(J,ICHAN)
                      ERROR(NDATA) = VSIG(J,ICHAN) / VEST(J,ICHAN)
                  ELSE IF ( YTYPE .EQ. 4 ) THEN
                      YDATA(NDATA) = VCAL(J,ICHAN) / VEST(J,ICHAN)
                      ERROR(NDATA) = VSIG(J,ICHAN) / 
     $                               ( SCAL(J,ICHAN) * VEST(J,ICHAN) )
                  END IF
                  ERROR(NDATA) = 
     $                  SQRT( ERROR(NDATA)**2 + MINERR(ICHAN)**2 )
                  CMEAN = CMEAN + YDATA(NDATA)
                  CRMS  = CRMS  + YDATA(NDATA)*YDATA(NDATA)
                  IMARK(NDATA) = J
              GO TO 450
              END IF
  440     CONTINUE
  450 CONTINUE
      IF ( NDATA .GT. 5 ) THEN
          CMEAN = CMEAN / FLOAT(NDATA)
          CRMS  = SQRT( CRMS/FLOAT(NDATA) - CMEAN*CMEAN )
      ELSE
          CMEAN = 0.
          CRMS  = 0.
      END IF
C-----------------------------------------------------------------------
C  Initialize the calibration arrays.
C
      DO I = 1, NFIT
          IF ( CALTYPE .EQ. 1 ) THEN
              TSYS(I) = TCAL0 + TCAL1*FLOAT(I)/FLOAT(NFIT)
          ELSE IF ( CALTYPE .EQ. 2 ) THEN
              HCAL0 = 0.D0
              HCAL1 = HMAX
              TSYS(I) = HCAL0 + HCAL1*FLOAT(I)/FLOAT(NFIT)
          ELSE IF ( CALTYPE .EQ. 3 ) THEN
              TSYS(I) = ZCAL0 + ZCAL1*FLOAT(I)/FLOAT(NFIT)
          ELSE IF ( CALTYPE .EQ. 4 ) THEN
              TSYS(I) = MCAL0 + MCAL1*FLOAT(I)/FLOAT(NFIT)
          ELSE IF ( CALTYPE .EQ. 5 ) THEN
              FCAL0 = 0.D0
              FCAL1 = FMAX
              TSYS(I) = FCAL0 + FCAL1*FLOAT(I)/FLOAT(NFIT)
          ELSE IF ( CALTYPE .EQ. 6 ) THEN
              NCAL0 = 0.D0
              NCAL1 = NMAX
              TSYS(I) = NCAL0 + NCAL1*FLOAT(I)/FLOAT(NFIT)
          ELSE
              WRITE(6,*) ' BOZO CITY, YOU CANNOT GET HERE '
          END IF
      END DO
C-----------------------------------------------------------------------
C      Sort the scans by the x axis.  This is to make the fits work 
C      properly and to make the plots look nice.
C
      IF ( NDATA .EQ. 0 ) GO TO 510
C
      CALL BUBBLE ( NDATA, XDATA, YDATA, ERROR, IMARK )
C
C   Determine the calibration curve.  The type of fit to use 
C   depends on FTYPE.
C         0 => no fit
C         1 => a gaussian smoothing of width ISMOOTH
C         2 => FFT smoothing, assuming the data is evenly spaced.
C         3 => a polynomial fit through the data.
C         4 => A symmetric polynomial, i.e. slope of zero at 0.
C
      IF ( CALTYPE .EQ. 1 ) THEN
          CAL0 = TCAL0
          CAL1 = TCAL1
      ELSE IF ( CALTYPE .EQ. 2 ) THEN
          CAL0 = HCAL0
          CAL1 = HCAL1
      ELSE IF ( CALTYPE .EQ. 3 ) THEN
          CAL0 = ZCAL0
          CAL1 = ZCAL1
      ELSE IF ( CALTYPE .EQ. 4 ) THEN
          CAL0 = MCAL0
          CAL1 = MCAL1
      ELSE IF ( CALTYPE .EQ. 5 ) THEN
          CAL0 = FCAL0
          CAL1 = FCAL1
      ELSE IF ( CALTYPE .EQ. 6 ) THEN
          CAL0 = NCAL0
          CAL1 = NCAL1
      ELSE
          WRITE(6,*) ' BOZO CITY, YOU CANNOT GET HERE '
          STOP
      END IF
      COEF(1) = 0.0
      COEF(2) = 0.0
      COEF(3) = 0.0
      IF ( FTYPE .EQ. 1 ) THEN
C
C  Type 1 fit to the data, a truncated, gaussian smoothing.
C  For a time cal, ISMOOTH is the half width in minutes
C                  R1      is the half width in hours
C
          ISMOOTH = MAX ( ISMOOTH, 10 )
          R1 = (CAL1-CAL0) * FLOAT(ISMOOTH) / 780.
          CALL GAUSS_SMOOTH ( XDATA, YDATA, YFIT, NDATA, R1,
     $                        TSYS, VSYS, NFIT, FITRMS )
C
      ELSE IF ( FTYPE .EQ. 2 ) THEN
          CALL MAKECAL ( XDATA, YDATA, YFIT, NDATA, TSYS, VSYS, NFIT,
     $                   ISMOOTH)
      ELSE IF ( FTYPE .EQ. 3 ) THEN
          CALL POLYFIT ( XDATA, YDATA, ERROR, YFIT, NDATA,
     $          0., 900., 0., 100., ORDER, TSYS, VSYS, NFIT,
     $          FITRMS  )
      ELSE IF ( FTYPE .EQ. 4 ) THEN
          CALL EVENPOLY ( XDATA, YDATA, ERROR, YFIT, NDATA,
     $           ORDER, COEF, TSYS, VSYS, NFIT, FITRMS, CHISQ )
      ELSE IF ( FTYPE .EQ. 5 ) THEN
          CALL EXPO_FIT ( XDATA, YDATA, ERROR, YFIT, NDATA, EXPO,
     $                      TSYS, VSYS, NFIT, COEF )
      ELSE
          DO J = 1, NFIT
              VSYS(J) = 1.0
          END DO
          DO J = 1, NDATA
              YFIT(J) = 1.
          END DO
      END IF
C-----------------------------------------------------------------------
C   quality of the fit statistics.
C   rms of the visibilities about the calibration curve
C   chi square of the fit.
C
      FITRMS  = 0.
      FITMEAN = 0.
      CHISQ   = 0.
      DO I = 1, NDATA
C
C   If the fitted data is too close to zero, it should not be used
C   in the rms, however, it should be included in the chi square.
C
          IF ( YFIT(I) .LE. 0.02 ) THEN
              WRITE(6,*) ' FIT IS TOO SMALL, DATA DELETED FROM RMS '
              RS = 0.
          ELSE
              RS = YDATA(I)/YFIT(I) - 1.
          END IF
          IF ( ERROR(I) .LT. .002 ) THEN
              WRITE(6,*) ' ERROR on point is awfully small ', error(i)
              DELTA = ( YDATA(I)-YFIT(I) ) / 0.002
          ELSE
              DELTA = ( YDATA(I)-YFIT(I) ) / ERROR(I)
          END IF
          FITMEAN = FITMEAN + RS
          FITRMS  = FITRMS  + RS * RS
          CHISQ   = CHISQ + DELTA*DELTA
      END DO
      IF ( NDATA .GE. 2 ) THEN
         FITMEAN = FITMEAN / FLOAT(NDATA)
         FITRMS  = SQRT (ABS( FITRMS/FLOAT(NDATA) - FITMEAN*FITMEAN ))
         CHISQ = CHISQ / FLOAT(NDATA)
      END IF
C
      CMIN = VSYS(1)
      CMAX = VSYS(1)
      DO J = 2, NFIT
          CMIN = MIN ( CMIN, VSYS(J) )
          CMAX = MAX ( CMAX, VSYS(J) )
      END DO
      CRANGE = 0.5 * ( CMAX - CMIN )
C
C   The fits end here
C-----------------------------------------------------------------------
C  Do the plotting.
C
  510 CONTINUE
      IF ( CALTYPE .EQ. 1 ) THEN
          TITLE(1) = 'TIME(HOURS)'
      ELSE IF ( CALTYPE .EQ. 2 ) THEN
          IF ( STYPE .EQ. 1 ) THEN
              TITLE(1) = 'LASER RMS'
          ELSE
              TITLE(1) = 'PHASE AUTOCORRELATION, 4 MS LAG'
          END IF
      ELSE IF ( CALTYPE .EQ. 3 ) THEN
          TITLE(1) = 'ZENITH DISTANCE (DEGREES)'
      ELSE IF ( CALTYPE .EQ. 4 ) THEN
          TITLE(1) = 'MIRROR ANGLE (DEGREES)'
      ELSE IF ( CALTYPE .EQ. 5 ) THEN
          TITLE(1) = 'FRINGE AMPLITUDE SIGNAL TO NOISE CHANNEL 1'
      ELSE IF ( CALTYPE .EQ. 6 ) THEN
          TITLE(1) = 'PHOTONS PER INTEGRATION THIS CHANNEL'
      ELSE
          TITLE(1) = ' '
      END IF
      IF ( YTYPE .EQ. 1 ) THEN
          TITLE(2) = 'OBSERVED VISIBILITY SQUARED'
      ELSE IF ( YTYPE .EQ. 2 ) THEN
          TITLE(2) = 'CALIBRATED VISIBILITY SQUARED'
      ELSE IF ( YTYPE .EQ. 3 ) THEN
          TITLE(2) = 'OBSERVED/ESTIMATED VIS SQUARED'
      ELSE IF ( YTYPE .EQ. 4 ) THEN
          TITLE(2) = 'CALIBRATED/ESTIMATED VIS SQUARED'
      END IF
      TITLE(3) = DATE
      IF ( FTYPE .EQ. 3 ) THEN
          WRITE(TITLE(4),1500) ICHAN, FILTERID(ICHAN), IBASE,
     $                         FTYPE, ORDER, FITRMS, CHISQ
      ELSE IF ( FTYPE .EQ. 4 ) THEN
          WRITE(TITLE(4),1500) ICHAN, FILTERID(ICHAN), IBASE,
     $                         FTYPE, ORDER, FITRMS, CHISQ
      ELSE IF ( FTYPE .EQ. 5 ) THEN
          WRITE(TITLE(4),1501) ICHAN, FILTERID(ICHAN), IBASE,
     $                         FTYPE, EXPO, FITRMS, CHISQ
      ELSE
          WRITE(TITLE(4),1500) ICHAN, FILTERID(ICHAN), IBASE, FTYPE,
     $                         MIN(ISMOOTH,999), FITRMS, CHISQ
      END IF
      WRITE(TITLE(5),1510) CMEAN, CRMS, CRANGE
      IF ( FTYPE .EQ. 4 ) THEN
          WRITE(TITLE(6),1515) (COEF(I),I=1,3), '  WEIGHTED FIT '
      ELSE
          WRITE(TITLE(6),1515) (COEF(I),I=1,3), '  EQUAL WT FIT '
      END IF
C
C Determine the range of the data, and use PGENV to set the range of the
C axes and to draw a box.  PGLABEL is used to label it.
C
      IF ( CALTYPE .EQ. 1 ) THEN
          XMIN = TCAL0
          XMAX = TCAL0 + TCAL1
      ELSE IF ( CALTYPE .EQ. 2 ) THEN
          XMIN = HCAL0
          XMAX = HCAL0 + HCAL1
      ELSE IF ( CALTYPE .EQ. 3 ) THEN
          XMIN = ZCAL0
          XMAX = ZCAL0 + ZCAL1
      ELSE IF ( CALTYPE .EQ. 4 ) THEN
          XMIN = MCAL0
          XMAX = MCAL0 + MCAL1
      ELSE IF ( CALTYPE .EQ. 5 ) THEN
          XMIN = FCAL0
          XMAX = FCAL0 + FCAL1
      ELSE IF ( CALTYPE .EQ. 6 ) THEN
          XMIN = NCAL0
          XMAX = NCAL0 + NCAL1
      ELSE
          WRITE(6,*) ' BOZO CITY.   YOU CANNOT GET HERE '
          STOP
      END IF
      IF ( (RTYPE .EQ. 1 ) .AND. (NDATA.GT.0) ) THEN
          YMIN = CMIN
          YMAX = CMAX
          DO J = 1, NDATA
              YMIN = MIN ( YMIN, YDATA(J) )
              YMAX = MAX ( YMAX, YDATA(J) )
          END DO
      ELSE IF ( RTYPE .EQ. 2 ) THEN
          YMIN = 0.
          YMAX = 1.1
      END IF
      CALL PGSCI ( LIGHT_BLUE )
C
C                            Set big letters so that PGENV will leave
C                            more room around the plot.
      CALL PGSCH ( 1.5 )
      CALL PGRNGE ( XMIN, XMAX, XLO, XHI )
      CALL PGRNGE ( YMIN, YMAX, YLO, YHI )
      CALL PGENV  ( XLO, XHI, YLO, YHI, 0, 1 )
C
C       In order, TITLE contains  X-axis label
C                                 Y-axis label
C                          The  big   top label.
C                          The second top label.
C                          The right hand label.
C                          Information line on bottom.
C
      CALL PGMTEXT ( 'T', 2.5, 0.5, 0.5, TITLE(3) )
      CALL PGSCH ( 1.0 )
      CALL PGMTEXT ( 'B', 3.5, 0.5, 0.5, TITLE(1) )
      CALL PGMTEXT ( 'L', 3.5, 0.5, 0.5, TITLE(2) )
      CALL PGMTEXT ( 'T', 1.5, 0.5, 0.5, TITLE(4) )
      CALL PGMTEXT ( 'R', 2.0, 0.5, 0.5, TITLE(5) )
      CALL PGMTEXT ( 'R', 4.0, 0.5, 0.5, TITLE(6) )
C
C    Write the star names with the colors in which they will be plotted.
C
      DO J = 1, NS
          WRITE(STARLIST(1:5),'(I5)') IPLOT(J)
          DISP = FLOAT(J)/FLOAT(NS+1)
C
C  Color 1 = white is reserved for the fitted line
C
          ICOLOR = MOD ( J+1, 15 ) + 1
          CALL PGSCI ( ICOLOR )
          CALL PGMTEXT ( 'B', 5.0, DISP, 0.5, STARLIST(1:5) )
      END DO
C
      IF ( NDATA .EQ. 0 ) GO TO 800
C
C   Plot the fitted line.
C
      IF ( FTYPE .NE. 0 ) THEN
          CALL PGSCI ( WHITE )
          CALL PGSLS ( 1 )
          CALL PGLINE ( NFIT, TSYS, VSYS )
      END IF
C
C    Plot the points, one star at a time
C
      DO 780 J = 1, NS
          NPTS = 0
          DO I = 1, NDATA
              IF ( STAR(IMARK(I)) .EQ. IPLOT(J) ) THEN
                  NPTS = NPTS + 1
                  XPLOT(NPTS) = XDATA(I)
                  YPLOT(NPTS) = YDATA(I)
                  ZPLOT(NPTS) = ERROR(I)
              END IF
          END DO
          IF ( NPTS .EQ. 0 ) GO TO 780
          ICOLOR = MOD( J + 1, 15 ) + 1
          CALL PGSCI( ICOLOR ) 
          IF ( LTYPE .EQ. 1 ) THEN
              CALL PGPOINT( NPTS, XPLOT, YPLOT, J+1 )
          ELSE IF ( LTYPE .EQ. 2 ) THEN
              CALL PGPOINT( NPTS, XPLOT, YPLOT, J+1 )
              CALL PGLINE(NPTS, XPLOT, YPLOT )
          ELSE IF ( LTYPE .EQ. 3 ) THEN
              CALL PGLINE(NPTS, XPLOT, YPLOT )
              DO I = 1, NPTS
                  VECX(1) = XPLOT(I)
                  VECX(2) = XPLOT(I)
                  VECY(1) = YPLOT(I) - ZPLOT(I)
                  VECY(2) = YPLOT(I) + ZPLOT(I)
                  CALL PGLINE(2, VECX, VECY )
              END DO
          END IF
  780 CONTINUE
  800 CONTINUE
C-----------------------------------------------------------------------
C   Prompt the user for what to do next.
C
      IF ( HCOPY ) THEN
          HCOPY = .FALSE. 
          CALL PGEND
          GO TO 145
      END IF
      CALL EGA_RESTORE_DEFAULT
  805 CONTINUE
      CALL LOCATE(0,0)
      IF ( NDATA .EQ. 0 ) THEN
          CALL WRITE_STRING ( 'No data.  Command:' )
      ELSE
          CALL WRITE_STRING ( 'Command:' )
      END IF
      CALL PGCURSE( RX, RY, CH )
      CALL CAPS ( CH )
      IF ( CH .EQ. '?' ) THEN
          CALL LOCATE(  0, 1 )
          CALL WRITE_STRING ( ' ?   Prints this help     ' )
          CALL LOCATE(  0, 2 )
          CALL WRITE_STRING ( ' A   Add a star           ' )
          CALL LOCATE(  0, 3 )
          CALL WRITE_STRING ( ' D   Delete a star        ' )
          CALL LOCATE(  0, 4 )
          CALL WRITE_STRING ( ' I   Identify a point     ' )
          CALL LOCATE(  0, 5 )
          CALL WRITE_STRING ( ' C   Toggle the channel   ' )
          CALL LOCATE(  0, 6 )
          CALL WRITE_STRING ( ' Z   Delete a point       ' )
          CALL LOCATE(  0, 7 )
          CALL WRITE_STRING ( ' R   Range of y axis      ' )
          CALL LOCATE(  0, 8 )
          CALL WRITE_STRING ( ' X   Toggle X axis type   ' )
          CALL LOCATE(  0, 9 )
          CALL WRITE_STRING ( ' Y   Toggle Y axis type   ' )
          CALL LOCATE(  0,10 )
          CALL WRITE_STRING ( ' L   Toggle line type     ' )
          CALL LOCATE(  0,11 )
          CALL WRITE_STRING ( ' S   Smooth the fit       ' )
          CALL LOCATE(  0,12 )
          CALL WRITE_STRING ( ' H   Make a hard copy     ' )
          CALL LOCATE(  0,13 )
          CALL WRITE_STRING ( ' K   Keep this chan cal   ' )
          CALL LOCATE(  0,14 )
          CALL WRITE_STRING ( ' F   Toggle type of fit   ' )
          CALL LOCATE(  0,15 )
          CALL WRITE_STRING ( ' G   Toggle seeing type   ' )
          CALL LOCATE(  0,16 )
          CALL WRITE_STRING ( ' M   Minimum error for fit' )
          CALL LOCATE(  0,17 )
          CALL WRITE_STRING ( ' E   Quit                 ' )
          CALL LOCATE(  0,18 )      
          CALL WRITE_STRING ( ' Q   Quit                 ' )
          CALL LOCATE(  0,19 )
          GO TO 805
C-----------------------------------------------------------------------
      ELSE IF ( CH .EQ. 'A' ) THEN
C
C    Add star to list.
C
          CALL LOCATE(0,1)
          L = 0
          STARLIST = ' '
          DO 808 K = 1, NSTAR
              DO 807 J = 1, NS
                  IF ( IPLOT(J) .EQ. SLIST(K) ) GO TO 808
  807         CONTINUE
              WRITE(STARLIST(5*L:5*(L+1)),'(I5)') SLIST(K)
              L = L + 1
              IF ( L .EQ. 15 ) THEN
                  CALL WRITE_STRING ( STARLIST )
                  CALL LOCATE(0,2)
                  L = 0
              END IF
  808     CONTINUE
          IF ( L .NE. 0 ) CALL WRITE_STRING( STARLIST(1:5*(L+1)) )
          CALL LOCATE (0,0)
          CALL WRITE_STRING ( 'Add star number:' )
          READ (5,*) ISTAR
          DO 810 J = 1, NSTAR
              IF ( ISTAR .NE. SLIST(J) ) GO TO 810
              NS = NS + 1
              NCALSTAR    = NS
              IPLOT(NS)   = ISTAR
              CALSTAR(NS) = ISTAR
              GO TO 815
  810     CONTINUE
  815     CONTINUE
C-----------------------------------------------------------------------
      ELSE IF ( CH .EQ. 'D' ) THEN
C
C   Delete star
C
          CALL LOCATE(0,0)
          DO 818 L = 1, NS, 15
              WRITE(STARLIST,'(15I5)') (IPLOT(J),J=L,MIN(NS,L+14))
              CALL WRITE_STRING ( STARLIST )
              CALL LOCATE(0,1+L/15)
  818     CONTINUE
          CALL WRITE_STRING ( 'Delete star number:' )
          READ (5,*) ISTAR
          J = 0
          DO 820 I = 1, NS
              IF ( IPLOT(I) .EQ. ISTAR ) J = I
  820     CONTINUE
          IF ( J .NE. 0 ) THEN
              DO 825 I = J+1, NS
                  IPLOT(I-1)   = IPLOT(I)
                  CALSTAR(I-1) = IPLOT(I)
  825         CONTINUE
              NS = NS - 1
              NCALSTAR = NS
              DO 830 I = 1, ISCAN
                  IF ( STAR(I) .EQ. ISTAR ) THEN
                      BAD(I) = 0
                  END IF
  830         CONTINUE
          END IF
C-----------------------------------------------------------------------
C    Toggle the Y-axis Scale format
C                          1   fixed scale, 0 to 1.2
C                          2   auto scale
C                          3   user input scale
C    
      ELSE IF ( CH .EQ. 'R' ) THEN
          IF ( (RTYPE .EQ. 2) .OR. (RTYPE .EQ. 3) ) THEN
              RTYPE = 3
              WRITE(6,*) ' INPUT Y AXIS RANGE '
              READ (5,*,END=100) YMIN,YMAX
              IF ( (YMIN.EQ.0.0) .AND. (YMAX.EQ.0.0) ) THEN
                  RTYPE = MOD(RTYPE,3) + 1
              END IF
          ELSE
              RTYPE = MOD(RTYPE,3) + 1
          END IF
C-----------------------------------------------------------------------
C   Toggle x axis type
C
      ELSE IF ( CH .EQ. 'X' ) THEN
          CALL LOCATE(  0, 1 )
          CALL WRITE_STRING   ( ' ENTER CALIBRATION TYPE: '//
     $        '1 TIME,  2 SEEING,  3 ZENITH ANGLE')
          CALL LOCATE(  0, 2 )
          CALL WRITE_STRING
     $      ( '4 MIRROR ANGLE,  5 FASN.,  6 N/INT' )
          CALL LOCATE(  0, 3 )
          READ ( 5,*,IOSTAT=IERR ) I
          IF ( IERR .NE. 0 ) I = 1
          IF ( ( I .LT. 1 ) .OR. ( I .GT. 6 ) ) THEN
              CALL WRITE_STRING ( ' BAD CALIBRATION TYPE ' )
          ELSE
              CALTYPE = I
          END IF
C-----------------------------------------------------------------------
      ELSE IF ( CH .EQ. 'Y' ) THEN
          YTYPE = MOD( YTYPE, 4 ) + 1
C-----------------------------------------------------------------------
      ELSE IF ( CH .EQ. 'L' ) THEN
          LTYPE = MOD(LTYPE,3) + 1
C-----------------------------------------------------------------------
      ELSE IF ( CH .EQ. 'F' ) THEN
C
          CALL LOCATE(  0, 1 )
          CALL WRITE_STRING
     $         ( ' ENTER FIT TYPE  0 NONE,   1 BOXCAR,    2 FFT     ' )
          CALL LOCATE(  17, 2 )
          CALL WRITE_STRING     ( '3 POLY,   4 EVEN POLY, 5 EXPONENT' )
          CALL LOCATE(  0, 3 )
          READ ( 5,*,IOSTAT=IERR ) I
          IF ( IERR .NE. 0 ) I = 0
          IF ( ( I .LT. 0 ) .OR. ( I .GT. 5 ) ) THEN
              CALL WRITE_STRING ( ' BAD FIT TYPE ' )
          ELSE
              FTYPE = I
          END IF
C-----------------------------------------------------------------------
C   Keep the calibration for this channel / type
C   NCALSCAN is the number of calibration scans for each channel.
C
      ELSE IF ( CH .EQ. 'K' ) THEN
C
C  Save the number of parameters in the fit (misnamed NFREE).
C
          IF ( ( FTYPE .EQ. 1 ) .OR. ( FTYPE .EQ. 2 ) ) THEN
              R1 = (CAL1-CAL0) * FLOAT(ISMOOTH) / 780.
              NFREE(ICHAN,CALTYPE) = ( XDATA(NDATA) - XDATA(1) ) / R1
          ELSE IF ( FTYPE .EQ. 3 ) THEN
              NFREE(ICHAN,CALTYPE) = ORDER + 1
          ELSE IF ( FTYPE .EQ. 4 ) THEN
              NFREE(ICHAN,CALTYPE) = ORDER/2 + 1
          ELSE IF ( FTYPE .EQ. 5 ) THEN
              NFREE(ICHAN,CALTYPE) = 2
          END IF
          NCALSCAN(ICHAN) = NDATA
C
C  Save the appropriate calibration.
C
          IF ( CALTYPE .EQ. 1 ) THEN
              DO I = 1, NFIT
                  TCAL(I,ICHAN) = VSYS(I)
              END DO
          ELSE IF ( CALTYPE .EQ. 2 ) THEN
C
C   For the seeing calibration, also keep the coefficients
C   of the seeing fit and the type of data being fit.
C
              DO I = 1, NFIT
                  HCAL(I,ICHAN) = VSYS(I)
              END DO
              IF ( FTYPE .EQ. 5 ) THEN
                  SEEING_COEF(1,ICHAN) = COEF(1)
                  SEEING_COEF(2,ICHAN) = COEF(2)
                  SEEING_COEF(3,ICHAN) = EXPO
              ELSE IF ( FTYPE .EQ. 4 ) THEN
                  IF ( ORDER .EQ. 2 ) THEN
                      SEEING_COEF(1,ICHAN) = COEF(1)
                      SEEING_COEF(2,ICHAN) = - COEF(2) / COEF(1)
                      SEEING_COEF(3,ICHAN) = ORDER
                  ELSE
                      SEEING_COEF(1,ICHAN) = 0.
                      SEEING_COEF(2,ICHAN) = 0.
                      SEEING_COEF(3,ICHAN) = ORDER
                  END IF
              ELSE
                  SEEING_COEF(1,ICHAN) = 0.
                  SEEING_COEF(2,ICHAN) = 0.
                  SEEING_COEF(3,ICHAN) = 0.
              END IF
              SEEING_TYPE = STYPE
          ELSE IF ( CALTYPE .EQ. 3 ) THEN
C
C    For zenith angle calibration, Record calibration and coefficients.
C
              DO I = 1, NFIT
                  ZCAL(I,ICHAN) = VSYS(I)
              END DO
              IF ( FTYPE .EQ. 4 ) THEN
                  ZD_COEF(0,ICHAN) = ORDER
                  ZD_COEF(1,ICHAN) = COEF(1)
                  ZD_COEF(2,ICHAN) = COEF(2)
                  ZD_COEF(3,ICHAN) = COEF(3)
              END IF
          ELSE IF ( CALTYPE .EQ. 4 ) THEN
              DO I = 1, NFIT
                  MCAL(I,ICHAN) = VSYS(I)
              END DO
          ELSE IF ( CALTYPE .EQ. 5 ) THEN
              DO I = 1, NFIT
                  FCAL(I,ICHAN) = VSYS(I)
              END DO
          ELSE IF ( CALTYPE .EQ. 6 ) THEN
              DO I = 1, NFIT
                  NCAL(I,ICHAN) = VSYS(I)
              END DO
          END IF
          CAL_CHISQ(ICHAN) = CHISQ
          CAL_RMS(ICHAN) = FITRMS
C-----------------------------------------------------------------------
C    Identify the star closest to the cursor
C
      ELSE IF ( CH .EQ. 'I' ) THEN
          CALL LOCATE (0,0)
C          WRITE(6,*) ' CURSOR POSITION = ', RX, RY
          INEAR = 1
          SCALE = (YMIN-YMAX)/(XMIN-XMAX)
          SCALE = SCALE * SCALE
          DISTANCE = SCALE*(RX-XDATA(1))*(RX-XDATA(1)) +
     $                     (RY-YDATA(1))*(RY-YDATA(1))
          DO 850 I = 1, NDATA
              D2 = SCALE*(RX-XDATA(I))*(RX-XDATA(I)) +
     $                   (RY-YDATA(I))*(RY-YDATA(I))
              IF ( D2 .LT. DISTANCE ) THEN
                  DISTANCE = D2
                  INEAR = I
              END IF
  850     CONTINUE
          CALL LOCATE ( 0,2 )
          WRITE(STARLIST,'(A10,I5,60X)')
     $                   ' STAR IS ', STAR(IMARK(INEAR))
          CALL WRITE_STRING ( STARLIST )
          CALL LOCATE (0,0)
          RX = XDATA(INEAR)
          RY = YDATA(INEAR)
          GO TO 805
C-----------------------------------------------------------------------
      ELSE IF ( CH .EQ. 'Z' ) THEN
          INEAR = 1
          SCALE = (YMIN-YMAX)/(XMIN-XMAX)
          SCALE = SCALE * SCALE
          DISTANCE = SQRT( SCALE*(RX-XDATA(1))*(RX-XDATA(1)) +
     $                           (RY-YDATA(1))*(RY-YDATA(1)) )
          DO 860 I = 1, NDATA
              D2 = SQRT( SCALE*(RX-XDATA(I))*(RX-XDATA(I)) +
     $                         (RY-YDATA(I))*(RY-YDATA(I)) )
              IF ( D2 .LT. DISTANCE ) THEN
                  DISTANCE = D2
                  INEAR = I
              END IF
  860     CONTINUE
          RX = XDATA(INEAR)
          RY = YDATA(INEAR)
          BAD(IMARK(INEAR)) = IOR( BAD(IMARK(INEAR)), MASK(ICHAN) )
C-----------------------------------------------------------------------
C  Toggle channel
C
      ELSE IF ( CH .EQ. 'C' ) THEN
          ICHAN = MOD(ICHAN,4) + 1
      ELSE IF ( CH .EQ. 'S' ) THEN
C
C  Smooth averge visibility
C
          CALL LOCATE (0,0)
          IF ( ( FTYPE .EQ. 3 ) .OR. ( FTYPE .EQ. 4 ) ) THEN
              CALL WRITE_STRING ( ' Order of fit:' )
              READ (5,*) ORDER
              ORDER = MIN ( MAX(ORDER,0), 10 )
          ELSE IF ( FTYPE .EQ. 5 ) THEN
              CALL WRITE_STRING ( ' Exponent of fit:' )
              READ (5,*) EXPO
          ELSE
              CALL WRITE_STRING ( 'Smoothing interval, minutes:')
              CALL LOCATE (29,0)
              READ (5,*) ISMOOTH
              ISMOOTH = MAX(ISMOOTH,0)
          END IF
C-----------------------------------------------------------------------
C   HARDCOPY
C
      ELSE IF ( CH .EQ. 'H' ) THEN
          HCOPY = .TRUE.
          CALL PGEND
          GO TO 145
C-----------------------------------------------------------------------
C  Toggle seeing type
C
      ELSE IF ( CH .EQ. 'G' ) THEN
          STYPE = MOD(STYPE, 2 ) + 1
C-----------------------------------------------------------------------
C   Input minimum allowed error 
C
      ELSE IF ( CH .EQ. 'M' ) THEN
          WRITE(6,1600) MINERR
          READ (5,*) MINERR
          DO J = 1, 4
              IF ( MINERR(J) .LT. .001 ) MINERR(J) = .001
          END DO
 1600 FORMAT ( ' Input minimum allowed error, all 4 channels '
     $        /' current values are ', 4F8.2 )
C-----------------------------------------------------------------------
      ELSE IF ( (CH.EQ.'Q') .OR. (CH.EQ.'E') ) THEN
C-----------------------------------------------------------------------
C   Finish up and save averages
C
          GO TO 870
      ELSE
C-----------------------------------------------------------------------
C  Bad command
C
          CALL LOCATE (0,1)
          CALL WRITE_STRING(' BAD COMMAND' )
          GO TO 805
      END IF
      CALL PGASK(.FALSE.)
      CALL PGADVANCE
      GO TO 150
  870 CONTINUE
      CALL PGEND
  900 CONTINUE
 1500 FORMAT ( 'CH=', I1, ' FILT=', A7, ' B=',I2, 2X,
     $        ' FIT=', I1, '/', I3,   '/', F5.3, 3X, 'CHI SQ =', F6.2 )
 1501 FORMAT ( 'CH=', I1, ' FILT=', A7, ' B=',I2, 2X,
     $        ' FIT=', I1, '/', F4.2, '/', F5.3, 3X, ' CHI SQ =', F6.2)
 1510 FORMAT ( 'MEAN V =', F6.3, ' RMS =' F6.3, ' FIT RNG =',F6.3, 1X )
 1515 FORMAT ( ' COEFFICIENTS = ', 3F8.3, A )
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE SMOOTH ( Y, N, ISMOOTH )
      IMPLICIT UNDEFINED (A-Z)
      SAVE
      INTEGER *4 N, ISMOOTH, I, J, NAVER
      REAL *4 Y(*), TEMP(1000)
      IF ( N .GT. 1000) THEN
          WRITE(6,*) ' Array is too big to smooth '
          GO TO 900
      END IF
      DO 120 I = 1, N
          TEMP(I) = 0.
          NAVER   = 0
          DO 110 J = MAX(I-ISMOOTH,1), MIN(I+ISMOOTH,N)
              IF ( Y(J) .EQ. 0. ) GO TO 110 
              TEMP(I) = TEMP(I) + Y(J)
              NAVER = NAVER + 1
  110     CONTINUE
          IF ( NAVER .GT. 0 ) THEN
              TEMP(I) = TEMP(I) / FLOAT( NAVER )
          END IF
  120 CONTINUE
      DO 150 I = 1, N
          Y(I) = TEMP(I)
  150 CONTINUE
  900 CONTINUE
      RETURN
      END
      SUBROUTINE MAKECAL (XDATA, YDATA, YFIT, NDATA, TSYS, VSYS,
     $                           N, ISMOOTH)
C
C  Make a calibration curve.  The input data is in XDATA, YDATA.  YFIT 
C  is the calibrated visibility for each observation.  The calibration
C  is also computed as VSYS, on a uniform grid, TSYS of length N.
C  ISMOOTH is the amount of smoothing.
C
C  The algorithm used is MODIFIED from numerical recipes.  It is an FFT
C  of the input data (ignoring the non-uniform data spacing), apodizes
C  it and transforms it back.  Then the data is interpolated onto a
C  uniform grid.
C
      IMPLICIT UNDEFINED (A-Z)
      SAVE
      REAL *4 XDATA(*), YDATA(*), YFIT(*), TSYS(*), VSYS(*), R1, RS
      INTEGER *4 NDATA, N, ISMOOTH, I, I1, I2, IDATA, J, K
      INTEGER *4 IFIRST, ILAST, ICNT
      REAL * 4 BOX
C
C  Interpolate the input data onto the grid TSYS.
C  Assume that the grig TSYS is uniform.
C
      BOX = 0.5 * ( TSYS(2) - TSYS(1) )
      IFIRST = N
      ILAST  = 1
      DO 120 J = 1, N
          VSYS(J) = 0.
          ICNT = 0
          DO 110 I = 1, NDATA
              IF ( ABS(XDATA(I)-TSYS(J)) .GE. BOX ) GO TO 110
              VSYS(J) = VSYS(J) + YDATA(I)
              ICNT = ICNT + 1
  110     CONTINUE
          IF ( ICNT .GT. 0 ) THEN
              VSYS(J) = VSYS(J) / FLOAT(ICNT)
              IFIRST = MIN(IFIRST,J)
              ILAST  = MAX(ILAST, J)
          ELSE
              VSYS(J) = -1.
          END IF
  120 CONTINUE
C
C   Linear interpolation for any missing data.
C
      J = IFIRST - 1
      DO 140 I = IFIRST, ILAST
          IF ( VSYS(I) .LT. 0. ) GO TO 140
          IF ( J .NE. I-1 ) THEN
              DO 130 K = J+1, I-1
                  VSYS(K) = VSYS(J) * ( TSYS(I)-TSYS(K) )
     $                    + VSYS(I) * ( TSYS(K)-TSYS(J) )
                  VSYS(K) = VSYS(K) / ( TSYS(I)-TSYS(J) )
  130         CONTINUE
          END IF
          J = I
  140 CONTINUE
C
C     Convert ISMOOTH from minutes to pixels
C
      RS = FLOAT( MAX(1, MIN(ISMOOTH, 780) ) ) / FLOAT(780)
      RS = 2 * N * RS
C
C     Smooth just the part of the curve that contains data
C
      I = ILAST - IFIRST + 1
      CALL SMOOFT( VSYS(IFIRST), I, RS )
C
C     Extrapolate the smooth curve by using the nearest point.
C
      DO 150 I = 1, IFIRST
          VSYS(I) = VSYS(IFIRST)
  150 CONTINUE
      DO 160 I = ILAST, N
          VSYS(I) = VSYS(ILAST)
  160 CONTINUE
C
C     Calculate the fitted visibility for each observation.
C
      BOX = TSYS(2) - TSYS(1)
      IF ( BOX .LE. 0. ) THEN
          DO I = 1, NDATA
              YFIT(I) = 1.
          END DO
      ELSE
          DO I = 1, NDATA
              RS = ( XDATA(I) - TSYS(1) ) / BOX
              I1 = RS
              RS = RS - FLOAT(I1)
              IF ( I1 .LT. 0 ) THEN
                  YFIT(I) = VSYS(1)
              ELSE IF ( I1 .LE. N-2 ) THEN
                  YFIT(I) = (1.-RS)*VSYS(I1+1) + RS * VSYS(I1+2)
              ELSE 
                  YFIT(I) = VSYS(N)
              END IF
          END DO
      END IF
C
      RETURN
      END
      REAL *4 FUNCTION GAUSS( X, D )
      IMPLICIT UNDEFINED (A-Z)
      SAVE
      REAL *4 X, D, ARG
      ARG = X*X/(0.5*D*D)
      IF ( ARG .LT. .1 ) THEN
          GAUSS = 1. - ARG * ( 1. - 0.5*ARG ) 
      ELSE IF ( ARG .GT. 5 ) THEN
          GAUSS = 0.
      ELSE
          GAUSS = EXP(-ARG)
      END IF
      RETURN
      END
      SUBROUTINE GAUSS_SMOOTH ( X, Y, YFIT, N, WIDTH,
     $                          XLINE, YLINE, NLINE, RMS )
C
C  Smooth the data using a truncated gaussian.
C
C  Input data                                 X(i), Y(i), i=1, N
C  Values fitted to the data                  YFIT(i),    i=1, N
C  1 sigma half width of smoothing (hours)    WIDTH
C
C  To facilitate plotting the function, the value of the fit, YLINE(i),
C  is evaluated on a grid of x values XLINE(i), i=1, NLINE input by
C  the user.
C
C  The quality of the fit is given by RMS, the rms of the difference
C  between the data and the fit.
C  
      IMPLICIT UNDEFINED (A-Z)
      SAVE
      INTEGER *4 N, NLINE
      REAL *4    X(*), Y(*), XLINE(*), YLINE(*), YFIT(*)
      REAL *4 RMS, WIDTH, MEAN, DELTA, WEIGHT, WT
      REAL *4 START_X, STOP_X, NSIG
      INTEGER *4 I, J, I0, I1, NPTS
      REAL *4 GAUSS
      EXTERNAL GAUSS
C.......................................................................
      DATA NSIG / 3.0 /     ! half width for truncation in sigma.
C.......................................................................
C    Convert full width in minutes to half width in hours.
      DELTA = WIDTH
C.......................................................................
C  I0 is the first point in range, I1 is the first point beyond range.
C
      I0 = 1
      I1 = 1
      DO I = 1, N
          START_X = X(I) - NSIG * DELTA
          STOP_X  = X(I) + NSIG * DELTA
          DO WHILE ( ( I0.LT.N ) .AND. (X(I0) .LT. START_X) )
              I0 = I0 + 1
          END DO
          DO WHILE ( (I1.LT.N+1) .AND. (X(I1).LT.STOP_X) )
              I1 = I1 + 1
          END DO
          NPTS = I1 - I0
          YFIT(I) = 0.
          WEIGHT = 0.
          IF ( NPTS .GT. 0 ) THEN
              DO J = I0, I1-1
                  WT = GAUSS( X(I)-X(J), DELTA )
                  YFIT(I) =  YFIT(I) + WT * Y(J)
                  WEIGHT  = WEIGHT + WT
              END DO
              YFIT(I) = YFIT(I) / WEIGHT
          END IF
      END DO
C.......................................................................
C  Evaluate the smoothed curve on a grid of points.
C  Extapolate data by assuming no change beyond the observations.
C
      I0 = 1
      I1 = 1
      DO I = 1, NLINE
          IF ( XLINE(I) .EQ. 0.0 ) THEN
              YLINE(I) = 0.
          ELSE IF ( XLINE(I) .LT. X(1) ) THEN
              YLINE(I) = YFIT(1)
          ELSE IF ( XLINE(I) .GT. X(N) ) THEN
              YLINE(I) = YFIT(N)
          ELSE
              START_X = XLINE(I) - NSIG * DELTA
              STOP_X  = XLINE(I) + NSIG * DELTA
              DO WHILE ( (I0.LT.N) .AND. (X(I0).LT.START_X) )
                  I0 = I0 + 1
              END DO
              DO WHILE ( (I1.LE.N) .AND. (X(I1).LT.STOP_X) )
                  I1 = I1 + 1
              END DO
              NPTS = I1 - I0
              YLINE(I) = 0.
              WEIGHT = 0.
              IF ( NPTS .GT. 0 ) THEN
                  DO J = I0, I1-1
                      WT = GAUSS( XLINE(I)-X(J), DELTA )
                      YLINE(I) =  YLINE(I) + WT * Y(J)
                      WEIGHT = WEIGHT + WT
                  END DO
                  IF ( WEIGHT .GT. .1 ) THEN
                      YLINE(I) = YLINE(I) / WEIGHT
                  ELSE
                      YLINE(I) = 0.
                  END IF
              END IF
          END IF
      END DO
 1200 FORMAT ( I5, 2F10.4, 4X, F10.4 )
C.......................................................................
C  Do a linear interpolation for missing data
C
      I0 = 1
      DO WHILE ( I0 .LT. NLINE-1 )
          IF ( YLINE(I0+1) .GT. 0.0 ) THEN
              I0 = I0 + 1
          ELSE
              I1 = I0
              DO WHILE ( YLINE(I1+1) .EQ. 0.0 )
                  I1 = I1 + 1
              END DO
              DO J = I0+1, I1-1
                  YLINE(J) =    YLINE(I0) + (YLINE(I1)-YLINE(I0))
     $               * (XLINE(J)-XLINE(I0))/(XLINE(I1)-XLINE(I0))
              END DO
              I0 = I1 + 1
          END IF
      END DO
C.......................................................................
      RETURN
      END
      SUBROUTINE EXPO_FIT ( X, Y, ERROR, YFIT, N, EXPO,
     $                      XLINE, YLINE, NLINE, COEF )
C
C  Do a weighted fit of the equation   Y = A + B x^expo
C  to a data set by varying the coefficients A and B.
C
C  Input data                 X(i), Y(i), i=1, N
C  Values fitted to the data  YFIT(i), i=1, N
C  Exponent of the fit        EXPO
C  coefficients of fit        COEF
C
C  This subroutine output the coefficients of the fit.
C  To facilitate plotting the function, it also outputs
C  the value of the fit, YLINE(i), on a grid of x values
C  XLINE(i), i=1, NLINE input by the user.
C
      IMPLICIT UNDEFINED (A-Z)
      INTEGER *4 N, NLINE, MAXPAR, I, NUNK
      PARAMETER (MAXPAR=5)
      REAL *4 X(*), Y(*), ERROR(*), XLINE(*), YLINE(*), YFIT(*)
      REAL *4 EXPO, COEF(2), RELERR, MINERR, WT
      REAL *8 A(MAXPAR,MAXPAR), AA(MAXPAR,MAXPAR), WEIGHT
      REAL *8 CC(MAXPAR), DD(MAXPAR), EE(MAXPAR), BB, VV
      SAVE
      DATA  RELERR / .01 /
C.......................................................................
      NUNK = 2
      IF ( NUNK+2 .GT. MAXPAR ) THEN
          WRITE(0,*) ' TOO many parameters.  Limit is ',MAXPAR
          NUNK = MAXPAR-2
      END IF
      IF ( ABS(EXPO) .LT. .1 ) THEN
          NUNK = 1.0
          DD(2) = 0.
          WRITE(6,*) ' EXPONENT IS TOO CLOSE TO ZERO   ',
     $               '   ONLY ONE COEF FIT '
      END IF
C               LS mode=1, initialize for least squares fit
      WEIGHT = 1.
      CALL LS(1,NUNK,BB,CC,DD,EE,VV, WEIGHT, A, AA, MAXPAR)
      DO I = 1, N
          BB    = Y(I)
          CC(1) = 1.D0
          CC(2) = X(I)**EXPO
          MINERR = RELERR * Y(I)
          WEIGHT = 1. / ( MAX(ERROR(I),MINERR) )**2
C
C              LS mode = 2, increment coefficient matrix for least sq fit.
C                 mode = 3, increment solution vector    for least sq fit.
C                           The weight is included as real*8
C
          CALL LS(2,NUNK,BB,CC,DD,EE,VV, WEIGHT, A, AA, MAXPAR )
          CALL LS(3,NUNK,BB,CC,DD,EE,VV, WEIGHT, A, AA, MAXPAR )
      END DO
C
C          LS mode = 4, do the least squares fit
C
      IF ( N .LT. 5 ) THEN
          COEF(1) = 0.
          COEF(2) = 0.
      ELSE 
          CALL LS(4, NUNK, BB, CC, DD, EE, VV, WEIGHT, A, AA, MAXPAR)
          COEF(1) = DD(1)
          COEF(2) = DD(2)
      END IF
C
C   Convert the coefficients so that they are of the form
C
C    V^2 = COEF(1) * ( 1. - COEF(2)*X^EXPO )
C
      COEF(2) = - COEF(2) / COEF(1)
C
      DO I = 1, N
          YFIT(I) = COEF(1) * ( 1. - COEF(2)*X(I)**EXPO )
      END DO
      DO I = 1, NLINE
          YLINE(I) = COEF(1) * ( 1. - COEF(2)*XLINE(I)**EXPO )
      END DO
      RETURN
      END