#!/usr/bin/python2 import sys import os import time import string from math import * import Numeric import MLab import Gnuplot from Tkinter import * import tkMessageBox g = Gnuplot.Gnuplot() g('set data style points') def myPrint(aString,logfile): print aString print >>logfile, aString return def getRootName(fileName): i = fileName.find('.') return fileName[0:i] def getRawFileNameFrom(daomasterMchFileName): i = daomasterMchFileName.find('.') newName = daomasterMchFileName[0:i] return newName + ".raw" def getTfrFileNameFrom(daomasterMchFileName): i = daomasterMchFileName.find('.') newName = daomasterMchFileName[0:i] return newName + ".tfr" def getMonsterFileNameFrom(daomasterMchFileName): i = daomasterMchFileName.find('.') newName = daomasterMchFileName[0:i] return newName + ".monster" def getMonFileNameFrom(daomasterMchFileName): # this is for the reduced version i = daomasterMchFileName.find('.') newName = daomasterMchFileName[0:i] return newName + ".mon" def getListOfApFileNamesFrom(mchFile): try: infile = open(mchFile,"rb") except: print "No file named %s exists" % fileName sys.exit(0) linesList = infile.readlines() infile.close nlines = len(linesList) cc = "\'" #print cc fileNameList = [] for aLine in linesList: i1 = aLine.find(cc) i2 = aLine.find(cc,i1+1) #print "i1 = %2d i2 = %2d" % (i1, i2,) fileName = aLine[i1+1:i2] fileName.strip() fileNameList.append(fileName) return(fileNameList) def getListOfImagesFrom(mchFile): imageNameList = getListOfApFileNamesFrom(mchFile) outputList = [] for image in imageNameList: i1 = image.find('_') i2 = image.find('_',i1+1) i3 = image.find('_',i2+1) i4 = image.find('_',i3+1) rootName = image[:i2] number = image[i2+1:i3] ext = image[i4-1:i4] outputList.append(rootName + '.' + number + '.' + 'fits' + '[' + ext + ']') return(tuple(outputList)) def convertApToAlsFileName(apFileName): i1 = apFileName.rfind('.') stripedName = apFileName[:i1] #print stripedName for i in range(1): i1 = stripedName.rfind('_') s1 = stripedName[:i1] s2 = stripedName[i1+1:] stripedName = s1 + '.' + s2 #print stripedName return stripedName def getBestMagnitude(m1,e1,sky1,sh1,chi1,m2,e2,sky2,sh2,chi2,m3,e3,sky3,sh3,chi3,badIdentifier): """getBestMagnitude(m1,e1,sky1,sh1,chi1,m2,e2,sky2,sh2,chi2,m3,e3,sky3,sh3,chi3,badIdentifier): Input: m1, m2, m3 : three magnitude values badIdentifier: number flagging bad mag values Output: the best magnitude estimate from those three values together with the associated values of msky, sharp, and chi Notice that the routine only checks the magnitude values for goodness. It assumes that if they are good, the associated merr, msky, sharp, and chi are also good. """ if good(m1) and good(m2) and good(m3): m = (m1 + m2 + m3 ) / 3.0 merr = sqrt(e1*e1 + e2*e2 + e3*e3) sky = (sky1 + sky2 + sky3) / 3.0 sharp = (sh1 + sh2 + sh3) / 3.0 chi = (chi1 + chi2 + chi3) / 3.0 dm = sqrt(((m1-m)**2 + (m2-m)**2 + (m3-m)**2)/2.0) n = 3 elif (not good(m1)) and good(m2) and good(m3): m = (m2 + m3) / 2.0 merr = sqrt(e2*e2 + e3*e3) sky = (sky2 + sky3) / 2.0 sharp = (sh2 + sh3) / 2.0 chi = (chi2 + chi3) / 2.0 dm = abs(m2-m3) n = 2 elif good(m1) and (not good(m2)) and good(m3): m = (m1 + m3) / 2.0 merr = sqrt(e1*e1 + e3*e3) sky = (sky1 + sky3) / 2.0 sharp = (sh1 + sh3) / 2.0 chi = (chi1 + chi3) / 2.0 dm = abs(m1-m3) n = 2 elif good(m1) and good(m2) and (not good(m3)): m = (m1 + m2) / 2.0 merr = sqrt(e1 + e2) sky = (sky1 + sky2) / 2.0 sharp = (sh1 + sh2) / 2.0 chi = (chi1 + chi2) / 2.0 dm = abs(m1-m2) n = 2 elif good(m1): m = m1 merr = e1 sky = sky1 sharp = sh1 chi = chi1 dm = 9.9999 n = 1 elif good(m2): m = m2 merr = e2 sky = sky2 sharp = sh2 chi = chi2 dm = 9.9999 n = 1 elif good(m3): m = m3 merr = e3 sky = sky3 sharp = sh3 chi = chi3 dm = 9.9999 n = 1 else: m = badIdentifier merr = 9.9999 sky = 99.9999 sharp = 99.9999 chi = 99.9999 dm = 9.9999 n = 0 return (m,merr,dm,sky,sharp,chi,n) def good(m): badDataIdentifier = 90.9999 if 0.0 < m and m < badDataIdentifier-1.0: return 1 else: return 0 def getMagnitudeOffset(mag1,mag2,dmagmax,mag3,mag3min,mag3max): # some checking of the input data # dmagmax = abs(dmagmax) if mag3min > mag3max: aux = mag3min mag3min = mag3max mag3max = aux # we loop creating arrays only with good data m1 = [] dm12 = [] m3 = [] for i in range(len(mag1)): if good(mag1[i]) and good(mag2[i]) and good(mag3[i]): m1.append(mag1[i]) dm12.append(mag1[i]-mag2[i]) m3.append(mag3[i]) if len(m1) > 0: m1 = Numeric.array(m1) dm12 = Numeric.array(dm12) m3 = Numeric.array(m3) else: print "one of the array into getMagnitudeOffsets had no good data" sys.exit(0) # we get a rough initial guess of the statistics npy,miny,maxy,meany,stdey,sigy,mediany = getStatsY(dm12,-1.9,1.9,m3,mag3min,mag3max) # we refine the guess here dmmin = mediany - 0.7 dmmax = mediany + 0.7 npy,miny,maxy,meany,stdey,sigy,mediany = getStatsY(dm12,dmmin,dmmax,m3,mag3min,mag3max) # and we go here for the final answer dmmin = mediany - dmagmax dmmax = mediany + dmagmax npy,miny,maxy,meany,stdey,sigy,mediany = getStatsY(dm12,dmmin,dmmax,m3,mag3min,mag3max) return(npy,miny,maxy,meany,stdey,sigy,mediany) def readAlsFile(fileName): """readAlsFile(fileName): Input: name of file containing allstar photometry file Output: id dictWithDataList This function reads an allstar file skipping the lines with comments. Then it reads the data, which it returns as a tuple containing a list with the ids, and a dictionary with lists of xc, yc, mag, merr, niter, sharp, chi, pier, and perror, indexed with id. id is a list of integers. """ try: infile = open(fileName,"rb") except: print "No file named %s exists" % fileName sys.exit(0) linesList = infile.readlines() infile.close nlines = len(linesList) # we first determine the number of comment lines # for counter in xrange(0,nlines-1,1): currentLine = linesList[counter] if currentLine[0] != '#': break numberOfCommentLines = counter # now we loop over the data lines extracting usefull fields # id = {} indice = {} dataListDict = {} starNumber = 0 for counter in xrange(numberOfCommentLines,nlines-1,2): starNumber = starNumber + 1 t1 = linesList[counter] t1 = t1[:-2] # getting rid of newline and \ characters t2 = linesList[counter+1] t2 = t2[:-1] t = t1 + t2 dataList = string.split(t) if dataList[3] == "INDEF": continue # we skip data which does not have good values floatList = map(string.atof,dataList[1:-1]) id[starNumber] = string.atoi(dataList[0]) dataListDict[string.atoi(dataList[0])] = floatList[0:9] + string.split(dataList[10]) indice[string.atoi(dataList[0])] = starNumber return(id,indice,dataListDict) def readDaoRaw(fileName): """readDaoRaw(fileName): Input: name of file containing daomaster output raw file Output: id dictWithDataList This function reads a daomaster raw file skipping the the first three lines. Then it reads the data, which it returns as a tuple containing a list with the ids, and a dictionary with lists of xc, yc, v1, ev1 v2, ev2, v3, ev3, vv, nvv, b1, eb1, b2, eb2, b3, eb3, bb, nbb, u1, eu1, u2, eu2, u3, eu3, uu, nuu, indexed with id. id is a list of integers. """ try: infile = open(fileName,"rb") except: print "No file named %s exists" % fileName sys.exit(0) linesList = infile.readlines() infile.close nlines = len(linesList) id = [] dataListDict = {} numberOfCommentLines = 3 for counter in xrange(numberOfCommentLines,nlines-1,2): t1 = linesList[counter] t1 = t1[:-1] # getting rid of newline character t2 = linesList[counter+1] t = t1 + t2 if t == '': break # here we read the data from the line # dataList = string.split(t) floatList = map(string.atof,dataList[1:]) id.append(string.atoi(dataList[0])) dataListDict[string.atoi(dataList[0])] = floatList #print floatList #sys.exit(0) return(id,dataListDict) def readDaoTfr(fileName): """readDaoTfr(fileName): Input: name of file containing daomaster output transfer table Output: n1, n2, n3, n4, n5, n6, n7, n8, n9 (dictionaries) This program returns nm(n1), nm(n2), etc, as associative arrays (Dictionaries). Each dictionary contains the id of the corresponding star in the master list. """ try: infile = open(fileName,"rb") except: print "No file named %s exists" % fileName sys.exit(0) linesList = infile.readlines() infile.close nlines = len(linesList) nm_1, nm_2, nm_3, nm_4, nm_5, nm_6, nm_7, nm_8, nm_9 = {}, {}, {}, {}, {}, {}, {}, {}, {} numberOfLinesToSkip = 10 # THIS SHOULD BE MADE MORE ROBUST for counter in xrange(numberOfLinesToSkip,nlines-1,1): t = linesList[counter] t = t[:-1] # getting rid of newline character if t == '': break # here we read the data from the line # dataList = string.split(t) intList = map(string.atoi,dataList[3:]) masterId = string.atoi(dataList[0]) nm_1[intList[0]] = masterId nm_2[intList[1]] = masterId nm_3[intList[2]] = masterId nm_4[intList[3]] = masterId nm_5[intList[4]] = masterId nm_6[intList[5]] = masterId nm_7[intList[6]] = masterId nm_8[intList[7]] = masterId nm_9[intList[8]] = masterId return(nm_1,nm_2,nm_3,nm_4,nm_5,nm_6,nm_7,nm_8,nm_9) def getTfrDict(tfrFile,column): try: tfrfile = open(tfrFile,"rb") except: print "problem opening %s" % (tfrFile,) lines = tfrfile.readlines() tfrfile.close() # now we fin the index for the first data line found = 0 counter = 0 while not found: if lines[counter][1] == '=': found = 1 counter = counter + 1 firstDataLine = counter # now we loop over the lines creating the dictionary # starNumber = {} for aLine in lines[firstDataLine:]: aLine = string.split(aLine[:-1]) starNumber[string.atoi(aLine[0])] = string.atof(aLine[column-1]) return starNumber def getValueFromTable(tableName,tableColumn,idArray,starNumberDict,badData): try: tablefile = open(tableName,'rb') except: print "problem reading %s" % (tableName,) lines = tablefile.readlines() tablefile.close() # we first create a dictionary with the values in # tableColumn indexed by row number columnValue = {} rownumber = 0 for aLine in lines: rownumber = rownumber + 1 aLine = string.split(aLine[:-1]) columnValue[rownumber] = string.atof(aLine[tableColumn-1]) # we the fill the new Float64 array with this values when # available, or badData otherwise # newArray = idArray*0.0 for i in range(len(idArray)): if starNumberDict.has_key(idArray[i]): newArray[i] = columnValue[starNumberDict[idArray[i]]] else: newArray[i] = badData return newArray def loadAlsArrays(indx, nm, alsFileName, xc, yc, mag, merr, msky, sharp, chi): id, starNumber, dataDict = readAlsFile(alsFileName) NNN = 0 for i in id.keys(): if nm.has_key(starNumber[id[i]]): if indx.has_key(nm[starNumber[id[i]]]): NNN = NNN + 1 zkey = indx[nm[starNumber[id[i]]]] xc[zkey] = dataDict[id[i]][0] yc[zkey] = dataDict[id[i]][1] mag[zkey] = dataDict[id[i]][2] merr[zkey] = dataDict[id[i]][3] msky[zkey] = dataDict[id[i]][4] sharp[zkey] = dataDict[id[i]][6] chi[zkey] = dataDict[id[i]][7] print "%6d stars processed - done" % (NNN,) del id del starNumber del dataDict return def plotit(x,x1,x2,y,y1,y2,*options): global g if len(options) == 0: d = Gnuplot.Data(x,y) elif len(options) == 1: d = Gnuplot.Data(x,y,with=options) else: print('inappropriate number of optional arguments') return if len(options) != 0: print "printing with options = %s" % options if x1 != x2: rangeXString = 'set xrange [' + str(x1) + ':' + str(x2) + ']' g(rangeXString) if y1 != y2: rangeYString = 'set yrange [' + str(y1) + ':' + str(y2) + ']' g(rangeYString) g.plot(d) return def getStatsY(y,ymin,ymax,x,xmin,xmax): filty = [] for i in range(len(y)): if ymin < y[i] and y[i] < ymax and xmin < x[i] and x[i] < xmax: filty.append(y[i]) filty = Numeric.array(filty) meany = MLab.mean(filty) mediany = MLab.median(filty) sigy = MLab.std(filty) miny = MLab.min(filty) maxy = MLab.max(filty) npointsy = len(filty) return(npointsy,miny,maxy,meany,sigy/sqrt(npointsy),sigy,mediany) class Monster: "Class to read and handle monster photometry files" def __init__(self,fileName): "Initialize the class by reading the data lists" self.infile = fileName print "I read from %s" % self.infile try: infile = open(fileName, "rb") except: print "Problem reading %s" % self.infile sys.exit(0) linesList = infile.readlines() infile.close self.nlines = len(linesList) # now we parse the first line # skipping the first gato char # self.headerString = linesList[0][1:-1] self.headerList = string.split(self.headerString) self.formatString = linesList[1][1:-1] self.formatList = string.split(self.formatString) if len(self.headerList) != len(self.formatList): print "incompatible length of data list and format list" return self.numberOfColumns = len(self.headerList) self.firstDataRow = 2 # first row that could contain data print "%6d columns: %s" % (self.numberOfColumns, self.headerString,) print "%6d columns: %s" % (self.numberOfColumns, self.formatString,) # now we read the data in each of the columns into our column list # # first we mapp column titles to column number (starting at 0) # counter = 0 self.columnNumber = {} for columnTitle in self.headerList: self.columnNumber[columnTitle] = counter counter = counter + 1 # now we fill the columns one row at a time # self.dataInColumn = {} for title in self.headerList: self.dataInColumn[title] = [] counter = 0 for aLine in linesList[self.firstDataRow:]: aLine = aLine[:-1] # remove ending new line char #print "aLine: %s" % (aLine,) t = string.split(aLine) newList = [] for jj in range(len(t)): if (self.formatList[jj]).find('d') >= 0: newList.append(string.atoi(t[jj])) elif (self.formatList[jj]).find('f') >= 0: newList.append(string.atof(t[jj])) t = newList if len(t) == 0: continue #print "length of t is %6d" % len(t) for title in self.headerList: (self.dataInColumn[title]).append(t[self.columnNumber[title]]) counter = counter + 1 for title in self.headerList: self.dataInColumn[title] = Numeric.array(self.dataInColumn[title]) self.numberOfDataRows = counter #print self.dataInColumn print "%6d data rows loaded" % counter del linesList return def save(self,fileName): #print 'procedure not ready' #return try: outfile = open(fileName, "wb") except: print "Problem opening %s" % fileName sys.exit(0) print >>outfile, "# %s" % self.headerString print >>outfile, "# %s" % self.formatString print self.formatList for counter in xrange(0,self.numberOfDataRows): newLine = "" dataList = [] for title in self.headerList: dataList.append(self.dataInColumn[title][counter]) #print title,self.columnNumber[title],self.formatList[self.columnNumber[title]],self.dataInColumn[title][counter] newLine = newLine + (self.formatList[self.columnNumber[title]] % self.dataInColumn[title][counter]) + " " print >>outfile, newLine outfile.close() return def show(self): n1 = len(string.split(self.headerString)) n2 = len(string.split(self.formatString)) print "%4d columns: %s" % (n1, self.headerString,) print "%4d columns: %s" % (n2, self.formatString,) return def append(self,colTitle,newArray,formatString): """append a new list as a new column. It checks that the length of the list equals the table length """ if len(newArray) != self.numberOfDataRows: print "length of list (%d) differ from table (%s)" % (len(newArray), self.numberOfDataRows,) return if formatString == "": formatString = "%12.4e" self.headerList = self.headerList + string.split(colTitle) self.formatList = self.formatList + string.split(formatString) self.headerString = self.headerString + " " + colTitle self.formatString = self.formatString + " " + formatString self.numberOfColumns = len(self.headerList) self.columnNumber[colTitle] = self.numberOfColumns - 1 # columns start at 0 print "new data in column %3d" % self.columnNumber[colTitle] self.dataInColumn[colTitle] = newArray print "%6d columns: %s" % (self.numberOfColumns, self.headerString,) print "%6d columns: %s" % (self.numberOfColumns, self.formatString,) return def overwrite(self,colTitle,newArray): # overwrites a new list over an existing column. # It checks that the length of the list # equals the table length # if len(newArray) != self.numberOfDataRows: print "length of list (%d) differ from table (%s)" % (len(newArray), self.numberOfDataRows,) return self.headerString = self.headerString self.numberOfColumns = self.numberOfColumns + 1 self.dataInColumn[colTitle] = newArray return def column(self,columnTitle): #print "returning data in column named %s" % columnTitle return(self.dataInColumn[columnTitle]) def hardcopy(self): g.hardcopy('monster.ps', enhanced=1, color=1) print ('\n******** Saved plot to postscript file "monster.ps" ********\n') return def plot(self,xcolTitle,x1,x2,ycolTitle,y1,y2,*options): global g x = self.column(xcolTitle) y = self.column(ycolTitle) if len(options) == 0: d = Gnuplot.Data(x,y) elif len(options) == 1: d = Gnuplot.Data(x,y,with=options) else: print('inappropriate number of optional arguments') return if len(options) != 0: print "printing with options = %s" % options #g = Gnuplot.Gnuplot() if x1 != x2: rangeXString = 'set xrange [' + str(x1) + ':' + str(x2) + ']' g(rangeXString) if y1 != y2: rangeYString = 'set yrange [' + str(y1) + ':' + str(y2) + ']' g(rangeYString) g.xlabel(xcolTitle) g.ylabel(ycolTitle) g.plot(d) #raw_input("Press any key to continue") return def tvmark(self,xcol,ycol): x = self.column(xcol) y = self.column(ycol) try: tmpfile = open("tvmark.tmp","wb") except: print "problem opening tmp file" return linesList = "" for i in xrange(len(x)): aLine = str(x[i]) + " " + str(y[i]) + "\n" linesList = linesList + aLine print >>tmpfile, linesList tmpfile.close() os.system("shiraf_tvmark 1 tvmark.tmp col=204 points=1") #os.system("rm tvmark.tmp") return def ds9_raise(): # we raise the display # os.system("xpaset -p ds9 raise") return def ds9_centerOn(xcen,ycen): # we center on the coordinates given # os.system("xpaset -p ds9 pan to " + str(xcen) + " " + str(ycen)) return def ds9_setZoomTo(zoom): # we zoom into this region # os.system("xpaset -p ds9 zoom to " + str(zoom)) def ds9_getNewXYAfterClick(): # we move to an unlikely place for the crosshair to be # os.system("xpaset -p ds9 crosshair -157 -132 physical") # we read the current position # comPipe = os.popen("xpaget ds9 crosshair") line = string.split(comPipe.readline()) comPipe.close() xPrev = line[0] yPrev = line[1] while 1: comPipe = os.popen("xpaget ds9 crosshair") line = string.split(comPipe.readline()) comPipe.close() xNew = line[0] yNew = line[1] if xNew != xPrev or yNew != yPrev: return(string.atof(xNew),string.atof(yNew)) time.sleep(0.05) def ds9_get2ClicksOffset(): (x1,y1) = ds9_getNewXYAfterClick() (x2,y2) = ds9_getNewXYAfterClick() return(x1-x2,y1-y2) def ds9_getOffsetInteractor(xc,yc): done = 0 while not done: ds9_raise() ds9_centerOn(xc,yc) tkMessageBox.showinfo("getCoordinates","After clicking OK on this window, " + \ "click on the star you will use " + \ "to calculate offset and then " + \ "click on the corresponding dot") (xOffset,yOffset) = ds9_get2ClicksOffset() print "applying offsets of %7.1f, %7.1f to coordinate list from frame 1" % \ (xOffset, yOffset,) message = "dx = " + str(xOffset) + " dy = " + str(yOffset) + "\n\n" + \ "Are these offsets OK?\n Click No to repeat" done = tkMessageBox.askyesno("Confirm",message) return(xOffset,yOffset) def Disp(frameToDisp,z1,z2): ds9_raise() comando = "shiraf_disp " + frameToDisp + " 1 zscale- zrange- z1=" + str(z1) + " z2=" + str(z2) print "Disp: %s" % comando os.system(comando) return def Find(frameToFind,outCoo,datapars,findpars): comando = "shiraf_daofind image=" + frameToFind + " output=" + outCoo + \ " datapars=" + datapars + " findpars=" + findpars + \ " skymap=" + "\\\"\\\"" + " starmap=" + "\\\"\\\"" + \ " boundary=\"nearest\"" + " interactive=no" print " -->%s" % comando os.system(comando) return def FindInteractive(frameToFind,outCoo,datapars,findpars): comando = "shiraf_daofind image=" + frameToFind + " output=" + outCoo + \ " datapars=" + datapars + " findpars=" + findpars + \ " skymap=" + "\\\"\\\"" + " starmap=" + "\\\"\\\"" + \ " boundary=\"nearest\"" + " interactive=yes" print " -->%s" % comando os.system(comando) return def Tvmark(cooFile,color,psize): ds9_raise() comando = "shiraf_tvmark 1 " + cooFile + " col=" + str(color) + " points=" + str(psize) os.system(comando) return def Phot(image,cooFile,outPhotFile,datapars,centerpars,fitskypars,photpars): comando = "shiraf_phot image=" + image + " coords=" + cooFile + \ " output=" + outPhotFile + \ " skyfile=" + "\\\"\\\"" + " plotfile=" + "\\\"\\\"" + \ " datapars=" + datapars + " centerpars=" + centerpars + \ " fitskypars=" + fitskypars + " photpars=" + photpars + \ " interactive=no radplots=no" os.system(comando) return def PhotInteractive(image,cooFile,outPhotFile,datapars,centerpars,fitskypars,photpars): comando = "shiraf_phot image=" + image + " coords=" + cooFile + \ " output=" + outPhotFile + \ " skyfile=" + "\\\"\\\"" + " plotfile=" + "\\\"\\\"" + \ " datapars=" + datapars + " centerpars=" + centerpars + \ " fitskypars=" + fitskypars + " photpars=" + photpars + \ " interactive=yes radplots=no" os.system(comando) return def Imcopy(inMEFImage,outFITSImage,chipN,irafSection): comando = "shiraf_imcopy input=" + inMEFImage + "[" + chipN + "]" + \ irafSection + " output=" + outFITSImage os.system(comando) return def writeDatapars(sigma): lines = """scale,r,h,1.,0.,,\"Image scale in units per pixel\" fwhmpsf,r,h,3.,0.,,\"FWHM of the PSF in scale units\" emission,b,h,yes,,,\"Features are positive ?\" sigma,r,h,""" + str(sigma) + """,0.,,\"Standard deviation of background in counts\" datamin,r,h,-30.,,,\"Minimum good data value\" datamax,r,h,60000.,,,\"Maximum good data value\" noise,s,h,\"poisson\",|poisson|,,\"Noise model\" ccdread,s,h,\"\",,,\"CCD readout noise image header keyword\" gain,s,h,\"OUTCONAD\",,,\"CCD gain image header keyword\" readnoise,r,h,4.,,,\"CCD readout noise in electrons\" epadu,r,h,0.,,,\"Gain in electrons per count\" exposure,s,h,\"EXPTIME\",,,\"Exposure time image header keyword\" airmass,s,h,\"AIRMSTAR\",,,\"Airmass image header keyword\" filter,s,h,\"FILT1NAM\",,,\"Filter image header keyword\" obstime,s,h,\"UTC\",,,\"Time of observation image header keyword\" itime,r,h,1.,,,\"Exposure time\" xairmass,r,h,INDEF,,,\"Airmass\" ifilter,s,h,INDEF,,,\"Filter\" otime,s,h,INDEF,,,\"Time of observation\" mode,s,h,\"ql\",,,""" datapfile = open (dataparsFileName(),"wb") print >>datapfile, lines datapfile.close() return def delDatapars(): os.system("rm " + dataparsFileName()) return def dataparsFileName(): return "zpdatapars.par" def writeFindpars(sigma): lines = "threshold,r,h," + str(sigma) + """,,,\"Threshold in sigma for feature detection\" nsigma,r,h,1.5,,,\"Width of convolution kernel in sigma\" ratio,r,h,1.,0.,1.,\"Ratio of minor to major axis of Gaussian kernel\" theta,r,h,0.,0.,180.,\"Position angle of major axis of Gaussian kernel\" sharplo,r,h,0.2,,,\"Lower bound on sharpness for feature detection \" sharphi,r,h,1.,,,\"Upper bound on sharpness for feature detection\" roundlo,r,h,-1.,,,\"Lower bound on roundness for feature detection\" roundhi,r,h,1.,,,\"Upper bound on roundness for feature detection\" mkdetections,b,h,no,,,\"Mark detections on the image display ?\" mode,s,h,\"ql\",,,""" findpfile = open(findparsFileName(),"wb") print >>findpfile, lines findpfile.close() return def delFindpars(): os.system("rm " + findparsFileName()) return def findparsFileName(): return "zpfindpars.par" def writeCenterpars(): lines = """calgorithm,s,h,\"centroid\",|centroid|gauss|none|ofilter|,,\"Centering algorithm\" cbox,r,h,7.,,,\"Centering box width in scale units\" cthreshold,r,h,0.,,,\"Centering threshold in sigma above background\" minsnratio,r,h,1.,0.,,\"Minimum signal-to-noise ratio for centering algorithm\" cmaxiter,i,h,10,,,\"Maximum number of iterations for centering algorithm\" maxshift,r,h,7.,,,\"Maximum center shift in scale units\" clean,b,h,no,,,\"Symmetry clean before centering ?\" rclean,r,h,1.,,,\"Cleaning radius in scale units\" rclip,r,h,2.,,,\"Clipping radius in scale units\" kclean,r,h,3.,,,\"Rejection limit in sigma\" mkcenter,b,h,yes,,,\"Mark the computed center on display ?\" mode,s,h,\"ql\",,,""" centpfile = open(centerparsFileName(),"wb") print >>centpfile, lines centpfile.close() return def delCenterpars(): os.system("rm " + centerparsFileName()) return def centerparsFileName(): return "zpcenterpars.par" def writeFitskypars(alg,innerRadius,width): lines = "salgorithm,s,h," + alg + """,|median|mode|centroid|gauss|crosscor|ofilter|histplot|radplot|constant|file|mean|,,\"Sky fitting algorithm\" annulus,r,h,""" + str(innerRadius) + """,,,\"Inner radius of sky annulus in scale units\" dannulus,r,h,""" + str(width) + """,,,\"Width of sky annulus in scale units\" skyvalue,r,h,0.,,,\"User sky value\" smaxiter,i,h,10,,,\"Maximum number of sky fitting iterations\" sloclip,r,h,0.,,,\"Lower clipping factor in percent\" shiclip,r,h,0.,,,\"Upper clipping factor in percent\" snreject,i,h,50,,,\"Maximum number of sky fitting rejection iterations\" sloreject,r,h,3.,,,\"Lower K-sigma rejection limit in sky sigma\" shireject,r,h,3.,,,\"Upper K-sigma rejection limit in sky sigma\" khist,r,h,3.,,,\"Half width of histogram in sky sigma\" binsize,r,h,0.1,,,\"Binsize of histogram in sky sigma\" smooth,b,h,no,,,\"Boxcar smooth the histogram\" rgrow,r,h,0.,,,\"Region growing radius in scale units\" mksky,b,h,no,,,\"Mark sky annuli on the display\" mode,s,h,\"ql\",,,""" fitsfile = open(fitskyparsFileName(),"wb") print >>fitsfile, lines fitsfile.close() return def delFitskypars(): os.system("rm " + fitskyparsFileName()) return def fitskyparsFileName(): return "zpfitskypars.par" def writePhotpars(apertureRadius): lines = """weighting,s,h,\"constant\",|constant|cone|gauss|,,\"Photometric weighting scheme for wphot\" apertures,s,h,\"""" + str(apertureRadius) + """\",,,\"List of aperture radii in scale units\" zmag,r,h,25.,,,\"Zero point of magnitude scale\" mkapert,b,h,no,,,\"Draw apertures on the display\" mode,s,h,\"ql\",,,""" photfile = open(photparsFileName(),"wb") print >>photfile, lines photfile.close() return def delPhotpars(): os.system("rm " + photparsFileName()) return def photparsFileName(): return "zpphotpars.par" def offsetCoo(cooIn,xOffset,yOffset,cooOut): infile = open(cooIn,"rb") outfile = open(cooOut,'wb') allLines = infile.readlines() for aLine in allLines: if aLine[0] != '#' and len(aLine) != 0: aLine = aLine[:-1] aLine = string.split(aLine) newX = string.atof(aLine[0]) + xOffset newY = string.atof(aLine[1]) + yOffset aLine[0] = "%10.3f" % newX aLine[1] = "%10.3f" % newY outLine = "" for word in aLine: outLine = outLine + " " + word print >>outfile, outLine outfile.close() infile.close() del allLines del aLine return def writeCoo(x,y,mag,id,cooOut): outfile = open(cooOut,'wb') for i in range(len(x)): print >>outfile, "%10.3f %10.3f %8.4f %6d" % (x[i], y[i], mag[i], id[i],) outfile.close() return def readMagFromPhot(photFile): # we start reading the lines from photFile comando = "shiraf_pdump " + photFile + " mag yes" comPipe = os.popen(comando) lines = comPipe.readlines() comPipe.close() # here we create the list to return the mags in mags = [] for j in range(len(lines)): line = string.split(lines[j][:-1]) if line[0] != 'INDEF': mags.append(string.atof(line[0])) else: mags.append(99.9999) # WE MUST RETURN THE WHOLE LIST return(mags) def calcPhotDiff(photFile1,photFile2,x1_Origin,y1_Origin,x2_Origin,y2_Origin): # we start by reading the lines from the two # phot files # comando1 = "shiraf_pdump " + photFile1 + " id,xc,yc,mag,merr yes" comPipe1 = os.popen(comando1) lines1 = comPipe1.readlines() comPipe1.close() comando2 = "shiraf_pdump " + photFile2 + " id,xc,yc,mag,merr yes" comPipe2 = os.popen(comando2) lines2 = comPipe2.readlines() comPipe2.close() print "lines1: %6d lines2: %6d " % (len(lines1),len(lines2),) # here we create the lists that we will return # xc1, yc1, xc2, yc2, mavg, dm21, merr, id = [], [], [], [], [], [], [], [] # now we extract the appropriate fields from the # lines and append them to the corresponding lists # if len(lines1) == len(lines2): for j in xrange(len(lines1)): line1 = lines1[j][:-1] line2 = lines2[j][:-1] # now we apply the appropriate chip offsets line1 = string.split(line1) line2 = string.split(line2) mag1 = line1[3] merr1 = line1[4] mag2 = line2[3] merr2 = line2[4] if mag1 != "INDEF" and merr1 != "INDEF" and mag2 != "INDEF" and merr2 != "INDEF": id1 = string.atoi(line1[0]) id2 = string.atoi(line2[0]) if id1 != id2: print "ids did not match" sys.exit(0) mag1 = string.atof(mag1) merr1 = string.atof(merr1) newX1 = string.atof(line1[1]) + x1_Origin newY1 = string.atof(line1[2]) + y1_Origin mag2 = string.atof(mag2) merr2 = string.atof(merr2) newX2 = string.atof(line2[1]) + x2_Origin newY2 = string.atof(line2[2]) + y2_Origin avgmag = (mag2 + mag1)/2.0 dm = mag2 - mag1 magerr = sqrt(merr1*merr1 + merr2*merr2) xc1.append(newX1) yc1.append(newY1) xc2.append(newX2) yc2.append(newY2) mavg.append(avgmag) dm21.append(dm) merr.append(magerr) id.append(id1) else: print "lists have different lengths" sys.exit(0) return(xc1,yc1,xc2,yc2,dm21,merr,mavg,id) def writePhotDiffTo(zpFileName,xc1,yc1,xc2,yc2,dm21,merr,mavg,id,chipN1,chipN2): if zpFileName == 'stdout': for i in xrange(len(xc1)): outLine = "%7.1f %7.1f %7.1f %7.1f %7.3f %6.3f %7.3f %7d" % \ (xc1[i], yc1[i], xc2[i], yc2[i], dm21[i], \ merr[i], mavg[i], chipN1*10000000 + chipN2*1000000 + id[i],) print outLine del outLine else: zpfile = open(zpFileName,'ab') for i in xrange(len(xc1)): outLine = "%7.1f %7.1f %7.1f %7.1f %7.3f %6.3f %7.3f %7d" % \ (xc1[i], yc1[i], xc2[i], yc2[i], dm21[i], \ merr[i], mavg[i], chipN1*10000000 + chipN2*1000000 + id[i],) print >>zpfile, outLine del outLine zpfile.close() return def readPhotDiffFrom(zpFileName): zpfile = open(zpFileName,'rb') lines = zpfile.readlines() zpfile.close() xc1,yc1,xc2,yc2,dm21,merr,mavg,id = [],[],[],[],[],[],[],[] for aLine in lines: aLine = aLine[:-1] # strip it aLine = string.split(aLine) # split it xc1.append(string.atof(aLine[0])) yc1.append(string.atof(aLine[1])) xc2.append(string.atof(aLine[2])) yc2.append(string.atof(aLine[3])) dm21.append(string.atof(aLine[4])) merr.append(string.atof(aLine[5])) mavg.append(string.atof(aLine[6])) id.append(string.atof(aLine[7])) del aLine del lines xc1 = Numeric.array(xc1, Numeric.Float64) yc1 = Numeric.array(yc1, Numeric.Float64) xc2 = Numeric.array(xc2, Numeric.Float64) yc2 = Numeric.array(yc2, Numeric.Float64) dm21 = Numeric.array(dm21, Numeric.Float64) merr = Numeric.array(merr, Numeric.Float64) mavg = Numeric.array(mavg, Numeric.Float64) id = Numeric.array(id, Numeric.Int32) return(xc1,yc1,xc2,yc2,dm21,merr,mavg,id)