##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009 #! /usr/bin/env python # #ColorMagDiagram.py # #PURPOSE: #Generate a color magnitude diagram for the different filters used. If #requested, the obseravtions will also be matched to the measurements #in a standard catalog. # #KEYWORDS: # TvFlag = 0 - Don't plot or display # 1 - Display and offer user imexamine # OverRide = 0 - Don't overwrite the existing .coo or .mag files # 1 - Overwrite the files # OverWrite = 0 - Read the paramters from the existing .par files # 1 - Rewrite the parameter files with the values in this file # SkipPSF = 1 - Skip the PSF fitting # SkipNStar = 1 - Skip the NStar # SkipDAOFind = 1 - Skip the DAOFind routine # SkipAllStar = 1 - Skip the AllStar routine # SkipPSF = 1 - Skip the first psf fitting # SkipPSF2 = 1 - Skip the fitting of the psf after subtraction of neighbors # Comment = ""- String describing the conditions used for the psf fit # #CALLING SEQUENCE: #ColorMagDiagram, GFitsFileName, IFitsFileName, YFitsFileName, [Keywords] # #NOTE: If there are less than 3 filters, set the filenames to '' # #Example Calling Sequence: #run ~/python/H4RGAnly/ColorMagDiagram.py '/nfs/slac/g/ki/ki09/lances/Reduced/07Apr26/M13/M13_G_SlopeFlatFieldedCenteredDitheredMeanWCS.fits' '/nfs/slac/g/ki/ki09/lances/Reduced/07Apr26/M13/M13_I_SlopeFlatFieldedCenteredDitheredMeanWCS.fits' '/nfs/slac/g/ki/ki09/lances/Reduced/07Apr26/M13/M13_Y_SlopeFlatFieldedCenteredDitheredMeanWCS.fits' --Catalog='/nfs/slac/g/ki/ki09/lances/Reduced/07Apr26/M13/Astrometry/M13UBVPhotometryAndProperMotionsKadla1966.txt' --TvFlag=1 import sys,os,getopt,re,string from pyraf import iraf import pyfits from pylab import * import numpy from numpy import * import numdisplay import matplotlib.pylab as mplot import TableIO #Set Up for plotting PlotFile='/afs/slac/u/ki/lances/python/SPIEPlotSettings.py' if os.path.exists(PlotFile) : execfile(PlotFile) #The filename will be specified as the first argument. These files have #the form [Object Name]_[Filter Used]_SlopeFlatFielded...fits UseG=1 ; UseI=1 ; UseY=1 GFitsFileName = sys.argv[1] if GFitsFileName == '' : UseG = 0 IFitsFileName = sys.argv[2] if IFitsFileName == '' : UseI = 0 YFitsFileName = sys.argv[3] if YFitsFileName == '' : UseY = 0 #Offsets between images GIXOffset = 10 #I X Coord - GIXOffset = G X Coord GIYOffset = 0 #I Y Coord - GIYOffset = G Y Coord GYXOffset = -20 #Y X Coord - GYXOffset = G X Coord GYYOffset = -14 #Y Y Coord - GYXOffset = G Y Coord IYXOffset = GYXOffset - GIXOffset IYYOffset = GYYOffset - GIYOffset #AIRMASSES GAirMass = 1.0039087 IAirMass = 1.0110843 YAirMass = 1.0150339 #Magntiude Zeropoints GZero = 0.25 IZero = 0.62 YZero = 2.173 #Threshold for matching pixels PixTh = 6 PltSc = 0.12607 GIMatchNum = 0 GYMatchNum = 0 IYMatchNum = 0 #Now parse the command line keyword arguments keywords = ['TvFlag=', 'SkipDAOFind=', 'SkipDAOPhot=', 'SkipPSF=', \ 'SkipNStar=', 'SkipSubStar=','SkipAllStar=','OverRide=', \ 'SkipPSF2=', 'OverWrite=', 'Comment=', 'Catalog=' ] OverWrite = 0 OverRide = 0 TvFlag = 0 SkipDAOFind = 0 SkipDAOPhot = 0 SkipPSF = 0 SkipPSF2 = 0 SkipNStar = 0 SkipSubStar = 0 SkipAllStar = 0 Catalog = '' Comment = '' opts, extraparams = getopt.getopt(sys.argv[4:],'',keywords) for o,p in opts: if o in ['--TvFlag']: TvFlag=int(p) elif o in ['--OverRide']: OverRide=int(p) elif o in ['--SkipDAOFind']: SkipDAOFind = int(p) elif o in ['--SkipDAOPhot']: SkipDAOPhot = int(p) elif o in ['--SkipPSF']: SkipPSF = int(p) elif o in ['--SkipNStar']: SkipNStar = int(p) elif o in ['--SkipSubStar']: SkipSubStar = int(p) elif o in ['--SkipAllStar']: SkipAllStar = int (p) elif o in ['--SkipPSF2']: SkipPSF2 = int(p) elif o in ['--OverWrite']: OverWrite = int(p) elif o in ['--Comment']: Comment = p elif o in ['--Catalog']: Catalog = p #To use the Catalog, run CatalogToTable to extract the physical coordinates #of the stars in the Catalog for each image in the '...XY.tab' files if Catalog != '' : BaseName = os.path.basename(Catalog) DirName = os.path.dirname(Catalog)+'/' BaseParts = BaseName.split('.') BaseKey = BaseParts[0] WCSTableName = DirName+BaseKey+'RADEC.tab' GPhyTableName = DirName+BaseKey+'_G_XY.tab' IPhyTableName = DirName+BaseKey+'_I_XY.tab' YPhyTableName = DirName+BaseKey+'_Y_XY.tab' #Extract various keys from the full path for each filter if UseG : GBaseName = os.path.basename(GFitsFileName) DirName = os.path.dirname(GFitsFileName) OutDir = DirName+'/Photometry/' GBaseParts = GBaseName.split('_') GObjectName = GBaseParts[0] GThisFilter = GBaseParts[1] GBaseParts = GBaseName.split('.') GBaseKey = GBaseParts[0] GDaoFindFile = OutDir+GBaseKey+'.coo' GDaoPhotFile = OutDir+GBaseKey+'.mag' GFStList = OutDir+GBaseKey+'.pst' GPSFGrList = OutDir+GBaseKey+'.psg' GPSFImHFile = OutDir+GBaseKey+'.psf.fits' GNStarFile = OutDir+GBaseKey+'.nst' GNRejFile = OutDir+GBaseKey+'.nrj' GAllStarFile = OutDir+GBaseKey+'.als' GRejFile = OutDir+GBaseKey+'.arj' GAllSubFile = OutDir+GBaseKey+'All.sub.fits' GIm = pyfits.getdata(GFitsFileName) GSubIm = pyfits.getdata(GAllSubFile) if UseI: IBaseName = os.path.basename(IFitsFileName) DirName = os.path.dirname(IFitsFileName) OutDir = DirName+'/Photometry/' IBaseParts = IBaseName.split('_') IObjectName = IBaseParts[0] IThisFilter = IBaseParts[1] IBaseParts = IBaseName.split('.') IBaseKey = IBaseParts[0] IIm = pyfits.getdata(IFitsFileName) IDaoFindFile = OutDir+IBaseKey+'.coo' IDaoPhotFile = OutDir+IBaseKey+'.mag' IFStList = OutDir+IBaseKey+'.pst' IPSFGrList = OutDir+IBaseKey+'.psg' IPSFImHFile = OutDir+IBaseKey+'.psf.fits' INStarFile = OutDir+IBaseKey+'.nst' INRejFile = OutDir+IBaseKey+'.nrj' IAllStarFile = OutDir+IBaseKey+'.als' IRejFile = OutDir+IBaseKey+'.arj' IAllSubFile = OutDir+IBaseKey+'All.sub.fits' IIm = pyfits.getdata(IFitsFileName) ISubIm = pyfits.getdata(IAllSubFile) if UseY: YBaseName = os.path.basename(YFitsFileName) DirName = os.path.dirname(YFitsFileName) YBaseParts = YBaseName.split('_') YObjectName = YBaseParts[0] YThisFilter = YBaseParts[1] YBaseParts = YBaseName.split('.') YBaseKey = YBaseParts[0] YDaoFindFile = OutDir+YBaseKey+'.coo' YDaoPhotFile = OutDir+YBaseKey+'.mag' YFStList = OutDir+YBaseKey+'.pst' YPSFGrList = OutDir+YBaseKey+'.psg' YPSFImHFile = OutDir+YBaseKey+'.psf.fits' YNStarFile = OutDir+YBaseKey+'.nst' YNRejFile = OutDir+YBaseKey+'.nrj' YAllStarFile = OutDir+YBaseKey+'.als' YRejFile = OutDir+YBaseKey+'.arj' YAllSubFile = OutDir+YBaseKey+'All.sub.fits' YIm = pyfits.getdata(YFitsFileName) YSubIm = pyfits.getdata(YAllSubFile) #Create a directory for Photometry output OutDir = DirName+'/Photometry/' if not os.path.exists(OutDir): os.mkdir(OutDir) #Get some of the keywords from the keywords file if os.path.isfile('/afs/slac/u/ki/lances/python/keywords.py'): execfile('/afs/slac/u/ki/lances/python/keywords.py'); #GET THE DATE Fullpath=os.getcwd() #Import the good packages from IRAF iraf.digiphot(_doprint=0) iraf.daophot(_doprint=0) iraf.apphot(_doprint=0) #KEY PARAMETERS Paramters to use FWHM = 10.5 #Full Width at Half Max GAIN = 5.0 #e-/ADU NUMFILES = 9.0 #Number of files used to make dither GAIN = NUMFILES*GAIN #If N images are averaged, GAIN=N*GAIN print 'Using FWHM of ' + str(FWHM) #PHOTPARS ############################################################# PhotParsFileName = OutDir+GObjectName+'_'+GThisFilter+'PhotPar.par' iraf.photpars.setParList(ParList = PhotParsFileName) #CENTERPARS ############################################################# CenterParsFileName = OutDir+GObjectName+'_'+GThisFilter+'CenterPar.par' iraf.centerpars.setParList(ParList=CenterParsFileName) #SKYPARS ################################################################ FitSkyParsFileName = OutDir+GObjectName+'_'+GThisFilter+'FitSkyPar.par' iraf.fitskypars.setParList(ParList=FitSkyParsFileName) #FINDPARS ############################################################## FindParsFileName = OutDir+GObjectName+'_'+GThisFilter+'FindPar.par' iraf.findpars.setParList(ParList = FindParsFileName) #DATAPARS - http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?datapars.hlp DataParsFileName = OutDir+GObjectName+'_'+GThisFilter+'DataPar.par' iraf.datapars.setParList(ParList=DataParsFileName) #DAOPARS ############################################################### DAOParsFileName = OutDir+GObjectName+'_'+GThisFilter+'DAOPar.par' iraf.daopars.setParList(ParList = DAOParsFileName) #Get the magnitudes, magnitude errors, CHI and SHARP from the .als file #G Image Data ######################## if UseG: GAllStarComps = iraf.pdump(GAllStarFile,\ 'xcenter, ycenter, mag, merr, chi, sharpness, msky', \ 'mag != INDEF', Stdout = 1) GAllStarLen = len(GAllStarComps) print 'There are : '+str(GAllStarLen) + ' stars in G' GXCen = zeros(GAllStarLen, dtype=float64) GYCen = zeros(GAllStarLen, dtype=float64) GMags = zeros(GAllStarLen, dtype=float64) GMErrs = zeros(GAllStarLen, dtype=float64) GChi = zeros(GAllStarLen, dtype=float64) GSharp = zeros(GAllStarLen, dtype=float64) GMSky = zeros(GAllStarLen, dtype=float64) for i, line in zip(range(GAllStarLen),GAllStarComps) : Pars = line.split() GXCen[i] = float(Pars[0]) GYCen[i] = float(Pars[1]) GMags[i] = float(Pars[2]) GMErrs[i] = float(Pars[3]) GChi[i] = float(Pars[4]) GSharp[i] = float(Pars[5]) GMSky[i] = float(Pars[6]) #I Image Data ##################################################### if UseI: IAllStarComps = iraf.pdump(IAllStarFile, \ 'xcenter, ycenter, mag, merr, chi, sharpness, msky', \ 'mag != INDEF', Stdout = 1) IAllStarLen = len(IAllStarComps) print 'There are: '+str(IAllStarLen)+ ' stars in I' IXCen = zeros(IAllStarLen, dtype=float64) IYCen = zeros(IAllStarLen, dtype=float64) IMags = zeros(IAllStarLen, dtype=float64) IMErrs = zeros(IAllStarLen, dtype=float64) IChi = zeros(IAllStarLen, dtype=float64) ISharp = zeros(IAllStarLen, dtype=float64) IMSky = zeros(IAllStarLen, dtype=float64) #Attempt to match the magnitudes GIGDupMags = [0] GIIDupMags = [0] GIDupMags = [0] GIGCutDupMags = [0] GIICutDupMags = [0] GICutDupMags = [0] #Get the offsets between the files GIMatchNum = 0 for i, line in zip(range(IAllStarLen),IAllStarComps) : Pars = line.split() IXCen[i] = float(Pars[0]) IYCen[i] = float(Pars[1]) IMags[i] = float(Pars[2]) IMErrs[i] = float(Pars[3]) IChi[i] = float(Pars[4]) ISharp[i] = float(Pars[5]) IMSky[i] = float(Pars[6]) if UseG: #Try to match the array XCenMatch = (where(abs((IXCen[i]-GIXOffset)-GXCen) < PixTh))[0] YCenMatch = (where(abs((IYCen[i]-GIYOffset)-GYCen[XCenMatch]) < PixTh))[0] if len(YCenMatch) == 1 : MatchInd = XCenMatch[YCenMatch] print 'Y Index'+str(MatchInd) print 'Found Match : '+str(GIMatchNum) + ' at star number ' +str(i) print 'XCen G: ' + str(GXCen[MatchInd]) + 'XCen I: '+str(IXCen[i]) print 'YCen G: ' + str(GYCen[MatchInd]) + 'YCen I: '+str(IYCen[i]) GIMatchNum +=1 GIGDupMags=append(GIGDupMags, GMags[MatchInd]-GZero) GIIDupMags=append(GIIDupMags, IMags[i]-IZero) GIDupMags =append(GIDupMags , GMags[MatchInd]-GZero-(IMags[i]-IZero)) if IChi[i] < 2 and GChi[MatchInd] < 2: GIGCutDupMags=append(GIGCutDupMags, GMags[MatchInd]-GZero) GIICutDupMags=append(GIICutDupMags, IMags[i]-IZero) GICutDupMags =append(GICutDupMags , GMags[MatchInd]-GZero-(IMags[i]-IZero)) #Y Image Data ##################################################### if UseY: YAllStarComps = iraf.pdump(YAllStarFile, \ 'xcenter, ycenter, mag, merr, chi, sharpness, msky', \ 'mag != INDEF', Stdout = 1) YAllStarLen = len(YAllStarComps) print 'There are: '+str(IAllStarLen)+ ' stars in I' YXCen = zeros(YAllStarLen, dtype=float64) YYCen = zeros(YAllStarLen, dtype=float64) YMags = zeros(YAllStarLen, dtype=float64) YMErrs = zeros(YAllStarLen, dtype=float64) YChi = zeros(YAllStarLen, dtype=float64) YSharp = zeros(YAllStarLen, dtype=float64) YMSky = zeros(YAllStarLen, dtype=float64) #Attempt to match the magnitudes GYGDupMags = [0] GYYDupMags = [0] GYDupMags = [0] IYIDupMags = [0] IYYDupMags = [0] IYDupMags = [0] GYGCutDupMags = [0] GYYCutDupMags = [0] GYCutDupMags = [0] for i, line in zip(range(YAllStarLen), YAllStarComps) : Pars = line.split() YXCen[i] = float(Pars[0]) YYCen[i] = float(Pars[1]) YMags[i] = float(Pars[2]) YMErrs[i] = float(Pars[3]) YChi[i] = float(Pars[4]) YSharp[i] = float(Pars[5]) YMSky[i] = float(Pars[6]) #Try to match the array against G and I if UseG : XCenMatch = (where(abs((YXCen[i]-GYXOffset)-GXCen) < PixTh))[0] YCenMatch = (where(abs((YYCen[i]-GYYOffset)-GYCen[XCenMatch]) < PixTh))[0] if len(YCenMatch) == 1 : MatchInd = XCenMatch[YCenMatch] print 'Y Index'+str(MatchInd) print 'Found Match : '+str(GYMatchNum) + ' at star number ' +str(i) print 'XCen G: ' + str(GXCen[MatchInd]) + 'XCen Y: '+str(YXCen[i]) print 'YCen G: ' + str(GYCen[MatchInd]) + 'YCen Y: '+str(YYCen[i]) GYMatchNum +=1 GYGDupMags=append(GYGDupMags,GMags[MatchInd]-GZero) GYYDupMags=append(GYYDupMags,YMags[i]-YZero) GYDupMags =append(GYDupMags, GMags[MatchInd]-GZero-(YMags[i]-YZero)) if YChi[i] < 2 : GYGCutDupMags=append(GYGCutDupMags,GMags[MatchInd]-GZero) GYYCutDupMags=append(GYYCutDupMags,YMags[i]-YZero) GYCutDupMags =append(GYCutDupMags, GMags[MatchInd]-GZero-(YMags[i]-YZero)) if UseI : XCenMatch = (where(abs((YXCen[i]-IYXOffset)-IXCen) < PixTh))[0] YCenMatch = (where(abs((YYCen[i]-IYYOffset)-IYCen[XCenMatch]) < PixTh))[0] if len(YCenMatch) == 1 : MatchInd = XCenMatch[YCenMatch] print 'Y Index'+str(MatchInd) print 'Found Match : '+str(IYMatchNum) + ' at star number ' +str(i) print 'XCen I: ' + str(IXCen[MatchInd]) + 'XCen Y: '+str(YXCen[i]) print 'YCen I: ' + str(IYCen[MatchInd]) + 'YCen Y: '+str(YYCen[i]) IYMatchNum +=1 IYIDupMags=append(IYIDupMags,IMags[MatchInd]-IZero) IYYDupMags=append(IYYDupMags,YMags[i]-YZero) IYDupMags =append(IYDupMags, IMags[MatchInd]-IZero-(YMags[i]-YZero)) print 'In G and I found : ' + str(GIMatchNum) print 'In G and Y found : ' + str(GYMatchNum) print 'In I and Y found : ' + str(IYMatchNum) if TvFlag > 0 : #Plot the magnitudes vs. error for G mplot.clf() mplot.semilogy(GMags,GMErrs,'b+') mplot.xlabel('Magnitude') mplot.ylabel('Magnitude Error') mplot.title('Magnitudes for '+GThisFilter+' Band') mplot.savefig(OutDir+'GMagVsMagErr.png') #Plot the magnitudes vs. error for I mplot.cla() mplot.semilogy(IMags,IMErrs,'r+') mplot.xlabel('Magnitude') mplot.ylabel('Magnitude Error') mplot.title('Magnitudes for '+IThisFilter+' Band') mplot.savefig(OutDir+'IMagVsMagErr.png') #Plot the magnitudes vs. error for I mplot.cla() mplot.plot(IChi,IMags,'r+') mplot.xlabel('Chi (Goodness of fit)') mplot.ylabel('Magnitude') mplot.title('Magnitudes for '+IThisFilter+' Band') mplot.savefig(OutDir+'IMagVsChi.png') mplot.cla() mplot.plot(GChi,GMags,'b+') mplot.xlabel('Chi (Goodness of fit)') mplot.ylabel('Magnitude') mplot.title('Magnitudes for '+GThisFilter+' Band') mplot.savefig(OutDir+'GMagVsChi.png') #Plot the median sky vs. magnitude mplot.cla() mplot.plot(GMSky ,GMags,'b+') mplot.xlabel('Median Sky') mplot.ylabel('Magnitude') mplot.title('Magnitudes for '+GThisFilter+' Band') mplot.savefig(OutDir+'GMagVsSky.png') mplot.cla() mplot.plot(IMSky ,IMags,'r+') mplot.xlabel('Median Sky') mplot.ylabel('Magnitude') mplot.title('Magnitudes for '+GThisFilter+' Band') mplot.savefig(OutDir+'IMagVsSky.png') #Plot the Color Magnitude Diagram for G and Y - 0th element is zero GYDupMags = GYDupMags[1:] GYYDupMags = GYYDupMags[1:] GYGDupMags = GYGDupMags[1:] mplot.clf() mplot.axis([GYDupMags.min(), GYDupMags.max(),\ GYYDupMags.max(), GYYDupMags.min()]) mplot.plot(GYDupMags, GYYDupMags, 'b+') mplot.title('Color Magnitude Diagram for M13') mplot.xlabel('g-y') mplot.ylabel('Magnitude in y') mplot.figtext(.15,.85, str(GAllStarLen) + ' stars in G') mplot.figtext(.15,.81, str(YAllStarLen) + ' stars in Y') mplot.figtext(.15,.77, str(GYMatchNum) + ' matched stars') mplot.savefig(OutDir+'M13ColorGYMagDiag.eps') #Plot the Color Magnitude Diagram for G and Y - 0th element is zero GYCutDupMags = GYCutDupMags[1:] GYYCutDupMags = GYYCutDupMags[1:] GYGCutDupMags = GYGCutDupMags[1:] mplot.clf() mplot.axis([GYCutDupMags.min(), GYCutDupMags.max(),\ GYYCutDupMags.max(), GYYCutDupMags.min()]) mplot.plot(GYCutDupMags, GYYCutDupMags, 'b+') mplot.title('Color Magnitude Diagram for M13') mplot.xlabel('g-y') mplot.ylabel('Magnitude in y') mplot.figtext(.15,.85, str(GAllStarLen) + ' stars in G') mplot.figtext(.15,.81, str(YAllStarLen) + ' stars in Y') mplot.figtext(.15,.77, str(GYMatchNum) + ' matched stars') mplot.savefig(OutDir+'M13ColorGYMagDiagCutChi2.eps') #Plot the Color Magnitude Diagram for G and Y - 0th element is zero IYDupMags = IYDupMags[1:] IYYDupMags = IYYDupMags[1:] IYIDupMags = IYIDupMags[1:] mplot.clf() mplot.axis([IYDupMags.min(), IYDupMags.max(),\ IYYDupMags.max(), IYYDupMags.min()]) mplot.plot(IYDupMags, IYYDupMags, 'b+') mplot.title('Color Magnitude Diagram for M13') mplot.xlabel('I-Y') mplot.ylabel('Magnitude in Y') mplot.figtext(.15,.85, str(IAllStarLen) + ' stars in I') mplot.figtext(.15,.81, str(YAllStarLen) + ' stars in Y') mplot.figtext(.15,.77, str(IYMatchNum) + ' matched stars') mplot.savefig(OutDir+'M13ColorIYMagDiag.eps') #Plot the Color Magnitude Diagram for G and I - 0th element is zero GIDupMags = GIDupMags[1:] GIIDupMags = GIIDupMags[1:] GIGDupMags = GIGDupMags[1:] mplot.clf() mplot.axis([GIDupMags.min(), GIDupMags.max(),\ GIIDupMags.max(), GIIDupMags.min()]) mplot.plot(GIDupMags, GIIDupMags, 'b+') mplot.title('Color Magnitude Diagram for M13') mplot.xlabel('g-i') mplot.ylabel('Magnitude in i') mplot.figtext(.15,.85, str(GAllStarLen) + ' stars in G') mplot.figtext(.15,.81, str(IAllStarLen) + ' stars in I') mplot.figtext(.15,.77, str(GIMatchNum) + ' matched stars') mplot.savefig(OutDir+'M13ColorGIMagDiag.eps') #Plot the Color Magnitude Diagram for G and I - 0th element is zero GICutDupMags = GICutDupMags[1:] GIICutDupMags = GIICutDupMags[1:] GIGCutDupMags = GIGCutDupMags[1:] mplot.clf() mplot.axis([GICutDupMags.min(), GICutDupMags.max(),\ GIICutDupMags.max(), GIICutDupMags.min()]) mplot.plot(GICutDupMags, GIICutDupMags, 'b+') mplot.title('Color Magnitude Diagram for M13') mplot.xlabel('g-i') mplot.ylabel('Magnitude in i') mplot.figtext(.15,.85, str(GAllStarLen) + ' stars in G') mplot.figtext(.15,.81, str(IAllStarLen) + ' stars in I') mplot.figtext(.15,.77, str(GIMatchNum) + ' matched stars') mplot.savefig(OutDir+'M13ColorGIMagDiagCutChi2.eps') #Tranformed to U #Print some values print 'There are : '+str(GAllStarLen) + ' stars in G' print 'There are : '+str(IAllStarLen) + ' stars in I' print 'There are : '+str(YAllStarLen) + ' stars in Y' print 'In G and I found : ' + str(GIMatchNum) print 'In G and Y found : ' + str(GYMatchNum) print 'In I and Y found : ' + str(IYMatchNum) #Now, if catalog was specified as a keyword, try to make matches with it if Catalog != '' : #The format of the catalog is described in the header CatalogEnt = TableIO.readTableAsArray(Catalog,'#') CIDs = CatalogEnt[:,0] #ID of Star CVMags = CatalogEnt[:,4] #Magnitude in V CBVMags = CatalogEnt[:,5] #B-V CBMags = CVMags+CBVMags #Magnitude in B CUBMags = CatalogEnt[:,6] #U-B CUMags = CBMags + CUBMags #Magnitude in U CPropX = CatalogEnt[:,7] #Proper motion in X-mas CPropY = CatalogEnt[:,8] #Proper motion in Y-mas CRAs = CatalogEnt[:,10] #Right Ascension CDECs = CatalogEnt[:,11] #Declination CRADECs = numpy.array([CRAs,CDECs]).transpose() CRADEClist = [] for i in range(len(CRAs)): #CRADEClist= append(CRADEClist, [str(CRADECs[i,0])+' '+str(CRADECs[i,1])]) CRADEClist.extend([str(CRADECs[i,0])+' '+str(CRADECs[i,1])]) #Get the physical coordinates in each image iraf.wcsctran.setParam('input', 'STDIN') iraf.wcsctran.setParam('output',"STDOUT") iraf.wcsctran.setParam('inwcs','world') iraf.wcsctran.setParam('outwcs','tv') #G, I, and Y slightly off iraf.wcsctran.setParam('image',GFitsFileName) GXY = iraf.wcsctran(Stdin=CRADEClist,Stdout=1,mode='h') iraf.wcsctran.setParam('image',IFitsFileName) IXY = iraf.wcsctran(Stdin=CRADEClist,Stdout=1,mode='h') iraf.wcsctran.setParam('image',YFitsFileName) YXY = iraf.wcsctran(Stdin=CRADEClist,Stdout=1,mode='h') XYLen = len(GXY) #Attempt to match the magnitudes CIDupMags = [0.] CGDupMags = [0.] CYDupMags = [0.] CatLen = len(CIDs) CGVMags = [0.] #Keep Track of the diff g-V CGBVMags = [0.] #And the corresponding B-V #Keep track of the pixel offsets between catalog and our observed centroids CGXDiffs = [0] ; CGYDiffs = [0] ; CGMatchInd = [0] CIXDiffs = [0] ; CIYDiffs = [0] ; CIMatchInd = [0] CYXDiffs = [0] ; CYYDiffs = [0] ; CYMatchInd = [0] #Go through the catalog entries and match them to the G,I,Y stars CGMatchNum = 0 ; CIMatchNum = 0 ; CYMatchNum = 0 ; CMatchNum=0 #Open a Obsfile for later input to the transformation algorithm ObsFileName = OutDir+GObjectName+'GIY'+'ObsFile.obs' ObsFile = open(ObsFileName,'w') StandFileName = OutDir+GObjectName+'Standard.std' StdFile = open(StandFileName,'w') #The catalog might have comments at the top #Valid lines should be indexed separately from the loop index ValidInd = 0 for i in range(CatLen): FoundOne = 0 ; FoundG = 0 ; FoundI = 0 ; FoundY = 0 XCenMatch= 0 ; YCenMatch = 0 if re.search('\d\.\d\d\d\d', GXY[i]): GRow = GXY[i].split() CatGX = float(GRow[0]) ; CatGY = float(GRow[1]) IRow = IXY[i].split() CatIX = float(IRow[0]) ; CatIY = float(IRow[1]) YRow = YXY[i].split() CatYX = float(YRow[0]) ; CatYY = float(YRow[1]) #Construct the output lines for the ObsFile if UseG : #Try to match the array XCenMatch = (where(abs(CatGX-GXCen) < PixTh))[0] YCenMatch = (where(abs(CatGY-GYCen[XCenMatch]) < PixTh))[0] if len(YCenMatch) == 1 : FoundOne = 1 ; FoundG = 1 MatchInd = XCenMatch[YCenMatch] print 'G Index'+str(MatchInd) print 'Found Match : '+str(CGMatchNum) + ' at star number ' +str(i) CGXDiffs = append(CGXDiffs, CatGX-GXCen[MatchInd][0]) CGYDiffs = append(CGYDiffs, CatGY-GYCen[MatchInd][0]) CGVMags = append(CGVMags, CVMags[ValidInd]-GMags[MatchInd]) print str(CVMags[ValidInd]-GMags[MatchInd]) CGBVMags = append(CGBVMags, CBVMags[ValidInd]) CGMatchInd = append(CGMatchInd, i) CGDupMags = append(CGDupMags, GMags[MatchInd]) GXStr = str(GXCen[MatchInd][0]); GYStr = str(GYCen[MatchInd][0]) GMagStr = str(GMags[MatchInd][0]); GMErrsStr= str(GMErrs[MatchInd][0]) CGMatchNum+=1 else : GXStr = 'INDEF' ; GYStr = 'INDEF' GMagStr = 'INDEF' ; GMErrsStr = 'INDEF' if UseI : XCenMatch = (where(abs(CatIX-IXCen) < PixTh))[0] YCenMatch = (where(abs(CatIY-IYCen[XCenMatch]) < PixTh))[0] if len(YCenMatch) == 1 : FoundOne = 1 FoundI = 1 MatchInd = XCenMatch[YCenMatch] print 'I Index'+str(MatchInd) print 'Found Match : '+str(CIMatchNum) + ' at star number ' +str(i) CIXDiffs = append(CIXDiffs, CatIX-IXCen[MatchInd][0]) CIYDiffs = append(CIYDiffs, CatIY-IYCen[MatchInd][0]) CIMatchInd = append(CIMatchInd, i) CIDupMags = append(CIDupMags, IMags[MatchInd]) IXStr = str(IXCen[MatchInd][0]); IYStr = str(IYCen[MatchInd][0]) IMagStr = str(IMags[MatchInd][0]); IMErrStr = str(IMErrs[MatchInd][0]) CIMatchNum +=1 else : IXStr = 'INDEF' ; IYStr = 'INDEF' IMagStr = 'INDEF' ; IMErrsStr = 'INDEF' if UseY : XCenMatch = (where(abs(CatYX-YXCen) < PixTh))[0] YCenMatch = (where(abs(CatYY-YYCen[XCenMatch]) < PixTh))[0] if len(YCenMatch) == 1 : Found1 = 1 FoundY = 1 MatchInd = XCenMatch[YCenMatch] print 'Y Index'+str(MatchInd) print 'Found Match : '+str(CYMatchNum) + ' at star number ' +str(i) CYXDiffs = append(CYXDiffs, CatYX-YXCen[MatchInd][0]) CYYDiffs = append(CYYDiffs, CatYY-YYCen[MatchInd][0]) CYMatchInd = append(CYMatchInd, i) CYDupMags = append(CYDupMags, YMags[MatchInd]) YXStr = str(YXCen[MatchInd][0]); YYStr = str(YYCen[MatchInd][0]) YMagStr = str(YMags[MatchInd][0]); YMErrsStr= str(YMErrs[MatchInd][0]) CYMatchNum +=1 else : YXStr = 'INDEF' ; YYStr = 'INDEF' YMagStr = 'INDEF' ; YMErrsStr = 'INDEF' #If at least one match was made, write this to the Observation file if FoundOne == 1 : #Print out ID, FILTER, AIRMASS, XCENTER, YCENTER, MAG, MERR ObsFile.write(str(CMatchNum).ljust(10)) ObsFile.write('00'.ljust(10)) ObsFile.write(str(GAirMass).ljust(10)) ObsFile.write(GXStr.ljust(10)) ObsFile.write(GYStr.ljust(10)) ObsFile.write(GMagStr.ljust(10)) ObsFile.write(GMErrsStr.ljust(10)) ObsFile.write('\n') ObsFile.write('*'.ljust(10)) ObsFile.write('10'.ljust(10)) ObsFile.write(str(IAirMass).ljust(10)) ObsFile.write(IXStr.ljust(10)) ObsFile.write(IYStr.ljust(10)) ObsFile.write(IMagStr.ljust(10)) ObsFile.write(IMErrsStr.ljust(10)) ObsFile.write('\n') ObsFile.write('*'.ljust(10)) ObsFile.write('20'.ljust(10)) ObsFile.write(str(YAirMass).ljust(10)) ObsFile.write(YXStr.ljust(10)) ObsFile.write(YYStr.ljust(10)) ObsFile.write(YMagStr.ljust(10)) ObsFile.write(YMErrsStr.ljust(10)) ObsFile.write('\n') #Write this to the Standard Observation File StdFile.write(str(CMatchNum).ljust(10)) StdFile.write(str(CVMags[ValidInd]).ljust(10)) StdFile.write(str(CBVMags[ValidInd]).ljust(10)) StdFile.write(str(CUBMags[ValidInd]).ljust(10)) StdFile.write(str(CatGX).ljust(15)) StdFile.write(str(CatGY).ljust(15)) #StdFile.write(str(RAs[i]).ljust(10)) StdFile.write('\n') CMatchNum+=1 #Increment this if a valid line was found ValidInd+=1 #Close the file ObsFile.close() StdFile.close() print 'Found ' + str(CGMatchNum) + 'stars from G' print 'Found ' + str(CIMatchNum) + 'stars from I' print 'Found ' + str(CYMatchNum) + 'stars from Y' if TvFlag > 0 : #Calculate the offset w/r to proper motions CGRDiff = sqrt(CGXDiffs**2+CGYDiffs**2) CGRProp = sqrt(CPropX[CGMatchInd]**2+CPropY[CGMatchInd]**2) CIRDiff = sqrt((CIXDiffs-1)**2+CIYDiffs**2) CIRProp = sqrt(CPropX[CIMatchInd]**2+CPropY[CIMatchInd]**2) CYRDiff = sqrt(CYXDiffs**2+CYYDiffs**2) CYRProp = sqrt(CPropX[CYMatchInd]**2+CPropY[CYMatchInd]**2) #Plot the proper motions and the centroid offset scatters: mplot.clf() mplot.plot(CGRDiff,CGRProp,'b+') mplot.axis([CGRDiff.min(), CGRDiff.max(),\ 0, 20]) mplot.plot(CIRDiff,CIRProp,'g+') mplot.axis([CIRDiff.min(), CIRDiff.max(),\ 0, 20]) mplot.plot(CYRDiff,CYRProp,'r+') mplot.axis([CYRDiff.min(), CYRDiff.max(),\ 0, 20]) mplot.clf() MeanCGXDiff = mean(CGXDiffs) MeanCGYDiff = mean(CGYDiffs) StdCGXDiff = std(CGXDiffs) StdCGYDiff = std(CGYDiffs) print 'Std. Dev of X Centroid Diffs in G : ' + str(PltSc*StdCGXDiff) print 'Std. Dev of Y Centroid Diffs in G : ' + str(PltSc*StdCGYDiff) SkyGXDiff = PltSc*(CGXDiffs-MeanCGXDiff) SkyGYDiff = PltSc*(CGYDiffs-MeanCGYDiff) mplot.plot(SkyGXDiff,SkyGYDiff,'b+') MeanCIXDiff = mean(CIXDiffs) MeanCIYDiff = mean(CIYDiffs) StdCIXDiff = std(CIXDiffs) StdCIYDiff = std(CIYDiffs) print 'Std. Dev of X Centroid Diffs in I : ' + str(PltSc*StdCIXDiff) print 'Std. Dev of Y Centroid Diffs in I : ' + str(PltSc*StdCIYDiff) SkyIXDiff = PltSc*(CIXDiffs-MeanCIXDiff) SkyIYDiff = PltSc*(CIYDiffs-MeanCIYDiff) mplot.plot(SkyIXDiff,SkyIYDiff,'g+') MeanCYXDiff = mean(CYXDiffs) MeanCYYDiff = mean(CYYDiffs) StdCYXDiff = std(CYXDiffs) StdCYYDiff = std(CYYDiffs) print 'Std. Dev of X Centroid Diffs in Y : ' + str(PltSc*StdCYXDiff) print 'Std. Dev of Y Centroid Diffs in Y : ' + str(PltSc*StdCYYDiff) SkyYXDiff = PltSc*(CYXDiffs-MeanCYXDiff) SkyYYDiff = PltSc*(CYYDiffs-MeanCYYDiff) mplot.plot(SkyYXDiff,SkyYYDiff,'r+') mplot.axis([-.7,.7,-.7,.7]) mplot.title('Difference Between Measured Centroids and Catalog Entries') mplot.xlabel('Difference in X (Arcseconds)') mplot.ylabel('Difference in Y (Arcseconds)') mplot.legend(('G','I','Y')) mplot.savefig(OutDir+'M13CentroidDifferences.png') mplot.clf() mplot.plot(CGBVMags,CGVMags,'r+')