{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "MiKvRj5u076V"
   },
   "source": [
    "## **ABOUT**\n",
    "This example illustrates the 2D multi-slice, Spin Echo (SE) acquisition using the `pypulseq` library. This sequence is typically used for T<sub>2</sub> weighted imaging. A 2D Fourier transform can be used to reconstruct images from this acquisition. Read more about SE [here](http://mriquestions.com/se-vs-multi-se-vs-fse.html).\n",
    "\n",
    "**Contact**: For issues, write to ks3621@columbia.edu\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Y98YDJr215fa"
   },
   "source": [
    "## **INSTALL** `pypulseq`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ogKNAZH3TmgA"
   },
   "outputs": [],
   "source": [
    "!pip install git+https://github.com/imr-framework/pypulseq.git@dev"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "UgqzEwle2xCd"
   },
   "source": [
    "## **IMPORT PACKAGES**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "3X7UsV832B6j"
   },
   "outputs": [],
   "source": [
    "from math import pi\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from pypulseq.Sequence.sequence import Sequence\n",
    "from pypulseq.calc_duration import calc_duration\n",
    "from pypulseq.make_adc import make_adc\n",
    "from pypulseq.make_delay import make_delay\n",
    "from pypulseq.make_sinc_pulse import make_sinc_pulse\n",
    "from pypulseq.make_trapezoid import make_trapezoid\n",
    "from pypulseq.opts import Opts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "UQ4AWw9l4et_"
   },
   "source": [
    "## **USER INPUTS**\n",
    "\n",
    "These parameters are typically on the user interface of the scanner computer console "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ssnNwiQH4q_0"
   },
   "outputs": [],
   "source": [
    "nsa = 1  # Number of averages\n",
    "n_slices = 3  # Number of slices\n",
    "Nx = 128\n",
    "Ny = 128\n",
    "fov = 220e-3  # mm\n",
    "slice_thickness = 5e-3  # s\n",
    "slice_gap = 15e-3  # s\n",
    "rf_flip = 90  # degrees\n",
    "rf_offset = 0\n",
    "print('User inputs setup')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "PeYeI0V45ZfD"
   },
   "source": [
    "## **SYSTEM LIMITS**\n",
    "Set the hardware limits and initialize sequence object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "XHs1LT965kqg"
   },
   "outputs": [],
   "source": [
    "system = Opts(max_grad=32, grad_unit='mT/m', max_slew=130, slew_unit='T/m/s', \n",
    "              grad_raster_time=10e-6, rf_ringdown_time=10e-6, \n",
    "              rf_dead_time=100e-6)\n",
    "seq = Sequence(system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ee-xBrpa7Zyn"
   },
   "source": [
    "## **TIME CONSTANTS**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "u2dW2nRf7obq"
   },
   "outputs": [],
   "source": [
    "TE = 100e-3  # s\n",
    "TR = 3  # s\n",
    "tau = TE / 2  # s\n",
    "readout_time = 6.4e-3\n",
    "pre_time = 8e-4  # s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "OTw7M03g79bH"
   },
   "source": [
    "## **RF**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "XDZyQrbL8I3Q"
   },
   "outputs": [],
   "source": [
    "flip90 = round(rf_flip * pi / 180, 3)\n",
    "flip180 = 180 * pi / 180\n",
    "rf90, gz90, _ = make_sinc_pulse(flip_angle=flip90, system=system, duration=4e-3, \n",
    "                                slice_thickness=slice_thickness, apodization=0.5, \n",
    "                                time_bw_product=4, return_gz = True)\n",
    "rf180, gz180, _ = make_sinc_pulse(flip_angle=flip180, system=system, \n",
    "                                  duration=2.5e-3, \n",
    "                                  slice_thickness=slice_thickness, \n",
    "                                  apodization=0.5, \n",
    "                                time_bw_product=4, phase_offset=90 * pi/180, return_gz = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "RFSHuUOG9LHK"
   },
   "source": [
    "## **READOUT**\n",
    "Readout gradients and related events"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Q8p-CttI9dk9"
   },
   "outputs": [],
   "source": [
    "delta_k = 1 / fov\n",
    "k_width = Nx * delta_k\n",
    "gx = make_trapezoid(channel='x', system=system, flat_area=k_width, \n",
    "                    flat_time=readout_time)\n",
    "adc = make_adc(num_samples=Nx, duration=gx.flat_time, delay=gx.rise_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "o829kzm8kVFB"
   },
   "source": [
    "## **PREPHASE AND REPHASE**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "m5zA1bMakTVs"
   },
   "outputs": [],
   "source": [
    "phase_areas = (np.arange(Ny) - (Ny / 2)) * delta_k\n",
    "gz_reph = make_trapezoid(channel='z', system=system, area=-gz90.area / 2,\n",
    "                         duration=2.5e-3)\n",
    "gx_pre = make_trapezoid(channel='x', system=system, flat_area=k_width / 2, \n",
    "                        flat_time=readout_time / 2)\n",
    "gy_pre = make_trapezoid(channel='y', system=system, area=phase_areas[-1], \n",
    "                        duration=2e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "5Css5esAkYHo"
   },
   "source": [
    "## **SPOILER**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "R1DOmoKKkawr"
   },
   "outputs": [],
   "source": [
    "gz_spoil = make_trapezoid(channel='z', system=system, area=gz90.area * 4,\n",
    "                          duration=pre_time * 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3F5JUpE9-4lo"
   },
   "source": [
    "## **DELAYS**\n",
    "Echo time (TE) and repetition time (TR). Here, TE is broken down into `delay1` and `delay2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "aOKRJclb_mDQ"
   },
   "outputs": [],
   "source": [
    "delay1 = tau - calc_duration(rf90) / 2 - calc_duration(gx_pre)\n",
    "delay1 -= calc_duration(gz_spoil) - calc_duration(rf180) / 2\n",
    "delay1 = make_delay(delay1)\n",
    "delay2 = tau - calc_duration(rf180) / 2 - calc_duration(gz_spoil)\n",
    "delay2 -= calc_duration(gx) / 2\n",
    "delay2 = make_delay(delay2)\n",
    "delay_TR = TR - calc_duration(rf90) / 2 - calc_duration(gx) / 2 - TE\n",
    "delay_TR -= calc_duration(gy_pre)\n",
    "delay_TR = make_delay(delay_TR)\n",
    "print(f'delay_1: {delay1}')\n",
    "print(f'delay_2: {delay1}')\n",
    "print(f'delay_TR: {delay_TR}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "6Dq4wT-UAEOR"
   },
   "source": [
    "## **CONSTRUCT SEQUENCE**\n",
    "Construct sequence for one phase encode and multiple slices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "B8ZmVkkrAXnK"
   },
   "outputs": [],
   "source": [
    "# Prepare RF offsets. This is required for multi-slice acquisition\n",
    "delta_z = n_slices * slice_gap\n",
    "z = np.linspace((-delta_z / 2), (delta_z / 2), n_slices) + rf_offset\n",
    "\n",
    "for k in range(nsa):  # Averages\n",
    "  for j in range(n_slices):  # Slices\n",
    "    # Apply RF offsets\n",
    "    freq_offset = gz90.amplitude * z[j]\n",
    "    rf90.freq_offset = freq_offset\n",
    "\n",
    "    freq_offset = gz180.amplitude * z[j]\n",
    "    rf180.freq_offset = freq_offset\n",
    "\n",
    "    for i in range(Ny):  # Phase encodes\n",
    "      seq.add_block(rf90, gz90)\n",
    "      gy_pre = make_trapezoid(channel='y', system=system, \n",
    "                              area=phase_areas[-i -1], duration=2e-3)\n",
    "      seq.add_block(gx_pre, gy_pre, gz_reph)\n",
    "      seq.add_block(delay1)\n",
    "      seq.add_block(gz_spoil)\n",
    "      seq.add_block(rf180, gz180)\n",
    "      seq.add_block(gz_spoil)\n",
    "      seq.add_block(delay2)\n",
    "      seq.add_block(gx, adc)\n",
    "      gy_pre = make_trapezoid(channel='y', system=system, \n",
    "                              area=-phase_areas[-j -1], duration=2e-3)\n",
    "      seq.add_block(gy_pre, gz_spoil)\n",
    "      seq.add_block(delay_TR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "l-YP9djBJCpC"
   },
   "source": [
    "## **PLOTTING TIMNG DIAGRAM**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "d_iCUR4nfoH9"
   },
   "outputs": [],
   "source": [
    "seq.plot(time_range=(0, 0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "fYNgdWc_KiK7"
   },
   "source": [
    "## **GENERATING `.SEQ` FILE**\n",
    "Uncomment the code in the cell below to generate a `.seq` file and download locally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "6iN0aeuuqKRe"
   },
   "outputs": [],
   "source": [
    "# seq.write('t2_se_pypulseq_colab.seq')  # Save to disk\n",
    "# from google.colab import files\n",
    "# files.download('t2_se_pypulseq_colab.seq')  # Download locally"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "4Q0b5w-lKtfP"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "write_t2_se.ipynb",
   "private_outputs": true,
   "provenance": []
  },
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}