/
simul_1647.py
106 lines (83 loc) · 2.55 KB
/
simul_1647.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
import numpy as np
import matplotlib.pyplot as plt
transit = 11.259
planet1 = (1 - 0.8016, -0.017, 0.016, 134.737) # (depth, drop_start, drop_end, start)
planet2_start = -0.46
planet2_end= -0.4348
planet2_mean = (planet2_end + planet2_start) / 2
planet2 = (1 - 0.8326, planet2_start, planet2_end, planet1[3] + (1+planet2_mean) * transit)
print(planet2_mean)
print(planet2[3])
def calc_phase(start, period, value):
diff = value - start
cycli = int(diff / period)
part = diff - cycli * period
phase = part / period
if phase < 0:
phase += 1
#print("HUH: ", start, period, value)
if phase > 0.5:
return -(1-phase)
else:
return phase
print(calc_phase(120,10,125) == 0.5)
print(calc_phase(120,10,126) == -0.4)
print(calc_phase(120,10,124) == 0.4)
print(calc_phase(120,10,115) == 0.5)
print(calc_phase(120,10,116) == -0.4)
print(abs(calc_phase(134.737, 11.259, 25.3852) - 0.287609912) < 0.001)
def calculate_planet_drop(planet, x):
phase = calc_phase(planet1[3], transit, x)
if planet[1] <= phase <= planet[2]:
mean = (planet[1] + planet[2]) / 2
return (1-(abs(mean-phase)/(mean-planet[1]))) * planet[0]
else:
return 0
step = 0.01
start = min(planet1[2], planet2[2])
end = 500
light = np.ones(int((end-start+step)/step))
light_err = np.zeros(int((end-start+step)/step))
time = np.arange(start, end, step, dtype=np.float64)
for i in range(len(light)):
x = start + i * step
light[i] -= calculate_planet_drop(planet1, x)
light[i] -= calculate_planet_drop(planet2, x)
import lightkurve as lk
lc = lk.LightCurve(time=time, flux=light, flux_err=light_err)
#plt.plot(light)
#plt.show()
"""
folded = lc.fold(period=transit, t0=plan²oet1[3])
folded.bin().scatter()
folded.plot_river()
plt.show()
"""
"""
folded = lc.fold(period=0.94, t0=260.360)
folded.bin().scatter()
folded.plot_river()
"""
folded = lc.fold(period=5.630, t0=264.200)
folded.bin().scatter()
folded.plot_river()
plt.show()
"""
def fold(light, fold_start, fold_period):
res_x = []
res_y = []
for i in range(len(light)):
x = start + i * step
phase = calc_phase(fold_start, fold_period, x)
res_x.append(phase)
res_y.append(light[i])
return (res_x, res_y)
def fold_plot(light, fold_start, fold_period):
folded = fold(light, fold_start, fold_period)
plt.scatter(folded[0], folded[1])
plt.show()
return folded
default_fold = fold_plot(light, planet1[3], transit)
#fold_plot(light, 260.950, 0.750)
fold_plot(light, 269.850, 10.540)
"""