FTN77,S
$FILES 3,3
      PROGRAM OBQCK(3,80)
*************************************************************************
*
*  PROGRAM TO OBSERVE ANY PLANET IN HIGH SPEED INTERLACED MODE AT BOTH
*  3.6CM AND 6.0CM SIMULTANEOUSLY UNTIL TERMINATED BY BIT 15.
*
*  ORIGINAL STRIPPED DOWN VERSION 1994/07/20 JFHQ
*
**************************************************************************
*
      IMPLICIT NONE
*
      REAL    HALIM
      PARAMETER (HALIM = 5.9)
*
      INTEGER NFRQ, FLUX, SDEV, EPCH, MXSMPL
      PARAMETER (NFRQ = 2, FLUX = 1, SDEV = 2, EPCH = 3, MXSMPL = 650)
*
      CHARACTER PLANETS(9)*8, RESFIL*9, CFILNAM*6
*
      LOGICAL HYPFLT, EXISTS
*
      INTEGER IPLN, IRETPM(5), LU, ISWIT, ITIME(5), IYEAR, IFRQ,
     *        IALT, IERR, IPRINT, FILNAM(3)
*
      REAL    WAVE(NFRQ), HAOFS(NFRQ), DEOFS(NFRQ), BMSEP(NFRQ),
     *        TUFLUX(NFRQ), RAP, DECP, TRA, TDEC, TMPFLX, HA, TIME
*
      DOUBLE PRECISION  OFFSRS(NFRQ,3), NDIODE(NFRQ,3), ONSRS(NFRQ,3),
     *                  MEDIAN, STDDEV, EPOCH
*
      EQUIVALENCE (CFILNAM,FILNAM(1))
*
*   WAVE   is the approximate standard wavelength in CM
*   HAOFS  is the Hour Angle offset of the dual beam centre
*   DEOFS  is the Declination offset of the dual beam centre
*   BMSEP  is the Hour Angle difference between the two beams
*   TUFLUX is the rough JANSKY calibration of the noise diode
*
      DATA    WAVE    /    3.6,      6.0    /,
     *        HAOFS   /   0.144,    0.084   /,
     *        DEOFS   /  -0.603,   -0.214   /,
     *        BMSEP   /   0.254,    0.254   /,
     *        TUFLUX  /   175.0,     24.0   /
*
      DATA    PLANETS / 'MERCURY', 'VENUS', 'EARTH', 'MARS',
     *                  'JUPITER', 'SATURN', 'URANUS', 'NEPTUNE',
     *                  'PLUTO'  /
*
      DATA    RESFIL  / 'QK000::15' /,  LU  / 7 /
*
*   Pick up command line parameters (ie. from SCHDL etc )
*
      CALL RMPAR(IRETPM)
      IPLN = IRETPM(1)
      IF (IRETPM(2) .NE. 0)  LU = IRETPM(2)
*
*   Set up console to receive output
*
      CALL FSYSU(LU,LU)
      PRINT '(A1)', CHAR(14B)   ! Form Feed
*
*   Reset the Switch Register
*
      ISWIT = 0
      CALL SWOUT(ISWIT)
*
*   Print suitable program header
*
      CALL EXEC(11,ITIME,IYEAR)
      WRITE (LU,1000)  ITIME(5), IYEAR
 1000 FORMAT(//3X,"SCHEDULED PLANETS PROGRAM (OBQCK) : DAY ",I3,",",I4
     *        /3X,41"="
     *       //3X," Version  <940721.1158>"
     *       //3X," Set S register bit 15 ON to END Program.",//)
*
*   Check out Hyperbola configuration
*
      CALL CHKHYP(HYPFLT)
      IF (HYPFLT) then
         PRINT '(1X,"*** Hyperbola is Tilted !!!! ******")'
         GOTO 999
      endif
*
*   Check a suitable planet was specified and set up filename
*
      IF (IPLN.LT.1 .OR. IPLN.GT.9 .OR. IPLN.EQ.3) then
         PRINT '(1X,"**** Illegal Planet specified !!! *****")'
         GOTO 999
      endif
      PRINT '(3X,"Commencing observation of ",A)', PLANETS(IPLN)
*
*   Set up appropriate sourse and results filenames
*
      CFILNAM = PLANETS(IPLN)(1:6)
      WRITE (RESFIL(3:5),'(I3.3)') ITIME(5)
      PRINT '(3X,"Results will be in file ",A/)', RESFIL
*
*   Check Hour Angle of source is appropriate
*
      CALL PLANET(LU,FILNAM,RAP,DECP,IERR)
      IF (IERR .NE. 0) GO TO 999
      HA = TIME(2)/3600.0 - RAP/15.0
      IF (-HALIM .GT. HA .OR. HA .GT. HALIM) then
         PRINT '(1X,"**** Planet is below HA Limits !!! *****")'
         GO TO 999
      endif
*
*   Print appropriate starting header and dummy previous offsource line
*
      WRITE (LU,2000)
 2000 FORMAT(1X," MJD Epoch      Source   Wave    Off-Source      ",
     *          " Noise-diode      On-Source        Flux (approx Jy)")
      PRINT '(1X,12X,2X,A10,F5.1,17X,17X,"_")',
     *    PLANETS(IPLN), WAVE(1)
      IPRINT = 12
*
*   Repeat Observations indefinitely till S register bit 15 is set
*
      DO WHILE (ISSW(15) .GE. 0)
*
*      Observe the two frequencies sequentially
*
         DO 900 IFRQ = 1,NFRQ
            IALT = MOD(IFRQ,2) + 1
*
*         Track Planet at appropriate frequency
*
            CALL PLANET(LU,FILNAM,RAP,DECP,IERR)
            IF (IERR.NE.0) GOTO 999
            TRA = RAP + (HAOFS(IFRQ)+BMSEP(IFRQ)/2)
     *                   / COS(DECP*0.017453293)
            TDEC = DECP + DEOFS(IFRQ)
*
*         Check present Hour Angle against limits
*
            HA = TIME(2)/3600.0 - TRA/15.0
            IF (HA .GT. HALIM) then
               PRINT '(/1X,A," has now set !")', PLANETS(IPLN)
               GO TO 999
            endif
*
            CALL TRACK (4, TRA, TDEC, .TRUE., .TRUE.)
*
*         Integrate for half the time with IALT noise diode off
*
            CALL SAMPLE(IFRQ,IALT,MXSMPL/2,1)
*
*         Save results as IALT off-source measurement
*
            OFFSRS(IALT,FLUX) = MEDIAN(IALT)
            OFFSRS(IALT,SDEV) = STDDEV(IALT)
            OFFSRS(IALT,EPCH) = EPOCH(IALT)
*
*         Integrate for rest of the time with IALT noise diode on
*
            CALL SAMPLE(IFRQ,IALT,MXSMPL/2,2)
*
*         Save results as IALT noise diode measurement
*
            NDIODE(IALT,FLUX) = MEDIAN(IALT)
            NDIODE(IALT,SDEV) = STDDEV(IALT)
            NDIODE(IALT,EPCH) = EPOCH(IALT)
*
*         And save full IFRQ integration as on-source measurement
*
            ONSRS(IFRQ,FLUX) = MEDIAN(IFRQ)
            ONSRS(IFRQ,SDEV) = STDDEV(IFRQ)
            ONSRS(IFRQ,EPCH) = EPOCH(IFRQ)
*
*         Estimate Current flux from On-Source where possible
*
            TMPFLX = 0.0
            IF (OFFSRS(IFRQ,EPCH) .NE. 0.0D0) then
               TMPFLX = -1.0 * TUFLUX(IFRQ)
     *                    * (ONSRS (IFRQ,FLUX) - OFFSRS(IFRQ,FLUX))
     *                    / (NDIODE(IFRQ,FLUX) - OFFSRS(IFRQ,FLUX))
            endif
*
*         Print appropriate results line in two pieces, allowing page-throw
*
            IF (TMPFLX.EQ.0.0) then
               PRINT '(E17.8)', ONSRS(IFRQ,FLUX)
            else
               PRINT '(E17.8,F6.1)', ONSRS(IFRQ,FLUX), TMPFLX
            endif
            IPRINT = IPRINT + 1
            IF (IPRINT.GE.54) then
               PRINT '(A1)', CHAR(14B)   ! Form Feed
               WRITE(LU,2000)
               IPRINT = 0
            endif
            PRINT '(1X,F12.6,2X,A10,F5.1,E17.8,E17.8,"_")',
     *          OFFSRS(IALT,EPCH), PLANETS(IPLN), WAVE(IALT),
     *          OFFSRS(IALT,FLUX), NDIODE(IALT,FLUX)
*
*         Finally write the results to the appropriate file
*
            CALL WRITE_TO(RESFIL,PLANETS(IPLN),IFRQ,IALT,
     *                     ONSRS,NDIODE,OFFSRS,WAVE)
  
  900    CONTINUE
      END DO
*
*   Exits here on ERROR or when S register bit 15 is set, with NDs off
*
  999 CALL EXEC(2,113B,0B,1,3,3)
      PRINT '(//,3X,"*** END ***")'
      END
*
**********************************************************************
*  Check hyperbola focus and tilt matches standards for this system  *
*   returns :                                                        *
*     FAULT       logical*2   false if focus & tilt within tolerance *
*                             otherwise set true                     *
**********************************************************************
*
      SUBROUTINE CHKHYP (FAULT)
*
      IMPLICIT NONE
*
      REAL*4    STDFOCUS, STDTILT, SUBFOCUS, SUBTILT, TOLERANCE
      LOGICAL   FAULT
      CHARACTER HYPMES*11, FOCMES*11, TILTMES*9, NOTMES*4
* 
      DATA  TOLERANCE / 0.8 /, STDFOCUS / 8.0 /, STDTILT / -0.67 /, 
     &      HYPMES  /'HYPERBOLA '/,
     &      FOCMES  /'FOCUSED TO '/,
     &      TILTMES /'TILTED TO'/,
     &      NOTMES  /' NOT'/
* 
      FAULT = .FALSE. 
      CALL RHYP (SUBFOCUS, SUBTILT) 
      IF (ABS(SUBFOCUS-STDFOCUS) .GT. TOLERANCE) THEN 
          FAULT = .TRUE.
          PRINT *,HYPMES//FOCMES,SUBFOCUS,NOTMES,STDFOCUS 
      END IF
      IF (ABS(SUBTILT-STDTILT) .GT. TOLERANCE) THEN 
          FAULT = .TRUE.
          PRINT *,HYPMES//TILTMES,SUBTILT,NOTMES,STDTILT
      END IF
*
      RETURN
      END
*
*********************************************************************
*
*       FUNCTION TO OBTAIN UNIVERSAL- OR SIDERIAL TIME.
*            Returns UT if L = 1, ST if L = 2 (in SECONDS)
*                    Time in decimal year for L = 3 (less 1900.)
*
*********************************************************************
*
      REAL FUNCTION TIME(L)
*
      IMPLICIT NONE
*
      INTEGER L, IT(5), MS, IS, M, IH, ID, LSTDAY, IYR
*
      REAL    BASEOF, RATE, OFFSET
*
      EQUIVALENCE (IT(1),MS), (IT(2),IS), (IT(3),M),
     *            (IT(4),IH), (IT(5),ID)
*
      DATA  BASEOF / 84 806.882 /, LSTDAY / 0 /,
     *      RATE / 0.002 737 909 3 /
*
*   Calculate UT
*
      CALL EXEC (11, IT, IYR)
      TIME = FLOAT(IH*60 + M) * 60.0 + FLOAT(IS) + FLOAT(MS)/100.
      IF (L .LT. 2)  GO TO 900
      IF (L .EQ. 3)  GO TO 300
*
*   Calculate ST.
*
      IF (ID .NE. LSTDAY)
     *   OFFSET = FLOAT((IYR-1979)*365 + (IYR-1977)/4 + ID - 258)
     *            * 236.55536 + BASEOF + 6644.4
      TIME = TIME + TIME * RATE + OFFSET
      TIME = TIME - FLOAT(IFIX(TIME/86400.0)) * 86400.0
      GO TO 900
*
*   Calculate time in decimal YEAR
*
  300 TIME = ((TIME/86400.) + FLOAT(ID))/365. + FLOAT(IYR)
      TIME = TIME - 1900.
*
  900 RETURN
      END
*
*********************************************************************
*
*  SOURCE TRACKING SUBROUTINE.   -   Uses DRIVE as a son program
*
*********************************************************************
*
      SUBROUTINE TRACK (INTYPE, INLONG, INLAT, INTRACK, INWAIT)
*
      IMPLICIT NONE
*
      INTEGER  INTYPE
      REAL     INLONG, INLAT
      LOGICAL  INTRACK, INWAIT
*
      INTEGER*2 TYPE, IBUF(7), DRIVE(3)
      REAL*4    ALONG, ALAT 
      LOGICAL*2 TRCK, WAIT
* 
      EQUIVALENCE   (IBUF(1),TYPE),   (IBUF(2),ALONG),
     *              (IBUF(4),ALAT),   (IBUF(6),TRCK), 
     *              (IBUF(7),WAIT)
* 
      DATA    DRIVE /'DRIVE '/
* 
      TYPE = INTYPE 
      ALONG = INLONG
      ALAT = INLAT
      TRCK = INTRACK
      WAIT = INWAIT
*
*   Call DRIVE with the parameters as sepcified in the call
*
           CALL EXEC (9,DRIVE,0,0,0,0,0,IBUF,7)
*
      RETURN
      END
*
**********************************************************************
*
*   Subroutine to enact DVM sampling NSAM times into 3.6cm/6cm arrays
*
***********************************************************************
*
      SUBROUTINE SAMPLE( IFRQ, IALT, NSAM, NRUN )
*
      IMPLICIT NONE
*
      INTEGER NFRQ, FLUX, SDEV, EPCH, MXSMPL
      PARAMETER (NFRQ = 2, FLUX = 1, SDEV = 2, EPCH = 3, MXSMPL = 650)
*
      INTEGER  NSMPL(NFRQ)
*
      REAL     DVM(NFRQ,MXSMPL)
*
      DOUBLE PRECISION MJD(NFRQ,2)
*
      COMMON /SMPLER/ NSMPL, DVM, MJD
*
      INTEGER  IFRQ, IALT, NSAM, NRUN, ITIME(5), IYEAR, ISAM
*
*   Choose standard DVM select setup ie. 3.6cm on DVM1, 6cm on DVM2
*
      CALL EXEC (2,113B,31B,1,5,5)
*
*   Always reset IALT sampler and record suitable starting EPOCH
*
      NSMPL(IALT) = 0
      CALL EXEC(11,ITIME,IYEAR)
      CALL MJDAY(IYEAR,ITIME,MJD(IALT,1))
*
      IF (NRUN .EQ. 1) then
*
*      On first run set all noise diodes off
*
         CALL EXEC (2,113B,0B,1,3,3)
*
*      Also reset IFRQ sampler and record same starting EPOCH
*
         NSMPL(IFRQ) = 0
         MJD(IFRQ,1) = MJD(IALT,1)
*
       else
*
*       On second run set IALT noise diode on
*
         CALL EXEC(2,113B,2*IALT-1,1,3,3)
*
      endif
*
*   Now pause to allow the radiometers to settle
*
      CALL EXEC(12,0,2,0,-2)
*
*   Setup both DVMs, default values, medium sampling rate
*
      WRITE (35,'("*S1")')
      WRITE (36,'("*S1")')
*
*   Now sample both DVMs NSAM times into sampler arrays
*
      DO 300 ISAM = 1, NSAM
*
         NSMPL(1) = NSMPL(1) + 1
  100    READ(35,*,ERR=100) DVM(1,NSMPL(1))
         NSMPL(2) = NSMPL(2) + 1
  200    READ(36,*,ERR=200) DVM(2,NSMPL(2))
*
  300 CONTINUE
*
*   Record completion epoch for both samplers
*
      CALL EXEC(11,ITIME,IYEAR)
      CALL MJDAY(IYEAR,ITIME,MJD(IALT,2))
      MJD(IFRQ,2) = MJD(IALT,2)
*
      RETURN
      END
*
***************************************************************************
*
*  Function to return the median value of corresponding DVM samples
*
***************************************************************************
*
      DOUBLE PRECISION FUNCTION MEDIAN(IFRQ)
*
      IMPLICIT NONE
*
      INTEGER  IFRQ
*
      INTEGER NFRQ, FLUX, SDEV, EPCH, MXSMPL
      PARAMETER (NFRQ = 2, FLUX = 1, SDEV = 2, EPCH = 3, MXSMPL = 650)
*
      INTEGER  NSMPL(NFRQ)
*
      REAL     DVM(NFRQ,MXSMPL)
*
      DOUBLE PRECISION MJD(NFRQ,2)
*
      COMMON /SMPLER/ NSMPL, DVM, MJD
*
      INTEGER  IPTR(MXSMPL), ICOUNT, INXT, IFND, IUP
  
      MEDIAN = 0.0D0
      IF (NSMPL(IFRQ) .LE. 0) GO TO 999
*
*   Sort into ascending order via pointer aided bubble sort
*
      ICOUNT = 0
      DO 400 INXT = 1, NSMPL(IFRQ)
*
*      Rise through pointer list looking for larger value
*
         DO 100 IFND = 1, ICOUNT
            IF (DVM(IFRQ,INXT) .LT. DVM(IFRQ,IPTR(IFND))) GO TO 200
  100    CONTINUE
*
*         Not found ie. this value largest - add to end of list
*
            ICOUNT = ICOUNT + 1
            IPTR(ICOUNT) = INXT
         GO TO 400
*
*         Larger value found - shuffle list up to make room and insert
*
  200       ICOUNT = ICOUNT + 1
            DO 300 IUP = ICOUNT-1, IFND, -1
               IPTR(IUP+1) = IPTR(IUP)
  300       CONTINUE
            IPTR(IFND) = INXT
*        GO TO 400
*
  400 CONTINUE
*
*   Median is then simply value at the middle of sorted list
*
      MEDIAN = DVM(IFRQ,IPTR(NSMPL(IFRQ)/2))
*
  999 RETURN
      END
*
***************************************************************************
*
*  Function to return the standard deviation of corresponding DVM samples
*
***************************************************************************
*
      DOUBLE PRECISION FUNCTION STDDEV(IFRQ)
*
      IMPLICIT NONE
*
      INTEGER  IFRQ, ISAM
*
      INTEGER NFRQ, FLUX, SDEV, EPCH, MXSMPL
      PARAMETER (NFRQ = 2, FLUX = 1, SDEV = 2, EPCH = 3, MXSMPL = 650)
*
      INTEGER  NSMPL(NFRQ)
*
      REAL     DVM(NFRQ,MXSMPL)
*
      DOUBLE PRECISION MJD(NFRQ,2)
*
      COMMON /SMPLER/ NSMPL, DVM, MJD
*
      DOUBLE PRECISION  AVG, AVSQ, DIFF
*
      STDDEV = 0.0D0
      IF (NSMPL(IFRQ) .LE. 0) GO TO 999
*
      AVG = 0.0D0
      AVSQ = 0.0D0
      DO 100 ISAM = 1, NSMPL(IFRQ)
         AVG = AVG + DVM(IFRQ,ISAM)
         AVSQ = AVSQ + DVM(IFRQ,ISAM)**2
  100 CONTINUE
      AVG = AVG / NSMPL(IFRQ)
      AVSQ = AVSQ / NSMPL(IFRQ)
*
      DIFF = AVSQ - AVG**2
      IF (DIFF .LT. 0.0D0) DIFF = 0.0D0
      STDDEV = SQRT ( DIFF * 2 / NSMPL(IFRQ))
*
  999 RETURN
      END
*
***************************************************************************
*
*  Function to return the average MJD epoch of corresponding DVM samples
*
***************************************************************************
*
      DOUBLE PRECISION FUNCTION EPOCH(IFRQ)
*
      IMPLICIT NONE
*
      INTEGER  IFRQ
*
      INTEGER NFRQ, FLUX, SDEV, EPCH, MXSMPL
      PARAMETER (NFRQ = 2, FLUX = 1, SDEV = 2, EPCH = 3, MXSMPL = 650)
*
      INTEGER  NSMPL(NFRQ)
*
      REAL     DVM(NFRQ,MXSMPL)
*
      DOUBLE PRECISION MJD(NFRQ,2)
*
      COMMON /SMPLER/ NSMPL, DVM, MJD
*
      EPOCH = 0.0D0
      IF (NSMPL(IFRQ) .LE. 0) GO TO 999
*
      EPOCH = (MJD(IFRQ,1) + MJD(IFRQ,2)) / 2.0D0
*
  999 RETURN
      END
*
************************************************************************
*
*   Subroutine to output results to file each loop.  Waits on busy file
*
************************************************************************
*
      SUBROUTINE WRITE_TO( RESFIL, SRSNAM, IFRQ, IALT, ONSRS, NDIODE,
     *                     OFFSRS, WAVE )
*
      IMPLICIT NONE
*
      INTEGER NFRQ, FLUX, SDEV, EPCH, MXSMPL
      PARAMETER (NFRQ = 2, FLUX = 1, SDEV = 2, EPCH = 3, MXSMPL = 650)
*
      CHARACTER RESFIL*(*), SRSNAM*(*), JUNK*2
*
      INTEGER   IFRQ, IALT, LUfile, IERR
*
      DOUBLE PRECISION ONSRS(NFRQ,3), NDIODE(NFRQ,3), OFFSRS(NFRQ,3)
*
      REAL   WAVE(NFRQ)
*
*   Keep trying to open results file ( if IERR = 508 ie. file locked )
*
      IERR = 508
      DO WHILE (IERR .EQ. 508)
         CALL Open_File(RESFIL,'UNKNOWN','END',LUfile,IERR)
         IF (IERR .EQ. 508) PRINT '("_")'
      END DO
      IF (IERR .NE. 0) then
         PRINT '(/1X,"I/O Error ",I4," opening ",A)', IERR, RESFIL
         GO TO 999
      endif
*
*   Write appropriate lines into results file
*
  300 WRITE(LUfile,1000) ONSRS(IFRQ,EPCH), SRSNAM, WAVE(IFRQ), 'ON',
     *              ONSRS(IFRQ,FLUX), ONSRS(IFRQ,SDEV)
      WRITE(LUfile,1000) OFFSRS(IALT,EPCH), SRSNAM, WAVE(IALT), 'OF',
     *              OFFSRS(IALT,FLUX), OFFSRS(IALT,SDEV)
      WRITE(LUfile,1000) NDIODE(IALT,EPCH), SRSNAM, WAVE(IALT), 'ND',
     *              NDIODE(IALT,FLUX), NDIODE(IALT,SDEV)
 1000 FORMAT(F12.6,2X,A10,1X,F4.1,2X,A2,E17.8,' +- ',E15.8)
*
*   Close results file so it can be looked at
*
      CLOSE (UNIT=LUfile, ERR=999)
*
  999 RETURN
      END
*
*********************************************************
*
      include &RHYP::16            ! Hyperbola readout
      include &PLNET::16           ! Planet coord reader
*
*********************************************************
