-
Notifications
You must be signed in to change notification settings - Fork 10
/
plot_radec_discovery_mollweide.python
112 lines (97 loc) · 4.51 KB
/
plot_radec_discovery_mollweide.python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#!/usr/bin/python
import xml.etree.ElementTree as ET
import subprocess, glob, os, 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("""
# gnuplot input
reset
set terminal svg
set key below
set output "plot_radec_discovery_mollweide.svg"
set title "Positions of Extrasolar planet in Equatorial coordinate system"
set label "Data taken from the Open Exoplanet Catalogue. Last updated on """+datetime.date.today().isoformat()+"""." at screen 0.31,0.12 front font ",8"
# Mollweide projection in gnuplot (ref: http://astrojct.blogspot.jp/)
bconv(b) = b*pi/180
lconv(l) = (l<180) ? l*pi/180 : (l-360)*pi/180
mwhigh(x) = 1.0 + -0.919061*(abs(1.0-x))**0.674635
mwmed(x) = mwhigh(x) -0.0807765 + 0.161136*x -0.0796311*x**2
mwlow(x) = mwmed(x) -3.53551e-05 + 0.000645749*x
mwst(x) = x < 0.2 ? mwlow(x) : (x < 0.9 ? mwmed(x) : mwhigh(x) )
mwt(b) = b>0 ? asin(mwst(sin(abs(b)))) : -asin(mwst(sin(abs(b))))
mwx(b,l) = 2.0*sqrt(2.0) * l * cos(mwt(b))
mwy(b,l) = sqrt(2.0) * sin(mwt(b))
set parametric
unset border
unset xtics
unset ytics
set pointsize 0.75
plot [-pi:pi] \
"-" using (mwx(bconv($2),lconv($1))):(mwy(bconv($2),lconv($1))) w p pt 7 lc 2 t "Radial velocity planets", \
"-" using (mwx(bconv($2),lconv($1))):(mwy(bconv($2),lconv($1))) w p pt 7 lc 3 t "Transiting planets", \
"-" using (mwx(bconv($2),lconv($1))):(mwy(bconv($2),lconv($1))) w p pt 7 lc 1 t "Directly imaged planets", \
"-" using (mwx(bconv($2),lconv($1))):(mwy(bconv($2),lconv($1))) w p pt 7 lc 4 t "Microlensing planets", \
"-" using (mwx($3/180*pi,$2/180*pi)):(mwy($3*pi/180,$2*pi/180)):1 with labels notit, \
mwx(pi*(-2)/6.,t),mwy(pi*(-2)/6.,t) w l lt 0 notit, \
mwx(pi*(-1)/6.,t),mwy(pi*(-1)/6.,t) w l lt 0 notit, \
mwx(pi*( 0)/6.,t),mwy(pi*( 0)/6.,t) w l lt 0 notit, \
mwx(pi*(+1)/6.,t),mwy(pi*(+1)/6.,t) w l lt 0 notit, \
mwx(pi*(+2)/6.,t),mwy(pi*(+2)/6.,t) w l lt 0 notit, \
mwx(t/2.,pi*(-3)/4.),mwy(t/2.,pi*(-3)/4.) w l lt 0 notit, \
mwx(t/2.,pi*(-2)/4.),mwy(t/2.,pi*(-2)/4.) w l lt 0 notit, \
mwx(t/2.,pi*(-1)/4.),mwy(t/2.,pi*(-1)/4.) w l lt 0 notit, \
mwx(t/2.,pi*( 0)/4.),mwy(t/2.,pi*( 0)/4.) w l lt 0 notit, \
mwx(t/2.,pi*(+1)/4.),mwy(t/2.,pi*(+1)/4.) w l lt 0 notit, \
mwx(t/2.,pi*(+2)/4.),mwy(t/2.,pi*(+2)/4.) w l lt 0 notit, \
mwx(t/2.,pi*(+3)/4.),mwy(t/2.,pi*(+3)/4.) w l lt 0 notit, \
mwx(t/2.,pi-2.*pi*0),mwy(t/2.,pi-2.*pi*0) w l lt -1 lw 2 notit, \
mwx(t/2.,pi-2.*pi*1),mwy(t/2.,pi-2.*pi*1) w l lt -1 lw 2 notit
""")
# retrieve data
def plotPosition(_discoverymethod):
available = 0
num = 0
for filename in glob.glob("open_exoplanet_catalogue/systems/*.xml"):
system = ET.parse(open(filename, 'r'))
planets = system.findall(".//planet")
for planet in planets:
discoverymethod = planet.findtext("./discoverymethod")
if discoverymethod == _discoverymethod:
ra = system.findtext("./rightascension")
dec = system.findtext("./declination")
#if ra is not None:
try:
ra = ra.split()
rightascension = (float(ra[0]) + float(ra[1])/60. + float(ra[2])/3600.)*15
dec = dec.split()
front = float(dec[0])
if front > 0.: # is it the only way?
declination = front + float(dec[1])/60. + float(dec[2])/3600.
else:
declination = front - float(dec[1])/60. - float(dec[2])/3600.
gnuplot.stdin.write("%e\t%e\n"%(rightascension, declination))
available += 1
except: # because "None", ra or dec not specified for this system. -_-?
pass
num +=1
gnuplot.stdin.write("\ne\n")
none = num - available
#print 'System =', num
#print 'None =', none
plotPosition("RV")
plotPosition("transit")
plotPosition("imaging")
plotPosition("microlensing")
# labels
gnuplot.stdin.write("60\t0\t60\n")
gnuplot.stdin.write("30\t0\t30\n")
gnuplot.stdin.write("60\t0\t60\n")
gnuplot.stdin.write("0h,0\t0\t0\n")
gnuplot.stdin.write("-30\t0\t-30\n")
gnuplot.stdin.write("-60\t0\t-60\n")
gnuplot.stdin.write("9h\t135\t0\n")
gnuplot.stdin.write("6h\t90\t0\n")
gnuplot.stdin.write("21h\t-45\t0\n")
gnuplot.stdin.write("18h\t-90\t0\n")
gnuplot.stdin.write("15h\t-135\t0\n\ne\n")
gnuplot.stdin.close()