# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + [markdown] id="D4yKHD9v9Sn_"
# # HIP 11915
#
# [The Solar Twin Planet Search II. A Jupiter twin around a solar twin](https://arxiv.org/abs/1507.03998)
#
# https://github.com/thiagofst/HIP11915
#
#

# + [markdown] id="2fWakgfvAVot"
# ### Google Drive
#

# + colab={"base_uri": "https://localhost:8080/"} id="cLPlr2yyocH6" executionInfo={"status": "ok", "timestamp": 1648820887102, "user_tz": 180, "elapsed": 2599, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="4085f57f-433e-4f3d-9f08-4f6be5f07866"
from google.colab import drive
drive.mount('/content/drive')

# + id="ucDInpLUpCOs" executionInfo={"status": "ok", "timestamp": 1648820887103, "user_tz": 180, "elapsed": 16, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
import os
os.chdir('/content/drive/MyDrive/SolarTwins/HIP11915')

# + colab={"base_uri": "https://localhost:8080/"} id="xvUlUzKoVxuV" executionInfo={"status": "ok", "timestamp": 1648820887105, "user_tz": 180, "elapsed": 16, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="cfa9f6f9-7276-4ee3-fd35-c51edaf05309"
# !ls

# + [markdown] id="U7tZO4SidHwT"
# ## Packages

# + id="_FAXq3TZbzpM" executionInfo={"status": "ok", "timestamp": 1648820892695, "user_tz": 180, "elapsed": 5600, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
# %matplotlib inline

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import arviz as az
import astropy.units as u

# + colab={"base_uri": "https://localhost:8080/"} id="0JjyVDm5dJI2" executionInfo={"status": "ok", "timestamp": 1648820902486, "user_tz": 180, "elapsed": 9798, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="a4d1a03b-eb47-440c-f091-e3b560816ba5"
# !pip install corner exoplanet astrobase

import corner
import exoplanet as xo
import astrobase.periodbase

# + id="79XWsB8nbXe3" executionInfo={"status": "ok", "timestamp": 1648820902489, "user_tz": 180, "elapsed": 17, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
from astropy.timeseries import LombScargle

# + colab={"base_uri": "https://localhost:8080/"} id="rwfU8uty-pqY" executionInfo={"status": "ok", "timestamp": 1648820907520, "user_tz": 180, "elapsed": 5045, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="6600dd53-a063-4160-f562-82e1b5d7fca3"
# !pip install radvel
import radvel
from radvel.plot import orbit_plots, mcmc_plots

# + id="LPTtFRyWl527" executionInfo={"status": "ok", "timestamp": 1648820907523, "user_tz": 180, "elapsed": 34, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
import copy
import pickle
import pdb
import configparser

# + id="lkXzrTDvu0gk" executionInfo={"status": "ok", "timestamp": 1648820907525, "user_tz": 180, "elapsed": 32, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
from astropy.timeseries import LombScargle

# + [markdown] id="DtMetGYodRmj"
# ## Data
#

# + colab={"base_uri": "https://localhost:8080/", "height": 427} id="8IDp00z6czWq" executionInfo={"status": "ok", "timestamp": 1648820907942, "user_tz": 180, "elapsed": 446, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="a5254439-ee84-482c-99f8-ab1b6c375973"
data = pd.read_csv('https://raw.githubusercontent.com/thiagofst/HIP11915/main/A/HIP11915_A.txt')
data.head(10)

# + id="QpMNzI4udUvm" executionInfo={"status": "ok", "timestamp": 1648820907943, "user_tz": 180, "elapsed": 11, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
x = data['time'].values
y = (data['mnvel'].values * (u.m/u.s))
yerr = (data['errvel'].values * (u.m/u.s))

# + id="UVjy_f5qc8he" executionInfo={"status": "ok", "timestamp": 1648820908350, "user_tz": 180, "elapsed": 415, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
fig, ax = plt.subplots(figsize = (12, 4))
ax.grid(alpha = 0.5)
ax.set(xlabel = 'Time / days', ylabel = 'Radial Velocity (%s)'%(y.unit))
ax.errorbar(x-x.min(), y.value - y.value.mean(), yerr = yerr.value, fmt = 'ok', ecolor = 'grey')

plt.tight_layout()
plt.savefig('HIP11915_rv.pdf')

# + colab={"base_uri": "https://localhost:8080/"} id="d0FZmKc2M4VK" executionInfo={"status": "ok", "timestamp": 1648820908351, "user_tz": 180, "elapsed": 14, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="2fc1c166-57a5-4801-e847-0bafbd38c365"
x.max()-x.min()

# + [markdown] id="q8oFDKGNfqaW"
# ## Lomb-Scargle

# + id="9s8MKk4MhdgE" executionInfo={"status": "ok", "timestamp": 1648820908352, "user_tz": 180, "elapsed": 12, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
phase = lambda t, period: np.remainder(t, period)/period

# + [markdown] id="SDS3y697utlW"
# ### `exoplanet`

# + id="_Rvns0FDfrox" executionInfo={"status": "ok", "timestamp": 1648820915712, "user_tz": 180, "elapsed": 7371, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
results = xo.estimators.lomb_scargle_estimator(
    x, 
    y.value,
    yerr.value, 
    max_peaks = 5, 
    min_period = 0.1, 
    max_period = 10000, 
    samples_per_peak = 50
)

# + colab={"base_uri": "https://localhost:8080/"} id="_TK7eRBlk8sO" executionInfo={"status": "ok", "timestamp": 1648820915713, "user_tz": 180, "elapsed": 24, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="3aa90f2a-bc77-4371-ce4e-b10f30271d80"
for i in range(5):
  peak = results['peaks'][i]
  print(peak['period'])

# + colab={"base_uri": "https://localhost:8080/"} id="ECQhfCk9mrxq" executionInfo={"status": "ok", "timestamp": 1648820915715, "user_tz": 180, "elapsed": 16, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="742a8ea5-3a72-4ae1-8d7b-7346c0dbcfd2"
peak = results['peaks'][1] # Main results "flag"

period_xo = peak['period'] # Orbital period
period_sigma_xo = peak['period_uncert'] # Orbital period uncertainty
print(f'P = {period_xo} +/- {period_sigma_xo} days')

freq_xo, power_xo = results['periodogram'] # Periodogram parameters
print(f'Number of frequencies: {len(freq_xo)}')

# + id="H0_VuklqmuVs" executionInfo={"status": "ok", "timestamp": 1648820917502, "user_tz": 180, "elapsed": 1799, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
fig, [ax1, ax2, ax3] = plt.subplots(nrows = 1, ncols = 3, figsize = (30, 5))

ax1.grid()
ax1.set(xlabel = 'Time / days', ylabel = 'Radial Velocity (%s)'%(y.unit), title = 'HIP 11915')
ax1.errorbar(x-x.min(), y.value - y.value.mean(), yerr = yerr.value, fmt = 'ok', ecolor = 'grey')

ax2.set_xscale('log')
ax2.set(xlabel = 'Frequency / days$^{-1}$', ylabel = r'$\hat{\Theta}_\mathrm{LS}$')
ax2.plot(freq_xo, power_xo, lw = 1, color = 'k')

ax3.grid()
ax3.set(xlabel = 'Phase (P = %.4f $\pm$ %.4f days)'%(period_xo, period_sigma_xo), ylabel = 'Radial Velocity (%s)'%(y.unit))
ax3.errorbar(np.concatenate([phase(x, period_xo), phase(x, period_xo)-1]), 
             np.concatenate([y.value - y.value.mean(), y.value - y.value.mean()]), 
             yerr = np.concatenate([yerr.value, yerr.value]), fmt = 'ok', ecolor = 'grey')

plt.tight_layout()
plt.savefig('HIP11915_exoplanet.pdf')

# + [markdown] id="1QT7J0w-uw-M"
# ### `astrobase`

# + colab={"base_uri": "https://localhost:8080/"} id="cFo8N9WiuwfI" executionInfo={"status": "ok", "timestamp": 1648821062872, "user_tz": 180, "elapsed": 145375, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="b992ff02-56a7-46e8-b68b-bd13558eeab1"
zgls = astrobase.periodbase.zgls.pgen_lsp(
    x, 
    y.value, 
    yerr.value, 
    magsarefluxes = False, 
    startp = 0.1, 
    endp = 10000, 
    stepsize = 1e-5, #1/(5 * (x.max() - x.min()))
    autofreq = False, 
    nbestpeaks = 4, 
    periodepsilon = 0.1, 
    sigclip = None, 
    nworkers = 4, 
    verbose = True
)

# + colab={"base_uri": "https://localhost:8080/"} id="PgEbmwIuvBDy" executionInfo={"status": "ok", "timestamp": 1648821062874, "user_tz": 180, "elapsed": 26, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="ca2281de-196f-4d19-c715-bc1e721cd74e"
print(zgls['nbestperiods'])
zgls

# + colab={"base_uri": "https://localhost:8080/"} id="ZGkEC9tExE7S" executionInfo={"status": "ok", "timestamp": 1648821062875, "user_tz": 180, "elapsed": 18, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="f10d092c-b3bc-4092-e415-3db4f36c41eb"
period_astrobase = zgls['nbestperiods'][1]
print(f'\nPeriod = {period_astrobase} days')

power_astrobase, periods_astrobase = zgls['lspvals'], zgls['periods'] # Periodogram

# + id="KzFUhbTHxKkD" executionInfo={"status": "ok", "timestamp": 1648821064159, "user_tz": 180, "elapsed": 1296, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
# Best result

fig, [ax1, ax2, ax3] = plt.subplots(nrows = 1, ncols = 3, figsize = (30, 5))

ax1.grid()
ax1.set(xlabel = 'Time / days', ylabel = 'Radial Velocity (%s)'%(y.unit), title = 'HIP 11915')
ax1.errorbar(x-x.min(), y.value - y.value.mean(), yerr = yerr.value, fmt = 'ok', ecolor = 'grey')

ax2.set_xscale('log')
ax2.set(xlabel = 'Frequency / days$^{-1}$', ylabel = r'$\hat{\Theta}_\mathrm{LS}$')
ax2.plot(1/periods_astrobase, power_astrobase, lw = 1, color = 'k')

# ax2_ = ax2.twinx()
# ax2_.set_ylim(tuple(lomb_scargle_BIC(ax1.get_ylim(), y.value, yerr.value)))
# ax2_.set_ylabel(r'$\Delta BIC$')

ax3.grid()
ax3.set(xlabel = 'Phase (P = %.4f days)'%(period_astrobase), ylabel = 'Radial Velocity (%s)'%(y.unit))
ax3.errorbar(np.concatenate([phase(x, period_astrobase), phase(x, period_astrobase)-1]), 
             np.concatenate([y.value - y.value.mean(), y.value - y.value.mean()]), 
             yerr = np.concatenate([yerr.value, yerr.value]), fmt = 'ok', ecolor = 'grey')

plt.tight_layout()
plt.savefig('HIP11915_astrobase.pdf')

# + [markdown] id="sAf_8RS5dGwY"
# ### `astropy` (low-resolution)

# + colab={"base_uri": "https://localhost:8080/"} id="XTL9e0bndIMM" executionInfo={"status": "ok", "timestamp": 1648821064161, "user_tz": 180, "elapsed": 14, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="6bb37e4e-5131-4feb-d48e-2ca4aea01ab5"
ls = LombScargle(x, y.value, yerr.value)
freq_astropy, power_astropy = ls.autopower(minimum_frequency = 1./10000, maximum_frequency = 1./100, samples_per_peak = 100, method = 'slow')
period_astropy = 1/freq_astropy[np.argmax(power_astropy)]
print(f'Period: {period_astropy} days')
print(f'Number of frequencies: {len(freq_astropy)}')

# + colab={"base_uri": "https://localhost:8080/"} id="pLKaS1W_nLFN" executionInfo={"status": "ok", "timestamp": 1648821066508, "user_tz": 180, "elapsed": 2356, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="596cbc0e-818a-4e49-df01-3f21dd9b795e"
print('The false-alarm probability (FAP) is the probability \n of measuring a peak of a given height (or higher)\n conditioned on the assumption that the data consists \n of Gaussian noise with no periodic component\n')

methods = ['davies', 'naive', 'baluev', 'bootstrap']

for k in np.arange(len(methods)):
  fap = ls.false_alarm_probability(power_astropy.max(), method = str(methods[k]))
  print(f'{str(methods[k])}: {fap}')

# + colab={"base_uri": "https://localhost:8080/"} id="80pN02hkdLN6" executionInfo={"status": "ok", "timestamp": 1648821070367, "user_tz": 180, "elapsed": 3862, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="c04cfe8b-12de-44bd-8ae6-f738bcd835c4"
probabilities = [0.1, 0.05, 0.01, 0.001]

for k in np.arange(len(methods)):
  probs = ls.false_alarm_level(probabilities, 
                               method = str(methods[k]),
                               samples_per_peak = 50, 
                               nyquist_factor = 5)
  print(f'{str(methods[k])}: {probs}')

# This tells us that to attain a 10% (=0.1) FAP requires the highest periodogram peak to be approximately 0.12; 
# 5% (=0.05) requires 0.13, and 1% (=0.01) requires 0.15.
# Values can change

# + id="axZxABwRdLuL" executionInfo={"status": "ok", "timestamp": 1648821079997, "user_tz": 180, "elapsed": 9636, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
fig, [ax1, ax2, ax3] = plt.subplots(nrows = 1, ncols = 3, figsize = (30, 5))

ax1.grid()
ax1.set(xlabel = 'Time / days', ylabel = 'Radial Velocity (%s)'%(y.unit), title = 'HIP 11915')
ax1.errorbar(x-x.min(), y.value - y.value.mean(), yerr = yerr.value, fmt = 'ok', ecolor = 'grey')

ax2.set_xscale('log')
ax2.set(xlabel = 'Frequency / days$^{-1}$', ylabel = r'$\hat{\Theta}_\mathrm{LS}$')
ax2.plot(freq_astropy, power_astropy, lw = 1, color = 'k')

methods = ['bootstrap']

for k in np.arange(len(methods)):
  for j in np.arange(len(probabilities)):
    probs = ls.false_alarm_level(probabilities, 
                                 method = str(methods[k]),
                                 samples_per_peak = 50, 
                                 nyquist_factor = 5)
    
    ax2.axhline(probs[j], lw = 0.6, ls = 'dashdot')#, label = f'{methods[k]}: {float(probabilities[j]*100)}%')
    # ax2.legend(loc = 'upper right', shadow = True)

ax3.grid()
ax3.set(xlabel = 'Phase (P = %.4f days)'%(period_astropy), ylabel = 'Radial Velocity (%s)'%(y.unit))
ax3.errorbar(np.concatenate([phase(x, period_astropy), phase(x, period_astropy)-1]), 
             np.concatenate([y.value - y.value.mean(), y.value - y.value.mean()]), 
             yerr = np.concatenate([yerr.value, yerr.value]), fmt = 'ok', ecolor = 'grey')

plt.tight_layout()
plt.savefig('HIP11915_astropy.pdf')

# + [markdown] id="4TWd0YKlxlOw"
# ### Results

# + colab={"base_uri": "https://localhost:8080/"} id="mKNN-F0DxtzQ" executionInfo={"status": "ok", "timestamp": 1648821080000, "user_tz": 180, "elapsed": 27, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="422122a5-dc1d-429a-c9bf-e67d3d8e4d6e"
print(f'Resulting orbital periods (days): \n\nastrobase: {period_astrobase}d\nexoplanet: {period_xo}+/-{period_sigma_xo}\nastropy: {period_astropy} (low-resolution)')

# + colab={"base_uri": "https://localhost:8080/"} id="sH9kShJExaZQ" executionInfo={"status": "ok", "timestamp": 1648821080002, "user_tz": 180, "elapsed": 23, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="a3f8ae43-018a-4d5e-df98-8995ef45a2c1"
# astrobase

periods = float('%.7f'%(period_astrobase))
# period_errs = period_sigma_xo
t0s = x[np.argmax(y)]
t0_errs = 0.3
Ks = np.abs(xo.estimate_semi_amplitude(periods, x, y, yerr, t0s=t0s))
print(Ks, "m/s")

msini_xo = np.abs(xo.estimate_minimum_mass(periods, x, y, yerr, t0s, m_star = 0.993))
print(msini_xo)

# exoplanet

periods = float('%.7f'%(period_xo))
# period_errs = period_sigma_xo
t0s = x[np.argmax(y)]
t0_errs = 0.3
Ks = np.abs(xo.estimate_semi_amplitude(periods, x, y, yerr, t0s=t0s))
print(Ks, "m/s")

msini_xo = np.abs(xo.estimate_minimum_mass(periods, x, y, yerr, t0s, m_star = 0.993))
print(msini_xo)

# + colab={"base_uri": "https://localhost:8080/"} id="rgnuW9iMoiXS" executionInfo={"status": "ok", "timestamp": 1648821080005, "user_tz": 180, "elapsed": 22, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="9a7c1386-8db4-45c9-d011-53dac986520b"
# astropy

periods = float('%.7f'%(period_astropy))
# period_errs = period_sigma_xo
t0s = x[np.argmax(y)]
t0_errs = 0.3
Ks = np.abs(xo.estimate_semi_amplitude(periods, x, y, yerr, t0s=t0s))
print(Ks, "m/s")

msini_xo = np.abs(xo.estimate_minimum_mass(periods, x, y, yerr, t0s, m_star = 0.993))
print(msini_xo)

# + [markdown] id="9NL9NIRCVVcK"
# ## Analysis of Variance

# + colab={"base_uri": "https://localhost:8080/"} id="qYg6RlpaVa7Q" executionInfo={"status": "ok", "timestamp": 1648821138672, "user_tz": 180, "elapsed": 58685, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="925c0601-0895-428f-c5b0-f997ee225beb"
aov = astrobase.periodbase.smav.aovhm_periodfind(
    x, 
    y.value, 
    yerr.value, 
    magsarefluxes = False, 
    startp = 1.1, 
    endp = 10000, 
    stepsize = 1e-5, #1/(5 * (x.max() - x.min()))
    autofreq = False, 
    normalize = True, 
    nharmonics = 6,
    nbestpeaks = 4, 
    periodepsilon = 0.1, 
    sigclip = None, 
    nworkers = 4, 
    verbose = True
)

# + colab={"base_uri": "https://localhost:8080/"} id="Hi3AQyYRVdxu" executionInfo={"status": "ok", "timestamp": 1648821138673, "user_tz": 180, "elapsed": 15, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="61419522-f5dc-4a39-f37c-7f8a5e50cd6b"
nbestperiods_aov = aov['nbestperiods']
aov

# + colab={"base_uri": "https://localhost:8080/"} id="GY2S16PRVfaS" executionInfo={"status": "ok", "timestamp": 1648821138674, "user_tz": 180, "elapsed": 9, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="55f9bfa3-f729-4db2-e390-0ee973c7a042"
period_aov = aov['bestperiod']
print(f'\nPeriod = {period_aov} days')
print(nbestperiods_aov)

power_aov, periods_aov = aov['lspvals'], aov['periods'] # Periodogram

# + id="y76m1RxPVhXC" executionInfo={"status": "ok", "timestamp": 1648821146788, "user_tz": 180, "elapsed": 8121, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
for i in range(len(nbestperiods_aov)):
  fig, [ax1, ax2, ax3] = plt.subplots(nrows = 1, ncols = 3, figsize = (30, 5))

  ax1.grid()
  ax1.set(xlabel = 'Time / days', ylabel = 'Radial Velocity (%s)'%(y.unit), title = 'HIP 11915')
  ax1.errorbar(x-x.min(), y.value - y.value.mean(), yerr = yerr.value, fmt = 'ok', ecolor = 'grey')

  ax2.set_xscale('log')
  ax2.set(xlabel = 'Frequencies / d$^{-1}$', ylabel = r'$\hat{\Theta}_\mathrm{MHAoV}$')
  ax2.plot(1./periods_aov, power_aov, lw = 1, color = 'k')

  ax3.grid()
  ax3.set(xlabel = 'Phase (P = %.4f days)'%(nbestperiods_aov[i]), ylabel = 'Radial Velocity (%s)'%(y.unit))
  ax3.errorbar(np.concatenate([phase(x, nbestperiods_aov[i]), phase(x, nbestperiods_aov[i])-1]), 
              np.concatenate([y.value - y.value.mean(), y.value - y.value.mean()]), 
              yerr = np.concatenate([yerr.value, yerr.value]), fmt = 'ok', ecolor = 'grey')

  plt.tight_layout()
  plt.savefig(f'HIP11915_aov_{nbestperiods_aov[i]}.pdf')

# + [markdown] id="IaB3YcuQVmBr"
# ## Phase Dispersion Minimisation

# + colab={"base_uri": "https://localhost:8080/"} id="-6F6ELG4VnzR" executionInfo={"status": "ok", "timestamp": 1648821172075, "user_tz": 180, "elapsed": 25306, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="25fc85e3-c479-43fb-8ef9-262cb8380d71"
pdm = astrobase.periodbase.spdm.stellingwerf_pdm(
    x, 
    y.value, 
    yerr.value, 
    magsarefluxes = False, 
    startp = 1.1, 
    endp = 10000, 
    stepsize = 1e-5, #1/(5 * (x.max() - x.min()))
    autofreq = False, 
    normalize = True,
    nbestpeaks = 4, 
    periodepsilon = 0.1, 
    sigclip = None, 
    nworkers = 4, 
    verbose = True
)

# + colab={"base_uri": "https://localhost:8080/"} id="5O8m0UEXVpHl" executionInfo={"status": "ok", "timestamp": 1648821172077, "user_tz": 180, "elapsed": 14, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="353914d7-b5b8-41eb-a0f4-1d46e01d942f"
nbestperiods_pdm = pdm['nbestperiods']
pdm

period_pdm = pdm['bestperiod']
print(f'\nPeriod = {period_pdm} days')
print(nbestperiods_pdm)
power_pdm, periods_pdm = pdm['lspvals'], pdm['periods'] # Periodogram

# + id="jIExwbttXEN4" executionInfo={"status": "ok", "timestamp": 1648821175525, "user_tz": 180, "elapsed": 3456, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
for i in range(len(nbestperiods_pdm)):
  fig, [ax1, ax2, ax3] = plt.subplots(nrows = 1, ncols = 3, figsize = (30, 5))

  ax1.grid()
  ax1.set(xlabel = 'Time / days', ylabel = 'Radial Velocity (%s)'%(y.unit), title = 'HIP 11915')
  ax1.errorbar(x-x.min(), y.value - y.value.mean(), yerr = yerr.value, fmt = 'ok', ecolor = 'grey')

  ax2.set_xscale('log')
  ax2.set(xlabel = 'Periods / days', ylabel = r'$\hat{\Theta}_\mathrm{PDM}$')
  ax2.plot(periods_pdm, power_pdm, lw = 1, color = 'k')

  ax3.grid()
  ax3.set(xlabel = 'Phase (P = %.4f days)'%(nbestperiods_pdm[i]), ylabel = 'Radial Velocity (%s)'%(y.unit))
  ax3.errorbar(np.concatenate([phase(x, nbestperiods_pdm[i]), phase(x, nbestperiods_pdm[i])-1]), 
              np.concatenate([y.value - y.value.mean(), y.value - y.value.mean()]), 
              yerr = np.concatenate([yerr.value, yerr.value]), fmt = 'ok', ecolor = 'grey')

  plt.tight_layout()
  plt.savefig(f'HIP11915_pdm_{nbestperiods_pdm[i]}.pdf')

# + [markdown] id="vBoCvoA8vwlI"
# ## RadVel
#
# Read the docs: https://radvel.readthedocs.io/en/latest/
#
# **Initial files**:
#
# - `.py` $-$ Keplerian fit configuration
# - `.txt` $-$ RV curve (It is alocated on GitHub)

# + colab={"base_uri": "https://localhost:8080/"} id="S3gxNzU_0IDt" executionInfo={"status": "ok", "timestamp": 1648821175527, "user_tz": 180, "elapsed": 19, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="259eb239-948e-4c19-8b4a-bd409ae99140"
os.chdir('/content/drive/MyDrive/SolarTwins/HIP11915')
# !ls

# + [markdown] id="lStV1lTP0_kf"
# #### Maximum likelihood fit

# + colab={"base_uri": "https://localhost:8080/"} id="QXu3A1LK0Kzb" executionInfo={"status": "ok", "timestamp": 1648821181521, "user_tz": 180, "elapsed": 6004, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="84e62b44-e40e-4cb1-b993-d7800bada7e7"
# !radvel fit -s HIP11915.py

# + [markdown] id="FoL_EMrz1FtQ"
# #### Plot the best-fit solution
#
# It creates a primary `_multipanel.pdf` file that will be upgraded with MCMC results

# + colab={"base_uri": "https://localhost:8080/"} id="SZV2O68A0TxJ" executionInfo={"status": "ok", "timestamp": 1648821186153, "user_tz": 180, "elapsed": 4640, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="ab0c73af-5b94-4e61-bf17-dacd6ea00fe7"
# !radvel plot -t rv -s HIP11915.py

# + [markdown] id="iyvF5MXA1NRX"
# #### MCMC

# + colab={"base_uri": "https://localhost:8080/"} id="q70hIkgz0ioz" executionInfo={"status": "ok", "timestamp": 1648821487033, "user_tz": 180, "elapsed": 300885, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="501f1f6b-3058-4f32-dd9c-a1de654fdefd"
# !radvel mcmc -s HIP11915.py

# + [markdown] id="Puuwj2yM-t9U"
# #### Update RV plots with MCMC

# + colab={"base_uri": "https://localhost:8080/"} id="hN9ODDip-S-W" executionInfo={"status": "ok", "timestamp": 1648821490290, "user_tz": 180, "elapsed": 3272, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="f42759c8-4d5e-4fb9-aaaf-9a3e0bdfa71b"
# !pip install celerite

# + colab={"base_uri": "https://localhost:8080/"} id="WA6-uyhQ-ydc" executionInfo={"status": "ok", "timestamp": 1648821559815, "user_tz": 180, "elapsed": 69539, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="1af732e7-97d0-40d3-ce57-8fd95a2d0818"
# !radvel plot -t rv corner trend -s HIP11915.py

# + [markdown] id="NZGhNJ8c_wST"
# #### Model comparison
#
# It performs model comparison testing to eliminate different sets of planets, their eccentricities, and RV trends. If this run well, a new table will be included in the summary report.

# + colab={"base_uri": "https://localhost:8080/"} id="VbLcUjqx_0yH" executionInfo={"status": "ok", "timestamp": 1648821566577, "user_tz": 180, "elapsed": 6774, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="486dec08-a1a3-4635-db60-33f565cf0c63"
# !radvel ic -t nplanets e trend -s HIP11915.py

# + [markdown] id="kMii0CzhA2aL"
# **Important** Derives physical parameters for the planetary system: Minimum mass (M$_{p}\sin{i}$) and orbital separation $a$. 

# + colab={"base_uri": "https://localhost:8080/"} id="nSj9YISxBELv" executionInfo={"status": "ok", "timestamp": 1648821573446, "user_tz": 180, "elapsed": 6876, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="dbee1bc0-9a11-4330-b44a-1514314758fd"
# !radvel derive -s HIP11915.py

# + [markdown] id="HkWivxJsBH1j"
# Corner plots for physical parameters

# + colab={"base_uri": "https://localhost:8080/"} id="fMElMVR3BJuQ" executionInfo={"status": "ok", "timestamp": 1648821579841, "user_tz": 180, "elapsed": 6401, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="a04438eb-04c0-404c-c44e-720564e35128"
# !radvel plot -t derived -s HIP11915.py

# + [markdown] id="bV9quDUpBODl"
# #### LaTeX files
#
# The main file is **_results.tex**
#
# `!radvel report...` is not working properly, it just creates a `_results.tex` file, but not the final `.pdf` file as expected. 
# (Not working $=$ It is taking so much time to run as I expected for just compile LaTeX files). 
#
# The next cell creates `.tex` table files: *parameters, priors, models comparisons* and *derived parameters*

# + colab={"base_uri": "https://localhost:8080/"} id="Ad7nOc8wBsW5" executionInfo={"status": "ok", "timestamp": 1648821586402, "user_tz": 180, "elapsed": 6571, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="0f8137bf-4aba-4a53-97eb-41e72742f4e3"
# !radvel table -t params priors ic_compare derived -s HIP11915.py

# + colab={"base_uri": "https://localhost:8080/", "height": 1000} id="6KKnNJurBy54" outputId="fb95d62d-82b3-49b9-c6c5-43c35049bdd2" executionInfo={"status": "ok", "timestamp": 1648821652728, "user_tz": 180, "elapsed": 66332, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
# !pip install folium==0.2.1
# !pip install pdflatex
# !apt install texlive-xetex texlive-fonts-recommended texlive-generic-recommended

# + id="nPk6FcCYCAcY" executionInfo={"status": "ok", "timestamp": 1648821652730, "user_tz": 180, "elapsed": 51, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
import pdflatex

# + id="Cwk03BwICCPg" executionInfo={"status": "ok", "timestamp": 1648821652733, "user_tz": 180, "elapsed": 49, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
# # !radvel report -s HIP11915.py

# + [markdown] id="ztfQHf3nl8Qm"
# ## Final results

# + [markdown] id="pWQQtfAmj3sK"
# ### Data
#
# - MCMC summary
# - Derived Keplerian quantities

# + id="Rn7dS0i5j878" executionInfo={"status": "ok", "timestamp": 1648821652735, "user_tz": 180, "elapsed": 48, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="927b9230-04c0-46ed-9181-fc4273368661" colab={"base_uri": "https://localhost:8080/", "height": 225}
mcmc_summary = pd.read_csv('HIP11915/HIP11915_post_summary.csv')
mcmc_summary.head()

# + [markdown] id="QQQd_d9Mwi6m"
# Multiply `k1` by 1000 and take its absolute value if originally the RV values were in km/s

# + id="AGOvEyoWmGAz" executionInfo={"status": "ok", "timestamp": 1648821652736, "user_tz": 180, "elapsed": 45, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="7c7e472d-2c1c-4a1d-baea-96b8e3fed046" colab={"base_uri": "https://localhost:8080/", "height": 225}
derived = pd.read_csv('HIP11915/HIP11915_derived_quantiles.csv')
derived.head()

# + [markdown] id="Dk8oj2lMwxV7"
# Convert `mpsini1` from Earth units to Jupiter units and multiply by 1000 if originally the RV values were in km/s

# + id="XpJZEVcRww1M" executionInfo={"status": "ok", "timestamp": 1648821652738, "user_tz": 180, "elapsed": 40, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
derived['mpsini1'] = np.abs((derived['mpsini1'])*(u.Mearth).to(u.Mjupiter))
derived['mpsini2'] = np.abs((derived['mpsini2'])*(u.Mearth).to(u.Mjupiter))

# + colab={"base_uri": "https://localhost:8080/", "height": 319} id="tMy5HOpsmuc0" executionInfo={"status": "ok", "timestamp": 1648821653829, "user_tz": 180, "elapsed": 1128, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="ae765769-759c-44ea-831e-97557437fd31"
derived_mcmc = pd.read_csv('HIP11915/HIP11915_derived.csv.bz2')

derived_mcmc['mpsini1'] = np.abs((derived_mcmc['mpsini1'])*(u.Mearth).to(u.Mjupiter))
derived_mcmc['mpsini2'] = np.abs((derived_mcmc['mpsini2'])*(u.Mearth).to(u.Mjupiter))

derived_mcmc.head(6)

# + colab={"base_uri": "https://localhost:8080/", "height": 288} id="0tsP2bZjq0gy" executionInfo={"status": "ok", "timestamp": 1648821655649, "user_tz": 180, "elapsed": 1851, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="7b85a3ba-0954-4694-a298-15e0c10a4a82"
mcmc_chains = pd.read_csv('HIP11915/HIP11915_chains.csv.bz2')
mcmc_chains.head(5)

# + colab={"base_uri": "https://localhost:8080/"} id="VUMgeNqz8d8C" executionInfo={"status": "ok", "timestamp": 1648821655650, "user_tz": 180, "elapsed": 17, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="4f234b33-7b20-43f8-e1f8-36a9d82367af"
print(mcmc_summary.per1[1])
print(mcmc_summary.per2[1])

# + [markdown] id="5CE5IVIGqYdj"
# ### Plots

# + [markdown] id="dr0mONfvqbSh"
# ##### Planet's mass and orbital separation

# + colab={"base_uri": "https://localhost:8080/"} id="esmE8paGmNtV" executionInfo={"status": "ok", "timestamp": 1648821658203, "user_tz": 180, "elapsed": 2564, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="7bf942de-7e5d-403a-91cd-039ec64fe9c7"
_ = corner.corner(derived_mcmc, 
              var_names = ['mpsini1', 'a1', 'mpsini2', 'a2'], 
              quantiles = [0.16, 0.5, 0.84], 
              truths = [derived.mpsini1[1], derived.a1[1], derived.mpsini2[1], derived.a2[1]], 
              show_titles = True, 
              title_fmt = '.3f', 
              labels = [r'$m_{b}\sin{i}$ [M$_{Jup}$]', r'a$_{b}$ [au]', r'$m_{c}\sin{i}$ [M$_{Jup}$]', r'a$_{c}$ [au]'],
              )

plt.savefig('HIP11915_corner_derived_pars.pdf')

# + [markdown] id="umwM5xclq92k"
# ##### Keplerian parameters

# + colab={"base_uri": "https://localhost:8080/"} id="8E3AqyTyk1kQ" executionInfo={"status": "ok", "timestamp": 1648821658204, "user_tz": 180, "elapsed": 20, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="b647aaa8-41f1-4d0e-f722-536cf9a3f30c"
mcmc_chains.keys()

# + colab={"base_uri": "https://localhost:8080/", "height": 73, "output_embedded_package_id": "1dCkfOvWpaSo-i_ZD6vSI13vk2Q_CeeWd"} id="hujNG5Vapnt0" executionInfo={"status": "ok", "timestamp": 1648821713566, "user_tz": 180, "elapsed": 55370, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="79f0ba65-0321-4ec1-d1ab-c9a22fc5513a"
# %matplotlib inline

_ = corner.corner(mcmc_chains, 
                  var_names = ['per1', 'k1', 'secosw1', 'sesinw1', 'tc1', 
                               'per2', 'k2', 'secosw2', 'sesinw2', 'tc2', 
                               'gamma_HARPS-B', 
                               'gamma_HARPS-A', 
                               'jit_HARPS-B', 
                               'jit_HARPS-A'], 
                  quantiles = [0.16, 0.50, 0.84], 
                  truths = [mcmc_summary.per1[1], 
                            mcmc_summary.k1[1],
                            mcmc_summary.secosw1[1],
                            mcmc_summary.sesinw1[1],
                            mcmc_summary.tc1[1], 
                            mcmc_summary.per2[1], 
                            mcmc_summary.k2[1],
                            mcmc_summary.secosw2[1],
                            mcmc_summary.sesinw2[1],
                            mcmc_summary.tc2[1], 
                            np.nan, np.nan, np.nan, np.nan], 
                  show_titles = True, 
                  title_fmt = '.5f',
                  labels = [r'P$_{b}$ [d]', 
                            r'K$_{b}$ (m/s)', 
                            r'[$\sqrt{e}\cos{\omega}$]$_{b}$', 
                            r'[$\sqrt{e}\sin{\omega}$]$_{b}$', 
                            r'T$_{b}$', 
                            r'P$_{c}$ [d]', 
                            r'K$_{c}$ (m/s)', 
                            r'[$\sqrt{e}\cos{\omega}$]$_{c}$', 
                            r'[$\sqrt{e}\sin{\omega}$]$_{c}$', 
                            r'T$_{c}$', 
                            r'$\gamma_{\rm H-B}$', 
                            r'$\gamma_{\rm H-A}$', 
                            r'$\sigma_{\rm H-b}$', 
                            r'$\sigma_{\rm H-A}$']
                  )

plt.tight_layout()
plt.savefig('HIP11915_corner.pdf')


# + [markdown] id="AfLd9B3z0FdS"
# ##### Phase diagram

# + id="IQ1lmVcRzlYZ" executionInfo={"status": "ok", "timestamp": 1648821713568, "user_tz": 180, "elapsed": 38, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
def load_status(statfile):
    config = configparser.RawConfigParser()
    gl = config.read(statfile)  

    return config

status = load_status('HIP11915/HIP11915_radvel.stat')

# + colab={"base_uri": "https://localhost:8080/"} id="UL7guJlEhjvu" executionInfo={"status": "ok", "timestamp": 1648821713573, "user_tz": 180, "elapsed": 40, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="bd0939c4-8e0f-4386-f81a-065e3e1f9a95"
post = radvel.posterior.load('HIP11915/HIP11915_post_obj.pkl')

print(post)

# + [markdown] id="H_Jk2V3i6u57"
# Positive K

# + id="cU8YHjdy2dLy" executionInfo={"status": "ok", "timestamp": 1648821713578, "user_tz": 180, "elapsed": 35, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
post.medparams['k1'] = np.abs(post.medparams['k1'])
post.uparams['k1'] = post.uparams['k1']

# + id="VEYzU4eQL9TN" executionInfo={"status": "ok", "timestamp": 1648821713581, "user_tz": 180, "elapsed": 36, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}}
post.medparams['k2'] = np.abs(post.medparams['k2'])
post.uparams['k2'] = post.uparams['k2']

# + colab={"base_uri": "https://localhost:8080/"} id="KTHFwfFj7dzO" executionInfo={"status": "ok", "timestamp": 1648821713587, "user_tz": 180, "elapsed": 40, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="9a41b3fe-3b0c-4dff-d9d3-90f76fe4ce06"
print(post.derived)

# + colab={"base_uri": "https://localhost:8080/", "height": 166} id="OcDJJr7NL6Dm" executionInfo={"status": "ok", "timestamp": 1648821715844, "user_tz": 180, "elapsed": 2294, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="b2c32bf5-befc-43e0-95db-0c361771a048"
RVPlot = orbit_plots.MultipanelPlot(post, 
                                    uparams = post.uparams,
                                    legend_kwargs={'loc': 'upper right'}, 
                                    status = status, 
                                    nobin = True, 
                                    highlight_last = False,
                                    show_rms = True,
                                    saveplot = 'HIP11915_rv_multipanel.png')

RVPlot.plot_multipanel()

# + colab={"base_uri": "https://localhost:8080/", "height": 166} id="Wq0LbaY7Ybjq" executionInfo={"status": "ok", "timestamp": 1648821719193, "user_tz": 180, "elapsed": 3356, "user": {"displayName": "<NAME>", "userId": "09108111230663562541"}} outputId="73f4635f-6eaa-430e-c39e-cd93e4b6d966"
RVPlot = orbit_plots.MultipanelPlot(post, 
                                    uparams = post.uparams,
                                    legend_kwargs={'loc': 'upper right'}, 
                                    status = status, 
                                    nobin = True, 
                                    highlight_last = False,
                                    show_rms = True,
                                    saveplot = 'HIP11915_rv_multipanel.pdf')

RVPlot.plot_multipanel()

StarcoderdataJupyterScriptsDedupFiltered