{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Two-station method\n",
    "Here we compute the phase velocity form seismograms of two stations lying on a common great circle. We make use of the fact that the phase of the Fourier tranform of the cross-correlation of the two seismograms is equal to the phase difference between the surface wave trains from which we calclate the phase velocity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Preparation: load packages, set some basic options  \n",
    "%matplotlib inline\n",
    "from obspy import UTCDateTime\n",
    "from obspy.clients.fdsn import Client\n",
    "from obspy.geodetics import locations2degrees\n",
    "import numpy as np\n",
    "import matplotlib.pylab as plt\n",
    "from multipleFilter import multipleFilterAnalysis, normalizeMFT \n",
    "plt.style.use('ggplot')\n",
    "plt.rcParams['figure.figsize'] = 15, 4\n",
    "plt.rcParams['lines.linewidth'] = 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Get data from earthquake (Gulf of Alaska, m=7.9, d = 24 km)\n",
    "client = Client(\"BGR\")\n",
    "t1 = UTCDateTime(\"2018-01-23T09:31:42.000\")\n",
    "st1raw = client.get_waveforms(\"GR\",\"TNS\",\"\",\"LHZ\",t1,t1 + 2 * 3600,       # data of station TNS of GRSN\n",
    "                          attach_response = True)\n",
    "inv1 = client.get_stations(network=\"GR\",station=\"TNS\",\n",
    "                           channel='LHZ',level=\"channel\",starttime=t1)    # inventory info about TNS\n",
    "client = Client(\"GFZ\")\n",
    "st2raw = client.get_waveforms(\"GE\",\"STU\",\"\",\"LHZ\",t1,t1 + 2 * 3600,       # data of station STU of GEOFON\n",
    "                          attach_response = True)\n",
    "inv2 = client.get_stations(network=\"GE\",station=\"STU\",\n",
    "                           channel='LHZ',level=\"channel\",starttime=t1)    # inventory info about STU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "st1 = st1raw.copy()                                                 # copy data\n",
    "st2 = st2raw.copy()\n",
    "st1.detrend('linear')                                               # detrend\n",
    "st2.detrend('linear')\n",
    "st1.remove_response(output='VEL',zero_mean = True,                  # remove response (pre-filter is important)\n",
    "                    pre_filt = (0.001, 0.005, 0.4, 0.5))\n",
    "st2.remove_response(output='VEL',zero_mean = True,                  # remove response(pre-filter is important)\n",
    "                    pre_filt = (0.001, 0.005, 0.4, 0.5))\n",
    "st1.plot()                                                          # plot data\n",
    "st2.plot()\n",
    "print(st1)\n",
    "print(st2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "co = inv1.get_coordinates(\"GR.TNS..LHZ\")                                # extract coordinates from inventory\n",
    "tns = (co[\"longitude\"],co[\"latitude\"])\n",
    "co = inv2.get_coordinates(\"GE.STU..LHZ\")\n",
    "stu = (co[\"longitude\"],co[\"latitude\"])\n",
    "event = (-149.09,56.05)                                                 # coordinates of event\n",
    "dis = locations2degrees(tns[1],tns[0],stu[1],stu[0])*np.pi/180.*6371.   # epicentral distance\n",
    "print(\"Distance in km: \",dis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"../images/phase_velocity_map.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "tw1 = UTCDateTime(\"2018-01-23T10:05:00.000\")          # cut seismograms to surface wave train\n",
    "tw2 = UTCDateTime(\"2018-01-23T10:16:00.000\")\n",
    "st1.trim(tw1,tw2)\n",
    "st2.trim(tw1,tw2)\n",
    "st1.plot()\n",
    "st2.plot()\n",
    "npts = st1[0].stats.npts\n",
    "dt = st1[0].stats.delta\n",
    "print(\"Number of samples: \",npts,\"Sampling interval: \",dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "cc = np.correlate(st2[0],st1[0],mode=\"full\")          # Cross correlation: cc(k)=sum_n a(n+k)b(n)\n",
    "tm = np.arange(-(npts-1),npts)*dt                     # includes negative lags as well\n",
    "plt.plot(tm,cc,'-k')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ns = 256                                       # limit length of cross correlation to first ns points\n",
    "ccp = cc[npts-1:npts-1+ns]                     # cross correlation for positive lag\n",
    "print(len(ccp))\n",
    "tm = np.arange(0,ns)*dt                        # associated time values\n",
    "plt.plot(tm,ccp,'-k')\n",
    "plt.xlabel(\"Time\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "nd = int(pow(2, np.ceil(np.log(ns)/np.log(2))))          # Fourier transfrom of the cross correlation\n",
    "fcp = np.fft.rfft(ccp,nd)\n",
    "freq = np.fft.rfftfreq(nd,dt)\n",
    "print(nd,dt,len(freq),freq[0:4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fmax = 0.1                                               # compute phase up to maximum frequency fmax\n",
    "fny = 1./(2.*dt)\n",
    "jmax = int(fmax/fny*(nd//2-1))\n",
    "print(jmax)\n",
    "spec = fcp[0:jmax]\n",
    "frs = freq[0:jmax]\n",
    "phase = np.unwrap(np.angle(spec))                        # unwrap the phase owing to 2*pi multiples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "plt.subplot(311)\n",
    "plt.plot(frs,np.abs(spec),\"-k\")                          # plot abs of spectrum\n",
    "plt.subplot(312)\n",
    "plt.plot(frs,np.angle(spec),\"-k\")                        # plot phase of spectrum\n",
    "plt.subplot(313)\n",
    "plt.plot(frs,phase,\"-k\")                                 # plot unwrapped phase\n",
    "plt.xlabel(\"Frequency [Hz]\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "j1 = int(jmax*0.1)\n",
    "j2 = int(jmax*0.6)\n",
    "phavel = -2.*np.pi*frs[j1:j2]*dis/phase[j1:j2]                # calculate phase velocity and plot\n",
    "plt.plot(frs[j1:j2],phavel,'-k')\n",
    "plt.xlabel(\"Frequency [Hz]\")\n",
    "plt.ylabel(\"Phase velocity [km/s]\")\n",
    "plt.show()"
   ]
  },
  {
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}