FTN77
$FILES 0,1
      PROGRAM DSCN8 (3,70)
******************************************************************
*  in file /MIKE/PROGS/DSCN8.FTN on the control computer
*  version <901109.1011>
*  For drift scans through extended sources at any frequency.
*  Uses disc source file for input data
*  Run string used for LU, file name, start source number
*
*  Programs called are :
*     DPRC2        precess 1950 RA,DEC to current epoch
*     DRIVE        produce drive command for COMND
*     DSCPL        plot drift scan
*     ALARM        set servo alarm on if fault occurs
*  files INCLUDED - see end of this listing
*
      PARAMETER (MAXI = 6, MAXDATA = 480, NSPBW = 12)
      LOGICAL   FILEND, START, FAULT, GFIT, IFBRK,
     &          PLOTREPEAT, DRIVEWAIT, TRACK, NOTEND
      CHARACTER NAME*20, FILENAME*10, CLU*2, CNSORCE*3,
     &          SYSTEM(MAXI)*4
      INTEGER*2 IANG(5), ITIME(5), IPRAM(5), ISTORE(1006), ANTBUF(7),
     &          IDVMSL(MAXI), NDSEL(MAXI), ISYNTH(MAXI), DRIVE(3)
      REAL      DATA(MAXDATA), AVRG(MAXDATA),
     &          RABMOFF(MAXI), DCBMOFF(MAXI), DNOISE(6), WAVE(MAXI),
     &          HPBW(MAXI), BWFN(MAXI), RADIOMETER(MAXI)
      REAL*8    DFREQ
*
      EQUIVALENCE
*         data buffer for DPRC2
     &            (IANG(1),RAP),           (IANG(3),DECP),
*         data buffer for DSCPL
     &            (ISTORE(1),DATA(1)),     (ISTORE(961),NAME(1:1)),
     &            (ISTORE(971),WAVE(1)),   (ISTORE(983),IYR),
     &            (ISTORE(984),ITIME(1)),  (ISTORE(989),DEC),
     &            (ISTORE(991),RAS50),     (ISTORE(993),RAE50),
     &            (ISTORE(995),LAST),      (ISTORE(996),IEPOCH),
     &            (ISTORE(997),INDEX),     (ISTORE(998),LUOUT),
     &            (ISTORE(999),LP),        (ISTORE(1000),IXPOS),
     &            (ISTORE(1001),IYPOS),    (ISTORE(1002),IXC),
     &            (ISTORE(1003),IYC),      (ISTORE(1004),IHITE),
     &            (ISTORE(1005),IWIDE),    (ISTORE(1006),MKNUM),
*         data buffer for DRIVE
     &            (ANTBUF(1), ITYPE),      (ANTBUF(2), RADRIVE),
     &            (ANTBUF(4), DECDRIVE),   (ANTBUF(6), TRACK),
     &            (ANTBUF(7), DRIVEWAIT)
*
      COMMON /LUS/ LU,MTLU,LU3330,LU3336,LU8662,LU8840A,LU8840B
*
*          index           1       2       3       4       5       6
      DATA SYSTEM    /'2.5S', '3.6D', '  6S', '  6D', ' 13S', ' 18S' /,
     &     WAVE      /   2.5,    3.6,    6.1,    6.2,   13.0,   18.0 /,
     &     HPBW      / 0.066,  0.096,  0.160,  0.160,  0.333,  0.493 /,
*    exptl hpbw's              0.106           0.163   0.338        1990 jun
     &     BWFN      / 0.150,  0.450,  0.470,  0.830,  0.620,  1.250 /,
     &     RABMOFF   /-0.405, +0.154, +0.555, +0.064, +0.057, +0.056 /,
     &     DCBMOFF   /-0.050, -0.565, -0.162, -0.181, +0.445, +1.700 /,
     &     DNOISE    / 148.0,   10.0,   14.0,   1.54,   3.96,   11.2 /,
     &     IDVMSL    /   77B,    11b,    33b,    33b,    55b,    77b /,
     &     NDSEL     /     2,      1,      2,      3,      4,      5 /,
     &     RADIOMETER/  18.0,    3.6,    6.0,    6.0,   13.0,   18.0 /,
     &     ISYNTH    /  8662,   3330,   8660,   8660,   5100,   8662 /,
*
     &     LP /12/, LUDVM/35/,  IFF/14b/, FILENAME /'          '/,
     &     DTUS /239.33/, DTSS /240.0/, DTR /0.017453293/,
     &     NTSETTLE /-3/,
     &     DRIVE / 'DRIVE '/, IPAGE / 1 /
*
*  WAVE          nominal operating wavelength (cm)
*  HPBW          half power beam width (degrees)
*  BWFN          beam width to first null (degrees)
*  RABMOFF       beam offset in right ascension (degrees)
*  DCBMOFF       beam offset from axis in declination (degrees)
*  DNOISE        noise tube calibration (Kelvins)
*  IDVMSL        dvm select code for unfiltered output
*  NDSEL         calibration noise diode selection code
*  RADIOMETER    radiometer selection
*  LU,LP         I/O and Plotter logical units
*  LUDVM         lu of dvm in use (35 or 36)
*  DTUS          degrees RA to universal seconds conversion
*  DTSS          degrees  RA to sidereal  seconds    "
*  DTR           degrees to radians                  "
*  MAXDATA       data array size for drift scan
*  NSPBW         number of data points per half power beamwidth
*  NTSETTLE      noise cal settling time (seconds)
*  IPAGE         single page chart advance on plotter
*
*     Initialise the program
*     ----------------------
*
*     pickup i/o lu, filename and start source number from run string
      CALL FPARM (CLU,FILENAME,CNSORCE)
*
*     set up user terminal LU
      READ (CLU,*,ERR=20) LU
   20 IF (LU .NE. 1 .AND. LU .NE. 7) LU = 7
*     set standard input, output device for HP 1000
      CALL FSYSU (LU,LU)
      LUOUT = LU
*     do a form feed on the Anadex printer
*     IF (LU .EQ. 7) PRINT *, CHAR(14b)
*
      READ (CNSORCE,*,ERR=30) NSORCE
   30 IF (NSORCE .LT. 1) NSORCE = 1
*
*     clear the switch register on the HP 1000
      CALL ISSR(0)
      START = .TRUE.
      FAULT = .FALSE.
*
      CALL EXEC(11,ITIME,IYR)
      WRITE (LU,1000) IYR, ITIME(5)
 1000 FORMAT(  3X,'DSCN8',I5,' DAY',I3/3X,17'='//
     *         3X,'Set S reg. bit 13 on to suspend'/
     *        18X,'14 on for next source'/
     *        18X,'15 on to end.'/)
*
*         start point for each source
*         ---------------------------
*
*         Read the source parameters from the disc file
  100     CALL DSCRD (LU,START,SYSTEM,FILENAME,NSORCE,FILEND,
     *                NAME,IEPOCH,RA,DEC,SCANL,INDEX,
     *                NSCANS,GFIT,PLOTREPEAT,DFREQ)
          IF (FILEND) THEN
              GO TO 900
          END IF
          START = .FALSE.
*
*         pick up latest beam offsets from SYPARM::16
          CALL RDSYS (LU,INDEX,SYSTEM,RABMOFF,DCBMOFF)
*
*         Set up parameters for the drift scan
*         ------------------------------------
*
*         default to precessed coordinates
          RAP  = RA
          DECP = DEC
          IF (IEPOCH .EQ. 1950) THEN
*             precess source coordinates
              CALL EXEC(9,5HDPRC2,IANG(1),IANG(2),IANG(3),IANG(4))
              CALL RMPAR(IANG)
          END IF
*
*         cosine of precessed declination
          COSDECP = COS (DECP * DTR)
*         Calculate beamwidths in degrees of RA at source declination
          HPBWDEC = HPBW(INDEX) / COSDECP
          BWFNDEC = BWFN(INDEX) / COSDECP
*
*         maximum scan length
          SCANMAX = MAXDATA / NSPBW * HPBWDEC
*         compare with wanted scan length
          SCANL = MIN (SCANL, SCANMAX)
*         minimum scan length
          IF (GFIT) THEN
*             must be long enough to fit baseline for Gaussian fit
              SCANMIN = BWFNDEC + 2 * HPBWDEC
          ELSE
*             must be long enough to include first minima of beam
              SCANMIN = BWFNDEC
          END IF
          SCANL = MAX (SCANL, SCANMIN)
*
*         period for NSAM samples (seconds)
          SAMINT = HPBWDEC * DTUS / NSPBW
*
*         Find current RA and DEC with beam offsets added
          ONRA  = RAP  + RABMOFF(INDEX) / COSDECP
          ONDEC = DECP + DCBMOFF(INDEX)
*
*         Find start and end RA in current and 1950 coordinates
          HALFSC = SCANL / 2.0
          RASTRT = POSRA (RAP - HALFSC)
          RAEND  = POSRA (RAP + HALFSC)
          RAS50  = POSRA (RA - HALFSC)
          RAE50  = POSRA (RA + HALFSC)
          DELTRA = RAEND - RAE50
*
*         plot number on current page & number of scans averaged
          NPLOT = 0
          ADAV = 0.0
*
          PRINT *, CHAR(14B)
          WRITE (LU,1700)
     &    NSORCE, NAME,
     &    RA, DEC, IEPOCH, SCANL, SYSTEM(INDEX), NSCANS, IPAGE,
     &    RABMOFF(INDEX), DCBMOFF(INDEX), RAP, DECP, DNOISE(INDEX),
     &        SAMINT
 1700     FORMAT
     &    (/3X,'SOURCE',I4,' : ',A,
     &     /6X,'RA',6X,'DEC',4X,'epoch',2X,'length',2X,'system',
     &            2X,'scans',2X,'pages'
     &     /3X,2F8.3,I6,F6.2,'deg',A7,4X,i3,4X,I2,
     &    //6X,'beam offsets',7X,'RAP',5X,'DECP',2X,'noise cal.',
     &            2X,'sample period',
     &     /3X,2F8.3,3X,2F8.3,F8.2,'K',6X,F6.2,'S')
*
          IF (DFREQ .GT. 0D0) THEN
              PRINT *,'   Set LO synthesizer',ISYNTH(INDEX),
     &                ' to',DFREQ,'MHz'
              CALL SETLO (ISYNTH,INDEX,DFREQ)
          ELSE
              PRINT *,'   LO synthesizer not set'
          END IF
*
*         Start of loop for each drift scan
*         ---------------------------------
*
          DO LOOP = 1,NSCANS
              WRITE (LU,'(/3x,''SCAN '',i3,'' : _'')') LOOP
*
*             reset the "computer failed" alarm
              CALL OBSOK
*
*             plot control
              NPLOT = NPLOT + 1
              IF (NPLOT .GT. 4) NPLOT = 1
*             select DVM inputs
              INSEL = IDVMSL(INDEX)
              CALL EXEC (2,113B,INSEL,1,5,5)
*
*             Find command HA, compare to east horizon limit
  410         COMHA = TIME(2) / DTSS - RASTRT
              COMHA = PMOD(COMHA)
              WRITE (LU,4100) COMHA
 4100         FORMAT(3X,'Command HA = ',F7.3)
*
*             West limit is west horizon
*             East limit is one hour below east horizon
              IF (COMHA .GT.  HALIM(+90.0,ONDEC) .OR.
     *            COMHA .LT. (HALIM(-90.0,ONDEC) -15.0)) THEN
*                 Source will not rise in less than 1 hour - try next source
                  WRITE (LU,4120)
 4120             FORMAT(6X,'Source > 1 hour below horizon')
                  GO TO 800
              END IF
*
*             Else calculate wait time until source rises, -ve for wait
  415         ISWAIT = (COMHA - HALIM(-90.0,ONDEC)) * DTUS
              IF (ISWAIT .LT. -1 .AND. COMHA .GT. 0.0) THEN
*                 Source below eastern horizon, wait until it has risen
                  MWAIT = -ISWAIT / 60
                  WRITE (LU,4150) MWAIT
 4150             FORMAT(3X,'WAIT',I5,' M')
                  CALL EXEC(12,0,2,0,ISWAIT)
              END IF
*
*             Drive to start coordinates
*             --------------------------
*             allow for calibration time, add beam offsets
  420         RADRIVE = POSRA (RASTRT - 0.30
     *                        + RABMOFF(INDEX) / COSDECP)
              DECDRIVE = DECP + DCBMOFF(INDEX)
              ITYPE = 4
              TRACK = .FALSE.
              DRIVEWAIT = .TRUE.
*             call drive program to send command to STEER
              CALL EXEC (9,DRIVE,0,0,0,0,0,ANTBUF,7)
              CALL EXEC (9,DRIVE,0,0,0,0,0,ANTBUF,7)
*
*             Do a noise tube calibration, return dvm units/K conversion
              CALL CALT2 (LU,LUDVM,NDSEL,INDEX,NTSETTLE,DNOISE,
     &                    GASCAL,FAULT)
              IF (FAULT) THEN
*                 abort if repeated calibration failures
                  GO TO 900
              END IF
*
*             find current antenna position
              CALL ANTPOS (RELHA,RELDEC)
*
*             Western horizon check, end scan if failed
  425         IF (RELHA .GT. (HALIM(+90.0,RELDEC)-0.5))  GO TO 700
*
*             Wait for beginning of drift scan
*
*             antenna HA + Beam offset + Pointing correction
              TRUHA = RELHA + RABMOFF(INDEX)/COS(RELDEC*DTR)
     &                      + HAPC(RELHA,RELDEC)
*             true current RA of telescope
              RANOW = RATEL(TRUHA)
              WRITE (LU,4200) RELHA, RELDEC, RAS50, RAE50, RASTRT,
     &                        RAEND, RANOW
 4200         FORMAT(/4X,'RELHA',5X,'RELDEC',4X,'RAS50',5X,'RAE50',
     &                5X,'RASTRT',4X,'RAEND',5X,'RANOW'/
     &                7(F10.3))
*
              RAENDMOD = RAEND
*             360 - 0 degree RA continuity
              IF (RASTRT .GT. 180.0 .AND. RAEND .LT. 180.0)
     &            RAENDMOD = RAEND + 360.0
*
*             time to wait in seconds
              ISWAIT = -((RASTRT - RANOW) * DTUS - SAMINT / 2)
              WRITE (LU,4300) -ISWAIT
 4300         FORMAT(3X,'wait',I5,' S_')
*
*             If start time has passed, start again
              IF (ISWAIT .GT. -1) THEN
                  GO TO 420
              ELSE
                  CALL EXEC (12,0,2,0,ISWAIT)
              END IF
*
*             Start of drift scan
*             -------------------
*
              CALL EXEC (11,ITIME,IYEAR)
              WRITE (LU,'(6x,''Start scan at '',4i5,'' _'')')
     &            (itime(i) , i = 5,2,-1)
*
              DO I = 1,MAXDATA
*                 initialize data point, data counter, value sum
                  DATA(I) = 0
                  COUNTS  = 0
                  END = TIME(1) + SAMINT
                  NOTEND = .TRUE.
*
                  DO WHILE (NOTEND)
*                     read from fluke multimeter
                      READ (LUDVM,*) VALUE
                      DATA(I) = DATA(I) + VALUE
                      COUNTS = COUNTS + 1
                      NOTEND = (TIME(1) .LT. END)
                  END DO
*
*                 take mean value
                  DATA(I) = DATA(I) / COUNTS
*
*                 increment last point counter
                  LAST = I
*
*                 end data capture if end RA reached
                  RANOW = RATEL(TRUHA)
                  IF (RASTRT .GT. 180.0 .AND. RAEND .LT. 180.0 .AND.
     *                RANOW .LT. 180.0) RANOW = RANOW + 360.0
                  IF (RANOW .GE. RAENDMOD) go to 600
              END DO
*
*             If data array full, find end RA
  600         RAEND = RATEL(TRUHA) - SAMINT / (DTSS * 2)
              RAE50 = POSRA (RAEND - DELTRA)
*
*             end of drift scan - scale data to Kelvins, add to average
*             ---------------------------------------------------------
              WRITE (LU,'(6X,''End scan''/)')
*
              DO I = 1,MAXDATA
                  IF (I .GT. LAST) DATA(I) = DATA(LAST)
                  DATA(I) = DATA(I) / GASCAL
                  AVRG(I) = (AVRG(I)*ADAV + DATA(I)) / (ADAV + 1.0)
              END DO
*             increment number of averaged spectra
              ADAV = ADAV + 1.0
              IADAV = ADAV
*
*             do a gaussian fit if wanted
              IF (GFIT) THEN
                  WRITE (LU,'(6x,''this scan :'')')
                  CALL DOGFT (LU,LAST,DATA,SYSTEM,INDEX)
                  IF (IADAV .GE. 2) THEN
                      WRITE (LU,'(6x,''average scan :'')')
                      CALL DOGFT (LU,LAST,AVRG,SYSTEM,INDEX)
                  END IF
              END IF
*
*             Check switch registers and number of scans for plot options
*             -----------------------------------------------------------
*
*             Plot individual scans if number done is greater than 1,
*             or bits 14 and 15 are clear and more than 1 scan wanted
              IF (IADAV .GT. 1 .OR. (ISSW(14) .GE. 0 .AND.
     *                               ISSW(15) .GE. 0 .AND.
     *                               NSCANS .gt. 1)) THEN
*                 set position and character size for individual plots
                  IXPOS = 1000
                  IYPOS = 2125 + (NPLOT - 1) * 2250
                  IXC = 100
                  IYC = 150
                  IHITE = 2100
                  IWIDE = 4000 - IYC * 3 / 2
                  MKNUM = 3
*
*                 Call DSCPL to plot individual scan
                  IPL = 1
                  IF (NPLOT .GT. 1) IPL = 2
C                 ICLAS = 0
C                 CALL CLRQ(1,ICLAS)
C                 CALL EXEC(20+140000B,0,ISTORE,1006,0,0,ICLAS,*890)
C                 CALL EXEC(24,5HDSCPL,ICLAS,IPL,IPAGE)
                  CALL EXEC(24,5HDSCPL,IPL,IPAGE,0,0,0,ISTORE,1006)
*
*                 plot average scan
                  DO IA = 1,MAXDATA
                      DATA(IA) = AVRG(IA)
                  END DO
                  IXPOS = 5990
C                 ICLAS = 0
C                 CALL CLRQ(1,ICLAS)
C                 CALL EXEC(20+140000B,0,ISTORE,1006,0,0,ICLAS,*890)
C                 CALL EXEC(24,5HDSCPL,ICLAS,IPL,0)
                  CALL EXEC(24,5HDSCPL,IPL,0,0,0,0,ISTORE,1006)
              END IF
*
*             If next source, or end bits set, or last scan done
              IF (ISSW(14) .LT. 0 .OR. ISSW(15) .LT. 0
     *                .OR. IADAV .EQ. NSCANS) THEN
                  IF (IADAV .EQ. 1
     *                    .OR. (IADAV .GT. 1 .AND. PLOTREPEAT)) THEN
*                     plot full page average scan
                      DO IA = 1,MAXDATA
                          DATA(IA) = AVRG(IA)
                      END DO
                      IXPOS = 1500
                      IYPOS = 5500
                      IXC = 150
                      IYC = 180
                      IHITE = 8000
                      IWIDE = 8500 - IYC * 3 / 2
                      MKNUM = 7
C                     ICLAS = 0
C                     CALL CLRQ(1,ICLAS)
C                     CALL EXEC(20+140000B,0,ISTORE,1006,0,0,ICLAS,*890)
C                     CALL EXEC(24,5HDSCPL,ICLAS,3,IPAGE)
                      CALL EXEC(24,5HDSCPL,3,IPAGE,0,0,0,ISTORE,1006)
                  END IF
*                 Print average scan
                  IF ( .NOT. GFIT) THEN
                      CALL DNORM (DATA,LAST,WAVE,INDEX)
                  END IF
                  CALL TABLE (LU,LAST,SCANL,DATA)
              endif
*
*             Check for program suspension
              IF (ISSW(13) .LT. 0 .OR. IFBRK(0)) THEN
                  PAUSE 'Suspended - type *GO,DSCN8 to restart'
                  CALL ISSR (0)
                  WRITE (LU,'(''Set S-reg. options now for DSCN8'')')
                  CALL EXEC (12,0,2,0,-30)
              END IF
*
*             Check for end
              IF (ISSW(15) .LT. 0 .OR. FAULT) GO TO 900
*
*             Check for next source
              IF (ISSW(14) .LT. 0) GO TO 800
*
  700     END DO
*
*         Go to next source
*         -----------------
*
  800     NSORCE = NSORCE + 1
          GO TO 100
*
*     fault calling DSCPL
C 890 WRITE (LU,'(''error return from class write to DSCPL'')')
C     FAULT = .TRUE.
*
*     end
*     ---
*
  900 IF (FAULT) CALL EXEC(10,5HALARM)
      STOP 'Bye from DSCN8'
      END
*
*
      SUBROUTINE DSCRD (LU,START,SYSTEM,FILENAME,NSORCE,  FILEND,
     &    NAME,IEPOCH,RA,DEC,SCANL,INDEX,NSCANS,GFIT,PLOTREPEAT,DFREQ)
************************************************************************
*  Read source data from disc file
*  inputs :
*     LU         integer*2     i/o lu
*     START      logical*2     at start of program or not
*     SYSTEM     character*4   array of system designations
*     FILENAME   character*10  source file name
*     NSORCE     integer*2     number of source in file to read
*  outputs :
*     FILEND     logical*2     true if end of source file encountered
*     NAME       character*20  name identifying the source
*     IEPOCH     integer*2     year of coordinates (e.g. 1950, 1990)
*     RA         real*4        RA in decimal degrees
*     DEC        real*4        DEC in decimal degrees
*     SCANL      real*4        drift scan length in degrees of RA
*     INDEX      integer*2     internal array index identifying the system
*     NSCANS     integer*2     number of drift scans to do on this source
*     GFIT       logical*2     true for Gaussian fit, else false
*     PLOTREPEAT logical*2     true for repeat of final average plot
*     DFREQ      real*8        LO syynthesizer frequency too be set
*  File format :
*     NAME     20 character    name, in quotes
*     SYSID    4  character    system identifier, in quotes
*     IEPOCH   integer         coordinate epoch - 1950 or current year
*     RA       real            Right Ascension in decimal degrees
*     DEC      real            Declination in decimal degrees
*     SCANL    real            Wanted scan length in degrees of RA
*     NSCANS   integer         Number of drift scans to be done on the source
*     IGFIT    integer         1 to do a gaussian fit to source, else 0
*     IREPEAT  integer         1 for repeat of final average plot, else 0
*     DFREQ    real*8          L O synthesizer frequency to be set
*
      PARAMETER (MAXI = 6)
      CHARACTER FILENAME*10, NAME*20, CBUFFER*80, SYSTEM(MAXI)*4,
     &          SYSID*4, SYSTEMP*4
      LOGICAL   START, FILEND, DATAERROR, GFIT, PLOTREPEAT
      REAL*8    DFREQ
*
*     Request file info if starting
*     -----------------------------
      IF (START .AND. (FILENAME(1:1) .EQ. ' ' .OR. NSORCE .LE. 0)) THEN
  110     WRITE (lu,'(3X,''Disc file name ? _'')')
          READ  (lu,'(a)',err=110) filename(1:10)
  120     WRITE (lu,'(/3X,''Start at source number ? _'')')
          READ  (lu,*,err=120) nsorce
      END IF
*
*     defaults :
      NSCANS  = 1
      IGFIT   = 0
      IREPEAT = 0
      FILEND  = .FALSE.
      DFREQ = 0D0
*
      OPEN (99,FILE=FILENAME,IOSTAT=IOS,ERR=700)
*
  200 REWIND (99)
      DO 210 I = 1,NSORCE
          READ (99,'(A)',IOSTAT=IOS,ERR=800,END=900) CBUFFER
  210 END DO
*
*     truncate read at last parameter in list in character buffer
      CBUFFER(79:79) = '/'
      READ (CBUFFER,*,IOSTAT=IOS,ERR=800,END=300)
     &    NAME,SYSID,IEPOCH,RA,DEC,SCANL,NSCANS,IGFIT,IREPEAT,DFREQ
*
      CLOSE (99)
*
*     Check correctable data
  300 IF (DEC .GT. 270.0) DEC = DEC - 360.0
      IF (NSCANS .LT. 1) NSCANS = 1
      GFIT = .FALSE.
      IF (IGFIT .EQ. 1) GFIT = .TRUE.
      PLOTREPEAT = .FALSE.
      IF (IREPEAT .EQ. 1) PLOTREPEAT = .TRUE.
      CALL UPCASE (NAME)
      CALL UPCASE (SYSID)
*
*     Check for fatal data errors
      DATAERROR = .FALSE.
*
*     right justify the input system index
      DO WHILE (SYSID(4:4) .EQ. ' ')
          SYSTEMP = '    '
          SYSTEMP(2:4) = SYSID(1:3)
          SYSID = SYSTEMP
      END DO
*
*     find the system index
      INDEX = 0
      DO I = 1, MAXI
          IF (SYSTEM(I) .EQ. SYSID) INDEX = I
      END DO
      IF (INDEX .EQ. 0) THEN
          DATAERROR = .TRUE.
      END IF
*
  320 IF (IEPOCH .NE. 1950 .AND. IEPOCH .LT. 1990
     *    .OR. RA .LT. 0.0 .OR. RA .GT. 360.0
     *    .OR. DEC .LT. -90.0 .OR. DEC .GT. 45.0
     *    .OR. SCANL .LE. 0.0)
     *        DATAERROR = .TRUE.
*
      IF ( .NOT. DATAERROR) THEN
          RETURN
      ELSE
*         Fatal error in data
          GO TO 800
      END IF
*
*         error opening file
  700     WRITE (lu,'(''Error '',i4,'' opening file '',a10)')
     *        ios, filename
          STOP 'BYE'
*
*         read error encountered, try next source
  800     WRITE (lu,'(''Error'',i5,'' at source'',i3/a)')
     &            ios, nsorce, cbuffer(1:80)
          NSORCE = NSORCE + 1
          GO TO 200
*
*         end of file
  900     WRITE (lu,'(''End of file at source'',i4)') nsorce
          FILEND = .TRUE.
          RETURN
          END
*
      SUBROUTINE DOGFT (LU,LAST,ARRAY,SYSTEM,INDEX)
************************************************************************
*     set up Gaussian fit simulating the beam
* inputs :
*     LU      integer*2       i/o lu
*     LAST    integer*2       last valid point in the array
*     ARRAY   real*4          data array
*     SYSTEM  character*4     array of system identifiers
*     INDEX   integer*2       index identifying system in use
*
      PARAMETER  (MAXI = 6)
      LOGICAL*2  DUALBEAM
      REAL*4     ARRAY(1), FLUXESTA(11), FLUXESTB(11)
      CHARACTER  SYSTEM(MAXI)*4
*
*     subtract linear fit to baseline of array
      CALL LNOFF (LAST,ARRAY,RMS)
      WRITE (lu,'(9x,''rms noise = '',1pg12.6,'' K'')') RMS
*
*     define the mid beam points in  the data array
      IF (SYSTEM(INDEX) .EQ. '  6D') THEN
          MIDBEAMA = LAST / 2 - 8
          MIDBEAMB = LAST / 2 + 8
          DUALBEAM = .TRUE.
      ELSE IF (SYSTEM(INDEX) .EQ. '3.6D') THEN
          MIDBEAMA = LAST / 2 - 13
          MIDBEAMB = LAST / 2 + 13
          DUALBEAM = .TRUE.
      ELSE
*         for single feeds
          MIDBEAM = LAST / 2
          DUALBEAM = .FALSE.
      END IF
*
      IF (DUALBEAM) THEN
*         make 11 estimates of flux for each beam
*         error in flux is rms/2 per beam owing to gaussian fit
          J = 0
          DO I = MIDBEAMA - 5, MIDBEAMA + 5
              J = J + 1
              CALL GFLUX (ARRAY,I,FLUXESTA(J))
          END DO
          J = 0
          DO I = MIDBEAMB - 5, MIDBEAMB + 5
              J = J + 1
              CALL GFLUX (ARRAY,I,FLUXESTB(J))
          END DO
*
*         find (signed) maximum flux in each beam
          FLUXA = FLUXESTA(1)
          FLUXB = FLUXESTB(1)
          DO J = 2, 11
              FLUXA = MAX (FLUXA, FLUXESTA(J))
              FLUXB = MIN (FLUXB, FLUXESTB(J))
          END DO
          FLUX = (FLUXA - FLUXB) / 2
          FLUXERROR = RMS / ( 2 * 1.414 )
          WRITE (lu,'(9x,''A beam flux = '',1pg12.6,''K'',
     &                3x,''B beam flux = '',1pg12.6,''K''/
     &                9x,''average flux= '',1pg12.6,'' +- _'')')
     &                FLUXA, FLUXB, FLUX
*
      ELSE
*         single beam flux fit
          J = 0
          DO I = MIDBEAM - 5, MIDBEAM + 5
              J = J + 1
              CALL GFLUX (ARRAY,I,FLUXESTA(J))
          END DO
          FLUX = FLUXESTA(1)
          DO J = 2, 11
              FLUX = MAX (FLUX,FLUXESTA(J))
          END DO
          FLUXERROR = RMS / 2
          WRITE (lu,'(9x,''flux = '',1pg12.6,'' +- _'')') FLUX
      END IF
*
      WRITE (lu,'(1pg12.6,'' K'')') FLUXERROR
      IF (FLUX .LT. 5 * FLUXERROR) THEN
          WRITE (lu,'(9x,''flux less than 5 * flux error :''/
     &        3x,''A beam flux estimates (K) :''/
     &        3x,6(1pg10.3,2x) / 3x,5(1pg10.3,2x))')
     &        (FLUXESTA(I), I = 1,11)
          IF (DUALBEAM) THEN
              WRITE (lu,'(3x,''B beam flux estimates (K) :''/
     &            3x,6(1pg10.3,2x) / 3x,5(1pg10.3,2x))')
     &            (FLUXESTB(I), I = 1,11)
          END IF
      END IF
      RETURN
      END
*
      SUBROUTINE LNOFF (N,DATA,RMS)
************************************************************************
*     remove slope between start and end points of data in array
* inputs :
*     N       integer*2       last valid point in array
*     DATA    real*4          data array
* outputs :
*     DATA    real*4          data array with slope removed
*     RMS     real*4          RMS scatter in first ten and last ten points
*
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*4     DATA(N), RMS
      INTEGER*2  ISTART(2), IEND(2)
*
*     use first ten and last ten points of data for slope
*     linear regression :
      MRANGE = 10
      SX  = 0.0
      SXX = 0.0
      SY  = 0.0
      SXY = 0.0
      DO K = 1,2
          IF (K .EQ. 1) THEN
              ISTART(K) = 1
              IEND(K)   = ISTART(K) + MRANGE - 1
          ELSE
              ISTART(K) = N - MRANGE + 1
              IEND(K)   = N
          END IF
*
          DO M = ISTART(K), IEND(K)
              X   = M
              SX  = SX + X
              SXX = SXX + X**2
              SY  = SY + DATA(M)
              SXY = SXY + X * DATA(M)
          END DO
      END DO
*
*     regression coefficients :
      M2 = MRANGE * 2
      SLOPE = ( SXY - SX*SY/M2 ) / ( SXX - SX*SX/M2 )
      CONSTANT = SY/M2 - SLOPE*SX/M2
*     subtract linear drift :
      DO I = 1, N
          DATA(I) = DATA(I) - ( SLOPE*I + CONSTANT )
      END DO
*
*     find rms noise over fitted region :
      VAR = 0.0
      DO K = 1,2
          DO M = ISTART(K), IEND(K)
              VAR = VAR + DATA(M) * DATA(M)
          END DO
      END DO
      RMS = SQRT (VAR/M2)
      RETURN
      END
*
      SUBROUTINE GFLUX (DATA,MIDBEAM,FLUX)
******************************************
*     calculate flux by fitting a gaussian
* inputs :
*     DATA        real*4      data array
*     MIDBEAM     integer*2   index of mid beam = center of Gaussian
* output :
*     FLUX        real*4      flux calculated from Gaussian fit
*
*     number of samples per half power beamwidth
      PARAMETER (NSPBW = 12)
*
      REAL DATA(1)
*
      IRANGE = NSPBW / 2
      SFLUX = 0.0
      SBMSQ = 0.0
      DO I = MIDBEAM-IRANGE, MIDBEAM+IRANGE
*         offset from center of beam in sampling units
          OFFSET  = FLOAT (I - MIDBEAM) / FLOAT (IRANGE)
*         normalised beam height at offset - use eqn for Gaussian
          BEAM = EXP (-0.693 * OFFSET**2)
          SFLUX = SFLUX + BEAM * DATA(I)
          SBMSQ = SBMSQ + BEAM**2
      END DO
*     flux in data input units
      FLUX = SFLUX / SBMSQ
      RETURN
      END
************************
      INCLUDE DSSUB.FTN
      INCLUDE HALIM.FTN
      INCLUDE RDSYS.FTN
      INCLUDE UPCASE.FTN
      INCLUDE COMMONLUS.SUB
      INCLUDE LOSET.SUB
      INCLUDE OBSOK.SUB
************************
