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$