Program TWOC C Version 2.03 C C by David Mozurkewich 9 Sep 1991 C C Does two color astrometry of the Mt Wilson Interferometer data. C It outputs one and two color astrometry. C IMPLICIT UNDEFINED (A-Z) REAL *8 PI, RAD, LON, LAT PARAMETER ( PI = 3.14159265358979D0 ) PARAMETER ( RAD = 180.D0 / PI ) PARAMETER ( LON = 118.059167D0 ) PARAMETER ( LAT = 34.2166944D0) INTEGER *4 MAXSCAN, nmax PARAMETER ( MAXSCAN = 4000 ) ! max number of points per scan CHARACTER *12 VERSN, VERSD PARAMETER( VERSN = 'VERSION 2.03' ) PARAMETER( VERSD = '1991 SEP 09' ) parameter (nmax=4000) C....................................................................... INTEGER *4 ISCAN, IPOINT, ICHAN, NPTS, N2C, n1ce, imark(nmax) INTEGER *4 I, J, NBODY, TDTOFF, NCOH, NSCAN INTEGER *4 STRTYR, STRTMON, STRTDAY, SPECT INTEGER *4 STRTSTAR, STRTSEC, STOPSEC, TEMP REAL *4 TCOH, TFIRST, TLAST, D1C, d1ce, NAVG(4), D2C(2), D2RMS(2) REAL *4 DRES(MAXSCAN), VSQ(MAXSCAN), PHASE(4,MAXSCAN) REAL *4 LAMBDA(4), FRINGE_OFFSET(3),dispf(4) REAL *8 DJUL, UT1, GAST, TJD, DRMS, drmse, VRAW(4) REAL *8 PDRES, PN(4), PNSQ(4), TAVG, RAA, DECA, RAM, DECM, MAG REAL *8 COHINT, ETIME, STIME, SCANTIME, DELAYC, VSUM(4), VRMS(4) REAL *8 PHIRMS, PHIAC, DELAY, LASER0(MAXSCAN), TIME(MAXSCAN) CHARACTER *1 ANS CHARACTER *3 MNTH(12) character scanid*11 logical noplot, hsearch C....................................................................... INCLUDE 'AVER.INC' C....................................................................... REAL *8 JD CHARACTER *8 HMS CHARACTER *15 STATUS EXTERNAL JD DATA MNTH / 'JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', $ 'JUL', 'AUG', 'SEP', 'OCT', 'NOV', 'DEC' / C----------------------------------------------------------------------- C Print out headers, initialize unit numbers. C Read the default baseline list. C CALL OPSTART ( versn, versd ) C----------------------------------------------------------------------- C Get the input parameters, open input and output files C STRTSTAR is the first star to look at. C STRTSEC and STOPSEC are time limits for the data C in UT seconds from start of the day. C After all the data is read in, control returns here to C process another file, if desired. 10 CONTINUE CALL OPINPUT( STATUS, STRTYR, STRTMON, STRTDAY, STRTSTAR, $ STRTSEC,STOPSEC,COHINT,SCANTIME,FRINGE_OFFSET,dispf ) IF (STATUS .EQ. 'QUIT' ) GO TO 900 noplot=.false. if(strtstar.lt.0)then noplot=.true. strtstar=abs(strtstar) if(strtstar.eq.9999)strtstar=0 endif STIME = DBLE (STRTSEC) ETIME = DBLE (STOPSEC) TDTOFF = 55 NCOH = 1000. * COHINT / TSAMPLE + 0.5 COHINT = 0.001 * TSAMPLE * NCOH NSCAN = SCANTIME / COHINT + 0.5D0 NSCAN = MIN ( NSCAN, MAXSCAN ) SCANTIME = COHINT * DBLE(NSCAN) WRITE(TERMM,*) ' Input file name = ', INDSN IF ( STRTMON .EQ. 0 ) THEN WRITE(TERMM,1020) ' Date = any ' ELSE WRITE(TERMM,1020) ' date = ', $ STRTYR, MNTH(STRTMON), STRTDAY END IF IF ( STRTSTAR .EQ. 0 ) THEN WRITE(TERMM,1030) ' Star = any ' ELSE WRITE(TERMM,1030) ' Star = ', STRTSTAR END IF WRITE(TERMM,*) ' Data averaged from ', HMS(STRTSEC), ' to ', $ HMS(STOPSEC) WRITE(TERMM,1040) ' Coherent int time = ', COHINT, NCOH WRITE(TERMM,1041) ' Maximum scan length = ', SCANTIME, NSCAN WRITE(6,*) ' Fringe offsets = ', FRINGE_OFFSET, ' microns' WRITE(TERMM,*) ' Write scan averages to output file' WRITE(6,*) ' Data offset from start of file = ', OFFSET C----------------------------------------------------------------------- C There are many bugs in the tape format. The filter identifications C are stored in the header records for each scan. But the first scan C on a tape often does not have a header record. The wavelengths of C the filters are needed to do the averaging. C C Therefore, we skip through the data until a header record is C encountered. We look up the filter wavelengths and rewind the file. C----------------------------------------------------------------------- hsearch=.true. DO I = 1, 4 FILTER(I) = -1 END DO WRITE(6,*) 'Search for filter wavelengths... ' 40 CONTINUE CALL GETADATA ( NCOH, TIME(1), DRES(1), PHASE(1,1), VSQ(1), $ PN, LASER0(1), PHIRMS, PHIAC, STATUS, hsearch ) IF ( STATUS .EQ. 'NO MORE DATA' ) THEN WRITE(6,*) ' NO HEADER RECORDS ON THIS TAPE !! ' GO TO 600 ELSE IF ( ( STATUS .EQ. 'END OF STAR' ) .OR. $ ( STATUS .EQ. 'NO FILTER ID' ) .OR. $ ( STATUS .EQ. 'BAD FILTER ID' ) ) THEN GO TO 40 END IF WRITE(6,'(A,4F8.1)') ' FILTER WAVELENGTHS ARE = ', WAVE DO I = 1, 4 LAMBDA(I) = WAVE(I) END DO C----------------------------------------------------------------------- WRITE(TERMP,1480) IDAY, MNTH(IMON), IYEAR C----------------------------------------------------------------------- C Skip to the first data to be processed. C hsearch=.false. IREC = 0 if(offset.ne.0)then status=' ' goto 4321 endif STATUS = 'FIRST' C WRITE(6,*) ' SEARCHING FOR FIRST REQUESTED SCAN. ' DO WHILE ( STATUS .NE. ' ' ) CALL GETADATA ( NCOH, TIME(1), DRES(1), PHASE(1,1), VSQ(1), $ PN, LASER0(1), PHIRMS, PHIAC, STATUS, hsearch ) c write(6,'(a,a)')' Search: status = ',status IF ( STATUS .EQ. 'NO MORE DATA' ) GO TO 600 IF ( ( TIME(1) .LT. STIME ) .OR. $ ( (ISTAR.NE.STRTSTAR).AND.(STRTSTAR.NE.0) ) ) THEN STATUS = 'EARLY' END IF END DO IF ((STRTYR .NE. IYEAR).AND.(STRTYR.NE.0)) THEN WRITE(TERME,1050) 'year', STRTYR, 'year', IYEAR GO TO 600 ELSE IF ((STRTMON.NE.IMON).AND.(STRTMON.NE.0)) THEN WRITE(TERME,1050) 'month', STRTMON, 'month', IMON GO TO 600 ELSE IF ((STRTDAY.NE.IDAY).AND.(STRTDAY.NE.0)) THEN WRITE(TERME,1050) 'day', STRTDAY, 'day', IDAY GO TO 600 END IF C----------------------------------------------------------------------- C The scans start here. C One SCAN is consecutive data on one star/baseline pair C with a maximum length of NSCAN coherent averages. C 4321 DO ISCAN = 1, 2048 IF ( TIME(1) .GT. ETIME ) THEN WRITE(6,*) ' Last time processed ' GO TO 600 END IF c write(6,*)' New scan: IREC + OFFSET = ',irec+offset c write(6,'(a,a)')' Status= ', status C----------------------------------------------------------------------- C Zero statistics variables for use by RDATA. C NLOCK = 1 NDATA = 0 DO I = 1, NFILT BDATA(I) = 0 END DO C----------------------------------------------------------------------- C Read the data for this scan. C NPTS = 0 DO ICHAN = 1, 4 NAVG(ICHAN) = 0 END DO IF ( STATUS .EQ. 'NO MORE DATA' ) GO TO 600 WRITE(6,'(/A,i5)') ' SCAN ', ISCAN DO I = 1, NSCAN CALL GETADATA ( NCOH, TIME(I), DRES(I), PHASE(1,I), $ VSQ(I), PN, LASER0(I), PHIRMS, PHIAC, STATUS, hsearch ) WRITE(6,'(A,I5,2A)') $ '+INTEGRATION ', NPTS, ' STATUS = ', STATUS IF ( STATUS .NE. ' ' ) THEN GO TO 120 END IF DO ICHAN = 1, 4 NAVG(ICHAN) = NAVG(ICHAN) + PN(ICHAN) END DO NPTS = I END DO 120 CONTINUE if ( npts .lt. 50 ) then WRITE(TERMM,1590) ISTAR, IBASE, NPTS WRITE(TERMP,1590) ISTAR, IBASE, NPTS WRITE(TERM2,1590) ISTAR, IBASE, NPTS GO TO 500 endif IF ( NPTS .EQ. 0 ) THEN WRITE(TERMM,1570) ISTAR, IBASE WRITE(TERMP,1570) ISTAR, IBASE WRITE(TERM2,1570) ISTAR, IBASE GO TO 500 END IF DO ICHAN = 1, 4 NAVG(ICHAN) = NAVG(ICHAN) / FLOAT(NPTS) END DO C----------------------------------------------------------------------- C Make the two color plots. C scanid='GOOD ID0000' 90 CALL CFI (DRES, PHASE, VSQ, MAXSCAN, NPTS, FRINGE_OFFSET, $ I,ISTAR,IBASE,LAMBDA,D2C,D2RMS,N2C,TERMP,imark,noplot,scanid, $ dispf) C----------------------------------------------------------------------- C Determine the averge time and 1 color residual delay. C TAVG = 0. D1C = 0. d1ce = 0. DRMS = 0. drmse= 0. n1ce = 0. DO I = 1, NPTS TAVG = TAVG + TIME(I) - TIME(1) D1C = D1C + DRES(I) DRMS = DRMS + DRES(I)**2 if(imark(i).eq.3)then d1ce = d1ce + dres(i) drmse= drmse + dres(i)**2 n1ce = n1ce + 1 endif END DO TAVG = TIME(1) + TAVG / FLOAT(NPTS) D1C = D1C / FLOAT(NPTS) if(n1ce.ne.0)then d1ce = d1ce / float(n1ce) drmse= sqrt ( drmse/float(n1ce) - d1ce**2 ) else d1ce=0 endif DRMS = SQRT ( DRMS/FLOAT(NPTS) - D1C**2 ) WEIGHT= NPTS C----------------------------------------------------------------------- C Calculate the delay for the mean time C DJUL = JD ( IDAY, IMON, IYEAR ) UT1 = TAVG / 86400.D0 CALL SIDTIM ( DJUL, UT1, 1, GAST ) TJD = DJUL + UT1 + DBLE(TDTOFF)/86400.D0 NBODY = 3 CALL APSTAR ( TJD, NBODY, RAM, DECM, PMRA, PMDEC, PAR, $ RV, RAA, DECA ) CALL DIURN ( TJD, GAST, RAA, DECA, RA, DEC ) HA = 15.D0 *( GAST - RA ) - LON IF ( HA .LT. -180.D0 ) HA = HA + 360.D0 DELAYC = CON + BX*COS(DEC/RAD)*COS(HA/RAD) $ - BY*COS(DEC/RAD)*SIN(HA/RAD) $ + BZ*SIN(DEC/RAD) ZA = RAD * ACOS( SIN(DEC/RAD)*SIN(LAT/RAD) + $ COS(HA/RAD)*COS(DEC/RAD)*COS(LAT/RAD) ) DELAY = DELAYC + D1C DO J = 1, 2 D2C(J) = D2C(J) END DO WRITE(6,*) ' DATE = ', IYEAR, IMON, IDAY WRITE(6,*) ' TIME = ', 24.D0*UT1 WRITE(6,*) ' SCAN, STAR, BASELINE ', ISCAN, ISTAR, IBASE WRITE(6,*) ' DELAY = ', DELAY WRITE(6,*) ' RES DELAYS 1-C, 2-C ', D1C, D2C WRITE(6,*) ' RES DELRMS 1-C, 2-C ', DRMS, D2RMS WRITE(6,*) ' NPTS, N2C = ', NPTS, N2C c write .two file here WRITE(TERM2,1580) MOD(IYEAR,100), IMON, IDAY, 24.*UT1, ISTAR, $ IBASE,DELAY,D1C,d1ce,D2C,DRMS,drmse,D2RMS,NPTS,N2C,scanid TEMP = INT ( 8.64D4*UT1 ) c write .prt file here WRITE(TERMP,1550) ISCAN, ISTAR, IBASE, NAME, $ HMS(TEMP), NLOCK, (4*NDATA)/NLOCK, $ DRES(1), DRMS, BDATA, VRAW(1), INT(NAVG(1)) c++++++++++++++++++++++++ c If the delay difference is off by more than 2 micron from c fringe_offset, we change the id! (call CFI one more time) c if(noplot.and.scanid(8:8).eq.'0')then if(abs(d2c(1)-d2c(2)-fringe_offset(ibase)).gt.2.0)then write(12,'(a,a)')' *** Delay offset too large, ', $ 'try other peak!' write(6,'(a,a)')' *** Delay offset too large, ', $ 'try other peak!' scanid(8:8)='1' scanid(10:11)='00' goto 90 endif endif c+++++++++++++++++++++++ c write .avg file here CALL WRITREC ( TERMO, UT1, DELAY, DRMS, ISCAN ) IF ( STATUS .EQ. 'NO MORE DATA' ) GO TO 600 500 CONTINUE END DO ! end of loop over scans (ISCAN). C----------------------------------------------------------------------- WRITE(TERMP,*) ' PROGRAM ENDS BECAUSE OF TOO MANY SCANS!!!!' 600 CONTINUE CALL OPCLOSE GO TO 10 900 CONTINUE STOP C----------------------------------------------------------------------- 1010 FORMAT ( A, 2F11.0 ) 1020 FORMAT ( A, I5, A4, I5 ) 1030 FORMAT ( A, I5 ) 1040 FORMAT ( A, F8.4, ' sec ', I6, ' points ' ) 1041 FORMAT ( A, F8.4, ' sec ', I6, ' coherent averages' ) 1045 FORMAT ( A, 4F8.6 ) 1050 FORMAT(' starting ', A, I5, ' is not equal to current ', A, I5) 1400 FORMAT ( I5, I2, F15.3, 4F8.1, 4F7.3 ) 1480 FORMAT ( ' DATE = ', I5, A4, I5 / $ ' NPTS Star B Residual Delay ', $ ' Std Deviation ', $ 'Fringe visibility squared points per fringe ', $ ' VED NPTS N2C MODE' / $ ' No. 1 C 2 C 2 C ', $ ' 1 C 2 C 2 C ', $ ' -2 -1 0 1 2 -2 -1 0 1 2 ', $ ' '/ $'------------------------------------------------------------', $'--------------------------------------------------') 1550 FORMAT ( I5, I5, I3, 1X, A8, 1X, A8, ' NLOCK=',I5, $ ' 4*NDATA/NLOCK=', I6,' DRES=', F10.1,' DRMS=', F7.1, $ ' BDATA=', I2,3I3,' VRAW(1)=', F5.2, 1X,' NAVG=', I5 ) 1570 FORMAT ( 4X, I5, I3, 5X, ' NO DATA ' ) 1580 FORMAT ( I3, '/', I2.2, '/', I2.2, F12.6, I5, I3, F14.4, 4F10.3, $ 4F7.2, 2I5 ,1x,a11) 1590 FORMAT ( 4X, I5, I3, 5X, ' TOO LITTLE DATA, NPTS = ',I5 ) END