plot_periodratio.python

#!/usr/bin/python
import xml.etree.ElementTree as ET
import subprocess, glob, os, math, datetime

# Open pipe gnuplot. You may want to change the terminal from svg to pdf for publication quality plots.
gnuplot = subprocess.Popen(['gnuplot',"-persist"], stdin=subprocess.PIPE, stdout=subprocess.PIPE)
gnuplot.stdin.write("""
set terminal svg 
set output "plot_periodratio.svg" 
set ylabel "Number of pairs"
set title "Period ratios in multiplanetary systems"
set xrange [1.2:4]
set lmargin 8
set yrange [0:]
set style fill solid 1.0

set label "Data taken from the Open Exoplanet Catalogue. Last updated on """+datetime.date.today().isoformat()+"""." at screen 0.31,0.15 front font ",8"
binwidth = 0.1
bin(x,width)=width*floor(x/width) + binwidth/2.0
set boxwidth binwidth

set multiplot 
""")


minratio = 1.2
maxratio = 4.0

# RV planets
data = []
gnuplot.stdin.write("""
set size 1,0.5
set origin 0,0.5
unset xtics
unset xlabel
set yrange [0:8]
set bmargin 0
plot \
"-" using (bin($1,binwidth)):2 smooth freq with boxes t "Radial velocity systems"
""")
for filename in glob.glob("open_exoplanet_catalogue/systems/*.xml"):
	system = ET.parse(open(filename, 'r'))
	planets = system.findall(".//planet")
	if len(planets)>1:
		for planet1 in planets:
			try:
				discoverymethod = planet1.findtext("./discoverymethod")
				period1 = float(planet1.findtext("./period"))
				if discoverymethod=="RV":
					for planet2 in planets:
						period2 = float(planet2.findtext("./period"))
						ratio = period1/period2
						if ratio>minratio and ratio<maxratio:
							data.append([ratio])
			except:
				pass

def compare(x, y):
	if x[0]>y[0]: return 1
	if x[0]<y[0]: return -1
	return 0

sorteddata = sorted(data, cmp=compare)

for ratio in sorteddata:
	gnuplot.stdin.write("%e\t%e\n"%(ratio[0],1.))

gnuplot.stdin.write("\ne\n")

gnuplot.stdin.write("""
unset title
set size 1,0.5
set origin 0,0.
set tmargin 0
set xtics
set bmargin 
set yrange [0:95]
set xlabel "Period ratio"
set ytics 10
plot \
"-" using (bin($1,binwidth)):2 smooth freq with boxes ls 2 t "Kepler objects of interest (KOI)" 
""")

# Kepler planets
data = []
for filename in glob.glob("open_exoplanet_catalogue/systems_kepler/*.xml"):
	system = ET.parse(open(filename, 'r'))
	planets = system.findall(".//planet")
	if len(planets)>1:
		for planet1 in planets:
			try:
				period1 = float(planet1.findtext("./period"))
				for planet2 in planets:
					period2 = float(planet2.findtext("./period"))
					ratio = period1/period2
					if ratio>minratio and ratio<maxratio:
						data.append([ratio])
			except:
				pass

def compare(x, y):
	if x[0]>y[0]: return 1
	if x[0]<y[0]: return -1
	return 0

sorteddata = sorted(data, cmp=compare)

for ratio in sorteddata:
	gnuplot.stdin.write("%e\t%e\n"%(ratio[0],1.))

gnuplot.stdin.write("\ne\n")
gnuplot.stdin.close()