SUBROUTINE SCANCALC C C Subroutine to calculate the scan averages, using the C calibration tables. C C UNIT 2 point calibration table C 3 scan calibration table C 4 visibility data file C C VERSION 1.0 23 March 1989 C IMPLICIT UNDEFINED (A-Z) SAVE REAL *4 PI, PISQ, RAD PARAMETER ( PI = 3.1415926535 ) PARAMETER (PISQ = PI*PI ) PARAMETER ( RAD = 180. / PI ) CHARACTER *11 VERSN, VERSD PARAMETER ( VERSN = 'VERSION 3.0' ) PARAMETER ( VERSD = '14 OCT 1988' ) INTEGER *2 I2SCAN, I2OK INTEGER *4 NSUM, IERR, N, I, J, K, L, VCNT, NPTS, ILAST, IOK INTEGER *4 NDATA(0:4), NDARK, NAME(4), SPECTRUM, ILEN, IFIRST INTEGER *4 MASK(4), IREC, ICOM REAL *4 INTTIME, BVAR(4), DEN, VAMP(4), VSQ, PAR, RV, PMRA, PMDEC REAL *4 COSMA, TEMP(16), WAVES(0:45) REAL *4 NEWMAG, NEWBV REAL *8 RAM, DECM, TEMP8 INCLUDE 'VPLOT.INC' CHARACTER * 1 ANS CHARACTER * 3 MNTH(12) CHARACTER * 7 FILTERS(0:45) CHARACTER * 8 ICMD LOGICAL NOPCAL, NOSCAL, INSCAN 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), IFILT(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 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), IFILT(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 MNTH / 'JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', $ 'JUL', 'AUG', 'SEP', 'OCT', 'NOV', 'DEC' / 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/01 ', '410/01 ', $ '435/01 ', '486/01 ', '656/01 ', $ '650/01 ', '750/999', '656/10 ', $ '037/0 ', '038/0 ', '039/0 ', $ '040/0 ', '623/10 ', '670/05 ', $ '701/05 ', '712/10 ', '754/05 ' / DATA WAVES / 700.0, $ 399.7, 400. , 400. , $ 451.8, 450.2, 450. , $ 501.2, 500. , 500. , $ 550.1, 550. , 550. , $ 600.9, 599.5, 600. , $ 650.6, 650. , 650. , $ 699.0, 700. , 700. , $ 750.1, 750. , 750. , $ 799.2, 800.0, 800. , $ 687. , 397. , 410. , $ 434. , 486. , 656. , $ 650. , 750. , 656. , $ 0. , 0. , 0. , $ 0. , 623. , 672. , $ 701. , 712. , 754. / DATA MASK / 1, 2, 4, 8 / C----------------------------------------------------------------------- C Read in data, using the calibration tables. C The visibility data file and the calibration files are C opened by VREAD and are not closed until a new file is input. C C The first record is a header record. C----------------------------------------------------------------------- READ (3,REC=1,IOSTAT=IERR) I2SCAN, I2OK, YEAR, MONTH, DAY, $ (FILT(J),J=1,4), FORMATID, (NDFRATIO(J),J=1,4), $ COHINT, (TEMP(J),J=1,2) IF ( IERR .NE. 0 ) THEN WRITE(6,*) ' CANNOT READ FIRST SCAN OF SCAN CALIBRATION ', $ 'TABLE' WRITE(6,*) ' RECALC goes down in flames ' STOP END IF IF ( COHINT .EQ. 0 ) COHINT = 4 WRITE(6,'(A,I5)') ' COHERENT INTEGRATION TIME = ', COHINT IF ( FORMATID .NE. 1 ) THEN WRITE(6,*) ' ***ERROR****** BAD FORMAT ID TYPE ', FORMATID WRITE(6,*) ' TYPE ''C'' TO CONTINUE ANYWAY, ' WRITE(6,*) ' ANYTHING ELSE TO ABORT NOW ' READ (5,'(A)') ANS CALL CAPS ( ANS ) IF ( ANS .NE. 'C' ) STOP FORMATID = 1 END IF ISCAN = I2SCAN WRITE(6,*) ISCAN, ' scans recorded in scan calibration table. ' IF ( (ISCAN.LE.0).OR.(YEAR.LE.0).OR.(DAY.LE.0) ) THEN WRITE(6,*) ' Corrupt header on scan calibration table' WRITE(6,*) ' Cannot load this data set ' GO TO 950 END IF WRITE(DATE,1100,IOSTAT=IERR) DAY, MNTH(MONTH), YEAR DO 135, J = 1, 4 FILTERID(J) = FILTERS(FILT(J)) WAVE(J) = WAVES (FILT(J)) 135 CONTINUE WRITE(6,*) ' FILTER NUMBERS ARE ', FILT WRITE(6,*) ' AT WAVELENGTHS ', WAVE WRITE(6,*) ' Data is for ', DATE 1100 FORMAT ( I2, 1X, A3, 1X, I4 ) IF ( IERR .NE. 0 ) THEN WRITE(6,*) ' Date in scan calibration table is corrupt' WRITE(6,*) ' Cannot load this data set' GO TO 950 END IF C----------------------------------------------------------------------- C Read the point calibration prescription from the scan calibration C table. C READ (3,REC=ISCAN+3,IOSTAT=IERR) (DVDP(J),J=1,4), PHIMIN, PHIMAX, $ PHIREF, PHIPEAK, (TEMP(J),J=1,7) IF ( IERR .NE. 0 ) THEN WRITE(6,*) ' CANNOT READ POINT CALIBRATION PARAMETERS FROM', $ ' SCAN CALIBRATION TABLE' WRITE(6,*) ' You are probably using an obsolete cal table. ' WRITE(6,*) ' Read attempt aborted ' GO TO 950 END IF IF ( PHIMAX .EQ. 0. ) THEN PHIMIN = 0. PHIMAX = 5. END IF WRITE(6,'(A,3F8.2)') ' PHIMIN, PHIMAX, PHIREF = ', $ PHIMIN, PHIMAX, PHIREF WRITE(6,'(A,4F8.3)') ' DVDP = ',DVDP WRITE(6,'(A,2F8.3)') ' PHIPEAK, LASER PEAK = ', PHIPEAK C----------------------------------------------------------------------- C Read the data file, according to the prescription in the C Scan calibration table. C The header/time record for scan I is in record REC0(I) C The data is in records REC0(I)+1 to REC0(I) C The dark count data is in DREC(I). if DREC(I) = 0, then no C dark count data was recorded C----------------------------------------------------------------------- WRITE(6,1110) DO 190 I = 1, ISCAN READ(3,REC=I+1,IOSTAT=IERR) I2SCAN, I2OK, $ REC0(I), REC1(I), DREC(I), (D(I,K),K=1,4), $ (SCAL(I,K),K=1,4), NDFSTATUS(I), (TEMP(K),K=1,3) J = I2SCAN OKSCAN(I) = I2OK IF ( IERR .NE. 0 ) THEN WRITE(6,*) ' Error ', IERR, ' reading scan calibration ', $ ' table record ', I ISCAN = I - 1 GO TO 950 ELSE IF ( J .NE. I ) THEN WRITE(6,*) ' Error in scan calibration table format ' WRITE(6,*) ' Scan ', I, ' is labeled ', J ISCAN = I - 1 GO TO 950 ELSE IF ( REC0(I) .EQ. REC1(I) ) THEN WRITE(6,*) ' ERROR scan ', I, ' has no data ' GO TO 950 END IF C----------------------------------------------------------------------- C Record the data stored in the scan and dark count records. C The dark counts themselves have already been read and are C stored in the scan calibration table. C----------------------------------------------------------------------- READ ( 4, REC=REC0(I), IOSTAT=IERR) IA IF ( IERR .NE. 0 ) THEN WRITE(6,*) ' Unexpected error ', IERR, $ ' encountered while reading data file ' WRITE(6,*) ' BOZO CITY ' STOP END IF INTTIME = FLOAT(TAU0) / 1000. STAR(I) = ISTAR BASE(I) = IBASE IF ( DREC(I) .NE. 0 ) THEN READ ( 4, REC=DREC(I), IOSTAT=IERR) IA DELAY(I) = 1.E-6 * DLAS ZD(I) = ZD0 HA(I) = HA0 / 15. DEC(I) = DEC0 COSMA = SIN(HA0/RAD) * COS(DEC0/RAD) IF ( BASE(I) .GT. 2 ) COSMA = - COSMA MA(I) = RAD * ACOS( COSMA ) DJITTER(I) = DJTR DUTYC(I) = TLAST - TFIRST IF ( DUTYC(I) .GT. 0. ) THEN NLOCK(I) = FLOAT(ILOCK) DUTYC(I) = INTTIME * FLOAT(REC1(I)-REC0(I)) ELSE NLOCK(I) = 0. DUTYC(I) = 0. END IF ELSE DELAY(I) = 0.0 ZD(I) = 0. HA(I) = 0. DEC(I) = 0. MA(I) = 0. DJITTER(I) = 0. DUTYC(I) = 0. NLOCK(I) = 0. DUTYC(I) = 0. END IF C----------------------------------------------------------------------- C Average the data stored in the data records C C NOTE: if any point data is deleted on a channel, that data is not C included in the average visibility for that channel. However, C the point has to be deleted on all four channels to be not C included in the average time and jitter calculations. That is, C the weak dependence time and jitter have on channel number are C ignored. C----------------------------------------------------------------------- HOURS (I) = 0. TJITTER(I) = 0. PHI(I,1) = 0. PHI(I,2) = 0. NDATA(0) = 0 DO J = 1, 4 NDATA(J) = 0 BDEN(I,J) = 0. BRMS(I,J) = 0. BVAR(J) = 0. VAMP(J) = 0. VSIG(I,J) = 0. END DO DO 160 K = REC0(I)+1, REC1(I) READ ( 4, REC=K, IOSTAT=IERR) IA IF ( IERR .NE. 0 ) THEN WRITE(6,*) ' Unexpected error ', IERR, $ ' encountered while reading data file ' WRITE(6,*) ' BOZO CITY ' STOP END IF C C Read point cal table to determine if data is usable. C READ ( 2, REC=K, IOSTAT=IERR) IOK IF ( IERR .NE. 0 ) THEN WRITE(6,*) ' Unexpected error ', IERR, $ ' encountered reading point cal file' END IF IF ( IOK .EQ. 15 ) THEN GO TO 160 END IF IF ( (PHI12.GT.PHIMAX) .OR. (PHI12.LT.PHIMIN) ) THEN GO TO 160 END IF NDATA(0) = NDATA(0) + 1 HOURS(I) = HOURS(I) + 24.*TIME0 TJITTER(I) = TJITTER(I) + JITTR IF ( JITTR .LT. 0 ) TJITTER(I) = TJITTER(I) + 256 PHI(I,1) = PHI(I,1) + PHIRMS PHI(I,2) = PHI(I,2) + PHI12 DO J = 1, 4 IF ( IAND(IOK,MASK(J)) .EQ. 0 ) THEN NDATA(J) = NDATA(J) + 1 BDEN(I,J) = BDEN(I,J) + ADEN(J) BRMS(I,J) = BRMS(I,J) + ADEN(J)*ADEN(J) BVAR(J) = BVAR(J) + AVAR(J) VAMP(J) = VAMP(J) + ANUM(J) VSIG(I,J) = VSIG(I,J) + ANUM(J)*ANUM(J) END IF END DO 160 CONTINUE IF ( NDATA(0) .EQ. 0 ) THEN WRITE(6,1115) I 1115 FORMAT ( '+ NO DATA FOR SCAN ', I5, / ) OKSCAN(I) = 15 HOURS(I) = 0. TJITTER(I) = 0. PHI(I,1) = 0. PHI(I,2) = 0. DO J = 1, 4 BDEN(I,J) = 0. BRMS(I,J) = 0. VIS(I,J) = 0. VCAL(I,J) = 1. VSIG(I,J) = 0. SCIN(I,J) = 0. END DO ELSE HOURS(I) = HOURS(I) / FLOAT(NDATA(0)) TJITTER(I) = TJITTER(I) / FLOAT(NDATA(0)) PHI(I,1) = PHI(I,1) / FLOAT(NDATA(0)) PHI(I,2) = PHI(I,2) / FLOAT(NDATA(0)) DO 170 J = 1, 4 IF ( NDATA(J) .LE. 1 ) THEN OKSCAN(I) = IOR( OKSCAN(I), MASK(J) ) BRMS(I,J) = 0. VIS(I,J) = 2.0 VSIG(I,J) = 0.0 ELSE BDEN(I,J) = BDEN(I,J) / FLOAT(NDATA(J)) BRMS(I,J) = BRMS(I,J)/FLOAT(NDATA(J)) - $ BDEN(I,J)*BDEN(I,J) BRMS(I,J) = SQRT( BRMS(I,J)/FLOAT(NDATA(J)-1) ) DEN = BDEN(I,J) - D(I,J) BVAR(J) = BVAR(J) / FLOAT(NDATA(J)) VAMP(J) = VAMP(J) / FLOAT(NDATA(J)) VSIG(I,J) = VSIG(I,J) / FLOAT(NDATA(J)) IF ( NDATA(J) .GT. 1 ) THEN VSIG(I,J) = SQRT( (VSIG(I,J)-VAMP(J)*VAMP(J)) $ / FLOAT(NDATA(J)-1) ) ELSE VSIG(I,J) = 0. END IF C C N.B. Throughout this entire program, what is called visibility C is, in reality, visibility squared. C C The reason for going to visibility squared everywhere, instead of C visibility, is so that the error estimate is correct. This is only C a problem when the visibilty is less than or about equal to the C estimated error, but these are very important data points. C C Note: the photon noise estimate, SIG_PHOT, is correct only for C 4 ms incoherent, half second coherent averages. C IF ( (DEN.GT..01).AND.(NDATA(J).GT.1) ) THEN VIS (I,J) = PISQ*VAMP(J) /(2.*DEN*DEN) VSIG(I,J) = PISQ*VSIG(I,J)/(2.*DEN*DEN) SIG_PHOT(I,J) = PISQ*PISQ/(500.*NDATA(J)) * $ (1./(DEN*DEN) + 4.*VIS(I,J)/(PISQ*DEN)) SIG_PHOT(I,J) = SQRT ( SIG_PHOT(I,J) ) ELSE WRITE(6,1050) BDEN(I,J), D(I,J) VIS (I,J) = 5. VSIG(I,J) = 0. OKSCAN(I) = IOR ( OKSCAN(I),MASK(J) ) SIG_PHOT(I,J) = 0. END IF 1050 FORMAT ( /,' bad ch, flux =', f10.2, ' dark = ', f10.2,/) C----------------------------------------------------------------------- C Correct the data for the correlation between the visibility C and the phase fluctuations. C----------------------------------------------------------------------- VIS(I,J) = VIS(I,J) / $ (1.+DVDP(J)*( PHI(I,2)-PHIREF)) IF (SCAL(I,J) .NE. 0. ) THEN VCAL(I,J) = VIS (I,J) / SCAL(I,J) ELSE VCAL(I,J) = VIS(I,J) END IF END IF 170 CONTINUE END IF WRITE(6,1120) I, STAR(I), BASE(I), HOURS(I), REC0(I), $ REC1(I), DREC(I), OKSCAN(I) 1110 FORMAT ( ' SCAN STAR B TIME REC0 REC1 DREC BAD ',/) 1120 FORMAT ( '+', 2I5, I3, F8.3, 2X, 4I6, $ ) 190 CONTINUE SCAN0 = 1 SCAN1 = ISCAN C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Read the calibration arrays. C----------------------------------------------------------------------- DO I = 1, NFIT READ (3,REC=I+5+ISCAN, IOSTAT=IERR) $ (TCAL(I,K),K=1,4), (HCAL(I,K),K=1,4), $ (ZCAL(I,K),K=1,4), (MCAL(I,K),K=1,4) IF ( IERR .NE. 0 ) THEN WRITE(6,*) ' FATAL ERROR ', IERR, ' READING SCAN ', $ ' CALIBRATION TABLE ' END IF END DO C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Read the auxiliary information records. C----------------------------------------------------------------------- IREC = ISCAN + NFIT + 5 DO I = 1, 4 D_ORDER(1,I) = 0 D_ORDER(2,I) = 0 D_RANGE(1,I) = 0. D_RANGE(2,I) = 0. D_RMS(I) = 0. D_CHISQ(I) = 0. NCAL_RMS(I) = 0. NCALSCAN(I) = 0 MAG0(I) = 0. CTERM(I) = 0. END DO ICOM = 0 IDUSK = 0 IDAWN = 0 NCALSTAR = 0 DO I = 1, 200 CALSTAR(I) = 0 END DO 200 CONTINUE IREC = IREC + 1 READ (3,REC=IREC, IOSTAT=IERR) TEMP8, (TEMP(I),I=3,16) IF ( IERR .NE. 0 ) GO TO 950 WRITE( ICMD, '(A8)', IOSTAT=IERR ) TEMP8 C----------------------------------------------------------------------- IF ( ICMD .EQ. 'D_ORDER ' ) THEN READ(3,REC=IREC, IOSTAT=IERR) TEMP(1), TEMP(2), $ D_ORDER, (TEMP(I),I=11,16) WRITE(6,'(2X,A,3X,8I5)') ICMD, D_ORDER ELSE IF ( ICMD .EQ. 'D_RANGE ' ) THEN READ(3,REC=IREC, IOSTAT=IERR) TEMP(1), TEMP(2), $ D_RANGE, IDUSK, IDAWN, (TEMP(I),I=11,16) WRITE(6,'(2X,A,3X,8F7.2)') ICMD, D_RANGE ELSE IF ( ICMD .EQ. 'D_STAT ' ) THEN READ(3,REC=IREC, IOSTAT=IERR) TEMP(1), TEMP(2), $ D_RMS, D_CHISQ, (TEMP(I),I=11,16) WRITE(6,'(2X,A,3X,8F7.2)') ICMD, D_RMS, D_CHISQ ELSE IF ( ICMD .EQ. 'NCAL_FIT' ) THEN READ(3,REC=IREC, IOSTAT=IERR) TEMP(1), TEMP(2), $ MAG0, CTERM, (TEMP(I),I=11,16) WRITE(6,'(2X,A,3X,8F7.2)') ICMD, NCAL_RMS ELSE IF ( ICMD .EQ. 'NCAL_RMS' ) THEN READ(3,REC=IREC, IOSTAT=IERR) TEMP(1), TEMP(2), $ NCAL_RMS, (TEMP(I),I=7,16) WRITE(6,'(2X,A,3X,8F7.2)') ICMD, NCAL_RMS ELSE IF ( ICMD .EQ. 'NCALSTAR' ) THEN READ(3,REC=IREC, IOSTAT=IERR) TEMP(1), TEMP(2), $ NCALSTAR, (CALSTAR(I),I=1,13) K = 13 DO WHILE ( K .LT. NCALSTAR ) IREC = IREC + 1 READ(3,REC=IREC, IOSTAT=IERR) (CALSTAR(I),I=K+1,K+16) K = K + 16 END DO WRITE(6,'(2X,A,3X,10(4I5,2X,4I5,/,13X)') ICMD, $ (CALSTAR(I),I=1,NCALSTAR) ELSE IF ( ICMD .EQ. 'NCALSCAN' ) THEN READ(3,REC=IREC, IOSTAT=IERR) TEMP(1), TEMP(2), $ NCALSCAN, (TEMP(I),I=7,16) WRITE(6,'(2X,A,3X,8F7.2)') ICMD, NCALSCAN ELSE IF ( ICMD .EQ. 'COMMENT ' ) THEN ICOM = ICOM + 1 IF ( ICOM .GT. 5 ) THEN WRITE(6,*) ' ERROR TOO MANY COMMENTS RECORDED ' ELSE READ (3,REC=IREC, IOSTAT=IERR) (TEMP(I),I=1,16) WRITE(COMMENT(ICOM), '(14A4)') (TEMP(I),I=3,16) END IF ELSE IF ( ICMD .EQ. 'END ' ) THEN WRITE(6,'(2X,A)') ICMD GO TO 950 ELSE WRITE(6,'(2X,2A)') ICMD, ' UNKNOWN COMMAND ' END IF GO TO 200 950 CONTINUE RETURN C----------------------------------------------------------------------- 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 ) 2020 FORMAT( I7, ' OBSERVED STARS' ) 2021 FORMAT( ' THE', I3, ' OBSERVED BASELINES ARE' ) 2025 FORMAT( 5I10 ) 2050 FORMAT ( 2I5, 2X, 2A8 ) END