Detectability_and_Probability.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyO743kg1I+mb0NubdW9GoYF",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/liamhall64/Habitable-Exomoons/blob/main/Detectability_and_Probability.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "collapsed": true,
        "id": "GK6yhXsRXupD",
        "outputId": "f94dd59b-f35b-4e29-cf0a-f65f61678519"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Collecting pytransit\n",
            "  Downloading PyTransit-2.6.8-py3-none-any.whl.metadata (5.7 kB)\n",
            "Collecting celerite\n",
            "  Downloading celerite-0.4.3-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (2.9 kB)\n",
            "Collecting emcee\n",
            "  Downloading emcee-3.1.6-py2.py3-none-any.whl.metadata (3.0 kB)\n",
            "Collecting corner\n",
            "  Downloading corner-2.2.2-py3-none-any.whl.metadata (2.2 kB)\n",
            "Requirement already satisfied: arviz in /usr/local/lib/python3.10/dist-packages (from pytransit) (0.19.0)\n",
            "Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (from pytransit) (1.26.4)\n",
            "Requirement already satisfied: scipy in /usr/local/lib/python3.10/dist-packages (from pytransit) (1.13.1)\n",
            "Requirement already satisfied: pandas in /usr/local/lib/python3.10/dist-packages (from pytransit) (2.1.4)\n",
            "Requirement already satisfied: xarray in /usr/local/lib/python3.10/dist-packages (from pytransit) (2024.9.0)\n",
            "Requirement already satisfied: tables in /usr/local/lib/python3.10/dist-packages (from pytransit) (3.8.0)\n",
            "Collecting uncertainties (from pytransit)\n",
            "  Downloading uncertainties-3.2.2-py3-none-any.whl.metadata (6.9 kB)\n",
            "Requirement already satisfied: numba in /usr/local/lib/python3.10/dist-packages (from pytransit) (0.60.0)\n",
            "Requirement already satisfied: astropy in /usr/local/lib/python3.10/dist-packages (from pytransit) (6.1.3)\n",
            "Requirement already satisfied: matplotlib in /usr/local/lib/python3.10/dist-packages (from pytransit) (3.7.1)\n",
            "Requirement already satisfied: tqdm in /usr/local/lib/python3.10/dist-packages (from pytransit) (4.66.5)\n",
            "Collecting semantic-version (from pytransit)\n",
            "  Downloading semantic_version-2.10.0-py2.py3-none-any.whl.metadata (9.7 kB)\n",
            "Requirement already satisfied: setuptools in /usr/local/lib/python3.10/dist-packages (from pytransit) (71.0.4)\n",
            "Collecting deprecated (from pytransit)\n",
            "  Downloading Deprecated-1.2.14-py2.py3-none-any.whl.metadata (5.4 kB)\n",
            "Requirement already satisfied: seaborn in /usr/local/lib/python3.10/dist-packages (from pytransit) (0.13.1)\n",
            "Collecting ldtk (from pytransit)\n",
            "  Downloading ldtk-1.8.4-py3-none-any.whl.metadata (992 bytes)\n",
            "Collecting pyopencl (from pytransit)\n",
            "  Downloading pyopencl-2024.2.7-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (4.7 kB)\n",
            "Collecting pyrr (from pytransit)\n",
            "  Downloading pyrr-0.10.3-py3-none-any.whl.metadata (832 bytes)\n",
            "Collecting meepmeep (from pytransit)\n",
            "  Downloading MeepMeep-0.6.3-py3-none-any.whl.metadata (714 bytes)\n",
            "Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib->pytransit) (1.3.0)\n",
            "Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.10/dist-packages (from matplotlib->pytransit) (0.12.1)\n",
            "Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib->pytransit) (4.53.1)\n",
            "Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib->pytransit) (1.4.7)\n",
            "Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib->pytransit) (24.1)\n",
            "Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib->pytransit) (10.4.0)\n",
            "Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib->pytransit) (3.1.4)\n",
            "Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.10/dist-packages (from matplotlib->pytransit) (2.8.2)\n",
            "Requirement already satisfied: dm-tree>=0.1.8 in /usr/local/lib/python3.10/dist-packages (from arviz->pytransit) (0.1.8)\n",
            "Requirement already satisfied: h5netcdf>=1.0.2 in /usr/local/lib/python3.10/dist-packages (from arviz->pytransit) (1.3.0)\n",
            "Requirement already satisfied: typing-extensions>=4.1.0 in /usr/local/lib/python3.10/dist-packages (from arviz->pytransit) (4.12.2)\n",
            "Requirement already satisfied: xarray-einstats>=0.3 in /usr/local/lib/python3.10/dist-packages (from arviz->pytransit) (0.8.0)\n",
            "Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas->pytransit) (2024.2)\n",
            "Requirement already satisfied: tzdata>=2022.1 in /usr/local/lib/python3.10/dist-packages (from pandas->pytransit) (2024.1)\n",
            "Requirement already satisfied: pyerfa>=2.0.1.1 in /usr/local/lib/python3.10/dist-packages (from astropy->pytransit) (2.0.1.4)\n",
            "Requirement already satisfied: astropy-iers-data>=0.2024.7.29.0.32.7 in /usr/local/lib/python3.10/dist-packages (from astropy->pytransit) (0.2024.9.16.0.32.21)\n",
            "Requirement already satisfied: PyYAML>=3.13 in /usr/local/lib/python3.10/dist-packages (from astropy->pytransit) (6.0.2)\n",
            "Requirement already satisfied: wrapt<2,>=1.10 in /usr/local/lib/python3.10/dist-packages (from deprecated->pytransit) (1.16.0)\n",
            "Requirement already satisfied: traitlets in /usr/local/lib/python3.10/dist-packages (from ldtk->pytransit) (5.7.1)\n",
            "Requirement already satisfied: tenacity in /usr/local/lib/python3.10/dist-packages (from ldtk->pytransit) (9.0.0)\n",
            "Requirement already satisfied: llvmlite<0.44,>=0.43.0dev0 in /usr/local/lib/python3.10/dist-packages (from numba->pytransit) (0.43.0)\n",
            "Requirement already satisfied: platformdirs>=2.2.0 in /usr/local/lib/python3.10/dist-packages (from pyopencl->pytransit) (4.3.6)\n",
            "Collecting pytools>=2024.1.5 (from pyopencl->pytransit)\n",
            "  Downloading pytools-2024.1.14-py3-none-any.whl.metadata (3.0 kB)\n",
            "Requirement already satisfied: multipledispatch in /usr/local/lib/python3.10/dist-packages (from pyrr->pytransit) (1.0.0)\n",
            "Requirement already satisfied: cython>=0.29.21 in /usr/local/lib/python3.10/dist-packages (from tables->pytransit) (3.0.11)\n",
            "Requirement already satisfied: numexpr>=2.6.2 in /usr/local/lib/python3.10/dist-packages (from tables->pytransit) (2.10.1)\n",
            "Requirement already satisfied: blosc2~=2.0.0 in /usr/local/lib/python3.10/dist-packages (from tables->pytransit) (2.0.0)\n",
            "Requirement already satisfied: py-cpuinfo in /usr/local/lib/python3.10/dist-packages (from tables->pytransit) (9.0.0)\n",
            "Requirement already satisfied: msgpack in /usr/local/lib/python3.10/dist-packages (from blosc2~=2.0.0->tables->pytransit) (1.0.8)\n",
            "Requirement already satisfied: h5py in /usr/local/lib/python3.10/dist-packages (from h5netcdf>=1.0.2->arviz->pytransit) (3.11.0)\n",
            "Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.7->matplotlib->pytransit) (1.16.0)\n",
            "Downloading PyTransit-2.6.8-py3-none-any.whl (10.9 MB)\n",
            "\u001b[2K   \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m10.9/10.9 MB\u001b[0m \u001b[31m18.4 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hDownloading celerite-0.4.3-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (307 kB)\n",
            "\u001b[2K   \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m308.0/308.0 kB\u001b[0m \u001b[31m11.4 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hDownloading emcee-3.1.6-py2.py3-none-any.whl (47 kB)\n",
            "\u001b[2K   \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m47.4/47.4 kB\u001b[0m \u001b[31m1.2 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hDownloading corner-2.2.2-py3-none-any.whl (15 kB)\n",
            "Downloading Deprecated-1.2.14-py2.py3-none-any.whl (9.6 kB)\n",
            "Downloading ldtk-1.8.4-py3-none-any.whl (38 kB)\n",
            "Downloading MeepMeep-0.6.3-py3-none-any.whl (33 kB)\n",
            "Downloading pyopencl-2024.2.7-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (698 kB)\n",
            "\u001b[2K   \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m698.1/698.1 kB\u001b[0m \u001b[31m19.4 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hDownloading pyrr-0.10.3-py3-none-any.whl (46 kB)\n",
            "\u001b[2K   \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m46.8/46.8 kB\u001b[0m \u001b[31m946.7 kB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hDownloading semantic_version-2.10.0-py2.py3-none-any.whl (15 kB)\n",
            "Downloading uncertainties-3.2.2-py3-none-any.whl (58 kB)\n",
            "\u001b[2K   \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m58.3/58.3 kB\u001b[0m \u001b[31m3.0 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hDownloading pytools-2024.1.14-py3-none-any.whl (89 kB)\n",
            "\u001b[2K   \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m89.9/89.9 kB\u001b[0m \u001b[31m2.5 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hInstalling collected packages: uncertainties, semantic-version, pytools, pyrr, emcee, deprecated, celerite, pyopencl, meepmeep, ldtk, corner, pytransit\n",
            "Successfully installed celerite-0.4.3 corner-2.2.2 deprecated-1.2.14 emcee-3.1.6 ldtk-1.8.4 meepmeep-0.6.3 pyopencl-2024.2.7 pyrr-0.10.3 pytools-2024.1.14 pytransit-2.6.8 semantic-version-2.10.0 uncertainties-3.2.2\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.10/dist-packages/pytools/persistent_dict.py:63: RecommendedHashNotFoundWarning: Unable to import recommended hash 'siphash24.siphash13', falling back to 'hashlib.sha256'. Run 'python3 -m pip install siphash24' to install the recommended hash.\n",
            "  warn(\"Unable to import recommended hash 'siphash24.siphash13', \"\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from astropy import units as u\n",
        "from astropy import constants as const\n",
        "\n",
        "!pip install pytransit celerite emcee corner\n",
        "from pytransit import QuadraticModel\n",
        "# from pytransit import UniformModel\n",
        "from scipy.optimize import minimize"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# TTV calculations"
      ],
      "metadata": {
        "id": "ce4IMPmaijKK"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def barycentre(Mp, Mm, dpm, moon_period):\n",
        "  dp = dpm/(1+(Mp/Mm)) # DIST BARYCENTRE TO PLANET\n",
        "  dm = dpm/(1+(Mm/Mp)) # DIST BARYCENTRE TO MOON\n",
        "  print(f'Planet-barycentre distance = {dp:.3f}\\nMoon-barycentre distance = {dm:.3f}')\n",
        "  return(dp, dm)\n",
        "\n",
        "# PARAMETERS\n",
        "Mp = 1*const.M_earth\n",
        "Mm = 0.0123*const.M_earth\n",
        "dpm = 0.3844e9*u.m\n",
        "moon_period = 27.322*u.day\n",
        "P = 1*u.year\n",
        "P = P.to(u.day)\n",
        "\n",
        "dp, dm = barycentre(Mp, Mm, dpm, moon_period)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "YRrTetmCkCvM",
        "outputId": "ef2a5e02-55c7-441a-fba9-103e3120474b"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Planet-barycentre distance = 4670670.750 m\n",
            "Moon-barycentre distance = 379729329.250 m\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "N_trials = 10 # NUMBER OF ANGLES\n",
        "N_obs = 10 # NUMBER OF OBSERVATIONS/TRANSITS (SHOULD BE WITHIN 5 YEARS AT LEAST)"
      ],
      "metadata": {
        "id": "Xe-pLTr_kWXn"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "t = np.linspace(0, 2*P, 1000)\n",
        "\n",
        "# BARYCENTRE\n",
        "ab = 1*const.au # SEPARATION\n",
        "Pb = P.to(u.s)  # PERIOD\n",
        "wb = (2*np.pi/Pb) *u.rad # ANGULAR FREQ\n",
        "alphab = wb*t\n",
        "xb = ab*np.sin(alphab)\n",
        "\n",
        "angle = []\n",
        "def position(N_trials):\n",
        "  for i in range(N_trials):\n",
        "    angle.append(np.random.uniform(0, 2*np.pi)*u.rad)\n",
        "    #print(f\"Start angle #{i+1} = {angle[i]:.3f}\") # IF WE WANT TO PRINT THE LIST OF ALL STARTING ANGLES OF 'N' REALISATIONS\n",
        "\n",
        "    # PLANET\n",
        "    wp = 2*np.pi/moon_period *u.rad\n",
        "    alphap = wp*t + angle[i]\n",
        "    xp = dp * np.sin(alphap) + xb\n",
        "\n",
        "    # MOON\n",
        "    wm = 2*np.pi/moon_period *u.rad\n",
        "    alpham = wm*t + np.pi*u.rad + angle[i]\n",
        "    xm = dm*np.sin(alpham) + xb\n",
        "  return angle, xb, xp, xm, wb, wp, wm, ab, alphap, alpham\n",
        "\n",
        "\n",
        "angle, xb, xp, xm, wb, wp, wm, ab, alphap, alpham = position(N_trials)\n",
        "print('Done')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "5nsaVyZbiqOA",
        "outputId": "5d12feb7-b279-4009-90ee-b305ff9e0583"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Done\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "t = np.arange(0,N_obs*P.value,P.value)*P.unit # TIME\n",
        "\n",
        "# BARYCENTRE\n",
        "ab = 1*const.au # SEPARATION\n",
        "Pb = P.to(u.s)  # PERIOD\n",
        "wb = (2*np.pi/Pb) *u.rad # ANGULAR FREQ\n",
        "alphab = wb*t\n",
        "xb = ab*np.sin(alphab)\n",
        "\n",
        "ttvp = []\n",
        "det = 0\n",
        "\n",
        "for i in range(N_trials):\n",
        "  angle = np.random.uniform(0, 2*np.pi)*u.rad\n",
        "\n",
        "  # PLANET\n",
        "  wp = 2*np.pi/moon_period *u.rad\n",
        "  alphap = wp*t + angle\n",
        "  xp = dp * np.sin(alphap) + xb\n",
        "  diffp = xb - xp\n",
        "\n",
        "  vb = (2*np.pi*ab)/P.to(u.s) # BARYCENTRE VELOCITY\n",
        "  #print(f'Velocity of system around star = {vb:.3f}')\n",
        "  ttvp.append(diffp/vb)\n",
        "  if N_trials <=10:\n",
        "    print(ttvp[i])\n",
        "else:\n",
        "  print('Too many to print!')\n",
        "ttvp = ttvp[i].to(u.day)\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "S94uWWAHuoSV",
        "outputId": "790253ef-03f3-4891-c1be-0b406407b80f"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[  95.62783773 -156.20803677  115.86007199   -0.6533191  -114.97555226\n",
            "  156.31713588  -96.66006525  -25.45040472  131.11701771 -152.06711983] s\n",
            "[-122.62587232  154.94947048  -87.15808183  -36.94734382  137.18058179\n",
            " -148.77955664   64.250011     61.7923516  -147.90982483  138.46072326] s\n",
            "[  32.99742725   90.49714804 -155.52026967  120.05960328   -7.02677377\n",
            " -110.54615104  156.69368525 -101.59927146  -19.13983082  127.51242481] s\n",
            "[  44.98616102 -141.02041564  145.93943761  -56.5649791   -69.35689141\n",
            "  150.46631921 -134.35738598   31.43829532   91.7935109  -155.71626557] s\n",
            "[  11.28864667 -122.76112586  154.91595099  -86.97744671  -37.15838372\n",
            "  137.28567069 -148.71079516   64.05182687   61.99190887 -147.98181855] s\n",
            "[ 147.8005812   -61.49080813  -64.54902292  148.8828413  -137.02140554\n",
            "   36.62855263   87.43051258 -154.99951953  122.42120229  -10.7448603 ] s\n",
            "[  92.89494158   30.10198974 -133.64961094  150.84437791  -70.57651492\n",
            "  -55.29180632  145.43533004 -141.61108411   46.2899665    78.93973596] s\n",
            "[ -92.88123919  -30.11868057  133.65850603 -150.83973002   70.56132712\n",
            "   55.30772099 -145.44168888  141.6037786   -46.27371682  -78.95443066] s\n",
            "[ 154.03491883 -125.89646148   16.41456655  103.67300589 -156.77602146\n",
            "  108.58389048    9.76578753 -121.80565573  155.14521228  -88.24331055] s\n",
            "[ -98.82414103  156.51236752 -113.0757979    -3.42060199  117.70690383\n",
            " -155.94115407   93.41967708   29.46153983 -133.30725012  151.02131022] s\n",
            "Too many to print!\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Light curve"
      ],
      "metadata": {
        "id": "nS8pCgcWiPK3"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# SWAPPED THE t0 FOR ttvp FOR THE TIME OF INFERIOR CONJUNCTION IN DAYS TO RUN IN A LOOP TO GET MULTIPLE LIGHT CURVES\n",
        "per = 365.25                         #orbital period in days\n",
        "rp_rs = const.R_earth/const.R_sun        #planet radius / stellar radius ratio\n",
        "ars =  const.au/const.R_sun              #semi-major axis / stellar radius ratio\n",
        "inc =  (90*u.deg).to(u.rad).value    #orbital inclination (in radians)\n",
        "ecc = 0.                             #eccentricity\n",
        "w = (90*u.deg).to(u.rad).value       #longitude of periastron (in radians)\n",
        "gamma = [0.231, 0.226]                 #limb darkening coefficients [u1, u2]\n",
        "\n",
        "texp = 100*u.s.to(u.day)             # the exposure time : we convert seconds to days\n",
        "t = np.arange(-0.5, 0.5, texp)       # new timegrid #note we use np.arange rather than np.linspace (DAYS) MANDEL AND........\n",
        "tm = QuadraticModel()                # a model that uses two limb-darkening coefficients\n",
        "tm.set_data(t)\n",
        "\n",
        "lc_Earth = []\n",
        "for j in range(len(ttvp)):\n",
        "  # ttvp[j] is a scalar Quantity, so we directly use it as t0_val\n",
        "  t0_val = ttvp[j]\n",
        "  lc_Earth.append(tm.evaluate(k=rp_rs, ldc=gamma, t0=t0_val, p=per, a=ars, i=inc, e=ecc, w=w))\n",
        "  #plt.plot(t, lc_Earth[j], '-o') # Plot the last added light curve\n",
        "print('Done')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "wLmFYK49iNY1",
        "outputId": "38ccaec2-00ee-47a1-9bf7-184d40c41a34"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Done\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "def planck(wl, T):\n",
        "  a = np.float64(1.191042768e8)*u.um**5 *u.W/ u.m**2 /u.sr/u.um\n",
        "  b = np.float64(14387.7516)*1*u.um * 1*u.K\n",
        "  try:\n",
        "    x = b/(wl*T)\n",
        "    bb = a/wl**5 / (np.exp(x) - 1.0)\n",
        "  except ArithmeticError:\n",
        "    bb = np.zeros(np.size(wl))\n",
        "  return bb\n",
        "\n",
        "s_lum, s_rad = 1*const.L_sun, 1*const.R_sun              # STELLAR VALUES\n",
        "T_s = (s_lum/(4*np.pi*s_rad**2*const.sigma_sb))**0.25    # TEMP-LUMINOSITY RELATIONSHIP\n",
        "\n",
        "wl = np.linspace(0.6, 5.3, 1000)*u.um                    # JWST NIRSPEC PRISM WAVELENGTH RANGE (OR HUBBLE 1.1-1.7)\n",
        "wl_full = np.linspace(0,6, 1000)*u.um\n",
        "BB_flux = np.pi*u.sr*planck(wl,T_s)                      # FLUX DENSITY - W/M^2/MICRONS\n",
        "BB_flux_full = np.pi*u.sr*planck(wl_full,T_s)\n",
        "\n",
        "####\n",
        "\n",
        "R_s = 1*u.Rsun.to(u.m) # host star radius in Rsun (code converts this to m)\n",
        "d = 10*u.pc.to(u.m) # distance to star in pc (code converts this to m)\n",
        "\n",
        "Flux = BB_flux*(R_s/d)**2\n",
        "\n",
        "D = 6.5*u.m; Atel = np.pi*(D/2)**2 # APERTURE OF SPECIFIC TELESCOPE\n",
        "trans = 0.5; QE = 0.8 # QUANTUM EFFICIENCY\n",
        "Power_per_micron = Flux*Atel*trans*QE # WATTS / MICRON\n",
        "Power = Power_per_micron*np.gradient(wl) # WATTS\n",
        "\n",
        "Photons_per_second = Power/(const.h*const.c/wl.to(u.m)) # CONVERT WATTS TO PHOTONS\n",
        "Photons_per_second = np.sum(Photons_per_second).value*1/u.s # TOTAL ELECTRONS PER SECOND\n",
        "\n",
        "exposure_time = 100*u.s # GENRALISED EXPOSURE\n",
        "Electrons = (Photons_per_second*exposure_time).value\n",
        "std = Electrons**0.5\n",
        "\n",
        "Noise = Electrons**0.5 # ASSUMING A NOISE LIMITED INSTRUMENT\n",
        "Noise = np.random.normal(0, Noise, len(lc_Earth[0]))\n",
        "\n",
        "'''\n",
        "for i in range(len(lc_Earth)):\n",
        "  lc = lc_Earth[0]*Electrons\n",
        "  lc = lc + Noise\n",
        "  plt.plot(t, lc, '-o', label = f'transit {i+1}')\n",
        "  plt.grid()\n",
        "  plt.title('Noise light curves for each ')\n",
        "  plt.ylabel('Relative signal')\n",
        "  plt.xlabel('Time (days)')\n",
        "\n",
        "  if N_obs <=10:\n",
        "    plt.legend()\n",
        "'''"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 53
        },
        "id": "b7EIeasS7nKW",
        "outputId": "eff21175-316c-43c4-93d5-37c9b55bc7f4"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "\"\\nfor i in range(len(lc_Earth)):\\n  lc = lc_Earth[0]*Electrons\\n  lc = lc + Noise\\n  plt.plot(t, lc, '-o', label = f'transit {i+1}')\\n  plt.grid()\\n  plt.title('Noise light curves for each ')\\n  plt.ylabel('Relative signal')\\n  plt.xlabel('Time (days)')\\n\\n  if N_obs <=10:\\n    plt.legend()\\n\""
            ],
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "string"
            }
          },
          "metadata": {},
          "execution_count": 11
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "def chi_squared(X, lc, Noise):\n",
        "    rp_rs = X[0]\n",
        "    S = X[1]\n",
        "    t_0 = X[2]\n",
        "    model =tm.evaluate(k=rp_rs, ldc=gamma, t0=t_0, p=per, a=ars, i=inc, e=ecc, w=w) *S\n",
        "    return np.sum(((model-lc))**2/Noise**2)\n",
        "\n",
        "t0 = ttvp\n",
        "data_lc = lc\n",
        "S= Electrons\n",
        "\n",
        "fit_init = [((lc.max()-lc.min())/ lc.max())**0.5, np.mean(lc[0:40]), t0[0].value+np.random.normal(0,0.001)] # these are initial guesses - we'll cheat a bit by putting in the known values but you can try starting with different initial values too\n",
        "# Access the first element of t0.value to ensure it's a scalar\n",
        "bounds =((rp_rs*0.99,rp_rs*1.1), (S*0.99,S*1.1), (np.min(ttvp.value),np.max(ttvp.value))) # the bounds over which the algorithm will vary the parameters\n",
        "# Use std instead of Noise\n",
        "fit  = minimize(chi_squared, fit_init, args=(data_lc, std), method='Nelder-Mead', jac=None, hess=None, hessp=None, bounds=bounds, constraints=(), tol=None, callback=None, options=None)\n",
        "final_fit = [fit['x'][0], fit['x'][1], fit['x'][2]]\n",
        "if N_obs <= 10:\n",
        "  print ('Fitted Rp/Rs', final_fit[0], 'Actual Rp/Rs', rp_rs)\n",
        "  print ('Fitted S', final_fit[1], 'Actual S', S)\n",
        "  print ('Fitted t0', final_fit[2], 'Actual t0', t0)\n",
        "\n",
        "model_fit = tm.evaluate(k=fit['x'][0], ldc=gamma, t0=fit['x'][2], p=per, a=ars, i=inc, e=ecc, w=w) *fit['x'][1]\n",
        "'''\n",
        "plt.figure('lc')\n",
        "plt.plot(t,data_lc, 'o')\n",
        "plt.plot(t,model_fit, '-')\n",
        "plt.show();\n",
        "'''\n",
        "print('Done')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 211
        },
        "id": "ZJC48-4y77F-",
        "outputId": "c8d8b84a-c5e6-495e-a2aa-3ed77e0089b3"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "error",
          "ename": "NameError",
          "evalue": "name 'lc' is not defined",
          "traceback": [
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
            "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
            "\u001b[0;32m<ipython-input-12-e155fb1376b4>\u001b[0m in \u001b[0;36m<cell line: 9>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0mt0\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mttvp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m \u001b[0mdata_lc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlc\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     10\u001b[0m \u001b[0mS\u001b[0m\u001b[0;34m=\u001b[0m \u001b[0mElectrons\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;31mNameError\u001b[0m: name 'lc' is not defined"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "t0_list = []\n",
        "\n",
        "for i in range(N_obs):\n",
        "  Noise = Electrons**0.5 # ASSUMING A NOISE LIMITED INSTRUMENT\n",
        "  Noise = (np.random.normal(0, Noise, len(lc_Earth[0])))\n",
        "\n",
        "  lc = lc_Earth[0]*Electrons\n",
        "  lc = lc + Noise\n",
        "\n",
        "  if N_obs <=10:\n",
        "    plt.plot(t, lc, '-o')\n",
        "    plt.grid()\n",
        "    plt.title('Noise light curves for each ttv')\n",
        "    plt.ylabel('Relative signal')\n",
        "    plt.xlabel('Time (days)')\n",
        "\n",
        "  t0 = ttvp\n",
        "  data_lc = lc\n",
        "  S = Electrons\n",
        "\n",
        "  fit_init = [((lc.max()-lc.min())/ lc.max())**0.5, np.mean(lc[0:40]), t0[0].value+np.random.normal(0,0.001)] # these are initial guesses - we'll cheat a bit by putting in the known values but you can try starting with different initial values too\n",
        "  # Access the first element of t0.value to ensure it's a scalar\n",
        "  bounds =((rp_rs*0.99,rp_rs*1.1), (S*0.99,S*1.1), (-0.01,0.01)) # the bounds over which the algorithm will vary the parameters\n",
        "  # Use std instead of Noise\n",
        "  fit  = minimize(chi_squared, fit_init, args=(data_lc, std), method='Nelder-Mead', jac=None, hess=None, hessp=None, bounds=bounds, constraints=(), tol=None, callback=None, options=None)\n",
        "  final_fit = [fit['x'][0], fit['x'][1], fit['x'][2]]\n",
        "  #print ('Fitted Rp/Rs', final_fit[0], 'Actual Rp/Rs', rp_rs_Earth)\n",
        "  #print ('Fitted S', final_fit[1], 'Actual S', S)\n",
        "  #print ('Fitted t0', final_fit[2], 'Actual t0', t0)\n",
        "\n",
        "  model_fit = tm.evaluate(k=fit['x'][0], ldc=gamma, t0=fit['x'][2], p=per, a=ars, i=inc, e=ecc, w=w) * fit['x'][1]\n",
        "\n",
        "  #plt.figure('lc')\n",
        "  #plt.plot(t,data_lc, 'o')\n",
        "  #plt.plot(t,model_fit, '-')\n",
        "\n",
        "  t0_list.append(final_fit[2])\n",
        "\n",
        "#plt.show()"
      ],
      "metadata": {
        "id": "-YBuKINOFcKb"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Probability\n",
        "\n",
        "This calculates the SNR of all light curves and we should expect that for Earth, each transit or observation there should be a detection of a transit timing variation due to size of the Moon and separation."
      ],
      "metadata": {
        "id": "-WzUFPNKbGUL"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "t0_list_sec = [t*u.day.to(u.s) for t in t0_list]\n",
        "stddev = np.std(t0_list_sec)*u.s\n",
        "\n",
        "det = 0\n",
        "\n",
        "for i in range(len(ttvp)):\n",
        "  SNR = np.abs(ttvp[i]/stddev)\n",
        "  #print(SNR) # IF ANY ARE GREATER THAN 5*SIGMA THEN SUCCESSFUL DETECTION\n",
        "\n",
        "  if SNR >= 5:\n",
        "    if N_obs <=10:\n",
        "      print('Detection')\n",
        "    det = det+1\n",
        "  else:\n",
        "    print(f'Obs {i+1} = No detection')\n",
        "\n",
        "prob = 100*det/N_obs\n",
        "print(f'Probability of detection out of {N_obs} transits = {prob:.1f}%')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "JX9GzPEx7Bhu",
        "outputId": "3c8e4181-8353-4694-a610-dbbd92cef92c"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Obs 10 = No detection\n",
            "Obs 29 = No detection\n",
            "Obs 48 = No detection\n",
            "Obs 67 = No detection\n",
            "Obs 86 = No detection\n",
            "Obs 215 = No detection\n",
            "Obs 234 = No detection\n",
            "Obs 253 = No detection\n",
            "Obs 272 = No detection\n",
            "Obs 291 = No detection\n",
            "Obs 310 = No detection\n",
            "Obs 329 = No detection\n",
            "Obs 348 = No detection\n",
            "Obs 367 = No detection\n",
            "Obs 386 = No detection\n",
            "Obs 405 = No detection\n",
            "Obs 424 = No detection\n",
            "Obs 443 = No detection\n",
            "Obs 462 = No detection\n",
            "Obs 572 = No detection\n",
            "Obs 591 = No detection\n",
            "Obs 610 = No detection\n",
            "Obs 629 = No detection\n",
            "Obs 648 = No detection\n",
            "Obs 667 = No detection\n",
            "Obs 686 = No detection\n",
            "Obs 705 = No detection\n",
            "Obs 724 = No detection\n",
            "Obs 743 = No detection\n",
            "Obs 762 = No detection\n",
            "Obs 781 = No detection\n",
            "Obs 800 = No detection\n",
            "Obs 819 = No detection\n",
            "Obs 929 = No detection\n",
            "Obs 948 = No detection\n",
            "Obs 967 = No detection\n",
            "Obs 986 = No detection\n",
            "Probability of detection out of 1000 transits = 96%\n"
          ]
        }
      ]
    }
  ]
}