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