timedelay.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run proof_setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, cosf, sinf, e, aoc, sinw, cosw, sini = sm.symbols(\"n, cosf, sinf, e, aoc, sinw, cosw, sini\", real=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "aoc*sini*(-e**2 + 1)*(cosf*sinw + cosw*sinf)/(cosf*e + 1)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# e * cos(f)\n",
    "ecosf = e * cosf\n",
    "\n",
    "# 1 - e^2\n",
    "ome2 = 1 - e ** 2\n",
    "\n",
    "# Partials\n",
    "dfdM = (1 + ecosf) ** 2 / (ome2 * sm.sqrt(ome2))\n",
    "dfde = (2 + ecosf) * sinf / ome2\n",
    "\n",
    "r = aoc * (1 - e**2) / (1 + e * cosf)\n",
    "dt = sini * r * (sinw * cosf + cosw * sinf)\n",
    "dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "aoc*n*sini*(cosf**2*cosw*e + cosf*cosw + cosw*e*sinf**2 - sinf*sinw)/sqrt(-e**2 + 1)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ddtdM = sm.simplify(sm.diff(dt, sinf)) * cosf * dfdM\n",
    "ddtdM -= sm.simplify(sm.diff(dt, cosf)) * sinf * dfdM\n",
    "ddtdte = sm.simplify(ddtdM * n)\n",
    "ddtdte"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "dtede = sm.simplify(sm.diff(dt, e) / (1 + ddtdte))\n",
    "dtedaoc = sm.simplify(sm.diff(dt, aoc) / (1 + ddtdte))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "aoc*sini*(cosf*sinw + cosw*sinf)*(cosf*(e**2 - 1) - 2*e*(cosf*e + 1))/(cosf*e + 1)**2"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sm.simplify(sm.diff(dt, e))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-sini*(e**2 - 1)*(cosf*sinw + cosw*sinf)/(cosf*e + 1)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sm.simplify(sm.diff(dt, aoc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-aoc*cosf*sini*(e**2 - 1)/(cosf*e + 1)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sm.simplify(sm.diff(dt, sinw))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-aoc*sinf*sini*(e**2 - 1)/(cosf*e + 1)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sm.simplify(sm.diff(dt, cosw))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-aoc*(e**2 - 1)*(cosf*sinw + cosw*sinf)/(cosf*e + 1)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sm.simplify(sm.diff(dt, sini))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import astropy.units as u\n",
    "import astropy.constants as c\n",
    "\n",
    "import theano\n",
    "import theano.tensor as tt\n",
    "import exoplanet as xo\n",
    "\n",
    "def compute_time_delay(n, tperi, aoc, sini, e, sinw, cosw, tobs):\n",
    "    \n",
    "    t = tt.dscalar()\n",
    "    M = n * (t - tperi)\n",
    "    sinf, cosf = xo.theano_ops.kepler.KeplerOp()(M, e)\n",
    "    tobs_mod = t + aoc*sini*(1-e**2)*(cosf*sinw + cosw*sinf)/(1 + e*cosf)\n",
    "    dtobsdt = 1 + n*aoc*sini*(cosf**2*cosw*e + cosf*cosw + cosw*e*sinf**2 - sinf*sinw)/np.sqrt(1-e**2)\n",
    "    \n",
    "    func = theano.function([t], [tobs_mod, dtobsdt])\n",
    "    \n",
    "    t0 = tobs - aoc\n",
    "    f, g = func(t0)\n",
    "    delta = tobs - f\n",
    "    while np.abs(delta) > 1e-12:\n",
    "        t0 += delta / g\n",
    "        f, g = func(t0)\n",
    "        delta = tobs - f\n",
    "    return t0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.00026858858915434537"
      ]
     },
     "execution_count": 150,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(10 * u.R_sun / c.c.to(u.R_sun / u.day)).value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.004028815929189739"
      ]
     },
     "execution_count": 156,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "aoc = (150 * u.R_sun / c.c.to(u.R_sun / u.day)).value\n",
    "w = -np.pi\n",
    "period = 10\n",
    "n = 2 * np.pi / period\n",
    "e = 0.0\n",
    "\n",
    "opsw = 1 + np.sin(w)\n",
    "E0 = 2 * np.arctan2(\n",
    "    np.sqrt(1 - e) * np.cos(w),\n",
    "    np.sqrt(1 + e) * opsw,\n",
    ")\n",
    "M0 = E0 - e * np.sin(E0)\n",
    "\n",
    "t0 = 0.0\n",
    "tref = t0 - M0 / n\n",
    "\n",
    "t = t0\n",
    "compute_time_delay(n, tref, aoc, 1.0, e, np.sin(w), np.cos(w), t) - t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.728056416212937"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "0.003977816955703428 * 24 * 60"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}