# -*- coding: utf-8 -*-
"""
Created on Wed Dec  6 11:19:43 2023

@author: nielbock
"""
import matplotlib.pyplot as plt
import numpy as np
import os
import matplotlib.ticker as ticker

# Get the user directory
wdir=os.getcwd()

LANG="DE"
#LANG="EN"

EARTHMASS_SI = 5.97217e24
EARTHRADIUS_SI = 6371e3

def sphere_radius(mass, bulk_density):
    mass_E = mass / EARTHMASS_SI
    volume = mass / bulk_density
    radius = (3 * volume / (4 * np.pi)) ** (1 / 3)
    radius_E = radius / EARTHRADIUS_SI
    return mass_E,radius_E

# Mass range for the plot in kg
mass = np.linspace(1e22, 1e27, 5000)  # Adjust the range as needed

# Bulk densities for different functions in kg/m^3
bulk_densities = [996.557, 4107, 5513, 7874]  # Add or modify densities as desired

# Custom labels for substances
if LANG == "DE":
    labels = ['Wasser (konstante Dichte)', 'Silikate (konstante Dichte)', 'Erde (konstante Dichte)', 'Eisen (konstante Dichte)']
else:
    labels = ['Water (incompressible)', 'Silicate (incompressible)', 'Earth (incompressible)', 'Iron (incompressible)']

# Custom colors for curves
colors = ['blue', 'grey', 'green','red']

# Create a figure and axes
fig, ax = plt.subplots()

# Plot the radius for each bulk density
for density, label, color in zip(bulk_densities, labels, colors):
    mass_E,radius_E = sphere_radius(mass, density)
    ax.plot(mass_E, radius_E, label=label, color=color)

# Plot Formula Mass vs. Radius for constant density
formula = r'$M = \frac{4}{3}\,\pi\,R^3\,\bar{\rho}$'
fig.text(0.1, 0.80, formula, transform=ax.transAxes, ha='left', va='bottom',fontsize=1.0 * plt.rcParams['font.size'])

# Set labels and title
if LANG == "DE":
    ax.set_xlabel('Masse (Erdmassen)')
    ax.set_ylabel('Radius (Erdradien)')
    ax.set_title('Kugelradius vs. Masse')
else:
    ax.set_xlabel('Mass (Earth mass)')
    ax.set_ylabel('Radius (Earth radius)')
    ax.set_title('Radius of Sphere vs. Mass')

# Set custom axis limits
x_min, x_max = 0, 10  # Custom x-axis limits
y_min, y_max = 0, 3.5  # Custom y-axis limits
ax.set_xlim(x_min, x_max)
ax.set_ylim(y_min, y_max)
ax.set_xticks(np.linspace(0, 10, 11))
ax.set_yticks(np.arange(0, 3.6, 0.5))

# Change x-axis tick format from point to comma for decimals
if LANG == "DE":
    #formatter = ticker.FuncFormatter(lambda x, pos: f'{x:.1f}'.replace('.', ','))
    #plt.gca().xaxis.set_major_formatter(formatter)
    formatter = ticker.FuncFormatter(lambda y, pos: f'{y:.1f}'.replace('.', ','))
    plt.gca().yaxis.set_major_formatter(formatter)

# Add legend with custom labels
ax.legend()

# Add legend with custom labels
legend_x = 0.98
legend_y = 0.05
ax.legend(loc='lower right', bbox_to_anchor=(legend_x, legend_y), fontsize=0.7 * plt.rcParams['font.size'])

# Export the plot to PNG and PDF with high resolution
dpi = 300  # Resolution in dots per inch

output_path = rf'{wdir}\Plots\Exoplanets_M_vs_R_naive'

if LANG == "DE":
    plt.savefig(f'{output_path}_de.png', dpi=dpi)
    plt.savefig(f'{output_path}_de.pdf', dpi=dpi)
else:
    plt.savefig(f'{output_path}_en.png', dpi=dpi)
    plt.savefig(f'{output_path}_en.pdf', dpi=dpi)

# Display the plot
plt.show()