Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add scintillometry walk through #84

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
389 changes: 389 additions & 0 deletions examples/Scintillometry-walk-through.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,389 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Baseband and Scintillometry Walk Through"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Welcome to the Scintillometry package walk through notebook. Here you will learn the knowledge of pulsar and the basic usage of baseband and scintillometry python packages."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Outline\n",
"* Introducation\n",
"* Read baseband data using baseband pacakge\n",
"* Dedispersion\n",
"* Integration\n",
" * Waterfall plot\n",
" * Folding\n",
" * Polarizations \n",
"* Output"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sample data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The sample data used in this notebook are stored at:\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Baseband data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Radio telescope receives the electromagnetic waves(EM waves) from far away sources. The EM wave field creates voltage in the telescope antenna. These voltage signals are amplified, filtered and than digitized by the backend electronics. Baseband data are fourier transformed voltage signals. Thus they are complex numbers. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Structure of baseband data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Polarizations \n",
"\n",
"The EM waves has polarizations. The radio antennas are generally designed to receive the polarized signal by two orthogonal linear polarization channels X and Y. From these two polarization channels, the [Stokes parameters]( https://en.wikipedia.org/wiki/Stokes_parameters) can be calculated. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Pulsars "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Useful links: <br> \n",
"For pulsar ephemeris (Pulsar motion parameters)<br> \n",
"http://www.atnf.csiro.au/research/pulsar/psrcat/ <br>\n",
"For pulse profile <br> \n",
"http://www.epta.eu.org/epndb/"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Pulsar data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Pulsar data are periodical signals "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Pulsar data process"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Baseband --> dedispersion --> Resample --> fold --> Dynamic spectrum --> Secondary Spectrum"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Required python package\n",
"1. [Baseband](https://github.com/mhvk/baseband) \n",
"2. [Numpy](http://www.numpy.org/) \n",
"3. [Astropy](http://www.astropy.org/) \n",
"4. [Scintillometry](https://github.com/mhvk/scintillometry) \n",
"#### Optional python package \n",
"5. [PINT](https://github.com/nanograv/PINT) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Read baseband data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"1. Decide your data types\n",
" * Supported data types\n",
" * [VDIF](https://baseband.readthedocs.io/en/stable/vdif/index.html)\n",
" * [MARK 5B](https://baseband.readthedocs.io/en/stable/mark5b/index.html)\n",
" * [MARK 4](https://baseband.readthedocs.io/en/stable/mark4/index.html)\n",
" * [DADA](https://baseband.readthedocs.io/en/stable/dada/index.html)\n",
" * [GUPPI](https://baseband.readthedocs.io/en/stable/guppi/index.html)\n",
" * [GSB](https://baseband.readthedocs.io/en/stable/gsb/index.html)\n",
"\n",
"\n",
"Note: There is a detailed example code in the link to each file types."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"2. Import the necessary module. For instance, we use GUPPI file type"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np \n",
"import astropy.units as u # import astropy units module, this helps you on the physical quantity unit convers\n",
"from baseband import guppi # import the guppi baseband reader \n",
"import os # This package will help you navigating the file system. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"3. Create the file handler for your files. \n",
"<br>3.1 Read a single file."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# I think the baseband part should go to the baseband package.\n",
"data_dir = '/home/luo/Projects/baseband/baseband/data'\n",
"data_file = os.path.join(data_dir, 'sample_puppi.raw')\n",
"fh = guppi.open(data_file, 'rb') # 'rb' means read binary file, since the data file is in binary format. "
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"ename": "UnicodeDecodeError",
"evalue": "'ascii' codec can't decode byte 0xf6 in position 0: ordinal not in range(128)",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mUnicodeDecodeError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-12-1d65c06704ee>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mfh\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread_header\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m~/.local/lib/python3.6/site-packages/baseband-2.0.0-py3.6.egg/baseband/guppi/base.py\u001b[0m in \u001b[0;36mread_header\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 98\u001b[0m \u001b[0mheader\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m`\u001b[0m\u001b[0;34m~\u001b[0m\u001b[0mbaseband\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mguppi\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mGUPPIHeader\u001b[0m\u001b[0;31m`\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 99\u001b[0m \"\"\"\n\u001b[0;32m--> 100\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mGUPPIHeader\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfromfile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfh_raw\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 101\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 102\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mread_frame\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmemmap\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mverify\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/.local/lib/python3.6/site-packages/baseband-2.0.0-py3.6.egg/baseband/guppi/header.py\u001b[0m in \u001b[0;36mfromfile\u001b[0;34m(cls, fh, verify)\u001b[0m\n\u001b[1;32m 116\u001b[0m \u001b[0mline\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'========='\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 117\u001b[0m \u001b[0;32mwhile\u001b[0m \u001b[0mline\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m8\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32min\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m'='\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m' '\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 118\u001b[0;31m \u001b[0mline\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfh\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m80\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdecode\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'ascii'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 119\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mline\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'END'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 120\u001b[0m \u001b[0;32mbreak\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mUnicodeDecodeError\u001b[0m: 'ascii' codec can't decode byte 0xf6 in position 0: ordinal not in range(128)"
]
}
],
"source": [
"fh.read_header()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>3.2 Read mulitple files together. \n",
"Note: Those files should be contiguous and sorted by time. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"4. Read and exam the header information"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"GUPPIFileReader(fh_raw=<_io.BufferedReader name='/home/luo/Projects/baseband/baseband/data/sample_puppi.raw'>)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fh."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\"\"\"Test of the scintillometry dedispersion routines on PSR B1937+21\"\"\"\n",
"import numpy as np\n",
"import astropy.units as u\n",
"from baseband import guppi\n",
"from scintillometry.dispersion import Dedisperse\n",
"from scintillometry.functions import Square\n",
"from scintillometry.fourier import get_fft_maker\n",
"import matplotlib.pylab as plt\n",
"\n",
"\n",
"FFT = get_fft_maker('pyfftw', flags=['FFTW_ESTIMATE'], threads=4)\n",
"\n",
"dm = 71.02227638 * u.pc / u.cm**3\n",
"f0 = (641.928224582360 + 2.78557046583437451/60.) * u.Hz # from polyco\n",
"\n",
"fh = guppi.open('puppi_58245_B1937+21_0799.0010.raw')\n",
"frequency = (fh.header0['OBSFREQ']\n",
" - fh.header0['OBSBW'] / 2. # Down to bottom end\n",
" + fh.header0['CHAN_BW'] / 2. # Up to center of channel 0\n",
" + fh.header0['CHAN_BW'] * np.arange(4)) * u.MHz\n",
"\n",
"dedisperse = Dedisperse(fh, dm, frequency=frequency, sideband=1,\n",
" FFT=FFT)\n",
"\n",
"square = Square(dedisperse)\n",
"\n",
"nbin = 64\n",
"lc = np.zeros(nbin)\n",
"for i in range(10):\n",
" print('At ', square.tell(unit=u.s))\n",
" pre = square.tell()\n",
" p = square.read(3125000).sum((1, 2))\n",
" post = square.tell()\n",
" ph = (np.arange(pre, post) / fh.sample_rate * f0).to_value(u.one)\n",
" ph = (ph - 0.35) % 1.\n",
" lc += np.bincount((ph * nbin).astype(int), weights=p)\n",
"\n",
"plt.ion()\n",
"plt.plot(np.arange(nbin)/nbin, lc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Useful links: <br>\n",
"&nbsp;&nbsp;&nbsp;&nbsp;Baseband documentation: https://baseband.readthedocs.io/en/stable/ \n",
"<br>\n",
"&nbsp;&nbsp;&nbsp;&nbsp;Baseband github: https://github.com/mhvk/baseband"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Dedispersion"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Integration"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Output"
]
},
{
"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.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}