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 * *********************************************************