FTN66 
$FILES 2,2
$ALIAS /DOUT1/,NOALLOCATE 
$ALIAS /POINT/,NOALLOCATE 
      PROGRAM SIXCM(3,80),  6cm Drift scan CR16 <910103.1337> 
* 
* 
******************************************************************
*                                                                *
*      SIX CENTIMETRE DRIFT CURVE OBSERVING AND ANALYSIS PROGRAM *
*                                                                *
*      ORIGINAL VERSION WRITTEN BY D CAMERON  JULY 1980          *
*                                                                *
*      REVISED BY G D NICOLSON    1983 MAY 25                    *
*                                                                *
******************************************************************
* 
      LOGICAL PSUSPD
* 
* 
      DIMENSION ABUF(400),IOBUF(20),ITIME(5),ATIME(5),ICIRX1(3),
     *          DSCAN(2),ISIZ(2),DECPNT(2),NAMDP(3),IANG(5),INAM(3),
     *          NPROG(3)
* 
      INTEGER NWAKE,CMD(12),SORSBF(40),FILNAM(3),FDCB(144), 
     *        SORSNM(4),DDCB(144),DATABF(61),FLNAM1(3),STEST1(4)
* 
      REAL COORDS(6)
* 
* 
      LOGICAL HASTOP,DECSTP,ENDFIL
* 
      EQUIVALENCE (CMD(1),RA0),(CMD(3),DEC0), 
     *            (CMD(5),IRAT),(CMD(6),STOPTM),
     *            (CMD(10),HASTOP),(CMD(11),DECSTP),
     *            (CMD(12),NWAKE),(COORDS(5),RHA),
     *            (COORDS(6),RDEC),(IANG(1),RAP),(IANG(3),DECP) 
* 
      COMMON/POINT/LOGCLS,NCCMD,J,COMCO/DOUT1/ISERVO
* 
* 
****** INITIALIZE PARAMETERS FOR STEER
* 
      DATA IRAT/4/,NAMDP/2HDP,2HRC,2H2 /, 
     *     ICIRX1 /2HCI, 2HR , 2HX1/,  INAM /2HOB, 2HCI, 2HR /, 
     *     IVELA /2HVE/,
     *     IPS   /2HPS/ 
     *     NPROG /2HPL,2HS7,1H7/
* 
      write (7,'(a)') char(14b) ! form feed 
      STOPTM=-1.
      ISIZ(1)=30
      ISIZ(2)=60
      FLAG=0.0
* 
* Scheduling routine allocates RN 
      CALL RNALC(IRN) 
* 
*     SELECT DVM FOR SIX CM FILTERED OUTPUT 
*     ------------------------------------- 
* 
      CALL EXEC (2,113B,33B,1,5,5)
* 
******CLEAR SWITCH REGISTER AND SET PLOTTER PAGE COUNT TO ZERO
* 
      CALL ISSR(0)
      IPLOT = 0 
* 
******READ TIME 
      CALL EXEC(11,ITIME,IYEAR) 
* 
* 
******WRITE HEADINGS AND OBSERVER INSTRUCTIONS
* 
      WRITE (7,5300) ITIME(5),IYEAR,ITIME(4), ITIME(3)
* 
      WRITE (7,5305)
* 
******CREATE AND OPEN DATA OUTPUT FILE
* 
  10  WRITE(7,5000) 
* 
  15  READ(7,5010)FLNAM1
* 
      CALL CREAT(DDCB,IERR,FLNAM1,ISIZ,3,0,-15) 
      IF(IERR.GE.0)GO TO 20 
      WRITE (7,5005) IERR 
      GO TO 15
* 
 20   CALL OPEN(DDCB,IERR,FLNAM1) 
      IF(IERR.LT.0)GO TO 10 
* 
      WRITE(7,5020) 
* 
      READ(7,*)GASF 
* 
50    DUMMY=0.
* 
70    WRITE(7,5030) 
* 
      READ(7,5010)FILNAM
* 
* 
******DETERMINE FIRST SOURCE TO BE OBSERVED 
* 
  71  CALL SIDTM(ST)
      ST = ST + 0.33
* 
****** READ SOURCE NO.  ZERO INDICATES SEARCH FOR SOURCE NEAREST MERIDIAN 
* 
      WRITE (7,5100)
      READ (7,*) NOSORS 
* 
****** WRITE HEADINGS 
* 
      WRITE (7,5120) GASF 
      WRITE (7,5130)
* 
******
* 
      SEARCH = 1
      IF (NOSORS.EQ.0)  SEARCH = 0
* 
****** OPEN FILE CONTAINING SOURCE LIST 
* 
      CALL OPEN(FDCB,IERR,FILNAM) 
      IF(IERR.GE.0) GO TO 145 
      WRITE(7,5140) IERR
      GO TO 70
* 
******  TEST WHETHER TO SEARCH FOR MERIDIAN SOURCE  
* 
 145  IF (SEARCH.GT.0) GO TO 155
      HAEAST = -5.0 
 150  NOSORS = 1
* 
****** READ SOURCE FROM FILE
* 
 155  CALL POSNT(FDCB,IERR,NOSORS,1)
 160  CALL READF(FDCB,IERR,SORSBF,40,LEN) 
      IF(IERR.GE.0) GO TO 190 
      WRITE(7,5150) IERR
 190  IF(LEN.GE.0) GO TO 200
           WRITE (7,5160) 
           IF(SEARCH.EQ.1) GO TO 150
                HAEAST = HAEAST - 5.0 
* 
                WRITE (7,5400) HAEAST 
 5400           FORMAT(/"MERIDIAN WINDOW SET TO ",F6.1,"   TO +5.0 DEG")
* 
                GO TO 150 
* 
 200  CALL CODE 
      READ(SORSBF,5170)SORSNM,RA,DEC,ON 
      IWAVE = ON
      IF(SEARCH.EQ.1) GO TO 210 
      HATEST = ST*15 - RA 
      IF (HATEST.GE.5..OR.HATEST.LE.HAEAST) GO TO 160 
      SEARCH = 1
* 
*     Test whether next source is Circinus X-1, Vela pulsar or P1641
* 
  210    IF(SORSNM(1).NE.ICIRX1(1)) GO TO 211 
* 
*           Schedule program OBCIR to make 3 scans of Circinus X-1
* 
             KFREQ = ON 
             CALL EXEC (9,INAM,KFREQ) 
*            do a form feed on printer after OBCIR finishes 
             write (7,'(a)') char(14b)
* 
*            Write headings 
*            ============== 
* 
*            Read time and year 
* 
            CALL EXEC (11,ITIME,IYEAR)
* 
             WRITE (7,7401) ITIME(5),IYEAR,ITIME(4),ITIME(3)
             WRITE (7,5130) 
             GO TO 213
* 
*            Test whether source is VELA PULSAR or P1641-45 
* 
  211      IF(SORSNM(1).NE.IVELA) GO TO 212 
* 
*            Schedule pulsar program
* 
             CALL EXEC ( 9,NPROG,3,IWAVE) 
* 
             GO TO 213
  212          IF(SORSNM(1).NE.IPS)  GO TO 215
* 
*                  Schedule pulsar program - option 4 (1641-45) 
* 
                   CALL EXEC ( 9,NPROG,4,IWAVE) 
* 
                   GO TO 213
* 
* 
* 
*     Re-initialise plotter page number 
* 
  213 IPLOT = 0 
* 
*     Advance printer to top of form
* 
      WRITE (7,'(a)') char(14b) 
      WRITE (7,5130)
* 
* 
****** Initialise parameters
* 
*     Read next source to be observed 
* 
      GO TO 160 
* 
 215  IPLOT = IPLOT + 1 
      IF (IPLOT.EQ.7)  IPLOT = 1
      SPMA = 0.0
      SPMB = 0.0
* 
****** COMPUTE RIGHT ASCENSION AND DECLINATION FOR STEER
* 
      RAP = RA
      DECP = DEC
* 
****** PRECESS 1950 CORDINATES TO CURRENT EPOCH 
* 
      CALL EXEC (9,NAMDP,IANG(1),IANG(2),IANG(3),IANG(4),IANG(5)) 
      CALL RMPAR (IANG) 
* 
*     Add beam offsets to right ascension and allow for filter delay
* 
      RAP = RAP + 0.070/(COS(DEC/57.3)) 
      RA0=RAP 
      DEC0=DECP 
      RA0=RA0-(.5+(.5/(COS(DEC0/57.3))))
* 
****** CORRECT FOR DECLINATION BEAM OFFSET
* 
      DEC0=DEC0-.203
* 
* 
      DECP0 = DEC0
      IF (DECP0.GE.270.) DECP0 = DECP0 -360.
* 
****** TEST WHETHER DECLINATION POINTING SCANS ARE REQUIRED 
*      AND APPLY HALF BEAMWIDTH CORRECTIONS 
* 
      IF(ON.NE.10) GO TO 220
      DEC0=DEC0-.08 
220   CONTINUE
* 
****** COMPUTE SCAN DURATION IN SECONDS 
* 
      TM = (240. + (240./(COS(DEC0/57.3)))) 
      ITM = TM
* 
******DRIVE TO FIRST SCAN 
* 
      DO 230 I=1,2
* 
*        Turn off noise tube at start of scan 
* 
         CALL EXEC (2,113B,0,1,3,3) 
* 
* 
*     SELECT DVM
* 
      CALL EXEC (2,113B,33B,1,5,5)
      HASTOP=.TRUE. 
      DECSTP=.TRUE. 
      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)
******Drive again after waiting 1 second
* 
      CALL EXEC (12,0,2,0,-1) 
      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)
      CALL EXEC(12,0,2,0,-10) 
* 
******END OF DRIVE
* 
****** WRITE SOURCE NAME
* 
       WRITE (7,5410) SORSNM
 5410  FORMAT (4A2,"_") 
* 
* 
******SAMPLE DATA FOR ITM SECONDS 
* 
      CALL SIDTM(ST)
      CALL DVMR(I,SORSNM,ITM,IOBUF,ABUF,RDEC,RHA,AVG1)
* 
******CALCULATE CENTRE OF SCAN
* 
      IF(RHA.LT.0.0)RHA=RHA+360.0 
      IF(RDEC.LT.0.0)RDEC=RDEC+360.0
      DECPNT(I) = RDEC
      IF (DECPNT(I).GE.270.) DECPNT(I) = DECPNT(I) - 360. 
      CENT1=(ST*15)-RHA 
      IF(CENT1.LT.0.0)CENT1=CENT1+360.
      CALL HAPC(RHA,RDEC,HAP) 
      RAPAP=RAP+HAP 
      CENT=(RAPAP-CENT1)
      IF(CENT.LT.0.0)CENT=CENT+360. 
      CENT=((CENT/15)*3600.)/5. +  1. 
      DHA = CENT
      CENT=ABS(CENT)
* 
****** ANALYZE SCAN BY USING GAUSSIAN BEAM FITTING PROCEDURE
* 
      CALL BEAMF(ITM,ABUF,RDEC,SPA,SPB,SIGMA,CENT,I,IPLOT,DHA 
     *,SORSNM,FILNAM) 
* 
******FIRE NOISE TUBE BETWEEN SCANS 
* 
      NTRY = -5 
      IF(I-1)240,240,250
 240  CALL EXEC(2,113B,0,1,3,3) 
      CALL EXEC (12,0,2,0,-2) 
      CALL EXEC(2,113B,3,1,3,3) 
      CALL NOISE(40,AVG,-2) 
      NTON = AVG - AVG1 - 250.
      NTRY = NTRY + 1 
      IF (NTRY.EQ.0) GO TO 250
      IF(NTON.LT.0) GO TO 240 
      CALL NOISE(400,AVG,-2)  
      AVG2=AVG
      AVG4=AVG2 - AVG1
      CALL EXEC (2,113B,0,1,3,3)
* 
******SCAN ANALYSIS 
* 
  250 FLUXF=GASF/AVG4 
      SIGMA=SIGMA*FLUXF 
      FLUXA=ABS(SPA)*FLUXF
      FLUXB=ABS(SPB)*FLUXF
      FLUX=(FLUXA+FLUXB)/2
      WRITE(7,5180) FLUXA,FLUXB,FLUX,SIGMA,DHA,RHA,RDEC 
      DSCAN(I) = ABS(SPA) + ABS(SPB)
      SPMA=SPA+SPMA 
      SPMB=SPB+SPMB 
      IF(ON.NE.10.)GO TO 230
      DEC0=DEC0+.16 
230   DUMMY=0.
      SPAM=(ABS(SPMA/2.)/AVG4)*GASF 
      SPBM=(ABS(SPMB/2.)/AVG4)*GASF 
 260  DUMMY=0.
* 
******WRITE OUT FLUX VALUES USING ASSUMED VALUE FOR NOISE TUBE
* 
 270  IF(ON.NE.10) GO TO 280
 290  DSTH = ALOG(DSCAN(1)) 
      DNTH = ALOG(DSCAN(2)) 
      PCONST =0.0289
      OFSET = PCONST*(DSTH-DNTH)
      WRITE(7,5260)SORSNM,
     *OFSET,RHA,RDEC
      GO TO 300 
 280  FLUX=(SPAM+SPBM)/2. 
      WRITE(7,5270)SORSNM,SPAM, 
     *SPBM,FLUX, AVG4 
 300  IF(ON.NE.10) GO TO 310
      CALL CODE 
      WRITE(DATABF,5190)SORSNM, 
     *SPAM,SPBM,OFSET,SIGMA,RHA,RDEC,ON 
      GO TO 320 
 310  CALL CODE 
      WRITE(DATABF,5190)SORSNM(1),SORSNM(2),SORSNM(3),SORSNM(4),
     *SPAM,SPBM,FLUX,SIGMA,RHA,RDEC,ON
320   WRITE (7,5210)
      CALL WRITF(DDCB,IERR,DATABF,61) 
      IF(IERR.LT.0) GO TO 10
      IF (ISSW(15).LT.0) GO TO 180
* 
*     Suspend program unconditionally 
* 
        IF (ISSW(13).GE.0) GO TO 7500 
          WRITE (7,7400)
 7400     FORMAT ("Program suspended indefinitely. To restart type" 
     *            " GO,SIXCM")
* 
*         Clear switch register 
          CALL ISSR(0)
          pause 
          write (7,'(a)') char(14b)    ! form feed after restart
* 
*     Re-initialise page number 
* 
      IPLOT = 0 
* 
*     Write headers 
* 
*     Read time 
* 
      CALL EXEC (11,ITIME,IYEAR)
* 
* 
      WRITE (7,7401) ITIME(5),IYEAR,ITIME(4),ITIME(3) 
 7401 FORMAT(/////"SIX CENTIMETRE OBSERVING PROGRAM DAY",I4,
     1       "  YEAR ",I4,I4,"h",I2,//) 
      WRITE (7,5130)
7500        CONTINUE
* 
*     Option to go to new source if bit 14 is set 
* 
        IF (ISSW(14).LT.0)  THEN
           CALL ISSR(0) 
           GO TO 71 
        END IF
* 
* Let SCHDL interrupt if it wants 
        CALL RNCLW(IRN,PSUSPD)
        IF (.NOT.PSUSPD)  GOTO 178
* 
*     Re-initialise page number 
* 
          IPLOT = 0 
* 
*     Read time 
* 
          CALL EXEC (11,ITIME,IYEAR)
* 
*     Write headers 
          WRITE (7,7401) ITIME(5),IYEAR,ITIME(4),ITIME(3) 
          WRITE (7,5130)
178     CONTINUE
* 
* 
* End of cycle
        GOTO 160
* 
 180  CALL CLOSE(DDCB,IERR) 
      WRITE(7,5220) 
      WRITE(7,5240) 
      CALL OPEN(DDCB,IERR,FLNAM1) 
 330  CALL READF(DDCB,IERR,DATABF,61) 
      IF(IERR.LT.0)GO TO 340
      CALL CODE 
      READ(DATABF,5190)SORSNM,FLUXA,FLUXB,FLUX,SIGMA,RHA,RDEC,ON
      IF(ON.LE.10) GO TO 330
      CALFA=ON/FLUX 
      CALL CLOSE(DDCB,IERR) 
      CALL OPEN(DDCB,IERR,FLNAM1) 
 360  CALL READF(DDCB,IERR,DATABF,61,LEN) 
      IF(LEN.LT.0)GO TO 350 
      CALL CODE 
      READ(DATABF,5190)SORSNM,SPAM,SPBM,FLUX,SIGMA,RHA,RDEC,ON
      IF(ON.EQ.10) GO TO 360
      SPAM=SPAM*CALFA 
      SPBM=SPBM*CALFA 
      FLUX=FLUX*CALFA 
 370  WRITE(7,5190)SORSNM,SPAM,SPBM,FLUX,SIGMA,RHA,RDEC,ON
      GO TO 360 
 340  WRITE(7,5250) 
 350  CALL CLOSE(DDCB,IERR) 
* 
*     Clear switch register 
* 
      CALL ISSR(0)
* 
         WRITE (7,5430) 
5430     FORMAT(//"*END*"//////)
* 
5000  FORMAT(///"ENTER NAME OF OUTPUT DATA FILE"/ 
     1 "USE SA FOLLOWED BY DAY NO. EG SA076"//) 
5005  FORMAT("FILE NAME ALREADY IN USE.  FMGR ERROR  ",I3,/ 
     1 "USE SB, SC, SD ETC FOLLOWED BY DAY NO.  EG SB076"//)
5010  FORMAT(3A2) 
5020  FORMAT("ENTER NOISE TUBE FLUX IN JANSKYS _")
5030  FORMAT("ENTER SOURCE FILE NAME _")
5050  FORMAT("TYPE 1 IF NEW MAGTAPE ON, OTHERWISE TYPE 0 _")
5060  FORMAT(3A2,2I6) 
5100  FORMAT("ENTER NO OF FIRST SOURCE IN FILE TO BE OBSERVED"/ 
     1 "OR FOR PROGRAM TO LOCATE SOURCE NEAREST THE MERIDIAN TYPE 0"
     2 //)
5300  FORMAT (///////"SIX CENTIMETRE OBSERVING PROGRAM DAY",I4, 
     1 "  YEAR ",I4,I4,"H",I2,/,60"*"/
     2 /"* THIS VERSION USES FLUKE 8840 DVM2" / 
     3 "*"/"*  Check TIME by typing *TI on system console VDU "/
     4 "  DVM1  must be in COMPUTER mode!"/"*"/"*") 
5305  FORMAT ("*  At end of run set SWITCH REGISTER bit 15 to ter"
     1 "minate"/
     2 "*  Bit 14 allows new source to be selected "/ 
     3 "*  Bit 13 suspends the program unconditionally."/ 
     4 "*  Use *GO,SIXCM to restart."/"*"/60"*")
5140  FORMAT("FILE NOT FOUND-FMP ERR")
5150  FORMAT("FILE ERROR-PROGRAM MUST BE TERMINATED") 
5160  FORMAT(/"END OF SOURCE FILE-SCAN FROM START") 
5170  FORMAT(4A2,F9.3,F9.3,F5.1)
5190  FORMAT(4A2,6F8.3,F5.1)
5210  FORMAT(/) 
5220  FORMAT("                   FINAL DATA ANALYSIS PRINT OUT")
5240  FORMAT(" SORSNM     FLUXA   FLUXB   FLUX    SIGMA   RHA    RDEC 
     *ON")
5250  FORMAT("NO CALIBRATOR OBSERVED-NO PRINT OUT") 
5120  FORMAT(20X,"PRELIMINARY FLUX FOR GAS FACTOR = ",F7.3) 
5130  FORMAT(6HSOURCE,2X,6HFLUX A,2X,6HFLUX B,2X,4HFLUX,2X, 
     *5HSIGMA,2X,10HNOISE CAL ,2X,6HOFFSET,3X,2HRA,4X,3HDEC)  
5260  FORMAT(4A2,40X,F6.3,2F8.3)
5270  FORMAT(4A2,F6.2,2X,F6.2,F6.2,8X,F5.0) 
5180  FORMAT(F6.2,2X,F6.2,F6.2,X,F6.3,14X,F5.3,2F8.3) 
      END 
      SUBROUTINE SIDTM(ST)
**************************
      DIMENSION ITIME(5),ATIME(5) 
      ICODE=11
      CALL EXEC(ICODE,ITIME,IYEAR)
      DO 10 J3=1,5
      ATIME(J3)=ITIME(J3) 
 10   DUMMY=0 
      UTIME=(ATIME(4)+ATIME(3)/60.+ATIME(2)/3600.+ATIME(1)/360000.) 
      IF (IYEAR.EQ.1990) ATIME(5) = ATIME(5) + 365
      ST=(UTIME*9.856)/3600.+UTIME+8.456556 + ATIME(5)*0.0657097  
  30  IF(ST.LT.24.) GO TO 40
      ST = ST - 24. 
      GO TO 30
  40  DUMMY=0.0 
      RETURN
      END 
$ALIAS /DOUT1/,NOALLOCATE 
$ALIAS /POINT/,NOALLOCATE 
      SUBROUTINE DVMR(I,SORSNM,ITM,IOBUF,ABUF,RDEC,RHA,AVG1)
************************************************************
      COMMON/DOUT1/ISERVO 
      DIMENSION  IOBUF(20),ABUF(400), SAMPLE(100) 
      INTEGER SORSNM(4) 
      ILENG=ITM/5+1 
      ICON=128
      MID = ILENG/2 
      DO 50  J1=1,ILENG 
*     
*     Sample once per second for five seconds and take median 
* 
             WRITE (35,'("*S1")') 
      CALL EXEC (2,113B,33B,1,5,5)
          DO 40 J2 = 1,98 
             READ (35,*) SAMPLE(J2) 
   40     CONTINUE
              IF (J1.NE.MID) GO TO 35 
* 
*     Read Hour Angle and Declination 
*   
              CALL EXEC (1,11,IOBUF(1),10,1,3)
* 
*             Mask off tachometer bits
* 
              IOBUF(1) = IOBUF(1).AND.171777B 
              IOBUF(4) = IOBUF(4).AND.171777B 
* 
                 CALL CODE  
                 READ (IOBUF,33) RHA,RDEC 
 33              FORMAT (2F6.3) 
* 
* 
* 
 35    CONTINUE 
* 
* 
*             Find median of 98 samples 
* 
                 NS = 98
                 MS = NS
 41              MS = MS/2
                 IF(MS.EQ.0) GO TO 45 
                    KS = NS - MS  
                    JS = 1
 42                 IS = JS 
 43                 LS = IS + MS
                       IF (SAMPLE(IS).LE.SAMPLE(LS)) GO TO 44 
                          BS = SAMPLE(IS) 
                          SAMPLE(IS) = SAMPLE(LS) 
                          SAMPLE(LS) = BS 
                          IS = IS - MS
                             IF(IS.GE.1) GO TO 43 
 44                             JS = JS + 1 
                                IF(JS.LE.KS) GO TO 42 
                                   GO TO 41 
* 
*         Sort complete 
* 
 45    ABUF(J1) = -SAMPLE(25)*100000
* 
 50   CONTINUE
* 
*      Average last 12 points of baseline for noise cal reference 
* 
      AVG1 = 0.0
      DO 60 J1 = ILENG-11, ILENG
      AVG1 = AVG1 + ABUF(J1)
 60   CONTINUE
      AVG1 = AVG1/12
      RETURN
      END 
      SUBROUTINE BEAMF(ITM,ABUF,RDEC,SPA,SPB,SIGMA,CENT,I,IPLOT,DHA 
     *                 ,SORSNM,FILNAM)
********************************************************************
      DIMENSION ABUF(400),ANEW(400),ITIME(5)
      INTEGER SORSNM(4),FILNAM(3) 
******DIFFERENT SLOPES TO BASE LINE 
   5  ILENG=ITM/5+1 
      ISP=1 
      ILP=24
      IC=ILENG
      N=48
      U=COS(RDEC/57.3)
******JMP TO LEAST SQ FIT FOR BASE LINE 
      CALL LSQB(ISP,ILP,ABUF,A,B,IC,N)
******SUBTRACT OFF SLOPE AND BIAS FROM CURVE
      A8=A
      B8=B
      DO 10 J1=1,ILENG
      AJ1=J1
      ANEW(J1)=ABUF(J1)-(A*AJ1+B) 
  10  CONTINUE
******DERIVE SIGMA VALUE FOR BASE LINE
      SIG=0.
      DO 20 J=1,24
      AJ1=J 
      SIG=(ANEW(J)*ANEW(J))+SIG 
  20  CONTINUE
      ILEN=ILENG-23 
      DO 30 J=ILEN,ILENG
      AJ1=J 
      SIG=(ANEW(J)*ANEW(J))+SIG 
  30  CONTINUE
      SIGMA=SQRT((SIG/48.)) 
******TEST FOR POINTS GREATER THAN SIGMA
      NREJ =0 
      DO 40 J=1,24
      IF(ABS(ANEW(J))-3.*SIGMA)40,40,50 
  50  ABUF(J)=A*J+B 
      NREJ = NREJ+1 
  40  CONTINUE
      ILEN=ILENG-23 
      DO  60 J=ILEN,ILENG 
      IF(ABS(ANEW(J))-3.*SIGMA)60,60,70 
  70  ABUF(J)=A*J+B 
      NREJ = NREJ+1 
  60  CONTINUE
      IF(NREJ.GT.0) GO TO 5 
******FIND CENTRE OF SCAN (ILENG/2) AND USE 5 PTS 
******ON EACH SIDE FOR LSQ FIT TO DTERMINE ZERO CROSSING
      ILE=CENT
      ISP=(ILE)-4./U
      ILP=(ILE)+4./U
      IC=-1 
      N = 1 + ILP - ISP 
      CALL LSQB(ISP,ILP,ANEW,A,B,IC,N)
******DETERMINE ZERO CROSSING 
      ZE1=CENT
      ZE=-B/A 
      DHA = (DHA - ZE)*0.0208 
******FIND CENTRE OF BOTH BEAMS 
      ABEAM=ZE-6.0000*(1./COS(RDEC/57.3)) 
      BBEAM=ZE+6.0000*(1./COS(RDEC/57.3)) 
******FIND ROUGH ESTIMATE FOR CENTRE OF BEAMS 
      DO 80  J=1,ILENG
      AJ=J
      IF(ABS(AJ-ABEAM)-.6)90,90,100 
  90  ABEM=J
 100  IF(ABS(AJ-BBEAM)-.6)110,110,80
 110  BBEM=J
  80  CONTINUE
******FIT CURVE TO A AND B BEAMS
      SYB=0.
      SB=0. 
      IABE=ABEM-5./U
      IABM=ABEM+5./U
      DO 120 J=IABE,IABM
      AJ=J
      YB=((AJ-ABEAM)/(3.84/U))
      SYB=SYB+ANEW(J)*EXP(-0.693*(YB*YB)) 
      SB=SB+(EXP(-0.693*(YB*YB)))*(EXP(-0.693*(YB*YB))) 
 120  CONTINUE
      SPA=SYB/SB
******B BEAM
      SYB=0.
      SB=0. 
      IABE=BBEM-5./U
      IABM=BBEM+5./U
      DO 130 J=IABE,IABM
      AJ=J
      BY=((AJ-BBEAM)/(3.84/U))
      SYB=SYB+ANEW(J)*EXP(-.693*(BY*BY))
      SB=SB+(EXP(-.693*(BY*BY)))*(EXP(-.693*(BY*BY))) 
 130  CONTINUE
      SPB=SYB/SB
******ANALYTICAL CURVE
      IABEM=ABEM-5./U 
      IABE=ABEM+5./U
      DO 220 J=IABEM,IABE 
      AJ=J
      AY=((AJ-ABEAM)/(3.84/U))
      ANEW(J)=EXP(-.693*(AY*AY))*SPA
 220  CONTINUE
      IABEM=BBEM-5./U 
      IABE=BBEM+5./U
      DO 140 J=IABEM,IABE 
      AJ=J
      YA=((AJ-BBEAM)/(3.84/U))
      ANEW(J)=+EXP(-.693*(YA*YA))*SPB 
 140  CONTINUE
******BEGIN TO PLOT 
      IF (IPLOT.GT.1) GO TO 145 
      IF (I.EQ.2) GO TO 145 
C     DO 143 IPAGE = 1,2
*     Write page marker at top right corner 
* 
      MX = 9990 
      MY = 9990 
* 
      WRITE (12) -1,1,MX,MY 
* 
          MY = 9875 
* 
          WRITE (12)  1,1,MX,MY 
          CALL EXEC (12,0,2,0,-2) 
          WRITE (12) -1,1,MX,MY 
      CALL EXEC(2,113B,7,1,3,3) 
      CALL EXEC(12,0,2,0,-1)
      CALL EXEC(2,113B,0,1,3,3) 
      CALL EXEC(12,0,2,0,-5)
 143  CONTINUE
      IXX =  0
      IXY =  200
      IYX = -160
      IYY =  0
      WRITE (12) -1,1,250,500 
      CALL EXEC (11,ITIME,IYEAR)
      WRITE (12,6000) IXX,IXY,IYX,IYY,ITIME(5),IYEAR,(FILNAM(M),M=1,3)
 6000 FORMAT ( 4I5,5HDAY  ,I3,8H  YEAR  ,I4, 15H  SOURCE FILE  ,3A2)
 145  X1=ILENG
      M1=-1 
      KX=(2*IPLOT+0.5)*666.7
      KY=(I-1)*5000 
      WRITE(12)M1,1,KX,KY 
      DO 150 K=1,2
      M1=1
      KY1=5000*I-1000 
      WRITE(12)M1,1,KX,KY1
      WRITE(12)-1,1,KX,KY 
      DO 160 J=1,ILENG
      AJ=J
      IF(K-1)170,180,170
 170  Y=ABUF(J)-(A8*AJ+B8)
      GO TO 190 
 180  Y=ANEW(J) 
 190  SCALE=(((ABS(SPA)+ABS(SPB))/2)*1.1) 
      Y=(((Y+SCALE)/(SCALE))+(IPLOT) - 0.75)*2000*.6667 
      X=(((AJ-1)/(X1-1))+((I-1)*1.25))*4000 
      IX=IFIX(Y)
      IY=IFIX(X)
      WRITE(12)1,1,IX,IY
 160  CONTINUE
      WRITE(12)-1,1,IX,IY 
      IF(K.EQ.2) GO TO  200 
      IF(I.EQ.2) GO TO 200
      IXX=0 
      IXY=100 
      IYX=-160
      IYY=0 
      KY2=KY1+100 
      WRITE(12)-1,1,KX,KY2
      WRITE(12,5000)IXX,IXY,IYX,IYY,(SORSNM(M),M=1,4) 
 200  DUMMY=0.0 
 150  CONTINUE
5000  FORMAT(4I5,4A2) 
      END 
      SUBROUTINE LSQB(ISP,ILP,ABUF,A,B,IC,N)
********************************************
      DIMENSION ABUF(400) 
      SSX=0.
      SX=0. 
      SXY=0.
      CI=0. 
      SY=0. 
      DO 10  K=1,2
      IF(CI)20,30,30
  30  DUMMY=0.
      DO 40  J1=ISP,ILP 
      AJ=J1 
      SSX=SSX+(AJ*AJ) 
      SX=SX+AJ
      SXY=SXY+AJ*ABUF(J1) 
      SY=SY+ABUF(J1)
  40  CONTINUE
  20  DUMMY=0.
      IF(IC)50,60,60
  50  CI=-1.
  60  ISP=IC-23 
      ILP=IC
  10  CONTINUE
******DERIVE LINEAR EQS FOR BEST FIT
      AN=N
      A=(SXY-((SX*SY)/AN))/(SSX-AN*((SX/AN)*(SX/AN))) 
      B=SY/AN-(A*(SX/AN)) 
      RETURN
      END 
$ALIAS /DOUT1/,NOALLOCATE 
$ALIAS /POINT/,NOALLOCATE 
      SUBROUTINE NOISE(NOPTS,AVG,DELAY) 
*************************************** 
******NOPTS=NO OF PTS TO BE SAMPLED 
C     AVG=AVG OF SAMPLES READ BY DVM
      COMMON/DOUT1/ISERVO 
      DIMENSION IOBUF(20) 
      DVM1=0
      CALL EXEC (2,113B,33B,1,5,5)
         WRITE (35,'("*S1")') 
      DO 10 I6=1,NOPTS
         READ  (35,*) RDVM
  10  DVM1 = RDVM+DVM1
      AVG = -(DVM1/NOPTS)*100000
5000  FORMAT(I1,I5) 
      RETURN
      END 
      SUBROUTINE HAPC(RHA,RDEC,HAP) 
*********************************** 
         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/
     *        -3.577006066E-01, 3.149323906E-01, 1.773597917E-03, 
     *         1.213048605E-04, 9.195761737E-08, 5.097878259E-09, 
     *        -1.196142822E-02,-3.325711682E-02, 1.367536276E-04, 
     *         2.956253983E-06,-4.676537927E-09, 8.220894774E-10, 
     *         6.843939369E-04, 6.105595596E-04,-9.183186131E-06, 
     *        -2.313902098E-07, 2.997697834E-10, 7.804750269E-11, 
     *         3.210506975E-05, 2.031257628E-05,-2.219498734E-07, 
     *        -1.971238918E-09, 1.246548062E-11, 3.854811155E-13, 
     *        -7.075766653E-07,-1.935502987E-07, 5.511228767E-09, 
     *         1.689931705E-10,-1.550025600E-13,-5.156162593E-14, 
     *        -1.690506921E-08,-5.670167102E-09, 9.960097836E-11, 
     *         2.191014772E-12,-3.998645622E-15,-6.508315344E-16/ 
******START OF FUNCTION 
         HA = RHA 
         DEC = RDEC 
         IF(HA.GT.270)HA=HA-360 
         IF(DEC.GT.270)DEC=DEC-360
         HAP = (((((((((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 
         HAP = HAP/1000.0 
         RETURN 
      END 
      INCLUDE &SSUB6::16
      END$
