# 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/forward_pop_synth.py
# def generate_population(dart, N, ra_in=None, dec_in=None):
# def get_v_kick(N, sigma=c.v_kick_sigma):
# def get_theta(N):
# def get_phi(N):
# def get_omega(N):
# def get_M1(N, alpha=c.alpha, M1_min=c.min_mass_M1, M1_max=c.max_mass_M1):
# def get_M2(M1, N, M2_min=c.min_mass_M2, M2_max=c.max_mass_M2):
# def get_a(N, a_min=c.min_a, a_max=c.max_a):
# def get_ecc(N):
# def get_t(N, min_time=c.min_t, max_time=c.max_t):
# def get_z(t_b, N, min_z=c.min_z, max_z=c.max_z):
# def get_new_ra_dec(ra, dec, theta_proj, pos_ang):
#     M1 = dart.generate_M1(N)
#     M2 = dart.generate_M2(M1, N)
#     A = (alpha+1.0) / (np.power(M1_max, alpha+1.0) - np.power(M1_min, alpha+1.0))
#     C = 1.0 / (np.log(a_max) - np.log(a_min))
#     C = 1.0 / (np.log(max_z) - np.log(min_z))


import numpy as np
import sys
import matplotlib
import corner
import matplotlib.pyplot as plt
from matplotlib import font_manager
from dart_board import posterior, forward_pop_synth
from scipy import stats
#           r"$m_f$"]


if sol == 'low':
    plt_range = ([1.2,2.6], [0.0,45.0], [24.0,30.0], [0.0, 0.25])
else:
    plt_range = ([1.2,2.6], [0.0,22.0], [24.0,30.0], [0.0, 0.25])

# Plot observational constraints
obs_Porb, obs_Porb_err = 27.405, 0.5
tmp_x = np.linspace(24.0, 30.0, 1000)
ax[2].plot(tmp_x, stats.norm.pdf(tmp_x, loc=obs_Porb, scale=obs_Porb_err), color='k', zorder=0.1)
ax[3].axvline(0.17, color='k', zorder=0.1)

# trad_Porb = posterior.A_to_P(trad_derived.T[0], trad_derived.T[1], trad_derived.T[2])
# Testing with no weights
# ax[0].hist(trad_derived.T[0], range=plt_range[0], bins=20, normed=True, color='C0', alpha=0.3, label='Traditional')
# ax[1].hist(trad_derived.T[1], range=plt_range[1], bins=20, normed=True, color='C0', alpha=0.3)
# ax[2].hist(trad_Porb, range=plt_range[2], bins=20, normed=True, color='C0', alpha=0.3)
# ax[3].hist(trad_derived.T[3], range=plt_range[3], bins=20, normed=True, color='C0', alpha=0.3)
# ax[0].hist(trad_derived.T[0], range=plt_range[0], bins=20, normed=True, weights=trad_likelihood, color='C0', alpha=0.3, label='Traditional')
# ax[1].hist(trad_derived.T[1], range=plt_range[1], bins=20, normed=True, weights=trad_likelihood, color='C0', alpha=0.3)
# ax[2].hist(trad_Porb, range=plt_range[2], bins=20, normed=True, weights=trad_likelihood, color='C0', alpha=0.3)
# ax[3].hist(trad_derived.T[3], range=plt_range[3], bins=20, normed=True, weights=trad_likelihood, color='C0', alpha=0.3)
# sin_i = np.sin(forward_pop_synth.get_theta(len(trad_Porb)))
# m_f = (trad_derived.T[1] * sin_i)**3 / (trad_derived.T[0]+trad_derived.T[1])**2
# ax[4].hist(m_f, range=plt_range[4], bins=20, normed=True, color='C0', alpha=0.3)


# Plot results from MCMC

# 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/forward_pop_synth.py
# def generate_population(dart, N, ra_in=None, dec_in=None):
# def get_v_kick(N, sigma=c.v_kick_sigma):
# def get_theta(N):
# def get_phi(N):
# def get_omega(N):
# def get_M1(N, alpha=c.alpha, M1_min=c.min_mass_M1, M1_max=c.max_mass_M1):
# def get_M2(M1, N, M2_min=c.min_mass_M2, M2_max=c.max_mass_M2):
# def get_a(N, a_min=c.min_a, a_max=c.max_a):
# def get_ecc(N):
# def get_t(N, min_time=c.min_t, max_time=c.max_t):
# def get_z(t_b, N, min_z=c.min_z, max_z=c.max_z):
# def get_new_ra_dec(ra, dec, theta_proj, pos_ang):
#     M1 = dart.generate_M1(N)
#     M2 = dart.generate_M2(M1, N)
#     A = (alpha+1.0) / (np.power(M1_max, alpha+1.0) - np.power(M1_min, alpha+1.0))
#     C = 1.0 / (np.log(a_max) - np.log(a_min))
#     C = 1.0 / (np.log(max_z) - np.log(min_z))


# Path: paper/scripts/J0513_plots.py
import numpy as np
import sys
import matplotlib
import corner
import matplotlib.pyplot as plt
from matplotlib import font_manager
from dart_board import posterior, forward_pop_synth
from scipy import stats

#           r"$m_f$"]


if sol == 'low':
    plt_range = ([1.2,2.6], [0.0,45.0], [24.0,30.0], [0.0, 0.25])
else:
    plt_range = ([1.2,2.6], [0.0,22.0], [24.0,30.0], [0.0, 0.25])

# Plot observational constraints
obs_Porb, obs_Porb_err = 27.405, 0.5
tmp_x = np.linspace(24.0, 30.0, 1000)
ax[2].plot(tmp_x, stats.norm.pdf(tmp_x, loc=obs_Porb, scale=obs_Porb_err), color='k', zorder=0.1)
ax[3].axvline(0.17, color='k', zorder=0.1)

# trad_Porb = posterior.A_to_P(trad_derived.T[0], trad_derived.T[1], trad_derived.T[2])
# Testing with no weights
# ax[0].hist(trad_derived.T[0], range=plt_range[0], bins=20, normed=True, color='C0', alpha=0.3, label='Traditional')
# ax[1].hist(trad_derived.T[1], range=plt_range[1], bins=20, normed=True, color='C0', alpha=0.3)
# ax[2].hist(trad_Porb, range=plt_range[2], bins=20, normed=True, color='C0', alpha=0.3)
# ax[3].hist(trad_derived.T[3], range=plt_range[3], bins=20, normed=True, color='C0', alpha=0.3)
# ax[0].hist(trad_derived.T[0], range=plt_range[0], bins=20, normed=True, weights=trad_likelihood, color='C0', alpha=0.3, label='Traditional')
# ax[1].hist(trad_derived.T[1], range=plt_range[1], bins=20, normed=True, weights=trad_likelihood, color='C0', alpha=0.3)
# ax[2].hist(trad_Porb, range=plt_range[2], bins=20, normed=True, weights=trad_likelihood, color='C0', alpha=0.3)
# ax[3].hist(trad_derived.T[3], range=plt_range[3], bins=20, normed=True, weights=trad_likelihood, color='C0', alpha=0.3)
# sin_i = np.sin(forward_pop_synth.get_theta(len(trad_Porb)))
# m_f = (trad_derived.T[1] * sin_i)**3 / (trad_derived.T[0]+trad_derived.T[1])**2
# ax[4].hist(m_f, range=plt_range[4], bins=20, normed=True, color='C0', alpha=0.3)


# Plot results from MCMC

MCMC_Porb = posterior.A_to_P(MCMC_derived.T[idx_M1], MCMC_derived.T[idx_M2], MCMC_derived.T[idx_a])