      SUBROUTINE PHOTOM
C
C   Does simple photometry on each star and filter.
C   28 April 1989
C
      IMPLICIT UNDEFINED(A-Z)
      SAVE
      INTEGER *4 MAXPAR
      PARAMETER (MAXPAR=20)
      REAL *4 PI, RAD
      PARAMETER ( PI = 3.1415926535 )
      PARAMETER ( RAD = 180. / PI )
      INCLUDE 'VPLOT.INC'
      INCLUDE 'HARDCOPY.INC'
      INTEGER *4 SYMBOL(5), SPTR, N, I, J, K, L, NPTS, NUNK
      INTEGER *4 ISTAR, ICOLOR, ICHAN, MASK(4), WHERE(CNT)
      INTEGER *4 N, IZERO, XTYPE, ICUR, COUNT
      REAL *4 XPLOT(CNT), YPLOT(CNT), ZPLOT(CNT), NPHOT
      REAL *4 YFIT(CNT), DELT, ERROR(CNT), DELTAM(CNT)
      REAL *4 RX, RY, DELTA, AMCO(4), WT
      REAL *4 XMIN, XMAX, YMIN, YMAX, XHI, XLO, YHI, YLO
      REAL *4 DX, DY, DIST, DIST2, AIRM(CNT), DMAX, AVG, RMS, SCALE
      CHARACTER * 1 CH
      CHARACTER *40 TITLE(6)
      LOGICAL HCOPY
C**********************************************************************
      REAL *8 A(MAXPAR,MAXPAR), AA(MAXPAR,MAXPAR), WEIGHT
      REAL *8 CC(MAXPAR), DD(MAXPAR), EE(MAXPAR), BB, VV
C----------------------------------------------------------------------
      DATA SYMBOL  / 16, 3, 5, 6, 13 /
      DATA MASK   / 1, 2, 4, 8 /
C.......................................................................
      DATA XTYPE, ICHAN    / 1, 1 /
      HCOPY  = .FALSE.
C-----------------------------------------------------------------------
  100 CONTINUE
      IF ( HCOPY ) THEN
          NPLOTS= NPLOTS + 1
          WRITE(PLOTFILE(6:8),'(I3.3)') NPLOTS
          CALL PGBEGIN(0, PLOTFILE//PLOT_DEV, 1, 1 )
          WRITE(6,*) ' CREATING FILE ', PLOTFILE
      ELSE
          CALL PGBEGIN(0, '/EGA', 1, 1)
      END IF
  110 CONTINUE
      YMAX = -1.E6
      YMIN =  1.E6
      XMAX = -1.E6
      XMIN =  1.E6
C-----------------------------------------------------------------------
C        Read the data into the arrays
C
      NUNK = 2
      NPTS = 0
      CALL LS(1,NUNK,BB,CC,DD,EE,VV, WEIGHT, A, AA, MAXPAR)
      DO 150 ISTAR = 1, NSTAR
          DO 140 I = 1, ISCAN
              IF ( STAR(I) .NE. SLIST(ISTAR)            ) GO TO 140
              IF ( IAND( OKSCAN(I),MASK(ICHAN) ) .NE. 0 ) GO TO 140
              IF ( BDEN(I,ICHAN) .LT. D(I,ICHAN)        ) GO TO 140
              IF ( BRMS(I,ICHAN) .LT. 0.                ) GO TO 140
              IF ( ZD(I)         .EQ. 0.                ) GO TO 140
              NPTS = NPTS + 1
              AIRM(NPTS) = 1. / COS ( PI * ZD(I) / 180. )
              NPHOT = BDEN(I,ICHAN)-D(I,ICHAN)
              IF ( IAND( NDFSTATUS(I), MASK(ICHAN) ) .NE. 0 ) THEN
                  NPHOT = NDFRATIO(ICHAN) * NPHOT
              END IF
              ERROR(NPTS) = 2.5*LOG10( 1. + BRMS(I,ICHAN)
     $                               / (BDEN(I,ICHAN)-D(I,ICHAN)) )
              IF ( ERROR(NPTS).LT.0 ) THEN
                  WRITE(6,*)
     $                   ' YOU SCREWED UP AGAIN, ERROR IS < ZERO!'
              END IF
              IF ( XTYPE .EQ. 1 ) THEN
                  XPLOT(NPTS) = HA(I)
              ELSE
                  XPLOT(NPTS) =   AIRM(NPTS)
              END IF
              YPLOT(NPTS) = MAG(ISTAR) + 2.5*LOG10( NPHOT )
     $                           -  BV(ISTAR) *  CTERM(ICHAN)
              ZPLOT(NPTS) = ISTAR
              WHERE(NPTS) = I
C              WRITE(6,'(2I5,2X,3F10.3)') I, SLIST(ISTAR),
C     $             XPLOT(NPTS), YPLOT(NPTS), ERROR(NPTS)
              IF ( NPTS .EQ. 1 ) THEN
                  XMIN = XPLOT(NPTS)
                  XMAX = XPLOT(NPTS)
                  YMIN = YPLOT(NPTS)
                  YMAX = YPLOT(NPTS)
              ELSE
                  XMIN = MIN ( XPLOT(NPTS), XMIN )
                  XMAX = MAX ( XPLOT(NPTS), XMAX )
              END IF
              CC(1) = 1.0
              CC(2) = AIRM (NPTS)
              BB    = YPLOT(NPTS)
C          
C      LS mode = 2, increment coefficient matrix for least sq fit.
C         mode = 3, increment solution vector    for least sq fit.
C                   The weight is included as real*8
C
              IF ( ERROR(NPTS) .LT. .02 ) THEN
                  WEIGHT = 1./(.02*.02)
              ELSE
                  WEIGHT = 1./(ERROR(NPTS)*ERROR(NPTS) )
              END IF
              CALL LS(2,NUNK,BB,CC,DD,EE,VV, WEIGHT, A, AA, MAXPAR)
              CALL LS(3,NUNK,BB,CC,DD,EE,VV, WEIGHT, A, AA, MAXPAR)
  140     CONTINUE
  150 CONTINUE
      IF ( NPTS .GT. NUNK ) THEN
          CALL LS(4,NUNK,BB,CC,DD,EE,VV, WEIGHT, A, AA, MAXPAR )
      ELSE
          WRITE(6,*) ' NOT ENOUGH DATA TO FIT THE PHOTOMETRY '
          DO 170 J = 1, NUNK 
              DD(J) = 0. 
  170     CONTINUE 
      END IF
C-----------------------------------------------------------------------
C   Save the fits
C
      MAG0(ICHAN) = DD(1)
      AMCO(ICHAN) = DD(2)
C-----------------------------------------------------------------------
C     Determine the mean residual magnitude for each star subtract
C     That amount from each scan.
C     ZPLOT is the star sequence number for each plotted scan.
      YMIN =  1.E6
      YMAX = -1.E6
      DO 200 ISTAR = 1, NSTAR
          COUNT = 0.
          DELTAM(ISTAR) = 0
          DO 185 I = 1, NPTS
              IF ( ZPLOT(I) .EQ. ISTAR ) THEN
                  COUNT = COUNT + 1
                  DELTAM(ISTAR) = DELTAM(ISTAR) + YPLOT(I)
     $                          - MAG0(ICHAN) - AMCO(ICHAN)*AIRM(I)
              END IF
  185     CONTINUE
          IF ( COUNT .EQ. 0 ) GO TO 200
          DELTAM(ISTAR) = DELTAM(ISTAR) / FLOAT(COUNT)
          DO 190 I = 1, NPTS
              IF ( ZPLOT(I) .EQ. ISTAR ) THEN
                  YPLOT(I) = YPLOT(I) - DELTAM(ISTAR)
                  YMIN = MIN ( YPLOT(I), YMIN )
                  YMAX = MAX ( YPLOT(I), YMAX )
              END IF
  190     CONTINUE
  200 CONTINUE
C-----------------------------------------------------------------------
C     Determine the residual for each scan and set the cursor on the
C     scan with the largest deviation.
C
      AVG  = 0.
      RMS  = 0.
      COUNT = 0.
      DMAX = 0
      DO 210 I = 1, NPTS
          YFIT(I) = MAG0(ICHAN) + AMCO(ICHAN)*AIRM(I)
          DELTA = ABS( YPLOT(I) - YFIT(I) )
          IF ( DELTA .GT. DMAX ) THEN
              ICUR = I
              DMAX = DELTA
          END IF
          AVG = AVG + YPLOT(I)
          RMS = RMS + YPLOT(I)*YPLOT(I)
          COUNT = COUNT + 1
  210 CONTINUE
      IF ( COUNT .GT. 1 ) THEN
          AVG = AVG / COUNT
          RMS = SQRT( RMS/COUNT - AVG*AVG )
      END IF
      CALL PGSCI ( 5 )
      CALL PGSCH ( 1.5 )
      CALL PGRNGE ( XMIN, XMAX, XLO, XHI )
      CALL PGRNGE ( YMIN, YMAX, YLO, YHI )
C                  Make y-axis plot backwards
      CALL PGENV  ( XLO,  XHI,  YHI, YLO, 0, 1 )
      IF ( XTYPE .EQ. 1 ) THEN
          TITLE(1) = 'HOUR ANGLE (HOURS)'
      ELSE
          TITLE(1) = 'AIRMASS'
      END IF
      TITLE(2) = 'SYSTEM MAGNITUDE FOR (B-V)=0'
      TITLE(3) = DATE
      WRITE(TITLE(4), 1305) ICHAN, FILTERID(ICHAN)
      TITLE(5) = 'COLORS DENOTE DIFFERENT STARS '
      WRITE(TITLE(6), 1310) MAG0(ICHAN), AMCO(ICHAN), RMS, COUNT
      CALL PGMTEXT ( 'T', 2.5, 0.5, 0.5, TITLE(3) )
      CALL PGSCH ( 1.0 )
      CALL PGMTEXT ( 'B', 3.5, 0.5, 0.5, TITLE(1) )
      CALL PGMTEXT ( 'L', 3.5, 0.5, 0.5, TITLE(2) )
      CALL PGMTEXT ( 'T', 1.5, 0.5, 0.5, TITLE(4) )
      CALL PGMTEXT ( 'R', 2.5, 0.5, 0.5, TITLE(5) )
      CALL PGMTEXT ( 'B', 5.0, 0.5, 0.5, TITLE(6) )
C-----------------------------------------------------------------------
C    Plot one line for each star
      IZERO = 1
      ZPLOT(NPTS+1) = -1234
      ICOLOR = 1
      DO 230 I = 2, NPTS+1
          IF ( ZPLOT(I-1) .NE. ZPLOT(I) ) THEN
              ICOLOR = MOD(ICOLOR,14) + 1
              CALL PGSCI ( ICOLOR )
              CALL PGLINE (I-IZERO, XPLOT(IZERO), YPLOT(IZERO)  )
              CALL PGPOINT(I-IZERO, XPLOT(IZERO), YPLOT(IZERO), 16 )
              IZERO = I
          END IF
  230 CONTINUE
C-----------------------------------------------------------------------
      IF ( HCOPY ) THEN
          CALL PGEND
          HCOPY = .FALSE.
          GO TO 100
      END IF
C-----------------------------------------------------------------------
C   Query the user for what to do next.
C
  300 CONTINUE
      RX = XPLOT(ICUR) 
      RY = YPLOT(ICUR) 
      CALL EGA_RESTORE_DEFAULT
  320 CONTINUE
      CALL LOCATE(0,0)
      CALL WRITE_STRING ( 'Command:' )
      CALL PGCURSE(RX, RY, CH )
      CALL CAPS ( CH )
C-----------------------------------------------------------------------
C   Ask for help.
      IF ( CH .EQ. '?' ) THEN
          CALL LOCATE  (0, 1)
          CALL WRITE_STRING ( ' ?       Displays this      ' )
          CALL LOCATE  (0, 2)
          CALL WRITE_STRING ( ' C       toggle channel     ' )
          CALL LOCATE  (0, 3)
          CALL WRITE_STRING ( ' X       toggle X axis      ' )
          CALL LOCATE  (0, 4)
          CALL WRITE_STRING ( ' H       make hardcopy      ' )
          CALL LOCATE  (0, 5)
          CALL WRITE_STRING ( ' I       identify a point   ' )
          CALL LOCATE  (0, 6)
          CALL WRITE_STRING ( ' D       delete a point     ' )
          CALL LOCATE  (0, 7)
          CALL WRITE_STRING ( ' Q       quit               ' )
          CALL LOCATE  (0, 8)
          GO TO 300
C-----------------------------------------------------------------------
C Toggle channel
      ELSE IF ( CH .EQ. 'C' ) THEN
          ICHAN = MOD ( ICHAN, 4 ) + 1
C-----------------------------------------------------------------------
C Toggle x axis
      ELSE IF ( CH .EQ. 'X' ) THEN
          XTYPE = MOD ( XTYPE, 2 ) + 1
C-----------------------------------------------------------------------
C Make a hard copy
      ELSE IF ( CH .EQ. 'H' ) THEN
          HCOPY = .TRUE.
          CALL PGEND
          GO TO 100
C-----------------------------------------------------------------------
C  Identify the plotted point closest to the cursor
      ELSE IF ( CH .EQ. 'I' ) THEN
          SCALE = (XMAX-XMIN)*(XMAX-XMIN) / ((YMAX-YMIN)*(YMAX-YMIN))
          DIST =       2.*(XMAX-XMIN)*(XMAX-XMIN)
          ICUR = 0
          DO 360 I = 1, NPTS
              DX = XPLOT(I) - RX
              DY = YPLOT(I) - RY
              DIST2 = DX*DX + SCALE*DY*DY
              IF ( DIST2 .LT. DIST ) THEN
                  DIST = DIST2
                  ICUR = I
              END IF
  360     CONTINUE
          IF ( ICUR .EQ. 0 ) THEN
              WRITE(6,*) ' INTERNAL ERROR IN DELETE ROUTINE '
          ELSE
              WRITE(6,*) ' STAR ', SLIST(ZPLOT(ICUR))
          END IF
          GO TO 300
C-----------------------------------------------------------------------
C  Delete the closest plotted point
C  Delete means mark the scan closest to the cursor as bad
      ELSE IF ( CH .EQ. 'D' ) THEN
          SCALE = (XMAX-XMIN)*(XMAX-XMIN) / ((YMAX-YMIN)*(YMAX-YMIN))
          DIST =       2.*(XMAX-XMIN)*(XMAX-XMIN)
          J = 0
          DO 400 I = 1, NPTS
              DX = XPLOT(I) - RX
              DY = YPLOT(I) - RY
              DIST2 = DX*DX + SCALE*DY*DY
              IF ( DIST2 .LT. DIST ) THEN
                  DIST = DIST2
                  J = I
              END IF
  400     CONTINUE
          IF ( J .EQ. 0 ) THEN
              WRITE(6,*) ' INTERNAL ERROR IN DELETE ROUTINE '
          ELSE
              OKSCAN(WHERE(J)) = IOR( OKSCAN(WHERE(J)),MASK(ICHAN))
              WRITE(6,*) ' DELETED STAR ', STAR(WHERE(J)),
     $                   '         SCAN ', WHERE(J)
          END IF
          GO TO 300
C-----------------------------------------------------------------------
C quit
      ELSE IF ( CH .EQ. 'Q' ) THEN
          GO TO 900
      END IF
      CALL PGASK ( .FALSE. )
      CALL PGADVANCE
      GO TO 110
  900 CONTINUE
      CALL PGEND
      RETURN
C-----------------------------------------------------------------------
 1305 FORMAT ( 'CHANNEL', I3, ' = ', A7, 20X )
 1310 FORMAT ( 'm0=', F6.2, ' AMCO=', F6.2, ' RMS=', F7.4, ' N=', I4 )
      END
