# Path: dart_board/posterior.py
# def ln_posterior(x, dart):
# def ln_likelihood(x, dart):
# def posterior_properties(x, output, dart):
#     def get_error_from_kwargs(param, **kwargs):
# def calculate_ADAF_efficiency(Mdot, Mdot_edd, delta=0.1):
# def calculate_L_x(M1, mdot, k1):
# def calc_prob_from_mass_function(M1, M2, f_obs, f_err):
#     def model_h(M1, M2, f):
#     def func_integrand(f, f_obs, f_err, M1, M2):
#     def calc_f(f_obs, f_err, M1, M2):
# def check_output(output, binary_type):
# def get_theta_proj(ra, dec, ra_b, dec_b):
# def get_dtheta_dalpha(alpha, delta, alpha_b, delta_b):
# def get_dtheta_ddelta(alpha, delta, alpha_b, delta_b):
# def get_domega_dalpha(alpha, delta, alpha_b, delta_b):
# def get_domega_ddelta(alpha, delta, alpha_b, delta_b):
# def get_J_coor(alpha, delta, alpha_b, delta_b):
#     M1 = np.exp(ln_M1)
#     M2 = np.exp(ln_M2)
#     M1 = np.exp(ln_M1)
#     M2 = np.exp(ln_M2)
#         N = 100000
# 
# Path: dart_board/posterior.py
# def ln_posterior(x, dart):
# def ln_likelihood(x, dart):
# def posterior_properties(x, output, dart):
#     def get_error_from_kwargs(param, **kwargs):
# def calculate_ADAF_efficiency(Mdot, Mdot_edd, delta=0.1):
# def calculate_L_x(M1, mdot, k1):
# def calc_prob_from_mass_function(M1, M2, f_obs, f_err):
#     def model_h(M1, M2, f):
#     def func_integrand(f, f_obs, f_err, M1, M2):
#     def calc_f(f_obs, f_err, M1, M2):
# def check_output(output, binary_type):
# def get_theta_proj(ra, dec, ra_b, dec_b):
# def get_dtheta_dalpha(alpha, delta, alpha_b, delta_b):
# def get_dtheta_ddelta(alpha, delta, alpha_b, delta_b):
# def get_domega_dalpha(alpha, delta, alpha_b, delta_b):
# def get_domega_ddelta(alpha, delta, alpha_b, delta_b):
# def get_J_coor(alpha, delta, alpha_b, delta_b):
#     M1 = np.exp(ln_M1)
#     M2 = np.exp(ln_M2)
#     M1 = np.exp(ln_M1)
#     M2 = np.exp(ln_M2)
#         N = 100000


import numpy as np
import pickle
import matplotlib
import corner
import matplotlib.pyplot as plt
from matplotlib import font_manager
from dart_board import posterior
from dart_board import posterior
from scipy import stats
          r"$v_{\rm k, i}\ ({\rm km}\ {\rm s}^{-1})$",
          r"$\theta_{\rm k}\ ({\rm rad.})$",
          r"$\phi_{\rm k}\ ({\rm rad.})$",
          r"$t_{\rm i}\ ({\rm Myr})$"]
# plt_range = ([13.9,13.95], [12,13], [0,4000], [0.7,1], [0,750], [0.0,np.pi], [0.0,np.pi], [15,22])
# plt_range = ([0,25], [0,20], [0,4000], [0.0,1], [0,750], [0.0,np.pi], [0.0,np.pi], [10,60])
plt_range = ([0,40], [0,30], [1,4], [0.0,1], [0,750], [0.0,np.pi], [0.0,np.pi], [0,60])


# Load traditional population synthesis results
trad_x_i = np.load("../data/HMXB_trad_chain.npy")
length, ndim = trad_x_i.shape

trad_likelihood = np.load("../data/HMXB_trad_lnprobability.npy")
trad_derived = np.load("../data/HMXB_trad_derived.npy")

length, ndim = trad_x_i.shape
trad_x_i = trad_x_i[int(length*frac):,:]
trad_likelihood = trad_likelihood[int(length*frac):]
trad_derived = trad_derived[int(length*frac):,:]

# trad_x_i = trad_x_i.reshape((len(trad_likelihood), 14))
# trad_derived = trad_derived.reshape((len(trad_likelihood), 9))
print(trad_x_i.shape)
print(trad_derived.shape)



# Make the orbital separation distribution in log-scale
chains.T[2] = np.log10(chains.T[2])

# Path: dart_board/posterior.py
# def ln_posterior(x, dart):
# def ln_likelihood(x, dart):
# def posterior_properties(x, output, dart):
#     def get_error_from_kwargs(param, **kwargs):
# def calculate_ADAF_efficiency(Mdot, Mdot_edd, delta=0.1):
# def calculate_L_x(M1, mdot, k1):
# def calc_prob_from_mass_function(M1, M2, f_obs, f_err):
#     def model_h(M1, M2, f):
#     def func_integrand(f, f_obs, f_err, M1, M2):
#     def calc_f(f_obs, f_err, M1, M2):
# def check_output(output, binary_type):
# def get_theta_proj(ra, dec, ra_b, dec_b):
# def get_dtheta_dalpha(alpha, delta, alpha_b, delta_b):
# def get_dtheta_ddelta(alpha, delta, alpha_b, delta_b):
# def get_domega_dalpha(alpha, delta, alpha_b, delta_b):
# def get_domega_ddelta(alpha, delta, alpha_b, delta_b):
# def get_J_coor(alpha, delta, alpha_b, delta_b):
#     M1 = np.exp(ln_M1)
#     M2 = np.exp(ln_M2)
#     M1 = np.exp(ln_M1)
#     M2 = np.exp(ln_M2)
#         N = 100000
# 
# Path: dart_board/posterior.py
# def ln_posterior(x, dart):
# def ln_likelihood(x, dart):
# def posterior_properties(x, output, dart):
#     def get_error_from_kwargs(param, **kwargs):
# def calculate_ADAF_efficiency(Mdot, Mdot_edd, delta=0.1):
# def calculate_L_x(M1, mdot, k1):
# def calc_prob_from_mass_function(M1, M2, f_obs, f_err):
#     def model_h(M1, M2, f):
#     def func_integrand(f, f_obs, f_err, M1, M2):
#     def calc_f(f_obs, f_err, M1, M2):
# def check_output(output, binary_type):
# def get_theta_proj(ra, dec, ra_b, dec_b):
# def get_dtheta_dalpha(alpha, delta, alpha_b, delta_b):
# def get_dtheta_ddelta(alpha, delta, alpha_b, delta_b):
# def get_domega_dalpha(alpha, delta, alpha_b, delta_b):
# def get_domega_ddelta(alpha, delta, alpha_b, delta_b):
# def get_J_coor(alpha, delta, alpha_b, delta_b):
#     M1 = np.exp(ln_M1)
#     M2 = np.exp(ln_M2)
#     M1 = np.exp(ln_M1)
#     M2 = np.exp(ln_M2)
#         N = 100000


# Path: paper/scripts/HMXB_plots.py
import numpy as np
import pickle
import matplotlib
import corner
import matplotlib.pyplot as plt
from matplotlib import font_manager
from dart_board import posterior
from dart_board import posterior
from scipy import stats

          r"$v_{\rm k, i}\ ({\rm km}\ {\rm s}^{-1})$",
          r"$\theta_{\rm k}\ ({\rm rad.})$",
          r"$\phi_{\rm k}\ ({\rm rad.})$",
          r"$t_{\rm i}\ ({\rm Myr})$"]
# plt_range = ([13.9,13.95], [12,13], [0,4000], [0.7,1], [0,750], [0.0,np.pi], [0.0,np.pi], [15,22])
# plt_range = ([0,25], [0,20], [0,4000], [0.0,1], [0,750], [0.0,np.pi], [0.0,np.pi], [10,60])
plt_range = ([0,40], [0,30], [1,4], [0.0,1], [0,750], [0.0,np.pi], [0.0,np.pi], [0,60])


# Load traditional population synthesis results
trad_x_i = np.load("../data/HMXB_trad_chain.npy")
length, ndim = trad_x_i.shape

trad_likelihood = np.load("../data/HMXB_trad_lnprobability.npy")
trad_derived = np.load("../data/HMXB_trad_derived.npy")

length, ndim = trad_x_i.shape
trad_x_i = trad_x_i[int(length*frac):,:]
trad_likelihood = trad_likelihood[int(length*frac):]
trad_derived = trad_derived[int(length*frac):,:]

# trad_x_i = trad_x_i.reshape((len(trad_likelihood), 14))
# trad_derived = trad_derived.reshape((len(trad_likelihood), 9))
print(trad_x_i.shape)
print(trad_derived.shape)



# Make the orbital separation distribution in log-scale
chains.T[2] = np.log10(chains.T[2])

trad_x_i.T[2] = posterior.P_to_A(trad_x_i.T[0], trad_x_i.T[1], trad_x_i.T[2])