FTN66 

$ALIAS /POINT/,NOALLOCATE 
$ALIAS /DOUT1/,NOALLOCATE 
      PROGRAM THRCM(3,80),  13cm Drift scan CR16 <910707.1552>
************************************************************************
*                                                                      *
*     Thirteen centimetre drift curve observing and analysis program   *
*     Adapted from 6 cm program "SIXCM"            1980 AUG 12         *
*     Revised .................................... 1983 MAY 05         *
*                                                                      *
************************************************************************
* 
      DIMENSION ABUF(400),IOBUF(20),ITIME(5),ATIME(5),ICIRX1(3),
     *          INAM(3), NPROG(3) 
      DIMENSION DSCAN(2),ISIZ(2),DECPNT(2),NAMDP(3),IANG(5) 
* 
      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
* 
******Initialise program parameters 
* 
      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,'('''')') !form feed
      STOPTM=-1.
      ISIZ(1)=30
      ISIZ(2)=60
      FLAG=0.0
      IPLOT=0 
* 
*     Select 13 cm filtered DVM input 
*     ------------------------------- 
* 
      CALL EXEC (2,113B,55B,1,5,5)
* 
*     Clear switch register 
* 
      CALL ISSR(0)
******READ TIME 
*     Read time 
* 
      CALL EXEC(11,ITIME,IYEAR) 
      YEAR=IYEAR-1950 
      A=ITIME(5)
      U=ITIME(4)
*     Write headings and operator instructions on Silent 700 terminal 
* 
      WRITE (7,5300) ITIME(5),IYEAR 
      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)
      GO TO 15
 20   CALL OPEN(DDCB,IERR,FLNAM1) 
      IF(IERR.LT.0)GO TO 10 
      WRITE(7,5020) 
      READ(7,*)GASF 
      ST2=6.5883
******Calculate sidereal time 
* 
      STO=ST2+A*0.0657
      IF(STO-24.)30,40,40 
40    STO=STO-24. 
30    ST1=STO+0.0657
      IF(ST1-24.)50,60,60 
60    ST1=ST1-24. 
50    DUMMY=0.
  70  WRITE(7,5030) 
      READ(7,5010)FILNAM
******POSITION MAG TAPE 
C     WRITE(7,5050) 
C     READ(7,*)NEWT 
C     IF(NEWT.EQ.1)GO TO 80 
C 90  CALL PTAPE(8,1,0) 
C100  CALL PTAPE(8,3,0) 
C     CALL PTAPE(8,0,1) 
C     IF(IEOF(8).GT.0) GO TO 100
C     CALL PTAPE(8,0,1) 
C     IF(IEOF(8).GT.0) GO TO 90 
C     CALL PTAPE(8,-1,0)
C 80  WRITE(8,5060)FILNAM,IYEAR,ITIME(5)
C     ENDFILE 8 
* 
*     Determine first source to be observed 
* 
 110  CALL SIDTM(UT,ST,UTIME,STO,ST1) 
      WRITE(7,5100) 
      READ(7,*)NUR
      POWERF = 1
      IF (NUR.EQ.0) POWERF = 0. 
      WRITE(7,5120)GASF 
      WRITE(7,5130) 
      GO TO 140 
 150  CALL POSNT(FDCB,IERR,1,1) 
      GO TO 160 
 140  CALL OPEN(FDCB,IERR,FILNAM) 
      IF(IERR.GE.0) GO TO 170 
      WRITE(7,5140)IERR 
 170  IF(IERR.LE.0) GO TO 70
* 
*     Read source cordinates from file
* 
      IF(POWERF.LE.0) GO TO 160 
      CALL POSNT(FDCB,IERR,NUR,1) 
      FLAG=1. 
 160  CALL READF(FDCB,IERR,SORSBF,40,LEN) 
      IF(ISSW(15).LT.0)GO TO 180
*   
      IPLOT=IPLOT+1 
      IF(IPLOT.EQ.8) IPLOT=1
      SPMA=0.0
      SPMB=0.0
      IF(IERR.GE.0) GO TO 190 
      WRITE(7,5150)IERR 
 190  IF(LEN.GE.0) GO TO 200
      IF(FLAG.GT.0) GO TO 180 
      ST=ST+0.33
      GO TO 140 
 195  WRITE(7,5160) 
 200  IF(LEN.LT.0) GO TO 150
      CALL CODE 
      READ(SORSBF,5170)SORSNM,RA,DEC,ON 
* 
*     Test if VELA is to be observed
* 
      IF (SORSNM(1). NE. IVELA) GO TO 211   
* 
      CALL EXEC (9,NPROG,3,181) 
* 
      IPLOT = 0 
      WRITE (7,'(a)') char(14b) 
      WRITE (7,5130)
      GO TO 160 
 211  IF(FLAG.GE.1)GO TO 210
      IF(ABS(RA-(ST*15.)).GE.5)GO TO 160
      FLAG=1. 
 210  RAP = RA
      DECP = DEC
      CALL EXEC (9,NAMDP,IANG(1),IANG(2),IANG(3),IANG(4),IANG(5)) 
      CALL RMPAR (IANG) 
* 
*     Add beam offsets and allow for filter delay 
* 
      RAP = RAP + 0.052/(COS(DEC/57.3)) + 0.010 
      RA0=RAP 
      DEC0=DECP 
      RA0=RA0-(0.5+(.5/(COS(DEC0/57.3)))) 
      TM=(240+(240/(COS(DEC0/57.3)))) 
      ITM=IFIX(TM)
      DEC0=DEC0 + 0.445 
      DECP0 = DEC0
      IF (DECP0.GE.270.) DECP0 = DECP0 -360.
      IF(ON.NE.10) GO TO 220
      DEC0=DEC0-.166
220   CONTINUE
* 
*     Drive to first scan 
* 
      DO 230 I=1,2
* 
*     Select DVM
* 
* 
      CALL EXEC (2,113B,55B,1,5,5)
* 
C      WRITE (7,217) RA0
C 217    FORMAT (55X,F7.3,"_")
      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 aftet waiting 10 seconds
* 
      CALL EXEC(12,0,2,2,-10) 
C       WRITE (7,219) RA0 
C219    FORMAT (1X,F7.3)
      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,-20) 
* 
* 
*     Telescope now at corect position
* 
*     Write source name 
* 
      WRITE (7,5410) SORSNM 
5410  FORMAT (4A2,"_")
* 
*     Sample data for ITM seconds 
* 
      CALL SIDTM(UT,ST,UTIME,STO,ST1) 
      CALL DVMR(I,SORSNM,ITM,IOBUF,ABUF,RDEC,RHA,AVG1)
*     Compute 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)
      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,4,1,3,3) 
      CALL NOISE(40,AVG,-10)
      NTON =ABS( AVG - AVG1)- 200.
      NTRY = NTRY + 1 
      IF (NTRY.EQ.0) GO TO 245
      IF(NTON.LT.0) GO TO 240 
      CALL NOISE(200,AVG,-5)
      AVG2=AVG
      AVG4=ABS(AVG2 - AVG1) 
 245  CALL EXEC (2,113B,0,1,3,3)
******SCAN ANALYSIS 
  250 FLUXF=ABS(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+.332
230   DUMMY=0.
      SPAM=(ABS(SPMA/2.)/AVG4)*GASF 
      SPBM=(ABS(SPMB/2.)/AVG4)*GASF 
 260  DUMMY=0.
******WRITE OUT FLUX VALUES USING PREVIOUS DAYS GASF
 270  IF(ON.NE.10) GO TO 280
 290  DSTH = ALOG(DSCAN(1)) 
      DNTH = ALOG(DSCAN(2)) 
      PCONST = 0.0599 
      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
      GO TO 160 
 180  CALL CLOSE(DDCB,IERR) 
C     ENDFILE 8 
C     ENDFILE 8 
      CALL EXEC (11,ITIME,IYEAR)
      WRITE (7,8000) ITIME(4),ITIME(3)
 8000 FORMAT (///"PROGRAM TERMINATED AT  ",I2,"H",I3,"M UT")
      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) 
         WRITE (7,5430) 
5430     FORMAT(//"*END*"//////)
5000  FORMAT(//"ENTER NAME OF OUTPUT DATA FILE"/
     1  "USE TA FOLLOWED BY DAY NO EG TA076 FOR DAY 076"//) 
5005  FORMAT("FILE NAME ALREADY IN USE."
     1/"USE TB, TC, TD ETC FOLLOWED BY DAY NO.  EG TB076"//)
5010  FORMAT(3A2) 
5020  FORMAT("ENTER NOISE TUBE FLUX IN JANSKYS  _") 
5030  FORMAT("ENTER SOURCE FILE NAME  _") 
5050  FORMAT("INPUT 1 IF NEW MAGTAPE ON,OTHERWISE TYPE 0  _") 
5060  FORMAT(3A2,2I6) 
5100  FORMAT("ENTER NO OF FIRST SOURCE IN SCHEDULE TO BE OBSERVED"/ 
     1 " OR FOR PROGRAM TO LOCATE SOURCE NEAREST THE MERIDIAN TYPE 0" 
     2 //)
5300  FORMAT(//////////"13 CENTIMETRE OBSERVING PROGRAM DAY",I4,
     1      "  YEAR ",I4,I2,"H",I2," UT",/,64"*"/,/ 
     2"*"/   "*CHECK TIME BY TYPING TI ON SYSTEM CONSOLE (VDU)"/
     3"*SET DVM1 TO 0.1 V RANGE, 0.1 SEC SAMPLE PERIOD",/ 
     4  "*AND SAMPLING RATE TO 'STOP'"/"*"/"*"
     5  "SELECT COMPUTER MODE FOR DVM1. "/"*"/"*")
5305  FORMAT("*"/"*AT END OF RUN SET SWITCH REGISTER BIT 15 TO TER" 
     1"MINATE"/"***************************************************") 
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(10X,"ASSUMED FLUX CALIBRATION FOR NOISE TUBE = ",F7.3) 
5130  FORMAT(//6HSOURCE,2X,6HFLUX A,2X,6HFLUX B,2X,4HFLUX,2X, 
/    *5HSIGMA,2X,4HNCAL,2X,4HNDEC,2X,6HOFFSET,3X,2HRA,4X,3HDEC) 
5260  FORMAT(4A2,40X,F6.3,2F8.3)
5270  FORMAT(4A2,F7.2,F7.2,F7.2,6X,F6.0)
5180  FORMAT(F7.2,F7.2,F7.2,F6.3,14X,F5.3,2F8.3)
      END 
      SUBROUTINE SIDTM(UT,ST,UTIME,STO,ST1) 
******************************************* 
      DIMENSION ITIME(5),ATIME(5) 
      ICODE=11
      CALL EXEC(ICODE,ITIME)
      DO 10 J3=1,5
      ATIME(J3)=ITIME(J3) 
 10   DUMMY=0 
      UTIME=(ATIME(4)+ATIME(3)/60.+ATIME(2)/3600.+ATIME(1)/360000.) 
      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 
      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 48  J1=1,ILENG 
*     
*     Sample once per second for five seconds and take median 
* 
             CALL EXEC (2,113B,55B,1,5,5) 
             WRITE (35,'("*S1")') 
          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
* 
 48   CONTINUE
* 
*      Average last 12 points of baseline for noise cal reference 
      AVG1 = 0.0
      DO 50 J1 = ILENG-11, ILENG
      AVG1 = AVG1 + ABUF(J1)
 50   CONTINUE
      AVG1 = AVG1/12
      RETURN
      END 
      SUBROUTINE BEAMF(ITM,ABUF,RDEC,SPA,SPB,SIGMA,CENT,I,IPLOT,DHA 
     *,SORSNM,FILNAM) 
******************************************************************* 
      DIMENSION ABUF(200),ANEW(200),SUM(2),ITIME(5) 
      INTEGER SORSNM(4),FILNAM(3) 
******FIT STRAIGHT LINE TO BASE LINE REGION TO REMOVE DRIFT 
   5  ILENG=ITM/5+1 
      ISP=1 
      ILP=24
      IC=ILENG
      N=48
      U=COS(RDEC/57.3)
      SECD = 1/U
      BMWTH = 0.10448*U 
      ALPHA = BMWTH*BMWTH 
******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 BY LOCATING CENTROID OF BEAM RESPONSE 
      NOPTS = 8*SECD
      MID = CENT - NOPTS
  75  SUM(1) = 0. 
      SUM(2) = 0. 
      MID = MID + 1 
      INC = 0 
      DO 85 J = 1,2 
      N1 = MID - NOPTS + INC
      N2 = MID - 1 + INC
* 
*     Test that N1 and N2 are within bounds for sampled data
* 
      DO 80 N = N1,N2 
      SUM(J) = SUM(J) + ANEW(N) 
  80  CONTINUE
      INC = NOPTS + 1 
  85  CONTINUE
      IF (SUM(1).GT.SUM(2)) GO TO 75
      ZE  = (SUM(1) - SUM(2))/(SUM(1) + SUM(2)) 
      SUM(1) = 0. 
      SUM(2) = 0. 
      DO 90 J = 1,2 
      DO 90 N = 1,NOPTS 
      N1 = N  -  1 + 2*(J-1)
      BMFN = EXP(-ALPHA*N1*N1)
      SUM(J) = SUM(J) + BMFN
  90  CONTINUE
      CONV = (SUM(1) - SUM(2))/(SUM(1) + SUM(2))
      DHA = (DHA -MID - ZE/CONV)*0.0208 
      BMCNT = MID - ZE/CONV 
      IBMCNT = BMCNT
******DETERMINE FLUX DENSITY BY MEANS OF LEAST SQUARES BEAM FIT 
      SPA = 0.
      BMFNSQ = 0. 
      N1 = IBMCNT - NOPTS 
      N2 = IBMCNT + NOPTS 
      DO 100 N = N1,N2
      BMFN = EXP(-ALPHA*(BMCNT-N)**2) 
      SPA = SPA + BMFN*ANEW(N)
      BMFNSQ = BMFNSQ + BMFN*BMFN 
 100  CONTINUE
      SPA =  SPA/BMFNSQ 
      SPB = SPA 
******ANALYTICAL CURVE
      DO 220 N = N1,N2
      AN = BMCNT - N
      BMFN = EXP(-ALPHA*AN*AN)
      ANEW(N) = SPA*BMFN
 220  CONTINUE
******BEGIN TO PLOT 
      IF ( IPLOT.GT.1) GO TO 145
      IF ( I.GT.1) GO TO 145
******ADVANCE CHART TWO PAGES   - ONE PAGE - MJG 1982 09 30 
      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,-3) 
 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
      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,"_") 
      KYMAX = 9999
      WRITE (12) -1,1,KX,KYMAX
      END 
      SUBROUTINE LSQB(ISP,ILP,ABUF,A,B,IC,N)
********************************************
      DIMENSION ABUF(200) 
      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 
      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(12,0,2,0,DELAY) 
      CALL EXEC (2,113B,55B,1,5,5)
      WRITE (35,'("*S1")')
      DO 10 I6=1,NOPTS
         READ (35,*) RDVM 
  10  DVM1 = DVM1 + RDVM
      AVG = -(DVM1/NOPTS)*100000
      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 
      END$