cc_from_many_sources.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise to simulate random signals at to receivers and to investigate how the source distribution affects the cross-correlation\n",
    "\n",
    "## Main Python packages used: \n",
    "\n",
    "- [**NumPy**](https://numpy.org) : for mathematical functions\n",
    "\n",
    "- [**Matplotlib**](https://matplotlib.org) : for plotting results\n",
    "\n",
    "- [**SciPy**](https://github.com/obspy/obspy/wiki) : to compute [**cross-correlations**](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.correlate.html)\n",
    "\n",
    "---\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d52ccebb8b75432280ac971395418734",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "30\n",
      "60\n",
      "90\n",
      "120\n",
      "150\n",
      "180\n",
      "210\n",
      "240\n",
      "270\n",
      "300\n",
      "330\n",
      "adding uncorrelated noise 1\n",
      "adding uncorrelated noise 2\n",
      "correlating\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "87eed3d40f4b40df83fa2ab5b6020d0d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d1be442814714161b690aad25d956eba",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy import signal\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#==============================================================\n",
    "#------------------------------------------------\n",
    "# parameters to modify\n",
    "#------------------------------------------------\n",
    "\n",
    "sigduration =1*24*60*60          # signal duration in sec\n",
    "\n",
    "noiselevel = .1                  # noise level\n",
    "\n",
    "#-------------------------- source distribution anisotropy\n",
    "anis = np.arange(0,180)\n",
    "az_to_source = 2*np.pi*np.arange(0,180)/180\n",
    "\n",
    "for i in range (0,180):\n",
    "    anis[i] = 1\n",
    "\n",
    "#for i in range (10,20):\n",
    "#    anis[i] = 10\n",
    "#\n",
    "#\n",
    "#anis = np.sin(az_to_source + np.pi*10/180)\n",
    "#idx = anis<0\n",
    "#anis[idx]=0\n",
    "\n",
    "#anis = np.abs(np.cos(az_to_source))\n",
    "#==============================================================\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#------------------ plotting mode\n",
    "%matplotlib widget\n",
    "#----------------------------\n",
    "\n",
    "plt.close('all')\n",
    "\n",
    "    \n",
    "\n",
    "#------------------------------------------------\n",
    "# plotting anisotropy\n",
    "#........................\n",
    "f1 = plt.figure()\n",
    "ax = f1.add_subplot(111, projection='polar')\n",
    "ax.plot(az_to_source,anis,'.')\n",
    "ax.set_theta_zero_location('N')\n",
    "ax.set_theta_direction(-1)\n",
    "plt.title('source strength azimuthal distribution')\n",
    "f1.show()\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#------------------------------------------------\n",
    "# subroutines for generating impulsive sgnal and noise\n",
    "#------------------------------------------------\n",
    "\n",
    "#-------------------------- single Ricker pulse\n",
    "\n",
    "def one_ricker(t,delay,f):\n",
    "    a = np.pi*np.pi*f*f*(t-delay)**2\n",
    "    return (1-2*a)*np.exp(-a)\n",
    "\n",
    "#-------------------------- random signal as sum of many Ricker pulses\n",
    "dt = 1.\n",
    "f = .07             # ricker frequency\n",
    "n_rick = 12\n",
    "\n",
    "trick = np.arange(-n_rick,n_rick)*dt\n",
    "rick0 = one_ricker(trick,0,f)\n",
    "\n",
    "def fast_random_rickers(t,nsig,delay,f,lam,nt,dt):\n",
    "    a = 2*np.random.rand(nsig)-1\n",
    "    r = np.random.rand(nsig)\n",
    "    dts = -np.log(r)/lam\n",
    "    ts = np.cumsum(dts)-250\n",
    "    tmax = nt*dt\n",
    "    drange = 1/f\n",
    "    u = 0*t\n",
    "    for i in range (0, nsig):\n",
    "        delI = ts[i]+delay\n",
    "        if delI > tmax:\n",
    "            break\n",
    "        dif_t = np.around((t-delI)/dt)\n",
    "        i_dif_t = dif_t.astype(int)\n",
    "        idx = (i_dif_t>-n_rick) & (i_dif_t<n_rick)\n",
    "        ii = np.around((t[idx]-delI)/dt)\n",
    "        iii = ii.astype(int)+n_rick\n",
    "        u[idx] = u[idx] + a[i]*rick0[iii]\n",
    "    return u\n",
    "\n",
    "\n",
    "#------------------------------------------------\n",
    "# defining parmeters\n",
    "\n",
    "dist = 500\n",
    "vel = 5\n",
    "\n",
    "#dt = .5\n",
    "npt = np.int(sigduration/dt)\n",
    "\n",
    "nsig = np.int(npt/5)\n",
    "nsig2 = 3*nsig\n",
    "\n",
    "lam = .0003\n",
    "lam2= lam\n",
    "\n",
    "\n",
    "\n",
    "#------------------------------------------------\n",
    "delay1 = 150;\n",
    "\n",
    "#..............................\n",
    "\n",
    "t = np.arange(npt)\n",
    "t = dt*t\n",
    "\n",
    "tcc = np.arange(-(npt-1),npt)\n",
    "tcc = dt*tcc\n",
    "\n",
    "\n",
    "sig1 = 0.*t\n",
    "sig2 = 0.*t\n",
    "\n",
    "#------------------------------------------------\n",
    "# computing signals from many sources\n",
    "#..............................\n",
    "\n",
    "for i in range (0,180):\n",
    "    ang1 = i*2\n",
    "    if ( int(ang1/30) == ang1/30) :\n",
    "        print(ang1)\n",
    "    angr1 = np.pi*ang1/180\n",
    "    delay2 = delay1 + dist*np.sin(angr1)/vel\n",
    "#\n",
    "    seed = np.random.randint(100000,size=1)\n",
    "    np.random.seed(seed)\n",
    "    sig1 = sig1 + anis[i]*fast_random_rickers(t,nsig,delay1,f,lam,npt,dt)\n",
    "    np.random.seed(seed)\n",
    "    sig2 = sig2 + anis[i]*fast_random_rickers(t,nsig,delay2,f,lam,npt,dt)\n",
    "\n",
    "print('adding uncorrelated noise 1')    \n",
    "sig1 = sig1 + noiselevel*fast_random_rickers(t,nsig,delay1,f,lam2,npt,dt)\n",
    "print('adding uncorrelated noise 2')    \n",
    "sig2 = sig2 + noiselevel*fast_random_rickers(t,nsig,delay1,f,lam2,npt,dt)\n",
    "\n",
    "print('correlating')    \n",
    "cc = signal.correlate(sig2,sig1,'full')/np.sqrt(np.sum(sig2**2)*np.sum(sig1**2))\n",
    "\n",
    "#..............................\n",
    "\n",
    "\n",
    "#------------------------------------------------\n",
    "# plotting stacked cross-correlation\n",
    "#..............................\n",
    "\n",
    "#------------------------------------------------\n",
    "# plotting signals\n",
    "#..............................\n",
    "\n",
    "plt.figure()\n",
    "\n",
    "plt.subplot(3,1,1)\n",
    "plt.plot(t,sig1,'r')\n",
    "plt.title('signal at receiver 1')\n",
    "\n",
    "plt.subplot(3,1,2)\n",
    "plt.plot(t,sig2,'b')\n",
    "plt.title('signal at receiver 2')\n",
    "\n",
    "plt.subplot(3,1,3)\n",
    "plt.plot(tcc,cc,'k')\n",
    "plt.xlim(-300,300)\n",
    "plt.title('cross-correlations between receivers 1 and 2')\n",
    "plt.xlabel('time (s)')\n",
    "\n",
    "plt.subplots_adjust(top=0.92, bottom=0.12, left=0.10, right=0.95, hspace=0.55, wspace=0.35)\n",
    "\n",
    "plt.show()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "#------------------------------------------------\n",
    "# plotting source-station geometry\n",
    "#..............................\n",
    "\n",
    "plt.figure()\n",
    "\n",
    "plt.axes().set_aspect('equal')\n",
    "\n",
    "y = 1250\n",
    "x = 1500\n",
    "plt.plot(x,y,'r^')\n",
    "x = 1000\n",
    "plt.plot(x,y,'b^')\n",
    "\n",
    "\n",
    "#..............................\n",
    "for i in range (0, 180):\n",
    "    ang = i*2\n",
    "    angr = np.pi*ang/180\n",
    "    x = 1250 + 750*np.sin(angr)\n",
    "    y = 1250 + 750*np.cos(angr)\n",
    "    plt.plot(x,y,'k*')\n",
    "    \n",
    "plt.ylim(0,2500)\n",
    "plt.xlim(0,2500)\n",
    "\n",
    "plt.axhline(y=1250,color='k',linestyle=':',linewidth=1)\n",
    "plt.axvline(x=1250,color='k',linestyle=':',linewidth=1)\n",
    "\n",
    "plt.title('sources (star) and receivers (triangle) positions')\n",
    "plt.xlabel('distance (km)')\n",
    "plt.ylabel('distance (km)')\n",
    "\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}