mr-plotter.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b468e3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "############################\n",
    "#@|++++++mr-plotter++++++#@|\n",
    "############################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a44b05f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pyvo\n",
    "import json\n",
    "import argparse\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import interpolate\n",
    "from astropy.constants import G\n",
    "import matplotlib.pyplot as plt\n",
    "import astropy.constants as const\n",
    "from astropy.io.votable import parse\n",
    "from configparser import ConfigParser"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "788c3728",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|+++++++++++++++++We load the dictionaries+++++++++++++++++++++++++\n",
    "zeng_models_colors = json.load(open('misc/dicts/zeng_models_colors.txt'))\n",
    "zeng_models_labels = json.load(open('misc/dicts/zeng_models_labels.txt'))\n",
    "color_coding_labels = json.load(open('misc/dicts/color_coding_labels.txt'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b218a084",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|++++Uncomment this for the mr-plotter.py version+++++++\n",
    "#parser = argparse.ArgumentParser()\n",
    "#parser.add_argument('config_file')\n",
    "#args = parser.parse_args()\n",
    "#config_file = args.config_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7ac21f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|++++Uncomment this for the mr.plotter.ipynb version+++++\n",
    "config_file = 'example5.ini'\n",
    "#@|++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de483a68",
   "metadata": {},
   "outputs": [],
   "source": [
    "path_models = 'theoretical_models/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "baac94f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "#|Turbet et al. (2020) theoretical M-R models\n",
    "def turbet2020_MR():\n",
    "    global R_pl_turb, M_pl_turb, m_pure_iron, m_earth_like, m_pure_rock, r_pure_iron, r_earth_like, r_pure_rock \n",
    "\n",
    "    #@|M-R relationships of the core of the planet based on Zeng et al. (2019)\n",
    "    #@|Earth-like, pure iron, and pure rock \n",
    "    pure_iron_df = pd.read_csv('theoretical_models/zeng_2019_pure_iron', sep = \"\\t\", header = None)\n",
    "    earth_like_df = pd.read_csv('theoretical_models/zeng_2019_earth_like', sep = \"\\t\", header = None)\n",
    "    pure_rock_df = pd.read_csv('theoretical_models/zeng_2019_pure_rock', sep = \"\\t\", header = None)\n",
    "\n",
    "    m_pure_iron, r_pure_iron = pure_iron_df[0].values,  pure_iron_df[1].values\n",
    "    m_earth_like, r_earth_like = earth_like_df[0].values,  earth_like_df[1].values\n",
    "    m_pure_rock, r_pure_rock = pure_rock_df[0].values,  pure_rock_df[1].values\n",
    "\n",
    "    Cores_db = np.array(['earth', 'rock', 'iron'])                          #@|Cores database\n",
    "    R_db = np.array([r_earth_like, r_pure_rock, r_pure_iron], dtype=object) #@|Corresponding radius database\n",
    "    M_db = np.array([m_earth_like, m_pure_rock, m_pure_iron], dtype=object) #@|Corresponding masses database\n",
    "\n",
    "    R_Core, M_Core = [], []\n",
    "    for core in Core_turb2020:\n",
    "        idx = np.where(core==Cores_db)[0]\n",
    "        R_Core.extend(R_db[idx])\n",
    "        M_Core.extend(M_db[idx])\n",
    "\n",
    "\n",
    "    R_pl_turb, M_pl_turb = [], []\n",
    "\n",
    "    for i,frac in enumerate(WMFs_turb2020):\n",
    "\n",
    "        r_core = R_Core[i]\n",
    "        m_core = M_Core[i]\n",
    "\n",
    "        R_cons,M_h20, Pt = 8.314, 1.8e-2, 0.1 # fixed\n",
    "        g_core = G.value * m_core * const.M_earth.value / (r_core*const.R_earth.value)**2\n",
    "        g = G.value * M_turb2020 * const.M_earth.value / (R_turb2020*const.R_earth.value)**2\n",
    "\n",
    "        alpha_1, alpha_2, alpha_3, alpha_4, alpha_5, alpha_6 = -3.550, 1.310, 1.099, 4.683e-1, 7.664e-1, 4.224e-1\n",
    "        beta_1, beta_2, beta_3, beta_4, beta_5 = 2.846, 1.555e-1, 8.777e-2, 6.045e-2, 1.143e-2\n",
    "        beta_6, beta_7, beta_8, beta_9, beta_10 = 1.736e-2, 1.859e-2, 4.314e-2, 3.393e-2, -1.034e-2\n",
    "\n",
    "        X = (np.log10(frac)-alpha_1) / alpha_2\n",
    "        Y = (np.log10(g)-alpha_3) / alpha_4\n",
    "        Z = (np.log10(Seff_turb2020)-alpha_5)/alpha_6\n",
    "\n",
    "        X = (np.log10(frac)-alpha_1) / alpha_2\n",
    "        Y = (np.log10(g)-alpha_3) / alpha_4\n",
    "        Z = (np.log10(Seff_turb2020)-alpha_5)/alpha_6\n",
    "\n",
    "        Teff = 10**(beta_1+beta_2*X+beta_3*Y+beta_4*Z+beta_5*X*Y+beta_6*Y**2+beta_7*X**3+beta_8*X**2*Y+beta_9*X*Y**2+beta_10*Y**4)\n",
    "\n",
    "        #@|z_atm \n",
    "        a = np.log(frac / (1-frac) * g_core**2 / 4 / np.pi / G.value / Pt) \n",
    "        b = R_cons * Teff / M_h20 / g_core\n",
    "        c = 1 / (r_core*const.R_earth.value)\n",
    "        z_atm = (1 / a / b - c)**(-1)\n",
    "\n",
    "        #r_new = (r_earth_like.values*const.R_earth.value + z_atm) / const.R_earth.value \n",
    "        #m_new = m_earth_like.values / (1-frac)\n",
    "\n",
    "        r_pl = (r_core*const.R_earth.value + z_atm) / const.R_earth.value \n",
    "        m_pl = m_core / (1-frac)\n",
    "\n",
    "        r_pl = np.array(r_pl)\n",
    "        \n",
    "        #@|####-Results-########\n",
    "        R_pl_turb.append(r_pl)\n",
    "        M_pl_turb.append(m_pl)\n",
    "        #@|####################\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e21b32da",
   "metadata": {},
   "outputs": [],
   "source": [
    "def aguich2021_MR():\n",
    "    \n",
    "    global R_aguich2021, M_aguich2021\n",
    "\n",
    "    R_aguich2021, M_aguich2021 = [], []\n",
    "    for i in range(len(x_core_aguich2021)):\n",
    "\n",
    "        idxs_aguich2021 = np.where((df_aguichine_2021['x_core'] == x_core_aguich2021[i]) & \\\n",
    "                                   (df_aguichine_2021['x_H2O'] == x_H2O_aguich2021[i]) & \\\n",
    "                                   (df_aguichine_2021['T_irr'] == Tirr_aguich2021[i]) & \\\n",
    "                                   (df_aguichine_2021['errcode'] == 0))[0] #@|only pick the errcode = 0 values\n",
    "\n",
    "        r_aguich2021 = (df_aguichine_2021['R_b'].values + df_aguichine_2021['R_a'].values)[idxs_aguich2021]\n",
    "        m_aguich2021 = (df_aguichine_2021['M_b'].values + df_aguichine_2021['M_a'].values)[idxs_aguich2021]\n",
    "\n",
    "\n",
    "        R_aguich2021.append(r_aguich2021)\n",
    "        M_aguich2021.append(m_aguich2021)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53b4952c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def lopez_fortney2014_MR():\n",
    "    \n",
    "    global R_lf2014, M_lf2014\n",
    "    \n",
    "    R_lf2014, M_lf2014 = [], []\n",
    "    \n",
    "    for i in range(len(age_lf2014)):\n",
    "        \n",
    "        df_lf2014 = pd.read_csv(f'{path_models}/Lopez&Fortney_2014_{age_lf2014[i]}_{opacity_lf2014[i]}', sep = \" \")\n",
    "        \n",
    "        idxs_flux_lf2014 = np.where(df_lf2014['Flux'].values == float(Seff_lf2014[1]))[0]\n",
    "        \n",
    "        \n",
    "        m_lf2014 = df_lf2014['Mass'][idxs_flux_lf2014].values\n",
    "        r_lf2014 = df_lf2014[f'{H_He[i]}'][idxs_flux_lf2014].values\n",
    "        \n",
    "        R_lf2014.append(r_lf2014)\n",
    "        M_lf2014.append(m_lf2014)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8e28b07",
   "metadata": {},
   "outputs": [],
   "source": [
    "def haldemann2024_MR():\n",
    "\n",
    "    global R_haldemann2024, M_haldemann2024\n",
    "    \n",
    "    R_haldemann2024, M_haldemann2024 = [], []\n",
    "\n",
    "    for i in range(len(models_haldemann2024)):\n",
    "        \n",
    "        df_haldemann2024 = pd.read_csv(f'{path_models}/Mass_Radius_{models_haldemann2024[i]}.dat', \\\n",
    "                               comment = '#', header = None)\n",
    "\n",
    "        if T_out_haldemann2024[i] == 50:\n",
    "            idx_R_haldemann2024 = 1\n",
    "        if T_out_haldemann2024[i] == 300:\n",
    "            idx_R_haldemann2024 = 2\n",
    "        if T_out_haldemann2024[i] == 800:\n",
    "            idx_R_haldemann2024 = 3\n",
    "        if T_out_haldemann2024[i] == 1500:\n",
    "            idx_R_haldemann2024 = 4\n",
    "        if T_out_haldemann2024[i] == 2000:\n",
    "            idx_R_haldemann2024 = 5\n",
    "\n",
    "\n",
    "        m_haldemann2024 = df_haldemann2024[0].values\n",
    "        r_haldemann2024 = df_haldemann2024[idx_R_haldemann2024].values\n",
    "        \n",
    "        R_haldemann2024.append(r_haldemann2024)\n",
    "        M_haldemann2024.append(m_haldemann2024)\n",
    "          "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74737fc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def seager2007_MR():\n",
    "    \n",
    "    global R_seager2007, M_seager2007\n",
    "    \n",
    "    R_seager2007, M_seager2007 = [], []\n",
    "    \n",
    "    m_seager2007 = np.linspace(0, 20, 1000)\n",
    "    \n",
    "    for i in range(len(models_seager2007)):\n",
    "    \n",
    "        if models_seager2007[i] == 'iron':\n",
    "                   \n",
    "            k1, k2, k3 = -0.209490, 0.0804, 0.394\n",
    "            m1, r1 = 5.80, 2.52\n",
    "        \n",
    "        if models_seager2007[i] == 'rock':\n",
    "                   \n",
    "            k1, k2, k3 = -0.209594, 0.0799, 0.413\n",
    "            m1, r1 = 10.55, 3.90\n",
    "\n",
    "        if models_seager2007[i] == 'water':\n",
    "                   \n",
    "            k1, k2, k3 = -0.209396, 0.0807, 0.375\n",
    "            m1, r1 = 5.52, 4.43\n",
    "                    \n",
    "        log_r_seager2007 = k1 + 1/3 * np.log10(m_seager2007/m1) - k2 * (m_seager2007/m1)**k3\n",
    "        r_seager2007 = 10**log_r_seager2007\n",
    "        r_seager2007 = r_seager2007 * r1\n",
    "\n",
    "        \n",
    "        M_seager2007.append(m_seager2007)\n",
    "        R_seager2007.append(r_seager2007)\n",
    "                    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd0181ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "def otegi2020_MR():\n",
    "    \n",
    "    global R_otegi_small, M_otegi_small, R_otegi_intermediate, M_otegi_intermediate\n",
    "    \n",
    "    M_otegi_small = np.linspace(0, 10, 100)\n",
    "    M_otegi_intermediate = np.linspace(6, 138, 100)\n",
    "    \n",
    "    R_otegi_small = 1.03 * M_otegi_small**0.29\n",
    "    R_otegi_intermediate = 0.70 * M_otegi_intermediate**0.63"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5aa87d39",
   "metadata": {},
   "outputs": [],
   "source": [
    "def parc2024_MR():\n",
    "    \n",
    "    global R_parc_small, M_parc_small, R_parc_intermediate,\\\n",
    "    M_parc_intermediate, R_parc_giant, M_parc_giant\n",
    "    \n",
    "    M_parc_small = np.linspace(0, 10, 100)\n",
    "    M_parc_intermediate = np.linspace(6, 138, 100)\n",
    "    M_parc_giant = np.linspace(138, 10000, 100)\n",
    "    \n",
    "    \n",
    "    R_parc_small = 1.02 * M_parc_small**0.28\n",
    "    R_parc_intermediate = 0.61 * M_parc_intermediate**0.67\n",
    "    R_parc_giant = 11.9 * M_parc_giant**0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbcee882",
   "metadata": {},
   "outputs": [],
   "source": [
    "def iso_density():\n",
    "    \n",
    "    global Radii_iso, Masses_iso\n",
    "    \n",
    "    #@|idodensity curves\n",
    "    radii = np.linspace(0, const.R_earth.value * 25, 10000)\n",
    "    Radii_iso, Masses_iso = [], []\n",
    "    for i in range(len(density)):\n",
    "        \n",
    "        masses = density[i] * 4 / 3 * np.pi *  radii**3\n",
    "        \n",
    "        Radii_iso.append(radii / const.R_earth.value)\n",
    "        Masses_iso.append(masses / const.M_earth.value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7b0aa69",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c91d45f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|Read the configuration file 'config.ini' into a config_object\n",
    "config_object = ConfigParser()\n",
    "config_object.read(f\"config/{config_file}\")\n",
    "\n",
    "#@|Sections\n",
    "try:\n",
    "    CATALOG_DATA = config_object['CATALOG_DATA']\n",
    "except:\n",
    "    pass\n",
    "try:\n",
    "    MY_DATA = config_object['MY_DATA']\n",
    "except:\n",
    "    pass\n",
    "try:\n",
    "    MODELS = config_object['MODELS']\n",
    "except:\n",
    "    pass\n",
    "try:\n",
    "    OPTIONAL_CONFIG = config_object['OPTIONAL_CONFIG']\n",
    "except:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8133f8a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4049843f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|++++++Exoplanet catalog++++++++++\n",
    "#@|NEA, Exoplanet.eu, or PlanetS\n",
    "#@|++++++++++++++++++++++++++++++++++\n",
    "catalog = CATALOG_DATA['catalog']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dc3813b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
    "#@|Load data (df) from the NASA Exoplanet Archive (Confirmed planets)\n",
    "#@|+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
    "\n",
    "if catalog == 'NEA':\n",
    "\n",
    "    try:\n",
    "        web_or_local = CATALOG_DATA['web_or_local']\n",
    "    except: \n",
    "        web_or_local = 'local'\n",
    "    \n",
    "    try:\n",
    "        ps_or_composite = CATALOG_DATA['ps_or_composite']\n",
    "    except:\n",
    "        ps_or_composite = 'ps'\n",
    "\n",
    "    if web_or_local == 'web':\n",
    "\n",
    "        if ps_or_composite == 'ps':\n",
    "            print('Downloading the \"Planetary Systems\" table from the NASA Exoplanet Archive ... (this might take a while)')\n",
    "\n",
    "        if ps_or_composite == 'composite':\n",
    "            print('Downloading the \"Planetary Systems Composite Data\" table from the NASA Exoplanet Archive ... (this might take a while)')   \n",
    "\n",
    "        request = pyvo.dal.TAPQuery('https://exoplanetarchive.ipac.caltech.edu/TAP', 'SELECT * FROM '+ps_or_composite)\n",
    "        table = pyvo.dal.TAPQuery.execute(request)\n",
    "\n",
    "        df = pd.DataFrame(table)\n",
    "\n",
    "    if web_or_local == 'local':\n",
    "\n",
    "        list_complete = sorted(os.listdir('catalog_data/NEA/'))\n",
    "        #list_complete\n",
    "\n",
    "        list_ps = []\n",
    "        list_composite = []\n",
    "        for item in list_complete:\n",
    "            if item[0:10] == 'PSCompPars':\n",
    "                list_composite.append(item)\n",
    "\n",
    "            else:\n",
    "                list_ps.append(item)\n",
    "\n",
    "        if ps_or_composite == 'ps':\n",
    "\n",
    "            df = pd.read_csv('catalog_data/NEA/'+list_ps[-1], comment = \"#\")\n",
    "\n",
    "        if ps_or_composite == 'composite':\n",
    "\n",
    "            df = pd.read_csv('catalog_data/NEA/'+list_composite[-1], comment = \"#\")\n",
    "            \n",
    "            \n",
    "#@|+++++++++++++++++++++++++++++++++++++++++++\n",
    "#@|Load data df from the Exoplanet.eu catalog \n",
    "#@|+++++++++++++++++++++++++++++++++++++++++++ \n",
    "    \n",
    "    \n",
    "if catalog == 'Exoplanet.eu':\n",
    "    \n",
    "    file = os.listdir('catalog_data/Exoplanet.eu/')[0]\n",
    "    df = pd.read_csv('catalog_data/Exoplanet.eu/'+file)\n",
    "    \n",
    "    df = df.rename(columns={'orbital_period':'pl_orbper', 'semi_major_axis':'pl_orbsmax',\\\n",
    "                           'eccentricity':'pl_orbeccen','discovered':'disc_year', \\\n",
    "                            'temp_calculated':'pl_eqt_Exoplanet.eu', 'temp_measured':'pl_eqt',\\\n",
    "                            'log_g':'st_logg','mag_v':'V Mag', 'mag_k':'K Mag','star_distance': 'sy_dist',\\\n",
    "                            'star_metallicity':'st_met', 'star_mass':'st_mass', 'star_radius': 'st_rad',\\\n",
    "                            'star_age': 'st_age', 'star_teff':'st_teff',\\\n",
    "                            'mass':'pl_bmasse','mass_error_min':'pl_bmasseerr2',\\\n",
    "                            'mass_error_max':'pl_bmasseerr1','radius':'pl_rade',\\\n",
    "                            'radius_error_min':'pl_radeerr2','radius_error_max':'pl_radeerr1' })            \n",
    "            \n",
    "            \n",
    "#@|++++++++++++++++++++++++++++++++++++++\n",
    "#@|Load data df from the PlanetS catalog \n",
    "#@|++++++++++++++++++++++++++++++++++++++\n",
    "\n",
    "if catalog == 'PlanetS':\n",
    " \n",
    "    votable_PlanetS = parse(\"catalog_data/PlanetS/PLANETS.vot\")\n",
    "    table_PlanetS = votable_PlanetS.get_first_table().to_table()\n",
    "    df = table_PlanetS.to_pandas()\n",
    "    \n",
    "    df = df.rename(columns={'Discovery Year':'disc_year', 'Discovery Method':'discoverymethod', \\\n",
    "                        'Discovery Facility':'disc_facility','Orbital Period [days]':'pl_orbper',\\\n",
    "                        'Orbit Semi-Major axis [au]':'pl_orbsmax','Eccentricity':'pl_orbeccen',\\\n",
    "                        'Insolation Flux [Earth Flux]':'pl_insol',\\\n",
    "                        'Insolation Flux [Earth Flux] - Computation':'pl_insol_PlanetS',\\\n",
    "                        'Equilibrium Temperature [K]':'pl_eqt',\\\n",
    "                        'Equilibrium Temperature [K] - Computation':'pl_eqt_PlanetS',\\\n",
    "                        'Stellar Effective Temperature [K]':'st_teff','Teff [K] (Gaia DR3)':'st_teff_gaia',\\\n",
    "                        'Stellar Radius [Rsun]':'st_rad','Stellar Mass [Msun]':'st_mass',\\\n",
    "                        'Stellar Metallicity [dex]':'st_met','[Fe/H] [dex] (Gaia DR3)':'st_met_gaia',\\\n",
    "                        'Stellar Surface Gravity [log(cm/s**2)]':'st_logg',\\\n",
    "                        'log(g) [log(cm/s**2)] (Gaia DR3)':'st_logg_gaia','Distance [pc]':'sy_dist',\\\n",
    "                        'V Mag':'sy_vmag','K Mag':'sy_kmag','Transmission Spectroscopy Metric (TSM) - Computation':'TSM',\\\n",
    "                        'Emission Spectroscopy Metric (ESM) - Computation':'ESM',\\\n",
    "                        'Measured Radial Velocity Semi-Amplitude [m/s]':'K_measured','Stellar Age [Gyr]':'st_age',\\\n",
    "                        'Stellar Age [Gyr] (Gaia DR3)':'st_age_gaia','Stellar Luminosity [log(Lsun)]':'st_lum',\\\n",
    "                        'Stellar Luminosity [log(Lsun)] (Gaia DR3)':'st_lum_gaia',\\\n",
    "                        'Planet Mass [Mjup]':'pl_bmasse','Planet Mass - Upper Unc [Mjup]':'pl_bmasseerr1',\\\n",
    "                        'Planet Mass - Lower Unc [Mjup]':'pl_bmasseerr2',\\\n",
    "                        'Planet Radius [Rjup]':'pl_rade','Planet Radius - Upper Unc [Rjup]':'pl_radeerr1',\\\n",
    "                        'Planet Radius - Lower Unc [Rjup]':'pl_radeerr2' })\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ef83a9a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49895a46",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
    "#@|We remove those planets with a precision in mass and radius lower than precision_radius and precision_mass\n",
    "#@|++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
    "\n",
    "precision_radius = float(CATALOG_DATA['precision_radius'])\n",
    "precision_mass = float(CATALOG_DATA['precision_mass'])\n",
    "\n",
    "\n",
    "#@|mass precision better than thres\n",
    "M_pl, M_pl_err = df['pl_bmasse'].values, \\\n",
    "(abs(df['pl_bmasseerr1'].values) + \\\n",
    " abs(df['pl_bmasseerr2'].values)) / 2\n",
    "idxs_M_pl_thres = np.where(M_pl_err / M_pl < (precision_mass / 100))[0]\n",
    "#print(str(len(idxs_M_pl_thres))+' of '+str(n_planets)+' planets ('+str(int(len(idxs_M_pl_thres) / n_planets * 100))+'%) have a mass measurement precision below '+str(int(thres*100))+'%.')\n",
    "df = df.iloc[idxs_M_pl_thres]\n",
    "\n",
    "#@|radius precision better than thres\n",
    "R_pl, R_pl_err = df['pl_rade'].values,\\\n",
    "(abs(df['pl_radeerr1'].values) +\\\n",
    " abs(df['pl_radeerr2'].values)) / 2\n",
    "idxs_R_pl_thres = np.where(R_pl_err / R_pl < (precision_radius / 100))[0]\n",
    "#print(str(len(idxs_R_pl_thres))+' of '+str(n_planets)+' planets ('+str(int(len(idxs_R_pl_thres) / n_planets * 100))+'%) have a mass and radius precision below '+str(int(thres*100))+'%.')\n",
    "df = df.iloc[idxs_R_pl_thres]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c38131a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c86404fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "color_coding = CATALOG_DATA['color_coding']\n",
    "\n",
    "#@|+++++Numerical and non-numerical color coding++++++++++\n",
    "if color_coding == 'disc_facility' or color_coding == 'disc_year' or color_coding == 'discoverymethod':\n",
    "    numerical_color_coding = False\n",
    "else:\n",
    "    numerical_color_coding = True\n",
    "\n",
    "\n",
    "try:\n",
    "    groups = [str(x.strip()) for x in CATALOG_DATA['groups'].split(',')]\n",
    "    \n",
    "    if color_coding == 'disc_year':\n",
    "        groups = np.array(groups).astype(int)\n",
    "    \n",
    "except:\n",
    "    pass\n",
    "\n",
    "try:\n",
    "    colors_groups = [str(x.strip()) for x in CATALOG_DATA['colors_groups'].split(',')]\n",
    "except:\n",
    "    pass\n",
    "\n",
    "try:\n",
    "    plot_all_planets = CATALOG_DATA['plot_all_planets'] == 'True'\n",
    "except:\n",
    "    plot_all_planets = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2ad8a05",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38ad8bac",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
    "#@|We remove the rows from the df that have no 'color_coding'\n",
    "#@|----------Only in the numerical color coding-------------\n",
    "#@|++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
    "\n",
    "cols = df.columns.values\n",
    "\n",
    "if numerical_color_coding:\n",
    "    \n",
    "\n",
    "    if color_coding in cols:\n",
    "        idxs_not_nan = np.where(~np.isnan(df[color_coding]))[0]\n",
    "        df = df.iloc[idxs_not_nan]\n",
    "\n",
    "    print('Your sample contains '+str(len(df))+' planets')\n",
    "\n",
    "        \n",
    "else:\n",
    "    \n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a765d31",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "089f5830",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|+++++++++++++++++++to be done++++++++++++++++++++++++++\n",
    "#@|idxs_without_ttvs = np.where(df['ttv_flag'] == 0)\n",
    "#@|df = df.iloc[idxs_without_ttvs]\n",
    "#@print('Your sample contains '+str(len(df))+' planets')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "505d327f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ddc76e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|+++++++++++++++++++++++++++++++++++++++++++++++++++\n",
    "#@|----Definition of the mass and radius arrays-----\n",
    "#@|+++++++++++++++++++++++++++++++++++++++++++++++++++\n",
    "\n",
    "#@|Note: The PlanetS catalog masses and radii are in Rjup. We convert them into Rearth. \n",
    "\n",
    "if catalog == 'PlanetS' or catalog == 'Exoplanet.eu':\n",
    "    \n",
    "    fac_R = const.R_jup.value / const.R_earth.value\n",
    "    fac_M = const.M_jup.value / const.M_earth.value\n",
    "    \n",
    "else:\n",
    "    \n",
    "    fac_M, fac_R = 1, 1\n",
    "    \n",
    "\n",
    "R_catalog, R_catalog_err_up, R_catalog_err_down = df['pl_rade'].values, df['pl_radeerr1'].values, df['pl_radeerr2'].values\n",
    "M_catalog, M_catalog_err_up, M_catalog_err_down = df['pl_bmasse'].values, df['pl_bmasseerr1'].values, df['pl_bmasseerr2'].values\n",
    "\n",
    "\n",
    "R_catalog, R_catalog_err_up, R_catalog_err_down = R_catalog*fac_R, R_catalog_err_up*fac_R, R_catalog_err_down*fac_R\n",
    "M_catalog, M_catalog_err_up, M_catalog_err_down = M_catalog*fac_M, M_catalog_err_up*fac_M, M_catalog_err_down*fac_M\n",
    "\n",
    "#@|Mass and radius upper and lower errors in the same array (for the plots)\n",
    "\n",
    "R_catalog_ERR = [(abs(R_catalog_err_down[i]), R_catalog_err_up[i]) for i in range(len(R_catalog))]\n",
    "M_catalog_ERR = [(abs(M_catalog_err_down[i]), M_catalog_err_up[i]) for i in range(len(M_catalog))]\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57f0e797",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c4e2c0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|#######--zeng et al. (2016,2019) config--##########\n",
    "try:\n",
    "    models_zeng = MODELS['models_zeng']\n",
    "    models_zeng = [x.strip() for x in models_zeng.split(',')]\n",
    "    zeng = True\n",
    "except:\n",
    "    zeng = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b926c682",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|#######--marcus et al. 2010 config--##########\n",
    "try:\n",
    "    models_marcus = MODELS['models_marcus']\n",
    "    models_marcus = [x.strip() for x in models_marcus.split(',')]\n",
    "    marcus = True\n",
    "except:\n",
    "    marcus = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af2a44a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|#######--turbet et al. 2020 config--##########\n",
    "try:\n",
    "    M_turb2020 = float(MODELS['M_turb2020'])\n",
    "    R_turb2020 = float(MODELS['R_turb2020'])\n",
    "    Seff_turb2020 =  float(MODELS['Seff_turb2020'])\n",
    "    WMFs_turb2020 = [float(x.strip()) for x in MODELS['WMFs_turb2020'].split(',')]\n",
    "    Core_turb2020 = [x.strip() for x in MODELS['Core_turb2020'].split(',')]\n",
    "    colors_turb2020 = [x.strip() for x in MODELS['colors_turb2020'].split(',')]\n",
    "    turb_2020 = True\n",
    "    turbet2020_MR()\n",
    "    \n",
    "except:\n",
    "    turb_2020 = False    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f040e0ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|#######--aguichine et al. 2020 config--##########\n",
    "try:\n",
    "    x_core_aguich2021 = [float(x.strip()) for x in MODELS['x_core_aguich2021'].split(',')]\n",
    "    x_H2O_aguich2021 = [float(x.strip()) for x in MODELS['x_H2O_aguich2021'].split(',')]\n",
    "    Tirr_aguich2021 = [float(x.strip()) for x in MODELS['Tirr_aguich2021'].split(',')]\n",
    "    colors_aguich2021 = [x.strip() for x in MODELS['colors_aguich2021'].split(',')]\n",
    "    df_aguichine_2021 = pd.read_csv('theoretical_models/aguichine_2021', comment = \"#\", sep = '\\t')\n",
    "    aguich2021 = True\n",
    "    aguich2021_MR()\n",
    "except:\n",
    "    aguich2021 = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "735f841e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|#######--lopez & fortney et al. 2014 config--##########\n",
    "try:\n",
    "    age_lf2014 = [x.strip() for x in MODELS['age_lf2014'].split(',')]\n",
    "    opacity_lf2014 = [x.strip() for x in MODELS['opacity_lf2014'].split(',')]\n",
    "    Seff_lf2014 = [x.strip() for x in MODELS['Seff_lf2014'].split(',')]\n",
    "    H_He = [x.strip() for x in MODELS['H_He'].split(',')]\n",
    "    colors_lf2014 = [x.strip() for x in MODELS['colors_lf2014'].split(',')]\n",
    "    lopez_fortney2014 = True\n",
    "    lopez_fortney2014_MR()\n",
    "except:\n",
    "    lopez_fortney2014 = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3726543",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|#######----------haldemann et al. (2024) config---------##########\n",
    "try:\n",
    "    models_haldemann2024 = [x.strip() for x in MODELS['models_haldemann2024'].split(',')]\n",
    "    T_out_haldemann2024 =  [float(x.strip()) for x in MODELS['T_out_haldemann2024'].split(',')]\n",
    "    colors_haldemann2024 = [x.strip() for x in MODELS['colors_haldemann2024'].split(',')]\n",
    "\n",
    "    haldemann2024 = True\n",
    "    haldemann2024_MR()\n",
    "except:\n",
    "    haldemann2024 = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85a9f1c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|#######----------seager et al. (2007) config---------##########\n",
    "try:\n",
    "    models_seager2007 = [x.strip() for x in MODELS['models_seager2007'].split(',')]\n",
    "    colors_seager2007 = [x.strip() for x in MODELS['colors_seager2007'].split(',')]\n",
    "\n",
    "    seager2007 = True\n",
    "    seager2007_MR()\n",
    "except:\n",
    "    seager2007 = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08f83f5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|#######------otegi et al. 2020 config------##########\n",
    "try:\n",
    "    models_otegi2020 = MODELS['relations_otegi2020']\n",
    "    models_otegi2020 = [x.strip() for x in models_otegi2020.split(',')]\n",
    "    otegi2020 = True\n",
    "    \n",
    "    try:\n",
    "        \n",
    "        colors_otegi2020 = MODELS['colors_otegi2020']\n",
    "        colors_otegi2020 = [x.strip() for x in colors_otegi2020.split(',')]\n",
    "        \n",
    "    except:\n",
    "        \n",
    "        colors_otegi2020[0] = 'k'\n",
    "        colors_otegi2020[1] = 'k'\n",
    "        \n",
    "    try:\n",
    "        \n",
    "        linestyles_otegi2020 = MODELS['linestyles_otegi2020']\n",
    "        linestyles_otegi2020 = [x.strip() for x in linestyles_otegi2020.split(',')]\n",
    "        \n",
    "    except:\n",
    "        \n",
    "        linestyles_otegi2020[0] = 'dashed'\n",
    "        linestyles_otegi2020[1] = 'dashed'\n",
    "         \n",
    "except:\n",
    "    \n",
    "    otegi2020 = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80de684a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|#######------parc et al. 2024 config------##########\n",
    "try:\n",
    "    models_parc2024 = MODELS['relations_parc2024']\n",
    "    models_parc2024 = [x.strip() for x in models_parc2024.split(',')]\n",
    "    parc2024 = True\n",
    "    \n",
    "    try:\n",
    "        \n",
    "        colors_parc2024 = MODELS['colors_parc2024']\n",
    "        colors_parc2024 = [x.strip() for x in colors_parc2024.split(',')]\n",
    "        \n",
    "    except:\n",
    "        \n",
    "        colors_parc2024[0] = 'k'\n",
    "        colors_parc2024[1] = 'k'\n",
    "        colors_parc2024[2] = 'k'\n",
    "        \n",
    "    try:\n",
    "        \n",
    "        linestyles_parc2024 = MODELS['linestyles_parc2024']\n",
    "        linestyles_parc2024 = [x.strip() for x in linestyles_parc2024.split(',')]\n",
    "        \n",
    "    except:\n",
    "        \n",
    "        linestyles_parc2024[0] = 'dashed'\n",
    "        linestyles_parc2024[1] = 'dashed'\n",
    "        linestyles_parc2024[2] = 'dashed'\n",
    "         \n",
    "except:\n",
    "    \n",
    "    parc2024 = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59977ca3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|#######--Isodensity config--##########\n",
    "try:\n",
    "    density = [float(x.strip()) for x in MODELS['density'].split(',')]\n",
    "    for i in range(len(density)):\n",
    "        density[i] = density[i] * 1000\n",
    "    colors_density = [x.strip() for x in MODELS['colors_density'].split(',')]\n",
    "    isodensity = True\n",
    "except:\n",
    "    isodensity = False "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56583c95",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97585460",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@|optional configuration\n",
    "\n",
    "#@|number of columns to make the plot (one or two)\n",
    "\n",
    "try:\n",
    "    n_cols = OPTIONAL_CONFIG['n_cols']\n",
    "except:\n",
    "    n_cols = 'one'\n",
    "\n",
    "#@|#####--color_max and color_min--###########\n",
    "\n",
    "if color_coding != 'none':\n",
    "    \n",
    "    if numerical_color_coding:\n",
    "        \n",
    "        try: \n",
    "            color_min = float(OPTIONAL_CONFIG['color_min'])\n",
    "\n",
    "        except:\n",
    "            \n",
    "            five_per_cent = int(len(df[color_coding].values)*0.05)\n",
    "            color_min = np.median(np.sort(df[color_coding].values)[:five_per_cent])\n",
    "\n",
    "\n",
    "        try:\n",
    "            color_max = float(OPTIONAL_CONFIG['color_max'])\n",
    "\n",
    "        except:\n",
    "\n",
    "            if catalog == 'NEA':\n",
    "\n",
    "                five_per_cent = int(len(df[color_coding].values)*0.05)\n",
    "                color_max = np.median(np.sort(df[color_coding].values)[-five_per_cent:])\n",
    "\n",
    "            if catalog == 'PlanetS':\n",
    "\n",
    "                five_per_cent = int(len(df[color_coding_S].values)*0.05)\n",
    "                color_max = np.median(np.sort(df[color_coding_S].values)[-five_per_cent:])\n",
    "    else:\n",
    "        \n",
    "        pass\n",
    "\n",
    "        \n",
    "    \n",
    "#@|############--log_x and log_y--##########\n",
    "\n",
    "try:\n",
    "    log_x = OPTIONAL_CONFIG['log_x'] == 'True'\n",
    "except:\n",
    "    log_x = True\n",
    "    \n",
    "try:\n",
    "    log_y = OPTIONAL_CONFIG['log_y'] == 'True'\n",
    "except:\n",
    "    log_y = True\n",
    "    \n",
    "\n",
    "#@|###########--xlim and ylim--###############\n",
    "\n",
    "try:\n",
    "    x_lims = [float(x.strip()) for x in OPTIONAL_CONFIG['xlim'].split(',')]\n",
    "except:\n",
    "    x_lims = [0.5, 21]\n",
    "    \n",
    "try:\n",
    "    y_lims = [float(x.strip()) for x in OPTIONAL_CONFIG['ylim'].split(',')]\n",
    "except:\n",
    "    y_lims  = [0.9, 2.8]\n",
    "    \n",
    "    \n",
    "#@|###########--legend location-###############\n",
    "\n",
    "try:\n",
    "    loc_legend = OPTIONAL_CONFIG['loc_legend']\n",
    "except:\n",
    "    loc_legend = 'upper right'\n",
    "    \n",
    "#@|############--plot the low-density super-Earths region--###############\n",
    "\n",
    "try:\n",
    "    low_d_sE = OPTIONAL_CONFIG['low_density_superEarths'] == 'True'\n",
    "except:\n",
    "    low_d_sE = False\n",
    "    \n",
    "    \n",
    "#@|############--Markersize of the NEA and my planets--###############  \n",
    "    \n",
    "try:\n",
    "    size_catalog_planets = int(OPTIONAL_CONFIG['size_catalog_planets'])\n",
    "except:\n",
    "    size_catalog_planets = 35\n",
    "    \n",
    "try:   \n",
    "    size_my_planets = int(OPTIONAL_CONFIG['size_my_planets'])\n",
    "except:\n",
    "    size_my_planets = 150\n",
    "    \n",
    "    \n",
    "#@|############--Grey shade below the 100% iron model by Zeng et al. (2019)--###############\n",
    "\n",
    "try:\n",
    "    shade_below_pure_iron = OPTIONAL_CONFIG['shade_below_pure_iron'] == 'True'\n",
    "except:\n",
    "    shade_below_pure_iron = True\n",
    "    \n",
    "    \n",
    "    \n",
    "#@|############--Dark background--##########\n",
    "\n",
    "try:\n",
    "    dark_background = OPTIONAL_CONFIG['dark_background'] == 'True'\n",
    "except:\n",
    "    dark_background = False\n",
    "    \n",
    "    \n",
    "#@|#########---linewidths of the models-----#############\n",
    "\n",
    "try:   \n",
    "    lw_models = float(OPTIONAL_CONFIG['lw_models'])\n",
    "except:\n",
    "    lw_models = 1.6\n",
    "    \n",
    "\n",
    "#@|########-------numbers of columns in the legend--------######\n",
    "\n",
    "try:   \n",
    "    n_cols_legend = int(OPTIONAL_CONFIG['n_cols_legend'])\n",
    "except:\n",
    "    n_cols_legend = 1\n",
    "    \n",
    "    \n",
    "#@|###########----capsizes-------#############\n",
    "\n",
    "try:\n",
    "    capsize = float(OPTIONAL_CONFIG['capsize'])\n",
    "except:\n",
    "    capsize = 0\n",
    "    \n",
    "    \n",
    "#@|###########----appearance-------#############  \n",
    "\n",
    "try:\n",
    "    appearance = str(OPTIONAL_CONFIG['appearance'])\n",
    "except:\n",
    "    appearance = 'standard'\n",
    "    \n",
    "    \n",
    "    \n",
    "#@|###########----ec_catalog for the plots with no color code-------#############  \n",
    "\n",
    "try:\n",
    "    ec_catalog = str(OPTIONAL_CONFIG['ec_catalog'])\n",
    "except:\n",
    "    ec_catalog = 'grey' \n",
    "    \n",
    "#@|###########----ec_catalog for the color coded plots-------#############  \n",
    "\n",
    "try:\n",
    "    ec_catalog_cc = str(OPTIONAL_CONFIG['ec_catalog_cc'])\n",
    "except:\n",
    "    ec_catalog_cc = 'whitesmoke'    \n",
    "    \n",
    "    \n",
    "#@|###########--My planets--###############   \n",
    "    \n",
    "text_boxes = False    \n",
    "\n",
    "mass_upper_limit = []\n",
    "\n",
    "for i in range(1, 21):\n",
    "    try:\n",
    "        globals()['m_p'+str(i)] = float(MY_DATA['m_p'+str(i)])\n",
    "        globals()['r_p'+str(i)] = float(MY_DATA['r_p'+str(i)])\n",
    "        globals()['c_p'+str(i)] = MY_DATA['c_p'+str(i)]\n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "    #|@uncertainties\n",
    "    \n",
    "    try:\n",
    "        globals()['m_p'+str(i)+'_err_up'] = float(MY_DATA['m_p'+str(i)+'_err_up'])\n",
    "        globals()['m_p'+str(i)+'_err_down'] = float(MY_DATA['m_p'+str(i)+'_err_down'])\n",
    "        globals()['r_p'+str(i)+'_err_up'] = float(MY_DATA['r_p'+str(i)+'_err_up'])\n",
    "        globals()['r_p'+str(i)+'_err_down'] = float(MY_DATA['r_p'+str(i)+'_err_down'])\n",
    "        mass_upper_limit.append(False)\n",
    "        \n",
    "    except:\n",
    "        \n",
    "        try:\n",
    "            globals()['m_p'+str(i)+'_err_down'] = float(MY_DATA['m_p'+str(i)+'_err_down'])\n",
    "            globals()['r_p'+str(i)+'_err_up'] = float(MY_DATA['r_p'+str(i)+'_err_up'])\n",
    "            globals()['r_p'+str(i)+'_err_down'] = float(MY_DATA['r_p'+str(i)+'_err_down'])\n",
    "            mass_upper_limit.append(True)\n",
    "            \n",
    "        except:\n",
    "            pass\n",
    "    \n",
    "    #@|text boxes\n",
    "    \n",
    "    try:\n",
    "        globals()['name_p'+str(i)] = MY_DATA['name_p'+str(i)]\n",
    "        #@|location of the text boxes\n",
    "        globals()['dis_x_p'+str(i)] = float(MY_DATA['dis_x_p'+str(i)])\n",
    "        globals()['dis_y_p'+str(i)] = float(MY_DATA['dis_y_p'+str(i)])\n",
    "        text_boxes = True  \n",
    "    except:\n",
    "        pass\n",
    "    \n",
    "mass_upper_limit = np.array(mass_upper_limit)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1521b461",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc98ffe4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_low_density_SE_region():\n",
    "    \n",
    "    x = [1.55, 1.58,  2., 2.5, 3.5, 2.6, 1.9]\n",
    "    y = [1.35, 1.5, 1.67, 1.68, 1.6, 1.38, 1.30]\n",
    "    #ax.scatter(x, y, zorder=3)\n",
    "\n",
    "    tck, u = interpolate.splprep([x + x[:1], y + y[:1]], s=0, per=True)\n",
    "    unew = np.linspace(0, 1, 100)\n",
    "    basic_form = interpolate.splev(unew, tck)\n",
    "    ax.plot(basic_form[0], basic_form[1], color='blue', lw=2,alpha = 0.7)\n",
    "    ax.fill(basic_form[0], basic_form[1], color='dodgerblue', alpha=0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e177243",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0925fe8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "###########################\n",
    "#@|--------Main---------#@|\n",
    "###########################\n",
    "\n",
    "#@|plot config_object\n",
    "zorder_lines = -100\n",
    "\n",
    "\n",
    "ec_my_planets = 'black'\n",
    "  \n",
    "    \n",
    "if dark_background:\n",
    "    plt.style.use(\"dark_background\")\n",
    "    ec_my_planets = 'white'\n",
    "    zeng_models_colors['zeng_2019_pure_iron'] = 'white'\n",
    "\n",
    "try:\n",
    "    cmap = OPTIONAL_CONFIG['cmap']\n",
    "except:\n",
    "    cmap = 'rainbow'\n",
    "\n",
    "\n",
    "#@|+++++++++++++++++++++++ \n",
    "#@|+++++Color coding++++++            \n",
    "#@|+++++++++++++++++++++++ \n",
    "\n",
    "    \n",
    "if color_coding in cols:\n",
    "    if n_cols == 'one':\n",
    "    \n",
    "        plt.figure(figsize = (10, 6.56))\n",
    "        \n",
    "    if n_cols == 'two':\n",
    "        \n",
    "        plt.figure(figsize = (20, 6.56))\n",
    "        \n",
    "    \n",
    "    #@|+++++++++++++++++++++++++++++++++ \n",
    "    #@|+++++Numerical Color coding++++++            \n",
    "    #@|+++++++++++++++++++++++++++++++++ \n",
    "    \n",
    "    if numerical_color_coding:\n",
    "\n",
    "    \n",
    "        plt.scatter(M_catalog, R_catalog, c = df[color_coding].values, cmap = cmap, s = size_catalog_planets, \\\n",
    "                    vmin = color_min, vmax = color_max, lw = 1.1, ec = ec_catalog_cc)\n",
    "        \n",
    "\n",
    "\n",
    "        if n_cols == 'one':\n",
    "            cbar = plt.colorbar(location = 'left', anchor=(2.95,0.95), aspect = 8, shrink=0.4)\n",
    "        if n_cols == 'two':\n",
    "            cbar = plt.colorbar(location = 'left', anchor=(2.3,0.95), aspect = 8, shrink=0.4)\n",
    "\n",
    "        \n",
    "        #@|plt.errorbar(M_catalog, R_catalog, yerr = np.array(R_catalog_ERR).T, xerr = np.array(M_catalog_ERR).T,\\\n",
    "                     #@|linestyle = \"None\", ecolor = 'grey', zorder = -1, lw = 1.8, capsize = capsize, alpha = 1)\n",
    "        \n",
    "        #@|++++++++++++++++++++++++++++++++++++++++++++\n",
    "        #@|Error bars of the same color of data points\n",
    "        #@|++++++++++++++++++++++++++++++++++++++++++++\n",
    "        #convert time to a color tuple using the colormap used for scatter\n",
    "        import matplotlib\n",
    "        import matplotlib.cm as cm\n",
    "        norm = matplotlib.colors.Normalize(vmin=color_min, vmax=color_max, clip=True)\n",
    "        mapper = cm.ScalarMappable(norm=norm, cmap=cmap)\n",
    "        time_color = np.array([(mapper.to_rgba(v)) for v in df[color_coding].values])\n",
    "\n",
    "        #loop over each data point to plot\n",
    "        for x, y, e1, e2, color in zip(M_catalog, R_catalog, \\\n",
    "                                       (abs(R_catalog_err_up) + abs(R_catalog_err_down)) / 2, \\\n",
    "                                       (abs(M_catalog_err_up) + abs(M_catalog_err_down)) / 2, time_color):\n",
    "            #plt.scatter(x, y, color=color)\n",
    "            plt.errorbar(x, y, yerr = e1, xerr=e2, lw = 1.8, capsize = capsize, color = color, zorder = -1)\n",
    "            \n",
    "            \n",
    "        #@|++++Color bar++++#@|    \n",
    "        cbar.set_label(color_coding_labels[color_coding], fontsize = 12)\n",
    "        cbar.ax.tick_params(labelsize=10)\n",
    "        #cbar.set_ticks_position('left')\n",
    "\n",
    "\n",
    "        #@|My planets\n",
    "        for i in range(1, 21):\n",
    "            try:\n",
    "                plt.scatter(globals()['m_p'+str(i)], globals()['r_p'+str(i)],\\\n",
    "                c = float(globals()['c_p'+str(i)]), ec = ec_my_planets, cmap = cmap, s = 250, \\\n",
    "                vmin = color_min, vmax = color_max, marker = 'h',lw = 1.7, zorder = 1000)\n",
    "                plt.errorbar(globals()['m_p'+str(i)], globals()['r_p'+str(i)],\\\n",
    "                 xerr = np.array([[globals()['m_p'+str(i)+'_err_down']],[globals()['m_p'+str(i)+'_err_up']]]), \\\n",
    "                 yerr = np.array([[globals()['r_p'+str(i)+'_err_down']],[globals()['r_p'+str(i)+'_err_up']]]), \\\n",
    "                 zorder = 999, lw = 1.7, capsize = capsize, color = ec_my_planets)\n",
    "\n",
    "\n",
    "                if text_boxes == True:\n",
    "\n",
    "                    plt.text(globals()['m_p'+str(i)] +  globals()['dis_x_p'+str(i)], \\\n",
    "                     globals()['r_p'+str(i)] +  globals()['dis_y_p'+str(i)], globals()['name_p'+str(i)], fontsize=12.5, \\\n",
    "                    verticalalignment = 'bottom', horizontalalignment = 'center', \\\n",
    "                     bbox = {\"boxstyle\": \"round\", 'facecolor': 'k', 'alpha': 0.7, 'pad': 0.4}, \\\n",
    "                     c = 'white', weight='bold', zorder = 10000)\n",
    "\n",
    "            except:\n",
    "                pass\n",
    "            \n",
    "            \n",
    "    #@|+++++++++++++++++++++++++++++++++++++ \n",
    "    #@|+++++Non numerical Color coding++++++            \n",
    "    #@|+++++++++++++++++++++++++++++++++++++         \n",
    "            \n",
    "    else:  \n",
    "        \n",
    "        if n_cols == 'one':\n",
    "    \n",
    "            plt.figure(figsize = (8, 7))\n",
    "\n",
    "        if n_cols == 'two':\n",
    "\n",
    "            plt.figure(figsize = (16, 7))\n",
    "            \n",
    "        \n",
    "        for i in range(len(groups)):\n",
    "\n",
    "            idxs_ac = np.where(df[color_coding].values==groups[i])[0] \n",
    "\n",
    "            \n",
    "            #@|change for the labels\n",
    "            if groups[i] == 'Transiting Exoplanet Survey Satellite (TESS)':\n",
    "                groups[i] = 'TESS'\n",
    "            \n",
    "            \n",
    "            plt.scatter(M_catalog[idxs_ac], R_catalog[idxs_ac], c = colors_groups[i], \\\n",
    "                        s = size_catalog_planets, lw = 1.1, ec = ec_catalog_cc, label = groups[i], zorder = i)\n",
    "            \n",
    "            plt.errorbar(M_catalog[idxs_ac], R_catalog[idxs_ac], yerr = np.array(R_catalog_ERR)[idxs_ac].T,\\\n",
    "                         xerr = np.array(M_catalog_ERR)[idxs_ac].T,\\\n",
    "                     linestyle = \"None\", ecolor = colors_groups[i], zorder = i-1, lw = 2, capsize = capsize)\n",
    "                     \n",
    "            \n",
    "        if plot_all_planets:\n",
    "            \n",
    "            plt.scatter(M_catalog, R_catalog, c = 'lightgrey', s = size_catalog_planets, \\\n",
    "                        ec = ec_catalog_cc, lw = 0.7, alpha = 0.6, zorder = -5, label  = 'Other')\n",
    "            plt.errorbar(M_catalog, R_catalog, yerr = np.array(R_catalog_ERR).T,\\\n",
    "                         xerr = np.array(M_catalog_ERR).T, \\\n",
    "                         linestyle = \"None\", ecolor = 'grey', \\\n",
    "                         zorder = -6, lw = 1.3, capsize = capsize, alpha = 0.6)\n",
    "            \n",
    "            \n",
    "        #@|My planets (no color coding)\n",
    "    \n",
    "        for i in range(1, 21):\n",
    "            try:\n",
    "\n",
    "                if text_boxes == True:\n",
    "\n",
    "                    plt.scatter(globals()['m_p'+str(i)], globals()['r_p'+str(i)],\\\n",
    "                                color = globals()['c_p'+str(i)], ec = ec_my_planets,\\\n",
    "                                marker = 'h', s = size_my_planets, lw = 1.7, zorder = 1000,)\n",
    "                    plt.errorbar(globals()['m_p'+str(i)], globals()['r_p'+str(i)],\\\n",
    "                                 xerr = np.array([[globals()['m_p'+str(i)+'_err_down']],\\\n",
    "                                                  [globals()['m_p'+str(i)+'_err_up']]]), \\\n",
    "                                 yerr = np.array([[globals()['r_p'+str(i)+'_err_down']],\\\n",
    "                                                  [globals()['r_p'+str(i)+'_err_up']]]), \\\n",
    "                                 zorder = 999, lw = 1.7, capsize = capsize, color = globals()['c_p'+str(i)])\n",
    "\n",
    "\n",
    "                    plt.text(globals()['m_p'+str(i)] +  globals()['dis_x_p'+str(i)], \\\n",
    "                             globals()['r_p'+str(i)] +  globals()['dis_y_p'+str(i)], \\\n",
    "                             globals()['name_p'+str(i)], fontsize = 12.5, \\\n",
    "                            verticalalignment = 'bottom', horizontalalignment = 'center', \\\n",
    "                             bbox = {\"boxstyle\": \"round\",'facecolor': globals()['c_p'+str(i)], \\\n",
    "                                     'alpha': 0.8, 'pad': 0.4}, c = 'white', weight='bold', zorder = 10000)\n",
    "\n",
    "                if text_boxes == False:\n",
    "\n",
    "                    plt.scatter(globals()['m_p'+str(i)], globals()['r_p'+str(i)],\\\n",
    "                                color = globals()['c_p'+str(i)], ec = ec_my_planets,\\\n",
    "                                marker = 'h', s = size_my_planets, lw = 1.7,\\\n",
    "                                zorder = 1000, label = globals()['name_p'+str(i)])\n",
    "                    plt.errorbar(globals()['m_p'+str(i)], globals()['r_p'+str(i)],\\\n",
    "                                 xerr = np.array([[globals()['m_p'+str(i)+'_err_down']],\\\n",
    "                                                  [globals()['m_p'+str(i)+'_err_up']]]), \\\n",
    "                                 yerr = np.array([[globals()['r_p'+str(i)+'_err_down']],\\\n",
    "                                                  [globals()['r_p'+str(i)+'_err_up']]]), \\\n",
    "                                 zorder = 999, lw = 1.7, capsize = capsize, color = globals()['c_p'+str(i)])   \n",
    "\n",
    "            except:\n",
    "            \n",
    "                pass\n",
    "            \n",
    "\n",
    "#@|+++++++++++++++++++++++ \n",
    "#@|++++No color coding++++            \n",
    "#@|+++++++++++++++++++++++ \n",
    "    \n",
    "elif color_coding == 'none':\n",
    "    if n_cols == 'one':\n",
    "    \n",
    "        plt.figure(figsize = (8, 7))\n",
    "        \n",
    "    if n_cols == 'two':\n",
    "        \n",
    "        plt.figure(figsize = (16, 7))\n",
    "        \n",
    "    #@|My planets \n",
    "    \n",
    "    for i in range(1, 21):\n",
    "        try:\n",
    "            \n",
    "            if text_boxes == True:\n",
    "                \n",
    "                plt.scatter(globals()['m_p'+str(i)], globals()['r_p'+str(i)],\\\n",
    "                            color = globals()['c_p'+str(i)], ec = ec_my_planets,\\\n",
    "                            marker = 'h', s = size_my_planets, lw = 1.7, zorder = 1000,)\n",
    "                plt.errorbar(globals()['m_p'+str(i)], globals()['r_p'+str(i)],\\\n",
    "                             xerr = np.array([[globals()['m_p'+str(i)+'_err_down']],\\\n",
    "                                              [globals()['m_p'+str(i)+'_err_up']]]), \\\n",
    "                             yerr = np.array([[globals()['r_p'+str(i)+'_err_down']],\\\n",
    "                                              [globals()['r_p'+str(i)+'_err_up']]]), \\\n",
    "                             zorder = 999, lw = 1.7, capsize = capsize, color = 'k')\n",
    "                \n",
    "                \n",
    "                plt.text(globals()['m_p'+str(i)] +  globals()['dis_x_p'+str(i)], \\\n",
    "                         globals()['r_p'+str(i)] +  globals()['dis_y_p'+str(i)], \\\n",
    "                         globals()['name_p'+str(i)], fontsize = 12.5, \\\n",
    "                        verticalalignment = 'bottom', horizontalalignment = 'center', \\\n",
    "                         bbox = {\"boxstyle\": \"round\",'facecolor': globals()['c_p'+str(i)], \\\n",
    "                                 'alpha': 0.8, 'pad': 0.4}, c = 'white', weight='bold', zorder = 10000)\n",
    "                \n",
    "            if text_boxes == False:\n",
    "                \n",
    "                plt.scatter(globals()['m_p'+str(i)], globals()['r_p'+str(i)],\\\n",
    "                            color = globals()['c_p'+str(i)], ec = ec_my_planets,\\\n",
    "                            marker = 'h', s = size_my_planets, lw = 1.7,\\\n",
    "                            zorder = 1000, label = globals()['name_p'+str(i)])\n",
    "                plt.errorbar(globals()['m_p'+str(i)], globals()['r_p'+str(i)],\\\n",
    "                             xerr = np.array([[globals()['m_p'+str(i)+'_err_down']],\\\n",
    "                                              [globals()['m_p'+str(i)+'_err_up']]]), \\\n",
    "                             yerr = np.array([[globals()['r_p'+str(i)+'_err_down']],\\\n",
    "                                              [globals()['r_p'+str(i)+'_err_up']]]), \\\n",
    "                             zorder = 999, lw = 1.7, capsize = capsize, color = 'k')\n",
    "                \n",
    "                \n",
    "   \n",
    " \n",
    "        except:\n",
    "            \n",
    "            pass\n",
    "        \n",
    "    #@|Catalog planets\n",
    "    \n",
    "    if appearance == 'standard':\n",
    "    \n",
    "        plt.scatter(M_catalog, R_catalog, c = 'lightgrey', s = size_catalog_planets, ec = ec_catalog, \\\n",
    "                    lw = 0.7, alpha = 1, zorder = 0)\n",
    "        plt.errorbar(M_catalog, R_catalog, yerr = np.array(R_catalog_ERR).T, xerr = np.array(M_catalog_ERR).T, \\\n",
    "                     linestyle = \"None\", ecolor = 'grey', zorder = -1, lw = 1.45, capsize = capsize, alpha = 1)\n",
    "     \n",
    "    if appearance == 'faint':\n",
    "        plt.scatter(M_catalog, R_catalog, c = 'lightgrey', s = size_catalog_planets, ec = ec_catalog, \\\n",
    "                    lw = 0.7, alpha = 0.5, zorder = 0)\n",
    "        plt.errorbar(M_catalog, R_catalog, yerr = np.array(R_catalog_ERR).T, xerr = np.array(M_catalog_ERR).T, \\\n",
    "                     linestyle = \"None\", ecolor = 'grey', zorder = -1, lw = 1.45, capsize = capsize, alpha = 0.2)\n",
    "        \n",
    "        \n",
    "        \n",
    "\n",
    "elif color_coding not in cols:\n",
    "    \n",
    "    raise Exception('The selected color code ('+color_coding+') is not availabe in the '+catalog+' catalog.')\n",
    "        \n",
    "#@|##########--Theoretical models--###############\n",
    "\n",
    "#@|Zeng et al. (2016, 2019)\n",
    "\n",
    "if zeng:\n",
    "\n",
    "    for model in models_zeng:\n",
    "        \n",
    "        if model == 'zeng_2019_pure_iron':\n",
    "            \n",
    "            plt.plot(pd.read_csv(path_models+model,  sep = '\\t', header = None)[0],\\\n",
    "                 pd.read_csv(path_models+model,  sep = '\\t', header = None)[1],\\\n",
    "                     lw = lw_models, linestyle = \"solid\" , zorder = -100,\\\n",
    "                     c = zeng_models_colors[model], label =  zeng_models_labels[model])\n",
    "            \n",
    "            \n",
    "            \n",
    "        elif model == 'zeng_2019_earth_like' or model == 'zeng_2016_20_Fe' or model == 'zeng_2019_pure_rock':\n",
    "            \n",
    "            plt.plot(pd.read_csv(path_models+model,  sep = '\\t', header = None)[0],\\\n",
    "                 pd.read_csv(path_models+model,  sep = '\\t', header = None)[1],\\\n",
    "                     lw = lw_models, linestyle = \"dashdot\" , zorder = -100,\\\n",
    "                     c = zeng_models_colors[model], label =  zeng_models_labels[model])\n",
    "            \n",
    "        elif model.find('H2O') != -1:\n",
    "            \n",
    "            plt.plot(pd.read_csv(path_models+model,  sep = '\\t', header = None)[0],\\\n",
    "                 pd.read_csv(path_models+model,  sep = '\\t', header = None)[1],\\\n",
    "                     lw = lw_models, linestyle = \"dashed\" , zorder = -100,\\\n",
    "                     c = zeng_models_colors[model], label =  zeng_models_labels[model])\n",
    "            \n",
    "           \n",
    "        else:\n",
    "        \n",
    "            plt.plot(pd.read_csv(path_models+model,  sep = '\\t', header = None)[0],\\\n",
    "                 pd.read_csv(path_models+model,  sep = '\\t', header = None)[1],\\\n",
    "                     lw = lw_models*1.4, linestyle = \"dotted\" , zorder = -100,\\\n",
    "                     c = zeng_models_colors[model], label =  zeng_models_labels[model])\n",
    "        \n",
    "    \n",
    "#@|Turbet et al. (2020)\n",
    "if turb_2020:\n",
    "    \n",
    "    for i in range(len(WMFs_turb2020)):\n",
    "\n",
    "        if Core_turb2020[i] == 'earth':\n",
    "            label_core = 'Earth'\n",
    "        elif Core_turb2020[i] == 'rock':\n",
    "            label_core = 'Rock'\n",
    "        elif Core_turb2020[i] == 'iron':\n",
    "            label_core = 'Iron'\n",
    "\n",
    "        label_wmf = WMFs_turb2020[i]\n",
    "\n",
    "        plt.plot(M_pl_turb[i], R_pl_turb[i], lw = lw_models*1.2, linestyle = 'dotted',\\\n",
    "                 zorder = -1000, c=colors_turb2020[i], label = label_core + \" + \"+ str(label_wmf*100)+\"% \"+r'$\\rm H_{2}O$ steam')\n",
    "\n",
    "        \n",
    "#@|Aguichine et al. (2021)\n",
    "if aguich2021:\n",
    "\n",
    "    for i in range(len(x_core_aguich2021)):\n",
    "\n",
    "        plt.plot(M_aguich2021[i], R_aguich2021[i], lw = lw_models*1.2, linestyle = 'dotted', zorder = -1000, \\\n",
    "                label = str(int(x_core_aguich2021[i]*100))+'% CMF & '+str(int(x_H2O_aguich2021[i]*100))+\\\n",
    "                 r'% WMF ('+str(int(Tirr_aguich2021[i]))+'K)', c = colors_aguich2021[i])\n",
    "        \n",
    "#@|Lopez & Fortney et al. (2014)\n",
    "if lopez_fortney2014:\n",
    "    \n",
    "    for i in range(len(age_lf2014)):\n",
    "        \n",
    "        plt.plot(M_lf2014[i], R_lf2014[i], lw = lw_models*1.2, linestyle = 'dotted', zorder = -1000, \\\n",
    "                label = H_He[i]+'% H/He, '+age_lf2014[i]+','+f' {Seff_lf2014[i]}'+r'$\\rm S_{\\oplus}$', \\\n",
    "                c = colors_lf2014[i])\n",
    "        \n",
    "        \n",
    "#@|Haldemannet al. (2014) \n",
    "if haldemann2024:\n",
    "    \n",
    "    for i in range(len(models_haldemann2024)):\n",
    "        \n",
    "        #@|solid lines for the 100% core planets\n",
    "        \n",
    "        if models_haldemann2024[i] == 'C0' or models_haldemann2024[i] == 'C1':\n",
    "            \n",
    "            plt.plot(M_haldemann2024[i], R_haldemann2024[i], lw = lw_models*1.2, linestyle = 'solid', zorder = -1000, \\\n",
    "                label = models_haldemann2024[i]+ \\\n",
    "                 r' ($T_{\\rm out}$ = ' +str(int(T_out_haldemann2024[i]))+' K)', \\\n",
    "                c = colors_haldemann2024[i])\n",
    "        \n",
    "        else:\n",
    "        \n",
    "            plt.plot(M_haldemann2024[i], R_haldemann2024[i], lw = lw_models, linestyle = 'dashed', zorder = -1000, \\\n",
    "                    label = models_haldemann2024[i]+ \\\n",
    "                     r' ($T_{\\rm out}$ = ' +str(int(T_out_haldemann2024[i]))+' K)', \\\n",
    "                    c = colors_haldemann2024[i])\n",
    "            \n",
    "#@|Seager et al. (2007)            \n",
    "if seager2007:\n",
    "    \n",
    "    for i in range(len(models_seager2007)):\n",
    "        \n",
    "        if models_seager2007[i] == 'iron':\n",
    "            \n",
    "            plt.plot(M_seager2007[i], R_seager2007[i], lw = lw_models*1.2, linestyle = 'solid', zorder = -1000, \\\n",
    "                label = r'100% Fe (0% $\\rm MgSiO_{3}$)', c = colors_seager2007[i])\n",
    "            \n",
    "        if models_seager2007[i] == 'rock':\n",
    "            \n",
    "            plt.plot(M_seager2007[i], R_seager2007[i], lw = lw_models*1.2, linestyle = 'dashdot', zorder = -1000, \\\n",
    "                label = r'0% Fe (100% $\\rm MgSiO_{3}$)', c = colors_seager2007[i])\n",
    "            \n",
    "        if models_seager2007[i] == 'water':\n",
    "            \n",
    "            plt.plot(M_seager2007[i], R_seager2007[i], lw = lw_models*1.2, linestyle = 'dashed', zorder = -1000, \\\n",
    "                label = r'100% $\\rm H_{2}O$', c = colors_seager2007[i])\n",
    "            \n",
    "            \n",
    "            \n",
    "#@|Marcus et al. (2010)\n",
    "if marcus:\n",
    "\n",
    "    for model in models_marcus:\n",
    "        plt.plot(pd.read_csv(path_models+model,  sep = '\\t', header = None)[0],\\\n",
    "                 pd.read_csv(path_models+model,  sep = '\\t', header = None)[1], lw = lw_models*1.4, linestyle = \"dotted\" ,\\\n",
    "                 zorder = -100, c = zeng_models_colors[model], label =  zeng_models_labels[model] )\n",
    "\n",
    "        \n",
    "#@|Isodensity lines\n",
    "if isodensity:\n",
    "\n",
    "    iso_density()\n",
    "\n",
    "    #for i in range(len(Masses_iso)):\n",
    "        #plt.plot(Masses_iso[i], Radii_iso[i], lw = lw_models, linestyle = '-.', zorder = -1000,\\\n",
    "                #label = fr'Isodensity ({density[i]/1000} $\\rm g / cm^3$)', c = colors_density[i])   \n",
    "        \n",
    "    for i in range(len(Masses_iso)):\n",
    "        plt.plot(Masses_iso[i], Radii_iso[i], lw = lw_models, linestyle = '-.', zorder = -1000,\\\n",
    "                label = fr'{density[i]/1000} $\\rm g / cm^3$ isodensity', c = colors_density[i]) \n",
    "\n",
    "        \n",
    "#@|Otegi et al. (2020)\n",
    "if otegi2020:\n",
    "    \n",
    "    otegi2020_MR()\n",
    "    \n",
    "    \n",
    "    if 'small' in models_otegi2020:\n",
    "        \n",
    "        idx = np.where(np.array(models_otegi2020)=='small')[0][0]\n",
    "        plt.plot(M_otegi_small, R_otegi_small, lw = 2.5, c = colors_otegi2020[idx], \\\n",
    "                 linestyle = linestyles_otegi2020[idx])\n",
    "        \n",
    "    if 'intermediate' in models_otegi2020:\n",
    "        \n",
    "        idx = np.where(np.array(models_otegi2020)=='intermediate')[0][0]\n",
    "        plt.plot(M_otegi_intermediate, R_otegi_intermediate, lw = 2.5, c = colors_otegi2020[idx], \\\n",
    "                 linestyle =  linestyles_otegi2020[idx])\n",
    "        \n",
    "        \n",
    "    #@|vertical lines (no! to remove)   \n",
    "    #plt.vlines(x = [10, 138],  ymin = [1.25, 6], ymax = [4.9, 25], linestyle = 'dotted', color = 'k', lw = 1.7)        \n",
    "        \n",
    "        \n",
    "#@|Parc et al. (2024)\n",
    "if parc2024:\n",
    "    \n",
    "    parc2024_MR()\n",
    "    \n",
    "    if 'small' in models_parc2024:\n",
    "        \n",
    "        idx = np.where(np.array(models_parc2024)=='small')[0][0]\n",
    "        plt.plot(M_parc_small, R_parc_small, lw = 2.5, c = colors_parc2024[idx], \\\n",
    "                 linestyle = linestyles_parc2024[idx])\n",
    "        \n",
    "    if 'intermediate' in models_parc2024:\n",
    "        \n",
    "        idx = np.where(np.array(models_parc2024)=='intermediate')[0][0]\n",
    "        plt.plot(M_parc_intermediate, R_parc_intermediate, lw = 2.5, c = colors_parc2024[idx], \\\n",
    "                 linestyle =  linestyles_parc2024[idx])\n",
    "        \n",
    "    if 'giant' in models_parc2024:\n",
    "        \n",
    "        idx = np.where(np.array(models_parc2024)=='giant')[0][0]\n",
    "        plt.plot(M_parc_giant, R_parc_giant, lw = 2.5, c = colors_parc2024[idx], \\\n",
    "                 linestyle = linestyles_parc2024[idx])\n",
    "        \n",
    "        \n",
    "    #@|vertical lines (only between intermediate and giant planets)\n",
    "    #plt.vlines(x = [10, 138],  ymin = [1.25, 6], ymax = [4.9, 25], linestyle = 'dotted', color = 'k', lw = 1.7)\n",
    "    plt.vlines(x = [138],  ymin = [6], ymax = [25], linestyle = 'dotted', color = 'k', lw = 1.7)\n",
    "    \n",
    "    \n",
    "#@configuration\n",
    "import matplotlib.ticker as mticker\n",
    "ax = plt.gca()\n",
    "if log_x == True:\n",
    "    ax.set_xscale('log')\n",
    "if log_y == True:\n",
    "    ax.set_yscale('log')\n",
    "    \n",
    "\n",
    "ax.xaxis.set_major_formatter(mticker.ScalarFormatter())\n",
    "ax.yaxis.set_minor_formatter(mticker.ScalarFormatter())\n",
    "ax.yaxis.set_major_formatter(mticker.ScalarFormatter())\n",
    "\n",
    "from matplotlib.ticker import FormatStrFormatter\n",
    "ax.yaxis.set_minor_formatter(FormatStrFormatter('%.1f'))\n",
    "ax.yaxis.set_major_formatter(FormatStrFormatter('%.1f'))\n",
    "\n",
    "\n",
    "ax.tick_params(which='major', length = 9, color = ec_my_planets, direction = 'in', width = 1.3)\n",
    "ax.tick_params(which='minor', length = 5, color = ec_my_planets, direction = 'in', width = 0.9)\n",
    "\n",
    "\n",
    "\n",
    "if color_coding == 'none':\n",
    "\n",
    "    ax.tick_params(axis='both', which='major', labelsize=14)\n",
    "    ax.tick_params(axis='both', which='minor', labelsize=14)\n",
    "\n",
    "    plt.ylabel(r'Planet radius $\\rm (R_{\\oplus})$', fontsize = 16)\n",
    "    plt.xlabel(r'Planet mass $\\rm (M_{\\oplus})$', fontsize = 16)    \n",
    "    \n",
    "else:\n",
    "    \n",
    "    ax.tick_params(axis='both', which='major', labelsize=13)\n",
    "    ax.tick_params(axis='both', which='minor', labelsize=13)\n",
    "\n",
    "    plt.ylabel(r'Planet radius $\\rm (R_{\\oplus})$', fontsize = 15)\n",
    "    plt.xlabel(r'Planet mass $\\rm (M_{\\oplus})$', fontsize = 15)   \n",
    "    \n",
    "plt.xlim(x_lims)\n",
    "plt.ylim(y_lims)\n",
    "\n",
    "\n",
    "if n_cols == 'one':\n",
    "    plt.legend(loc = loc_legend, fontsize = 10, markerscale = 1, ncol = n_cols_legend)\n",
    "if n_cols == 'two':\n",
    "    plt.legend(loc = loc_legend, fontsize = 12.5, ncol = n_cols_legend)\n",
    "\n",
    "#@|low density super-Earths region\n",
    "\n",
    "if low_d_sE:\n",
    "    plot_low_density_SE_region()\n",
    "    \n",
    "#@| Region below the 100% iron model by Zeng et al. (2019)\n",
    "    \n",
    "if shade_below_pure_iron:\n",
    "    \n",
    "    if haldemann2024:\n",
    "        \n",
    "        if 'C0' in models_haldemann2024:\n",
    "\n",
    "            idxs_C0 = np.where(np.array(models_haldemann2024) == 'C0')[0]\n",
    "            T_out_min = np.min(np.array(T_out_haldemann2024)[idxs_C0])\n",
    "            \n",
    "            df_haldemann2024 = pd.read_csv(f'{path_models}/Mass_Radius_C0.dat', \\\n",
    "                               comment = '#', header = None)\n",
    "            \n",
    "        if 'C1' in models_haldemann2024 and 'C0' not in models_haldemann2024:\n",
    "            \n",
    "            idxs_C1 = np.where(np.array(models_haldemann2024) == 'C1')[0]\n",
    "            T_out_min = np.min(np.array(T_out_haldemann2024)[idxs_C1])\n",
    "            \n",
    "            df_haldemann2024 = pd.read_csv(f'{path_models}/Mass_Radius_C1.dat', \\\n",
    "                               comment = '#', header = None)\n",
    "            \n",
    "            \n",
    "        if T_out_min == 50:\n",
    "            idx_R_haldemann2024 = 1\n",
    "\n",
    "        if T_out_min == 300:\n",
    "            idx_R_haldemann2024 = 2\n",
    "\n",
    "        if T_out_min == 800:\n",
    "            idx_R_haldemann2024 = 3\n",
    "\n",
    "        if T_out_min == 1500:\n",
    "            idx_R_haldemann2024 = 4\n",
    "\n",
    "        if T_out_min == 2000:\n",
    "            idx_R_haldemann2024 = 5\n",
    "\n",
    "        m_pure_iron = df_haldemann2024[0].values\n",
    "        r_pure_iron = df_haldemann2024[idx_R_haldemann2024].values\n",
    "\n",
    "        plt.fill_between(m_pure_iron, r_pure_iron, 0,  alpha = 0.6, color = 'grey', zorder = -1000)\n",
    "            \n",
    "    \n",
    "    else:\n",
    "        \n",
    "        #@|Zeng et al. (2019)\n",
    "        r_pure_iron = pd.read_csv(path_models+'zeng_2019_pure_iron',  sep = '\\t', header = None)[1].values\n",
    "        m_pure_iron = pd.read_csv(path_models+'zeng_2019_pure_iron',  sep = '\\t', header = None)[0].values\n",
    "        plt.fill_between(m_pure_iron, r_pure_iron, 0,  alpha = 0.6, color = 'grey', zorder = -1000)\n",
    "    \n",
    "\n",
    "\n",
    "plt.savefig(f'output/{config_file[:-4]}.pdf', bbox_inches = 'tight', pad_inches = 0.1)\n",
    "plt.savefig(f'output/{config_file[:-4]}.png', bbox_inches = 'tight', pad_inches = 0.1, dpi = 400)\n",
    "\n",
    "\n",
    "print('')\n",
    "print('\\033[1m' + f'Your diagram {config_file[:-4]}.pdf/png has been successfully generated and saved')\n",
    "print('\\033[0m')\n",
    "print('     ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')\n",
    "print('      Thank you for using Mister plotter :-D We hope to see you again soon!')\n",
    "print('     ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')\n",
    "print('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9275869a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3861c9d9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8a0848a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}