Photometry+light-curve.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Module: do_aperture_photometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function uses the photutils packages to do aperture photometry of the aligned & reduced data.\n",
    "\n",
    "INPUT:\n",
    "\n",
    "     Filename of the the data image and mean FWHM.\n",
    "     \n",
    "OUTPUT:\n",
    "\n",
    "     Magnitude of sources of different apertures and their corresponding errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def do_aperture_photometry(filename, fwhm):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now select a mode for input images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "choice=input(\"Enter the mode. Please use the keywords\\n 'single' for single image, 'multiple'  for multiple images: \\n\\n\")\n",
    "\n",
    "if (choice=='single'):\n",
    "    print('single mode')\n",
    "elif (choice ==' multiple'):\n",
    "    print(' multiple image mode')\n",
    "else:\n",
    "    print('Enter valid choice!!!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now    \n",
    "       1. Identify the sources, their position. \n",
    "       2. Create the multiple apertures around the sources.\n",
    "       3. Get the magnitude and magnitude error for different apertures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data,header=fits.getdata(filename,header=True)\n",
    "data=2*data\n",
    "exposure=header['EXPOSURE']\n",
    "print('Exposure is',exposure)\n",
    "\n",
    "sigma_clip = SigmaClip(sigma=3., maxiters=10)\n",
    "bkg_estimator = MedianBackground()\n",
    "bkg = Background2D(data, (10,10), filter_size=(3, 3),sigma_clip=sigma_clip, bkg_estimator=bkg_estimator)\n",
    "back=bkg.background # this is the background we need for the background subtraction.\n",
    "back2=np.median(bkg.background)\n",
    "\n",
    "\n",
    "mask = data == 0\n",
    "unit = u.electron / u.s\n",
    "\n",
    "\n",
    "xdf_image = CCDData(data, unit=unit, meta=header, mask=mask)\n",
    "norm_image = ImageNormalize(vmin=1e-4, vmax=5e-2, stretch=LogStretch(), clip=False)\n",
    "xdf_image_clipped = np.clip(xdf_image, 1e-4, None)\n",
    "\n",
    "mean, median, std = sigma_clipped_stats(xdf_image.data, sigma=3.0, maxiters=20, mask=xdf_image.mask)\n",
    "\n",
    "print('Finding the sources')\n",
    "\n",
    "daofind = DAOStarFinder(fwhm=fwhm, threshold=10*std) # 3 sigma above the background.\n",
    "sources = daofind(data-back)\n",
    "\n",
    "#sources_findpeaks = find_peaks(xdf_image.data, mask=xdf_image.mask, threshold=30.*std, box_size=30, centroid_func=centroid_2dg)\n",
    "\n",
    "print('We have found:',len(sources),' sources')\n",
    "#print(sources)\n",
    "\n",
    "#print(sources['xcentroid'], sources['ycentroid'],sources['fwhm'])\n",
    "#positions=sources['xcentroid'], sources['ycentroid']\n",
    "positions=np.genfromtxt('co_ordinates_list.txt',unpack=True,usecols=(0,1))\n",
    "#print(positions)\n",
    "radii=[ fwhm,2*fwhm, 3*fwhm]\n",
    "\n",
    "#positions=(sources['xcentroid'], sources['ycentroid'])\n",
    "apertures = [CircularAperture(positions, r=r) for r in radii]\n",
    "\n",
    "an_ap = CircularAnnulus(positions, r_in=5*fwhm, r_out=7*fwhm)\n",
    "#apers = [apertures, annulus_apertures]\n",
    "\n",
    "\n",
    "#bkg_sigma=mad_std(data)\n",
    "effective_gain=exposure\n",
    "error=calc_total_error(data,back,effective_gain)\n",
    "\n",
    "\n",
    "#error=0.1*data\n",
    "phot_table = aperture_photometry(data-back, apertures,error=error)\n",
    "phot_table2=aperture_photometry(data-back,an_ap)\n",
    "\n",
    "\n",
    "bkg_mean = phot_table2['aperture_sum'] / an_ap.area()\n",
    "bkg_sum = bkg_mean * an_ap.area()\n",
    "\n",
    "\n",
    "final_sum0=phot_table['aperture_sum_0']-bkg_sum\n",
    "final_sum1=phot_table['aperture_sum_1']-bkg_sum\n",
    "final_sum2=phot_table['aperture_sum_2']-bkg_sum\n",
    "\n",
    "\n",
    "mag_back=-2.5*np.log10(bkg_mean/exposure)+25\n",
    "mag_0=-2.5*np.log10(final_sum0/exposure)+25\n",
    "mag_1=-2.5*np.log10(final_sum1/exposure)+25\n",
    "mag_2=-2.5*np.log10(final_sum2/exposure)+25\n",
    "\n",
    "flux_err_0=phot_table['aperture_sum_err_0']\n",
    "mag_err_0=1.09*flux_err_0/final_sum0\n",
    "\n",
    "flux_err_1=phot_table['aperture_sum_err_1']\n",
    "mag_err_1=1.09*flux_err_1/final_sum1\n",
    "\n",
    "flux_err_2=phot_table['aperture_sum_err_2']\n",
    "mag_err_2=1.09*flux_err_2/final_sum2\n",
    "\n",
    "fig=plt.figure()\n",
    "plt.imshow(data,cmap='gray',origin='lower',vmin=500,vmax=900)\n",
    "colors=['red','green','yellow']\n",
    "for i in range(len(apertures)):\n",
    "    apertures[i].plot(color=colors[i], alpha=0.7)\n",
    "\n",
    "an_ap.plot(color='green', alpha=0.7)\n",
    "plt.show()\n",
    "\n",
    "with open ('results.dat','w') as r:\n",
    "    for i in range (len(phot_table)):\n",
    "        #print(final_sum0[i],final_sum1[i],final_sum2[i],final_sum3[i],final_sum4[i],final_sum5[i],file=r)\n",
    "        print(mag_back[i],mag_0[i],mag_err_0[i],mag_1[i],mag_err_1[i],mag_2[i],mag_err_2[i],file=r)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Photometry is completed & now we are ready to plot the light curve...!!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Module: plot_light_curve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function is used to plot the the light curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_light_curve(filename,list):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#~~~~~~~~~~~~~~~~~~~~~~for plotting purposes~~~~~~~~#\n",
    "    mp.rcParams['font.family']='serif'\n",
    "    mp.rcParams['xtick.major.size']=10\n",
    "    mp.rcParams['xtick.major.width']=2\n",
    "    mp.rcParams['xtick.minor.size']=7\n",
    "    mp.rcParams['xtick.minor.width']=2\n",
    "    mp.rcParams['ytick.major.size']=10\n",
    "    mp.rcParams['ytick.major.width']=2\n",
    "    mp.rcParams['ytick.minor.size']=7\n",
    "    mp.rcParams['ytick.minor.width']=2\n",
    "    mp.rcParams['axes.linewidth']=1.5\n",
    "    mp.rcParams['xtick.labelsize']=36\n",
    "    mp.rcParams['ytick.labelsize']=36\n",
    "#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#\n",
    "\n",
    "    filelist=list\n",
    "    files=sorted(glob(filelist))\n",
    "    julian_date=np.zeros(len(files))\n",
    "    for i in range(0,len(files)):\n",
    "        data=fits.open(files[i])\n",
    "        header=data[0].header\n",
    "        julian_date[i]=header['JD']\n",
    "\n",
    "    # This is a dummy for flux\n",
    "    julian_date-=2458870\n",
    "    fluxd=np.zeros(len(files))\n",
    "    for i in range(len(files)):\n",
    "        fluxd[i]=random.randint(5,8)\n",
    "\n",
    "    flux_errd=0.05*fluxd\n",
    "    #################\n",
    "\n",
    "    fig,ax=plt.subplots()\n",
    "    plt.plot(julian_date,fluxd,'ko',markersize=25,label='I band light curve')\n",
    "    plt.errorbar(julian_date,fluxd,yerr=flux_errd,capsize=5,fmt=' ')\n",
    "    plt.ylabel(\"Magnitudes (instrumental)\",fontsize=36)\n",
    "    plt.xlabel(\"Julian Dates (2458870+)\",fontsize=36)\n",
    "    #plt.title(\"Deviation in lag values with respect to the variability\",fontsize=36)\n",
    "    ax.tick_params(axis=\"both\",which='minor',direction=\"in\")\n",
    "    ax.tick_params(axis=\"both\",which='major',direction=\"in\")\n",
    "    ax.yaxis.set_ticks_position('both')\n",
    "    ax.xaxis.set_ticks_position('both')\n",
    "    ax.minorticks_on()\n",
    "    ax.legend(fontsize=30)\n",
    "    plt.show()\n"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}