Program ObsNterp
c
c Program ObsNterp evaluates the Comet Ephemeris given the JED Observation
c Time. The user must input the Interpolation Time for the ephemeris to be
c evaluated.
c The computations are performed using the Lagrange method of interpolatio
c
c The number of known points is fixed at seven(7).
c
c_____________________________________________________________________________
c
c .. Designed and Implemented by:___ Ravenel N. Wimberly ___
c Sterling Software
c Jet Propulsion Laboratory
c Astrometry Network
c International Halley Watch
c_____________________________________________________________________________
c
c .. Declare Variables
c ..
Character Flag*1, Month*6, Filenm*14
Double Precision Tobs, Values(11), Ramn, Decmn
Integer Year, Imonth, Id, Ih, Key
Real Day
c
c ... Clear Screen
Write(*,'(//////////////////////////)')
c_____________________________________________________________________________
c ... Prompt User for Comet Ephemeris File Name
1 write(*,'(a)') ' What is the name of your Comet Ephemeris File ? '
read (*,'(a)',end=20) Filenm
c ... Open Comet Ephemeris file for processing
open(10,file=filenm,status='old',access='sequential',
* form='formatted')
c_____________________________________________________________________________
c______________________PROCESS DATA___________________________________________
c ... Clear Screen
10 Write(*,'(//////////////////////////)')
c
Write(*,*) ' At the Prompt Please Input Epoch for Interpolation'
c
write(*,'(a)') ' Please Enter Year of Interest (1985) --> '
read (*,*) Year
write(*,'(a)') ' Please Enter Month of Interest (09) --> '
read (*,*) Imonth
write(*,'(a)') ' Please Enter Day of Interest (23.47) --> '
read (*,*) Day
c ____________Set flag for Gregorian Calendar___________
Flag='G'
Tobs=0.0d0
c ----------Call Jdate to Obtain Julian Date for Interpolation---------
call jdate(Flag,Year,Imonth,Day,Month,Tobs)
c __________Display Date Results__________
write(*,'(/,10h Date Is: ,i4,a,f10.5,3h = ,f16.7,//)')
* Year,Month,Day,Tobs
c_____________________________________________________________________________
c ... Evaluate Comet Ephemeris at Time TOBS
c
call Interp(Tobs,Values)
c
c ... Convert Ra & Dec from radians to Hours and Degrees
key= -1
call argch2(key,Values(1),Values(2),Ih,Ramn,Flag,Id,Decmn)
c_____________________________________________________________________________
c ... Display Results of Interpolation
write(*,'(1x,3hRA=,i4,f7.3,5h DEC=,a,i3,f6.2,7h Delta=,f8.4,8h DEL
*DOT=,f8.4,3h R=,f7.4,/,6h RDOT=,f8.4,7h THETA=,f5.1,6h BETA=,f5.1,
*6h MOON=,f6.1,7h PSANG=,f6.1,7h PSAMV=,f6.1,//)') Ih,Ramn,Flag,Id,
* Decmn,(values(i),i=3,11)
c_____________________________________________________________________________
c ... Determine if user wants to interpolate with this set of data
c
Write(*,*) ' Interpolate More Points With This Data Set ?? -'
Write(*,'(a)') ' (1 = YES, 0 = NO) -----> '
Read (*,*,err=10,end=10) ians
If(ians.eq.1) go to 10
c
c ... Determine if user wants to interpolate with a new set of data
c
Close(10,status='keep')
Write(*,*) ' Interpolate With A New Set of Data ?? -'
Write(*,'(a)') ' (1 = YES, 0 = NO) -----> '
Read (*,*,err=1,end=1) ians
If(ians.eq.1) go to 1
c
c ... If not.. Clear screen and terminate program
c
20 Write(*,'(//////////////////////////)')
Close(10,status='keep')
End
SUBROUTINE JDATE(FLAG, JAHR,MONTH,DAY,CNAME, DJD)
C SUBROUTINE DATE COMPUTES JULIAN DATE FROM CIVIL DATE,OR VICE VERSA.
C JAHR IS YEAR AD OR BC(INTEGER),DJD IS DOUBLE PRECISION JULIAN DATE.
C JAHR AND MONTH ARE INTEGERS, DAY A SINGLE PRECISION VARIABLE.
C IF DJD IS TO BE COMPUTED, IT MUST BE ENTERED AS ZERO.
C CIVIL DATE MAY BE EITHER JULIAN OR GREGORIAN CALENDAR,DEPENDING ON FLAG.
C FLAG IS 'G' FOR GREGORIAN CIVIL DATE, 'J' FOR JULIAN CIVIL DATE.
C CNAME IS THE NAME OF THE MONTH,SUITABLE FOR WRITING OUT (A6).
DIMENSION NO(12),DEPOCH(2)
CHARACTER*6 VNAME(12), CNAME
CHARACTER*1 FLAG
INTEGER Y, JAHRJ, JAHRG
DOUBLE PRECISION DJD,DEPOCH,DMJD
DATA VNAME/ ' JAN ',' FEB ',' MAR ',' APR ',' MAY ',' JUNE',
1 ' JULY',' AUG ',' SEPT',' OCT ',' NOV ',' DEC '/
C JD 2159716.5 IS 1 JAN 1201AD (GREGORIAN CALENDAR).
DATA DEPOCH / 260423.5D0, 2159716.5D0/
DATA NO/31,28,31,30,31,30,31,31,30,31,30,31/
C
IF ( DJD.GT.1.D-13 ) GO TO 50
Y = JAHR
M = MONTH
IF ( FLAG.EQ.'G' ) GO TO 10
C THE MULLER-WIMBERLY STATEMENT FOR JULIAN CALENDAR FOLLOWS.
J=367*Y-7*(Y+5001+(M-9)/7)/4+275*M/9+1729777
C DBLE(FLOAT(J)) MAY NOT GIVE FULL PRECISION ON SOME MACHINES.
DJD = DBLE(DAY+FLOAT(J) ) -.5D0
C
CNAME = VNAME(MONTH)
RETURN
10 CONTINUE
C THE MULLER-WIMBERLY STATEMENT FOR GREGORIAN CALENDAR FOLLOWS.
J=367*Y-7*(Y+(M+9)/12)/4-3*((Y+(M-9)/7)/100+1)/4+275*M/9+1721029
DJD = DAY + DBLE(J) - .5D0
C
CNAME = VNAME(MONTH)
RETURN
50 CONTINUE
JUMP = 0
MODE = 2
IF ( FLAG.EQ.'J' ) MODE = 1
DMJD=DJD-DEPOCH(MODE)
L=DMJD
FL=FLOAT(L)
REMAIN=DMJD-DBLE(FL)
LEAP = 0
GO TO (51,52) , MODE
51 JAHR = JAHRJ(L)
IF ( MOD(JAHR,4).EQ.0) JUMP = 1
GO TO 60
52 N400 = L/146097
IF ( L.EQ.146097 ) LEAP = 1
M400 = MOD(L,146097) - LEAP
N100 = M400/36524
L = MOD( M400,36524 )
JAHR = JAHRG( L,N100,N400 )
IF (MOD(JAHR,4).NE.0) GO TO 60
IF ( MOD(JAHR,100).NE.0) JUMP = 1
IF ( MOD(JAHR,400).EQ.0 ) JUMP = 1
60 CONTINUE
L = L + LEAP
L = MOD (L,1461)
L = (L/1460)*365 + MOD(L,365) + 1
DO 55 MONTH = 1,12
N = NO(MONTH)
IF (MONTH.EQ.2) N = 28+JUMP
L = L-N
IF (L.LE.0) GO TO 56
55 CONTINUE
56 CONTINUE
CNAME = VNAME(MONTH)
DAY = FLOAT(L+N) + REMAIN
RETURN
END
INTEGER FUNCTION JAHRJ(IL)
INTEGER IL
C
JAHRJ = (IL -IL/1460)/365 - 3999
C
RETURN
END
INTEGER FUNCTION JAHRG(IL,I100,I400)
INTEGER IL,I100,I400
C
JAHRG = (IL-IL/1460)/365+100*I100+400*I400+1201
C
RETURN
END
Subroutine Interp(tobs,v)
c************************************************************
c..
c Interp Reads and Interpolates Comet Ephemeris File at Time
c Tobs and Uses a Seven Point LAGRANGIAN INTERPOLATION
c************************************************************
c..INPUT
c Tobs = D.P. JED Observation Time When COMET EPHEMERIS
c FILE is to be Interpolated.
c
c..OUTPUT
c V = (V(I),I=1,11)=D.P. Interpolated Comet Values
c***********************************************************
Implicit Double Precision(A-H,O-Z)
Character*1 Isign
Dimension Val(11,7), T(7), V(11)
Integer key, Ih, Id, Iread
c
Data Iread/1/, Oldtob/0.0d0/, rad/0.0174532925199D0/
c _________Set Output Vector To Zero___________
Do 1 mz=1,11
V(mz)= 0.0d0
1 Continue
c ___________________________
if(Oldtob.gt.Tobs) Iread= 1
if(Iread.ne.1) go to 10
Rewind 10
c ___________________________
c ____Store 7 Records of COMET Ephemeris into val____
c
DO 5 I=1,7
Read(10,'(16x,f10.1,(i3,f7.3),(1x,a,i2,f6.2),2(1x,f7.4,f9.4),
* f6.1,2(1x,f5.1),2f6.1)',end=80)
* T(I),Ih,Ramn,Isign,Id,Decmn,(val(M,I),M=3,11)
c ___________________Convert Ra & Dec to Radians____________________
key= +1
call argch2(key,Val(1,I),Val(2,I),Ih,Ramn,Isign,Id,Decmn)
c ------------------------------------------------------------------
5 CONTINUE
c
Iread=2
c
10 continue
c ________________________
IF(TOBS.EQ.T(4)) then
do 20 m= 1,11
v(m) = val(m,4)
20 continue
go to 70
ENDIF
c ________________________________________________
IF((TOBS-T(4))*(TOBS-T(3)).LT.0.0D0 .or.
* (TOBS-T(5))*(TOBS-T(4)).LT.0.0D0) go to 30
c ________________________________________________
c
c TOBS IS NOT BETWEEN REC 3 AND 4, OR 4 AND 5.. SO MOVE 6 RECORDS UP IN STACK
Do 28 i=1,6
T(i)=T(i+1)
do 24 m=1,11
Val(m,i)=Val(m,i+1)
24 continue
28 Continue
c
c _____________ NOW READ RECORD INTO 7TH AREA OF STACK_________________
Read(10,'(16x,f10.1,(i3,f7.3),(1x,a,i2,f6.2),2(1x,f7.4,f9.4),
* f6.1,2(1x,f5.1),2f6.1)',end=80)
* T(7),Ih,Ramn,Isign,Id,Decmn,(val(M,7),M=3,11)
c ___________________Convert Ra & Dec to Radians____________________
key= +1
call argch2(key,Val(1,7),Val(2,7),Ih,Ramn,Isign,Id,Decmn)
c ------------------------------------------------------------------
go to 10
c _______________________________________________
30 Do 60 iio= 1,7
x= 1
do 40 iin= 1,7
if(iin.eq.iio) go to 40
x= x * (Tobs - T(iin)) / (T(iio) - T(iin))
40 continue
do 50 m= 1,11
if(m.eq.1) then
if(((Val(1,7)/rad)-(Val(1,iio)/rad)).gt.345.0d0) then
v(m) = v(m) + x * (val(1,iio)+(360.0d0*rad))
go to 50
endif
endif
v(m) = v(m) + x * val(m,iio)
50 continue
60 Continue
c _______________________________________________
70 Oldtob = Tobs
return
80 Oldtob = 9999999.0d0
write(*,'(///,46h Observation Time is not in this Ephemeris Set)')
return
End
SUBROUTINE ARGCH2(KEY,RA,DEC,IH,AM,ISGN,ID,AMIN)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CHARACTER*1 ISGN,IMINUS,IPLUS
DATA IPLUS,IMINUS /'+','-'/
DATA RAD/0.0174532925199D0/
C****
C**** KEY .NE. -1 ---- HOURS AND DEGREES TO RADIANS
C**** KEY .EQ. -1 ---- RADIANS TO HOURS AND DEGREES
C****
IF(KEY.EQ.-1) GO TO 10
H=IH
D=ID
RA=15.D0*(H+AM/6.D1)* RAD
DEC=(D+AMIN/6.D1)*RAD
IF(ISGN.EQ.'-') DEC=-DEC
RETURN
10 HOURS=RA/RAD/15.D0
HOURS=DMOD(HOURS,24.D0)
IH=HOURS
AM=(HOURS-IH)*6.D1
AM=AM+1.D-6
IF(AM.LT.6.D1) GO TO 12
AM=AM-6.D1
IH=IH+1
12 IF(IH.GE.24) IH=IH-24
IF(DEC) 15,20,20
15 ISGN=IMINUS
GO TO 25
20 ISGN=IPLUS
25 DEG=DABS(DEC)/RAD
ID=DEG
AMIN=(DEG-ID)*6.D1
RETURN
END