#----------------------------------------------------------- #Utility script to read exoplanet and parent star #data from the California Carnegie Planet Search tables #(http://exoplanets.org), and compute various derived #quantities and simplified tables for plotting. This #database is quite well quality-controlled, but there is #some bad data in there. For example, one of OGLE stars #in the April 2010 database is listed with a much too large #distance, which makes its luminosity come out absurdly high #(about a million times that of the Sun). This may be corrected #in future updates, but be wary of strange numbers. # #Note that information on extrasolar planets can also be found #at exoplanet.eu. The site exoplanets.org, maintained by Geoff #Marcy, produces output that is generally easier to work with, though. # # #The input tables contain non-numeric data, #so the readTable utility can't be used. The data #crunching to extract the data from the format in the #database, and needed to link planet data with star data #involves some Python constructs that the student is not #expected to have mastered at this point. The constructs #in question involve file input/output and dictionaries. #This script is nonetheless included as a chapter script #so that the student can do problems involving the most #current list of planets. # #To use this script: # (a) Download the a file of star data and a file # of planet data from exoplanets.org. The 2010 # files of star data and planet data are already located # in the Chapter 1 datasets directory, so you need # not download them. A processed output file # made from this data is also in that directory, # and can be plotted using ReadExoplanets.py. You # only need to run this script if you download updated # data, or if you want to modify it to save some information # into the output file that isn't saved by default. # # (b) Edit the file names for the files in the # script, so that the script can find them # # (c) Run the script. It will dump a selection # of useful planet data in column format to # the file PlanetOut.temp.txt . This can be # read into another script using readTable(...) # from ClimateUtilities. # # You can also modify the output by modifying # what is stuffed into the Curve() object for output. # Also, exoplanets.org gives you a lot of control over # what data you want, so if you want to play with # other kinds of data (like the "metalicity" of the star) # you can select to get those columns, and then edit # the loops so that they save the data you want to a list. # To map the column name to the column number, just follow # the example given for 'Vmag', but use the column name # (or add a "standardized" synonym to the dictionary of # standard names). #------------------------------------------------------------------ #ToDo: # # **Add an output option which compiles this data # into an exoplanets.py data handbook, # analogous to planets.py # # **Compile a dictionary that allows the user to associate # a line in the processed data table with the name of # the planet (e.g. Gliese581c). import string,math from ClimateUtilities import * #Path to Workbook datasets datapath = '/Users/rtp1/Havsornen/PlanetaryClimateBook/WorkbookDatasets/' delimiter = ',' #Column separator for the data files Missing = 'NaN' #Missing data code #Edit the following to put in the file names you are using. #Remember to put in the path name to the directory if the #file is not in the same directory as this script. starFileName = 'stars2010.txt' #Star file name here. planetFileName = 'planets2010.txt' #Planet file name here. #Tack on the directory path starFileName = datapath+'Chapter1Data/exoplanets/' + starFileName planetFileName = datapath+'Chapter1Data/exoplanets/' + planetFileName #Functions needed for magnitude and T calculation def T(BV): logT = (14.551 - BV)/3.684 if logT > 3.961: a,b,c = [.344,-3.402,8.037] logT = (-b - (b*b-4*a*c)**.5)/(2.*a) return 10.**logT def Bolo(T): x = math.log10(T) - 4. return -8.499 * x**4 + 13.421* x**3- 8.131*x**2 - 3.901 *x - 0.438 #----------------Main Program Starts Here-------------------- #Open the files starFile = open(starFileName) planetFile = open(planetFileName) #Read in and process stellar data lines = starFile.readlines() #The following assumes the header is in line istart and #the units are in line istart+1. istart is hard-coded for now, #but we can add a scan routine later to locate the start of the #data. You may need to edit these if there is some extra stuff #at the top of the file. istart = 0 header = lines[istart].split(delimiter) header[-1] = header[-1][:-1] #Gets rid of the newline character at the end # #Now make up the dictionary mapping star related variables #to the corresponding columns, based on the header. This happens #automatically. You shouldn't need to change this column = 0 starDict = {} for h in header: starDict[h] = column column += 1 #Dictionary mapping standard names to column names. You will only #need to edit this if the column names have changed since the 2010 #version of the data base, or if there is additional data not listed #here, which you wish to look at. Even then, you don't need to edit #this dictionary if you are willing to refer to the data by its #column name. # #Dictionary of star table position of useful information # name Star name # Vmag V magnitude # BminusV B magnitude minus V magnitude (color index) # dist Distance in parsecs # Teff Effective temperature # FeToH Fe to H ratio (index of metallicity) # M Mass in units of Solar masses # RHK Magnetic activity index (based on Calcium H and K lines) #These are the ones used in the main calculation. Feel free to add more starDictStandard = {'name':'STAR','BminusV':'BMV','Vmag':'V',\ 'Teff':'TEFF','dist':'DISTANCE','M':'MSTAR'} #Now add synonyms to starDict for key in starDictStandard.keys(): starDict[key] = starDict[starDictStandard[key]] #Empty lists to hold the data we want to save ColorIndex = [] TList = [] IBoloList = [] starNames = [] MstarList = [] #List of star masses starNameIndex={} #Dictionaries linking star ID to index RHKList = [] #ToDo: Check for missing stellar data and put in #a missing data code if necessary. The following just #skips stars with missing data needed for luminosity calculation count = 0 for line in lines[istart+2:]: buff = line.split(',') #Check for missing data if 'NaN' in buff[starDict['Vmag']] or 'NaN' in buff[starDict['dist']]: continue if 'NaN' in buff[starDict['Teff']] and 'NaN' in buff[starDict['BminusV']]: continue #First get a star identifier so we can #link the planet with its star starName = buff[starDict['name']] #Put in the dictionary entry if len(starName.strip())>0: starNameIndex[starName.strip()] = count starNames.append(starName) else: print "Star name missing from row",count # #Now look for temperature data. If it's missing, #recompute temperature from B-V #Note: Is there a way to do this with a try ... except #block? Problem is that atof('nan') doesn't raise an exception BminusV = string.atof(buff[starDict['BminusV']]) #Fix the bad BmV in the database for GJ 176 if 'GJ 176' in starName: BminusV = 1.51 if 'NaN' in buff[starDict['Teff']]: Teff = T(BminusV) else: Teff = string.atof(buff[starDict['Teff']]) #There are a few glitches in the exoplanets.org #database, where Teff is given a low value instead #of NaN when the data is missing. Deal with that here if Teff < 10.: Teff = T(BminusV) TList.append(Teff) #Convert magnitude to absolute magnitude d = string.atof(buff[starDict['dist']]) Vmag = string.atof(buff[starDict['Vmag']]) AbsVmag = Vmag - 5.*math.log10(d/10.) #Add in bolometric correction AbsBoloMag = AbsVmag + Bolo(Teff) #Convert to V luminosity (relative to Sun) I = 10.**((4.72-AbsBoloMag)/2.5) IBoloList.append(I) # #Save other miscellaneous data MstarList.append(atof(buff[starDict['M']])) RHKList.append(atof(buff[starDict['RHK']])) # count += 1 #Put the star data in a separate Curve() object and #write it out separately from planets if you want. #--------Now read in the planet data----------------------- lines = planetFile.readlines() #The following assumes the header is in line istart and #the units are in line istart+1. istart is hard-coded for now, #but we can add a scan routine later to locate the start of the #data istart = 0 header = lines[istart].split(delimiter) header[-1] = header[-1][:-1] # #Now make up the dictionary mapping star related variables #to the corresponding columns, based on the header column = 0 planetDict = {} for h in header: planetDict[h] = column column += 1 #Dictionary mapping standard names to column names. You will only #need to edit this if the column names have changed since the 2010 #version of the data base, or if there is additional data not listed #here, which you wish to look at. Even then, you don't need to edit #this dictionary if you are willing to refer to the data by its #column name. # # # name Planet name # period orbital period in Earth days # ecc eccentricity # mass Planet maximum mass in Jupiter masses # rsm semi-major axis in A.U. # planetDictStandard = {'name':'NAME','period':'PER','ecc':'ECC',\ 'mass':'MSINI','rsm':'A'} #Now add synonyms to planetDict for key in planetDictStandard.keys(): planetDict[key] = planetDict[planetDictStandard[key]] #Empty lists to hold the data we want to save MPlanetList = [] rsmList = [] eccList = [] yearList = [] IPlanetList = [] TstarPlanetList = [] planetNames = [] count = 0 for line in lines[istart+2:]: buff = line.split(delimiter) #Extract the star name from the planet name planet = buff[planetDict['name']].strip() star = planet[:-1] star = star.strip() #Find index of the host star. In the #exoplanets.org database, as of 2010, every star #has a standard name. Note that some of the missing #stars are missing because luminosity data couldn't #be computed. Instead, we should put in a missing data #code for luminosity, and index the star name, so that #other characteristics of the planet can be plotted. #Since we are mostly interested in climate, we take the #view that if we don't know luminosity we can't do anything #much about the climate, so we leave out the planet try: starIndex = starNameIndex[star] except: print "Couldn't link planet ",planet, "to its star" continue #ToDo: Check for missing data ecc = string.atof(buff[planetDict['ecc']]) m = string.atof(buff[planetDict['mass']]) rsm = string.atof(buff[planetDict['rsm']]) per = string.atof(buff[planetDict['period']]) # planetNames.append(buff[planetDict['name']].strip()) MPlanetList.append(m) rsmList.append(rsm) eccList.append(ecc) yearList.append(per) IPlanetList.append(IBoloList[starIndex]) TstarPlanetList.append(TList[starIndex]) # count += 1 #Convert lists to numpy arrays so we can do arithmetic #on them IPlanet = numpy.array(IPlanetList) rsm = numpy.array(rsmList) ecc = numpy.array(eccList) # #Compute a few derived quantities axrat = (1.+ecc)/(1.-ecc) #Ratio of apastron to periastron distance LstarPeri = IPlanet/(rsm*rsm*(1-ecc)**2) #Stellar constant at periastron LstarAp = IPlanet/(rsm*rsm*(1+ecc)**2) #Stellar constant at apastron LstarSM = IPlanet/(rsm*rsm) #Stellar constant at semi-major axis. #For eccentricities up to .25 this is nearly #equal to the annual mean. When e = .5 , #the annual mean exceeds this by a factor #of 1.16. (see Chapter 7) #Put some data in a Curve object and write out c = Curve() c.addCurve(MPlanetList,'mPlanet') c.addCurve(rsmList,'rsm') c.addCurve(eccList,'ecc') c.addCurve(axrat,'PeriApRat') c.addCurve(yearList,'yearLength') c.addCurve(IPlanetList,'Luminosity') c.addCurve(LstarPeri,'PeriSolarConstant') c.addCurve(LstarAp,'ApSolarConstant') c.addCurve(LstarSM,'SMSolarConstant') c.addCurve(TstarPlanetList,'TStar') c.dump('exoplanetsProcessed.txt') #Puts the output in the current directory