procedure lickast (input) ######################################################################### # This is a procedure to modify images originally acquired at Lick # # so that the astrometric keywords are fully accurate. The problem # # is that the RA, DEC, and HA cards in a Lick image are produced from # # the data supplied by TELCO. Users of Lick telescopes know that # # TELCO is not the best of systems, it is setup and calibrated by # # humans, and the encoders on the telescopes have scale and hysteresis # # errors which cause the recorded position to drift during a night. # # # # Anyway, the bottom line is that if you are working with something # # that wants accurate values of HJD, you need a better position. # # # # In this script we take as input some accurate catalog values of # # RA, DEC, and the EPOCH of these. (Note that EPOCH is a misnomer from # # the early days of FITS and should be EQUINOX. This is combined with # # the (more likely to be correct) values of the DATE-OBS and TIME in # # the image header. The original values are saved under new keywords # # named RATELCO, DECTELCO, and the catalog values are inserted. # # The task noao.astutil.asttimes, along with the OBSERVAT keyword, # # is used to get an accurate value of the local sidereal time to put # # into the ST keyword and to calculate a better value of HA. # # Finally, the noao.astutil tasks setairmass and setjd are called to # # finish the insertion of the other keywords which are useful when # # reducing photometry, spectra, etc. # # # # As one final nice thing to do, the script also looks for a ccdtype # # keyword in each entry in the catalog file. This allows the procedure # # to add in the IRAF imagetyp card with the value specified in the # # fifth column of the object catalog. # # # # Steve Allen (sla@lick.ucsc.edu) 1993 August 30 # # # # This IRAF cl script is Copyright (c) 1993,1994 Steven L. Allen # ######################################################################### string input {prompt = 'Input images'} file catalog {prompt = 'Catalog of object positions'} string observat="" {prompt = 'Default Observatory ID'} #real dra=-50. {prompt = 'Default Right Ascension'} #real ddec=-100. {prompt = 'Default Declination'} #real dequinox=-10000. {prompt = 'Default Precession Equinox'} struct *ilist="" # file structure pointers for reading lists struct *catlist="" # file structure pointer for reading catalog struct line="" {length=160} # a line read from the object catalog file begin string l_input,l_cat # internal copies of input and output string ifile, ofile # temporary files containing expanded lists string in # loop variables for stepping thru file lists string valstr # a keyword value real ha, lmst # accurate value of HA and LMST int year, month, day # from the DATE-OBS header card real time # from the TIME header card int idx # index into DATE-OBS string string curobs # observatory for the current index struct title, otitle # string sra,sdec,seq,styp # strings read from the object catalog file real ra,dec,equinox int i,j bool incat # we make use of several noao.astutil tasks as a part of this script noao astutil # fix so script will work in the background cache ("asttimes") cache ("imgets") cache ("observatory") # Pick names for the temporary files # Note that your variable tmp$ had better point to a big enough place ifile = mktemp ("tmp$lickast") ofile = mktemp ("tmp$lickcat") # prompt the user for the query parameters if not on command line l_input = input l_cat = catalog ra=0;dec=0;equinox=0 # expand the input and output image templates # (i.e., the given names may have contained wildcard characters) sections (l_input, > ifile) otitle="" # set file pointer to look at our input list, and loop ilist = ifile while (fscan (ilist, in) != EOF) { ################################################################ # we want to find this object in the catalog and look up its # ra, dec, and precession equinox therein title="" hselect(in, 'title', yes) | scan(title) if (title == otitle) { # Same object as was in the previous image. We don't need # to search thru the catalog file again because # the variables are already set. } else { # (title != otitle) # Not the same object as was in the previous image. # Save this object name for next time. otitle = title j = strlen(title) # in case we do not find the object in the catalog incat = no # Rewind catalog file pointer in preparation for reading it. catlist = l_cat # Scan thru each entry in the catalog while (fscan (catlist, line) != EOF) { if (substr(line,1,j) == title) { # we have found a matching object name incat = yes #print('match ',line) i=stridx("'",line)+1 j=strlen(line) sra=""; sdec=""; seq=""; styp=""; valstr=substr(line,i,j) # extract the astrometrically accurate positions #print('valstr|',valstr,'|') i=fscan(valstr,sra,sdec,seq,styp) #print('i ',i,' sra ',sra,' sdec ',sdec,' seq ',seq,' styp|',styp,'|') # break out of the catalog scanning while loop break } } catlist = "" } # endif (title == otitle) ################################################################ if ( ! incat ) { print ('Object ', title, ' is not in catalog ', l_cat ) # do not replace any of the position cards in this image # proceed on to the next one next } else { # add it to the output list for setjd and setairmass print ( in, >> ofile ) } ################################################################ # we presume without checking that the image has RA, DEC, and HA # and we make backup copies to show what TELCO originally said # For reference, the Lick data acquisition queries TELCO for # RA, DEC, and HA at the end of every exposure just before readout. ######## valstr = "" hselect( in, 'ratelco', yes) | scan (valstr) if (valstr == "") { # We have just determined that RATELCO does not exist. We add it. imgets( in, 'ra') # Note an interesting bug in cl syntax. If you want to add # a sexagesimal constant to a header, you have to do this: # hedit(in, 'ratelco', "(str('21:51:59.0'))", add+, show-, ver-) # But if you want to add a sexagesimal variable to a header, # you must do it in two steps, like this: hedit(in, 'ratelco', 'foo' , add+ , show-, ver-) hedit(in, 'ratelco', ""//imgets.value//"" , update+, show-, ver-) } if (sra == "TELCO") { # leave the TELCO position untouched hselect (in, 'ra', yes) | words | scan(ra) } else if (fscan(sra, ra) == 1) { # the catalog has a legitimate value. We insert it. hedit(in, 'ra', ""//sra//"", update+, show-, ver-) } ######## valstr = "" hselect( in, 'dectelco', yes) | scan (valstr) if (valstr == "") { # We have just determined that DECTELCO does not exist. We add it. imgets( in, 'dec') hedit(in, 'dectelco' , 'foo' , add+, show-, ver-) hedit(in, 'dectelco' , ""//imgets.value//"", update+, show-, ver-) } if (fscan(sdec, dec) == 1) { # the catalog has a legitimate value. We insert it. hedit(in, 'dec', ""//sdec//"", update+, show-, ver-) } ######## valstr = "" hselect( in, 'hatelco', yes) | scan (valstr) if (valstr == "") { # We have just determined that HATELCO does not exist. We add it. imgets( in, 'ha') hedit(in, 'hatelco' , 'foo' , add+, show-, ver-) hedit(in, 'hatelco' , ""//imgets.value//"" , update+, show-, ver-) } # we don't have a value for HA to put in here (yet) ################################ # parse header cards to get inputs for asttimes # Lick stores the UT of the start of the observation in TIME imgets( in, 'time') time=real(imgets.value) hedit(in, 'ut', time, add+, show-, ver-) # Lick stores the UT date of the start of the observation in DATE-OBS imgets( in, 'date-obs') idx=stridx("01234567898",imgets.value) valstr=substr(imgets.value,idx,strlen(imgets.value)) day =int(substr(valstr,1,2)) month=int(substr(valstr,4,5)) year =int(substr(valstr,7,8)) ################################ # values of the observatory task need to be set here!!!!! valstr = "" hselect( in, 'observat', yes) | scan (valstr) if (valstr != "") { # any OBSERVAT card in the header always takes precedence curobs=valstr } else { # there was no OBSERVAT card in the image if (observat != "") { # the user specified one, so we use and insert it here curobs=observat hedit(in, 'observat', observat, add+, show-, ver-) } else { # the user said nothing, so we assume lick by default # but notice that this can screw up asttimes if you # happen to hit a month boundary. You really want to # have that OBSERVAT card in there pointing to # a site with timezone=0. curobs="lick" } } if (observatory.observatory != curobs) { observatory( "set", curobs) } ################################ # call asttimes to get accurate value of LMST # Notice that this is assuming the Lick convention where the TIME # card was given as UT, and it is correcting for whatever timezone # is actually contained in the observatory database. time=time-observatory.timezone if (time < 0) { time = time + 24. # this code fails at month start day = day - 1 } asttimes(year=year , month=month , day=day , time=time, observatory=")_.observatory", header=no, >> "dev$null" ) ha = asttimes.lmst - ra if (ha > 12.) ha = ha - 24. if (ha <= -12.) ha = ha + 24. ################################ # we now have a value for EPOCH either from user or from asttimes # (EPOCH is a misnomer from early FITS days and should be EQUINOX.) valstr = "" hselect( in, 'epoch', yes) | scan (valstr) if (valstr != "") { # There is a pre-existing EPOCH card in the image. # (EPOCH is a misnomer from early FITS days and should be EQUINOX.) # We presume that the existing card is correct, and leave it. } else if (fscan(seq, equinox) == 1) { # the catalog has a legitimate value. We insert it. #hedit(in, 'epoch', 'foo' , add+, show-, ver-) #hedit(in, 'epoch', ""//seq//"", update+, show-, ver-) hedit(in, 'epoch', equinox, add+, show-, ver-) } else { # we assume that the precession equinox is "of date" hedit(in, 'epoch', asttimes.epoch , add+, show-, ver-) } ################################ # the object catalog file may also specify a CCD image type card valstr = "" hselect( in, 'imagetyp', yes) | scan (valstr) if (valstr != "") { #print (' valstr was not blank') # There is a pre-existing IMAGETYP card in the image. # We presume that the existing card is correct, and leave it. } else if (styp == "zero") { # If the object catalog has an image type of "zero" then we # look at the exposure time card. i = 0 hselect (in, 'exposure', yes) | scan (i) print (' styp is zero, exp is ',i) # By arbitrary fiat we assert that 1 second exposures at Lick # may be "zero" frames, but that longer ones are "dark". if ( i > 1 ) styp = "dark" hedit(in, 'imagetyp', styp, add+, show-, ver-) # put it back, because we may not reread the catalog styp = "zero" } else if (styp != "") { #print (' vstyp was not blank') # The catalog has a legitimate value. We insert it. hedit(in, 'imagetyp', styp, add+, show-, ver-) } else { #print (' we are lostd') # We have no idea what the CCD image type should be. # So we do nothing. } ################################ # We need to put LMST into the ST keyword for noao.astutil.setairmass # Note that we used the Lick TIME card (which is begin time) directly, # so we get ST at the beginning of the exposure. hedit(in, 'st' , asttimes.lmst, add+, show-, ver-) ######## # finally, we know enough to add the accurate HA keyword hedit(in, 'ha' , ha , update+, show-, ver-) } # end while (fscan (ilist, in) != EOF) # call setjd to set accurate values of JD, HJD, LJD for obs midpoint # this uses OBSERVAT, DATE-OBS, TIME, EXPOSURE, RA, DEC, and EPOCH # (EPOCH is a misnomer from early FITS days and should be EQUINOX.) setjd( "@"//ofile, date="date-obs", time="time", exposure="exposure", ra="ra", dec="dec", epoch="epoch", jd="jd", hjd="hjd", ljd="ljd", utdate=yes, uttime=yes, listonly=no) # call setairmass to set an accurate value of AIRMASS and UTMIDDLE # this uses OBSERVAT, DATE-OBS, EXPOSURE, RA, DEC, EPOCH, ST and UT # (EPOCH is a misnomer from early FITS days and should be EQUINOX.) setairmass( "@"//ofile, intype="beginning", outtype="effective", date="date-obs", exposure="exposure", airmass="airmass", utmiddle="utmiddle", show-, update=yes, override=yes) # clear this list so we don't remember its content in the par file ilist = "" line = "" # get rid of the temporary files we created as workspace delete (ifile , ver-) delete (ofile , ver-) end