{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "from pycheops.models import FactorModel, EclipseModel \n",
    "from lmfit import report_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "T_0 = 0.0\n",
    "P = 2.0\n",
    "D = 0.02\n",
    "W = 0.04\n",
    "S = 0.50\n",
    "F = 0.005  # Planet-star flux ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Parameters([('s_c',\n",
       "             <Parameter 's_c', 0.9950248756218907, bounds=[-inf:inf], expr='1/(1+F)'>),\n",
       "            ('p_c',\n",
       "             <Parameter 'p_c', 0.004975124378109453, bounds=[-inf:inf], expr='1/(1+1/F)'>),\n",
       "            ('T_0', <Parameter 'T_0', 0.0, bounds=[-inf:inf]>),\n",
       "            ('P', <Parameter 'P', value=2.0 (fixed), bounds=[1e-15:inf]>),\n",
       "            ('D', <Parameter 'D', value=0.02 (fixed), bounds=[0:1]>),\n",
       "            ('W', <Parameter 'W', 0.04, bounds=[0:0.3]>),\n",
       "            ('S', <Parameter 'S', 0.5, bounds=[0:1]>),\n",
       "            ('f_c', <Parameter 'f_c', value=0 (fixed), bounds=[-1:1]>),\n",
       "            ('f_s', <Parameter 'f_s', value=0 (fixed), bounds=[-1:1]>),\n",
       "            ('a_c', <Parameter 'a_c', value=0.001 (fixed), bounds=[0:inf]>),\n",
       "            ('k',\n",
       "             <Parameter 'k', 0.1414213562373095, bounds=[0:1], expr='sqrt(D)'>),\n",
       "            ('aR',\n",
       "             <Parameter 'aR', 8.464330808241401, bounds=[1:inf], expr='2/(pi*W*sqrt((1-S)/k))'>),\n",
       "            ('rho',\n",
       "             <Parameter 'rho', 2.0342563549033525, bounds=[0:inf], expr='0.013418*aR**3/P**2'>),\n",
       "            ('b',\n",
       "             <Parameter 'b', 0.414091611332738, bounds=[-inf:1.3], expr='sqrt(((1-k)**2-S*(1+k)**2)/(1-S))'>),\n",
       "            ('F', <Parameter 'F', 0.005, bounds=[0:inf]>)])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = FactorModel(prefix='s_') + FactorModel(prefix='p_')*EclipseModel()\n",
    "pars = m.make_params(T_0=T_0, W=W, S=S)\n",
    "pars['P'].set(value=P, vary=False)\n",
    "pars['D'].set(value=D, vary=False)\n",
    "pars.add(name='F', value=F, min=0)\n",
    "pars['a_c'].set(value=0.001, vary=False)\n",
    "pars['s_c'].set(expr='1/(1+F)')  # Flux from star \n",
    "pars['p_c'].set(expr='1/(1+1/F)')  # Flux from planet\n",
    "pars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.9925, 1.0025)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 12500 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "t = np.arange(T_0 + (0.5-W)*P,T_0 + (0.5+W)*P, 0.001)\n",
    "yerr = 50e-6\n",
    "flux_true = m.eval(params=pars, t=t)\n",
    "flux = flux_true + normal(0.0, yerr, size=len(t))\n",
    "fig,ax=subplots()\n",
    "ax.scatter(t, flux,s=4)\n",
    "ax.plot(t, flux_true,'c')\n",
    "ax.set_xlabel('Time')\n",
    "ax.set_ylabel('Flux')\n",
    "ax.set_xlim(min(t),max(t))\n",
    "ax.set_ylim(1-1.5*F,1+0.5*F)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[Fit Statistics]]\n",
      "    # fitting method   = leastsq\n",
      "    # function evals   = 21\n",
      "    # data points      = 161\n",
      "    # variables        = 4\n",
      "    chi-square         = 150.570525\n",
      "    reduced chi-square = 0.95904793\n",
      "    Akaike info crit   = -2.78261832\n",
      "    Bayesian info crit = 9.54299914\n",
      "[[Variables]]\n",
      "    s_c:  0.99502521 +/- 5.8485e-06 (0.00%) == '1/(1+F)'\n",
      "    p_c:  0.00497479 +/- 5.8485e-06 (0.12%) == '1/(1+1/F)'\n",
      "    T_0: -5.6934e-05 +/- 2.5433e-05 (44.67%) (init = 0)\n",
      "    P:    2 (fixed)\n",
      "    D:    0.02 (fixed)\n",
      "    W:    0.03996047 +/- 8.0580e-05 (0.20%) (init = 0.04)\n",
      "    S:    0.50167037 +/- 0.00290831 (0.58%) (init = 0.5)\n",
      "    f_c:  0 (fixed)\n",
      "    f_s:  0 (fixed)\n",
      "    a_c:  0.001 (fixed)\n",
      "    k:    0.14142136 +/- 0.00000000 (0.00%) == 'sqrt(D)'\n",
      "    aR:   8.48689157 +/- 0.04138093 (0.49%) == '2/(pi*W*sqrt((1-S)/k))'\n",
      "    rho:  2.05056602 +/- 0.00000000 (0.00%) == '0.013418*aR**3/P**2'\n",
      "    b:    0.40948697 +/- 0.00808932 (1.98%) == 'sqrt(((1-k)**2-S*(1+k)**2)/(1-S))'\n",
      "    F:    0.00499966 +/- 5.9071e-06 (0.12%) (init = 0.005)\n",
      "[[Correlations]] (unreported correlations are < 0.100)\n",
      "    C(W, S) = -0.951\n",
      "    C(S, F) = -0.222\n",
      "    C(W, F) =  0.128\n"
     ]
    }
   ],
   "source": [
    "result = m.fit(flux, pars, t=t, weights=1/yerr)\n",
    "report_fit(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<Figure size 432x432 with 2 Axes>, GridSpec(2, 1, height_ratios=[1, 4]))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 33596 characters]\n",
      "text/plain": [
       "<Figure size 432x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "result.plot(xlabel='Time',ylabel='Flux')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pyvenv",
   "language": "python",
   "name": "pyvenv"
  },
  "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
}
