Viewing contents of file '../idllib/iuedac/iuelib/pro/findelt.pro'
;******************************************************************************
;+
;*NAME:
;
; FINDELT
;
;*CLASS:
;
;*PURPOSE:
;
; This program will search the file ELEMENTS.LIS for the set of orbital
; elements appropriate for the given observation date. If the image is
; too recent to have an appropriate set of elements in the list, or if
; there is an error opening the list file, the procedure ORBEL is called
; to extract the elements from the science image header, if an image was
; given in the first place. Failing all else, the user will be prompted
; for the elements.
;
;*CALLING SEQUENCE:
;
; FINDELT,IMAGET,RJD,JD,A,E,I,O,W,M
;
;*PARAMETERS:
;
; IMAGET (REQ) (I) (0) (S) - Image file name (cam+isn+type)
;
; RJD (REQ) (I) (0) (D) - Julian date of observation
;
; JD (REQ) (O) (0) (D) - Julian date - epoch of orbital elements
;
; A (REQ) (O) (0) (D) - Semi-major axis length of orbit in km
;
; E (REQ) (O) (0) (D) - Eccentricity of the orbit
;
; I (REQ) (O) (0) (D) - Inclination of the orbit
;
; O (REQ) (O) (0) (D) - Longitude of the ascending node
;
; W (REQ) (O) (0) (D) - Argument of pericenter
;
; M (REQ) (O) (0) (D) - Mean anomaly
;
;*SYSTEM VARIABLES:
;
; !err
; !iuer.dat
;
;*SUBROUTINES CALLED:
;
; PARCHECK
; ORBEL
;
;*FILES USED:
;
; ELEMENTS.LIS
; TEMP.LIS
;
;*NOTES:
;
; This program chooses the element set on or nearest in time to the
; observation date where there is no intervening delta V manuever.
; Note that the element set chosen may not be consistent with the
; element set found in the science image header, which is the most
; recent previous set to the observation date, regardless of possible
; intervening delta V manuevers.
; Images may be taken immediately following a delta V maneuver, when the
; old (and now invalid) set of elements is still loaded on the operations
; computer (and thus in the header).
;
; tested with IDL Version 2.1.0 (sunos sparc) 26 Jul 91
; tested with IDL Version 2.1.0 (ultrix mispel) N/A
; tested with IDL Version 2.1.0 (vms vax) 25 July 1991 by LLT
;
;*MODIFICATION HISTORY:
;
; Written 25 January 1990 by LLT.
; 4/9/90 LLT Prolog
; 10/90 LLT Modified
; 11/5/90 LLT Read elements from ELEMENTS.LIS
; 7-dec-90 LLT Change the way delta Vs are handled.
; 19-dec-90 LLT Add call to ORBEL if date is too recent. Remove option
; of using RADCOR's elements.
; 4-Jan-91 LLT Choose the nearest set instead of previous set.
; 12-Mar-91 LLT Use temp file to read element set chosen; thus avoiding
; tab characters and blanks in different versions of the
; ELEMENTS.LIS file.
; 4-11-91 GRA added branch for file name processing on vms
; 25July91 LLT clean up, add parcheck, update prolog, change to new
; logicals, fix branch to check science image header.
; tested on VAX
; 26July91 PJL tested on SUN; updated prolog
; 2 Sep 94 LLT replace IUER_DAT with !iuer.dat
;-
;******************************************************************************
pro findelt,imaget,rjd,jd,a,e,i,o,w,m
npar=n_params(0)
if npar eq 0 then begin
print,'FINDELT,IMAGET,RJD,JD,A,E,I,O,W,M'
retall
endif ;npar eq 0
parcheck,npar,9,'FINDELT'
jd=0.d0 & epoch=0.d0 & a=0.d0 & e=0.d0 & i=0.d0 & o=0.d0 & w=0.d0 & m=0.d0
; Below there are some options that we are no longer offering to the user...
;...........................................................................
answer=' '
;print,'Default orbital elements are those used by RADCOR (22 Nov 1979).'
;read,'Do you want to use the default elements [def=no]',answer
;if strupcase(strmid(answer,0,1)) eq 'Y' then begin
; jd=44199.5d & a=42163.2d & e=0.2359693d
; o=193.96197d & w=270.9129979d & m=246.56d & i=28.27283731d
; return & endif
yn=' '
;if imaget ne '' then $
;read,'Would you prefer to use the elements in the label [y/n; def=no]',yn
;if strpos(strupcase(strmid(yn,0,1)),'Y') ge 0 then goto,sih
;...........................................................................
; Get ready to search the elements.lis file for the orbital elements
print,'Ready to check elements.lis'
dn=[0,31,59,90,120,151,181,212,243,273,304,334] ;Day number starting each month
yjd0=0.d0
yjd=0.d0
get_lun,u
on_ioerror,sih ;Check science image header if can't open file
openr,u,!iuer.dat + 'elements.lis' ;File orbital elements
on_ioerror,label ;On eof or io error, do some tests or check label
line0=string(bytarr(132))
line1=strarr(1)
;repeat readf,u,line1 until strmid(line1,2,4) eq '----' ;Skip column headings
;on_matherror,continue
repeat begin
line0=line1(0) & yjd0=yjd ;Save line and date
continue:
;Read a new line,
repeat readf,u,line1 until $ ;skipping comments
((strmid(line1(0),1,1) ge '0') and (strmid(line1(0),1,1) le '9'))
year=fix(strmid(line1(0),1,2)) ;Extract the date
day=double(strmid(line1(0),5,2)) ;from line1, get
mon=fix(strmid(line1(0),3,2)) ;day number.
day=day+dn(mon-1)
if mon gt 2 then begin ;If year is a leap year,
day=day+(year/4.0 eq fix(year/4.0)) ;and if month is March
if (fix(year/100.) eq year/100.) then $ ;or later, add a day to
day=day-(year/400. ne fix(year/400.)) ;the day number.
endif ;later than February
yjd=43508.5d0+365.d0*(year-78.d0)+ $ ;Calculate the Julian
1.d0*fix((year-77.d0)/4.d0+1.d-04)+day ;date. Repeat if date
endrep until ((yjd gt rjd)) ;is le obs. date.
; Determine which set of elements is closer to the observation date, and use
; that set if there is no intervening delta-V maneuver.
dyjd=yjd-rjd & dyjd0=rjd-yjd0
if strpos(strupcase(line0),'DELTA') ge 0 then dyjd0=1000
if strpos(strupcase(line1(0)),'DELTA') ge 0 then dyjd=1000
if dyjd0 gt dyjd then begin
line0=line1(0)
yjd0=yjd
endif ;dyjd0 gt dyjd
useline0:
print,'Epoch of element set used: ',strmid(line0,1,6),' ',string(yjd0)
jd=yjd0 ;Julian date of epoch
; Get rid of blanks, tab characters, etc. Different versions of the element
; table have had different combinations of blanks and tab characters, making
; it difficult to extract the elements directly from the string.
openw,tun,'temp.tmp',/get_lun ;Open temporary file for writing
printf,tun,line0 ;Write line0 to this temporary file
close,tun ;Close the file
openr,tun,'temp.tmp' ;Open the file again - this time to read
readf,tun,temp,a,e,i,o,w,m ;Read the orbital elements as numbers
close,tun ;Close the file
free_lun,tun
;a=double(strmid(line0,9,7)) ;Semimajor axis (km)
;e=double(strmid(line0,18,9)) ;Eccentricity
;i=double(strmid(line0,28,7)) ;Inclination (deg)
;o=double(strmid(line0,37,7)) ;R.A. of ascending node (deg)
;w=double(strmid(line0,46,7)) ;Argument of perigee (deg)
;m=double(strmid(line0,55,7)) ;Mean anomaly at epoch (deg)
free_lun,u
return
;Error opening ELEMENTS.LIS, or obs. date is too recent
label:
if eof(u) then print,'End-of-file' else begin
print,strmessage(-!err) ,!err
print,line0
print,line1
print,year,mon,day
endelse ;not the end-of-file
; If the last element set is less than two weeks earlier than the observation
; date, then just go ahead and use that last set.
if (eof(u) and ((rjd-yjd0) lt 15.0) and strpos(strupcase(line0),'DELTA') lt 0) $
then goto,useline0
on_ioerror,ask ;Ask user for elements if all else fails
if yjd eq rjd then begin
line0=line1
yjd0=yjd
goto,useline0
endif ; yjd
sih:
if (yjd0 eq 0.d0) and (yjd ne rjd) then print,'Could not open ELEMENTS.LIS'
if imaget ne '' then begin ;Look in the science image header for the
orbel,imaget,orr,ord ;elements.
jd=ord(0)+ord(1)/86400.d0 ;Julian date of epoch
if jd ne 0 then begin ;There were elements in the label
a=ord(2) ;Semimajor axis (km)
e=ord(3) ;Eccentricity
i=ord(4) ;Inclination (deg)
o=ord(5) ;R.A. of ascending node (deg)
w=ord(6) ;Argument of pericenter (deg)
m=ord(7) ;Mean anomaly at epoch (deg)
return
endif ;jd ne 0
endif ; imaget
; Ask the user for the elements.
ask:
print,'Enter the orbital elements of IUE at the time of observation:'
ymd=' '
print,'If you have the epoch in the form of a Julian date, hit return.'
read,'What is the epoch [yymmdd]? ',ymd
if ymd eq '' then read,'What is the epoch [Julian date]? ',jd
read,'What is the length of the semimajor axis [km]? ',a
read,'What is the eccentricity? ',e
read,'What is the inclination [deg]? ',i
read,'What is the RA of the ascending node [deg]? ',o
read,'What is the argument of pericenter [deg]? ',w
read,'What is the mean anomaly [deg]? ',m
if ymd ne '' then begin ;Convert year-month-day to
year=fix(strmid(ymd,0,2)) ;Julian date.
day=double(strmid(ymd,4,2))
mon=fix(strmid(ymd,2,2))
day=day+dn([mon-1]) ;Calculate day number.
if mon gt 2 then begin ;Add a day if leap year and month
day=day+(year/4.0 eq fix(year/4.0)) ;is March or later.
if (fix(year/100.) eq year/100.) then day=day-(year/400. ne fix(year/400.))
endif ; mon
jd=43508.5d0+365.d0*(year-78.d0)+ $ ;Calculate Julian date.
1.d0*fix((year-77.d0)/4.d0+1.d-04)+day
endif ;ymd ne ''
return
end ;findelt.pro