gprMax/tools/Jupyter_notebooks/example_Ascan_2D.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A-scan from a metal cylinder (2D)\n",
    "\n",
    "This example is the GPR modelling equivalent of 'Hello World'! It demonstrates how to simulate a single trace (A-scan) from a metal cylinder buried in a dielectric half-space. The input needed to create the model is:\n",
    "\n",
    "### my_cylinder_Ascan_2D.in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%writefile ../../user_models/cylinder_Ascan_2D.in\n",
    "#title: A-scan from a metal cylinder buried in a dielectric half-space\n",
    "#domain: 0.240 0.210 0.002\n",
    "#dx_dy_dz: 0.002 0.002 0.002\n",
    "#time_window: 3e-9\n",
    "\n",
    "#material: 6 0 1 0 half_space\n",
    "\n",
    "#waveform: ricker 1 1.5e9 my_ricker\n",
    "#hertzian_dipole: z 0.100 0.170 0 my_ricker\n",
    "#rx: 0.140 0.170 0\n",
    "\n",
    "#box: 0 0 0 0.240 0.170 0.002 half_space\n",
    "#cylinder: 0.120 0.080 0 0.120 0.080 0.002 0.010 pec\n",
    "\n",
    "#geometry_view: 0 0 0 0.240 0.210 0.002 0.002 0.002 0.002 cylinder_half_space n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geometry of a metal cylinder buried in a dielectric half-space\n",
    "<img style=\"float: left\" src=\"cylinder_half_space_geo.png\" width=\"50%\" height=\"50%\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The geometry of the scenario is straightforward - the transparent area around the boundary of the domain represents the PML region. The red cell is the source and the blue cell is the receiver.\n",
    "\n",
    "For this initial example a detailed description of what each command in the input file does and why each command was used is given. The following steps explain the process of building the input file:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Determine the constitutive parameters for the materials\n",
    "\n",
    "There will be three different materials in the model representing air, the dielectric half-space, and the metal cylinder. Air (free space) already exists as a built-in material in gprMax which can be accessed using the ``free_space`` identifier. The metal cylinder will be modelled as a Perfect Electric Conductor, which again exists as a built-in material in gprMax and can be accessed using the ``pec`` identifier. So the only material which has to be defined is for the dielectric half-space. It is a non-magnetic material, i.e. $\\mu_r=1$ and $\\sigma_*=0$ and with a relative permittivity of six, $\\epsilon_r=6$, and zero conductivity, $\\sigma=0$. The identifier ``half_space`` will be used.\n",
    "\n",
    "    #material: 6 0 1 0 half_space\n",
    "\n",
    "### 2. Determine the source type and excitation frequency\n",
    "\n",
    "These should generally be known, often based on the GPR system or scenario being modelled. Low frequencies are used where significant penetration depth is important, whereas high frequencies are used where less penetration and better resolution are required. In this case a theoretical Hertzian dipole source fed with a Ricker waveform with a centre frequency of $f_c=1.5~\\textrm{GHz}$ will be used to simulate the GPR antenna (see the bowtie antenna example model for how to include a model of the actual GPR antenna in the simulation).\n",
    "\n",
    "    #waveform: ricker 1 1.5e9 my_ricker\n",
    "    #hertzian_dipole: z 0.100 0.170 0 my_ricker\n",
    "\n",
    "The Ricker waveform is created with the ``#waveform`` command, specifying an amplitude of one, centre frequency of 1.5 GHz and picking an arbitrary identifier of ``my_ricker``. The Hertzian dipole source is created using the ``#hertzian_dipole`` command, specifying a z direction polarisation (perpendicular to the survey direction if a B-scan were being created), location on the surface of the slab, and using the Ricker waveform already created.\n",
    "\n",
    "### 3. Calculate a spatial resolution and domain size\n",
    "\n",
    "In the guidance section it was stated that a good *rule-of-thumb* was that the spatial resolution should be one tenth of the smallest wavelength present in the model. To determine the smallest wavelength, the highest frequency and lowest velocity present in the model are required. The highest frequency is not the centre frequency of the Ricker waveform!  \n",
    "\n",
    "You can use the following code to plot builtin waveforms and their FFTs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from gprMax.waveforms import Waveform\n",
    "from tools.plot_source_wave import check_timewindow, mpl_plot\n",
    "\n",
    "w = Waveform()\n",
    "w.type = 'ricker'\n",
    "w.amp = 1\n",
    "w.freq = 1.5e9\n",
    "timewindow = 3e-9\n",
    "dt = 1.926e-12\n",
    "\n",
    "timewindow, iterations = check_timewindow(timewindow, dt)\n",
    "plt = mpl_plot(w, timewindow, dt, iterations, fft=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By examining the spectrum of a Ricker waveform it is evident much higher frequencies are present, i.e. at a level -40dB from the centre frequency, frequencies 2-3 times as high are present. In this case the highest significant frequency present in the model is likely to be around 4 GHz. To calculate the wavelength at 4 GHz in the half-space (which has the lowest velocity) use:\n",
    "\n",
    "$$\\lambda = \\frac{c}{f \\sqrt{\\epsilon_r}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from math import sqrt\n",
    "\n",
    "# Speed of light in vacuum (m/s)\n",
    "c = 299792458\n",
    "\n",
    "# Highest relative permittivity present in model\n",
    "er = 6\n",
    "\n",
    "# Maximum frequency present in model\n",
    "fmax = 4e9\n",
    "\n",
    "# Minimum wavelength\n",
    "wmin = c / (fmax * sqrt(er))\n",
    "\n",
    "# Maximum spatial resolution\n",
    "dmin = wmin / 10\n",
    "\n",
    "print('Minimum wavelength: {:g} m'.format(wmin))\n",
    "print('Maximum spatial resolution: {:g} m'.format(dmin))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This would give a minimum spatial resolution of 3 mm. However, the diameter of the cylinder is 20 mm so would be resolved to 7 cells. Therefore a better choice would be 2 mm which resolves the diameter of the rebar to 10 cells.\n",
    "\n",
    "    #dx_dy_dz: 0.002 0.002 0.002\n",
    "\n",
    "The domain size should be enough to enclose the volume of interest, plus allow 10 cells (if using the default value) for the PML absorbing boundary conditions and approximately another 10 cells of between the PML and any objects of interest. In this case the plan is to take a B-scan of the scenario (in the next example) so the domain should be large enough to do that. Although this is a 2D model one cell must be specified in the infinite direction (in this case the z direction) of the domain.\n",
    "\n",
    "    #domain: 0.240 0.210 0.002\n",
    "\n",
    "### 4. Choose a time window\n",
    "\n",
    "It is desired to see the reflection from the cylinder, therefore the time window must be long enough to allow the electromagnetic waves to propagate from the source through the half-space to the cylinder and be reflected back to the receiver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d = 0.090\n",
    "t = (2 * d) / (c / sqrt(6))\n",
    "print('Minimum time window: {:g} s'.format(t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the minimum time required, but the source waveform has a width of 1.2 ns, to allow for the entire source waveform to be reflected back to the receiver an initial time window of 3 ns will be tested.\n",
    "\n",
    "    #time_window: 3e-9\n",
    "\n",
    "The time step required for the model is automatically calculated using the CFL condition (for this case in 2D).\n",
    "\n",
    "### 5. Create the objects\n",
    "\n",
    "Now physical objects can be created for the half-space and the cylinder. First the ``#box`` command will be used to create the half-space and then the ``#cylinder`` command will be given which will overwrite the properties of the half-space with those of the cylinder at the location of the cylinder.\n",
    "\n",
    "    #box: 0 0 0 0.240 0.170 0.002 half_space\n",
    "    #cylinder: 0.120 0.080 0 0.120 0.080 0.002 0.010 pec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the model\n",
    "\n",
    "You can now run the model using:\n",
    "    \n",
    "    python -m gprMax user_models/cylinder_Ascan_2D.in\n",
    "\n",
    "**Tip**: You can use the ``--geometry-only`` command line argument to build a model and produce any geometry views but not run the simulation. This option is useful for checking the geometry of the model is correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "from gprMax.gprMax import api\n",
    "\n",
    "filename = os.path.join(os.pardir, os.pardir, 'user_models', 'cylinder_Ascan_2D.in')\n",
    "api(filename, n=1, geometry_only=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## View the results\n",
    "\n",
    "### Plot the A-scan\n",
    "\n",
    "You should have produced an output file ``cylinder_Ascan_2D.out``. You can view the results using:\n",
    "\n",
    "    python -m tools.plot_Ascan user_models/cylinder_Ascan_2D.out\n",
    "    \n",
    "You can use the following code to experiment with plotting different field/current components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "from gprMax.receivers import Rx\n",
    "from tools.plot_Ascan import mpl_plot\n",
    "\n",
    "filename = os.path.join(os.pardir, os.pardir, 'user_models', 'cylinder_Ascan_2D.out')\n",
    "outputs = Rx.defaultoutputs\n",
    "#outputs = ['Ez']\n",
    "plt = mpl_plot(filename, outputs, fft=False)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}