-
Notifications
You must be signed in to change notification settings - Fork 0
/
plot_individual_singles.py
129 lines (88 loc) · 3.02 KB
/
plot_individual_singles.py
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
113
114
115
116
117
118
119
120
121
## Plot the eccentricity preference for each single planet system
##
import numpy as np
import matplotlib.pyplot as plt
import random
import scipy
from scipy import stats
import os
from scipy import interpolate
from cksecc_helper import *
cks_data, koi_data, koi_errs = load_data()
koilist, propertymatrix = divide_into_singles_and_multis(cks_data, koi_data, koi_errs, singles=True, more_properties=True)
e, inc, like, emaxval = read_in_singles()
propertymatrix.append(like)
koilist, propertymatrix = clean_sample(koilist, propertymatrix, customremove=True)
msun = propertymatrix[0]
rsun = propertymatrix[1]
period = propertymatrix[2]
radius = propertymatrix[3]
teff = propertymatrix[4]
metal = propertymatrix[5]
age = propertymatrix[6]
dilution = propertymatrix[7]
durations = propertymatrix[8]
rpors = propertymatrix[9]
like = propertymatrix[-1]
alltotal = np.sum(np.log(like), axis=0)
allmaxllike = max(alltotal)
allllike = alltotal-allmaxllike
def get_aors(period, rho_star):
# [universal gravitational constant]^(1/3) * (day)^(2/3) * [solar density]^(1/3)* ( 1/(3*pi) )^(1/3)
aorsfactor = 4.206
aors = aorsfactor * period**(2./3.) * (rho_star*(4./3.*np.pi))**(1./3.)
if (aors < 1):
aors = np.nan
return aors
def compute_circular_edgeon_duration(period, rpors, aors):
prefactor = 1./np.pi
T0 = prefactor * period * np.arcsin( (1.+rpors)/aors )
return T0
######
kois=koilist
print "Nsystems=", len(kois)
print "koi loglike dur dur_e dur_circ Rearth"
preflist = []
d0list = []
dlist = []
radlist = []
for i in range(len(like)):
l = like[i]
ll = np.log(l)
ll = ll-max(ll)
sigma = np.sqrt(-ll*2)
#print koi[i]
plt.figure(1)
plt.semilogy(e, sigma)#, color='k')
plt.figure(2)
plt.plot(e, sigma)#, color='k')
#plt.plot(e, ll, color='k')
#if kois[i] % 1 != 0.01:
# print kois[i]
#if min(ll) < -2:
if np.sqrt(-min(ll)*2) > 3:
rho_stari = msun[i] / (4./3.*np.pi*rsun[i]**3.)
aorsi = get_aors(period[i], rho_stari)
print kois[i], np.sqrt(-min(ll)*2), durations[i], compute_circular_edgeon_duration(period[i], rpors[i], aorsi), radius[i]
preflist.append(ll[-1]-ll[0])
d0list.append(durations[i])
radlist.append(radius[i])
rho_stari = msun[i] / (4./3.*np.pi*rsun[i]**3.)
aorsi = get_aors(period[i], rho_stari)
dlist.append(compute_circular_edgeon_duration(period[i], rpors[i], aorsi))
prefsort = list(reversed(np.argsort(preflist)))
np.savetxt('koipreflist.txt', np.transpose([np.array(kois)[prefsort], np.array(preflist)[prefsort],
np.array(d0list)[prefsort], np.array(dlist)[prefsort], np.array(radlist)[prefsort]]), fmt=['%.2f', '%15.7f', '%15.7f', '%15.7f', '%15.7f'])
plt.figure(1)
plt.ylabel('Disfavored at $\sigma$ Level', fontsize=12)
plt.xlabel('$\sigma_e$', fontsize=12)
plt.xlim((0.,0.7))
plt.ylim((0.1, 8.))
plt.savefig('figs4/individual_singles.png')
plt.close()
plt.figure(2)
plt.ylabel('Disfavored at $\sigma$ Level', fontsize=12)
plt.xlabel('$\sigma_e$', fontsize=12)
plt.xlim((0.,0.7))
plt.ylim((0., 40.))
plt.savefig('figs4/individual_singles2.png')