get_harmonic_phase_offset.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import fftpack\n",
    "import os\n",
    "import subprocess\n",
    "from astropy.io import fits\n",
    "from astropy.table import Table, Column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "################################################################################\n",
    "def fits_in(in_file, meta_dict, test=False):\n",
    "    \"\"\"\n",
    "    Read in an eventlist in .fits format to make the cross spectrum. Read\n",
    "    in a clock-corrected GTI'd event list, populate the light curves, compute\n",
    "    cross spectrum per energy channel and keep running sum (later to be an\n",
    "    average) of the cross spectra.\n",
    "\n",
    "    Assumes that the reference band times are the center of the bins, and the\n",
    "    CI times are at the front of the bins.\n",
    "\n",
    "    I take the approach: start time <= segment < end_time, to avoid double-\n",
    "    counting and/or skipping events.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    in_file : str\n",
    "        The full path of the FITS data file being analyzed.\n",
    "\n",
    "    ref_band_file : str\n",
    "        Name of FITS optical or IR data file for reference band. This one file\n",
    "        has the reference band for the whole data set. Gaps are ok.\n",
    "\n",
    "    meta_dict : dict\n",
    "        Dictionary of necessary meta-parameters for data analysis.\n",
    "\n",
    "    test : boolean\n",
    "        True if only running one segment of data for testing, False if analyzing\n",
    "        the whole data file. Default=False\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    cross_spec:  np.array of floats\n",
    "        3-D array of the raw cross spectrum, per segment.\n",
    "        Dimensions: [n_bins, detchans, n_seg]\n",
    "\n",
    "    ci_whole : ccf_lc.Lightcurve object\n",
    "        Channel of interest for this data file.\n",
    "\n",
    "    ref_whole : ccf_lc.Lightcurve object\n",
    "        Reference band for this data file.\n",
    "\n",
    "    n_seg : int\n",
    "        Number of segments in this data file.\n",
    "\n",
    "    dt_whole : np.array of floats\n",
    "        1-D array of timestep between light curve bins for each segment. These\n",
    "        will be different if artificially adjusting the QPO frequency in Fourier\n",
    "        space (currently only doing that per data file, not per segment).\n",
    "\n",
    "    df_whole : np.array of floats\n",
    "        1-D array of frequency step between Fourier bins for each segment. These\n",
    "        will be different if artificially adjusting the QPO frequency in Fourier\n",
    "        space (currently only doing that per data file, not per segment).\n",
    "\n",
    "    exposure : float\n",
    "        The total (used) exposure of the data file.\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    assert tools.power_of_two(meta_dict['n_bins']), \"ERROR: n_bins must be a \"\\\n",
    "            \"power of 2.\"\n",
    "    meta_dict['obs_epoch'] = tools.obs_epoch_rxte(in_file)\n",
    "\n",
    "    print(\"Input file: %s\" % in_file)\n",
    "\n",
    "    ## Determining print iterator for segments\n",
    "    if meta_dict['n_bins'] == 32768:\n",
    "        print_iterator = int(10)\n",
    "    elif meta_dict['n_bins'] < 32768:\n",
    "        print_iterator = int(10)   #TODO: change back to 20\n",
    "    else:\n",
    "        print_iterator = int(1)\n",
    "\n",
    "    #######################################################\n",
    "    ## Check if the FITS file exists; if so, load the data\n",
    "    #######################################################\n",
    "\n",
    "    time = np.asarray([])\n",
    "    channel = np.asarray([])\n",
    "    pcuid = np.asarray([])\n",
    "\n",
    "    ## Reading an event list from an astropy table FITS file\n",
    "    # try:\n",
    "    #     data_table = Table.read(in_file)\n",
    "    #     time = data_table['TIME']\n",
    "    #     channel = data_table['CHANNEL']\n",
    "    #     pcuid = data_table['PCUID']\n",
    "    # except IOError:\n",
    "    #     print(\"\\tERROR: File does not exist: %s\" % in_file)\n",
    "    #     exit()\n",
    "\n",
    "    ## Reading an event list from a normal FITS table\n",
    "    try:\n",
    "        fits_hdu = fits.open(in_file)\n",
    "        time = fits_hdu[1].data.field('TIME')  ## Data is in ext 1\n",
    "        channel = fits_hdu[1].data.field('CHANNEL')\n",
    "        pcuid = fits_hdu[1].data.field('PCUID')\n",
    "        fits_hdu.close()\n",
    "    except IOError:\n",
    "        print(\"\\tERROR: File does not exist: %s\" % in_file)\n",
    "        exit()\n",
    "\n",
    "    # try:\n",
    "    #     fits_hdu = fits.open(in_file)\n",
    "    # except IOError:\n",
    "    #     print(\"\\tERROR: File does not exist: %s\" % in_file)\n",
    "    #     sys.exit()\n",
    "    #\n",
    "    # data = fits_hdu[1].data\n",
    "    # fits_hdu.close()\n",
    "\n",
    "    ###################\n",
    "    ## Initializations\n",
    "    ###################\n",
    "\n",
    "    n_seg = 0\n",
    "    ci_whole = ccf_lc.Lightcurve(n_bins=meta_dict['n_bins'],\n",
    "            detchans=meta_dict['detchans'], type='ci')\n",
    "    ref_whole = ccf_lc.Lightcurve(n_bins=meta_dict['n_bins'],\n",
    "            detchans=meta_dict['detchans'], type='ref')\n",
    "    cs_whole = np.zeros((meta_dict['n_bins'], meta_dict['detchans'], 1),\n",
    "            dtype=np.complex128)\n",
    "    dt_whole = np.array([])\n",
    "    df_whole = np.array([])\n",
    "    exposure = 0\n",
    "    # print(set(pcuid))\n",
    "    # exit()\n",
    "    start_time = time[0]\n",
    "    final_time = time[-1]\n",
    "\n",
    "    ###################################\n",
    "    ## Selecting PCU for interest band\n",
    "    ###################################\n",
    "\n",
    "    PCU2_mask = pcuid == 2\n",
    "    time_pcu2 = time[PCU2_mask]\n",
    "    chan_pcu2 = channel[PCU2_mask]\n",
    "\n",
    "    all_time_ci = np.asarray(time_pcu2, dtype=np.float64)\n",
    "    all_energy_ci = np.asarray(chan_pcu2, dtype=np.float64)\n",
    "\n",
    "    ######################################\n",
    "    ## Getting reference band light curve\n",
    "    ######################################\n",
    "\n",
    "    refpcu_mask = pcuid != 2\n",
    "    all_time_ref = np.asarray(time[refpcu_mask], dtype=np.float64)\n",
    "    all_energy_ref = np.asarray(channel[refpcu_mask], dtype=np.float64)\n",
    "    all_rate_ref = None\n",
    "    all_err_ref = None\n",
    "\n",
    "    seg_end_time = start_time + meta_dict['n_seconds']\n",
    "\n",
    "    ############################\n",
    "    ## Looping through segments\n",
    "    ############################\n",
    "\n",
    "    print(\"Segments computed:\")\n",
    "\n",
    "    while (seg_end_time + (meta_dict['adjust_seg'] * meta_dict['dt'])) <= \\\n",
    "            final_time:\n",
    "\n",
    "        ## Adjusting segment length to artificially line up the QPOs\n",
    "        seg_end_time += (meta_dict['adjust_seg'] * meta_dict['dt'])\n",
    "\n",
    "        ## Get events for channels of interest\n",
    "        time_ci = all_time_ci[np.where(all_time_ci < seg_end_time)]\n",
    "        energy_ci = all_energy_ci[np.where(all_time_ci < seg_end_time)]\n",
    "\n",
    "        ## Chop current segment off the rest of the list\n",
    "        for_next_iteration_ci = np.where(all_time_ci >= seg_end_time)\n",
    "        all_time_ci = all_time_ci[for_next_iteration_ci]\n",
    "        all_energy_ci = all_energy_ci[for_next_iteration_ci]\n",
    "\n",
    "        ## Get events for reference band\n",
    "        time_ref = all_time_ref[np.where(all_time_ref < seg_end_time)]\n",
    "        if not meta_dict['ref_file']:\n",
    "            energy_ref = all_energy_ref[np.where(all_time_ref < seg_end_time)]\n",
    "            rate_ref = [0]\n",
    "        else:\n",
    "            rate_ref = all_rate_ref[np.where(all_time_ref < seg_end_time)]\n",
    "            err_ref = all_err_ref[np.where(all_time_ref < seg_end_time)]\n",
    "\n",
    "        ## Chop current segment off the rest of the list\n",
    "        for_next_iteration_ref = np.where(all_time_ref >= seg_end_time)\n",
    "        all_time_ref = all_time_ref[for_next_iteration_ref]\n",
    "        if not meta_dict['ref_file']:\n",
    "            all_energy_ref = all_energy_ref[for_next_iteration_ref]\n",
    "        else:\n",
    "            all_rate_ref = all_rate_ref[for_next_iteration_ref]\n",
    "            all_err_ref = all_err_ref[for_next_iteration_ref]\n",
    "\n",
    "        ########################################################################\n",
    "        ## At the end of a segment, populate light curve and make cross spectrum\n",
    "        ########################################################################\n",
    "\n",
    "        if len(time_ci) > 0 and \\\n",
    "                (len(time_ref) > 0 or\n",
    "                (meta_dict['ref_file'] and \\\n",
    "                len(rate_ref) == meta_dict['n_bins'])):\n",
    "\n",
    "            ##############################################################\n",
    "            ## Populate the light curves for interest and reference bands\n",
    "            ##############################################################\n",
    "\n",
    "            rate_ci_2d = tools.make_2Dlightcurve(np.asarray(time_ci),\n",
    "                    np.asarray(energy_ci), meta_dict['n_bins'],\n",
    "                    meta_dict['detchans'], start_time, seg_end_time)\n",
    "\n",
    "            if not meta_dict['ref_file']:\n",
    "                rate_ref_2d = tools.make_2Dlightcurve( np.asarray(time_ref),\n",
    "                        np.asarray(energy_ref), meta_dict['n_bins'],\n",
    "                        meta_dict['detchans'], start_time, seg_end_time)\n",
    "                ## Stack the reference band\n",
    "                rate_ref = stack_reference_band(rate_ref_2d, instrument=\"PCA\",\n",
    "                        obs_epoch=meta_dict['obs_epoch'])\n",
    "\n",
    "            ###########################\n",
    "            ## Make the cross spectrum\n",
    "            ###########################\n",
    "\n",
    "            cs_seg, ci_seg, ref_seg = make_cs(rate_ci_2d, rate_ref, meta_dict)\n",
    "\n",
    "            ######################################################\n",
    "            ## Only keep and use segments where the variance > 0.\n",
    "            ######################################################\n",
    "\n",
    "            ## Append segment to arrays\n",
    "            fft_ci_whole = np.dstack((fft_ci_whole, fft_ci_seg))\n",
    "            fft_ref_whole = np.dstack((fft_ref_whole, fft_ref_seg))\n",
    "\n",
    "            ## Sum across segments -- arrays, so it adds by index\n",
    "            exposure += (seg_end_time - start_time)\n",
    "            n_seg += 1\n",
    "            ci_whole.mean_rate += ci_seg.mean_rate\n",
    "            ref_whole.mean_rate += ref_seg.mean_rate\n",
    "\n",
    "            if n_seg % print_iterator == 0:\n",
    "                print(\"\\t\", n_seg)\n",
    "\n",
    "            if test is True and n_seg == 1:  # For testing\n",
    "                break\n",
    "\n",
    "            start_time = seg_end_time\n",
    "            seg_end_time += meta_dict['n_seconds']\n",
    "\n",
    "        ## This next bit deals with gappy data\n",
    "        else:\n",
    "            start_time = max(all_time_ci[0], all_time_ref[0])\n",
    "            seg_end_time = start_time + meta_dict['n_seconds']\n",
    "\n",
    "        ## End of 'if there are counts in this segment'\n",
    "\n",
    "    ## End of while-loop\n",
    "\n",
    "    fft_ci_whole = fft_ci_whole[:,:,1:]\n",
    "    fft_ref_whole = fft_ref_whole[:,:,1:]\n",
    "\n",
    "\n",
    "    return fft_ci_whole, fft_ref_whole, n_seg"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}