{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "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, ThermalPhaseModel\n",
    "from lmfit import report_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "T_0 = 0.5\n",
    "P = 2.0\n",
    "a_th = 0.012\n",
    "b_th = -0.003\n",
    "c_th = 0.0027\n",
    "D = 0.02\n",
    "W = 0.04\n",
    "S = 0.28"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Parameters([('_T_0', <Parameter '_T_0', 0.5, bounds=[-inf:inf], expr='T_0'>),\n",
       "            ('_P', <Parameter '_P', 2.0, bounds=[1e-15:inf], expr='P'>),\n",
       "            ('_a_th', <Parameter '_a_th', 0.012, bounds=[-inf:inf]>),\n",
       "            ('_b_th', <Parameter '_b_th', -0.003, bounds=[-inf:inf]>),\n",
       "            ('_c_th', <Parameter '_c_th', 0.0027, bounds=[0:inf]>),\n",
       "            ('T_0', <Parameter 'T_0', 0.5, bounds=[-inf:inf]>),\n",
       "            ('P', <Parameter 'P', value=2.0 (fixed), bounds=[1e-15:inf]>),\n",
       "            ('D', <Parameter 'D', 0.02, bounds=[0:1]>),\n",
       "            ('W', <Parameter 'W', 0.04, bounds=[0:0.3]>),\n",
       "            ('S', <Parameter 'S', 0.28, 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 (fixed), bounds=[0:inf]>),\n",
       "            ('c', <Parameter 'c', 1.0, bounds=[-inf:inf]>),\n",
       "            ('_A',\n",
       "             <Parameter '_A', 0.012369316876852983, bounds=[-inf:inf], expr='hypot(_a_th,_b_th)'>),\n",
       "            ('_Fmax',\n",
       "             <Parameter '_Fmax', 0.013384658438426494, bounds=[0:inf], expr='_c_th+(_a_th+_b_th+_A)/2'>),\n",
       "            ('_Fmin',\n",
       "             <Parameter '_Fmin', 0.0010153415615735111, bounds=[0:inf], expr='_Fmax - _A'>),\n",
       "            ('_ph_max',\n",
       "             <Parameter '_ph_max', -0.4610104348113153, bounds=[-inf:inf], expr='arctan2(_b_th,-_a_th)/(2*pi)'>),\n",
       "            ('k',\n",
       "             <Parameter 'k', 0.1414213562373095, bounds=[0:1], expr='sqrt(D)'>),\n",
       "            ('aR',\n",
       "             <Parameter 'aR', 7.053609006867835, bounds=[1:inf], expr='2/(pi*W*sqrt((1-S)/k))'>),\n",
       "            ('rho',\n",
       "             <Parameter 'rho', 1.1772316868653663, bounds=[0:inf], expr='0.013418*aR**3/P**2'>),\n",
       "            ('b',\n",
       "             <Parameter 'b', 0.8475196745964672, bounds=[-inf:1.3], expr='sqrt((1-k)**2-S*(1+k)**2)/(1-S)'>)])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = ThermalPhaseModel(prefix=\"_\")*EclipseModel()*FactorModel()\n",
    "pars = m.make_params(T_0=T_0, \n",
    "                     _a_th=a_th, _b_th=b_th, _c_th=c_th,\n",
    "                     D=D, W=W, S=S)\n",
    "pars['P'].set(value=P, vary=False)\n",
    "pars['_T_0'].set(expr='T_0')\n",
    "pars['_P'].set(expr='P')\n",
    "pars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.5, 2.499000000000002)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 22104 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "t = np.arange(T_0,T_0 + P, 0.001)\n",
    "yerr = 450e-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(0,1.5*pars['_F'].value)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Max flux from model = 0.0134, F = 0.0134\n",
      "Phase of maximum flux from model = 0.5390, ph_max = 0.5390\n",
      "Peak-to-trough amplitude from model = 0.013, A = 0.012\n"
     ]
    }
   ],
   "source": [
    "print('Max flux from model = {:0.4f}, F = {:0.4f}'.format(max(flux_true),pars['_Fmax'].value))\n",
    "print('Phase of maximum flux from model = {:0.4f}, ph_max = {:0.4f}'\n",
    "      .format((t[argmax(flux_true)]-T_0)/P,pars['_ph_max'].value % 1))\n",
    "print('Peak-to-trough amplitude from model = {:0.3f}, A = {:0.3f}'\n",
    "      .format(max(flux_true)-min(flux_true),pars['_A'].value))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[Fit Statistics]]\n",
      "    # fitting method   = leastsq\n",
      "    # function evals   = 88\n",
      "    # data points      = 2000\n",
      "    # variables        = 8\n",
      "    chi-square         = 2011.39108\n",
      "    reduced chi-square = 1.00973448\n",
      "    Akaike info crit   = 27.3587657\n",
      "    Bayesian info crit = 72.1659854\n",
      "[[Variables]]\n",
      "    _T_0:     0.49990573 == 'T_0'\n",
      "    _P:       2.00000000 == 'P'\n",
      "    _a_th:    0.01258269 (init = 0.012)\n",
      "    _b_th:   -0.00319169 (init = -0.003)\n",
      "    _c_th:    0.00288370 (init = 0.0027)\n",
      "    T_0:      0.49990573 (init = 0.5)\n",
      "    P:        2 (fixed)\n",
      "    D:        0.02720969 (init = 0.02)\n",
      "    W:        0.04062058 (init = 0.04)\n",
      "    S:        0.28994368 (init = 0.28)\n",
      "    f_c:      0 (fixed)\n",
      "    f_s:      0 (fixed)\n",
      "    a_c:      0 (fixed)\n",
      "    c:        0.95113742 (init = 1)\n",
      "    _A:       0.01298117 == 'hypot(_a_th,_b_th)'\n",
      "    _Fmax:    0.01406978 == '_c_th+(_a_th+_b_th+_A)/2'\n",
      "    _Fmin:    0.00108861 == '_Fmax - _A'\n",
      "    _ph_max: -0.46046311 == 'arctan2(_b_th,-_a_th)/(2*pi)'\n",
      "    k:        0.16495360 == 'sqrt(D)'\n",
      "    aR:       7.55385346 == '2/(pi*W*sqrt((1-S)/k))'\n",
      "    rho:      1.44588404 == '0.013418*aR**3/P**2'\n",
      "    b:        0.77626826 == 'sqrt((1-k)**2-S*(1+k)**2)/(1-S)'\n"
     ]
    }
   ],
   "source": [
    "result = m.fit(flux, pars, t=t, weights=1/yerr)\n",
    "report_fit(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<Figure size 432x432 with 2 Axes>, GridSpec(2, 1, height_ratios=[1, 4]))"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 53076 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',data_kws={'ms':1},ax_fit_kws={'xlim':[min(t),max(t)]})"
   ]
  },
  {
   "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
}
