Program DELAYS C C Program to calculate geometric delay for running the NPOI. C Input is the FK% catalog star number, and the current date C and time. C C Output is a list of delays and delay rates for C every 5 minutes for the next hour. C C C Note: This program and Buscher use different sign conventions C for longitudes. C C Modified to do delays for only one star, but for the full night. C C IMPLICIT UNDEFINED (A-Z) Real *8 b1, b2, b3, con, bx, by, bz Real *8 ra, ha, dec, za, delay, delay_rate Real fdl1(2), fdl2(2), fdl1_dot(2) integer idate, iyear, imon, iday, hour, minute, utsec(2), istar integer starlist(5), icol, ihour REAL *8 PI, RAD, longitude, latitude PARAMETER ( PI = 3.14159265358979D0 ) PARAMETER ( RAD = 180.D0 / PI ) PARAMETER ( longitude= 111. + 32./60. + 6.0/3600. ) PARAMETER ( latitude = 35. + 5./60. + 48.0/3600. ) INTEGER *4 NBODY, TDTOFF, SPECT REAL *4 PMRA4, PMDEC4, PAR4, RV4, MAG4, BV INTEGER *4 I REAL *4 ANGDIA REAL *8 DJUL, UT1, GAST, TJD, sideral REAL *8 RAM, DECM, PMRA, PMDEC, PAR, RV REAL *8 appRA(5), appDEC(5) REAL *8 JD CHARACTER *8 HMS, name(5) character *11 ctime character *17 cline(5) EXTERNAL JD, HMS C----------------------------------------------------------------------- C TDTOFF, the cummulative number of leap seconds (=55 on 1 Jan 1988) TDTOFF = 60 C----------------------------------------------------------------------- C Input star number and time for observations. C C write(6,'(A,$)') 'Input FK5 star number:' read (5,*) starlist(1) write(6,'(A,$)') ' Input UT date YYYYMMDD:' read (5,*) idate C Parse the date and time iyear = idate / 10000 imon = mod ( idate/100, 100 ) iday = mod ( idate, 100 ) C----------------------------------------------------------------------- C The baseline for this observation - Center from East C This is the delay that has to be applied to the east FDL C C This is in local coordinates, with C C B1 toward the east C B2 toward the north C B3 vertical B1 = -17.489 B2 = 7.247 B3 = 0.0 CON= -17.9841 C....................................................................... C Convert the baseline to "Polar coordinates" C C Bx is toward (hour angle, dec ) = (0,0) C By is toward the East point on the horizon C Bz is toward the Earth's north pole Bx = -B2 * sin ( latitude / RAD ) + B3 * cos ( latitude / RAD ) By = B1 Bz = B2 * cos ( latitude / RAD ) + B3 * sin ( latitude / RAD ) c Bx = -4.1626934 c By = -17.490000 c Bz = 5.9236459 write(6,*) ' bx = ', bx write(6,*) ' by = ', by write(6,*) ' bz = ', bz write(6,*) ' con= ', con open(7,file='dave.out') C----------------------------------------------------------------------- C DJUL is the julian date at 0 hours UT in fractional days C UT1 is the universal time, in fractional days. Since it is C only used in the precession calculation, it need not be C calculated for every observation. Use the value for the C middle of the night. C TJD is the terrestial dynamical time, which differs from UT by C by the inclusion of leap seconds. DJUL = JD ( IDAY, IMON, IYEAR ) hour = 7 UT1 = real(hour) / 24.0 TJD = DJUL + UT1 + DBLE(TDTOFF)/86400.D0 C----------------------------------------------------------------------- C For each star, get C RAM, DECM the mean star position C PMRA, PMDEC, the proper motion C RV, PAR radial velocity and parallax C MAG4, SPECT, the V magnitude, spectral class C ANGDIA the angular diameter in arcseconds istar = 1 CALL CATALOG( starlist(istar), name(istar), $ ram, decm, $ pmra4, pmdec4, PAR4, RV4, MAG4, BV, SPECT, ANGDIA ) write(7,'(A,i5,2x,A,A, f6.1)') $ ' star ', starlist(istar), name(istar), ' V mag = ', mag4 pmra = PMRA4 pmdec = PMDEC4 par = PAR4 rv = RV4 C Convert the catalog positions to apparant places C C APSTAR computes the apparant place of the star (appRA, appDEC) C at the current time corrected for everything except diurnal C abberation. NBODY = 3 CALL APSTAR (TJD, NBODY, RAM, DECM, PMRA, PMDEC, PAR, $ RV, appRA(istar), appDEC(istar) ) C----------------------------------------------------------------------- write ( 7,*) ' year = ', iyear write ( 7,*) ' month = ', imon write ( 7,*) ' day = ', iday write ( 7,*) ' Delay is in millimeters ' write ( 7,*) ' Delay rate is in millimeters / second ' c write(7,1100) ( starlist(istar), name(istar), istar= 1, 5 ) c c 1100 format ( ' time ', 3x, 5(' delay rate '), c $ / ' ut1 sid ', 3x, 5(' meters um/sec '), c $ / ' sec hrs ', 3x, 5( 2x, i5, 2x, A8, 3x ) ) c write (7,1150) 1150 format ( 2('HH:MM utsec fdl1 delay fdl1 rate fdl2 delay ') $ / 2(' mm mm/sec mm ')) do ihour = 0, 6, 1 do minute = 0, 59, 10 do icol = 1,2 hour = ihour + 6 * (icol-1) ut1 = ( hour + minute/60. ) / 24. CALL SIDTIM( DJUL, UT1, 1, GAST ) sideral = gast - longitude/15.0 if ( sideral .lt. 0. ) then sideral = sideral + 24.0 else if ( sideral .gt. 24. ) then sideral = sideral - 24.0 end if write( ctime, '(i6,f5.1)' ) 60*(60*hour+minute), sideral C write(6,*) ' sideral = ', sideral C....................................................................... C SIDTIM determines the sideral time in England (GAST) C DIURN computes the diurnal abberation and outputs the star C position to use for the rest of the calculation. C----------------------------------------------------------------------- c do istar = 1, 5 C CALL DIURN ( TJD, GAST, appRA(istar), appDEC(istar), C $ RA, DEC ) ra = appRA (istar) dec = appDEC(istar) ha = 15.D0 * ( sideral - ra ) IF ( HA .LT. -180.D0 ) then HA = HA + 360.D0 else if ( ha .gt. 180.d0 ) then ha = ha - 360.d0 end if C if ( istar .eq. 1 ) write(6,*) ' ra = ', ra C if ( istar .eq. 1 ) write(6,*) ' dec = ', dec C if ( istar .eq. 1 ) write(6,*) ' ha = ', ha C Calculate delay in millimeters delay = CON + BX*COS(DEC/RAD)*COS(HA/RAD) $ - BY*COS(DEC/RAD)*SIN(HA/RAD) $ + BZ*SIN(DEC/RAD) delay = 1000. * delay C if ( istar .eq. 1 ) write(6,*) ' delay = ', delay C Calculate delay rate in millimeters/second ha = ha + 15.*( 1/3600.) delay_rate = CON + BX*COS(DEC/RAD)*COS(HA/RAD) $ - BY*COS(DEC/RAD)*SIN(HA/RAD) $ + BZ*SIN(DEC/RAD) delay_rate = 1000.*delay_rate - delay ZA = RAD * ACOS( SIN(DEC/RAD)*SIN(latitude/RAD) + $ COS(DEC/RAD)*COS(latitude/RAD)*COS(HA/RAD) ) C C Figure out where to set the delay lines. Only FDL 1 moves. C C Note: on this baseline, the delay is always < 0 C and the delay rate is always > 0 C C C Center hut fdl1(icol) = - delay + 1000. fdl1_dot(icol) = -delay_rate C East hut fdl2(icol) = 1000. utsec(icol) = 3600*hour + 60*minute end do if ( minute .eq. 0 ) write(7,*) ' ' write (7,1200) ihour,minute, utsec(1),fdl1(1),fdl1_dot(1),fdl2(1), $ ihour+6,minute, utsec(2),fdl1(2),fdl1_dot(2),fdl2(2) 1200 format ( 2(i2, ':', i2.2, i6, 3f12.3, 5x) ) end do end do STOP END REAL *8 FUNCTION JD ( IDAY, IMON, IYEAR ) C C Returns the julian day for 0 hours ut. IMPLICIT UNDEFINED(A-Z) SAVE INTEGER *4 IDAY, IMON, IYEAR, YEAR0, YEAR1, DOY PARAMETER (YEAR0=1986) PARAMETER (YEAR1=1995) REAL *8 JD0(YEAR0:YEAR1) EXTERNAL DOY DATA JD0 / 2446430.5D0, 2446795.5D0, 2447160.5, 2447526.5, $ 2447891.5D0, 2448256.5D0, 2448621.5, 2448987.5, $ 2449352.5D0, 2449717.5D0 / IF (( IYEAR .LT. YEAR0 ).OR.( IYEAR .GT. YEAR1 ) ) THEN WRITE(0,'(A,I5)') ' INVALID YEAR FOR COMPUTING JD:', IYEAR ELSE JD = JD0(IYEAR) + DOY( IDAY, IMON, IYEAR ) END IF RETURN END INTEGER *4 FUNCTION DOY( IDAY, IMON, IYEAR ) C C Returns Day of Year for date in arguements ( 1 Jan => doy=1 ) IMPLICIT UNDEFINED(A-Z) SAVE INTEGER *4 IDAY, IMON, IYEAR, LEAP, I, NDAY(12,0:1) DATA NDAY / 0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, $ 0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30 / C LEAP=0 for leap lear LEAP = MIN( 1, MOD( IYEAR, 4 ) ) DOY = IDAY DO 100 I = 1, IMON DOY = DOY + NDAY(I, LEAP) 100 CONTINUE RETURN END SUBROUTINE NEXTDAY ( IDAY, IMON, IYEAR ) C C Increments the date by one C IMPLICIT UNDEFINED(A-Z) INTEGER *4 IDAY, IMON, IYEAR, DAYS(12,0:1), LEAP DATA DAYS / 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31, $ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 / LEAP = MIN ( 1, MOD( IYEAR, 4 ) ) IDAY = IDAY + 1 IF ( IDAY .GT. DAYS(IMON,LEAP) ) THEN IMON = IMON + 1 IDAY = 1 IF ( IMON .EQ. 13 ) THEN IMON = 1 IYEAR = IYEAR + 1 END IF END IF RETURN END SUBROUTINE DIURN( TJD, GAST, RAA, DCA, RAT, DCT ) C C Calculates diurnal abberation C TJD Julian date ( in Terestial dynamical time ) C GAST Greenwich aparent sideral time C RAA Apparent geocentric right ascen. C DCA Apparent geocentric decl. C RAT Topocentric right ascen. ( output ) C DCT Topocentric decl. ( output ) C IMPLICIT UNDEFINED(A-Z) SAVE REAL *8 PI, TWOPI, RADCON, C, ERAD, AU, VROT, DLAT, DLON, HT REAL *8 SINLAT, COSLAT, RHO, ROTVEL, ABRCON REAL *8 TJD, RAA, DCA, RAT, DCT, GAST, HA REAL *8 SINHA, COSHA, SINDEC, COSDEC, ALT, SINALT PARAMETER ( PI = 3.141592653589793D0 ) PARAMETER ( TWOPI = 2.D0*PI ) PARAMETER ( RADCON = PI / 180.D0 ) PARAMETER ( C = 2.99792538D8 ) PARAMETER ( ERAD = 6378140.D0 ) PARAMETER ( AU = 1.49597870D11 ) PARAMETER ( VROT = 7.2921151467D-5 ) C Mt Wilson geocentric coordinates. PARAMETER ( DLAT = 35. + 5./60. + 48.0/3600. ) PARAMETER ( DLON = 111. + 32./60. + 6.0/3600. ) PARAMETER ( HT = 1742.0D0 ) LOGICAL FIRST DATA FIRST / .TRUE. / C IF ( FIRST ) THEN FIRST = .FALSE. SINLAT = DSIN ( DLAT * RADCON ) COSLAT = DCOS ( DLAT * RADCON ) RHO = ( ERAD + HT - 21400.D0 * SINLAT*SINLAT ) / ERAD ROTVEL = VROT * RHO * ERAD ABRCON = ROTVEL * COSLAT / C END IF HA = ( (GAST - RAA ) * 15.D0 - DLON ) * RADCON SINHA = DSIN( HA ) COSHA = DCOS( HA ) SINDEC = DSIN( DCA * RADCON ) COSDEC = DCOS( DCA * RADCON ) C C Test for object above the hroizon C SINALT = SINLAT * SINDEC + COSLAT * COSDEC * COSHA ALT = DASIN ( SINALT ) / RADCON IF ( ALT .LE. -1.D0 ) THEN RAT = 99.9D0 DCT = 99.9D0 ELSE C Apply correction for dirunal abberation RAT = RAA + ( ABRCON * COSHA / COSDEC ) /(15.D0*RADCON) DCT = DCA + ( ABRCON * SINHA * SINDEC ) / RADCON END IF RETURN END SUBROUTINE CATALOG( ISTAR, NAME, RAM, DECM, PMRA, PMDEC, PAR, $ RV, MAG, BV, SPECTRUM, ANGDIA ) C C CATALOG finds a star in the FK5 catalog. Returns star name ( NAME) C position ( RAM, DECM ) in fractional hours and degrees C proper motion ( PMRA, PMDEC ) parallax (PAR), C radial velocity ( RV), V magnitude (MAG) C C written and tested, 7 Dec, 1986 David Mozurkewich C FK5 catalog added, 6 Mar, 1987 C IMPLICIT UNDEFINED (A-Z) SAVE INTEGER*4 ISTAR, IREC, LCAT, I, IERR, INAME(2), STARNO, SPECTRUM REAL*8 RAM, DECM REAL*4 PAR, RV, MAG, PMRA, PMDEC, RAD, ANGDIA, BV CHARACTER *8 NAME LOGICAL FIRST DATA FIRST / .TRUE. / DATA LCAT, RAD / 53, 57.29577951D0 / C----------------------------------------------------------------------- IF ( FIRST ) THEN OPEN ( UNIT = LCAT, FILE='X:\\FK5.BIN', RECL=60, $ FORM='UNFORMATTED', ACCESS='DIRECT', STATUS='OLD', $ IOSTAT=IERR ) IF ( IERR .NE. 0 ) THEN WRITE(0,*) ' ERROR', IERR, ' OPENING CATALOG FILE ' WRITE(0,*) ' CATALOG DIES ' STOP END IF FIRST = .FALSE. END IF C IF ( (ISTAR .LE. 0) .OR. (ISTAR.GT.2000) ) THEN WRITE(0,*) ' Bad FK5 number ', ISTAR RAM = 0.D0 DECM = 0.D0 ELSE READ( LCAT,REC=ISTAR) STARNO, INAME(1), INAME(2), $ RAM, DECM, PMRA, PMDEC, PAR, RV, MAG, BV, SPECTRUM, ANGDIA IF ( STARNO .NE. ISTAR ) THEN WRITE(0,*) ' CATALOG is corrupt! ' WRITE(0,*) ' Record ', ISTAR, ' contains star ', STARNO WRITE(0,*) ' CATALOG bombs ' STOP END IF WRITE(NAME,'(2A4)') INAME(1), INAME(2) END IF RETURN 200 CONTINUE WRITE(0,1021) I, IERR STOP 1021 FORMAT ( ' fatal error while reading FK5.bin catalog '/ $ ' error number:', I5, ' record number:', I5/ ' too bad !') END