PROGRAM VPLOT C C This is the second, of three steps of data reduction for the C amplitude data from the Mount Wilson Optical Interferometer. C C FRITZ averages the raw data files and produces 'reasonable' sized C files, a few megabytes per night. C C VPLOT reads the visibility files output by FRITZ, makes several plots C of parameters related to seeing and system performance, makes C plots of the visibilities, allows some simple averaging, editing C and calibration. The output is a small file for use by OICAL. C C OI2CIT reads the output of VPLOT, re-calibrates the visibilities and C enters the data into the Cal-tek model fitting package. C C VERSION 4.02 19 July 1989 C C 4.01 Added automatic editing C 4.02 Added phase rms calibration C IMPLICIT UNDEFINED(A-Z) REAL *4 PI, PISQ PARAMETER ( PI = 3.1415926535 ) PARAMETER (PISQ = PI*PI ) CHARACTER *12 VERSN, VERSD PARAMETER ( VERSN = 'VERSION 4.02' ) PARAMETER ( VERSD = ' 1989 Jul 12' ) INTEGER *4 NSUM, IERR, IREC, ISYS INTEGER *4 N, I, J, K, L, MONTH, DAY, YEAR, VCNT, NPTS, JBASE INTEGER *4 NDATA, NDARK, ICHAN, IPLOT(10), FITTED(2), LINES(5) INTEGER *4 POINT(5), IFILT(4), BEGIN, NAME(4), SPECTRUM, NC INTEGER *4 MASK(4), NSTAT, NPAR(4) REAL *4 INTTIME, CALAVG, CALRMS, DEN, NPHOT, VAMP(4), VSQ REAL *4 CONST, SIGMA, FLUX(4), NAVG, NRMS, VAVG, VRMS REAL *4 VMEAN, MEAN, PAR, RV, PMRA, PMDEC, MAGLIM REAL *4 TEMP(9), NEWDIAM REAL *8 RAM, DECM C Global variables. INCLUDE 'VPLOT.INC' C Filenames for laserjet plots INCLUDE 'HARDCOPY.INC' REAL *4 XPLOT(CNT), YPLOT(CNT,7) C CHARACTER * 1 ANS CHARACTER * 3 MNTH(12) CHARACTER * 8 ICMD, FILTERS(0:32) CHARACTER *10 DEV CHARACTER *12 OUTNAME CHARACTER *40 TITLE(6) CHARACTER *80 STRING C LOGICAL HARD C C The input data format C INTEGER *2 IA(38), ICODE, ISTAR, IBASE, TAU0, DTAU0 INTEGER *2 TSAMPLE, IMON, IDAY, IYEAR, JITTR, ILOCK INTEGER *4 IDARK(4), FILT(4), DJTR REAL *8 TIME0, DLAS REAL *4 ANUM(6), STROKE, ADEN(4), AVAR(4), HA0, DEC0, ZD0 REAL *4 TFIRST, TLAST, PHI12, PHIRMS C DATA REC TIME REC DARK REC EQUIVALENCE (IA( 1), TIME0 ), $ (IA( 5), ANUM(1), ICODE ), $ (IA( 6), ISTAR ), $ (IA( 7), IBASE ), $ (IA( 8), TAU0, DTAU0 ), $ (IA( 9), TSAMPLE, DLAS ), $ (IA(10), IMON ), $ (IA(11), IDAY ), $ (IA(12), IYEAR ), $ (IA(13), STROKE, IDARK(1)), $ (IA(15), FILT(1) ) EQUIVALENCE (IA(17), ADEN(1) ), $ (IA(21), HA0 ), $ (IA(23), DEC0 ), $ (IA(25), AVAR(1), ZD0 ), $ (IA(27), DJTR ), $ (IA(29), TFIRST ), $ (IA(31), TLAST ), $ (IA(33), JITTR, ILOCK ), $ (IA(35), PHIRMS ), $ (IA(37), PHI12 ) C********************************************************************** DATA NPLOTS, NSTAT / 0, 0 / DATA PLOTFILE, PLOT_DEV / 'PLOT.NNN', '/LPM' / DATA MNTH / 'JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', $ 'JUL', 'AUG', 'SEP', 'OCT', 'NOV', 'DEC' / DATA FITTED / 1, 10 / DATA LINES / 11, 22, 33, 44, 55 / DATA POINT / 1, 2, 3, 4, 5 / DATA DEV / '/EGA ' / C DATA HARD / .FALSE. / DATA MASK / 1, 2, 4, 8 / DATA FILTERS / '700/999', $ '400/10 ', '400/25 ', '400/40 ', $ '450/10 ', '450/25 ', '450/40 ', $ '500/10 ', '500/25 ', '500/40 ', $ '550/10 ', '550/25 ', '550/40 ', $ '600/10 ', '600/25 ', '600/40 ', $ '650/10 ', '650/25 ', '650/40 ', $ '700/10 ', '700/25 ', '700/40 ', $ '750/10 ', '750/25 ', '750/40 ', $ '800/10 ', '800/25 ', '800/40 ', $ '687/40 ', '397/1 ', '410/1 ', $ '434/1 ', '486/1 ' / C---------------------------------------------------------------------- WRITE(6,*) ' VPLOT averages, edits, and plots visibilities' WRITE(6,*) ' ', VERSN, ' ', VERSD ISCAN = 0. JBASE = 0 NSTAR = 0 NBASE = 0 FNAME = ' ' 10 CONTINUE IF ( FNAME .EQ. ' ' ) THEN WRITE(6,*) ' No data, you must enter a data file. ' ICMD = 'READ' ELSE WRITE (6,5000) READ (5,'(A)',END=10) ICMD CALL CAPS ( ICMD ) END IF 5000 FORMAT( / ' ? for help, command :', $ ) IF ( (ICMD(1:5) .EQ. 'HELP ').OR.(ICMD(1:1).EQ.'?') ) THEN WRITE(6,'(A)') $' Usual sequence of commands ', $' READ inport data from a file ', $' SEEING to get some idea of data quality ', $' AUTOPED point calibration ', $' NCAL check alignment of fibers and sky transparancy ', $' DARK background editor ', $' PHOTOM check for vignetting ', $' SYSVIS determine system calibration ', $' VCAL apply calibration to the data ', $' EXPORT output data for model fitting ', $' COMMENT add comment to summary file ', $' SUMMARY print/display summary file ', $' ', $' INDEX prints an index of other commands ' C ELSE IF ( ICMD(1:5) .EQ. 'INDEX' ) THEN WRITE(6,'(A)') $' INDEX this listing ? for help ', $' DELAY plots vs delay PED manual point editor ', $' SED manual scan editor STAT gives mean vis and rms ', $' DIAM assign star diams UNCAL resets cal table ', $' EXAM examine some data SIGMA examines uncertainties ', $' NPLOT plot photon level VPLOT nice vis plots ', $' VIS plot visibilities OUTCAL Output cal tables ', $' HARD set hardcopy device type ', $' SAVE Save the scan tables FITDIA Help fitting diameters ', $' Q quit ', $' ' C ELSE IF ( ICMD(1:5) .EQ. 'READ ' ) THEN IF ( FNAME .NE. ' ' ) THEN WRITE(6,*) ' The current data will be lost unless you ', $ ' save it NOW! ' WRITE(6,'(A,$)') ' Do you want to save the current data? ' READ (5,'(A)') ANS CALL CAPS ( ANS ) IF ( ANS .EQ. 'Y' ) THEN CALL VSAVE END IF END IF CALL VREAD C ELSE IF ( ICMD(1:6) .EQ. 'PHOTOM' ) THEN CALL PHOTOM C C Determine mean PHOTON FLUX FOR EACH STAR C ELSE IF ( ICMD(1:5) .EQ. 'NCAL ' ) THEN CALL NCALIB C C Apply calibration curves to the observed visibility C ELSE IF ( ICMD(1:5) .EQ. 'VCAL ' ) THEN CALL APPLY_CAL ( NPAR ) CALL GET_SYS ( NPAR ) WRITE(6,'(A,4F6.3)') ' Calibration errors are ', SYS_CAL ELSE IF ( ICMD(1:5) .EQ. 'PCAL ' ) THEN WRITE(6,*) ' INPUT MAXIMUM ALLOWED PHASE RMS ' READ (5,*,ERR=10) PHIMAX WRITE(6,*) ' INPUT MINIMUM ALLOWED PHASE RMS ' READ (5,*,ERR=10) PHIMIN WRITE(6,*) ' INPUT PHASE RMS DATA SHOULD BE CORRECTED TO ' READ (5,*,ERR=10) PHIREF WRITE(6,*) ' INPUT MAGNITUDE OF CORRELATION, dV/dPHI ', $ ' for all 4 channels ' READ (5,*,ERR=10) (DVDP(J),J=1,4) WRITE(6,*) ' THESE VALUES WILL BE STORED WHEN YOU SAVE', $ ' THE CALIBRATION TABLES ' WRITE(6,*) ' YOU MUST RE-READ THE DATA FOR THEM TO HAVE', $ ' AN AFFECT ' ELSE IF ( ICMD(1:6) .EQ. 'UNCAL ') THEN C C Determine which calibrations are set C WRITE(6,*) ' CHAN TIME SEEING ZD MA' DO J = 1, 4 TITLE(1) = ' ' WRITE(TITLE(1)(1:5),'(I5)') J CALL CHECKCAL( TCAL(1,J), NFIT, TITLE(1)( 6:12) ) CALL CHECKCAL( HCAL(1,J), NFIT, TITLE(1)(13:19) ) CALL CHECKCAL( ZCAL(1,J), NFIT, TITLE(1)(20:26) ) CALL CHECKCAL( MCAL(1,J), NFIT, TITLE(1)(27:33) ) WRITE(6,*) TITLE(1) END DO WRITE(6,*) ' WHICH CALIBRATION DO YOU WANT TO RESET ?' WRITE(6,*) ' T for time ' WRITE(6,*) ' S for seeing ' WRITE(6,*) ' Z for zenith angle ' WRITE(6,*) ' M for mirror angle ' WRITE(6,*) ' A for all ' READ (5,'(A)') ANS CALL CAPS ( ANS ) IF ( (ANS .EQ. 'T').OR.(ANS .EQ. 'A') ) THEN DO ICHAN = 1, 4 NFREE(ICHAN,1) = 0 DO I = 1, NFIT TCAL(I,ICHAN) = 1. END DO END DO WRITE(6,*) ' TIME calibration reset ' END IF IF ( (ANS .EQ. 'S').OR.(ANS .EQ. 'A') ) THEN DO ICHAN = 1, 4 NFREE(ICHAN,2) = 0 DO I = 1, NFIT HCAL(I,ICHAN) = 1. END DO END DO WRITE(6,*) ' SEEING calibration reset ' END IF IF ( (ANS .EQ. 'Z').OR.(ANS .EQ. 'A') ) THEN DO ICHAN = 1, 4 NFREE(ICHAN,3) = 0 DO I = 1, NFIT ZCAL(I,ICHAN) = 1. END DO END DO WRITE(6,*) ' ZENITH DISTANCE calibration reset ' END IF IF ( (ANS .EQ. 'M').OR.(ANS .EQ. 'A') ) THEN DO ICHAN = 1, 4 NFREE(ICHAN,4) = 0 DO I = 1, NFIT MCAL(I,ICHAN) = 1. END DO END DO WRITE(6,*) ' MIRROR ANGLE calibration reset ' END IF WRITE(6,*) ' Remember, You must do a VCAL for this ', $ ' to have an effect ' ELSE IF ( ICMD(1:6) .EQ. 'OUTCAL' ) THEN OPEN ( UNIT=18, FILE=FNAME(1:IROOT)//'.TAB', STATUS='UNKNOWN', $ IOSTAT=IERR ) IF ( IERR .NE. 0 ) THEN WRITE(6,*) ' Error opening calibration table file' WRITE(6,*) ' TOO BAD ' ELSE WRITE(18,1800) DO 250 I = 1, NFIT WRITE(18,1810) I, $ (TCAL(I,J),J=1,4), (HCAL(I,J),J=1,4), $ (ZCAL(I,J),J=1,4), (MCAL(I,J),J=1,4) 250 CONTINUE CLOSE(UNIT=18,STATUS='KEEP' ) WRITE(6,*) ' Calibration tables written to file ', $ FNAME(1:IROOT)//'.TAB' END IF 1800 FORMAT ( ' I TIME HA ZD MIRROR ANGLE ' ) 1810 FORMAT ( I5, 4(2X, 4F6.3) ) ELSE IF ( ICMD(1:5) .EQ. 'STAT ' ) THEN WRITE(6,*) ' INPUT CHANNEL NUMBER:' READ (5,*,ERR=10) ICHAN IF (( ICHAN .LT. 1 ).OR.(ICHAN.GT.4)) THEN WRITE(6,*) ' BAD CHANNEL NUMBER ' GO TO 10 END IF NSTAT = NSTAT + 1 WRITE(OUTNAME, '(A5,I3.3,4X)') 'STAT.', NSTAT OPEN ( 33, FILE = OUTNAME ) WRITE( 6,2220) ICHAN, FILTERID(ICHAN), DATE WRITE(33,2220) ICHAN, FILTERID(ICHAN), DATE DO 300 J = 1, NBASE DO 290 I = 1, NSTAR VCNT = 0 VMEAN = 0. MEAN = 0. VRMS = 0. NPHOT = 0. CALAVG= 0. CALRMS= 0. NC = 0 DO 280 K = SCAN0, SCAN1 IF ( (BASE(K) .EQ. BLIST(J) ) .AND. $ (STAR(K) .EQ. SLIST(I) ) .AND. $ (IAND(OKSCAN(K),MASK(ICHAN)) .EQ. 0 ) ) THEN VCNT = VCNT + 1 VMEAN = VMEAN + VIS(K,ICHAN) VRMS = VRMS + VIS(K,ICHAN)*VIS(K,ICHAN) NPHOT = NPHOT + BDEN(K,ICHAN) - D(K,ICHAN) MEAN = MEAN + VEST(K,ICHAN) IF ( VCAL(K,ICHAN) .NE. 0. ) THEN CALAVG = CALAVG + VCAL(K,ICHAN) CALRMS = CALRMS + VCAL(K,ICHAN)*VCAL(K,ICHAN) NC = NC + 1 END IF END IF 280 CONTINUE IF ( VCNT .EQ. 0 ) THEN XPLOT(I) = 0. YPLOT(I,1) = 0. ELSE VMEAN = VMEAN / FLOAT(VCNT) MEAN = MEAN / FLOAT(VCNT) VRMS = SQRT ( VRMS / FLOAT(VCNT) - VMEAN*VMEAN ) NPHOT = NPHOT / FLOAT(VCNT) END IF IF ( NC .GT. 0 ) THEN CALAVG = CALAVG / FLOAT(NC) CALRMS = SQRT ( CALRMS/FLOAT(NC) - CALAVG*CALAVG ) END IF WRITE( 6,2230) SLIST(I), BLIST(J), VMEAN, VRMS, VCNT, $ MEAN, CALAVG, CALRMS, NC, NPHOT WRITE(33,2230) SLIST(I), BLIST(J), VMEAN, VRMS, VCNT, $ MEAN, CALAVG, CALRMS, NC, NPHOT 2220 FORMAT ( ' SCAN AVERAGES FOR CHANNEL ', I3, ' = ', A, 2X, A $ / 10X, ' STAR B Vrms #OBS', $ ' Vest rms #CAL #PHOT' ) 2230 FORMAT ( I15, I3, 2F7.3, I6, 3F7.3, I6, F7.1 ) IF ( VCNT .LT. 2 ) THEN XPLOT(I) = 0. YPLOT(I,1) = 0. ELSE XPLOT(I) = ANGDIA(I) YPLOT(I,1) = VMEAN END IF 290 CONTINUE CLOSE (UNIT=33 ) WRITE(6,*) ' WAITING TO PLOT ' READ (5,'(A)',END=295) ANS 295 CONTINUE C TITLE(1) = 'ESTIMATED ANGULAR DIAMETER (MAS)' C TITLE(2) = 'MEAN VISIBILITY' C TITLE(3) = DATE C IF ( BLIST(J) .EQ. 1 ) THEN C TITLE(4) = '12 METER N-S BASELINE' C ELSE IF ( BLIST(J) .EQ. 2 ) THEN C TITLE(4) = '12 METER E-S BASELINE' C ELSE C WRITE(TITLE(4),'(A,I3,24X)') 'VB BASELINE #', BLIST(J) C END IF C WRITE(TITLE(5),2202) 'CHANNEL', ICHAN, FILTERID(ICHAN) C 2202 FORMAT ( A, I4, 3X, A ) C TITLE(6) = 'ONE POINT PER STAR' C CALL EZMXY ( XPLOT, YPLOT, NSTAR, CNT, 1, POINT, DEV, TITLE) 300 CONTINUE ELSE IF ( ICMD(1:5) .EQ. 'SIGMA' ) THEN WRITE(6,'(A,$)') ' INPUT CHANNEL NUMBER:' READ (5,*,ERR=10) ICHAN IF (( ICHAN .LT. 1 ).OR.(ICHAN.GT.4)) THEN WRITE(6,*) ' BAD CHANNEL NUMBER ' GO TO 10 END IF WRITE(OUTNAME, '(A6,I3.3)') 'SIGMA.', ICHAN OPEN ( 33, FILE = OUTNAME ) WRITE( 6,2250) ICHAN, FILTERID(ICHAN), DATE WRITE(33,2250) ICHAN, FILTERID(ICHAN), DATE DO K = 1, ISCAN IF ( IAND(OKSCAN(K),MASK(ICHAN) ) .NE. 0 ) THEN WRITE( 6,2260) K, STAR(K), VIS(K,ICHAN), $ VCAL(K,ICHAN), VEST(K,ICHAN), VSIG(K,ICHAN), $ SIG_PHOT(K,ICHAN) WRITE(33,2260) K, STAR(K), VIS(K,ICHAN), $ VCAL(K,ICHAN), VEST(K,ICHAN), VSIG(K,ICHAN), $ SIG_PHOT(K,ICHAN) END IF END DO 2250 FORMAT ( ' SCAN SUMMARY FOR CHANNEL ', I3, ' = ', A, 2X, A,/ $ / ' SCAN STAR VISIBILITY SQUARED UNCERTAINTY ', $ / ' NO NO OBS CAL EST STAT PHOT ' ) 2260 FORMAT ( 2I5, 2X, 3F7.3, 2X, 2F7.3 ) ELSE IF ( ICMD(1:7) .EQ. 'COMMENT' ) THEN WRITE(6,'(I3,A)') (I,COMMENT(I),I=1,5) WRITE(6,'(A,$)') ' Modify which comment (1...5)?' READ (5,*) I IF ( (I.LT.1).OR.(I.GT.5) ) THEN WRITE(6,*) ' BAD COMMENT LINE ' ELSE WRITE(6,*) ' ENTER TEXT ' WRITE(6,'(I3,A)') I, COMMENT(I) WRITE(6,'(I3,$)' ) I READ (5,'(A)') COMMENT(I) END IF ELSE IF ( ICMD(1:7) .EQ. 'SUMMARY' ) THEN C C Gives summary of data reduction done to date, the calibration C finished and some data quality indicators. C CALL GET_SYS ( NPAR ) STRING = FNAME(1:IROOT) // '.SUM' OPEN ( UNIT=34, FILE=STRING, STATUS='UNKNOWN' ) WRITE(34,2310) DATE, FNAME, ISCAN, HOURS(1), HOURS(ISCAN), $ FORMATID WRITE(34,2320) (I, FILTERID(I), MAG0(I), CTERM(I), $ NCAL_RMS(I), I=1,4) WRITE(34,2325) (I, (D_ORDER(J,I),J=1,2), D_RMS(I), $ D_CHISQ(I), (D_RANGE(J,I),J=1,2),I=1,4) WRITE(34,2330) PHIPEAK WRITE(34,2335) (I, (ZD_COEF(J,I),J=0,3),I=1,4) WRITE(34,2340) NCALSTAR, (CALSTAR(I),I=1,NCALSTAR) WRITE(34,2360) SEEING_TYPE, FILTERID, $ ((SEEING_COEF(J,I),I=1,4),J=1,3), NCALSCAN, $ CAL_CHISQ, CAL_RMS, SYS_CAL, CAL_ERR WRITE(34,2370) ((NFREE(I,J),I=1,4),J=1,4), NPAR WRITE(34,2380) C WRITE( 6,2310) DATE, FNAME, ISCAN, HOURS(1), HOURS(ISCAN), $ FORMATID WRITE( 6,2320) (I, FILTERID(I), MAG0(I), CTERM(I), $ NCAL_RMS(I), I=1,4) WRITE( 6,2325) (I, (D_ORDER(J,I),J=1,2), D_RMS(I), $ D_CHISQ(I), (D_RANGE(J,I),J=1,2),I=1,4) WRITE( 6,2330) PHIPEAK WRITE( 6,2335) (I, (ZD_COEF(J,I),J=0,3),I=1,4) WRITE( 6,2340) NCALSTAR, (CALSTAR(I),I=1,NCALSTAR) WRITE( 6,2360) SEEING_TYPE, FILTERID, $ ((SEEING_COEF(J,I),I=1,4),J=1,3), NCALSCAN, $ CAL_CHISQ, CAL_RMS, SYS_CAL, CAL_ERR WRITE( 6,2370) ((NFREE(I,J),I=1,4),J=1,4), NPAR WRITE( 6,2380) C DO J = 1, 4 TITLE(1) = ' ' WRITE(TITLE(1)(1:5),'(I5)') J CALL CHECKCAL( TCAL(1,J), NFIT, TITLE(1)( 6:12) ) CALL CHECKCAL( HCAL(1,J), NFIT, TITLE(1)(13:19) ) CALL CHECKCAL( ZCAL(1,J), NFIT, TITLE(1)(20:26) ) CALL CHECKCAL( MCAL(1,J), NFIT, TITLE(1)(27:33) ) WRITE(34,'(A)') TITLE(1) WRITE( 6,'(A)') TITLE(1) END DO COMMENT WRITE( 6,'(I3,A70)') ( I, COMMENT(I), I=1,5 ) WRITE(34,'(I3,A70)') ( I, COMMENT(I), I=1,5 ) CLOSE ( UNIT = 34 ) 2310 FORMAT ( 20X, 'DATA SUMMARY FOR ', A, $ / 20X, 'DATA FILE NAME ', A, $ // ' Number of scans ', I5, $ / ' range of ut times =', F10.3, ' - ', F10.3, $ / ' format id number = ', I5 // ) 2320 FORMAT ( // 20X, ' PHOTOMETRIC CALIBRATION ' $ / ' CHAN FILTER MAG0 C TERM rms' $ / (I5, 2X, A8, 2X, F8.3, 1X, F8.3, 1X, F8.3) ) 2325 FORMAT ( // 20X, ' DARK COUNT CALIBRATION ' $ / ' CHAN ORDER RMS CHISQ RANGE ' $ / (I5, I6, ',', I3, F8.3, 1X, F4.1, F7.2, ' -', F5.2) ) 2330 FORMAT ( / 5X, ' PHI RMS PEAK = ', F8.2, $ / 5X, ' LASER RMS PEAK = ', F8.2 ) 2335 FORMAT ( // 5X, ' Zenith angle calibration ' $ / ' CHAN ORDER A0 A1 A2 ', $ 4(/ I5, 4F8.5) ) 2340 FORMAT ( / 5X, ' Number of calibration stars = ', I5, $ / ' calibration star list ', 10(4I5,2x,4I5,/,25X) ) 2360 FORMAT ( // 20X, ' VISIBILITY CALIBRATION ' $ / ' Seeing calibration type = ', I5 $ / ' CHAN1 CHAN2 CHAN3 CHAN4 ', $ / ' FILTER ', 4(1X,A8,1X), $ / ' V0 ', 4(F9.5,1X), $ / ' A ', 4(F9.5,1X), $ / ' EXP ', 4(F9.2,1X), $ / ' No. scans ', 4(I9, 1X), $ / ' Chi square', 4(F9.2,1X), $ / ' rms of fit', 4(F9.3,1X), $ / ' system rms', 4(F9.3,1X), $ / ' Export rms', 4(F9.3,1X) ) 2370 FORMAT ( 10X, ' Number of parameters in the fits ' $ / ' CAL TYPE 1 2 3 4' $ / ' time ', 4I5, $ / ' seeing ', 4I5, $ / ' zenith distance ', 4I5, $ / ' mirror angle ', 4I5, $ / ' total parameters', 4I5 ) 2380 FORMAT (/ ' CHAN TIME SEEING ZD MA' ) ELSE IF ( ICMD(1:5) .EQ. 'DARK ' ) THEN CALL DARKPLOT ELSE IF ( ICMD(1:6) .EQ. 'FITDIA' ) THEN CALL FITDIA ELSE IF ( ICMD(1:4) .EQ. 'DIAM' ) THEN WRITE(6,*) ' ENTER STAR NUMBER ' READ (5,*,END=10) I DO 330 J = 1, NSTAR IF ( I .EQ. SLIST(J) ) THEN K = J GO TO 335 END IF 330 CONTINUE WRITE(6,*) ' STAR IS NOT IN LIST ' GO TO 10 335 CONTINUE WRITE(6,*) ' CURRENT DIAMETER IS ', ANGDIA(K) WRITE(6,*) ' INPUT NEW DIAMETER: ' READ (5,*,END=10) NEWDIAM ANGDIA(K) = NEWDIAM WRITE(6,1700) WAVE DO 350 I = 1, ISCAN IF ( STAR(I) .NE. SLIST(K) ) GO TO 350 DO J = 1, 4 CALL QCAL ( BASE(I), HA(I), DEC(I), WAVE(J), $ ANGDIA(K), VEST(I,J) ) END DO WRITE(6,1705) STAR(I), BASE(I), HA(I), DEC(I), $ ANGDIA(J), (VEST(I,L),L=1,4) 350 CONTINUE 1700 FORMAT ( 10X, ' wavelengths = ', 4F8.0, $ / ' STAR B HA DEC DIA CH 1 CH 2 CH 3', $ ' CH 4' ) 1705 FORMAT ( I5, I3, 2F6.2, F6.1, 4F8.3 ) ELSE IF ( ICMD(1:6) .EQ. 'NPLOT ' ) THEN WRITE(6,*) ' INPUT CHANNEL NUMBER: ' READ (5,*,END=450) ICHAN IF ( (ICHAN.LT.1).OR.(ICHAN.GT.4) ) GO TO 450 WRITE(6,*) ' PLOT DARK COUNT TOO (Y/N)?' READ (5,'(A)') ANS CALL CAPS (ANS) IF ( ANS .EQ. 'Y' ) THEN J = 2 ELSE IF ( ANS .EQ. 'N' ) THEN J = 1 ELSE GO TO 450 END IF TITLE(1) = 'TIME (HOURS)' WRITE( TITLE(2),1304) ICHAN, FILTERID(ICHAN) DO 440 K = 1, NSTAR TITLE(3) = DATE WRITE(TITLE(4),1300) SLIST(K), STARNAME(K) TITLE(5) = 'DASHED LINE IS 10 TIMES DARK COUNT' NPTS = 0 NAVG = 0. NRMS = 0. DO 420 I = 1, ISCAN IF ( STAR(I) .NE. SLIST(K) ) GO TO 420 IF ( IAND(OKSCAN(I),MASK(ICHAN)).NE. 0 ) GO TO 420 NPTS = NPTS + 1 XPLOT(NPTS) = HOURS(I) YPLOT(NPTS,1) = BDEN(I,ICHAN) - D(I,ICHAN) YPLOT(NPTS,2) = 10. * D(I,ICHAN) NAVG = NAVG + YPLOT(NPTS,1) NRMS = NRMS + YPLOT(NPTS,1)*YPLOT(NPTS,1) 420 CONTINUE IF ( NPTS .LT. 3 ) GO TO 440 NAVG = NAVG / FLOAT(NPTS) NRMS = SQRT ( NRMS/FLOAT(NPTS) - NAVG*NAVG ) WRITE(TITLE(6),1306) NAVG, NRMS CALL EZMXY ( XPLOT, YPLOT, NPTS, CNT, J, $ LINES, DEV, TITLE) 440 CONTINUE 450 CONTINUE 1300 FORMAT ( 'FK5', I5, ' = ', A8, 21X ) 1304 FORMAT ( 'PHOTONS / 4 ms CHANNEL',I2, ' = ', A7 ) 1305 FORMAT ( 'CHANNEL', I3, ' = ', A7, 20X ) 1306 FORMAT ( ' = ', F7.2, 5X, ' RMS = ', F7.2, 8X ) 1307 FORMAT ( 'VIS SQUARED CHANNEL',I2, ' = ', A7 ) 1308 FORMAT ( ' = ', F7.2, 2X, ' VRMS = ', F7.2, 7X ) ELSE IF ( ICMD(1:4) .EQ. 'VIS ' ) THEN WRITE(6,*) ' INPUT CHANNEL NUMBER: ' READ (5,*,END=469) J IF ( (J.LT.1).OR.(J.GT.4) ) GO TO 450 TITLE(1) = 'TIME (HOURS)' WRITE( TITLE(2),1307) J, FILTERID(J) DO 465 K = 1, NSTAR TITLE(3) = DATE WRITE(TITLE(4),1300) SLIST(K), STARNAME(K) TITLE(5) = ' ' NPTS = 0 VAVG = 0. VRMS = 0. DO 460 I = 1, ISCAN C IF ( STAR(I) .EQ. SLIST(K) ) THEN NPTS = NPTS + 1 XPLOT(NPTS) = HOURS(I) YPLOT(NPTS,1) = VIS(I,J) VAVG = VAVG + VIS(I,J) VRMS = VRMS + VIS(I,J)*VIS(I,J) END IF 460 CONTINUE IF ( NPTS .LT. 3 ) GO TO 465 VAVG = VAVG / FLOAT(NPTS) VRMS = SQRT ( VRMS/FLOAT(NPTS) - VAVG*VAVG ) WRITE(6,1308) VAVG, VRMS WRITE(TITLE(6),1308) VAVG, VRMS CALL EZMXY ( XPLOT, YPLOT, NPTS, CNT, 1, $ LINES, DEV, TITLE) 465 CONTINUE 469 CONTINUE ELSE IF ( ICMD(1:6) .EQ. 'SEEING' ) THEN CALL SEEING ELSE IF ( ICMD(1:4) .EQ. 'EXAM' ) THEN C C Examine data for an individual star C WRITE(6,*) ' INPUT STAR NUMBER: ' READ (5,*,END=10) J WRITE(6,1500) DO 600 I = 1, ISCAN IF ( STAR(I) .NE. J ) GO TO 600 WRITE(6,1510) I, $ ( K, BDEN(I,K), D(I,K), VIS(I,K), K=1,4 ) 600 CONTINUE 1500 FORMAT (' SCAN CHAN',12X,'N ',10X,'DARK ',10X,'VIS ' ) 1510 FORMAT ( I5, 4(I5, 3F15.5, /, 5X ) ) ELSE IF ( ICMD(1:6) .EQ. 'DELAY ' ) THEN TITLE(1) = 'DELAY' TITLE(2) = 'DELAY LINE JITTER' TITLE(3) = DATE TITLE(4) = ' ' TITLE(5) = ' ' TITLE(6) = ' ' CALL EZMXY ( DELAY, DJITTER, ISCAN, CNT, 1, POINT, $ DEV, TITLE ) TITLE(2) = 'TRACKING JITTER' CALL EZMXY ( DELAY, TJITTER, ISCAN, CNT, 1, POINT, $ DEV, TITLE ) TITLE(2) = 'PHASE RMS (MICRONS)' CALL EZMXY ( DELAY, PHI(1,1), ISCAN, CNT, 1, POINT, $ DEV, TITLE ) TITLE(2) = 'PHASE AUTOCORRELATION, 4ms LAG' CALL EZMXY ( DELAY, PHI(1,2), ISCAN, CNT, 1, POINT, $ DEV, TITLE ) TITLE(2) = 'SCINTILLATION INDEX' CALL EZMXY ( DELAY, SCIN(1,1), ISCAN, CNT, 4, POINT, $ DEV, TITLE ) ELSE IF ( ICMD(1:6) .EQ. 'VPLOT ' ) THEN C Make nice, multi star, visibility plots CALL STARPLOT ( title, DEV ) ELSE IF ( ICMD(1:6) .EQ. 'SYSVIS' ) THEN CALL SYSVIS ( DEV ) C ELSE IF ( ICMD(1:5) .EQ. 'SOFT ' ) THEN C DEV = '/EGA ' C HARD = .FALSE. ELSE IF ( ICMD(1:5) .EQ. 'HARD ' ) THEN WRITE(6,*) ' If files PLOTS.nnn exist in current directory' WRITE(6,*) ' Plots will not be generated.' WRITE(6,*) ' Set hardcopy device. Input 1 to 6 ' WRITE(6,*) ' 1 /LJL 2 /LJM 3 /LJH ' WRITE(6,*) ' 4 /LPM 5 /LPH ' READ (5,*,END=10) I IF ( I .EQ. 1 ) THEN DEV = '/LJL ' ELSE IF ( I .EQ. 2 ) THEN DEV = '/LJM ' ELSE IF ( I .EQ. 3 ) THEN DEV = '/LJH ' ELSE IF ( I .EQ. 4 ) THEN DEV = '/LPM ' ELSE IF ( I .EQ. 5 ) THEN DEV = '/LPH ' ELSE WRITE(6,*) ' BAD DEVICE TYPE ' DEV = '/LJM ' END IF COMMENT HARD = .TRUE. PLOT_DEV = DEV ELSE IF ( ICMD(1:3) .EQ. 'SED' ) THEN CALL SEDIT ELSE IF ( ICMD(1:4) .EQ. 'PED' ) THEN CALL PEDIT 650 CONTINUE WRITE(6,'(A,$)') $ ' DO YOU WANT TO RE-CALCULATE THE SCAN AVERAGES ?' READ (5,'(A)') ANS CALL CAPS ( ANS ) IF ( ANS .EQ. 'Y' ) THEN CALL VSAVE CALL SCANCALC ELSE IF ( ANS .NE. 'N' ) THEN GO TO 650 END IF ELSE IF ( ICMD(1:7) .EQ. 'AUTOPED' ) THEN CALL AUTOPED WRITE(6,*) $ ' Please wait while the scan averages are re-computed' C C Note that the scan talbes have to be saved first, so that the C allowed range of the phase autocorrelation will be saved. C CALL VSAVE CALL SCANCALC ELSE IF ( ICMD(1:3) .EQ. 'PEX' ) THEN CALL PEXPRT ELSE IF ( ICMD(1:6) .EQ. 'DEBUG ' ) THEN WRITE(6,1600) DO 790 I = 1, ISCAN WRITE(6,1610) I, STAR(I), (SCAL(I,J),J=1,4), $ (VSIG(I,J),J=1,4) 790 CONTINUE 1600 FORMAT ( ' SCAN STAR SCAL1 SCAL2 SCAL3 SCAL4', $ ' VSIG1 VSIG2 VSIG3 VSIG4' ) 1610 FORMAT ( 2I5, 1X, 4F6.3, 1X, 4F6.3 ) ELSE IF ( ICMD(1:7) .EQ. 'EXPORT ' ) THEN C CALL GET_SYS ( NPAR ) WRITE(6,'(A,4F6.3)') ' Calibration errors are ', SYS_CAL WRITE(6,'(A,$)') ' Input values for export' READ (5,*,END=10) CAL_ERR WRITE(6,*) ' Scans will be saved in file ''', FNAME(1:IROOT) $ //'.CAL''' OPEN ( UNIT=33, FILE=FNAME(1:IROOT)//'.CAL', $ STATUS='UNKNOWN' ) OPEN ( UNIT=34, FILE=FNAME(1:IROOT)//'.CL5', $ STATUS='UNKNOWN' ) DO 890 J = 1, NBASE WRITE(33,1900) MOD(YEAR,100), MONTH, DAY, BLIST(J), $ (WAVE(I),I=2,4) WRITE(34,1901) MOD(YEAR,100), MONTH, DAY, BLIST(J), $ (WAVE(I),I=2,4) WRITE( 6,1900) MOD(YEAR,100), MONTH, DAY, BLIST(J), $ (WAVE(I),I=2,4) DO 880 I = 1, ISCAN IF ( OKSCAN(I) .EQ. 15 ) GO TO 880 IF ( BLIST(J).NE. BASE(I) ) GO TO 880 DO 885 ICHAN = 2, 4 TEMP(3*ICHAN-5) = VIS(I,ICHAN) IF ( IAND(OKSCAN(I),MASK(ICHAN)) .EQ. 0 ) THEN TEMP(3*ICHAN-4) = VCAL(I,ICHAN) IF ( ABS(SCAL(I,ICHAN)) .GT. .005 ) THEN C C Calibration error is a fraction of the calibrated vis C TEMP(3*ICHAN-3) = SQRT ( $ ( CAL_ERR(ICHAN)*VCAL(I,ICHAN) )**2 + $ ( VSIG(I, ICHAN)/SCAL(I,ICHAN) )**2 ) ELSE TEMP(3*ICHAN-3) = 0. END IF ELSE TEMP(3*ICHAN-4) = 2.00 TEMP(3*ICHAN-3) = 0.00 END IF 885 CONTINUE Comment WRITE( 6,1950) HOURS(I), STAR(I), TEMP WRITE(33,1950) HOURS(I), STAR(I), TEMP WRITE(34,1951) HOURS(I), STAR(I), TEMP 880 CONTINUE 890 CONTINUE CLOSE(UNIT=33) CLOSE(UNIT=34) 1900 FORMAT ( 3X, 3I2.2, 1H., I5, 3(F7.1,15X) ) 1901 FORMAT ( 3X, 3I2.2, 1H., I5, 3(F7.1,21X) ) 1950 FORMAT ( F10.6, I5, 3(3F7.3,1X) ) 1951 FORMAT ( F10.6, I5, 3(3F9.5,1X) ) ELSE IF ( ICMD(1:5) .EQ. 'SAVE ' ) THEN IF ( FNAME .EQ. ' ' ) THEN WRITE(6,*) ' There is no data to save ' ELSE WRITE(6,*) ' Calibration tables are being saved' CALL VSAVE END IF ELSE IF ( ICMD(1:2) .EQ. 'Q ' ) THEN IF ( FNAME .NE. ' ' ) THEN 895 CONTINUE WRITE(6,*) ' Don''t forget to save the summary file! ' WRITE(6,*) ' ' WRITE(6,*) ' Do you want to save the calibration tables?' READ (5,'(A)') ANS CALL CAPS ( ANS ) IF ( ANS .EQ. 'Y' ) THEN CALL VSAVE ELSE IF ( ANS .NE. 'N' ) THEN GO TO 895 END IF END IF 897 CONTINUE WRITE(6,*) ' Are you sure you want to quit (Y/N)?' READ (5,'(A)') ANS CALL CAPS ( ANS ) IF ( ANS .EQ. 'Y' ) THEN GO TO 950 ELSE IF ( ANS .NE. 'N' ) THEN GO TO 897 END IF ELSE WRITE(6,*) ' INVALID COMMAND, ', ICMD END IF GO TO 10 950 CONTINUE STOP END SUBROUTINE CAPS ( STRING ) CHARACTER *(*) STRING INTEGER *4 I, K, OFF OFF = ICHAR('A') - ICHAR('a') K = LEN(STRING) DO 100 I = 1, K IF ((STRING(I:I).LE.'z') .AND. (STRING(I:I).GE.'a')) $ STRING(I:I) = CHAR( ICHAR(STRING(I:I))+OFF ) 100 CONTINUE RETURN END SUBROUTINE LINFIT ( X0, Y0, YFIT, N, SLOPE, CONST, SIGMA ) C C Fit a streight line through the points X0(I), Y0(I), I = 1, N C The fitted points are in YFIT(I) C SLOPE is the fitted slope C CONST is the fitted y-intersept C SIGMA is the standard dev of residuals IMPLICIT UNDEFINED(A-Z) REAL *4 X0(*), Y0(*), YFIT(*), SLOPE, CONST, SIGMA REAL *4 XX, XY, YY, X, Y, RMS, AVG, RES INTEGER *4 N, I C C Check for illegal input C IF ( N .LE. 0 ) THEN WRITE(6,*) ' NO DATA TO FIT, BOZO ' SLOPE = 0. CONST = 0. SIGMA = 0. GO TO 900 ELSE IF ( N .EQ. 1 ) THEN WRITE(6,*) ' ONE POINT DOESN''T MAKE MUCH OF A LINE ' SLOPE = 0. CONST = Y0(1) YFIT(1) = Y0(1) SIGMA = 0. GO TO 900 END IF C C Zero the accumulators C X = 0. Y = 0. XY = 0. XX = 0. YY = 0. C C Increment the accumulators C DO 50 I = 1, N X = X + X0(I) Y = Y + Y0(I) XX = XX + X0(I)*X0(I) XY = XY + X0(I)*Y0(I) YY = YY + Y0(I)*Y0(I) 50 CONTINUE C C Determine the fit C X = X / FLOAT(N) Y = Y / FLOAT(N) XY = XY / FLOAT(N) XX = XX / FLOAT(N) YY = YY / FLOAT(N) SLOPE = ( XY - X*Y ) / ( XX - X*X ) CONST = Y - SLOPE*X SIGMA = YY -Y*Y + SLOPE * (XY -X*Y) SIGMA = SQRT ( MAX(0., SIGMA) ) C C Evalulate the fitted line and residual delays C RMS = 0. AVG = 0. DO 100 I = 1, N YFIT(I) = CONST + SLOPE*X0(I) RES = Y0(I) - YFIT(I) AVG = AVG + RES RMS = RMS + RES * RES 100 CONTINUE AVG = AVG / FLOAT(N) RMS = SQRT ( RMS / FLOAT(N) - AVG*AVG ) WRITE(6,*) ' LINFIT: RMS = ', RMS, ' EVALUATED AS ', SIGMA SIGMA = RMS 900 CONTINUE RETURN END SUBROUTINE APPLY_CAL ( NPAR ) C C Apply the calibration curves derived in SYSVIS to all of the data. C (including deleted points ) C 30 october 1988 C IMPLICIT UNDEFINED(A-Z) INTEGER *4 I, ICHAN, I0, NPAR(4) REAL *4 RI, TEMP C Global variables. INCLUDE 'VPLOT.INC' C C Divide observed visibilities by C an estimated system visibility C C The calibrations arrays have x-axes defined by C ( For 0 < I < NFIT ) C time = TCAL0 + TCAL1 * I / NFIT C H.A. = HCAL0 + HCAL1 * I / NFIT C Z.D. = ZCAL0 + ZCAL1 * I / NFIT C M.A. = MCAL0 + MCAL1 * I / NFIT C C The other calibrations can be calculated in SYSVIS but cannot be C but cannot be applied to the data. C C Calculate the number of parameters fit to the data. C DO ICHAN = 1, 4 NPAR(ICHAN) = 0 DO I = 1, 4 NPAR(ICHAN) = NPAR(ICHAN) + NFREE(ICHAN,I) END DO END DO DO I = 1, ISCAN DO ICHAN = 1, 4 VCAL(I,ICHAN) = VIS(I,ICHAN) SCAL(I,ICHAN) = 1. END DO END DO DO 400 I = 1, ISCAN C C Time calibration C IF ( HOURS(I) .LT. TCAL0 ) THEN DO ICHAN = 1, 4 VCAL(I,ICHAN) = VCAL(I,ICHAN) / TCAL(1,ICHAN) SCAL(I,ICHAN) = SCAL(I,ICHAN) * TCAL(1,ICHAN) END DO ELSE IF ( HOURS(I) .GT. TCAL0+TCAL1 ) THEN DO ICHAN = 1, 4 VCAL(I,ICHAN) = VCAL(I,ICHAN) / TCAL(NFIT,ICHAN) SCAL(I,ICHAN) = SCAL(I,ICHAN) * TCAL(NFIT,ICHAN) END DO ELSE RI = NFIT * ( HOURS(I)-TCAL0 ) / TCAL1 I0 = INT ( RI ) RI = RI - FLOAT(I0) DO ICHAN = 1, 4 TEMP = (1.-RI)*TCAL(I0,ICHAN) + RI*TCAL(I0+1,ICHAN) VCAL(I,ICHAN) = VCAL(I,ICHAN) / TEMP SCAL(I,ICHAN) = SCAL(I,ICHAN) * TEMP END DO END IF C C Seeing calibration, formerly hour angle calibration. C Calibration can be done versus either laser(1) or phase(2) rms. C IF ( PHI(I,SEEING_TYPE) .LT. HCAL0 ) THEN DO ICHAN = 1, 4 VCAL(I,ICHAN) = VCAL(I,ICHAN) / HCAL(1,ICHAN) SCAL(I,ICHAN) = SCAL(I,ICHAN) * HCAL(1,ICHAN) END DO ELSE IF ( PHI(I,SEEING_TYPE) .GT. HCAL0+HCAL1 ) THEN DO ICHAN = 1, 4 VCAL(I,ICHAN) = VCAL(I,ICHAN) / HCAL(NFIT,ICHAN) SCAL(I,ICHAN) = SCAL(I,ICHAN) * HCAL(NFIT,ICHAN) END DO ELSE RI = NFIT * ( PHI(I,SEEING_TYPE)-HCAL0 ) / HCAL1 I0 = INT ( RI ) RI = RI - FLOAT(I0) DO ICHAN = 1, 4 TEMP = (1.-RI)*HCAL(I0,ICHAN) + RI*HCAL(I0+1,ICHAN) VCAL(I,ICHAN) = VCAL(I,ICHAN) / TEMP SCAL(I,ICHAN) = SCAL(I,ICHAN) * TEMP END DO END IF C C Zenith distance calibration. C IF ( ZD(I) .LT. ZCAL0 ) THEN DO ICHAN = 1, 4 VCAL(I,ICHAN) = VCAL(I,ICHAN) / ZCAL(1,ICHAN) SCAL(I,ICHAN) = SCAL(I,ICHAN) * ZCAL(1,ICHAN) END DO ELSE IF ( ZD(I) .GT. ZCAL0+ZCAL1 ) THEN DO ICHAN = 1, 4 VCAL(I,ICHAN) = VCAL(I,ICHAN) / ZCAL(NFIT,ICHAN) SCAL(I,ICHAN) = SCAL(I,ICHAN) * ZCAL(NFIT,ICHAN) END DO ELSE RI = NFIT * ( ZD(I)-ZCAL0 ) / ZCAL1 I0 = INT ( RI ) RI = RI - FLOAT(I0) DO ICHAN = 1, 4 TEMP = (1.-RI)*ZCAL(I0,ICHAN) + RI*ZCAL(I0+1,ICHAN) VCAL(I,ICHAN) = VCAL(I,ICHAN) / TEMP SCAL(I,ICHAN) = SCAL(I,ICHAN) * TEMP END DO C C Mirror angle calibration. C END IF IF ( MA(I) .LT. MCAL0 ) THEN DO ICHAN = 1, 4 VCAL(I,ICHAN) = VCAL(I,ICHAN) / MCAL(1,ICHAN) SCAL(I,ICHAN) = SCAL(I,ICHAN) * MCAL(1,ICHAN) END DO ELSE IF ( ZD(I) .GT. MCAL0+MCAL1 ) THEN DO ICHAN = 1, 4 VCAL(I,ICHAN) = VCAL(I,ICHAN) / MCAL(NFIT,ICHAN) SCAL(I,ICHAN) = SCAL(I,ICHAN) * MCAL(NFIT,ICHAN) END DO ELSE RI = NFIT * ( MA(I)-MCAL0 ) / MCAL1 I0 = INT ( RI ) RI = RI - FLOAT(I0) DO ICHAN = 1, 4 TEMP = (1.-RI)*MCAL(I0,ICHAN) + RI*MCAL(I0+1,ICHAN) VCAL(I,ICHAN) = VCAL(I,ICHAN) / TEMP SCAL(I,ICHAN) = SCAL(I,ICHAN) * TEMP END DO END IF 400 CONTINUE WRITE(6,*) ' VISIBILITIES CALIBRATED ' RETURN END SUBROUTINE CHECKCAL( CAL, NCAL, RESULT ) C C Check to see if this calibration is set C IMPLICIT UNDEFINED (A-Z) INTEGER *4 NCAL, I REAL *4 CAL(*), YMAX, YMIN, MEAN, DELTA CHARACTER *7 RESULT IF ( NCAL .LT. 1 ) THEN YMAX = 0. YMIN = 0. ELSE YMAX = CAL(1) YMIN = CAL(1) DO I = 1, NCAL YMAX = MAX( CAL(I), YMAX ) YMIN = MIN( CAL(I), YMIN ) END DO END IF DELTA = YMAX - YMIN MEAN = ABS( 0.5*(YMAX+YMIN) - 1. ) IF (( DELTA .GT. .01 ) .OR. (MEAN .GT. .01) ) THEN RESULT = ' SET ' ELSE RESULT = ' ' END IF RETURN END