C*************************************************************************** C C &PLNET::16 C C This series of subroutines is designed to read a MICA ephemeris from a C file FILNAM on cartridge 16 and interpolate ( either linearly or quad- C ratically where possible ) to the current time. The file is assumed to C have ephemerides in strict increasing time order only. This ephemeris C file should be generated using MICA ( Multi-year Interactive Computer C Almanac ) for each suitable planet using C C .. .. <(Planet)> .. .. C .. C C for a suitable set of dates and UT1 (**NB**) times. Note that all the C header lines etc will be skipped. Note also that Calendar not Julian C dates should be used ( option chosen for User readability ). C C The use of Geocentirc/Topocentric co-ordinates is dependent on the C planet concerned. Topocentric effects for Jupiter are less that 1mdeg C but for Venus at closest approach, one may be as much as 9 mdeg off. C Note that the time spacing of Topocentric co-ordinates should be C fine enough to resolve the diurnal topocentric effect. Geocentric C co-ordinates should be spaced so as to suitably resolve the planet C motion against the stars. Topocentric is recommended for at least C Venus and Mars and perhaps Mercury though the Sun would probably C render such obervations impossible, C C*************************************************************************** C C Subroutine PLANET (LU, FILNAM, RA, DEC, IERR) C ================= C To retrieve current interpolated position of any planet from C MICA output lists contained in file FILNAM::16 C C inputs : LU Integer LU of error display C FILNAM Integer(3) Source file name C C outputs: RA Single Precision Topocentric RA of date C DEC Single Precision Topocentric DEC of date C IERR Integer IOSTAT IERR of I/O C Implicit None C Integer TabSize C Parameter ( TabSize = 4 ) C Real Interpolate C Integer LU, FILNAM(3), IERR C Real RA, DEC C Logical LOOPED C Character Filename*10, Ephemeris*10 C Double Precision MJDtab(TabSize), MJDnow, MJDepoch C Real RAtab(TabSize), DECtab(TabSize), RAcoeff(3), DECcoeff(3) C Integer Order, I, IOF, IOL, INX, IPR, LUfile C Integer*2 IYEAR, ITIME(5) C Common /MEMORY/ Ephemeris,MJDepoch,RAcoeff,DECcoeff,Order C C Initialise C IERR = 0 C C Determine current MJD from computer clock C Call EXEC(11,ITIME,IYEAR) Call MJDAY(IYEAR,ITIME,MJDnow) C C Check Validity of current ephemeris coefficients C Write (Filename,'(3A2,"::16")') (FILNAM(I), I=1,3) If (Filename .ne. Ephemeris .or. MJDnow .gt. MJDepoch) then C C Load suitable table of values for interpolation from file C IOF = 0 LOOPED = .FALSE. Call Open_File(Filename,'OLD','START',LUfile,IERR) If (IERR .ne. 0) then Write (LU,'(1X,"Could not open ",A)') Filename Go to 999 endif C C Loop, reading next ephemeris entry until a future entry is found C 100 If (IOF .eq. TabSize) LOOPED = .TRUE. IOF = mod(IOF,TabSize) + 1 Call Get_Entry(LUfile,MJDtab(IOF),RAtab(IOF), * DECtab(IOF),IERR) If (IERR .gt. 0) * Write (LU,'(1X,"I/O Error",I5," on ",A)') IERR, Filename If (IERR .ne. 0) Go to 999 If (MJDtab(IOF) .lt. MJDnow) Go to 100 If (IOF .lt. 2 .and. .not. LOOPED) Go to 990 C C Try see if next entry is better suited than old historical entry C IOL = mod(IOF-3+TabSize,TabSize) + 1 IPR = mod(IOF-2+TabSize,TabSize) + 1 INX = mod(IOF ,TabSize) + 1 Call Get_Entry(LUfile,MJDtab(INX),RAtab(INX),DECtab(INX),IERR) C C Close Ephemeris File again to minimise competition for resource C Close ( UNIT=LUfile, ERR=150 ) C If (IERR .eq. 0) then If ( (IOF .lt. 3 .and. .not. LOOPED) .or. * ((MJDtab(INX)-MJDnow) .le. (MJDnow - MJDtab(IOL))) ) then IF (IOF .eq. TabSize) LOOPED = .TRUE. IOL = IPR IPR = IOL IOF = INX endif endif IERR = 0 C C Check times of values ( within 2 days reqd ) according to Order C 150 If (IOF .lt. 3 .and. .not. LOOPED) Go to 200 C If ((MJDtab(IOF)-MJDnow) .gt. 2.0) then C C Try as linear with other two values where possible C If (MJDtab(IPR) .lt. MJDnow) Go to 990 IOF = IPR IPR = IOL Go to 200 C endif If ((MJDnow-MJDtab(IOL)) .gt. 2.0) Go to 200 Order = 2 C Go to 300 C 200 If ((MJDtab(IOF)-MJDnow) .gt. 2.0 .or. * (MJDnow-MJDtab(IPR)) .gt. 2.0 ) Go to 990 Order = 1 C C Calculate suitable interpolation coefficients C 300 Call Find_Coeff(Order,MJDtab,RAtab, TabSize,IOF,RAcoeff ) Call Find_Coeff(Order,MJDtab,DECtab,TabSize,IOF,DECcoeff) MJDepoch = MJDtab(IOF) C C Ephemeris is Complete ! C Ephemeris = Filename endif C C Interpolate topocentric RA and DEC of date to present time C RA = Interpolate(MJDnow,RAcoeff, Order,MJDepoch) DEC = Interpolate(MJDnow,DECcoeff,Order,MJDepoch) Go to 999 C C Error trap for unuseable ephemeris values C 990 Write (LU,'(1X,"Unable to interpolate Ephemeris ",A)') Filename IERR = -999 C 999 Return End C C Subroutine MJDAY (IYR, ITIME, MJD) C ================ C To calculate the Modified Julian Day, as defined in the C Astronomical Almanac, for a given year, day & time C C inputs : IYR Integer*2 Year ( from 1983 - ) C ITIME(5) Integer*2 Day No, Hr, Min, Sec, 100th C C outputs: MJD Double Precision Modified Julian Day Number C C Implicit None C Integer*2 IYR, ITIME(5) C Integer JYR C Double Precision MJD C Logical LEAPYR C MJD = ( ( (DBLE(ITIME(1))/100D0 + DBLE(ITIME(2)) ) / 60.0D0 * + DBLE(ITIME(3)) ) / 60.0D0 * + DBLE(ITIME(4)) ) / 24.0D0 * + DBLE(ITIME(5)) C If (IYR .lt. 1983) then C Print '(I4," - year not catered for in MJD routine MJDAY")', * IYR MJD = 0.0D0 C else C MJD = MJD + 45334.0D0 ! 1983 offset DO 100 JYR = 1984, IYR MJD = MJD + 365.0D0 If (LEAPYR(JYR-1)) MJD = MJD + 1.0D0 100 CONTINUE C endif C Return End C C Logical Function LEAPYR (IYEAR) C ======================= C To determine if any given year is a leap year C C inputs : IYEAR Integer Year to be assessed C C outputs: LEAPYR Logical .TRUE. if IYEAR is a leap year C Implicit None C Integer IYEAR C LEAPYR = (((MOD(IYEAR,4) .eq. 0) .AND. (MOD(IYEAR,100) .ne. 0)) * .OR. (MOD(IYEAR,400) .eq. 0)) C Return End C C Subroutine Open_File ( Filename, FileStat, Position, LU, IERR ) C ==================== C To open any sequential file if not already opened and to return C the LU number associated with the file. New LU's are allocated C sequentially from 65 upwards as necessary. Note that the file C must match the status FileStat or the open will fail. C C inputs : Filename Character*(*) Name of file to be opened C FileStat Character*(*) 'OLD', 'NEW' or 'UNKNOWN' C Position Character*(*) 'START' or 'END' as reqd C C outputs: LU Integer Logical Unit attatched C IERR Integer I/O IOSTAT error number C Implicit None C Character Filename*(*), FileStat*(*), Position*(*), Junk C Integer LU, Default, IERR C Logical Exists, OpenNow, Used C Data Default / 65 / C C Initialise C IERR = 0 Call Upcase (Position) C C first we must determine if the required file is open already C Inquire ( FILE=Filename, IOSTAT=IERR, ERR=999, * EXIST=Exists, OPENED=OpenNow, NUMBER=LU ) C If (.not. OpenNow) then C C if file is not yet opened, determine free LU to use for file C Used = .TRUE. LU = Default - 1 Do While (Used) LU = LU + 1 Inquire ( UNIT=LU, IOSTAT=IERR, ERR=999, OPENED=Used ) enddo C C and open required file FILENAME, ensuring it exists already C Open ( FILE=Filename, UNIT=LU, IOSTAT=IERR, ERR=999, * STATUS=FileStat, ACCESS='SEQUENTIAL' ) C endif C C if Position='START' then simply rewind file and Exit C If (Position .eq. 'START') then Rewind ( UNIT=LU, IOSTAT=IERR, ERR=999 ) goto 999 endif C C if Position='END' rewind, then locate EOF C If (Position .eq. 'END') then Rewind ( UNIT=LU, IOSTAT=IERR, ERR=999 ) Do While (.TRUE.) Read ( LU, '(A1)', IOSTAT=IERR, ERR=999, END=900 ) Junk enddo 900 IERR = 0 Go to 999 endif C 990 Print '("Unknown Position ",A," in Open_File")', Position IERR = -999 C 999 Return End C C Subroutine Get_Entry ( FileLU, MJD, RA, DEC, IERR) C =================== C To get the next ephemeris value from the file attatched to logical C unit FileLU. The next values of RA and DEC found in the file will C be returned, along with the corresponding Modified Julian Date. C C inputs : FileLU Integer LU associated with file C C outputs: MJD Double Precision MJD of ephemeris retrieved C RA Real Right Ascension Ephemeris C DEC Real Declination Ephemeris C IERR Integer IOSTAT ierr value if any C Implicit None C Logical LEAPYR C Character Buffer*80, Month(12)*3, MON*3, DSGN*1 C Integer FileLU, IERR, I, IDY, RAH, RAM, DCD, DCM, Length(12) C Integer*2 IYEAR, ITIME(5) C Double Precision MJD C Real RA, DEC, RAS, DCS C Data Month / 'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', * 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec' /, * Length / 31, 28, 31, 30, 31, 30, * 31, 31, 30, 31, 30, 31 / C C Initialise C IERR = 0 C C Read next line in file till decodable or EOF/error detected C 100 Read (FileLU, '(A80)', END=990, IOSTAT=IERR, ERR=999) Buffer IYEAR = 0 Read (Buffer, 1010, ERR=100) IYEAR, MON, IDY, * (ITIME(I), I=4,1,-1), RAH, RAM, RAS, DSGN, DCD, DCM, DCS 1010 Format(I4,1X,A3,4(1X,I2),1X,I1,5X,2(I2,1X),F6.3, * 5X,A1,1X,2(I2,1X),F5.2) If (IYEAR .EQ. 0) Go to 100 C C Correct 10ms and DOY components of ITIME, then find MJD C ITIME(1) = ITIME(1)*10 ITIME(5) = 0 DO 200 I = 1,12 If (MON .eq. Month(I)) Go to 300 ITIME(5) = ITIME(5) + Length(I) IF (I .eq. 2 .and. LEAPYR(IYEAR)) ITIME(5) = ITIME(5) + 1 200 CONTINUE Print '(A3," - illegal month descriptor in Get_Entry")',MON IERR = -111 Go to 999 300 ITIME(5) = ITIME(5) + IDY Call MJDAY(IYEAR,ITIME,MJD) C C Decimalise and return RA and DEC values C RA = 15.0*(RAH + RAM/60.0 + RAS/3600.0) DEC = DCD + DCM/60.0 + DCS/3600.0 IF (DSGN .eq. '-') DEC = -DEC Go to 999 C 990 Print '("End of file encountered on Ephemeris")' 999 Return End C C C Subroutine Find_Coeff(Order, Ttab, Xtab, Size, IOF, Xcoef) C ===================== C To determine the coefficients of a polynomial fit of order Order C ( only 1 or 2 allowed ) to independent variable in Xtab as a C function of dependent variable Ttab. IOF is the index of the C last value to fit, and the polynomial is constructed about this C point. Both Ttab, Xtab are taken to be circular lists of length C Size. No checking is performed. C C inputs : Order Integer Order of Polynomial to fit C Ttab(*) Double Precision Dependent Var (circ. list) C Xtab(*) Real Independent Var ( " ) C Size Integer Modulo of Circular Lists C IOF Integer Index of last value in list C C outputs: Xcoef(*) Real Polynomial Coefficients C C Implicit None C Integer Order, Size, IOF, IOL, IPR C Double Precision Ttab(*) C Real Xtab(*), Xcoef(*) C C Setup suitable pointers to list values IOL,IPR,IOF in order C IOL = mod(IOF-3+Size,Size) + 1 IPR = mod(IOF-2+Size,Size) + 1 C C Zero order coefficient is always value at IOF C XCoef(1) = Xtab(IOF) C C Linear Fit consists simply of slope C If (Order .eq. 1) then Xcoef(2) = (Xtab(IPR)-Xtab(IOF)) / (Ttab(IPR)-Ttab(IOF)) C C Quadratic fit is more complicated C elseif (Order .eq. 2) then Xcoef(2) = ((Xtab(IOL)-Xtab(IOF))*(Ttab(IPR)-Ttab(IOF))**2 * -(Xtab(IPR)-Xtab(IOF))*(Ttab(IOL)-Ttab(IOF))**2 ) * /(Ttab(IPR)-Ttab(IOL))/(Ttab(IOL)-Ttab(IOF)) * /(Ttab(IPR)-Ttab(IOF)) Xcoef(3) = ((Xtab(IOL)-Xtab(IOF))*(Ttab(IPR)-Ttab(IOF)) * -(Xtab(IPR)-Xtab(IOF))*(Ttab(IOL)-Ttab(IOF)) ) * /(Ttab(IOL)-Ttab(IPR))/(Ttab(IOL)-Ttab(IOF)) * /(Ttab(IPR)-Ttab(IOF)) C C Anything beyond that is not necessary here C else Print '(I2," - Order not catered for Find_Coeff routine")', * Order endif C Return End C C Real Function Interpolate (Tnow, Xcoef, Order, Tepoch) C ======================= C To evaluate a set of polynomial coefficients of order Order C at a specific value. C C inputs : Tnow Double Precision Value of Independent variable C Xcoef(*) Real Polynomial Coefficients in T C Order Integer Order of polynomial expansion C Tepoch Double Precision Value about which expanded C C outputs: Inte... Real Dependent value at Tnow C Implicit None C Double Precision Tnow, Tepoch C Real Xcoef(*) C Integer Order, I C Interpolate = Xcoef(1) Do 100 I = 1, Order Interpolate = Interpolate + Xcoef(I+1)*(Tnow-Tepoch)**I 100 Continue C Return End C======================= INCLUDE UPCASE::16 C=======================