C @(#)redsubs2.for	19.1 (ES0-DMD) 02/25/03 13:28:47
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
      SUBROUTINE POLYGON(ARRAY,IFLAG,LABEL)
C  @(#)redsubs2.for	19.1  (ESO-IPG)  02/25/03  13:28:47
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT (c) European Southern Observatory 1992
C.IDENT		redsubs2.for
C.MODULE
C.AUTHOR	Andrew T. Young
C.KEYWORD
C.LANGUAGE	FORTRAN 77
C.PURPOSE       more subroutines for REDUCE
C.COMMENTS
C.VERSION	5.9
C.RETURNS
C.ENVIRONMENT
C.
C-----------------------------------------------------------------------------
C
C
C       scans array for items with NSTR=iflag, and does polygon interpolation.
C
C
      IMPLICIT NONE
      REAL ARRAY, XBAR, YBAR, SLOPE, FRACT, YMAX, YMIN
C
      INTEGER IFLAG, LEN, LWORD, JBGN, NIGHT, K, J,
     1         JMIN, JMAX, JEND, K1, K2
C
      INCLUDE 'MID_REL_INCL:obs.inc'
C
C
      DIMENSION ARRAY(MXOBS)
C
      CHARACTER A, LABEL*(*), CARD*79
C
      LOGICAL HELP
      EXTERNAL HELP
C
C
C
      LEN=LWORD(LABEL)
      JBGN=1
C
      DO 100 NIGHT=1,NIGHTS
C
        CALL SPACE2
        CALL NEED(24)
        CARD(:LEN)=LABEL(:LEN)
        CARD(LEN+1:LEN+6)='  for '
        CALL JD2DAT(REAL(2400000.5D0+DTZERO(NIGHT)),CARD(LEN+7:))
        CALL TV(CARD)
        K=0
        YMAX=-3.E33
        YMIN=3.E33
C
        DO 10 J=JBGN,NDATA
C
C         find data for current night.
          IF (NITE(J).NE.NIGHT) GO TO 20
          IF (NSTR(J).EQ.IFLAG)THEN
C           add to XS, YS.
            K=K+1
            XS(K)=DJOBS(J)
            YS(K)=ARRAY(J)
            YMAX=MAX(YMAX,YS(K))
            YMIN=MIN(YMIN,YS(K))
            IF (K.EQ.1) JMIN=J
            JMAX=J
          END IF
C
   10   CONTINUE
        J=NDATA+1
C
   20   CONTINUE
        JEND=J-1
C
C    JBGN, JEND index the first and last data of ANY kind for the night.
C    JMIN, JMAX index the first and last IFLAG data for the night.
C
C
C       K is count of data for night N.
C
        IF (YMAX.GT.YMIN) THEN
C           K > 1.
   22       CALL NEED(24)
            CALL PLOT(K,XS,YS,'*')
            CALL TVN('     Time (decimal days) -->')
            CALL ASKN('OK to interpolate with a polygon?',A)
            IF (A.EQ.'N') THEN
   24           CALL TV('Enter C to replace data with a Constant,')
                IF (K.GE.6)
     1               CALL TVN('      L to smooth with a Linear fit,')
                CALL ASKN('   or Q to Quit.',A)
                IF (A.EQ.'C') THEN
C                   be lazy: fit line, but force it horizontal.
                    CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE)
                    SLOPE=0.
                ELSE IF (A.EQ.'L') THEN
                    IF (K.LT.6) THEN
                        CALL TV('Not enough data to fit a line safely.')
                        CALL TVN('Try a constant instead.')
                        GO TO 24
                    END IF
                    CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE)
                ELSE
                   GO TO 24
                END IF
C
                DO 28 J=JBGN,JEND
                    ARRAY(J)=(DJOBS(J)-XBAR)*SLOPE+YBAR
   28           CONTINUE
C
                GO TO 41
            ELSE IF (A.EQ.'Y' .OR. A.EQ.'O') THEN
C               go on.
            ELSE IF (HELP(A)) THEN
                CALL TV('Reply NO to smooth data.')
                GO TO 22
            ELSE
                GO TO 22
            END IF
        ELSE IF (YMAX.EQ.YMIN) THEN
            CALL TV('Only 1 value available for this night!')
            CALL SPACE2
            DO 30 J=JBGN,JEND
                ARRAY(J)=ARRAY(JMIN)
   30       CONTINUE
            GO TO 99
        ELSE
            CALL TV('NO DATA AVAILABLE FOR THIS NIGHT')
            CALL STSEPI
        END IF
C
C       Now interpolate.
C
        K1=1
        K2=2
C
        DO 40 J=JBGN,JEND
           IF (J.LT.JMIN) THEN
               ARRAY(J)=ARRAY(JMIN)
           ELSE IF (J.GT.JMAX) THEN
               ARRAY(J)=ARRAY(JMAX)
           ELSE
               IF (DJOBS(J).GT.XS(K2)) THEN
                   K1=K2
                   K2=K2+1
               END IF
C
C               linear interp...
               FRACT=(DJOBS(J)-XS(K1))/(XS(K2)-XS(K1))
               ARRAY(J)=YS(K1)+(YS(K2)-YS(K1))*FRACT
           END IF
   40   CONTINUE
C
   41   CARD(66:78)=CARD(53:65)
        CARD(53:65)=CARD(40:52)
        CARD(40:52)=CARD(27:39)
        CARD(27:39)=CARD(14:26)
        CARD(14:26)=CARD(:13)
        CARD(:13)='Interpolated '
        CALL NEED(24)
        CALL TV(CARD)
C       fudge to get right plot limits....
        XS(1)=DJOBS(JBGN)
        XS(K)=DJOBS(JEND)
        CALL PLOT(-K,XS,YS,' ')
        DO 50 J=JBGN,JEND
        CALL PLOT(-1,DJOBS(J),ARRAY(J),'+')
   50   CONTINUE
C       now undo fudge.
        XS(1)=DJOBS(JMIN)
        XS(K)=DJOBS(JMAX)
        CALL PLOT(K,XS,YS,'*')
        CALL RTNCON('     Time (decimal days) -->',28)
C
   99   JBGN=JEND+1
C       IF (JBGN.GT.NDATA) RETURN
  100 CONTINUE
C
      RETURN
      END
      SUBROUTINE DARKFIT(NDET,NDUSED,LABEL)
C
C       scans SIGNAL for items with NSTR=-1, models dark level vs. DTMP,
C       and subtracts dark current for detector number NDET.
C
C
C     C A U T I O N :  This routine has *** NOT *** been debugged!
C     ===============
C
      IMPLICIT NONE
C
      REAL XLIMS, YLIMS, SIGMIN, SIGMAX, YMAX, 
     1   YMIN, XBAR, YBAR, SLOPE, XBAR1, YBAR1, SLOPE1, TRY, XBAR2, 
     2         YBAR2, SLOPE2, FIT
      INTEGER NDET, LEN, LWORD, K, J, JMIN,
     1   JMAX, K3, K1, I, K2, NIX
C
      INCLUDE 'MID_REL_INCL:mbands.inc'
C     PARAMETER (MBANDS=9)
      INTEGER MXCOLR
      PARAMETER (MXCOLR=3*MBANDS)
      INTEGER NDUSED(MXCOLR)
C
      INCLUDE 'MID_REL_INCL:obs.inc'
C     This defines SIGNAL, NSTR, KBND, etc.
C
C     Includes for dlsq:
      INCLUDE 'MID_REL_INCL:kingfit.inc'
      INCLUDE 'MID_REL_INCL:dlsq.inc'
C
      CHARACTER A1, LABEL*16, CARD*79
C
      DIMENSION XLIMS(2),YLIMS(2)
C
      LOGICAL HELP
      EXTERNAL HELP
C
      EXTERNAL YP2EXP
C
C
C       Begin execution:
C
      LEN=LWORD(LABEL)
C
    2 K=0
      SIGMAX=-3.E33
      SIGMIN=3.E33
      CALL SPACE2
      CARD(:LEN)=LABEL(:LEN)
      CARD(LEN+1:)='  Log(DARK) vs. Temperature'
C
      DO 10 J=1,NDATA
C
C       find  dark data    for  current detector.
        IF (NSTR(J).EQ.-1 .AND. KBND(J).EQ.NDET)THEN
C           add to XS, YS.
            K=K+1
            XS(K)=DTMP(J)
            IF (SIGNAL(J).GT.0.) THEN
               YS(K)=LOG(SIGNAL(J))
            ELSE
               CARD='Negative dark reading at MJD ='
               WRITE(CARD(35:),'(F12.5)') DTZERO(NITE(J))+DJOBS(J)
               CARD(48:)='in channel'
               WRITE(CARD(55:),'(I3)') NDET
               CALL TV(CARD)
               K=K-1
               GO TO 10
            END IF
            SIGMAX=MAX(SIGMAX,SIGNAL(J))
            SIGMIN=MIN(SIGMIN,SIGNAL(J))
C
C           SYMB is night label:
            IF (NITE(J).LT.10) THEN
                SYMB(K)=CHAR(ICHAR('0')+NITE(J))
            ELSE
                SYMB(K)=CHAR(ICHAR('A')+NITE(J)-1)
            END IF
C
            IF (K.EQ.1) THEN
               JMIN=J
            END IF
C               do this every time, in case only 1 datum.
            JMAX=J
        END IF
C
   10 CONTINUE
C
C       K is count of data for current detector.
C
   21 YMAX=-3.E33
      YMIN=3.E33
C
      IF (SIGMAX.GT.SIGMIN) THEN
C           K > 1.
C           set page size:
            CALL PLOT(0,78.,22.,'P')
C           set different-symbol mode:
   22       CALL PLOT(0,78.,22.,'D')
            CALL NEED(24)
            CALL TV(CARD)
C           plot each night with its own symbol.
            CALL PLOT(K,XS,YS,SYMB)
C           set same-symbol mode:
            CALL PLOT(0,78.,22.,'S')
            IF (K.GT.10) THEN
                CALL TVN('     Temperature -->')
                CALL ASKN('OK to fit with 2 lines?',A1)
            ELSE IF (HELP(A1)) THEN
                CALL TV('Reply YES only if concave upward')
                CALL TVN('(high at ends, lower in middle)')
                GO TO 22
            ELSE
                A1='N'
                CALL RTNCON('     Temperature -->', 20)
            END IF
C
            IF (A1.EQ.'N') THEN
   24           CALL TV('Enter C to replace data with a Constant,')
                CALL TVN('      L to smooth with a Linear fit,')
                CALL ASKN('   or Q to Quit.',A1)
C
                IF (A1.EQ.'C') THEN
C                   be lazy: fit line, but force it horizontal.
                    CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE)
                    SLOPE=0.
                    NIX=1
                    IX(1)=2
                ELSE IF (A1.EQ.'L') THEN
                    IF (K.LT.6) THEN
                        CALL TV('Not enough data to fit 1 line safely.')
                        CALL TVN('Try a constant instead.')
                        GO TO 24
                    END IF
C                       fit a simple exponential in temperature:
                    CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE)
                    NIX=0
                ELSE
                   GO TO 24
                END IF
C
C
C               Refine one-exponential fit:
C
                IPR=0
                KIP=2
                PG(1)=EXP(YBAR)
                PG(2)=SLOPE
                DO 26 J=1,K
                  X(1,J)=XS(J)-XBAR
                  Y(J)=EXP(YS(J))
C                 use expected ordinates to weight:
                  W(J)=1./EXP(YBAR+X(1,J)*SLOPE)
   26           CONTINUE
C               refine fit:
C      subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
                CALL DLSQ(YP2EXP,K,2,   1,    2,  NIX,  1,    0,  1.E-3)
                SLOPE=LOG(P(2))
                YBAR=P(1)
C
C
C               now subtract fit from ALL data with this detector:
C
                DO 28 J=1,NDATA
                   IF ((NSTR(J).GT.0 .AND. NDUSED(KBND(J)).NE.NDET) .OR.
     1                 (NSTR(J).EQ.-1 .AND. KBND(J).NE.NDET)) GO TO 28
C                  subtract dark:
                   SIGNAL(J)=SIGNAL(J)-EXP((DTMP(J)-XBAR)*SLOPE+YBAR)
                   IF (NSTR(J).EQ.-1) THEN
                      YMAX=MAX(YMAX,SIGNAL(J))
                      YMIN=MIN(YMIN,SIGNAL(J))
                   END IF
   28           CONTINUE
C
                GO TO 61
C
C
            ELSE IF (A1.EQ.'Y' .OR. A1.EQ.'O') THEN
C               go on to do two-exponential fit at *****.
C
            ELSE IF (A1.EQ.'H')THEN
                CALL TV('If the plot is generally concave upward,')
                CALL TVN('and the different nights (different symbols)')
                CALL TVN('appear to follow the same curve, reply YES.')
                CALL TV('If the plot is straight or concave downward,')
                CALL TVN('or some nights are offset, reply NO.')
                CALL TV('Overlapping points are denoted by "$".')
                CALL SPACE2
                CALL RTNCON(' ',1)
                GO TO 21
C
            ELSE
                GO TO 22
C
            END IF
C
      ELSE IF (SIGMAX.EQ.SIGMIN) THEN
            CALL TV('Only 1 dark value available!')
            CALL SPACE2
            DO 30 J=1,NDATA
C               subtract dark:
                IF (NSTR(J).GT.0 .AND. NDUSED(KBND(J)).EQ.NDET)
     1             SIGNAL(J)=SIGNAL(J)-SIGMIN
   30       CONTINUE
C           skip plot.
            GO TO 99
      ELSE
            CALL TV('NO DATA AVAILABLE -- cannot subtract dark')
            RETURN
      END IF
C
C     *****
C
C       Here to fit 2 lines (exponentials), and remove their sum.
C
      K3=K/3
C
      IF (K3.LT.9) THEN
          CALL TV('CAUTION: probably too few points to get a fit.')
          CALL ASKN('Were there any wild points in the plot above?',A1)
C
          IF (A1.EQ.'Y')THEN
              CALL TV('Then you should try to fit just one line.')
              CALL RTNCON(' ',1)
              GO TO 21
          ELSE IF (A1.EQ.'N') THEN
C               go ahead.
          ELSE IF (A1.EQ.'H') THEN
C               get help: show plot again.
                GO TO 21
          ELSE
          END IF
C
      END IF
C
      CALL SORT2(XS,YS,K)
C     fit trial line to bottom 1/3 of data.
      CALL ROBLIN(XS,YS,K3, XBAR1,YBAR1,SLOPE1)
C
C     remove lower trial line, and fit the upper line.
C
      K1=0
C
      DO 32 I=1,K3
C        move upper third to lower third:
         K2=K-K3+I
         TRY=EXP(YS(K2))-EXP(YBAR1+SLOPE1*(XS(K2)-XBAR1))
         IF (TRY.LE.0.) THEN
C           skip this point.
         ELSE
            K1=K1+1
C           save K1 in K1+K.
            XS(K1+K)=XS(K1)
            YS(K1+K)=YS(K1)
C           move K2 to K1.
            YS(K1)=LOG(TRY)
            XS(K1)=XS(K2)
         END IF
   32 CONTINUE
C
      IF (K1.LT.6)THEN
          CALL TV('Sorry; can''t fit 2 lines.  Try 1 line.')
          CALL RTNCON(' ',1)
          GO TO 2
      ELSE
          CALL ROBLIN(XS,YS,K1, XBAR2,YBAR2,SLOPE2)
          IF (SLOPE2.LT.0. .OR. SLOPE2.LT.SLOPE1) THEN
C               unphysical result; complain.
              CALL TV('Sorry; attempted fit is nonsense.')
              CALL RTNCON('Try again, using 1 line.',24)
              GO TO 2
          ELSE
          END IF
      END IF
C
C     Now restore saved data to lower K1:
      DO 35 I=1,K1
         XS(K1)=XS(K1+K)
         YS(K1)=YS(K1+K)
   35 CONTINUE
C
C     *****  Polish the result.  *****
C
C     Set up for DLSQ:
C
      IPR=0
      NIX=0
C     Set up data:
      XBAR=(XBAR1+XBAR2)*0.5
      DO 36 I=1,K
         X(1,I)=XS(I)-XBAR
         Y(I)=YS(I)
         W(I)=1./(EXP(YBAR1+SLOPE1*(XS(I)-XBAR1)) +
     1            EXP(YBAR2+SLOPE2*(XS(I)-XBAR2)) )
   36 CONTINUE
      KIP=4
      PG(1)=EXP(YBAR1+SLOPE1*(XBAR-XBAR1))
      PG(2)=SLOPE1
      PG(3)=EXP(YBAR2+SLOPE2*(XBAR-XBAR2))
      PG(4)=SLOPE2
C               refine fit:
      DO 40 I=1,3
C     Subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
        CALL DLSQ(YP2EXP, K,    2,    1,    2,   NIX, 1,    0,   1.E-3)
C
        IF (P(1).LE.0.D0 .OR. P(3).LE.0.D0) THEN
          CALL TV('Sorry -- unable to fit two lines.  Try one.')
          CALL RTNCON(' ',1)
          GO TO 2
        END IF
C       Update values:
        DO 38 J=1,4
   38   PG(J)=P(J)
        SLOPE1=P(2)
        YBAR1=LOG(P(1))-SLOPE1*(XBAR-XBAR1)
        SLOPE2=P(4)
        YBAR2=LOG(P(3))-SLOPE2*(XBAR-XBAR2)
C       Revise weights:
        DO 39 J=1,K
          W(J)=1./(EXP(YBAR1+SLOPE1*(XS(I)-XBAR1)) +
     1             EXP(YBAR2+SLOPE2*(XS(I)-XBAR2)) )
   39   CONTINUE
   40 CONTINUE
C
C
C
C     Subtract the interpolated values from ALL data for this det.
C
      DO 60 J=1,NDATA
         IF ((NSTR(J).GT.0 .AND. NDUSED(KBND(J)).NE.NDET) .OR.
     1            (NSTR(J).EQ.-1 .AND. KBND(J).NE.NDET)) GO TO 60
         FIT=EXP(YBAR1+SLOPE1*(DTMP(J)-XBAR1)) +
     1        EXP(YBAR2+SLOPE2*(DTMP(J)-XBAR2))
         SIGNAL(J)=SIGNAL(J)-FIT
         IF (NSTR(J).NE.-1) GO TO 60
C        use only dark data to set limits for residual plot.
         YMAX=MAX(YMAX,SIGNAL(J))
         YMIN=MIN(YMIN,SIGNAL(J))
   60 CONTINUE
C
C       Now inspect resids. vs. time:
C
   61 CARD(:LEN)=LABEL(:LEN)
      CARD(LEN+1:)=' Dark RESIDUALS  vs. Time'
      CALL NEED(24)
      CALL TV(CARD)
C     set plot limits....
      XLIMS(1)=DJOBS(1)
      XLIMS(2)=DJOBS(NDATA)
      YLIMS(1)=YMIN
      YLIMS(2)=YMAX
      CALL PLOT(0,XLIMS,YLIMS,'L')
      DO 80 J=1,NDATA
          IF (NSTR(J).NE.-1 .OR. KBND(J).NE.NDET) GO TO 80
          IF (NITE(J).LT.10) THEN
             A1=CHAR(ICHAR('0')+NITE(J))
          ELSE
             A1=CHAR(ICHAR('A')+NITE(J)-1)
          END IF
          CALL PLOT(-1,DJOBS(J),SIGNAL(J),A1)
   80 CONTINUE
C     force plot:
      CALL PLOT(1,0.,0.,'+')
      CALL RTNCON('     Time (decimal days) -->',28)
C
   99 CONTINUE
C
      RETURN
      END
      SUBROUTINE YP2EXP
C
C   Fits  Y = P(1)*EXP(P(2)*Z(1)) + P(3)*EXP(P(4)*Z(1))
C
C
      IMPLICIT NONE
C
C     Includes for dlsq:
      INCLUDE 'MID_REL_INCL:kingfit.inc'
      INCLUDE 'MID_REL_INCL:dlsq.inc'
C
      DOUBLE PRECISION TERM
C
C
C   Start with first term:
C
C
C     KIP is set by DARKFIT.
C
      IP(1)=1
      IP(2)=2
C
      TERM=P(2)*Z(1)
      TERM=EXP(MIN(TERM,30.D0))
      PART(1)=TERM
      PART(2)=P(1)*TERM*Z(1)
C
      YT = P(1)*TERM
C
      IF (KIP.EQ.2) RETURN
C
C   Here for second term:
C
      IP(3)=3
      IP(4)=4
C
      TERM=P(4)*Z(1)
      TERM=EXP(MIN(TERM,30.D0))
      PART(3)=TERM
      PART(4)=P(3)*TERM*Z(1)
C
      YT = YT + P(3)*TERM
C
C
      RETURN
      END
      SUBROUTINE DARKPOLY(NDET,NDUSED,LABEL)
C
C       scans SIGNAL for items with NSTR=-1, and does polygon interpolation
C       and subtraction of dark current for detector number NDET.
C
C
      IMPLICIT NONE
C
      REAL XLIMS, YLIMS, SIGMIN, SIGMAX, SIGA, SIGZ, YMAX, 
     1         YMIN, XBAR, YBAR, SLOPE, FRACT
      INTEGER NDET, LEN, LWORD, JBGN, NIGHT, K, J,
     1         JMIN, JMAX, JEND, K1, K2, NIX
C
      INCLUDE 'MID_REL_INCL:mbands.inc'
C     PARAMETER (MBANDS=9)
      INTEGER MXCOLR
      PARAMETER (MXCOLR=3*MBANDS)
      INTEGER NDUSED(MXCOLR)
C
      INCLUDE 'MID_REL_INCL:obs.inc'
C     This defines NITE, SIGNAL, NSTR, KBND, etc.
C
C     Includes for dlsq:
      INCLUDE 'MID_REL_INCL:kingfit.inc'
      INCLUDE 'MID_REL_INCL:dlsq.inc'
C
      CHARACTER A1, LABEL*16, CARD*79
C
      DIMENSION XLIMS(2),YLIMS(2)
C
      LOGICAL HELP
      EXTERNAL HELP
C
      EXTERNAL YPLIN
C
C
C       Begin execution:
C
      LEN=LWORD(LABEL)
      JBGN=1
C
      DO 100 NIGHT=1,NIGHTS
C
        CALL SPACE2
        CALL NEED(24)
        CARD(:LEN)=LABEL(:LEN)
        CARD(LEN+1:LEN+10)='  DARK on '
        CALL JD2DAT(REAL(2400000.5D0+DTZERO(NIGHT)),CARD(LEN+11:))
        CARD(LWORD(CARD)+3:)='=  MJD'
        WRITE (CARD(LWORD(CARD)+1:),'(F7.0)') DTZERO(NIGHT)
        CALL TV(CARD)
        SIGMAX=-3.E33
        SIGMIN=3.E33
        K=0
C
        DO 10 J=JBGN,NDATA
C
C         find dark data for current night and detector.
          IF (NITE(J).NE.NIGHT) GO TO 20
c         KBND holds NDET for DARK data, remember.
          IF (NSTR(J).EQ.-1 .AND. KBND(J).EQ.NDET)THEN
C           add to XS, YS.
            K=K+1
            XS(K)=DJOBS(J)
            YS(K)=SIGNAL(J)
            IF (K.EQ.1) JMIN=J
            IF (K.EQ.1) SIGA=SIGNAL(J)
            JMAX=J
            SIGZ=SIGNAL(J)
            SIGMAX=MAX(SIGMAX,SIGNAL(J))
            SIGMIN=MIN(SIGMIN,SIGNAL(J))
          END IF
C
   10   CONTINUE
        J=NDATA+1
C
   20   CONTINUE
        JEND=J-1
C
C    JBGN, JEND index the first and last data of ANY kind for the night.
C    JMIN, JMAX index the first and last DARK data for the night.
C    SIGA, SIGZ are the first and last dark LEVELS for the night.
C    SIGMIN, SIGMAX are the lowest and highest dark LEVELS for the night.
C
C
        YMAX=-3.E33
        YMIN=3.E33
C
C       K is count of data for night N.
C
        IF (SIGMAX.GT.SIGMIN) THEN
   22       CALL NEED(24)
            CALL PLOT(K,XS,YS,'*')
            CALL TVN('     TIME (decimal days) -->')
            CALL ASKN('OK to interpolate with a polygon?',A1)
            IF (A1.EQ.'N') THEN
   24           CALL TV('Enter C to replace data with a Constant,')
                IF (K.GE.6)
     1               CALL TVN('      L to smooth with a Linear fit,')
                CALL ASKN('   or Q to Quit.',A1)
                IF (A1.EQ.'C') THEN
C                   be lazy: fit line, but force it horizontal.
                    CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE)
                    SLOPE=0.
                    NIX=1
                    IX(1)=2
                ELSE IF (A1.EQ.'L') THEN
                    IF (K.LT.6) THEN
                        CALL TV('Not enough data to fit a line safely.')
                        CALL TVN('Try a constant instead.')
                        GO TO 24
                    END IF
                    CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE)
                    NIX=0
                ELSE IF (HELP(A1)) THEN
                    CALL TV('Reply NO to smooth data.')
                    GO TO 22
                ELSE
                   GO TO 24
                END IF
C
C
C               Refine fit:
C
                IPR=0
                PG(1)=YBAR
                PG(2)=SLOPE
                DO 26 J=1,K
                  X(1,J)=XS(J)-XBAR
                  Y(J)=YS(J)
                  W(J)=1.
   26           CONTINUE
C               refine fit:
C     Subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
                CALL DLSQ(YPLIN,K,2,  1,    2,   NIX, 1,    0,   1.E-3)
                SLOPE=P(2)
                YBAR=P(1)
C
C
C               Subtract fit:
C
                DO 28 J=JBGN,JEND
                    IF((NSTR(J).GT.0 .AND. NDUSED(KBND(J)).NE.NDET) .OR.
C                       real star, not observed with current det...
     1                 (NSTR(J).LE.0 .AND. KBND(J).NE.NDET)) GO TO 28
C                       ...or dark, not observed with current det.
C                   OK, subtract interpolated dark line:
                    SIGNAL(J)=SIGNAL(J)-((DJOBS(J)-XBAR)*SLOPE+YBAR)
                    IF (NSTR(J).EQ.-1) THEN
C                       reset limits for residual plot.
                        YMAX=MAX(YMAX,SIGNAL(J))
                        YMIN=MIN(YMIN,SIGNAL(J))
                    END IF
   28           CONTINUE
C
C               plot residuals.
C
                CALL NEED(24)
                CARD(66:78)=CARD(53:65)
                CARD(53:65)=CARD(40:52)
                CARD(40:52)=CARD(27:39)
                CARD(27:39)=CARD(14:26)
                CARD(14:26)=CARD(:13)
                CARD(:13)='RESIDUALS of '
                CALL TV(CARD)
C               set plot limits....
                XLIMS(1)=DJOBS(JBGN)
                XLIMS(2)=DJOBS(JEND)
                YLIMS(1)=YMIN
                YLIMS(2)=YMAX
                CALL PLOT(0,XLIMS,YLIMS,'L')
                CALL XAXIS(XLIMS)
                DO 29 J=JBGN,JEND
                   IF (NSTR(J).NE.-1 .OR. KBND(J).NE.NDET) GO TO 29
C                  plot only dark resids., not star data.
                   CALL PLOT(-1,DJOBS(J),SIGNAL(J),'*')
   29           CONTINUE
C               force plot:
                CALL PLOT(1,0.,0.,'+')
                CALL RTNCON('     Time (decimal days) -->',28)
                GO TO 99
C
            ELSE IF (A1.EQ.'Y' .OR. A1.EQ.'O') THEN
C               go on.
            ELSE
                GO TO 22
            END IF
        ELSE IF (K.EQ.1) THEN
            CALL TV('Only 1 datum available for this night!')
            CALL SPACE2
            DO 30 J=JBGN,JEND
C                subtract dark:
                 IF((NSTR(J).GT.0 .AND. NDUSED(KBND(J)).NE.NDET) .OR.
     1                 (NSTR(J).LE.0 .AND. KBND(J).NE.NDET)) THEN
c                   skip it.
                 ELSE
                    SIGNAL(J)=SIGNAL(J)-SIGA
                 END IF
   30       CONTINUE
C           skip plot.
            GO TO 99
        ELSE
            CALLTV('NO DATA AVAILABLE FOR THIS NIGHT -- no subtraction')
            GO TO 99
        END IF
C
C       Do polygon interpolation here.
C
        K1=1
        K2=2
C
        DO 40 J=JBGN,JEND
           IF (NDUSED(KBND(J)).NE.NDET) GO TO 40
C          skip observations with other detectors.
           IF (J.LT.JMIN) THEN
C              Before 1st dark (at JMIN); use it.
               SIGNAL(J)=SIGNAL(J)-SIGA
           ELSE IF (J.GT.JMAX) THEN
C              After last dark (at JMAX); use it.
               SIGNAL(J)=SIGNAL(J)-SIGZ
           ELSE
C              Can interpolate.
               IF (DJOBS(J).GT.XS(K2)) THEN
C                  Step to next interval.
                   K1=K2
                   K2=K2+1
               END IF
C
C               linear interp...
                FRACT=(DJOBS(J)-XS(K1))/(XS(K2)-XS(K1))
                SIGNAL(J)=SIGNAL(J)-(YS(K1)+(YS(K2)-YS(K1))*FRACT)
           END IF
C          note that we have no plot of residuals at the DARK data points,
C          as they are negligible for polygon fit.
   40   CONTINUE
C
   99   JBGN=JEND+1
C       IF (JBGN.GT.NDATA) RETURN
  100 CONTINUE
C
      RETURN
      END
      SUBROUTINE ROBLIN (X, Y, N,  XBAR, YBAR, SLOPE)
C
C       Fits a robust line to N pairs (X,Y), having medians (Xbar, Ybar)
C       and SLOPE.  Method: see Hoaglin, Mosteller & Tukey, Ch.5.
C
C       Does NOT mess up input data; but leaves them sorted on X in XDUM,YDUM.
C
C
C       Input arguments:
C
      INTEGER N
      REAL X(N), Y(N)
C
C       Output arguments:
C
      REAL XBAR,YBAR,SLOPE
C
C   Scratch arrays:
C
      INCLUDE 'MID_REL_INCL:obs.inc'
C     (defines MXOBS.)
      REAL XDUM,YDUM,ZDUM
      COMMON /SCRAT/ XDUM(MXOBS), YDUM(MXOBS), ZDUM(MXOBS)
C
C   Local variables:
C
      INTEGER I,J,N3,MID,MIDL,MIDR
      REAL XL,XR,YL,YR,B,B0,B1,B2,RL,RR,DELTA1,DELTAX,DR0,DR1
C
C
C   *****  BEGIN EXECUTION  *****
C
      IF (N.GT.MXOBS)THEN
            CALL TV('Too many data to sort in available slots.')
            CALL TV('Expand MXOBS and recompile.')
            CALL STSEPI
      ELSE IF (N.EQ.2)THEN
            XBAR=(X(1)+X(2))/2.
            YBAR=(Y(1)+Y(2))/2.
            IF (X(2).NE.X(1)) THEN
              SLOPE=(Y(2)-Y(1))/(X(2)-X(1))
            ELSE
              SLOPE=0.
            END IF
            RETURN
      ELSE IF (N.EQ.1)THEN
            XBAR=X(1)
            YBAR=Y(1)
            SLOPE=0.
            RETURN
      ELSE IF (N.LE.0)THEN
            CALL TV('ROBLIN can''t fit a line to 0 points!')
            CALL TVN('BUG in program')
            CALL STSEPI
      END IF
C
      DO 1 I=1,N
      XDUM(I)=X(I)
      YDUM(I)=Y(I)
      ZDUM(I)=Y(I)
    1 CONTINUE
C
      CALL SORT1(ZDUM,N)
      CALL SORT2(XDUM,YDUM,N)
      MID=(N+1)/2
C
      IF (MID.EQ.N/2) THEN
C                             N EVEN
            XBAR=(XDUM(MID)+XDUM(MID+1))*0.5
            YBAR=(ZDUM(MID)+ZDUM(MID+1))*0.5
      ELSE
C                             N ODD
            XBAR=XDUM(MID)
            YBAR=ZDUM(MID)
      END IF
C
C             Pick size of end thirds:
C
      N3=(N+1)/3
C                     Pick end thirds:
      MIDL=(N3+1)/2
      MIDR=N+1-MIDL
C
      IF (MIDL.EQ.N3/2) THEN
C                             N3 EVEN
            XL=(XDUM(MIDL)+XDUM(MIDL+1))*0.5
            XR=(XDUM(MIDR)+XDUM(MIDR-1))*0.5
      ELSE
C                             N3 ODD
            XL=XDUM(MIDL)
            XR=XDUM(MIDR)
      END IF
      DELTAX=(XR-XL)
C
C       Now find y medians for left and right groups.
C                             LEFT group:
      DO 10 I=1,N3
   10 ZDUM(I)=YDUM(I)
      CALL SORT1(ZDUM,N3)
C
      IF (MIDL.EQ.N3/2) THEN
C                     N3 EVEN
            YL=(ZDUM(MIDL)+ZDUM(MIDL+1))*0.5
      ELSE
C                     N3 ODD
            YL=ZDUM(MIDL)
      END IF
C                             RIGHT group:
      DO 20 I=1,N3
   20 ZDUM(I)=YDUM(N+1-I)
      CALL SORT1(ZDUM,N3)
C
      IF (MIDL.EQ.N3/2) THEN
C                     N3 EVEN
            YR=(ZDUM(MIDL)+ZDUM(MIDL+1))*0.5
      ELSE
C                     N3 ODD
            YR=ZDUM(MIDL)
      END IF
C
      B0=(YR-YL)/DELTAX
      B=B0
C
C       B0 is initial slope.
C
C                     *****  BEGIN ITERATIONS  *****
      DO 50 J=1,5
C             Put RESIDS. in zdum:
C                             LEFT:
        DO 30 I=1,N3
   30   ZDUM(I)=YDUM(I)-(YBAR+B*(XDUM(I)-XBAR))
        CALL SORT1(ZDUM,N3)
C
        IF (MIDL.EQ.N3/2) THEN
C                     N3 EVEN
            RL=(ZDUM(MIDL)+ZDUM(MIDL+1))*0.5
        ELSE
C                     N3 ODD
            RL=ZDUM(MIDL)
        END IF
C                             RIGHT:
        DO 40 I=1,N3
   40   ZDUM(I)=YDUM(N+1-I)-(YBAR+B*(XDUM(N+1-I)-XBAR))
        CALL SORT1(ZDUM,N3)
C
        IF (MIDL.EQ.N3/2) THEN
C                     N3 EVEN
            RR=(ZDUM(MIDL)+ZDUM(MIDL+1))*0.5
        ELSE
C                     N3 ODD
            RR=ZDUM(MIDL)
        END IF
C
        IF (J.EQ.1) THEN
            DR0=RR-RL
            DELTA1=DR0/DELTAX
            B1=B0+DELTA1
            IF (B1.EQ.B0) GOTO 100
            B=B1
        ELSE
            DR1=RR-RL
            IF (DR1.EQ.0.) GOTO 100
            IF (ABS((DR1-DR0)/DR1).LT.1.E-4) GOTO 100
            B2=B1-DR1*(B1-B0)/(DR1-DR0)
            B=B2
            IF (ABS(B2-B1).LT.1.E-6*ABS(B1)) GO TO 100
C                     prepare for next round:
            B0=B1
            B1=B2
            DR0=DR1
        END IF
   50 CONTINUE
C
  100 CONTINUE
      SLOPE=B
C       Current value of b is converged.  Recalculate all resids.:
C
      DO 110 I=1,N
  110 ZDUM(I)=YDUM(I)-(YBAR+B*(XDUM(I)-XBAR))
      CALL SORT1(ZDUM,N)
C
      IF (MID.EQ.N/2) THEN
C                             N EVEN
            YBAR=(ZDUM(MID)+ZDUM(MID+1))*0.5+YBAR
      ELSE
C                             N ODD
            YBAR=ZDUM(MID)+YBAR
      END IF
C
C
      RETURN
      END
      SUBROUTINE SORT1(X,N)
C
C       Sorts X, in place.
C
C  BOOTHROYD'S SHELL-SORT:  COMM.ACM 6,445(1963) ALG.201.
C
C  SEE H.LORIN,"SORTING & SORT SYSTEMS" (ADDISON-WESLEY,1975)
C  CH.3 & 11 FOR THEORY AND TESTS.
C
C
      IMPLICIT NONE
C
      INTEGER N
      REAL X(N)
C
C
C  Local variables:
C
      INTEGER I,J,M
      REAL DUM
C
C  **** BEGIN ****
C
      IF(N.EQ.1)RETURN
      I=1
    3 I=I+I
      IF(I.LE.N)GO TO 3
      M=I-1
    5 M=M/2
      IF(M.EQ.0)RETURN
      DO 20 J=1,N-M
      DO 10 I=J,1,-M
      IF(X(I+M).GE.X(I)) GO TO 20
      DUM=X(I)
      X(I)=X(I+M)
      X(I+M)=DUM
   10 CONTINUE
   20 CONTINUE
      GO TO 5
C
      END
      SUBROUTINE SORT2(X,Y,N)
C
C       Sorts (X,Y) pairs on X, in place.
C
C  BOOTHROYD'S SHELL-SORT:  COMM.ACM 6,445(1963) ALG.201.
C
C  SEE H.LORIN,"SORTING & SORT SYSTEMS" (ADDISON-WESLEY,1975)
C  CH.3 & 11 FOR THEORY AND TESTS.
C
C
      IMPLICIT NONE
C
      INTEGER N
      REAL X(N),Y(N)
C
C
C  Local variables:
C
      INTEGER I,J,M
      REAL DUM
C
C  **** BEGIN ****
C
      IF(N.EQ.1)RETURN
      I=1
    3 I=I+I
      IF(I.LE.N)GO TO 3
      M=I-1
    5 M=M/2
      IF(M.EQ.0)RETURN
      DO 20 J=1,N-M
      DO 10 I=J,1,-M
      IF(X(I+M).GE.X(I)) GO TO 20
      DUM=X(I)
      X(I)=X(I+M)
      X(I+M)=DUM
      DUM=Y(I)
      Y(I)=Y(I+M)
      Y(I+M)=DUM
   10 CONTINUE
   20 CONTINUE
      GO TO 5
C
      END
      SUBROUTINE YPLIN
C
C  YP for use by DLSQ in linear fits.
C
C
C  fits:   Y = P(1) + P(2)*Z(1)
C
C
      IMPLICIT NONE
C
C     Includes for dlsq:
      INCLUDE 'MID_REL_INCL:kingfit.inc'
      INCLUDE 'MID_REL_INCL:dlsq.inc'
C
C
C
      YT = P(1) + P(2)*Z(1)
C
      PART(1)=1.D0
      PART(2)=Z(1)
C
      IP(1)=1
      IP(2)=2
      KIP=2
C
C
      RETURN
      END
      SUBROUTINE SKYSUB(BANDS,KMIN,KMAX,DTYPE,KDIAPHRAGM,EXTIN)
C
C       scans SIGNAL for items with DTYPE='Y', and does modelling,
C       interpolation, and subtraction of sky for each passband.
C
C
      IMPLICIT NONE
C
C
      INCLUDE 'MID_REL_INCL:mbands.inc'
C     PARAMETER (MBANDS=9)
C
      INCLUDE 'MID_REL_INCL:obs.inc'
      REAL XDUM,YDUM,ZDUM
      COMMON /SCRAT/ XDUM(MXOBS), YDUM(MXOBS), ZDUM(MXOBS)
C
C
C   Args:
C
      INTEGER KMIN, KMAX, KDIAPHRAGM
      CHARACTER DTYPE(MXOBS)
      REAL EXTIN(MBANDS)
C
C
      REAL COSOB, SINOB, RASUN, DESUN, HASUN, SOLONG, ELMOON, BMOON, 
     1    SINPHI, COSPHI, ELONG, ELROT, UTROT, STUTZ, TWOPI, ST2UT, 
     2    TNOON1, TNOON2, PI, DEGRAD, ALAT, XLIMS, YLIMS, ALTS, T, 
     3    YMAX, YMIN, AIRMAX, SALTS, OBLIQ, SIGMIN, 
     4    SIGMAX, RISE, SET, DAWN, DUSK, UTRAD, TLST, HAMOON, COSDEL, 
     5    SINDEL, COSHA, XM, YM, ZM, AZMOON, DELAZ, COSEL, XBAR, 
     6    YBAR, SLOPE, OLDVAR, ELONGMAX, ELONGMIN, ZENITH
C
      INTEGER NB, LEN, LWORD, JBGN, NIGHT, NTWI,
     1       K, J, JMIN, JMAX, JEND, MOONUP, NERROR, NIX, NPTS, NFIT
C
C   COMMONS FOR SPHERICAL TRIG.:
C
      COMMON /CSUN/ COSOB,SINOB,RASUN,DESUN,HASUN,SOLONG,ELMOON,BMOON
      COMMON /SPHERE/ SINPHI,COSPHI,ELONG,ELROT,UTROT,STUTZ,TWOPI,ST2UT,
     1 TNOON1,TNOON2,PI,DEGRAD,ALAT
      INTEGER MXCOLR
      PARAMETER (MXCOLR=3*MBANDS)
C
C       commons for star catalog:
      INCLUDE 'MID_REL_INCL:mstars.inc'
C     PARAMETER (MSTARS=1650)
      CHARACTER *32 STARS
      COMMON /SCATA/ STARS(MSTARS)
C
C     Includes for dlsq:
      INCLUDE 'MID_REL_INCL:kingfit.inc'
      INCLUDE 'MID_REL_INCL:dlsq.inc'
C
C
C   Commons for 12 Diaphragms:
C
      INTEGER NP2NDIA(12*MBANDS)
      INTEGER NDIA2NP(12,MBANDS), NDIAS,NDIA
C          indices: NDIA,K
      REAL      XNDIA(12,MBANDS)
      COMMON /DIACOM/XNDIA,NP2NDIA,NDIA2NP,NDIAS,NDIA
C
      CHARACTER DIANAM(12)*4
      COMMON /DIACH/DIANAM
C
C
C   LOCAL variables:
C
      CHARACTER A1, CARD*80, BANDS(MXCOLR)*8, DATSTR*16
C
      REAL ERRSAVE(12)
      INTEGER LORDER(12),L,NN, JDUSK,JDAWN, NPARMS,NREFDIA
C
      LOGICAL TWILIT, SHORT
C
C
      EXTERNAL YPLIN, YPSKY
C
C
      DIMENSION XLIMS(2),YLIMS(2)
      DIMENSION ALTS(2), SALTS(5)
C
C       true altitude of Moon at rise & set:
      DATA ALTS/-.0145,-.0145/
C
C
C       true altitude of Sun at limit of astronomical twilight:
      DATA SALTS/-.309,0.,0.,0.,-.309/
C
C
C
C       Begin execution:
C
      CALL SPACE2
      CALL TV('               SKY subtraction:')
      CALL TV('You will see graphs of sky brightness, plotted against')
      CALL TVN('time, airmass, and lunar elongation, for each band,')
      CALL TVN('on each night separately.')
      CALL TV('Then you can choose the method of fitting or')
      CALL TVN('interpolating sky brightness, for each case.')
      CALL SPACE
      CALL RTNCON(' ',1)
      CALL PLOT(0,78.,23.,'P')
C
      DO 1000 NB=KMIN,KMAX
        LEN=LWORD(BANDS(NB))
        JBGN=1
      DO 1000 NIGHT=1,NIGHTS
        CALL SPACE2
        CALL NEED(24)
    2   CARD(:LEN)=BANDS(NB)
        CARD(LEN+1:LEN+10)='  SKY on'
        CALL JD2DAT(REAL(2400000.5D0+DTZERO(NIGHT)),CARD(LEN+11:))
        DATSTR=CARD(LEN+6:)
        SHORT=.FALSE.
C
    3   K=0
        DO 5 NDIA=1,NDIAS
           NDIA2NP(NDIA,1)=0
    5   CONTINUE
        YMAX=-3.E33
        YMIN=3.E33
        AIRMAX=1.
C
        DO 10 J=JBGN,NDATA
C
C         find range of sky data for current night and passband.
          IF (NITE(J).NE.NIGHT) GO TO 12
          IF (KBND(J).NE.NB) GO TO 10
          IF (DTYPE(J).EQ.'Y')THEN
C           add to XS, YS.
            K=K+1
            XS(K)=DJOBS(J)
            YS(K)=SIGNAL(J)
            YMAX=MAX(YMAX,SIGNAL(J))
            YMIN=MIN(YMIN,SIGNAL(J))
            AIRMAX=MAX(AIRMAX,AIRM(J))
            IF (K.EQ.1) JMIN=J
            IF (K.EQ.1) SIGMIN=SIGNAL(J)
            JMAX=J
            SIGMAX=SIGNAL(J)
          END IF
C
   10   CONTINUE
        J=NDATA+1
C
   12   CONTINUE
        JEND=J-1
C
C    JBGN, JEND index the first and last data of ANY kind for the night.
C    JMIN, JMAX index the first and last SKY data for the night & band.
C    SIGMIN, SIGMAX are the first and last sky LEVELS for the night.
C
C       K is count of data for night N.
C
        T=(DTZERO(NIGHT)+INT(DJOBS((JBGN+JEND)/2))-51544.5)/36525.
        CALL STUTZR(T)
C       STUTZR sets sidereal time for later use by UTMOON & UTSUN.
C       Set OBLIQ for SUN (called by UTSUN).
        OBLIQ=(23.43929-T*(1.30042E-2+T*(1.639E-7-T*5.036E-7)))*DEGRAD
        COSOB=COS(OBLIQ)
        SINOB=SIN(OBLIQ)
C
C     Adjust YMAX for possible outliers:
C
        CALL SORT2(YS,XS,K)
C       (Note reversal of Y and X.)
        YM=(YS((K+1)/2)+YS(K/2+1))*0.5
C       YM is median.
        NN=(K+2)/4
        ZM=YS(K-NN)-YS(NN)
C       ZM is interquartile range.
        YMAX=MIN(YMAX,YM+4.*ZM)
C
        IF (K.GT.1) THEN
C
C           Take shortcut if required:
            IF (SHORT) THEN
               A1='T'
               GO TO 102
            END IF
C
C           set up plots; first, TIME:
C
C           set plot limits:
            XLIMS(1)=DJOBS(JMIN)
            XLIMS(2)=DJOBS(JMAX)
            YLIMS(1)=YMIN
            YLIMS(2)=YMAX
            CALL PLOT(0,XLIMS,YLIMS,'L')
C
C           plot moonrise & moonset.  First, guarantee RISE < SET:
            RISE=DJOBS(JMIN)
            SET=DJOBS(JMAX)
C           UTMOON needs T in Julian centuries from J2000.
            T=((2400000.5D0-2451545.D0)+DTZERO(NIGHT)+
     1                   INT(DJOBS((JBGN+JEND)/2)))/36525.
            CALL UTMOON(T,ALTS,1,*21)
C           UTMOON returns result in UTROT, in /SPHERE/.
            RISE=UTROT+INT(DJOBS((JBGN+JEND)/2))
C           we have at most 18 intervals between YMIN & YMAX.
            CALL PLOT(-1,RISE,(YMIN+14.*YMAX)/15.,'|')
            CALL PLOT(-1,RISE,(YMIN*2.+13.*YMAX)/15.,'|')
            CALL PLOT(-1,RISE,YMAX,'R')
            IF (DJOBS(JMIN).LT.RISE .AND. DJOBS(JMAX).GT.RISE)
     1         CARD(LWORD(CARD)+5:)='R = moonrise'
   21       CALL UTMOON(T,ALTS,2,*22)
            SET=UTROT+INT(DJOBS((JBGN+JEND)/2))
            CALL PLOT(-1,SET,(YMIN+14.*YMAX)/15.,'|')
            CALL PLOT(-1,SET,(YMIN*2.+13.*YMAX)/15.,'|')
            CALL PLOT(-1,SET,YMAX,'S')
            IF (DJOBS(JMIN).LT.SET .AND. DJOBS(JMAX).GT.SET)
     1         CARD(LWORD(CARD)+5:)='S = moonset'
   22       CONTINUE
C
C           RISE & SET now are correct for Moon.  Do SUN:
            DUSK=DJOBS(JMIN)
            DAWN=DJOBS(JMAX)
C           Here the argument is DAYS from J2000.
            CALL UTSUN(T*36525.,SALTS,1,*23)
            DUSK=UTROT+INT(DJOBS((JBGN+JEND)/2))
   23       CALL UTSUN(T*36525.,SALTS,5,*24)
            DAWN=UTROT+INT(DJOBS((JBGN+JEND)/2))
   24       CONTINUE
C           DUSK & DAWN are now limits of astronomical twilight.
C
            TWILIT=.FALSE.
            JDUSK=JEND
            JDAWN=JBGN
            DO 26 J=JMIN,JMAX
                IF (KBND(J).NE.NB) GO TO 26
                IF ( (DAWN.LT.DUSK  .AND.
     1               DJOBS(J).GT.DAWN .AND. DJOBS(J).LT.DUSK) .OR.
     2               (DAWN .GT.DUSK  .AND.
     3               (DJOBS(J).GT.DAWN .OR. DJOBS(J).LT.DUSK) )  ) THEN
C                   Twilight.
                    A1='t'
                    TWILIT=.TRUE.
                ELSE
C                   Not twilit.
                    A1='*'
                    JDUSK=MIN(JDUSK,J)
                    JDAWN=MAX(JDAWN,J)
                END IF
                IF (DTYPE(J).EQ.'Y')THEN
                    CALL PLOT(-1,DJOBS(J),SIGNAL(J),A1)
                END IF
   26       CONTINUE
            IF (TWILIT) CARD(LWORD(CARD)+5:)='t = twilight sky'
            CALL TV(CARD)
            CALL SPACE
            CALL PLOT(1,0.,0.,'+')
            WRITE(CARD,'(A,F7.0)')'    TIME (MJD =',DTZERO(NIGHT)
            CARD(LWORD(CARD):)=' + decimal of day) -->'
            CALL RTNCON(CARD,LWORD(CARD))
C
C               End of TIME plot; begin AIRMASS plot.
C
            CALL NEED(24)
            CARD(:32)='SKY brightness vs. AIRMASS, for'
            CARD(33:LEN+33)=BANDS(NB)
            CARD(LEN+34:)=DATSTR
            CARD(LWORD(CARD)+5:)='M = Moon up'
            CALL TV(CARD)
            CALL SPACE
            XLIMS(1)=1.
            XLIMS(2)=AIRMAX
            CALL PLOT(0,XLIMS,YLIMS,'L')
            K=0
            MOONUP=0
            NTWI=0
            DO 30 J=JMIN,JMAX
              IF (KBND(J).NE.NB) GO TO 30
              IF (DTYPE(J).EQ.'Y')THEN
                K=K+1
                IF ( (DAWN.LT.DUSK  .AND.
     1               DJOBS(J).GT.DAWN .AND. DJOBS(J).LT.DUSK) .OR.
     2               (DAWN .GT.DUSK  .AND.
     3               (DJOBS(J).GT.DAWN .OR. DJOBS(J).LT.DUSK) )  ) THEN
C                   Twilight.
                    A1='t'
                    NTWI=NTWI+1
                    KX(1,K)=-1
                ELSE IF ( (RISE.LT.SET  .AND.
     1               DJOBS(J).GT.RISE .AND. DJOBS(J).LT.SET) .OR.
     2               (RISE .GT.SET  .AND.
     3               (DJOBS(J).GT.RISE .OR. DJOBS(J).LT.SET) )  ) THEN
C                   Moon is up.  Flag in KX.
                    A1='M'
                    MOONUP=MOONUP+1
C                   Flag for model solution.
                    KX(1,K)=1
                ELSE
C                   Moon is down.
                    A1='-'
                    KX(1,K)=0
                END IF
                SYMB(K)=A1
                CALL PLOT(-1,AIRM(J),SIGNAL(J),A1)
C               Set KX(2,*) to NDIA, or 0 if no diaphragm:
                DO 28 NDIA=1,NDIAS
                  IF (DIAPH(J).EQ.DIANAM(NDIA)) THEN
                     KX(2,K)=NDIA
                     NDIA2NP(NDIA,1)=NDIA2NP(NDIA,1)+1
                     GO TO 29
                  END IF
   28           CONTINUE
C               Here only if no dias.
                KX(2,K)=0
C                          Debugging:
C                          IF (KDIAPHRAGM.NE.-1) CALL TV('Dia.error!')
   29           CONTINUE
C               Set KX(3,*) to J for use in later plots.
                KX(3,K)=J
C               Set X(1,*) to airmass.
                X(1,K)=AIRM(J)
C               Set X(4,*) to nearest star signal.
                CALL HALO(NB,DTYPE,KDIAPHRAGM,J,JBGN,JEND, X(4,K))
                Y(K)=SIGNAL(J)
                W(K)=1./ABS(SIGNAL(J))
              END IF
   30       CONTINUE
            CALL PLOT(1,50.,-50.,' ')
            CALL RTNCON('     AIRMASS -->',16)
C
C               End of AIRMASS plot; begin LUNAR ELONGATION plot.
C
            IF(MOONUP.LT.2) GO TO 100
            ELONGMAX=0.
            ELONGMIN=50.
            CALL NEED(24)
            CARD(:41)='SKY brightness vs. LUNAR ELONGATION, for'
            CARD(42:LEN+42)=BANDS(NB)
            CARD(LEN+43:)=DATSTR
            CARD(LWORD(CARD)+5:)='M = Moon up'
            CALL TV(CARD)
            CALL SPACE
            XLIMS(1)=0.
            XLIMS(2)=180.
            CALL PLOT(0,XLIMS,YLIMS,'L')
            K=0
            DO 40 J=JMIN,JMAX
              IF (KBND(J).NE.NB) GO TO 40
              IF (DTYPE(J).EQ.'Y')THEN
                K=K+1
                T=(DTZERO(NIGHT)-51544.5D0+DJOBS(J))/36525.
                UTRAD=MOD(DJOBS(J),1.)*TWOPI
                TLST=STUTZ+(UTRAD*1.00273791)+ELONG
C               MOON uses T for lunar position, TLST for topocentric ephemeris.
                CALL MOON(T, TLST)
C               MOON returns lunar coords. in RASUN, DESUN in /CSUN/.
                HAMOON=TLST-RASUN
                COSDEL=COS(DESUN)
                SINDEL=SIN(DESUN)
                COSHA=COS(HAMOON)
                XM= -COSDEL*SIN(HAMOON)
                YM=COSPHI*SINDEL-SINPHI*COSDEL*COSHA
                ZM=SINPHI*SINDEL+COSPHI*COSDEL*COSHA
                IF (ZM.GT.ALTS(1)) THEN
C                  Set X(2,*) to LUNAR airmass.  (ZM is cosZ(Moon).)
                   X(2,K)=(.0096467+ZM*(.148386+ZM*1.0024342))/
     1                    (.000303978+ZM*(.0102963+ZM*(.149864+ZM)))
                ELSE
C                  Moon below horizon.
                   X(2,K)=50.
C                                  WRITE(CARD,*)'ZM=',ZM
C                                  IF(KX(1,K).EQ.1)CALL TV(CARD)
                   IF(KX(1,K).EQ.1) SYMB(K)='m'
                END IF
                AZMOON=ATAN2(YM,XM)
                DELAZ=AZMOON-AZIM(J)
C               approximate cosz(sky) as 1/airmass.
                COSEL=ZM/AIRM(J) +
     1               SQRT((1.-ZM*ZM)*(1.-1./AIRM(J)**2))*COS(DELAZ)
C               Set X(3,*) to lunar elongation.
                X(3,K)=ACOS(COSEL)
                IF (KX(1,K).EQ.1) THEN
                   ELONGMAX=MAX(ELONGMAX,X(3,K))
                   ELONGMIN=MIN(ELONGMIN,X(3,K))
                END IF
                XS(K)=ACOS(COSEL)/DEGRAD
                YS(K)=SIGNAL(J)
              END IF
   40       CONTINUE
            CALL PLOT(0,0.,0.,'D')
            CALL PLOT(K,XS,YS,SYMB)
            CALL PLOT(0,0.,0.,'S')
            CALL RTNCON('     LUNAR ELONGATION (deg.) -->',32)
C
C               End of plots.
C
  100       CALL SPACE
            CALL TVN('Do you want to interpolate in TIME alone, or')
            CALL TVN('do you want to try fitting sky with a MODEL?')
            CALL TV('(Enter T for Time; M to Model the sky; or)')
            CALL ASKN('(        R to re-display the plots.)',A1)
            CALL SPACE
C
  102       IF (A1.EQ.'T') THEN
                K=0
                DO 105 J=JMIN,JMAX
C
C                  restore sky data for current night and passband.
                   IF (KBND(J).NE.NB) GO TO 105
                   IF (DTYPE(J).EQ.'Y')THEN
C                    add to XS, YS.
                     K=K+1
                     XS(K)=DJOBS(J)
                     YS(K)=SIGNAL(J)
                   END IF
C
  105           CONTINUE
C
  110           CALL TV('Enter C to replace sky data with a Constant,')
                CALL TVN('      L to smooth with a Linear fit,')
                CALL TVN('      N to use the sky Nearest the star,')
                CALL TVN('      P to use the sky Preceding the star,')
                CALL TVN('      F to use the sky Following the star,')
                CALL TVN('      R to Re-display plots,')
                CALL TVN('   or Q to Quit.')
                CALL ASKN('?',A1)
                IF (A1.EQ.'C') THEN
C                   be lazy: fit line, but force it horizontal.
                    CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE)
                    SLOPE=0.
                    PG(2)=0.
                    NIX=1
                    IX(1)=2
                ELSE IF (A1.EQ.'L') THEN
                    IF (K.LT.6) THEN
                        CALL TV('Not enough data to fit a line safely.')
                        CALL TVN('Try a constant instead.')
                        GO TO 110
                    END IF
                    CALL ROBLIN(XS,YS,K, XBAR,YBAR,SLOPE)
                    PG(2)=SLOPE
                    NIX=0
                ELSE
C                  for neighboring-sample subtraction:
C
                   IF (A1.EQ.'N') THEN
C                     use NEAREST neighbor.
                      CALL TV('CAUTION: This method is not robust!')
                      CALL SKYNBR(NB,DTYPE,KDIAPHRAGM,JBGN,JEND)
                   ELSE IF (A1.EQ.'P') THEN
C                     SKY comes first.
                      CALL TV('CAUTION: This method is not robust!')
                      CALL SKYADJ(NB,DTYPE,KDIAPHRAGM,JBGN,JEND,NERROR)
                      IF (NERROR.GT.0) GO TO 2
                   ELSE IF (A1.EQ.'F') THEN
C                     STAR comes first, so scan in reverse order.
                      CALL TV('CAUTION: This method is not robust!')
                      CALL SKYADJ(NB,DTYPE,KDIAPHRAGM,JEND,JBGN,NERROR)
                      IF (NERROR.GT.0) GO TO 2
                   ELSE IF (A1.EQ.'R') THEN
                      GO TO 2
                   ELSE
                      CALL TV('Usually N is your best choice; but')
                      CALL TVN('if you have been careful to take sky')
                      CALL TVN('ALWAYS before or after every star, you')
                      CALL TVN('can choose P for sky-then-star, or')
                      CALL TVN('           F for star-then-sky.')
                      GO TO 110
                   END IF
C                  skip residual plots if sky is subtracted exactly.
                   CALL TV('(Sky samples were subtracted exactly; so')
                   CALL TVN(' there are no residuals to plot.)')
                   GO TO 999
                END IF
C
C               here when linear model was used.
C
                IPR=0
                PG(1)=YBAR
C               PG(2) was set above.
                DO 120 J=1,K
                  X(1,J)=XS(J)-XBAR
                  Y(J)=YS(J)
                  W(J)=1.
  120           CONTINUE
C               refine fit:
C     Subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
                CALL DLSQ(YPLIN,K,2,  1,    2,   NIX, 1,    0,   1.E-3)
                SLOPE=P(2)
                YBAR=P(1)
C
                YMAX=-3.E33
                YMIN=3.E33
C
                DO 130 J=JBGN,JEND
                    IF (KBND(J).NE.NB) GO TO 130
C                   subtract fitted sky:
                    SIGNAL(J)=SIGNAL(J)-((DJOBS(J)-XBAR)*SLOPE+YBAR)
                    IF (DTYPE(J).EQ.'Y') THEN
                        YMAX=MAX(YMAX,SIGNAL(J))
                        YMIN=MIN(YMIN,SIGNAL(J))
                    END IF
  130           CONTINUE
C
C               plot sky resids.
C
                CALL NEED(24)
                CARD(:26)='Sky resids.  vs. TIME for'
                CARD(27:)=BANDS(NB)
                CALL TV(CARD)
C               set plot limits....
                XLIMS(1)=DJOBS(JMIN)
                XLIMS(2)=DJOBS(JMAX)
                YLIMS(1)=YMIN
                YLIMS(2)=YMAX
                CALL PLOT(0,XLIMS,YLIMS,'L')
                DO 140 J=JBGN,JEND
                   IF (KBND(J).NE.NB) GO TO 140
                   CALL PLOT(-1,DJOBS(J),SIGNAL(J),'*')
  140           CONTINUE
C               force plot:
                CALL PLOT(1,0.,0.,'+')
                CALL RTNCON('     Time (decimal days) -->',28)
                GO TO 999
C
            ELSE IF (A1.EQ.'M') THEN
C               go on.
C
            ELSE IF (A1.EQ.'R') THEN
                GO TO 2
C
            ELSE IF (A1.EQ.'H') THEN
                CALL SPACE2
                CALL TV('Reply M to detect bad sky measurements.')
                CALL TV('Reply T to use a neighboring sky sample,')
                CALL TVN('either before, after, or nearest the star.')
                CALL SPACE2
                GO TO 100
C
            ELSE
                GO TO 100
            END IF
        ELSE IF (K.EQ.1) THEN
            CALL TV('Only 1 datum available for this night!')
            CALL SPACE2
            DO 150 J=JBGN,JEND
C               subtract sky:
                IF (KBND(J).EQ.NB)
     1              SIGNAL(J)=SIGNAL(J)-SIGMIN
  150       CONTINUE
C           skip plot.
            GO TO 999
        ELSE
            CALL TV('NO SKY DATA AVAILABLE for this night')
            Call ASK('Do you want to continue?',A1)
            IF (A1.EQ.'Y') THEN
               RETURN
            ELSE
               CALL STSEPI
            END IF
        END IF
C
C
C       Now  *****  MODEL  *****  and interpolate sky.
C
        IPR=0
        NPTS=K
C       remove const. term:
        NIX=1
        IX(1)=1
        PG(1)=0.
C       set up diaphragm factors:
        NPARMS=11
        NN=0
        NREFDIA=0
        DO 157 K=12,11+NDIAS
           PG(K)=1.
           NDIA=K-11
               WRITE(CARD,*)NDIA2NP(NDIA,1),' data in dia.',DIANAM(NDIA)
               CALL TV(CARD)
C          Set NPARMS to largest one actually used.
           IF (NDIA2NP(NDIA,1).GT.0) NPARMS=K
C          Find leading diaphragm...
           IF (NDIA2NP(NDIA,1).GT.NN) THEN
C             flag NDIA as possible NREFDIA.
              NN=NDIA2NP(NDIA,1)
              NREFDIA=NDIA
           END IF
  157   CONTINUE
C
C
C       Plot star haloes:
        K=0
        DO 160 J=1,NPTS
           IF (X(4,J).GT.0.) THEN
             K=K+1
             XS(K)=X(4,J)
             YS(K)=Y(J)
           END IF
  160   CONTINUE
        CALL ROBLIN(XS,YS,K,  XBAR, YBAR, SLOPE)
C       Fraction of star included:
        IF (SLOPE.GT.0. .AND. SLOPE.LT.0.01) THEN
          CALL SORT2(XS,YS,K)
          ZM=XS(K)
C         ZM is max. value of XS.
          CALL NEED(24)
          CALL CENTER('Sky vs. corresponding Star signal:')
          CALL TVN('    Sky')
          CALL PLOT(-K,XS,YS,' ')
          CALL PLOT(0,0.,0.,'O')
          DO 165 J=1,20
             XM=J*ZM/20.
             YM=YBAR+SLOPE*(XM-XBAR)
             CALL PLOT(-1,XM,YM,'+')
  165     CONTINUE
          CALL PLOT(0,0.,1.,'O')
          DO 166 J=1,NPTS
            IF (X(4,J).EQ.0.) GO TO 166
            CALL PLOT(-1,X(4,J),Y(J),SYMB(J))
  166     CONTINUE
          CALL PLOT(1,0.,0.,'+')
          CARD='   Star signal -->      (Sky is 0.xxxx of Star)'
          WRITE (CARD(32:38),'(F7.4)') SLOPE
          CALL RTNCON(CARD,49)
        ELSE
          CALL TV('Contribution of star to sky measures not detected.')
          SLOPE=0.
        END IF
C
C       Set default order:
        DO 167 J=1,10
           LORDER(J)=J
  167   CONTINUE
C
C       estimate PG's:
C
        PG(11)=SLOPE
C
C       Start for dark sky:
C         Fit rough line to dark sky:
          K=0
          DO 170 J=1,NPTS
            IF (KX(1,J).EQ.0) THEN
C             Add to new list:
              K=K+1
              XS(K)=X(1,J)
              YS(K)=Y(J) - PG(11)*X(4,J)
C             set initial wt.:
              W(J)=1.
            ELSE
C             zero weights if Moon is up, or twilight:
              W(J)=0.
            END IF
  170     CONTINUE
          NN=K
C
        IF (NPTS-MOONUP-NTWI.GT.6 .AND. AIRMAX.GT.1.5) THEN
C             Debugging:
C             CALL TV('DATA:')
C             CALL PLOT(-NN,XS,YS,'*')
C         NN data now set for moonless part.
C
          OLDVAR=3.E33
          NFIT=1
C         Set order of relaxation:
          LORDER(3)=4
          LORDER(4)=5
          LORDER(5)=3
C
C         Fit line to data vs. airmass:
          CALL ROBLIN(XS,YS,NN, XBAR,YBAR,SLOPE)
C             Debugging:
C             DO 174 K=1,NN
C               CALL PLOT(-1, XS(K), YBAR+SLOPE*(XS(K)-XBAR),'-')
C 174         CONTINUE
C             CALL PLOT(1,0.,0.,'+')
C             CALL RTNCON('  AIRMASS -->',13)
          PG(2)=YBAR/XBAR
C
          IF (YBAR.LE.XBAR*SLOPE) THEN
C           Negative or zero intercept; form can't be fit.
C                             CALL TV('Negative intercept.')
            PG(3)=0.D0
            PG(4)=0.D0
            PG(5)=0.D0
          ELSE
C           Positive intercept; OK.
            IF (SLOPE.GT.0.5*YBAR/XBAR) THEN
C             P(4) is small.
C                             CALL TV('P(4) small.')
              PG(4)=(SQRT(YBAR/(XBAR*SLOPE))-1.)/XBAR
              PG(5)=PG(4)**2
            ELSE
C             P(4) is big.
C                             CALL TV('P(4) big; using P(5).')
              PG(4)=1./XBAR
              PG(5)=(YBAR/XBAR-SLOPE)/((YBAR+SLOPE*XBAR)*XBAR)
              PG(2)=(YBAR/XBAR+YBAR*PG(4)+YBAR*XBAR*PG(5))
              LORDER(3)=5
              LORDER(4)=4
            END IF
          END IF
C                     Override trial line with a priori values:
                      PG(4)=EXTIN(NB)
                      PG(5)=0.5*PG(4)*PG(4)
                      PG(2)=YBAR/(1.+(PG(4)+PG(5)*XBAR)*XBAR)
        ELSE
C         Not enough data or range to fit.
                              CALL TV('Can''t fit dark sky.')
          PG(2)=YMIN*2.
          PG(4)=EXTIN(NB)
          PG(5)=0.5*PG(4)*PG(4)
        END IF
        PG(3)=PG(5)*PG(2)*0.1
        ZENITH=PG(1)+(PG(2)+PG(3))/(1.+PG(4)+PG(5))
        IF (NPTS-MOONUP-NTWI.LT.6) GO TO 200
C
C       PG's now set.
C
C       hold all fixed but P(2):
        DO 180 J=2,11
          CALL FIXP(J,PG(J),NIX)
  180   CONTINUE
        CALL UNFIXP(2,NIX)
C
C       Relax successive parameters:
        DO 190 L=3,6
           CARD='Solving for x sky parameters...'
C          IF(L.EQ.3) CARD(24:24)=' '
           WRITE(CARD(12:13),'(I2)') L-1
           CALL TV(CARD)
           IF (NPTS+NIX.LT.NPARMS+4) GO TO 191
C
C          Fix all sky parms.:
           DO 181 K=1,11
             IX(K)=K
  181      CONTINUE
           NIX=11
           IF (NDIAS.GT.0) THEN
C             Fix reference diaphragm:
              K=NREFDIA+11
              CALL FIXP(K,PG(K),NIX)
C                    Debugging:
C                    WRITE(CARD,*)'Holding PG(',K,') at',PG(K)
C                    CALL TV(CARD)
           END IF
C
C          See if we can release P(11), the star fraction:
           IF (NPTS.GT.2*MOONUP) CALL UNFIXP(11,NIX)
C
C         First, vary ONLY the new parameter.
          CALL UNFIXP(LORDER(L-1),NIX)
C
C      Subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
          CALL DLSQ(YPSKY,NPTS,NPARMS, 1,  NPARMS,NIX, 4,    2,   1.E-3)
C
C          Next, see if result is acceptable:
           IF (WVAR.GT.0.D0) PG(LORDER(L-1))=P(LORDER(L-1))
C          Check each case individually:
           IF (P(2).LT.0.D0) CALL FIXP(2,0.,NIX)
           IF (P(3).LT.0.D0) CALL FIXP(3,0.,NIX)
           IF (P(5).LT.0.D0) CALL FIXP(5,0.,NIX)
           IF (P(11).LT.0.D0) CALL FIXP(11,0.,NIX)
C
           DO 182 K=3,L
              CALL UNFIXP(LORDER(K-1),NIX)
  182      CONTINUE
C
C      Subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
          CALL DLSQ(YPSKY,NPTS,NPARMS, 1,  NPARMS,NIX, 4,    2,   1.E-3)
C
C         Give up if solution failed:
          IF (WVAR.EQ.0.D0 .OR. WVAR.GT.1.2*OLDVAR) GO TO 189
C
C         Debugging:
C             CALL TV('DATA:')
C             CALL PLOT(-NN,XS,YS,'*')
C             DO 184 K=1,NPTS
C               remember that YC is DP...
C               IF (W(K).GT.0.) CALL PLOT(-1,X(1,K),REAL(YC(K)),'C')
C 184         CONTINUE
C             CALL PLOT(1,0.,0.,'+')
C             CALL RTNCON('  AIRMASS -->',13)
C
C         Check solution for plausibility:
          IF (P(2).LE.0.D0 .OR. P(3).LT.0.D0 .OR. P(5).LT.0.D0 .OR.
     1       (P(4).LT.0.D0 .AND. P(4)*P(4).GE.4.D0*P(5)) .OR.
     2       (P(4).GT.100.D0 .OR. P(5).GE.100.D0) .OR.
     3       (P(11).GT.0.01D0 .OR. P(11).LT.0.D0) ) THEN
C            Trouble.
             GO TO 189
          ELSE
C            OK.
          END IF
C
C         Adjust weights of moonless data:
          ZENITH=PG(1)+(PG(2)+PG(3))/(1.+PG(4)+PG(5))
          DO 185 K=1,NPTS
            IF (KX(1,K).EQ.0) W(K)=EXPTIM(KX(3,K))/(ABS(YC(K)) + ZENITH)
  185     CONTINUE
C         Update parameters:
          DO 187 K=1,5
             PG(K)=P(K)
             ERRSAVE(K)=SP(K)
C             WRITE(CARD,'(''P('',I2,'') ='',G9.2)') K,P(K)
C             CALL TV(CARD)
C             WRITE(CARD,'(4X,''+/-'',G9.2)') SP(K)
C             CALL TVN(CARD)
  187     CONTINUE
C         Star wings:
          PG(11)=P(11)
C          IF (SP(11).GT.0.) THEN
C             WRITE(CARD,'(''P(11) ='',G9.2)') P(11)
C             CALL TV(CARD)
C             WRITE(CARD,'(4X,''+/-'',G9.2)') SP(11)
C             CALL TVN(CARD)
C          END IF
C         Diaphragms:
          DO 188 K=12,NPARMS
             IF (P(K).GT.0.) PG(K)=P(K)
  188     CONTINUE
          NFIT=NPARMS-NIX
          WRITE(CARD,'(I5,A)') NFIT,'-parameter fit accepted.'
          CALL TV(CARD)
          OLDVAR=WVAR
C         Release next parameter and try again.
  189     CALL UNFIXP(LORDER(L),NIX)
  190   CONTINUE
C
  191   CONTINUE
C
C       Plot fit:
        CALL SPACE
        CALL NEED(24)
        WRITE (CARD,'(A,I3,A,I4,2A)')
     1    'Fit',NFIT,' parameters to',NN,' dark-sky data in ',BANDS(NB)
        CALL CENTER(CARD)
        CALL TVN('    Sky')
        ZENITH=PG(1)+(PG(2)+PG(3))/(1.+PG(4)+PG(5))
C       Reserve space:
        XLIMS(1)=1.0
        XLIMS(2)=AIRMAX
        YLIMS(1)=0.7*ZENITH
        XM=SQRT(AIRMAX)
        YLIMS(2)=1.4*MAX(ZENITH, PG(1)+((PG(2)+AIRMAX*PG(3))*AIRMAX)/
     1                             (1.+AIRMAX*(PG(4)+AIRMAX*PG(5))),
     2         PG(1)+((PG(2)+XM*PG(3))*XM)/(1.+XM*(PG(4)+XM*PG(5))))
        CALL PLOT(0,XLIMS,YLIMS,'L')
        CALL PLOT(0,0.,0.,'O')
C       Show fitted fcn.:
        DO 194 K=1,40
           XS(1)=1.+(AIRMAX-1.)*.025*K
           YS(1)=PG(1)+((PG(2)+XS(1)*PG(3))*XS(1))/
     1                (1.+XS(1)*(PG(4)+XS(1)*PG(5)))
           CALL PLOT(-1,XS,YS,'+')
  194   CONTINUE
        CALL PLOT(0,0.,1.,'O')
C       Show data:
        DO 195 K=1,NPTS
           IF(KX(1,K).EQ.0) CALL PLOT(-1,X(1,K),Y(K)-PG(11)*X(4,K), '*')
  195   CONTINUE
        CALL PLOT(1,0.,0.,'+')
        CALL RTNCON('   Airmass -->',14)
C
C     Start for MOONLIGHT:
C
C       Rough estimates for large-airmass parameters:
  200   PG(9)=0.35*EXTIN(NB)
        PG(10)=EXP(-EXTIN(NB))
C
C         Fit rough fcn. to lunar aureole:
          K=0
          DO 210 J=1,NPTS
            IF (KX(1,J).EQ.1) THEN
C             Add to new list:
              K=K+1
              XM=X(1,J)+X(2,J)
              ZM=X(1,J)*X(2,J)
              XS(K)=X(3,J)
C                        remove dark-sky component...
              YS(K)=(Y(J) - (PG(1)+((PG(2)+X(1,J)*PG(3))*X(1,J))/
     1                        (1.+X(1,J)*(PG(4)+X(1,J)*PG(5))))
C                     ... remove star halo ...
     2                    - PG(11)*X(4,J)
C                        ...and scale to unit airmass.
     3                  ) * X(3,J)/(X(1,J)*(EXP(-PG(9)*XM)+PG(10)/ZM))
C             partial allowance for photon noise:
              W(J)=EXPTIM(KX(3,J))/(ABS(Y(J))+ZENITH)
            ELSE
C             Keep old wt. for dark sky.
            END IF
  210     CONTINUE
C
        IF (MOONUP.GT.6 .AND. ELONGMAX-ELONGMIN.GT.0.2) THEN
C         Sort on elongation & use bottom of list for aureole:
          CALL SORT2(XS,YS,K)
          DO 211 J=1,K
C            Select points in aureole:
             IF (XS(J).GT.1.) GO TO 212
  211     CONTINUE
C         Whole set is in aureole.  Use 1st half:
          J=NPTS/2
  212     IF (J.EQ.1) THEN
C            No aureole data. Pick something:
             J=3
          END IF
C                      Debugging:
C                      WRITE(CARD,*) J,' points in initial fit.'
C                      CALL TV(card)
C         First J now suitable to estimate PG(6).
          CALL ROBLIN(XS,YS,J, XBAR,YBAR,SLOPE)
          PG(6)=YBAR-XBAR*SLOPE
          PG(7)=SLOPE
          IF (PG(6).LT.0.D0) THEN
C            Aureole not detectable.
             PG(6)=0.D0
             PG(7)=YBAR/XBAR
          END IF
C         Now move down *last* third...
          DO 216 J=1,K/3
             XS(J)=XS(K-K/3+J)
C            ...and divide by elongation.
             YS(J)=YS(K-K/3+J)/XS(J)
  216     CONTINUE
          CALL ROBLIN(XS,YS,K/3, XBAR,YBAR,SLOPE)
          PG(8)=MAX(SLOPE,0.)
        ELSE
C         Not enough data to fit with Moon.
          IF (K.GT.0) THEN
C           Use median Y:
            CALL SORT2(YS,XS,K)
            PG(7)=(YS((K+1)/2)+YS(K/2+1))*0.5
          ELSE
            PG(7)=0.
          END IF
          PG(6)=0.
          PG(8)=0.
        END IF
C
C       PG's now set for moonlight terms.
C
        IF (ELONGMAX-ELONGMIN.LT.2.) THEN
C          Rearrange order of solving:
           IF (ELONGMIN.GT.1.5) THEN
              LORDER(6)=8
              LORDER(8)=10
              LORDER(10)=6
              IF (ELONGMIN.GT.2.) THEN
                 LORDER(7)=9
                 LORDER(9)=7
              END IF
           ELSE IF (ELONGMIN.GT.0.5) THEN
              LORDER(6)=7
              LORDER(7)=6
              IF (ELONGMAX.LT.1.) THEN
                 LORDER(7)=9
                 LORDER(8)=10
                 LORDER(9)=6
                 LORDER(10)=8
              END IF
           ELSE IF (ELONGMAX.LT.1.0) THEN
              LORDER(7)=9
              LORDER(8)=7
              LORDER(9)=10
              LORDER(10)=8
           END IF
        END IF
C
        OLDVAR=3.E33
        NFIT=0
        IF (MOONUP.LT.3) GO TO 241
        DO 240 J=6,10
C          Relax parameters one by one:
           DO 221 K=1,11
              IX(K)=K
  221      CONTINUE
           NIX=11
           IF (NDIAS.GT.0) THEN
C             Fix reference diaphragm:
              K=NREFDIA+11
              CALL FIXP(K,PG(K),NIX)
           END IF
C
C          Release P(11), the star fraction:
           CALL UNFIXP(11,NIX)
C
           CARD='Solving for x moonlight parameters...'
           IF(J.EQ.6) CARD(24:24)=' '
           WRITE(CARD(12:13),'(I2)') J-5
           CALL TV(CARD)
C
C          First, vary ONLY the new parameter.
           CALL UNFIXP(LORDER(J),NIX)
           IF (MOONUP.LT.NPARMS+4-NIX) GO TO 241
C
C     subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
           CALL DLSQ(YPSKY,NPTS,NPARMS, 1,NPARMS,NIX, 4,    2,   1.E-3)
C
C          Next, see if result is acceptable:
           IF (WVAR.GT.0.D0) PG(LORDER(J))=P(LORDER(J))
C          Check each case individually:
           IF (P(6).LT.0.D0) CALL FIXP(6,0.,NIX)
           IF (P(8).LT.0.D0) CALL FIXP(8,0.,NIX)
           IF (P(9).LT.0.D0) CALL FIXP(9,0.,NIX)
           IF (P(10).LT.0.D0) CALL FIXP(10,0.,NIX)
           IF (P(9).GT.9.D0 .AND. P(10).EQ.0.D0) CALL FIXP(9,9.,NIX)
           IF (P(9).GT.10.D0) CALL FIXP(9,1.,NIX)
           IF (P(11).LT.0.D0) CALL FIXP(11,0.,NIX)
C
C          Then, vary parameters 6 to J.
           DO 222 K=6,J
              CALL UNFIXP(LORDER(K),NIX)
  222      CONTINUE
           IF (MOONUP.LT.NPARMS+4-NIX) GO TO 241
C
C     subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
           CALL DLSQ(YPSKY,NPTS,NPARMS, 1,NPARMS,NIX, 4,    2,   1.E-3)
C
C          Re-extract data for plots:
           YMAX=-3.E33
           YMIN=3.E33
           K=0
           DO 226 L=1,NPTS
             IF (KX(1,L).EQ.1) THEN
C              Add to new list:
               K=K+1
               XS(K)=X(3,L)
               XM=X(1,L)+X(2,L)
               ZM=X(1,L)*X(2,L)
C                         remove dark-sky component...
               YS(K)=(Y(L) - (P(1)+((P(2)+X(1,L)*P(3))*X(1,L))/
     1                         (1.+X(1,L)*(P(4)+X(1,L)*P(5))))
C                ...subtract star halo, and scale to unit airmass.
     2           - PG(11)*X(4,L) ) / (X(1,L)*(EXP(-P(9)*XM) + P(10)/ZM))
               IF(NDIAS.GT.1) YS(K)=YS(K)/P(KX(2,L)+11)
               YMAX=MAX(YMAX,YS(K))
               YMIN=MIN(YMIN,YS(K))
             END IF
  226      CONTINUE
           NN=K
C
C          Skip plotting:
C
           IF (NIX.GT.0) GO TO 231
C
C          Fix up plot limits:
C
           XLIMS(1)=0.1
           XLIMS(2)=ELONGMAX
C          Adjust for possible outliers:
           CALL SORT2(YS,XS,K)
C          (Note reversal of Y and X.)
           YM=(YS((K+1)/2)+YS(K/2+1))*0.5
C          YM is median.
           L=(K+2)/4
           ZM=YS(K-L)-YS(L)
C          ZM is interquartile range.
           YMAX=MIN(YMAX,YM+4.*ZM)
           YMIN=MAX(YMIN,YM-4.*ZM)
           YLIMS(1)=YMIN
           YLIMS(2)=YMAX
C
C          Show result:
           CALL SPACE
           CALL NEED(24)
           WRITE(CARD,'(I5,A,I4,1X,2A)')
     1           NPARMS-NIX,'-parameter fit ( + + + ) to',MOONUP,
     2           BANDS(NB)(:LWORD(BANDS(NB))),' data ( * ) with Moon'
           CALL CENTER(CARD)
           CALL TVN('Normalized Sky')
           CALL PLOT(0,0.,0.,'O')
C          Reserve space:
           CALL PLOT(0,XLIMS,YLIMS,'L')
           IF (YLIMS(1).LT.0.) CALL XAXIS(XLIMS)
C          Plot fit to aureole terms:
           DO 230 K=1,40
              XM=0.1+(.02+.003*K)*K
              YM=(P(6)/XM + P(7) + P(8)*XM)
              CALL PLOT(-1,XM,YM,'+')
  230      CONTINUE
           CALL PLOT(0,0.,1.,'O')
           CALL PLOT(NN,XS,YS,'*')
           CALL RTNCON('     LUNAR ELONGATION (rad.) -->',32)
C
C          Give up if solution failed:
  231      IF (WVAR.EQ.0.D0 .OR. WVAR.GT.1.2*OLDVAR) GO TO 240
C          Tests for plausibility:
           IF (P(6).LT.0.D0 .OR. P(8).LT.0.D0 .OR.
     1         P(9).LT.0.D0 .OR. P(10).LT.0.D0 .OR.
     2        (P(9).GT.10.D0 .AND. P(10).EQ.0.D0) .OR.
     3       (P(11).GT.0.01D0 .OR. P(11).LT.0.D0) ) GO TO 240
C
C          Result accepted.
C
C          Revise weights:
           DO 235 K=1,NPTS
C             partial allowance for photon noise:
              IF(KX(1,K).EQ.1) W(K)=EXPTIM(KX(3,K))/(ABS(Y(K))+ZENITH)
  235      CONTINUE
C          Update parameters:
           DO 237 K=6,10
             PG(K)=P(K)
             ERRSAVE(K)=SP(K)
C             WRITE(CARD,'(''P('',I2,'') ='',G9.2)') K,P(K)
C             CALL TV(CARD)
C             WRITE(CARD,'(4X,''+/-'',G9.2)') SP(K)
C             CALL TVN(CARD)
  237      CONTINUE
C          Wings of star image:
           IF (P(11).GT.0.D0) PG(11)=P(11)
           IF (SP(11).GT.0.) THEN
C             WRITE(CARD,'(''P(11) ='',G9.2)') P(11)
C             CALL TV(CARD)
C             WRITE(CARD,'(4X,''+/-'',G9.2)') SP(11)
C             CALL TVN(CARD)
           END IF
C         Diaphragms:
          DO 238 K=12,NPARMS
             IF (P(K).GT.0.) PG(K)=P(K)
  238     CONTINUE
C
           NFIT=NPARMS-NIX
           WRITE(CARD,'(I5,A)') NFIT,'-parameter fit accepted.'
           CALL TV(CARD)
           OLDVAR=WVAR
C
  240   CONTINUE
C
C       Re-extract data for final plot:
  241   YMAX=-3.E33
        YMIN=3.E33
        K=0
        DO 244 L=1,NPTS
          IF (KX(1,L).EQ.1) THEN
C           Add to new list:
            K=K+1
            XS(K)=X(3,L)
            XM=X(1,L)+X(2,L)
            ZM=X(1,L)*X(2,L)
C                      remove dark-sky component...
            YS(K)=(Y(L) - (PG(1)+((PG(2)+X(1,L)*PG(3))*X(1,L))/
     1                       (1.+X(1,L)*(PG(4)+X(1,L)*PG(5))))
C           ...subtract star halo, and scale to unit airmass.
     2         - PG(11)*X(4,L) ) / (X(1,L)*(EXP(-PG(9)*XM) + PG(10)/ZM))
            YMAX=MAX(YMAX,YS(K))
            YMIN=MIN(YMIN,YS(K))
          END IF
  244   CONTINUE
C
C       Fix up plot limits:
        XLIMS(1)=ELONGMIN
        XLIMS(2)=ELONGMAX
        IF (K.EQ.0) GO TO 300
C       Adjust for possible outliers:
        CALL SORT2(YS,XS,K)
C       (Note reversal of Y and X.)
        YM=(YS((K+1)/2)+YS(K/2+1))*0.5
C       YM is median.
        L=(K+2)/4
        ZM=YS(K-L)-YS(L)
C       ZM is interquartile range.
        YMAX=MIN(YMAX,YM+3.*ZM)
        YMIN=MAX(YMIN,YM-3.*ZM)
        IF (YMAX.EQ.YMIN) GO TO 300
        YLIMS(1)=YMIN
        YLIMS(2)=YMAX
C
C       Show result:
        CALL SPACE
        CALL NEED(24)
        WRITE(CARD,'(I5,A,I4,1X,2A)')
     1           NFIT,'-parameter fit ( + + + ) to',MOONUP,
     2           BANDS(NB)(:LWORD(BANDS(NB))),' data ( * ) with Moon'
        CALL CENTER(CARD)
        CALL TVN('Normalized Sky')
        CALL PLOT(0,0.,0.,'O')
C       Reserve space:
        CALL PLOT(0,XLIMS,YLIMS,'L')
        IF (YLIMS(1).LT.0.) CALL XAXIS(XLIMS)
C       Plot fit to aureole terms:
        DO 250 K=1,40
           XM=ELONGMIN+(ELONGMAX-ELONGMIN)*(.02+.003*K)*K
           YM=(PG(6)/XM + PG(7) + PG(8)*XM)
           CALL PLOT(-1,XM,YM,'+')
  250   CONTINUE
        CALL PLOT(0,0.,1.,'O')
        CALL PLOT(NN,XS,YS,'*')
        CALL RTNCON('     LUNAR ELONGATION (rad.) -->',32)
C
  300   CONTINUE
C
C
C   Now display fitted equation:
C
        DO 310 J=1,10
          SP(J)=ERRSAVE(J)
  310   CONTINUE
        CALL SKYEQN(PG,SP,CARD)
C
C
C   Don't forget diaphragm scaling:
C
        IF (NPARMS.GT.11) THEN
          CALL SPACE
          CALL TV('Diaphragm scaling:')
        END IF
        DO 320 J=12,NPARMS
          WRITE(CARD,'(A,F7.3,'' +/-'',F7.3)') DIANAM(J-11),P(J),SP(J)
          CALL TV(CARD)
  320   CONTINUE
C
C
C   Now display overall fit:
C
C       Hold all fixed; evaluate YC's:
        DO 355 K=1,NPARMS
           IX(K)=K
  355   CONTINUE
        NIX=NPARMS
C
C     subroutine DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
           CALL DLSQ(YPSKY,NPTS,NPARMS, 1, NPARMS,NIX, 4,    2,   1.E-1)
C
        CALL SPACE
        CALL NEED(24)
        CARD='Measured vs. model ( + + ) sky brightness in '//BANDS(NB)
        CARD(LWORD(CARD)+1:)=':'
        CALL CENTER(CARD)
        CALL TVN('Observed sky')
        XLIMS(1)=3.E33
        YLIMS(1)=0.
        YLIMS(2)=0.
        K=0
        DO 360 J=1,NPTS
C          omit twilight.
           IF (KX(1,J).LT.0) GO TO 360
           K=K+1
           XS(K)=YC(J)
           YS(K)=Y(J)
           XLIMS(1)=MIN(XLIMS(1),XS(K))
           IF (KX(1,J).EQ.0) THEN
             SYMB(K)='*'
           ELSE IF (KX(1,J).EQ.1) THEN
             SYMB(K)='M'
C          ELSE IF (KX(1,J).EQ.-1) THEN
C            SYMB(K)='t'
           ELSE
             SYMB(K)='?'
           END IF
  360   CONTINUE
        CALL SORT2(XS,YS,K)
        IF (XS(K).LT.2.*XS(K-1)) THEN
          XLIMS(2)=XS(K)
        ELSE
          XLIMS(2)=XS(K-1)
        END IF
        YLIMS(1)=0.8*XLIMS(1)
        YLIMS(2)=1.2*XLIMS(2)
        CALL PLOT(0,XLIMS,YLIMS,'L')
C       Draw unit diagonal:
        DO 370 J=1,20
           XM=J*XLIMS(2)/20.
           CALL PLOT(-1,XM,XM,'+')
  370   CONTINUE
        CALL PLOT(0,0.,0.,'D')
        CALL PLOT(K,XS,YS,SYMB)
        CALL RTNCON('     Calculated sky -->',23)
C
C       Plot star halo resids.:
C
        CALL SPACE
        CALL NEED(24)
        CARD='Residuals vs. Star signal in '//BANDS(NB)
        CARD(LWORD(CARD)+1:)=':'
        CALL CENTER(CARD)
        CALL TVN('  Residual')
        XLIMS(1)=0.
        YLIMS(1)=0.
        YLIMS(1)=0.
        YLIMS(2)=0.
        K=0
        DO 380 J=1,NPTS
           IF (KX(1,J).LT.0) GO TO 380
           K=K+1
           XS(K)=X(4,J)
           YS(K)=Y(J)-YC(J)
           XLIMS(2)=MAX(XLIMS(2),XS(K))
           YLIMS(1)=MIN(YLIMS(1),YS(K))
           YLIMS(2)=MAX(YLIMS(2),YS(K))
  380   CONTINUE
        CALL PLOT(0,XLIMS,YLIMS,'L')
        CALL XAXIS(XLIMS)
        CALL PLOT(K,XS,YS,SYMB)
        CALL PLOT(0,0.,0.,'S')
        CARD='   Star signal -->      (Sky is 0.xxxx of Star)'
        WRITE (CARD(32:38),'(F7.4)') PG(11)
        CALL RTNCON(CARD,49)
C
C       Plot resids. vs. time:
C
        CALL SPACE
        CALL NEED(24)
        CARD='O/C Ratio vs. Time in '//BANDS(NB)
        CARD(LWORD(CARD)+1:)=':'
        CALL CENTER(CARD)
        CALL TVN('   Ratio')
        XLIMS(1)=DJOBS(JDUSK)
        XLIMS(2)=DJOBS(JDAWN)
        YLIMS(1)=0.6
        YLIMS(2)=1.5
        CALL PLOT(0,XLIMS,YLIMS,'L')
C       Draw "axis" at 1.0:
        DO 388 J=1,80
           XM=XLIMS(1)+(XLIMS(2)-XLIMS(1))*J/80.
           CALL PLOT(-1,XM,1.,'-')
  388   CONTINUE
        DO 390 J=1,NPTS
           IF (KX(1,J).LT.0) GO TO 390
           IF (YC(J).EQ.0.D0) YC(J)=1.E-6
           CALL PLOT(-1,DJOBS(KX(3,J)),Y(J)/REAL(YC(J)),SYMB(J))
  390   CONTINUE
C       force plot:
        CALL PLOT(1,10.,0.,'+')
        WRITE(CARD,'(A,F7.0)')'    TIME (MJD =',DTZERO(NIGHT)
        CARD(LWORD(CARD):)=' + decimal of day) -->'
        CALL RTNCON(CARD,LWORD(CARD))
C
C
C    List outliers:
C
        CALL SPACE2
        CALL NEED(10)
        CARD='Sky anomalies in '//BANDS(NB)
        WRITE(CARD(LWORD(CARD)+2:),'(''on MJD'',F7.0)')
     1       DTZERO(NIGHT)+INT(DJOBS((JBGN+JEND)/2))
        CALL TV(CARD)
        CARD='DECIMAL    SKY      LUNAR     LUNAR    STELLAR'
        CARD(52:)='SKY       SKY       O/C'
        CALL TV(CARD)
        CARD='  DAY    AIRMASS   AIRMASS  ELONGATION  SIGNAL  OBSERVED'
        CARD(59:)='CALCULATED  RATIO'
        CALL TVN(CARD)
        WRITE(CARD,'(''------'',7(3X,''-------''))')
        CALL TVN(CARD)
        NN=0
        DO 450 J=1,NPTS
          IF (KX(1,J).GE.0 .AND.
     1        W(J)*(Y(J)-REAL(YC(J)))**2 .GT. SCALE*36.) THEN
            WRITE(CARD,'(F6.4,F9.3,2F10.3,1PE12.2,0P2F10.3,F8.3)')
     1           DJOBS(KX(3,J)),(X(K,J),K=1,4),Y(J),YC(J),Y(J)/YC(J)
            IF (KX(1,J).NE.1) THEN
               CARD(20:35)='   -         -'
            END IF
            CALL TVN(CARD)
            NN=NN+1
          END IF
  450   CONTINUE
        IF (NN.EQ.0) THEN
           CALL TV('   (none)')
        ELSE
           WRITE(CARD,'(I3,'' anomalous sky data'')') NN
           IF (NN.EQ.1) CARD(22:23)='um'
           CALL TV(CARD)
        END IF
C
C
C    Ask for approval:
C
  500   CARD=
     1    'Do you want to subtract the fitted model from stellar data?'
        CALL ASK(CARD,A1)
        IF (A1.EQ.'Y' .OR. A1.EQ.'O') THEN
C         Go on.
        ELSE IF (A1.EQ.'N') THEN
C         Take short-cut to 110.
          SHORT=.TRUE.
          GO TO 3
        ELSE IF (A1.EQ.'R') THEN
          GO TO 2
        ELSE IF (A1.EQ.'H') THEN
          IF (NPTS.GT.20*NN .AND. SCALE.LT.0.1*ZENITH) THEN
C            Clearly a good fit.
             CALL TV('     YES, that''s a good fit.')
          ELSE IF (NPTS.LT.10*NN .OR. SCALE.GT.0.2*ZENITH) THEN
C            Clearly a bad fit.
             CALL TV('     NO, that''s a POOR fit.')
          ELSE
C            Not clear.
             CALL SPACE
             WRITE(CARD,'(I5,''% rejects; typical residual is'',F5.2)')
     1              (NN*100+NPTS/2)/NPTS,SCALE/ZENITH
             CARD(LWORD(CARD)+2:)='of zenith brightness.'
             CALL TV(CARD)
             CALL TV('Hard to say.  You may want to try it both ways,')
             CALL TVN('and compare the results.')
          END IF
          GO TO 500
        ELSE
          CALL SPACE2
          CALL TV('Enter either YES or NO, or')
          CALL TVN('H for Help, or R to Re-display plots.')
          GO TO 500
        END IF
C
C
C     Now subtract sky fit from data:
C
C    JBGN, JEND index the first and last data of ANY kind for the night.
C    JMIN, JMAX index the first and last SKY data for the night & band.
C
C
        IF (NTWI.GT.0) THEN
          CALL TV('Using nearest-neighbor subtraction in twilight...')
          CALL SKYNBR(NB,DTYPE,KDIAPHRAGM,JBGN,JDUSK)
          CALL SKYNBR(NB,DTYPE,KDIAPHRAGM,JDAWN,JEND)
        END IF
C
C       Subtract fitted fcn. from rest:
        CALL TV('Using fitted model...')
        DO 600 J=JDUSK,JDAWN
           IF (NB.NE.KBND(J)) GO TO 600
           Z(1)=AIRM(J)
           IF ( (RISE.LT.SET  .AND.
     1               DJOBS(J).GT.RISE .AND. DJOBS(J).LT.SET) .OR.
     2               (RISE .GT.SET  .AND.
     3               (DJOBS(J).GT.RISE .OR. DJOBS(J).LT.SET) )  ) THEN
C               Moon is up.
                KZ(1)=1
                T=(DTZERO(NIGHT)-51544.5D0+DJOBS(J))/36525.
                UTRAD=MOD(DJOBS(J),1.)*TWOPI
                TLST=STUTZ+(UTRAD*1.00273791)+ELONG
C               MOON uses T for lunar position, TLST for topocentric ephemeris.
                CALL MOON(T, TLST)
C               MOON returns lunar coords. in RASUN, DESUN in /CSUN/.
                HAMOON=TLST-RASUN
                COSDEL=COS(DESUN)
                SINDEL=SIN(DESUN)
                COSHA=COS(HAMOON)
                XM= -COSDEL*SIN(HAMOON)
                YM=COSPHI*SINDEL-SINPHI*COSDEL*COSHA
                ZM=SINPHI*SINDEL+COSPHI*COSDEL*COSHA
                IF (ZM.GT.0.) THEN
C                  Set Z(2) to LUNAR airmass.
                   Z(2)=1./ZM
                ELSE
C                  Moon below horizon.
                   Z(2)=100.
                END IF
                AZMOON=ATAN2(YM,XM)
                DELAZ=AZMOON-AZIM(J)
C               approximate cosz(sky) as 1/airmass.
                COSEL=ZM/AIRM(J) +
     1               SQRT((1.-ZM*ZM)*(1.-1./AIRM(J)**2))*COS(DELAZ)
C               Set Z(3) to lunar elongation.
                Z(3)=ACOS(COSEL)
           ELSE
C               Moon is down.
                KZ(1)=0
                Z(2)=0.
           END IF
C
           IF (DTYPE(J).EQ.'S') THEN
C             Star.
              Z(4)=SIGNAL(J)
           ELSE IF (DTYPE(J).EQ.'Y') THEN
C             Sky.
              CALL HALO(NB,DTYPE,KDIAPHRAGM,J,JBGN,JEND, Z(4))
           ELSE
C             Dark, etc.
              GO TO 600
           END IF
C
C          Put NDIA in KZ(2):
           DO 580 NDIA=1,NDIAS
              IF (DIAPH(J).EQ.DIANAM(NDIA)) THEN
                 KZ(2)=NDIA
                 GO TO 581
              END IF
  580      CONTINUE
           KZ(2)=0
  581      CONTINUE
C
C          Z's and KZ(1) are now set for YPSKY:
           CALL YPSKY
C          Sky brightness now returned in YT.
           SIGNAL(J)=SIGNAL(J)-YT
  600   CONTINUE
C
C
C       Prepare for next night:
C
  999   JBGN=JEND+1
C
 1000 CONTINUE
C               END of loops over bands and nights.
C
      RETURN
      END
      SUBROUTINE YPSKY
C
C   Fits  Y = [Y1 + Y2 + Y3 + P(1)] * P(diaphragm),  where
C
C
C                                 2
C            P(2)*Z(1) + P(3)*Z(1)
C   Y1 = ---------------------------   is the sky background, and
C                                 2
C        1 + P(4)*Z(1) + P(5)*Z(1)
C
C
C
C   Y2 = Z(1) * [P(6)/Z(3) + P(7) + P(8)*Z(3)] * {exp[-P(9)*v1] + P(10)/v3}
C
C                                      is the scattered moonlight.
C
C
C   Y3 = Z(4) * P(11)       is the part of the star included in sky measures.
C
C
C
C     P(1) is a possible zero-offset, if dark was not subtracted.
C
C
C     Independent variables:
C     ---------------------
C
C     Z(1) = observed airmass           v1 = Z(1) + Z(2) = sum of airmasses
C
C     Z(2) = lunar airmass              v3 = Z(1)*Z(2) = product of airmasses
C
C     Z(3) = lunar elongation
C
C     Z(4) = nearest star signal, or zero if none found.
C
C
C     KZ(1) = Moon flag (1 if Moon up)
C
C     KZ(2) = diaphragm number (KDIA)
C
C
C
C
      IMPLICIT NONE
C
C     Includes for dlsq:
      INCLUDE 'MID_REL_INCL:kingfit.inc'
      INCLUDE 'MID_REL_INCL:dlsq.inc'
C
      DOUBLE PRECISION DEN,  V1, V3, FIRST, SECOND, FZ1, EXPFAC, Y2
C
      INTEGER I, KDIA
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C
      IF (KZ(1).LT.0) THEN
C        Twilight:  Not in model.
         KIP=0
         YT=0.
         RETURN
      END IF
C
C
C                                 2
C            P(2)*Z(1) + P(3)*Z(1)
C   Y1 = ---------------------------   is the sky background...
C                                 2
C        1 + P(4)*Z(1) + P(5)*Z(1)
C
C
      KIP=6
      IP(1)=1
      IP(2)=2
      IP(3)=3
      IP(4)=4
      IP(5)=5
      IP(6)=11
C
      DEN = 1.D0 + (P(4) + P(5)*Z(1))*Z(1)
      YT = (P(2) + P(3)*Z(1))*Z(1)/DEN
C
      PART(1)=1.D0
      PART(2)=Z(1)/DEN
      PART(3)=Z(1)*PART(2)
      PART(4)= -YT*PART(2)
      PART(5)= -YT*PART(3)
      PART(11)=Z(4)
C
      YT = YT + P(1) + P(11)*Z(4)
C
      IF (KZ(1).EQ.0) THEN
C        Don't need moonlight, so finish up:
         IF (KZ(2).GT.0) THEN
C           Scale old partials by diaphragm factor...
            KDIA=KZ(2)+11
            DO 30 I=1,KIP
               PART(IP(I))=PART(IP(I))*P(KDIA)
   30       CONTINUE
C           Scale by diaphragm factor, and add to KIP:
            KIP=7
            IP(7)=KDIA
            PART(KDIA)=YT
            YT=P(KDIA)*YT
         END IF
         RETURN
      END IF
C
C
C     Here if we need to add moonlight:
C
C
C   Y2 = Z(1) * [P(6)/Z(3) + P(7) + P(8)*Z(3)] * {exp[-P(9)*v1] + P(10)/v3}
C
C                                      is the scattered moonlight.
C
      KIP=11
      IP(7)=6
      IP(8)=7
      IP(9)=8
      IP(10)=9
      IP(11)=10
C
      V1=Z(1)+Z(2)
      V3=Z(1)*Z(2)
C
      FIRST=P(6)/Z(3) + P(7) + P(8)*Z(3)
      FZ1=FIRST*Z(1)
      IF (-P(9)*V1.GT.30.D0) THEN
C        try to prevent floating-point exceptions.
         EXPFAC=EXP(30.D0)
      ELSE
         EXPFAC=EXP(-P(9)*V1)
      END IF
      SECOND=EXPFAC + P(10)/V3
C
      Y2 = FZ1*SECOND
C
      PART(7)=SECOND*Z(1)
      PART(6)=PART(7)/Z(3)
      PART(8)=Z(3)*PART(7)
      PART(9)= -V1*EXPFAC*FZ1
      PART(10)=FZ1/V3
C
      YT = YT + Y2
C
C        Diaphragm scaling:
         IF (KZ(2).GT.0) THEN
C           Scale old partials by diaphragm factor...
            KDIA=KZ(2)+11
            DO 50 I=1,KIP
               PART(IP(I))=PART(IP(I))*P(KDIA)
   50       CONTINUE
C           Scale by diaphragm factor, and add to KIP:
            KIP=12
            IP(12)=KDIA
            PART(KDIA)=YT
            YT=P(KDIA)*YT
         END IF
C
      RETURN
      END
      SUBROUTINE SKYADJ(NB,DTYPE,KDIAPHRAGM,JBGN,JEND,NERROR)
C
C       Picks an ADJACENT sky sample, and subtracts it from star.
C
C
      IMPLICIT NONE
C
      REAL SKY, PHNOISE, ERREST
      INTEGER NB, KDIAPHRAGM, JBGN, JEND, NERROR, J3, J,
     1         LWORD, JJ
C
      INCLUDE 'MID_REL_INCL:obs.inc'
C
      INCLUDE 'MID_REL_INCL:mstars.inc'
C     PARAMETER (MSTARS=1650)
C       MSTARS is max.number of star-catalog entries.
C       commons for star catalog:
      CHARACTER *32 STARS
      COMMON /SCATA/ STARS(MSTARS)
C
      CHARACTER DTYPE(MXOBS)
C
      CHARACTER A, CARD*48, PRECEDED*8, FOLLOWING*9, DSKY*4
C
C
      SKY=-3.E33
      DSKY=' '
      NERROR=0
C
C     Decide which way to scan:
C
      IF (JBGN.LT.JEND) THEN
          J3=1
          PRECEDED='preceded'
          FOLLOWING='FOLLOWING'
      ELSE
          J3=-1
          PRECEDED='followed'
          FOLLOWING='PRECEDING'
      END IF
C
      DO 134 J=JBGN,JEND,J3
C
         IF (KBND(J).NE.NB .OR. NSTR(J).LT.0) GO TO 134
C        here if bands match.
         IF (DTYPE(J).EQ.'Y') THEN
C           this is sky; reset it.
            SKY=SIGNAL(J)
            PHNOISE=PHVAR(J)
            ERREST=ESTER(J)
            IF (KDIAPHRAGM.NE.-1) THEN
C              flag sky diaphragm code; expect star to match.
               DSKY=DIAPH(J)
            END IF
         ELSE
C           this is star; subtract sky.
C
            IF (SKY.EQ.-3.E33) THEN
C              ERROR: no previous sky!
               CALL SPACE2
               CARD=STARS(NSTR(J))
               CARD(LWORD(CARD)+2:LWORD(CARD)+6)='at'
               WRITE(CARD(LWORD(CARD)+2:),'(F7.4)')DJOBS(J)
               CALL TV(CARD)
               CARD='!! Star not '//PRECEDED//' by sky !!'
               CALL TVN(CARD)
               CARD='Do you want to use the '//FOLLOWING//' sky?'
               CALL ASK(CARD,A)
C
               IF (A.EQ.'Y' .OR. A.EQ.'O')THEN
C                 look for following sky:
                  DO 132 JJ=J+J3,JEND,J3
                     IF (KBND(JJ).NE.NB .OR.
     1                   NSTR(JJ).LT.0) GO TO 132
                     IF(DTYPE(JJ).EQ.'Y') THEN
                        IF(KDIAPHRAGM.NE.-1 .AND.
     1                    DIAPH(JJ).NE.DIAPH(J)) GOTO 132
C                       fix this one, but do not set SKY.
                        SIGNAL(J)=SIGNAL(J)-SIGNAL(JJ)
C                       accumulate photon variances.
                        PHVAR(J)=PHVAR(J)+PHVAR(JJ)
                        IF (ESTER(J).NE.0.)
     1                      ESTER(J)=SQRT(ESTER(J)**2 + ESTER(JJ)**2)
                        GO TO 134
                     END IF
  132             CONTINUE
                  CALL TV('No sky found! Quitting.')
                  CALL STSEPI
               ELSE
                  GO TO 999
               END IF
C
            ELSE
C              probably normal case.  Check DIAPH:
               IF (KDIAPHRAGM.NE.-1 .AND.
     1                DIAPH(J).NE.DSKY) THEN
                  CALL SPACE2
                  CARD=STARS(NSTR(J))
                  CARD(LWORD(CARD)+2:LWORD(CARD)+6)='at'
                  WRITE(CARD(LWORD(CARD)+2:),
     1                '(F7.4)')DJOBS(J)
                  CALL TV(CARD)
                  CALL TV(
     1              'Sky and Star taken through DIFFERENT diaphragms.')
                  CALL TVN('Sky subtraction impossible.')
                  GO TO 999
               END IF
C              NORMAL case. Subtract sky:
               SIGNAL(J)=SIGNAL(J)-SKY
               PHVAR(J)=PHVAR(J)+PHNOISE
               IF(ESTER(J).NE.0.) ESTER(J)=SQRT(ESTER(J)**2 + ERREST**2)
            END IF
C
         END IF
C
  134 CONTINUE
C
C       Normal return.
C
      RETURN
C
C
C       error return:
C
  999 CALL ASK('Do you want to try another sky-subtraction method?',A)
      IF (A.EQ.'Y' .OR. A.EQ.'O') THEN
          NERROR=1
          RETURN
      ELSE IF (A.EQ.'H') THEN
                   CALL TV('The photometric system changes with')
                   CALL TVN('diaphragm size, so you cannot combine')
                   CALL TVN('data taken through different holes.')
                   CALL RTNCON('Try the "Nearest" option.',25)
                   NERROR=1
                   RETURN
      ELSE
          CALL TV('Quitting.')
          CALL STSEPI
      END IF
C
      END
      SUBROUTINE HALO(NB,DTYPE,KDIAPHRAGM,J,JBGN,JEND, SKYUSE)
C
C       Picks the NEAREST STAR sample to sky(J), and returns it in SKYUSE.
C
C  Arguments:
C
C	NB	is band number
C	DTYPE   is 'Y' for sky, 'S' for star
C	KDIAPHRAGM is column pointer for diaphragm data
C	JBGN	is observation number for start of night
C	JEND	is observation number for end of night
C
C
      IMPLICIT NONE
C
C
      INTEGER NB, KDIAPHRAGM, JBGN, JEND
C
      INCLUDE 'MID_REL_INCL:obs.inc'
C
      CHARACTER DTYPE(MXOBS)
C
C
C
      REAL SKY, SKY2, SKYAIR, SKYAZ, SKYTIM,   PI,TWOPI,
     1   SKYAIR2, SKYAZ2, SKYTIM2, S1, S2, SKYUSE
C
      INTEGER J, JJ
C
C   Meaning of local variables:
C
C		Previous	Following
C		--------	---------
C     Star        SKY		  SKY2
C     Airmass	SKYAIR		SKYAIR2
C     Azimuth	SKYAZ		SKYAZ2
C     Time	SKYTIM		SKYTIM2
C     Diaphragm  DSKY		 DSKY2
C
C
      INCLUDE 'MID_REL_INCL:mstars.inc'
C     PARAMETER (MSTARS=1650)
C       MSTARS is max.number of star-catalog entries.
C       commons for star catalog:
      CHARACTER *32 STARS
      COMMON /SCATA/ STARS(MSTARS)
C
      CHARACTER DSKY*4, DSKY2*4
C
      DATA PI/3.1415926535/
C
C
C     *****  BEGIN  *****
C
      TWOPI=PI+PI
      SKY=-3.E33
      SKY2=-3.E33
      DSKY=' '
      DSKY2=' '
      SKYAIR2 = 0.0
      SKYAZ2  = 0.0
      SKYTIM2 = 0.0
C
C     J is current sky datum in band NB.
C
      IF (DTYPE(J).EQ.'Y') THEN
C        this is sky.
         IF (KDIAPHRAGM.NE.-1) THEN
C            save new diaphragm.
             DSKY=DIAPH(J)
         ELSE
C            no diaphragm info.
         END IF
      ELSE
         CALL TV('Called HALO with non-sky datum.')
         RETURN
      END IF
C
C     This is sky; look for PREVIOUS star.
   20       IF (SKY.LT.0.) THEN
C                No previous star.  Find one.
                 DO 28 JJ=J-1,JBGN,-1
                     IF (KBND(JJ).NE.NB .OR. NSTR(JJ).LT.0) GO TO 28
                     IF(DTYPE(JJ).EQ.'S') THEN
                        IF(KDIAPHRAGM.NE.-1 .AND.
     1                    DIAPH(JJ).NE.DIAPH(J)) GOTO 28
C                       OK; JJ is the previous star we seek.
                        DSKY=DIAPH(JJ)
                        SKY=SIGNAL(JJ)
                        SKYAIR=AIRM(JJ)
                        SKYAZ=AZIM(JJ)
                        SKYTIM=DJOBS(JJ)
                        GO TO 40
                     END IF
   28            CONTINUE
C                No previous star.
                 SKY=-3.E33
                 DSKY=' '
            ELSE
C               probably normal case.  Check DIAPH:
                IF (KDIAPHRAGM.NE.-1 .AND. DIAPH(J).NE.DSKY) THEN
C                  Mismatch; force search for previous star.
                   SKY=-3.E33
                   GO TO 20
                ELSE
C                  NORMAL case. Existing star OK.
                END IF
            END IF
C
C        SKY is now valid previous value. Look for FOLLOWING value:
C
   40   IF (SKY2.LT.0.) THEN
C           Look for next star.
            DO 48 JJ=J+1,JEND
                  IF (KBND(JJ).NE.NB .OR. NSTR(JJ).LT.0) GO TO 48
                  IF(DTYPE(JJ).EQ.'S') THEN
                        IF(KDIAPHRAGM.NE.-1 .AND.
     1                    DIAPH(JJ).NE.DIAPH(J)) GOTO 48
C                       OK; JJ is the following star we seek.
                        DSKY2=DIAPH(JJ)
                        SKY2=SIGNAL(JJ)
                        SKYAIR2=AIRM(JJ)
                        SKYAZ2=AZIM(JJ)
                        SKYTIM2=DJOBS(JJ)
                        GO TO 60
                  END IF
   48       CONTINUE
C           No following star.
            SKY2=-3.E33
            DSKY2=' '
        ELSE
C           existing sky2 might still be OK. Check it:
            IF (SKYTIM2.GT.DJOBS(J)) THEN
C               Time is OK.  Check diaphragm:
                IF (KDIAPHRAGM.NE.-1) THEN
                    IF (DIAPH(J).EQ.DSKY2) THEN
C                       we are OK; go on.
                    ELSE
C                       need a new SKY2.
                        SKY2=-3.E33
                        GO TO 40
                    END IF
                ELSE
C                   No problem.  Go on.
                END IF
            ELSE
C               We are past the old second star.  Move it back & get new one:
                SKY=SKY2
                DSKY=DSKY2
                SKYAIR=SKYAIR2
                SKYAZ=SKYAZ2
                SKYTIM=SKYTIM2
                SKY2=-3.E33
                GO TO 40
            END IF
        END IF
C
C       Now compare to find "nearer" star:
C
   60   IF (SKY.GT.0. .AND. SKY2.GT.0.) THEN
C           Everything is hunky-dory.  Find closer sample:
            S1=ABS(SKYTIM-DJOBS(J))*20. + ABS(SKYAIR-AIRM(J))
            IF (ABS(SKYAZ-AZIM(J)).GT.PI) THEN
               S1=S1+ (SKYAIR+AIRM(J))*(TWOPI-ABS(SKYAZ-AZIM(J)))
            ELSE
               S1=S1+ (SKYAIR+AIRM(J))*ABS(SKYAZ-AZIM(J))
            END IF
            S2=ABS(SKYTIM2-DJOBS(J))*20. + ABS(SKYAIR2-AIRM(J))
            IF (ABS(SKYAZ2-AZIM(J)).GT.PI) THEN
               S2=S2+ (SKYAIR2+AIRM(J))*(TWOPI-ABS(SKYAZ2-AZIM(J)))
            ELSE
               S2=S2+ (SKYAIR2+AIRM(J))*ABS(SKYAZ2-AZIM(J))
            END IF
C
            IF (S2.LT.S1) THEN
                SKYUSE=SKY2
                IF (S2.GT.0.1) THEN
C                  complain.
C                  CALL SKYFAR()
                END IF
            ELSE
                SKYUSE=SKY
                IF (S1.GT.0.1) THEN
C                  complain.
                END IF
            END IF
        ELSE IF (SKY.GT.0.) THEN
C           SKY is the only good one.  Use it.
            SKYUSE=SKY
        ELSE IF (SKY2.GT.0.) THEN
C           SKY2 is the only good one.  Use it.
            SKYUSE=SKY2
        ELSE
C           NEITHER star is any good.
C           set value to zero and return:
            SKYUSE=0.
            RETURN
        END IF
C
C       OK, we have now rigged SKYUSE to be the value to use.  Report it.
C
C
C       Normal return.
C
      RETURN
C
      END
      SUBROUTINE SKYNBR(NB,DTYPE,KDIAPHRAGM,JBGN,JEND)
C
C       Picks the NEAREST sky sample, and subtracts it from star.
C
C  Arguments:
C
C	NB	is band number
C	DTYPE   is 'Y' for sky
C	KDIAPHRAGM is column pointer for diaphragm data
C	JBGN	is observation number for start of night
C	JEND	is observation number for end of night
C
C
      IMPLICIT NONE
C
C
      INTEGER NB, KDIAPHRAGM, JBGN, JEND
C
      INCLUDE 'MID_REL_INCL:obs.inc'
C
      CHARACTER DTYPE(MXOBS)
C
C
C
      REAL SKY, SKY2, SKYAIR, SKYAZ, SKYTIM,   PI,TWOPI,
     1   PHNOISE, ERREST, SKYAIR2, SKYAZ2, SKYTIM2, 
     2         PHN2, ERR2, S1, S2, SKYUSE, PHNUSE, ERRUSE
C
      INTEGER J, JJ
C
C   Meaning of local variables:
C
C		Previous	Following
C		--------	---------
C     Sky         SKY		  SKY2
C     Airmass	SKYAIR		SKYAIR2
C     Azimuth	SKYAZ		SKYAZ2
C     Time	SKYTIM		SKYTIM2
C     Phot.noisePHNOISE		 PHN2
C     Est.error ERREST		 ERR2
C     Diaphragm  DSKY		 DSKY2
C
C
      INCLUDE 'MID_REL_INCL:mstars.inc'
C     PARAMETER (MSTARS=1650)
C       MSTARS is max.number of star-catalog entries.
C       commons for star catalog:
      CHARACTER *32 STARS
      COMMON /SCATA/ STARS(MSTARS)
C
      CHARACTER CARD*48, DSKY*4, DSKY2*4
C
      DATA PI/3.1415926535/
C
C
      TWOPI=PI+PI
      SKY=-3.E33
      SKY2=-3.E33
      DSKY=' '
      DSKY2=' '
C
      DO 100 J=JBGN,JEND
C
         IF (KBND(J).NE.NB .OR. NSTR(J).LT.0) GO TO 100
C        here if bands match.
C
         IF (DTYPE(J).EQ.'Y') THEN
C           this is sky.
            IF (KDIAPHRAGM.NE.-1) THEN
C               save new diaphragm.
                DSKY=DIAPH(J)
            ELSE
C               no diaphragm info.
            END IF
            SKY=SIGNAL(J)
            SKYAIR=AIRM(J)
            SKYAZ=AZIM(J)
            SKYTIM=DJOBS(J)
            PHNOISE=PHVAR(J)
            ERREST=ESTER(J)
            GO TO 100
         ELSE
C           this is star; look for PREVIOUS sky.
   20       IF (SKY.LT.0.) THEN
C                No previous sky.  Find one.
                 DO 28 JJ=J-1,JBGN,-1
                     IF (KBND(JJ).NE.NB .OR. NSTR(JJ).LT.0) GO TO 28
                     IF(DTYPE(JJ).EQ.'Y') THEN
                        IF(KDIAPHRAGM.NE.-1 .AND.
     1                    DIAPH(JJ).NE.DIAPH(J)) GOTO 28
C                       OK; JJ is the previous sky we seek.
                        DSKY=DIAPH(JJ)
                        SKY=SIGNAL(JJ)
                        SKYAIR=AIRM(JJ)
                        SKYAZ=AZIM(JJ)
                        SKYTIM=DJOBS(JJ)
                        PHNOISE=PHVAR(J)
                        ERREST=ESTER(J)
                        GO TO 40
                     END IF
   28            CONTINUE
C                No previous sky.
                 SKY=-3.E33
                 DSKY=' '
            ELSE
C               probably normal case.  Check DIAPH:
                IF (KDIAPHRAGM.NE.-1 .AND. DIAPH(J).NE.DSKY) THEN
C                  Mismatch; force search for previous sky.
                   SKY=-3.E33
                   GO TO 20
                ELSE
C                  NORMAL case. Existing sky OK.
                END IF
            END IF
         END IF
C
C        SKY is now valid previous value. Look for FOLLOWING value:
C
   40   IF (SKY2.LT.0.) THEN
C           Look for next sky.
            DO 48 JJ=J+1,JEND
                  IF (KBND(JJ).NE.NB .OR. NSTR(JJ).LT.0) GO TO 48
                  IF(DTYPE(JJ).EQ.'Y') THEN
                        IF(KDIAPHRAGM.NE.-1 .AND.
     1                    DIAPH(JJ).NE.DIAPH(J)) GOTO 48
C                       OK; JJ is the following sky we seek.
                        DSKY2=DIAPH(JJ)
                        SKY2=SIGNAL(JJ)
                        SKYAIR2=AIRM(JJ)
                        SKYAZ2=AZIM(JJ)
                        SKYTIM2=DJOBS(JJ)
                        PHN2=PHVAR(JJ)
                        ERR2=ESTER(JJ)
                        GO TO 60
                  END IF
   48       CONTINUE
C           No following sky.
            SKY2=-3.E33
            DSKY2=' '
        ELSE
C           existing sky2 might still be OK. Check it:
            IF (SKYTIM2.GT.DJOBS(J)) THEN
C               Time is OK.  Check diaphragm:
                IF (KDIAPHRAGM.NE.-1) THEN
                    IF (DIAPH(J).EQ.DSKY2) THEN
C                       we are OK; go on.
                    ELSE
C                       need a new SKY2.
                        SKY2=-3.E33
                        GO TO 40
                    END IF
                ELSE
C                   No problem.  Go on.
                END IF
            ELSE
C               We are past the old second sky.  Move it back & get new one:
                SKY=SKY2
                DSKY=DSKY2
                SKYAIR=SKYAIR2
                SKYAZ=SKYAZ2
                SKYTIM=SKYTIM2
                PHNOISE=PHN2
                ERREST=ERR2
                SKY2=-3.E33
                GO TO 40
            END IF
        END IF
C
C       Now compare to find "nearer" sky:
C
   60   IF (SKY.GT.0. .AND. SKY2.GT.0.) THEN
C           Everything is hunky-dory.  Find closer sample:
            S1=ABS(SKYTIM-DJOBS(J))*20. + ABS(SKYAIR-AIRM(J))
            IF (ABS(SKYAZ-AZIM(J)).GT.PI) THEN
               S1=S1+ (SKYAIR+AIRM(J))*(TWOPI-ABS(SKYAZ-AZIM(J)))
            ELSE
               S1=S1+ (SKYAIR+AIRM(J))*ABS(SKYAZ-AZIM(J))
            END IF
            S2=ABS(SKYTIM2-DJOBS(J))*20. + ABS(SKYAIR2-AIRM(J))
            IF (ABS(SKYAZ2-AZIM(J)).GT.PI) THEN
               S2=S2+ (SKYAIR2+AIRM(J))*(TWOPI-ABS(SKYAZ2-AZIM(J)))
            ELSE
               S2=S2+ (SKYAIR2+AIRM(J))*ABS(SKYAZ2-AZIM(J))
            END IF
C
            IF (S2.LT.S1) THEN
                SKYUSE=SKY2
                PHNUSE=PHN2
                ERRUSE=ERR2
                IF (S2.GT.0.1) THEN
C                  complain.
C                  CALL SKYFAR()
                END IF
            ELSE
                SKYUSE=SKY
                PHNUSE=PHNOISE
                ERRUSE=ERREST
                IF (S1.GT.0.1) THEN
C                  complain.
                END IF
            END IF
        ELSE IF (SKY.GT.0.) THEN
C           SKY is the only good one.  Use it.
            SKYUSE=SKY
            PHNUSE=PHNOISE
            ERRUSE=ERREST
        ELSE IF (SKY2.GT.0.) THEN
C           SKY2 is the only good one.  Use it.
            SKYUSE=SKY2
            PHNUSE=PHN2
            ERRUSE=ERR2
        ELSE
C           AARRRGHHH!!  NEITHER sky is any good!
            CARD='Can''t find a sky sample to match '//STARS(NSTR(J))
            CALL TV(CARD)
            WRITE(CARD,'(''at '',F7.4)') DJOBS(J)
            IF (KDIAPHRAGM.NE.-1) CARD(30:)='in diaphragm '//DIAPH(J)
            CALL TVN(CARD)
            CALL TV('This value will be REJECTED.')
            CALL SPACE
C           set value to zero and continue:
            SKYUSE=SIGNAL(J)
            PHVAR(J)=1.E10
            ERRUSE=1.E10
        END IF
C
C       OK, we have now rigged SKYUSE to be the value to subtract.  Do it.
C
        SIGNAL(J)=SIGNAL(J)-SKYUSE
        PHVAR(J)=PHVAR(J)+PHNUSE
        IF (ESTER(J).NE.0.) ESTER(J)=SQRT(ESTER(J)**2 +ERRUSE**2)
C
  100 CONTINUE
C
C       Normal return.
C
      RETURN
C
      END
      SUBROUTINE REDSUB(BANDS,KMIN,KMAX,LASTK,MAXK,DTYPE)
C
C       Subtracts red-leak data from corresponding bands.
C       Picks the NEAREST red-leak sample, and subtracts it from star.
C
C
      IMPLICIT NONE
C
      REAL RLFACT, RL, RLVAR, ERVAR
      INTEGER KMIN, KMAX, LASTK, MAXK, NBRL, K, KRL, JBGN, NIGHT, 
     1         KK, J, JMIN, JMAX, JEND, ISTR, JOFF, J1, J2
C
      INCLUDE 'MID_REL_INCL:obs.inc'
C
      INCLUDE 'MID_REL_INCL:mstars.inc'
C     PARAMETER (MSTARS=1650)
C       MSTARS is max.number of star-catalog entries.
C       commons for star catalog:
      CHARACTER *32 STARS
      COMMON /SCATA/ STARS(MSTARS)
C
      CHARACTER DTYPE(MXOBS)
C
      INCLUDE 'MID_REL_INCL:mbands.inc'
C     PARAMETER (MBANDS=9)
      INTEGER MXCOLR
      PARAMETER (MXCOLR=3*MBANDS)
      CHARACTER *8 BANDS(MXCOLR)
C
C       Common for red leaks:
C
      COMMON /REDLK/ NBRL(MXCOLR), RLFACT(MXCOLR)
C
C
C       Useful   bands run from KMIN to LASTK.
C       Measured bands run from KMIN to KMAX, including red leaks.
C       MAXK  marks end of expected red-leak bands, starting at NBANDS+1.
C
C
C    EXAMPLE of red-leak storage:
C
C       In UBVRI system, NBANDS = 5.  Suppose we measure only B,V,R, & BRL.
C
C               K     BAND(K)    NBRL(K)
C               -     -------    -------
C               1        U          0
C               2        B          6    <-- KMIN = 2
C               3        V          0
C               4        R          0    <-- LASTK = 4
C               5        I          0
C               6       BRL         0    <-- KMAX = MAXK = 6
C
C
C
C
      DO 1000 K=KMIN,LASTK
          IF (NBRL(K).EQ.0) GO TO 1000
         KRL=NBRL(K)
         JBGN=1
C
      DO 900 NIGHT=1,NIGHTS
         KK=0
         DO 10 J=JBGN,NDATA
C           get limits for current night & passband.
            IF (NITE(J).NE.NIGHT) GO TO 20
            IF (DTYPE(J).EQ.'Y') GO TO 10
            KK=KK+1
            IF (KK.EQ.1) JMIN=J
            JMAX=J
   10    CONTINUE
         J=NDATA+1
   20    JEND=J-1
C
C       JBGN, JEND index the first & last data of ANY kind for the night.
C       JMIN, JMAX index the first & last STAR data for the night & band.
C       KK counts star data for the night.
C
C
      DO 100 J=JMIN,JMAX
C
         IF (DTYPE(J).EQ.'Y' .OR. NSTR(J).LT.0) GO TO 100
C        here if datum is stellar.
C
         IF (KBND(J).EQ.K) THEN
C
C           this observation requires red-leak correction.
            ISTR=NSTR(J)
C           scan forward & back for nearest red-leak datum for same star:
            DO 50 JOFF=1,50
C               look ahead...
                J1=J+JOFF
                IF (J1.GT.JMAX) GO TO 40
                IF (KBND(J1).NE.KRL .OR. NSTR(J1).NE.ISTR) GO TO 40
                RL=SIGNAL(J1)
                RLVAR=PHVAR(J1)
                ERVAR=ESTER(J1)
                GO TO 60
   40           CONTINUE
C               look back...
                J2=J-JOFF
                IF (J2.LT.JMIN) GO TO 50
                IF (KBND(J2).NE.KRL .OR. NSTR(J2).NE.ISTR) GO TO 50
                RL=SIGNAL(J2)
                RLVAR=PHVAR(J2)
                ERVAR=ESTER(J2)
                GO TO 60
   50       CONTINUE
            CALL TV('No red-leak datum found for')
            CALL TVN(STARS(ISTR))
C
   60       CONTINUE
            IF (RL.LT.0.) THEN
               RL=0.
               CALL TV('Negative red-leak found for')
               CALL TVN(STARS(NSTR(J)))
               GO TO 100
            END IF
C
C           Correct for red leak:
            SIGNAL(J)=SIGNAL(J)-RL*RLFACT(K)
            PHVAR(J)=PHVAR(J)+RLVAR*RLFACT(K)
            IF(ESTER(J).NE.0.)ESTER(J)=SQRT(ESTER(J)**2+ERVAR*RLFACT(K))
C
        END IF
C
  100 CONTINUE
C
C     prepare for next pass.
      JBGN=JEND+1
  900 CONTINUE
C
 1000 CONTINUE
C
C       Normal return.
C
      RETURN
C
      END
      SUBROUTINE NEED(LINES)
C
C       Tell it how many lines you need for a plot.  It reads the LOG
C       keyword and adds spaces to start a new page if necessary.
C
C
      IMPLICIT NONE
C
      INTEGER LINES, NGOT, LINE, IUNIT, NULLS, ISTAT, NLPP, I
C
C
      CALL STKRDI('LOG',6,1, NGOT,LINE,IUNIT,NULLS,ISTAT)
      IF (ISTAT.NE.0) CALL TERROR(9999,'Could not read Keyword LOG')
      CALL STKRDI('LOG',7,1, NGOT,NLPP,IUNIT,NULLS,ISTAT)
      IF (ISTAT.NE.0) CALL TERROR(9999,'Could not read Keyword LOG')
      IF(LINE+LINES.GT.NLPP)THEN
C        we must space to a new page.
         DO 1 I=1,NLPP-LINE
            CALL TVN(' ')
    1    CONTINUE
      ELSE
C        we are OK.
      END IF
C
      RETURN
      END
