FTN66 
$ALIAS /DOUT1/,NOALLOCATE 
$ALIAS /POINT/,NOALLOCATE 
      PROGRAM OBCIR(3,80),  Multi-freq Cir X-1 CR16 <990122.1505> 
C********************************************************************** 
C                                                                     * 
C     PROGRAM TO OBSERVE CIRCINUS X-1 AT MULTIPLE FREQUENCIES         * 
C     ADAPTED FROM PROGRAM OBCIR                                      * 
C     WRITTEN BY GDN 1980 06 27                                       * 
C     NEW 3.6 CM BACKGROUND FOR WIDER FEED SEPARATION  1994-12-07     * 
C     CORRECTED BUG RATEL CALCULATION   AFFECTING ONLY CALIBRATOR 
C                                       1994-12-14                    * 
*********************************************************************** 
C 
C 
      DIMENSION SCAN(4,100), GDVM(2), IOBUF(20), BACRND(96),
     *          BCRND1(96), BCRND2(90), BCRND3(90), SAMPLE(200) 
      DIMENSION NAMDP(3),IANG(5),FLUX(3),IHEAD(8),IPAGE(2),IWVLN(2) 
      DIMENSION FREQ(3), HASTRT(3), DECOFF(3), NODVM(3), NOSCAN(3), 
     *          NTUBE(3), NOMNT(3), HFSCAN(3), NREAD(3), CALJY(3),
     *          KSTART(3), NSAMPL(3), IHBMPT(3), IBEAMA(3),RDVM(400), 
     *          IBMSEP(3), NBMPTS(3), BMWDTH(3), THRMS(3), RAOFF(3),
     *          IPRAM(5), RACAL(5), DECCAL(5) 
C 
      INTEGER IDVM1 
      INTEGER NWAKE,CMD(12),ITIME(5)
C 
      LOGICAL HASTOP,DECSTP 
C 
      REAL ATIME(5),COORDS(6) 
C 
      EQUIVALENCE (CMD(1),RA0),(CMD(3),DEC0), 
     *(CMD(5),IRAT),(CMD(6),STOPTM),
     *(CMD(10),HASTOP),(CMD(11),DECSTP),
     *(CMD(12),NWAKE),(COORDS(1),RHA),
     *(COORDS(2),RDEC), 
     *(IANG(1),RA50),(IANG(3),DEC50)
C 
      COMMON/POINT/LOGCLS,NCCMD,J,COMCO 
      COMMON/DOUT1/ISERVO 
C 
      DATA IRAT/4/
      DATA NAMDP/2HDP,2HRC,2H2 /,IHEAD/2HCI,2HRC,2HIN,2HUS,2H X,
     *2H-1,2H D,2HAY/,IPAGE/2HPA,2HGE/,IWVLN/2H C,2HM. /
* 
* 
      DATA BCRND1/0.056,0.104,0.107,0.113,0.090,0.083,0.067,0.075,
     *            0.053,0.056,0.041,0.037,0.030,0.044,0.051,0.035,
     *            0.014,-.003,-.005,-.006,-.021,-.023,-.029,-.017,
     *            -.011,-.004,-.011,-.031,-.064,-.077,-.092,-.112,
     *            -.153,-.186,-.209,-.223,-.229,-.225,-.198,-.175,
     *            -.137,-.099,-.068,-.067,-.066,-.073,-.068,-.085,
     *            -.082,-.066,-.043,-.015,0.006,0.037,0.053,0.060,
     *            0.078,0.074,0.081,0.056,0.042,-.003,-.049,-.081,
     *            -.112,-.132,-.169,-.208,-.260,-.287,-.307,-.312,
     *            -.323,-.290,-.245,-.210,-.171,-.128,-.074,-.044,
     *            -.026,-.010,0.003,0.013,0.022,0.034,0.053,0.072,
     *            0.077,0.081,0.072,0.058,0.042,-.000,-.059,-.083/
* 
      DATA BCRND2/-.14,-.06,0.04,0.10,0.16,0.28,0.38,0.44,0.48,0.52,
     *            0.54,0.58,0.58,0.56,0.52,0.48,0.46,0.42,0.38,0.32,
     *            0.28,0.24,0.20,0.16,0.14,0.12,0.12,0.10,0.08,0.08,
     *            0.00,-.02,-.08,-.12,-.18,-.22,-.24,-.28,-.30,-.28,
     *            -.28,-.24,-.20,-.16,-.10,-.01,0.12,0.21,0.28,0.35,
     *            0.41,0.42,0.41,0.38,0.31,0.26,0.20,0.13,0.06,0.04,
     *            0.02,0.04,0.08,0.14,0.22,0.32,0.42,0.52,0.60,0.62,
     *            0.60,0.56,0.52,0.44,0.38,0.30,0.20,0.12,0.08,0.04,
     *            0.02,0.02,0.06,0.14,0.22,0.34,0.44,0.62,0.76,0.88 / 
C 
* 
      DATA BCRND3/0.12,0.17,0.22,.028,0.32,0.35,0.38,0.37,0.37,0.38,
     *            0.36,0.33,0.32,0.27,0.22,0.17,0.14,0.12,0.07,0.04,
     *            0.01,0.00,0.00,0.01,0.02,0.04,0.08,0.14,0.20,0.26,
     *            0.32,0.40,0.50,0.62,0.75,0.86,0.98,1.07,1.17,1.32,
     *            1.44,1.53,1.70,1.84,2.00,2.04,2.10,2.12,2.11,2.11,
     *            2.09,2.08,2.05,2.01,1.96,1.87,1.82,1.72,1.60,1.48,
     *            1.40,1.29,1.21,1.13,1.07,0.97,0.91,0.80,0.75,0.68,
     *            0.60,0.50,0.44,0.36,0.33,0.27,0.23,0.18,0.12,0.09,
     *            0.04,-.02,-.05,-.09,-.14,-.18,-.21,-.27,-.31,-.24 / 
* 
      DATA  FREQ    /  3.6 ,    6.25,      13.1   /,
     *      HASTRT  /  1.082,   1.85,       1.85  /,
     *      RAOFF   /  0.184,   0.208,      0.0   /,
     *      DECOFF  / -0.556,  -0.220,      0.423 /,
     *      NODVM   / 11B   ,  33B   ,     55B    /,
     *      NOSCAN  /  1    ,   1    ,      1     /,
     *      NTUBE   /  1    ,   3    ,      4     /,
     *      NOMNT   / 5600  ,  1150  ,    2000    /,
     *      HFSCAN  /  1.00 ,   1.50 ,      1.50  /,
     *      NREAD   /  5    ,    8   ,       8    /,
     *      CALJY   / 175.0 ,   21.7 ,      38.5  /,
     *      KSTART  /  71   ,   65   ,      65    /,
     *      NSAMPL  /  96   ,   90   ,      90    /,
     *      IBEAMA  /  33   ,   35   ,      35    /,
     *      IBMSEP  /  22   ,   14   ,      14    /,
     *      IHBMPT  /   4   ,    4   ,       4    /,
     *      NBMPTS  /   8   ,    9   ,       9    /,
     *      BMWDTH  / 0.0433,  0.0357,    0.1428  /,
     *      THRMS   /   450.,    150.,      200.  / 
* 
* 
      DATA T0 /299.37/, P0 /16.5768/, PDOT /-0.0000352/ 
* 
* 
* 
C     INITIALISE PARAMETERS 
C     ********************* 
* 
* 
*      READ PARAMETERS PASSED IN RUN COMMAND
* 
       CALL RMPAR (IPRAM) 
* 
*      SET SWITCH REGISTER
* 
       CALL SWOUT (IPRAM(1))
* 
*     Reset computer alarm
* 
C     CALL OBSOK
* 
      STOPTM=-1.
  1   CALL EXEC(11,ITIME,IYR) 
      UT = (ITIME(2)/3600. + ITIME(3)/60. + ITIME(4)) 
      ST = UT*1.002738  + 8.46022  + ITIME(5)*0.0657097 
 2999 IF(ST.LT.24.0) GO TO 3000 
      ST = ST - 24.0
      GO TO 2999
 3000 WRITE (7,100) ITIME(5),IYR,ST 
 100  FORMAT (/////,10X,"CIRCINUS X-1 OBSERVATIONS DAY ",I3,2X,I4,
     *       "   ST  =  ",F5.2,"Hrs") 
* 
*     TEST WHETHER SOURCE IS ABOVE HORIZON  
*     *********************************** 
 3010 SORSUP =  21.0 - ST 
C     IF(SORSUP.LT.0.0.OR.SORSUP.GT.12.0) GO TO 3050
* 
* 
* 
      RA50 =  229.201 
      DEC50 = 303.013 
      IF (IPRAM(5).EQ.2) THEN 
         RA50 =  263.366
         DEC50 = 303.468
         ELSE 
         CONTINUE 
         END IF 
      IPAGNO = 1
C     PRECESS COORDINATES TO CURRENT EPOCH USING DPRC2
C     ************************************************
      CALL EXEC (9,NAMDP,IANG(1),IANG(2),IANG(3),IANG(4),IANG(5)) 
      CALL RMPAR(IANG)
      RANOW = RA50
      WRITE (7,110) RA50,DEC50
 110  FORMAT (/,10X"CURRENT COORDINATES ARE RA = ",F8.3," DEC = ",F8.3,/
     1 10X,52"*",////,
/    2 17X,       "NOISE WAIT TRANSIT  A BEAM  B BEAM AVERAGE    RMS"/
     3 10X," H.A.   CAL  TIME  TIME     FLUX    FLUX   FLUX     NOISE 
     4WAVE   DEC          PHASE CYCLE"/   
     5 10X," ****   ***  ****  ****     ****    ****   ****     *****") 
C 
C 
C     ENTER MAIN  LOOP OF OBSERVING PROGRAM TO SELECT FREQUENCY 
C     **********************************************************
C 
 112  DO 1000 IFREQ = 1,3 
* 
*     Check SWITCH REGISTER positions..... 0 skip 3.6 cm
*                                          1 skip 6.3 cm
*                                          2 skip 13. cm
* 
* 
*     Check if all switches have been set and if so write suitable
*     message.
* 
* 
      IF (ISSW(0).LT.0.AND.ISSW(1).LT.0.AND.ISSW(2).LT.0) GO TO 108 
          GO TO 109 
* 
 108              WRITE (7,115) 
 115               FORMAT ("At least ONE wavelength must be in use! ",
     *                     " RESET switch register as follows..."//,
     *                     "  BIT 0 set will skip 3.6 cm"/, 
     *                     "  BIT 1 set will skip 6.3 cm"/, 
     *                     "  BIT 2 set will skip 13  cm"/, 
     *                     " Program will suspend for 30 seconds "//) 
* 
* 
      CALL EXEC (12,0,2,0,-30)
* 
 109  IF (ISSW(0).LT.0.AND.IFREQ.EQ.1) GO TO 1000 
* 
         IF (ISSW(1).LT.0.AND.IFREQ.EQ.2) GO TO 1000
* 
            IF (ISSW(2).LT.0.AND.IFREQ.EQ.3) GO TO 1000 
C 
C     Apply offset corrections to source coordinates
C 
*     Set declination offset
* 
      DCOFF = DECOFF(IFREQ) 
* 
*     Check whether day time dec correction is needed 
* 
C     CALL EXEC (11,ITIME, IYR) 
* 
* 
C          IF (ITIME(4).GE.07.AND.ITIME(4).LT.15) 
                  DCOFF = DECOFF(IFREQ) 
      RA0 = RA50 - HASTRT(IFREQ)
      DEC0 = DEC50 + DCOFF
      DECON = DEC0
C 
C     Select correct radiometer for current frequency 
C 
      NDVM = NODVM(IFREQ) 
      CALL EXEC (2,113B,NDVM,1,5,5) 
C 
C 
C     Drive telescope to start of scan
C     ********************************* 
 10   DUMMY = 0.
      IF (IPRAM(3).GT.0) THEN 
         NSCAN = IPRAM(3) 
         ELSE 
         NSCAN = NOSCAN(IFREQ)
         END IF 
      NBAD = 0
      DO 660 N= 1,NSCAN 
      DEC0 = DECON - 0.010*IPRAM(4)*((IPRAM(3)-1)*0.5 - (N-1))  
      AN = N      
      HASTOP = .TRUE. 
      DECSTP = .TRUE. 
C 
C     first drive call
      CALL RNRQ (12B,NWAKE,ISTAT) 
      NCMD = NCCMD
      CALL EXEC (20,0,CMD,12,0,0,NCMD)
      CALL RNRQ (5,NWAKE,ISTAT) 
      CALL RNRQ (40B,NWAKE,ISTAT) 
C 
C     DRIVE AGAIN IMMEDIATELY 
      CALL RNRQ (12B,NWAKE,ISTAT) 
      NCMD = NCCMD
      CALL EXEC (20,0,CMD,12,0,0,NCMD)
      CALL RNRQ (5,NWAKE,ISTAT) 
      CALL RNRQ (40B,NWAKE,ISTAT) 
C 
C     DO NOISE TUBE CALIBRATION 
C     ************************* 
C 
      CALL EXEC (2,113B,NDVM,1,5,5) 
      CALL EXEC (12,0,2,0,-1) 
 310  GDVM(1) =0.0
      GDVM(2) =0.0
      NTRY = 0
C        Read Fluke 8840 20 times per second
  309 DO 350 K = 1,2
         CALL EXEC (2,113B,NDVM,1,5,5)
         WRITE (35,'("*S1")') 
      DO 330 L = 1,400
         READ  (35,*) RDVM(L) 
 330     CONTINUE 
         DO 333 NZ = 1,400  
            GDVM(K) =  (GDVM(K) + RDVM(NZ)) 
 333        CONTINUE
            GDVM(K) = -GDVM(K)/400.   
            IF (K.EQ.2) GO TO 350 
C     FIRE NOISE TUBE AND WAIT TO SETTLE
      NTNO = NTUBE(IFREQ) 
      CALL EXEC (2,113B,NTNO,1,3,3) 
      CALL EXEC (12,0,2,0,-5) 
 350  CONTINUE
C     TURN OFF NOISE TUBE 
      CALL EXEC (2,113B,0,1,3,3)
      GASCAL = ABS((GDVM(2) - GDVM(1))*100000)
C     IF(GASCAL.GT.300.) GO TO 365
C     GASCAL = NOMNT (IFREQ)
 365  NSECAL = GASCAL*0.1 
C     WRITE (7,5000) GDVM(1), GDVM(2), NSECAL 
C5000 FORMAT( "Off = ",F8.2, " On = ",F8.2," Noise cal = ",I5)
* 
*     Read Hour Angle 
* 
      CALL EXEC (1,11,IOBUF(1),10,1,2)
*   
*     Mask off tachometer bits
* 
      IOBUF(1) = IOBUF(1).AND.171777B 
* 
      CALL CODE 
         READ (IOBUF,367) HA
 367     FORMAT (F6.3)
* 
      WRITE (7,335) HA,NSECAL 
 335  FORMAT (  8X,F7.3,1X,I5 "_")
C     SAMPLE DVM FOR 12 MINUTES, PUTTING MARKER ON AFTER 6 MINUTES
C     ************************************************************
C     WAIT FOR BEGINNING OF SCAN
      CALL EXEC (11,ITIME,IYR)
      YR = IYR
      HRS = ITIME(4)
      AMINS = ITIME(3)
      SECS =  ITIME(2)
      AMSECS = ITIME(1) 
      UT = HRS + (AMINS +  SECS/60. + AMSECS/6000.)/60. 
      DAY = ITIME(5)
      STO = 8.46022  + DAY*0.0657097  
      IF (STO.GE.24.) STO = STO - 24. 
      ST = STO + UT*1.002738
      IF (ST.GE.24.) ST = ST - 24.
      IF (HA.GT.270.) HA = HA - 360.
         HAIN = HA
         DECIN = DEC0 
           CALL PTCOR (HAIN, DECIN, DHA, DDEC)
           DHA = DHA/1000.
      RATEL = ST*15 - HA - DHA -0.040 
      IF (RATEL.LT.0.0) RATEL = RATEL + 360 
      IF (RATEL.GT.360.0) RATEL = RATEL - 360.0 
      IWAIT = (RATEL - RANOW + HFSCAN(IFREQ) - RAOFF(IFREQ))*240  
      WRITE (7,370) IWAIT 
 370  FORMAT (I4"_")
      IF(IWAIT.GE.0) IWAIT = -1 
      CALL EXEC (12,0,2,0,IWAIT)
      HAEND = 87.95 
* 
*     If source has set then test whether program should restart
*     when source next rises
* 
      IF (HA.GT.HAEND.AND.ISSW(14).LT.0) GO TO 3050 
* 
* 
*     Test whether source has set and if so terminate program 
* 
 369  IF(HA.LE.HAEND) GO TO 373 
 3050 SUSTIM = 9.0 - ST 
      IF(SUSTIM.LT.0) SUSTIM = SUSTIM + 24.0
      WRITE (7,371)  SUSTIM 
 371  FORMAT(/ 4X,"SOURCE BEYOND LIMITS. RISES AGAIN IN ",F6.2, 
     * " Hrs."//) 
* 
* 
      GO TO 1000
 373  DO 550 K = 1,2
        IEND = NSAMPL(IFREQ)
        JEND = IEND/2 
        KREAD = (NREAD(IFREQ))*20 - 2 
        DO 545 L = 1,JEND 
         CALL EXEC (2,113B,NDVM,1,5,5)
            WRITE (35,'("*S1")')
         DO 540 M = 1,KREAD 
            READ  (35,*) SAMPLE(M)
 540     CONTINUE 
* 
*   Sort data into ascending order
* 
         NS = KREAD 
          MS = NS 
  541     MS = MS/2 
          IF(MS.EQ.0) GO TO 1544
             KS = NS - MS 
             JS = 1 
  542        IS = JS
  543        LS = IS + MS 
               IF(SAMPLE(IS).LE.SAMPLE(LS)) GO TO 544 
                  BS = SAMPLE(IS) 
                  SAMPLE(IS) = SAMPLE(LS) 
                  SAMPLE(LS) = BS 
                  IS = IS - MS
                    IF(IS.GE.1) GO TO 543 
  544                 JS = JS + 1 
                      IF(JS.LE.KS) GO TO 542
                        GO TO 541 
* 
*     Sort complete 
* 
 1544    L1 = (K-1)*JEND + L
         SCAN (1,L1) = -(SAMPLE(KREAD/2))*100000
* 
 545     CONTINUE 
         IF (K.EQ.2) GO TO 550
         CALL EXEC (11,ITIME,IYR) 
         WRITE (7,548) ITIME(4),ITIME(3),ITIME(2) 
 548     FORMAT (1X,3I3"_") 
 550     CONTINUE 
* 
* 
*     CALCULATE PHASE AT SOURCE TRANSIT 
* 
* 
      NOYRS = IYR - 1976
      NODAYS = NOYRS*365 + (IYR -1973)/4 + ITIME(5) 
      TRNSTM = (ITIME(4))/24.0 + (ITIME(3))/(60.0*24.0) + NODAYS
* 
*     CALCULATE NO OF CYCLES FROM DAY T0, 1976
* 
      DO 552 NCYCL = 1,660
         EPOCH = T0 + (P0 + PDOT*NCYCL)*NCYCL 
            IF (EPOCH.GT.TRNSTM)  GO TO 553 
  552         CONTINUE
       NCYCL = NCYCL - 1
  553  PERIOD = P0 + 2*PDOT*NCYCL 
       PHASE = (TRNSTM - (EPOCH - PERIOD))/PERIOD 
* 
       DO 555 M = 1,10
      N1 = (M-1)*9 + 1
 555  CONTINUE
C     SCALE SCANS TO JANSKYS AND SUBTRACT BACKGROUND
C     **********************************************
C 
      SCALE = CALJY(IFREQ)/GASCAL 
      DO 602 L = 1,IEND 
C 
C     SELECT CORRECT BACKGROUND 
C 
         IF(IFREQ.EQ.1) BACRND(L) = BCRND1(L) 
         IF(IFREQ.EQ.2) BACRND(L) = BCRND2(L) 
         IF(IFREQ.EQ.3) BACRND(L) = BCRND3(L) 
C 
C 
      TEMP = -SCAN(1,L) 
      SCAN(1,L) = SCALE*TEMP
*     
*     SKIP BACKGROUND SUBTRACTION IF IPRAM(5) GREATER THAN  0 
* 
      IF (IPRAM(5).GT.0) THEN 
         SCAN(2,L) = SCAN(1,L)
         ELSE 
         SCAN(2,L) = SCAN(1,L)  - BACRND(L) 
         END IF 
602   CONTINUE
C     REMOVE LINEAR DRIFT FROM BACKGROUND SUBTRACTED SCAN 
C     *************************************************** 
      ISTART = 1
      SX  = 0.0 
      SXX = 0.0 
      SY  = 0.0 
      SXY = 0.0 
      DO 604 K=1,2
      DO 603 M = ISTART, ISTART+24
      SX =  SX + M
      SXX = SXX+ M*M
      SY  = SY + SCAN(2,M)
      SXY = SXY + M*SCAN(2,M) 
 603  CONTINUE
      ISTART = ISTART + KSTART(IFREQ) 
 604  CONTINUE
      A = (SXY - SX*SY/50.)/(SXX - SX*SX/50.) 
      B = SY/50. - A*(SX/50.) 
C     SUBTRACT OFF DRIFT FROM RAW DATA, AND BACKGROUND CORRECTED SCAN 
C     *************************************************************** 
      DO 605 L = 1,IEND 
      TEMP1 = SCAN(1,L) 
      TEMP2 = SCAN(2,L) 
      DRIFT = A*L + B 
      SCAN(1,L) = TEMP1 - DRIFT 
      SCAN(2,L) = TEMP2 - DRIFT 
 605  CONTINUE
C     AVERAGE SCANS 
C     ************* 
* 
*     Compute rms noise in each scan
* 
      ISTART = 1
      VAR = 0.
      DO 625 K = 1,2
        DO 620 M = ISTART, ISTART+24
          VAR = VAR + (SCAN(2,M))**2
 620      CONTINUE
        ISTART = KSTART(IFREQ)
 625    CONTINUE
      RMS = 1000.*SQRT(VAR/48.) 
          IF(RMS.LT.THRMS(IFREQ)) GO TO 609 
             NBAD = NBAD + 1
                GO TO 609 
* 
*     Average scans 
* 
* 
 609  DO 608 L=1,IEND 
          SCAN(3,L) = (SCAN(3,L)*(N-1-NBAD) + SCAN(2,L))/(N-NBAD) 
 608  CONTINUE
* 
* 
C     PLOT SCANS ON SCALE 1 JY = 1 INCH 
C     ********************************* 
C 
C         WRITE HEADINGS ON PLOTTER before plotting first scan
C         *************************         MJG 1987 sept 14
C         do a chart advance before first plot on page
C ****
******    OLD PLOTTING CODE DISCARDED HERE
******************************************* 
******
C     CALCULATE  FLUX DENSITY FOR A AND B BEAMS, AND AVERAGE FLUX 
C     *********************************************************** 
C 
      IF(IFREQ.EQ.3) GO TO 800
* 
* 
      DO 720 M=1,2
      SBEAM = 0.
      SBMSQ = 0.
      SFLUX = 0.
      K1 = IBEAMA(IFREQ) + IBMSEP(IFREQ)*(M-1)
      K2 = K1 + NBMPTS(IFREQ) 
      KC = K1+IHBMPT(IFREQ) 
      DO 710 KS= K1,K2
      K3 = KC-KS
      HPBW = BMWDTH(IFREQ)
      BEAM = EXP(-1.0*HPBW*K3*K3) 
      SFLUX = SFLUX + BEAM*SCAN(2,KS) 
      SBMSQ = SBMSQ + BEAM*BEAM 
 710  CONTINUE
      FLUX(M)  = 1000.*SFLUX/SBMSQ
  720 CONTINUE
      FLUX(3)  = (FLUX(2) - FLUX(1))/2. 
      WRITE (7,730) FLUX,RMS,FREQ(IFREQ),DEC0,PHASE,NCYCL 
 730  FORMAT (4F8.0,F7.1,2F10.3,I5) 
* 
* 
 800  IF(IFREQ.LT.3) GO TO 659
* 
      SFLUX = 0.0 
      SBMSQ = 0.0 
        DO 820  M = 37,54 
                HPBW = BMWDTH(IFREQ)
                BMPT = M - 45.5 
                BEAM = EXP(-1.0*HPBW*BMPT*BMPT) 
                SFLUX = SFLUX + BEAM*SCAN(2,M)
                SBMSQ = SBMSQ + BEAM*BEAM 
 820            CONTINUE
* 
      FLUX(1) = 1000.*SFLUX/SBMSQ 
* 
      WRITE (7,830) FLUX(1),RMS,FREQ(IFREQ) 
 830  FORMAT (16X,2F8.0,F7.1) 
* 
*     Check whether BIT 15 is set. If so then terminate program.
* 
 659  IF (ISSW(15).LT.0) GO TO 1201 
* 
* 
 660  CONTINUE
      WRITE (7,770) 
 770  FORMAT( ) 
* 
*     Print scan data 
* 
      DO 779 IDATA = 1,IEND 
      WRITE (7,783) SCAN(1,IDATA) 
 783  FORMAT(F8.3,"_")
      IREM = IDATA/10 
      REM  = IREM 
      ADATA = IDATA 
      TEST = ADATA/10 - REM 
        IF (TEST.EQ.0.00000)  WRITE (7,784) 
 784    FORMAT (/)
 779  CONTINUE  
      WRITE (7,784) 
      IPAGNO = IPAGNO + 1 
 1000 CONTINUE
 1201 WRITE (7,1300)
 1300 FORMAT (//"** END** ")
      IF (IPRAM(2).EQ.1) THEN 
          CALL EXEC (12,0,2,0,-10)
          CALL EXEC (2,113B,0,1,3,3)
          CALL EXEC (2,113B,7,1,3,3)
          CALL EXEC (12,0,2,0,-2) 
          CALL EXEC (2,113B,0,1,3,3)
        ELSE
        END IF
      END 
      SUBROUTINE PTCOR (HAIN,DECIN,DHA,DDEC)
* 
********************************************************
*                                                      *
*     ROUTINE TO COMPUTE TELESCOPE POINTING ERROS      *
*     ABSTRACTED FROM COMND.                           *
*                                                      *
*     HA ERRORS COMPUTED FROM 5TH ORDER 2D POLYNOMIAL  *
*     DEC ERRORS FROM A 3RD ORDER 2D POLYNOMIAL        *
* 
*     ADAPTED BY GD NICOLSON   1982 MAY 18             *
*     REVISED AS A SUBROUTINE  1985 NOV 18             *
* 
********************************************************
* 
* 
      DATA LU /1/, IDEC /2HDE/, IHA /2HHA/
*     Calculate HA error for range of hour angle -90 to +90 
*     in 10 degree increments 
* 
             IF (HAIN.GE.360) HAIN = HAIN - 360 
               DHA  = HAPC (HAIN,DECIN) 
               DDEC = DECPC (HAIN,DECIN)  
* 
* 
* 
      END 
* 
* 
      FUNCTION HAPC(HAIN,DECIN) 
C     ================
C 
C     POINTING CORRECTION IN DEG TO BE ADDED TO TELESCOPE RA
C 
* 
         DATA H0D0,H1D0,H2D0,H3D0,H4D0,H5D0,
     *        H0D1,H1D1,H2D1,H3D1,H4D1,H5D1,
     *        H0D2,H1D2,H2D2,H3D2,H4D2,H5D2,
     *        H0D3,H1D3,H2D3,H3D3,H4D3,H5D3,
     *        H0D4,H1D4,H2D4,H3D4,H4D4,H5D4,
     *        H0D5,H1D5,H2D5,H3D5,H4D5,H5D5/
     *          -.151191063D+02,   .699575742D+00,   .138305730D-02,
     *           .106305573D-03,  -.157377145D-07,  -.789691772D-09,
     *          -.233881793D+00,  -.129692110D-01,  -.374238417D-04,
     *           .488644625D-05,   .317452014D-08,  -.121006399D-09,
     *           .704822802D-02,   .184092354D-03,  -.208731803D-05,
     *          -.128614688D-07,  -.575064514D-10,   .776776250D-11,
     *           .182036606D-05,   .544097068D-05,   .596358955D-07,
     *          -.131267735D-08,  -.174680497D-11,   .889230070D-13,
     *          -.669984226D-05,   .234148979D-07,   .144874677D-08,
     *           .116459191D-10,   .210315949D-13,  -.349541451D-14,
     *          -.764562857D-07,  -.111917596D-08,   .584015733D-11,
     *           .243451974D-12,   .352199776D-15,  -.401941555D-16 
     *   /
C 
         HA = PMOD(HAIN)
         DEC = PMOD(DECIN)
         HAPC = (((((((((H5D5*HA+H4D5)*HA+H3D5)*HA+ 
     *             H2D5)*HA + H1D5)*HA + H0D5)*DEC +
     *             ((((H5D4*HA + H4D4)*HA + H3D4)*HA +
     *             H2D4)*HA + H1D4)*HA + H0D4)*DEC +
     *             ((((H5D3*HA + H4D3)*HA + H3D3)*HA +
     *             H2D3)*HA + H1D3)*HA + H0D3)*DEC +
     *             ((((H5D2*HA + H4D2)*HA + H3D2)*HA +
     *             H2D2)*HA + H1D2)*HA + H0D2)*DEC +
     *             ((((H5D1*HA + H4D1)*HA + H3D1)*HA +
     *             H2D1)*HA + H1D1)*HA + H0D1)*DEC +
     *             ((((H5D0*HA + H4D0)*HA + H3D0)*HA +
     *             H2D0)*HA + H1D0)*HA + H0D0 
         RETURN 
      END 
      FUNCTION DECPC(HAIN,DECIN)
C     ================= 
C 
C     POINTING CORRECTION IN DEG TO BE ADDED TO TELESCOPE DEC 
* 
* 
C 
         DATA H0D0,H1D0,H2D0,H0D1,H1D1,H2D1,H0D2,H1D2,H2D2/ 
     *           .109778425D-01,  -.481926553D+00,  -.129761141D-01,
     *           .174368342D-02,  -.196094415D-02,  -.121406856D-03,
     *          -.535791318D-02,   .268791977D-04,  -.145048381D-06 
     *   /
C 
         HA = PMOD(HAIN)
         DEC = PMOD(DECIN)
         DECPC = (((H2D2*HA + H1D2)*HA + H0D2)*DEC +
     *               ((H2D1*HA + H1D1)*HA + H0D1))*DEC +
     *                (H2D0*HA + H1D0)*HA + H0D0
         RETURN 
      END 
      FUNCTION PMOD(ANGLE)
C     ====================
C 
C     RETURNS A VALUE OF ANGLE IN PMOD, -180.0<PMOD<180.0 
C 
       PMOD = ANGLE 
         IF (PMOD .GT.90.0) PMOD = PMOD - 360.0 
         RETURN 
      END 
      SUBROUTINE OBSOK
**********************
C     subroutine to send "computer okay" alarm reset signal <990122.1505> 
c     to be called at the start of each observing cycle, at least once / hr 
c     See alarm panel for legal CPU OK alarm settings.
c 
c     ensure DIO channel 4 output is zero (+4V) 
      CALL EXEC (2,113b,0,1,4,4)
c     wait one second 
      CALL EXEC (12,0,2,0,-1) 
c     set DIO channel 4 output line 1 high  (0V)
      CALL EXEC (2,113b,1,1,4,4)
c     wait one second 
      CALL EXEC (12,0,2,0,-1) 
c     set output to zero again
      CALL EXEC (2,113b,0,1,4,4)
      RETURN
      END 
********* 
      END$
      END$
