fitchisqperiod.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# coding: utf-8\n",
    "import numpy as np\n",
    "from matplotlib.pyplot import subplots, show\n",
    "from kapteyn import kmpfit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.load(\"chisqs.B0329+54.npz\")\n",
    "period = data['period']\n",
    "chisq = data['chisq']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss(p,x):\n",
    "    return p[0]*np.exp(-(x-p[2])**2/p[1]**2)+p[3]\n",
    "\n",
    "def absexp(p,x):\n",
    "    return p[0]*np.exp(-np.abs(x-p[2])/p[1])+p[3]\n",
    "\n",
    "def lorentz(p,x):\n",
    "    return p[0]*p[1]**2/((x-p[2])**2 + p[1]**2)+p[3]\n",
    "\n",
    "model = absexp\n",
    "\n",
    "fitobj=kmpfit.simplefit(model,[6e-10, 1e-6, period.mean()-1e-5, 0.4e-10], period, chisq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"\\nFit status kmpfit:\")\n",
    "print(\"====================\")\n",
    "print(\"Best-fit parameters:        \", fitobj.params)\n",
    "print(\"Asymptotic error:           \", fitobj.xerror)\n",
    "print(\"Error assuming red.chi^2=1: \", fitobj.stderr)\n",
    "print(\"Chi^2 min:                  \", fitobj.chi2_min)\n",
    "print(\"Reduced Chi^2:              \", fitobj.rchi2_min)\n",
    "print(\"Iterations:                 \", fitobj.niter)\n",
    "print(\"Number of free pars.:       \", fitobj.nfree)\n",
    "print(\"Degrees of freedom:         \", fitobj.dof)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.interpolate import interp1d\n",
    "from scipy.optimize import differential_evolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interp = interp1d(period,chisq,kind='cubic')\n",
    "bestperiod = differential_evolution(lambda x: -interp(x), np.array([period.min(),period.max()])[None,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bestperiod.x - period.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = subplots()\n",
    "ax.plot(period,chisq)\n",
    "ax.plot(period, model(fitobj.params, period))\n",
    "ax.axvline(period.mean())\n",
    "ax.axvline(fitobj.params[2])\n",
    "ax.axvline(bestperiod.x)\n",
    "ax.set_xlim(fitobj.params[2]-1e-6, fitobj.params[2]+1e-6)\n",
    "#ax.set_ylim(1.2e-7, 1.35e-7)\n",
    "#ax.set_xlim(56.5,57)\n",
    "show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(period.mean()-fitobj.params[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from observation import Observation\n",
    "obs = Observation(\"data/obs-10-04-2018/B2111+46_10-04-2018.fits.gz\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "import numpy as np\n",
    "\n",
    "hdulist = fits.open(\"data/obs-10-04-2018/B2310+42-04-2018.fits.gz\")\n",
    "col_period = fits.Column(name='period', array=data['period'], format='D')\n",
    "col_DM = fits.Column(name='DM', array=data['DM'], format='D')\n",
    "col_chisq = fits.Column(name='chisq', array=data['chisq'], format='D')\n",
    "t = fits.BinTableHDU.from_columns([col_period, col_DM, col_chisq])\n",
    "t.header['fitpars'] = str(fitobj.params)\n",
    "t.header['bestp'] = fitobj.params[2]\n",
    "t.header['fitfunc'] = model.__name__\n",
    "t.header['fitpars'] = \"-\"\n",
    "t.header['bestp'] = bestperiod.x[0]\n",
    "t.header['fitfunc'] = 'maximum'\n",
    "hdulist.append(t)\n",
    "#hdulist.writeto(\"data/obs-10-04-2018/B2310+42-04-2018-withP.fits\")"
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}