Sim_Two/sim_two.py at main · StellarTeamUFC/Sim_Two · GitHub

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
# -*- coding: utf-8 -*-
"""Sim_Two.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1EasJjhr9O1DjzrBQhB2D4wsnG9UX-ysJ

# Libraries
"""

!pip install fbm

from scipy.optimize import fsolve
import numpy as np
from fbm import fgn

import matplotlib.pyplot as plt

"""# Variables"""

''' Variables Declaration '''
# Solar parameters
M_sun = 1.989 * 10**30                          # Sun's mass in kg
M_earth = 5.972 * 10**24
M_jup = 1.899 * 10**27

# Input Parameters
P_orb = 50                                      # Orbital period in days
ai = 90                                         # Inclination angle of star (ranged over 0 and 90 degrees)
omega = (np.pi/180)*(340)                       # Aargument of pericenter in true orbital plane (rad)
excent = 0.3                                   # Eccentricity
tp = 10

np_ = 0.5                                         # Planet mass fraction (in Jupter mass unit)
M_p = np_ * M_jup                               # Planet mass

ns = 1                                          # Stellar mass fraction (in Sun mass unit)
M_star = ns * M_sun                             # Stellar mass

"""# Execution"""

# Eq Torres
K = 203.29 * ((M_p * np.sin(np.pi/180*ai))/M_jup) * ((M_sun/(M_star + M_p))**(2/3)) * ((P_orb)**(-1/3)) * (1/np.sqrt(1-excent**2))

cadence = 8
times = np.linspace(1, 100, 100)
range_ = P_orb/len(times)
times *= range_

def transcendental(Et, excent, eqLeft):
  return Et - excent*np.cos(Et) - eqLeft

EtS = []

for time in times:
  eqLeft = (2*np.pi/P_orb) * (time-tp)                                  # Mass function
  res = fsolve(lambda x: transcendental(x, excent, eqLeft), x0=-10)
  if time == 20:
    print(res[0])
  #print(res[0], '\n')
  EtS.append(res[0])

ni_t = np.arccos((np.cos(EtS) - excent) / (1 - excent*np.cos(EtS)))

#times = times + 1
#times = P_orb/times

RV = K * (np.cos(omega + ni_t) + excent*np.cos(omega))

tam = 16

def plotar(times, y):
  plt.rcParams["figure.figsize"] = (16,8)
  plt.plot(times, y, '.k')
  plt.xlabel('Time (days)', fontsize = tam+5)
  plt.xticks(fontsize = tam)
  plt.ylabel('Radial Velocity (m/s)', fontsize = tam+5)
  plt.yticks(fontsize = tam)
  plt.show()

plotar(times, RV)

"""# Noise"""

''' Noisy time series '''

stdnoise = 5
fs = len(times)
noise = fgn(n=fs, hurst=0.0003)

RV_noise = RV + noise

plotar(times, RV_noise)

# Remove periodically points (2 in each 3 points)

RV_noise_red = RV_noise[0::3]
times_red = times[0::3]/24

plotar(times_red, RV_noise_red)

x = np.linspace(1,10,10)
print(x)
print(x[0::3])