etas_linear_density_figs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reevaluation of linear aftershock distributions\n",
    "\n",
    "Reevaluate the linear spatial distributions. Are we calculating the initial densities correctly? Probably. Looking back at:\n",
    " - ETASmf/yodabass.py.calcEqDists()\n",
    "\n",
    "The main thing, i think, is that in order to accuratel validate the initial rates, we have to carefully select only the aftershocks. We do this in calcEqDists() by selecting approximately N_omori events after the mainshock... so it gets a bit complicated.\n",
    "\n",
    "The older ETASmf, yodabass.py code suffers some Python2 and numpy.datetime64... and other incompatibilities. Some of the salient bits probaly need to be ported out. Probably the best strategy is to leave behind some of the helper functions, like the eqcatalog() class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "failed while loading urllib and/or urllib. maybe python 3.x?\n"
     ]
    }
   ],
   "source": [
    "#\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "\n",
    "import datetime as dtm\n",
    "import matplotlib.dates as mpd\n",
    "import pytz\n",
    "tzutc = pytz.timezone('UTC')\n",
    "\n",
    "#import operator\n",
    "import math\n",
    "import random\n",
    "import numpy\n",
    "import scipy\n",
    "import scipy.optimize as spo\n",
    "import itertools\n",
    "import sys\n",
    "#import scipy.optimize as spo\n",
    "import os\n",
    "import operator\n",
    "#from PIL import Image as ipp\n",
    "import multiprocessing as mpp\n",
    "#\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "#import functools\n",
    "#\n",
    "#import shapely.geometry as sgp\n",
    "#\n",
    "from mpl_toolkits.basemap import Basemap as Basemap\n",
    "from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas\n",
    "from geographiclib.geodesic import Geodesic as ggp\n",
    "#\n",
    "#\n",
    "#import ANSStools as atp\n",
    "from yodiipy import ANSStools as atp\n",
    "#import bindex\n",
    "import contours2kml\n",
    "import globalETAS as gep\n",
    "from eq_params import *\n",
    "#import roc_generic            # we'll eventually want to move to a new library of roc tools.\n",
    "import random\n",
    "import geopy\n",
    "#\n",
    "#colors_ =  mpl.rcParams['axes.color_cycle']\n",
    "colors_ = ['b', 'g', 'r', 'c', 'm', 'y', 'k']\t\t# make sure these are correct...\n",
    "#\n",
    "emc = {'lat':32.13, 'lon':-115.30, 'event_date':dtm.datetime(2010,4,4,22,40,41, tzinfo=pytz.timezone('UTC'))}\n",
    "nepal_epi_lon = 84.698\n",
    "nepal_epi_lat = 28.175\n",
    "\n",
    "#\n",
    "# set up dict objects for different earthquakes. maybe from eq_params.py ?\n",
    "param_keys = ['lat_center', 'lon_center', 'to_dt', 'mc', 'mc_auto', 'cat_len_plus']\n",
    "nepal_params = {key:val for key,val in zip(param_keys, [nepal_epi_lat, nepal_epi_lon, \n",
    "                                                        dtm.datetime(2015,4,30, tzinfo=pytz.timezone('UTC')),\n",
    "                                                        4.0, 4.5, 220])}\n",
    "emc_params =   {key:val for key,val in zip(param_keys, [emc['lat'], emc['lon'],\n",
    "                                                      emc['event_date'] + dtm.timedelta(days=5), 2.5, 5.5, 220])}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class LD_figure(object):\n",
    "    def __init__(self, lat_center=nepal_epi_lat, lon_center=nepal_epi_lon, d_lat_auto=.5, d_lon_auto=.5,\n",
    "                       to_dt=dtm.datetime(2015,4,30, tzinfo=pytz.timezone('UTC')), dt_0=30, mc=2.5,\n",
    "                 mc_auto=4.5, d_lambda=1.76, dm=1.0, Ds=[1.0, 2.0], cat_len_plus=250.):\n",
    "        '''\n",
    "        # we want a liner density vs ETAS estimated linear densit figure.\n",
    "        # first, get a catalog. find the biggest event in the catalog, plot linear distance from that event.\n",
    "        # then, figure out the ETAS paramaters and see what we get and show for the domain 1. < D < 2.\n",
    "        '''\n",
    "        #\n",
    "        #eqcat = atp.catfromANSS(lon=lons, lat=lats, minMag=cm, dates0=dates, Nmax=None, fout=None, rec_array=True)\n",
    "        #eqcat.sort(order='event_date')\n",
    "        # auto_cat(lon_center=None, lat_center=None, d_lat_0=.25, d_lon_0=.5, dt_0=10,  mc=2.5, mc_0=4.5, to_dt=None, catlen=5.0*365.0 range_factor=5., rec_array=True, **kwargs)\n",
    "        eq_cat_params = atp.auto_cat_params(lat_center=lat_center, lon_center=lon_center,\n",
    "                                            d_lat=d_lat_auto, d_lon=d_lon_auto, dt_0=dt_0, mc_0=mc_auto, \n",
    "                                            to_dt=to_dt, range_factor=5.)\n",
    "        print('eq_cat_params: ', eq_cat_params)\n",
    "        # r_val = mpd.num2date(mpd.datestr2num(date_in.astype(str)), tz=tz_out)\n",
    "        eqcat = atp.catfromANSS(lat=eq_cat_params['lat'], lon=eq_cat_params['lon'], \n",
    "                                minMag=mc, dates0=[dtm64_to_datetime(eq_cat_params['mainshock_date']), \n",
    "                                                   dtm64_to_datetime(eq_cat_params['mainshock_date']) + dtm.timedelta(days=cat_len_plus)],\n",
    "                                fout=None, rec_array=True)\n",
    "        #\n",
    "        # find the biggest event; dump everything before that.\n",
    "        #eqcat.sort(order='mag')   # ... except that we want the index too. so let's just spin it:\n",
    "        mainshock = eqcat[0]\n",
    "        mainshock_index = 0\n",
    "        for j,rw in enumerate(eqcat):\n",
    "            if rw['mag']>mainshock['mag']: mainshock\n",
    "            #\n",
    "        #\n",
    "        eqcat.sort(order='event_date')\n",
    "        #\n",
    "        # get a measured dm:\n",
    "        try:\n",
    "            m1 = sorted(eqcat['mag'])[-2]\n",
    "            dm_meas = mainshock['mag']-m1\n",
    "        except:\n",
    "            dm_meas=None\n",
    "        dm_meas = (dm_meas or dm)\n",
    "        #\n",
    "        # now, we need a distance formula. this can be using geographiclib or spherical dists. we could hook up some\n",
    "        # of the class structures in the various etas codes...\n",
    "        #dists = []\n",
    "        #for j,ev in enumerate(eqcat):\n",
    "        #    g1=ggp.WGS84.Inverse(self.mainshock['lat'], self.mainshock['lon'], rw['y'], rw['x'])\n",
    "        #    r=g1['s12']/1000.)\n",
    "        #\n",
    "        dists = [[j,ggp.WGS84.Inverse(mainshock['lat'], mainshock['lon'], ev['lat'], ev['lon'])['s12']/1000.] \n",
    "                 for j,ev in enumerate(eqcat)]\n",
    "        #\n",
    "        dists2 =[[k,j,r] for k, (j,r) in enumerate(sorted(dists, key=lambda rw: rw[1]))]\n",
    "        #\n",
    "        self.__dict__.update(locals())\n",
    "        #\n",
    "    def plot_mean_N_prime(self, fignum=0, d_lambda=None, dm=None, dm_meas=None, mc=None):\n",
    "        #\n",
    "        f_size = (12,6)\n",
    "        mainshock=self.mainshock\n",
    "        d_lambda = (d_lambda or self.d_lambda)\n",
    "        d_lambda = (d_lambda or 1.76)\n",
    "        \n",
    "        dm = (dm or self.dm)\n",
    "        mc = (mc or self.mc)\n",
    "        #\n",
    "        dm_meas = (dm_meas or self.dm_meas)\n",
    "        #\n",
    "        fg = plt.figure(figsize=f_size)\n",
    "        plt.clf()\n",
    "        ax1 = fg.add_axes([.1,.1,.35,.8])\n",
    "        ax2 = fg.add_axes([.5,.1,.35,.8])\n",
    "        #\n",
    "        ax1.plot(*zip(*[[r,k] for k,j,r in self.dists2]), '.-')\n",
    "        ax1.set_ylabel('Number of events, $N$')\n",
    "        ax1.set_xlabel('Distance $r$')\n",
    "        #\n",
    "        #plt.figure(figsize=f_size)\n",
    "        #plt.clf()\n",
    "        #ax2=plt.gca()\n",
    "        X,Y = zip(*[[r,k/r] for k,j,r in self.dists2 if r!=0])\n",
    "        ax2.plot(X,Y, marker='.', ls='-')\n",
    "        #\n",
    "        # mark the (approximate) rupture boundary:\n",
    "        ax2.plot(numpy.ones(2)*.5*10**(.5*mainshock['mag'] - d_lambda), [min(Y), max(Y)], ls='--', lw=2., marker='')\n",
    "        ax2.plot(numpy.ones(2)*1.0*10**(0.5*mainshock['mag'] - d_lambda), [min(Y), max(Y)], ls='-', lw=2., marker='')\n",
    "        #ax2.plot(numpy.ones(2)*(1.0*10**(.5*mainshock['mag'] - d_lambda) + 1.0*10**(.5*6.3 - d_lambda)), [min(Y), max(Y)], ls='-', lw=2., marker='')\n",
    "        #\n",
    "        Y_D2   = 10.**self.log_N_prime_max(m=None, dm=None, D=2.0, mc=mc)\n",
    "        Y_D1   = 10.**self.log_N_prime_max(m=None, dm=None, D=1.0, mc=mc)\n",
    "        Y_D2_b = 10.**self.log_N_prime_max(m=None, dm=dm_meas, D=2.0, mc=mc)\n",
    "        Y_D1_b = 10.**self.log_N_prime_max(m=None, dm=dm_meas, D=1.0, mc=mc)\n",
    "\n",
    "        ax2.plot([min(X), max(X)], [Y_D2, Y_D2], 'm-' )\n",
    "        ax2.plot([min(X), max(X)], [Y_D1, Y_D1], 'm--' )\n",
    "        ax2.plot([min(X), max(X)], [Y_D2_b, Y_D2_b], 'c-' )\n",
    "        ax2.plot([min(X), max(X)], [Y_D1_b, Y_D1_b], 'c--' )\n",
    "\n",
    "        #\n",
    "        #ax2.set_yscale('log')\n",
    "        #ax2.set_xscale('log')\n",
    "        #\n",
    "        #self.__dict__.update(locals())\n",
    "\n",
    "    def plot_mean_N_prime(self, fignum=0, d_lambda=None, dm=None, dm_meas=None, mc=None):\n",
    "        #\n",
    "        f_size = (12,6)\n",
    "        mainshock=self.mainshock\n",
    "        d_lambda = (d_lambda or self.d_lambda)\n",
    "        d_lambda = (d_lambda or 1.76)\n",
    "        \n",
    "        dm = (dm or self.dm)\n",
    "        mc = (mc or self.mc)\n",
    "        #\n",
    "        dm_meas = (dm_meas or self.dm_meas)\n",
    "        #\n",
    "        fg = plt.figure(figsize=f_size)\n",
    "        plt.clf()\n",
    "        ax1 = fg.add_axes([.1,.1,.35,.8])\n",
    "        ax2 = fg.add_axes([.5,.1,.35,.8])\n",
    "        #\n",
    "        ax1.plot(*zip(*[[r,k] for k,j,r in self.dists2]), '.-')\n",
    "        ax1.set_ylabel('Number of events, $N$')\n",
    "        ax1.set_xlabel('Distance $r$')\n",
    "        #\n",
    "        #plt.figure(figsize=f_size)\n",
    "        #plt.clf()\n",
    "        #ax2=plt.gca()\n",
    "        X,Y = zip(*[[r,k/r] for k,j,r in self.dists2 if r!=0])\n",
    "        ax2.plot(X,Y, marker='.', ls='-')\n",
    "        #\n",
    "        # mark the (approximate) rupture boundary:\n",
    "        ax2.plot(numpy.ones(2)*.5*10**(.5*mainshock['mag'] - d_lambda), [min(Y), max(Y)], ls='--', lw=2., marker='')\n",
    "        ax2.plot(numpy.ones(2)*1.0*10**(0.5*mainshock['mag'] - d_lambda), [min(Y), max(Y)], ls='-', lw=2., marker='')\n",
    "        #ax2.plot(numpy.ones(2)*(1.0*10**(.5*mainshock['mag'] - d_lambda) + 1.0*10**(.5*6.3 - d_lambda)), [min(Y), max(Y)], ls='-', lw=2., marker='')\n",
    "        #\n",
    "        Y_D2   = 10.**self.log_N_prime_max(m=None, dm=None, D=2.0, mc=mc)\n",
    "        Y_D1   = 10.**self.log_N_prime_max(m=None, dm=None, D=1.0, mc=mc)\n",
    "        Y_D2_b = 10.**self.log_N_prime_max(m=None, dm=dm_meas, D=2.0, mc=mc)\n",
    "        Y_D1_b = 10.**self.log_N_prime_max(m=None, dm=dm_meas, D=1.0, mc=mc)\n",
    "\n",
    "        ax2.plot([min(X), max(X)], [Y_D2, Y_D2], 'm-' )\n",
    "        ax2.plot([min(X), max(X)], [Y_D1, Y_D1], 'm--' )\n",
    "        ax2.plot([min(X), max(X)], [Y_D2_b, Y_D2_b], 'c-' )\n",
    "        ax2.plot([min(X), max(X)], [Y_D1_b, Y_D1_b], 'c--' )\n",
    "        \n",
    "        #ax2.set_yscale('log')\n",
    "        #ax2.set_xscale('log')\n",
    "\n",
    "    def plot_N_prime(self, fignum=0, d_lambda=None, dm=None, dm_meas=None, mc=None):\n",
    "        #\n",
    "        f_size = (12,6)\n",
    "        mainshock=self.mainshock\n",
    "        d_lambda = (d_lambda or self.d_lambda)\n",
    "        d_lambda = (d_lambda or 1.76)\n",
    "        \n",
    "        dm = (dm or self.dm)\n",
    "        mc = (mc or self.mc)\n",
    "        #\n",
    "        dm_meas = (dm_meas or self.dm_meas)\n",
    "        #\n",
    "        fg = plt.figure(figsize=f_size)\n",
    "        plt.clf()\n",
    "        ax1 = fg.add_axes([.1,.1,.35,.8])\n",
    "        ax2 = fg.add_axes([.5,.1,.35,.8])\n",
    "        #\n",
    "        N = len(self.dists2)\n",
    "        ax1.plot(*zip(*[[r,k] for k,j,r in self.dists2]), '.-')\n",
    "        #ax1.plot(*zip(*[[r,N-k] for k,j,r in self.dists2]), '.-')\n",
    "        ax1.set_ylabel('Number of events, $N$')\n",
    "        ax1.set_xlabel('Distance $r$')\n",
    "        #\n",
    "        #plt.figure(figsize=f_size)\n",
    "        #plt.clf()\n",
    "        #ax2=plt.gca()\n",
    "        X,Y = zip(*[[r,1./(r-self.dists2[j][2])] for k,j,r in self.dists2[1:] if (r-self.dists2[j][2]>0)])\n",
    "        ax2.plot(X,Y, marker='.', ls='-')\n",
    "        #\n",
    "        # mark the (approximate) rupture boundary:\n",
    "        ax2.plot(numpy.ones(2)*.5*10**(.5*mainshock['mag'] - d_lambda), [min(Y), max(Y)], ls='--', lw=2., marker='')\n",
    "        ax2.plot(numpy.ones(2)*1.0*10**(0.5*mainshock['mag'] - d_lambda), [min(Y), max(Y)], ls='-', lw=2., marker='')\n",
    "        #ax2.plot(numpy.ones(2)*(1.0*10**(.5*mainshock['mag'] - d_lambda) + 1.0*10**(.5*6.3 - d_lambda)), [min(Y), max(Y)], ls='-', lw=2., marker='')\n",
    "        #\n",
    "        Y_D2   = 10.**self.log_N_prime_max(m=None, dm=None, D=2.0, mc=mc)\n",
    "        Y_D1   = 10.**self.log_N_prime_max(m=None, dm=None, D=1.0, mc=mc)\n",
    "        Y_D2_b = 10.**self.log_N_prime_max(m=None, dm=dm_meas, D=2.0, mc=mc)\n",
    "        Y_D1_b = 10.**self.log_N_prime_max(m=None, dm=dm_meas, D=1.0, mc=mc)\n",
    "\n",
    "        ax2.plot([min(X), max(X)], [Y_D2, Y_D2], 'm-' )\n",
    "        ax2.plot([min(X), max(X)], [Y_D1, Y_D1], 'm--' )\n",
    "        ax2.plot([min(X), max(X)], [Y_D2_b, Y_D2_b], 'c-' )\n",
    "        ax2.plot([min(X), max(X)], [Y_D1_b, Y_D1_b], 'c--' )\n",
    "\n",
    "        #\n",
    "        ax2.set_yscale('log')\n",
    "        ax2.set_xscale('log')\n",
    "        #\n",
    "        #self.__dict__.update(locals())\n",
    "        \n",
    "    def GR_dist(self, fignum=None):\n",
    "        N = len(self.eqcat)\n",
    "        GR = [[m, N-j] for j,m in enumerate(sorted(self.eqcat['mag']))]\n",
    "        #\n",
    "        if fignum is not None:\n",
    "            plt.figure(fignum)\n",
    "            plt.clf()\n",
    "            ax=plt.gca()\n",
    "            ax.plot(*zip(*GR), ls='-', marker='.')\n",
    "            #\n",
    "            ax.set_yscale('log')\n",
    "    def log_N_prime_max(self, m=None, mc=None, dm=1.0, D=1.5, d_lambda=1.76):\n",
    "        m = (m or self.mainshock['mag'])\n",
    "        mc = (mc or self.mc)\n",
    "        d_lambda = (d_lambda or self.d_lambda)\n",
    "        dm = (dm or self.dm)\n",
    "        #print('***', m, mc, dm)\n",
    "        #\n",
    "        return numpy.log10(2) - (.5*m - d_lambda) + (2./(2.+D))*numpy.log10(1.+D/2.) + (D/(2.+D))*(m - dm - mc)\n",
    "\n",
    "# eventually, consolidate the generalized datetime-wrangler code.\n",
    "def dtm64_to_datetime(x, tz_out=pytz.timezone('UTC')):\n",
    "    if isinstance(x, dtm.datetime): return x\n",
    "    return  mpd.num2date(mpd.datestr2num(x.astype(str)), tz=tz_out)\n",
    "\n",
    "def log_N_prime_max(m, mc, dm=1.0, D=1.5, d_lambda=1.76):\n",
    "    return numpy.log10(2) - (.5*m - d_lambda) + (2./(2.+D))*numpy.log10(1.+D/2.) + (D/(2.+D))*(m - dm - mc)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data handle fetched...\n",
      " no file.\n"
     ]
    },
    {
     "ename": "IndexError",
     "evalue": "list index out of range",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-3-7a013ae678e0>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mA\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mLD_figure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mnepal_params\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mA\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot_mean_N_prime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-2-fa84b214d27d>\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, lat_center, lon_center, d_lat_auto, d_lon_auto, to_dt, dt_0, mc, mc_auto, d_lambda, dm, Ds, cat_len_plus)\u001b[0m\n\u001b[1;32m     14\u001b[0m         eq_cat_params = atp.auto_cat_params(lat_center=lat_center, lon_center=lon_center,\n\u001b[1;32m     15\u001b[0m                                             \u001b[0md_lat\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0md_lat_auto\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0md_lon\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0md_lon_auto\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdt_0\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdt_0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmc_0\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmc_auto\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 16\u001b[0;31m                                             to_dt=to_dt, range_factor=5.)\n\u001b[0m\u001b[1;32m     17\u001b[0m         \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'eq_cat_params: '\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meq_cat_params\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     18\u001b[0m         \u001b[0;31m# r_val = mpd.num2date(mpd.datestr2num(date_in.astype(str)), tz=tz_out)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/home/myoder/Dropbox/Research/globalETAS/yodiipy/ANSStools.py\u001b[0m in \u001b[0;36mauto_cat_params\u001b[0;34m(lon_center, lat_center, d_lat_0, d_lon_0, dt_0, mc_0, to_dt, range_factor, **kwargs)\u001b[0m\n\u001b[1;32m    593\u001b[0m         \u001b[0;31m#\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    594\u001b[0m         \u001b[0;31m#biggest_earthquake = filter(lambda x: x['mag']==max(cat_0['mag']), cat_0)[0]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 595\u001b[0;31m         \u001b[0mmainshock\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0mcat_0\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnames\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mx\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mj\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mx\u001b[0m \u001b[0;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'mag'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcat_0\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'mag'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcat_0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\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[0m\u001b[1;32m    596\u001b[0m         \u001b[0;31m#\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    597\u001b[0m         \u001b[0;31m# now, get new map domain based on rupture length, etc.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: list index out of range"
     ]
    }
   ],
   "source": [
    "A = LD_figure(**nepal_params)\n",
    "A.plot_mean_N_prime()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A.plot_N_prime()\n",
    "# lats=[28.175-1., 28.175+1.], lons = [84.7-1., 84.7+1], to_dt=dtm.datetime(2015,5,1, tzinfo=pytz.timezone('UTC')) \n",
    "#\n",
    "A = LD_figure(**nepal_params)\n",
    "A.plot_mean_N_prime()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "A = LD_figure(**emc_params)\n",
    "\n",
    "A.plot_mean_N_prime()\n",
    "A.plot_N_prime()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A.GR_dist(fignum=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "***  -0.021050267773593537 -0.002108153484397598 -0.0002108184790551193\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmUAAAHVCAYAAACwrYzPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4nFd5/vH7mRnt++rdlnfHjmMn\ncZyVJGQPIQv7EoKBhNCWLhRoWUpbKNBCC4W2UNoQQxMIJCEs8Q8IIXHiLECcONiOV3mNLdmyFlub\ntWvm/P6YkbzNjCSPNO9I8/1cl67Z3tE80WSk2+ec9znmnBMAAAC85fO6AAAAABDKAAAAUgKhDAAA\nIAUQygAAAFIAoQwAACAFEMoAAABSAKEMAAAgBRDKAAAAUgChDAAAIAUEvC7gbJSXl7uqqiqvywAA\nABjSq6++2uScqxjquHEZyqqqqrRhwwavywAAABiSmR0YznFMXwIAAKQAQhkAAEAKIJQBAACkAEIZ\nAABACiCUAQAApABCGQAAQAoglAEAAKQAQhkAAEAKIJQBAACkAEIZAABACiCUAQAApABCGQAAQAog\nlAEAAKQAQhkAAEAKGJVQZmY3mVm1me0xs09HeTzLzB6JPL7ezKpOeuwzkfurzezG0agHAABgvEk4\nlJmZX9K3Jd0sabGk95jZ4tMOu1tSs3NunqRvSPpq5LmLJb1b0hJJN0n678j3AwAASCujMVK2UtIe\n59w+51yvpIcl3X7aMbdLeiBy/TFJ15qZRe5/2DnX45zbL2lP5Pt5qrsvqNauPq/LAAAAaWQ0Qtk0\nSTUn3a6N3Bf1GOdcv6RWSWXDfG7SfeWJnbriq894XQYAAEgjoxHKLMp9bpjHDOe54W9gdq+ZbTCz\nDY2NjSMscWR8ZjGqAAAAGBujEcpqJc046fZ0SYdjHWNmAUlFko4N87mSJOfcfc65Fc65FRUVFaNQ\ndmxmUsiRygAAQPKMRih7RdJ8M5ttZpkKL9xfc9oxayStilx/u6RnnHMucv+7I2dnzpY0X9LLo1BT\nQnwMlAEAgCQLJPoNnHP9Zvbnkp6U5Jf0PefcNjP7J0kbnHNrJK2W9AMz26PwCNm7I8/dZmaPStou\nqV/SR51zwURrSpSZMVIGAACSKuFQJknOuV9L+vVp9/3DSde7Jb0jxnO/LOnLo1HHaDGTyGQAACCZ\n6OgfhckIZQAAIKkIZVGE15SRygAAQPIQyqIIn33pdRUAACCdEMqi8JnJMX8JAACSiFAWhYmRMgAA\nkFyEsijC23KK0TIAAJA0hLIoIpmMMzABAEDSEMqi8A2MlHlcBwAASB+EsigGdkmnqz8AAEgWQlkU\nPt/AmjKPCwEAAGmDUBYHI2UAACBZCGVRDKwpAwAASBZCWRQDmYyRMgAAkCyEsih8tMQAAABJRiiL\nwiLnXzJSBgAAkoVQFsVg81hvywAAAGmEUBbFYPPYkMeFAACAtEEoi4KF/gAAINkIZVGwzRIAAEg2\nQlkUjJQBAIBkI5RFYcY2SwAAILkIZVEM9PN3pDIAAJAkhLIoWFMGAACSjVAWBWvKAABAshHKomCb\nJQAAkGyEsijYZgkAACQboSwKY6QMAAAkGaEsClpiAACAZCOURTG4pozzLwEAQJIQyqI4cfalt3UA\nAID0QSiLYrBPGfOXAAAgSQhlcTBSBgAAkoVQFsXASBk9/QEAQLIQyqJgTRkAAEg2QlkUPlpiAACA\nJCOURTEweUlHfwAAkCyEsihoHgsAAJKNUBbFiTVlpDIAAJAchLIoTpx9CQAAkByEsihYUwYAAJKN\nUBaFL/JTIZMBAIBkIZRFYZGxMkbKAABAshDKoqB5LAAASDZCWRRsswQAAJKNUBYFI2UAACDZCGVR\nsM0SAABINkJZFLTEAAAAyUYoi4JtlgAAQLIRyqIYWFPmSGUAACBJCGVRDK4p87gOAACQPghlUbAh\nOQAASDZCWRS+welLb+sAAADpg1AWFdssAQCA5CKURTE4UuZtGQAAII0QyqI40RKDWAYAAJKDUBYF\na8oAAECyEcqisME1ZR4XAgAA0gahLAqaxwIAgGQjlEVxok+Zt3UAAID0QSiLYqCjP+dfAgCAZCGU\nRcFIGQAASDZCWRSDe18SygAAQJIQyqIYmLykoz8AAEgWQlkUg81jPa4DAACkD0JZFLTEAAAAyUYo\ni4I1ZQAAINkIZVGwpgwAACQboSwKRsoAAECyEcqiONGnjFQGAACSI6FQZmalZvaUme2OXJbEOG5V\n5JjdZrbqpPvXmVm1mW2KfFUmUs9oObHQ39s6AABA+kh0pOzTktY65+ZLWhu5fQozK5X0j5IulrRS\n0j+eFt7udM4tj3w1JFjPqBicvqQpBgAASJJEQ9ntkh6IXH9A0h1RjrlR0lPOuWPOuWZJT0m6KcHX\nHVNsswQAAJIt0VA2yTlXJ0mRy2jTj9Mk1Zx0uzZy34DvR6Yu/95scCfwM5jZvWa2wcw2NDY2Jlh2\nfCz0BwAAyRYY6gAze1rS5CgP/d0wXyNa0BqIO3c65w6ZWYGkn0q6S9KD0b6Jc+4+SfdJ0ooVK8Y0\nLtESAwAAJNuQocw5d12sx8ys3symOOfqzGyKpGhrwmolXX3S7emS1kW+96HIZbuZ/UjhNWdRQ1ky\nsc0SAABItkSnL9dIGjibcpWkx6Mc86SkG8ysJLLA/wZJT5pZwMzKJcnMMiS9WdLWBOsZFWyzBAAA\nki3RUPYVSdeb2W5J10duy8xWmNn9kuScOybpi5JeiXz9U+S+LIXD2WuSNkk6JOm7CdYzKlhTBgAA\nkm3I6ct4nHNHJV0b5f4Nku456fb3JH3vtGM6JF2YyOuPFdaUAQCAZKOjfxSMlAEAgGQjlEXDNksA\nACDJCGVR+GJ2SwMAABgbhLIoBlpiMFIGAACShVAWhY8NyQEAQJIRyqIwDYyUeVwIAABIG4SyKAab\nx9LTHwAAJAmhLApj+hIAACQZoSyKE33KSGUAACA5CGVRnOjo72kZAAAgjRDKoqCjPwAASDZCWRRG\nR38AAJBkhLIoBprHEskAAECyEMpiMGOhPwAASB5CWQw+M9aUAQCApCGUxWBiTRkAAEgeQlkMPjPW\nlAEAgKQJeF1AqjJjpAwAgPEqFHJq7uxVfVuP6tu71dDWHb7e1q2G9h41tPeosa1b3//gSi2cXOB1\nuZIIZTGFF/p7XQUAADiZc05tXf2qb+9W/clBa+B6e7ca2nrU0N6tvuCZf8hLcjNUWZCtysIsza0o\nU2YgdSYNCWUxhBf6k8oAAEiWjp7+waDVcEbo6hkMYt19oTOeW5gd0KTCbE0qzNbFc/LC1wuyNKkw\nW5WF2ZpUmKXKguyUCmGnI5TFEF7o73UVAACMf845NXf2qa61S0dau1XX2q261i7VtXbrSGu3jkRC\n1/Ge/jOem5Ph1+SibFUWZGnZ9GJNKswaDF+TTgpbOZl+D/7LRhehLAZaYgAAMLRQyOlYZ6+OtHbr\ncEuXjrR1D4atk0NYT/+po1t+n2lSQZYmF2Vr0eQCXbWgYjBkTSo4MbqVnxUYbOo+0RHKYmGhPwAg\nzYVCTk0dPaprGQhaXapriwSulm7VtXWpvrVHvcFTA1eG3zSpMFtTirK1dHqxbliSrcmF2ZpanK3J\nRTmaUpSt8vws+X3pEbaGi1AWgy9NUjkAIH119wV1qKVLh1u6dKg5fFk7cLslPMp1+mL5TL9Pk4uy\nNbkoWxfOLBkMWZOLsgcvy/Oy5CNwjRihLAZaYgAAxrOBdVyHmsMB65Tw1Rq+PNrRe8pzfKbIiFaO\nLphZoqnFOZpalK0pRTmDoas0LzNtphOTjVAWA2vKAACprD8Y0pG2btU2nxq2Bm4fbulWV1/wlOfk\nZPg1rSRHU4tztGRqkaYVZ4dvF+VoWkmOJhVmK8OfumcnTnSEshjYZgkA4CXnnBqP96jmWJdqmztV\nc6xTNce6VNPcqZrmTtW1dKv/tDYB5fmZmlacowWTCvTGhZWaWhwOW9OKw1/FuRmMcqUwQlkMZkZL\nDADAmGrt6lPNsc5I6IoErmOdqmkOB7HT+3GV52dqekmuls8o0a3n5WhGaa6mR0LX1OIcZWeM/7YQ\n6YxQFkN4fSKpDABw9nr7Q6pt7tSBYwMjXaeGr7buU/tyFWQHNKMkV3Mr8nT1ggrNKM3VjNIcTS8J\nh6/cTP5sT2S8uzH4zBQ6s2EwAACn6OoN6sCxDh042qkDRwcuO/X60Q4dbuk6ZdYlK+DT9JLwCNcF\nM0s0ozRHM0pyw+GrJFdFuRne/YfAc4SyGHycfQkAiGjt6jspcHXo9aOdOhgJXg3tPaccW5yboVll\nebpgZoneev40zSzL06yyXM0qzVV5Pq0iEBuhLAYzU5BQBgBpo7WrT/saj2t/U4debwoHrwPHwiGs\npbPvlGMrC7JUVZanKxdUqKosVzPL8lRVlqtZpXmMduGsEcpi8PtoiQEAE01Pf1AHj3ZqX1OH9jV2\naH/T8chlxyk9u3wmTS3OUVVZnm5ZOiU80hUZ8ZpZmsvaLowJ/q+KgelLABifQiGnI23dJ0LXYADr\nUG1z5ylrvCoKsjS7PE/XL56kORV5ml2er9nleZpZmqvMAP26kFyEshh8tMQAgJTW1RvU3sbj2tMQ\n/trf1KF9TeEgdnIridxMv2aX52nZjGLdcf40zSnP05yKPFWV56kwm6lGpA5CWQxsswQAqaG1s097\nGtsHw9fuyOWhlq7BZSZ+n2lmaa5ml+fpsrllkVGvPM2tyFdlQRYNUzEuEMpiCLfEIJQBQDI459TY\n3hMOXpHRr9314euNJ53dmBXwaU5Fvi6YWaJ3rpiheZX5ml+Zr1lleUw3YtwjlMXg9xkjZQAwygbC\nV3V9u6qPtGtX/YkRsJMbqRZkBTS3Ml9XL6jQvMr8SPgq0LSSHPlpKYEJilAWA9ssAUBi2rr7tLu+\nXTuPtGvXkchlfbuaT2ovUZaXqXmV+bp12VTNr8zXvMoCzZ/ElCPSE6EsBp+F/0UHAIivpz+ovQ0d\n2hUJYNVH2rSrPrzma0Bepl8LJhfopnMna8GkAi2cXKCFkwpUlp/lYeVAaiGUxcDZlwBwKufCrSa2\nH27T9sNt2hmZgtzf1KFg5Bdmht80tyJfK6pK9N5JM7VocoEWTCrQ9JIcRr6AIRDKYvCZBn/JAEC6\n6QuGtLfx+GAA214X/jq5s/3M0lwtmFSgm5ZMDo98TS7Q7PI8ZfhZcA+cDUJZDD4W+gNIE23dfdpx\nuE076k6Er11Hjqs3GO71lRXwadHkAt187mQtnlKoxVMLtXByofKz+BMCjCY+UTH4jG2WAEw8De3d\n2nqoVVtq27S9rlXb69pUc+zE2q+yvEwtnlqoD15epcVTC7V4SqFml+cpwOgXMOYIZTGwzRKA8a6x\nvUdbD7XqtdpWbTnUqi2HWlTfFu75ZSbNLsvTedOL9e6LZg4GMM56BLxDKIsh3BKDUAZgfGg63qMt\nh1q1tbZVrx1q1dZDrapr7ZYUDmBzyvN06ZwyLZ1erKXTirR4KtOPQKrhExmDz6RQaOjjACDZWrv6\ntLmmRa/VtoRHwGpbdTgSwKRwAFs5u1RLpxVp6bQiLZlWRAADxgE+pTH4fab+IKkMgLf6giFVH2nX\nxpoWbTzYrE01LdrX2DH4+OzyPK2oigSw6UVaMrVQBWyyDYxLhLIYfExfAkgy55wOt3Zr08ETAWzL\noVb19If/gVien6nlM4r11vOn6fyZJVo6vUiFBDBgwiCUxcA2SwDGWmdvvzbVtGhTTYs2HgxfDmy+\nnRnw6dyphbrz4llaPrNY588opgErMMERymJgmyUAo62+rVsbXm/WhgPH9OqBZm073DbYpHp2eZ6u\nmFeu5TOKtXxGsc6ZUqjMAG0ogHRCKIvBZ6YgoQzAWQqFnHY3HA8HsNeb9cqBY4P9wLIzfFo2vVh/\netVcXTirRMtnFKskL9PjigF4jVAWg8+Msy8BDFt3X1Cba1q04UCzNrweHglr6+6XJJXnZ2nFrBKt\nurRKK6pKtZhRMABREMpioHksgHi6eoP648FmvbTvqF7ad1SbalrUFwz/zphXma9bzpuiC2eVasWs\nEs0qy2UtGIAhEcpiYJslACfr7O3XqweatX7fMb2076g214ZDmN9nOndakT50+WxdVFWqC2eVMBUJ\n4KwQymLw+cSaMiCNdfSEQ9hL+45q/f5j2lzTov5QOISdN71Id18xR5fMKdWKqlIaswIYFfwmiYFt\nloD00tsf0saDzfrdnia9uKdJr9W2qj/kFIiEsHuvnKOL55TpwlklhDAAY4LfLDH4mb4EJjTnnKrr\n2/Xi7nAIW7/vmLr6goMjYR+5ao4uiYSw3Ex+VQIYe/ymiYGF/sDEc7ilSy/uadLv9jTpd3uOqul4\nuFHr3Io8vXPFdF0+r1yXzC2jSz4ATxDKYmCbJWD86+zt1+/3HNXzuxv14p6mwT0jy/OzdMW8Ml0+\nr1yXzyvX1OIcjysFAEJZTEafMmDccc5pb+Nxratu1LrqRr28/5h6gyHlZvp18exSvXflTL1hfoUW\nTMqnRQWAlEMoi8HvY/oSGA8GRsOerW7Qc7saVdsc7po/vzJfqy6bpasXVmpFVYmyAn6PKwWA+Ahl\nMTB9CaSm8GhYh9ZVN5wxGnbZ3HL9yVVzdfXCCk0vyfW6VAAYEUJZDOGWGF5XAUCS+oMhvXqgWU/v\nqNfTOxq0vym8NozRMAATCaEsBp+F/0UOwBsdPf16flejntpRr2d3Nqi5s08ZftOlc8v1oStm642M\nhgGYYAhlMfjMFGSoDEiq+rZuPbW9Xk/vqNfv9xxVbzCkopwMXbOoUtcvnqQ3zC9XAe0qAExQhLIY\n/D6mL4FkeL2pQ09sPaLfbK3T5tpWSdLM0lzddeksXb94klbMKlHA7/O4SgAYe4SyGIzmscCY2dNw\nXE9sqdMTW49oe12bJGnZ9CL9zY0LdcPiSZpXScsKAOknoVBmZqWSHpFUJel1Se90zjVHOe43ki6R\n9KJz7s0n3T9b0sOSSiX9UdJdzrneRGoaLT62WQJGzcCWRk9sOaInttZpV/1xSdKFs0r0uVvO0U3n\nTmZ9GIC0l+hI2aclrXXOfcXMPh25/akox/2bpFxJHznt/q9K+oZz7mEz+x9Jd0v6ToI1jQq2WQIS\n45zTjrp2/WrLYT2x5Yj2NXXIZ9JFVaX6wm1LdOOSyZpclO11mQCQMhINZbdLujpy/QFJ6xQllDnn\n1prZ1SffZ+G5iWskvfek539eKRPKWOgPnI2DRzu1ZvMhPb7psHY3HJffZ7psbpnufsNs3bB4sioK\nsrwuEQBSUqKhbJJzrk6SnHN1ZlY5gueWSWpxzvVHbtdKmhbrYDO7V9K9kjRz5syzLHf4fD6mL4Hh\namjr1i9fq9Pjmw9rc02LJGllVam+dMe5etPSKSrNy/S4QgBIfUOGMjN7WtLkKA/9XYKvHW0Vb8wY\n5Jy7T9J9krRixYoxj0tMXwLxtXb26Tfb6rRm82H9Ye9RhZy0ZGqhPnPzIt26bCqbfAPACA0Zypxz\n18V6zMzqzWxKZJRsiqSGEbx2k6RiMwtERsumSzo8guePKbZZAs7UFwzp+V2NeuzVWq3d0aDeYEhV\nZbn682vm67ZlUzWvMt/rEgFg3Ep0+nKNpFWSvhK5fHy4T3TOOTN7VtLbFT4Dc0TPH2tsswScsPNI\nmx7bUKtfbDqspuM9KsvL1J2XzNRbzp+mpdOKaF8BAKMg0VD2FUmPmtndkg5KeockmdkKSX/inLsn\ncvsFSYsk5ZtZraS7nXNPKnxSwMNm9iVJGyWtTrCeUeOL/I0JhZx8Pv7gIP0c6+jVmk2H9Ngfa7X1\nUJsy/KZrFlXq7RfO0NULK5RBQ1cAGFUJhTLn3FFJ10a5f4Oke066/YYYz98naWUiNYwVf+Rf/iHn\n5Iu6/A2YePqDIT23q1E/2VCrtTvr1Rd0OndaoT5/62LdtnwaC/YBYAzR0T+GgdExpjCRDg63dOmR\nV2r06IYa1bV2qzw/U6surdLbLpyuc6YUel0eAKQFQlkMA0tkWOyPiao/GNKz1Y368csHta66QU7S\nlfMr9I+3LtG151QyPQkASUYoi8EXSWVkMkw0h1q69MjLB/XIhhrVt/WosiBLH33jPL1zxQzNKGWr\nIwDwCqEshoG1/UFSGSaAUMjphT1NeuD3r+vZ6nDnmqsWVOiLt8/UNYsqFWBUDAA8RyiLwXfSQn9g\nvDre06+fvlqrB37/uvY1dag8P0t//sZ5etdFM9gAHABSDKEshsHpy5DHhQBnYV/jcT34hwN67NVa\nHe/p1/IZxfrmu5br5qWTlRXwe10eACAKQlkMPhb6Y5wJhZye29Wo//v963puV6My/KY3nzdVqy6r\n0vIZxV6XBwAYAqEshhMtMQhlSG3dfUE9vumQvvvCfu1pOK7Kgix9/PoFes/KmaooyPK6PADAMBHK\nYhjYNoaF/khVLZ29emj9QX3/d6+r6XiPFk8p1DfftVxvWjpFmQEW7gPAeEMoi8FPSwykqJpjnVr9\n4n49uqFGnb1BXbWgQvdeOUeXzS1jD0oAGMcIZTGwpgypZtvhVn1n3V79ekud/D7Tbcum6cNXztai\nyXTcB4CJgFAWw4mWGB4XgrS38WCzvvXMHq3d2aCCrIA+fOUcffCy2ZpclO11aQCAUUQoi2FwmyVS\nGTyyft9RfevZPXphd5OKczP0yRsW6K5Lq1SUk+F1aQCAMUAoi4HmsfCCc04v7mnSf63do5dfP6by\n/Cx99k2LdOfFs5SXxccVACYyfsvH4PcxfYnkcc7phd1N+vendmlTTYsmF2br87cu1rtXzlR2Bs1e\nASAdEMpiMBb6I0le3n9MX3uyWi+/fkzTinP0z29ZqrddOI3O+wCQZghlMQxus0QowxjZXNOirz+1\nS8/valRFQZb+6fYletdFMwhjAJCmCGUxcPYlxkr1kXZ9/bfV+u32epXkZuizb1qkuy6pUk4mYQwA\n0hmhLIaBPmVBUhlGyaGWLn3tyWr9YtMh5WcG9PHrF+iDl1epIJuzKQEAhLKY2PsSo6Wtu0/fWbdX\nq1/cL5P0kSvn6k+umqPi3EyvSwMApBBCWQw+tllCgvqCIf1o/UH9x9rdOtbRq7eeP02fuHGhphXn\neF0aACAFEcpiYJslnC3nnJ7cVq+v/man9jd16LK5Zfrsm87RudOKvC4NAJDCCGUxsNAfZ2NHXZs+\nv2ab1u8/pvmV+fr+By7S1Qsr2CgcADAkQlkMxkJ/jEBrZ5/+/alq/eClAyrKydCX7jhX775ohgJ+\nn9elAQDGCUJZDAMd/elThniCIadHN9To356sVktnr953ySx9/PoFLOIHAIwYoSwGpi8xlFcPNOvz\na7Zpy6FWrZxdqs/fukSLpxZ6XRYAYJwilMXANkuIpaWzV//y6516ZEONJhdm6z/fc75uPW8K68YA\nAAkhlMVwYqSMUIYw55we33RYX/zldrV09ekjV87RX147X3lZfIwAAInjr0kMg6Es5HEhSAkHjnbo\nc7/Yqhd2N2nZjGL94C1LmaoEAIwqQlkMAyfNMVKW3vqCIX33hX36j6d3K8Pv0xduW6L3XTJr8EQQ\nAABGC6EsBmP6Mu1tPdSqT/5ks3YeadeNSybp87ct0ZQiuvEDAMYGoSwGtllKX739IX372T369rN7\nVJKXqf+960LduGSy12UBACY4QlkMbLOUnrYfbtMnf7JZ2+vadMfyqfr8bUvoOQYASApCWQwDI2V0\n9E8PfcGQvrNur/5z7W4V52YwOgYASDpCWQy0xEgfexqO62OPbNTWQ226ddlUfeG2JSrNY3QMAJBc\nhLIYBs6uC9ISY8JyzulHLx/UF3+5XTkZfn3nzgt089IpXpcFAEhThLIYBkJZP43KJqTmjl596qev\n6bfb63XFvHJ9/Z3LNKkw2+uyAABpjFAWw0AoY/py4nlxd5M+8ZNNOtbRq7970zm6+4rZ8tF3DADg\nMUJZDH5j+nKi6QuG9LUnq/W/z+/T3Io8rV51kc6dVuR1WQAASCKUxeQb6OjP2ZcTQl1rl/7iRxu1\n4UCz3nvxTP39LYuVk+n3uiwAAAYRymIYXOjP9OW49/yuRn3skU3q6Qvqv95zvm5dNtXrkgAAOAOh\nLAY/fcrGvWDI6T/W7tZ/PbNbCyoL9N/vu0BzK/K9LgsAgKgIZTGw0H98azreo489vEkv7mnS2y6Y\nri/dcS7TlQCAlEYoi2GwJUaQUDbebKlt1b0/2KBjHb3617edp3esmD64wTwAAKmKUBaDj5Gycenx\nTYf0t4+9pvL8LP30Ty/j7EoAwLhBKIuBNWXjSzDk9K+/2an/fX6fVs4u1X/feYHK87O8LgsAgGEj\nlMXA2ZfjR2tnn/7y4Y16blej7rpklv7h1sXK8Pu8LgsAgBEhlMUwuCE5I2UpbX9Thz70f6+otrlT\n//LWpXrPyplelwQAwFkhlMXAhuSp75XXj+nDD26Qz0w/+vAluqiq1OuSAAA4a4SyGAa2QmT6MjWt\n2XxYn3x0s6aX5Oj7H7xIs8ryvC4JAICEEMpiMDP5faZgiKGyVOKc03+v26t/e7JaK6tKdd/7L1Rx\nbqbXZQEAkDBCWRx+M6YvU0hfMKS//8VWPfxKjW5fPlX/+vbzlBWgISwAYGIglMXh89GnLFV09Qb1\n0R/9Uc/sbNBfXDNPH79+AQ1hAQATCqEsjvBIGaHMa61dfbrngVe04UCzvvyWc3XnxbO8LgkAgFFH\nKIvD5yOUea2hvVvvX/2y9jYe17fec4FuOW+K1yUBADAmCGVx+H3G9KWHDh7t1F3fW6+Gth6tXnWR\nrlxQ4XVJAACMGUJZHAGfqZ+RMk9UH2nXXavXqzcY0kMfvlgXzCzxuiQAAMYUoSwOnxkd/T2w/XCb\n3rd6vTL8pkc/cqkWTCrwuiQAAMYcoSwOP2vKkm7roVa9b/V65WT49eMPX6KqcprCAgDSA7s2x+Ez\no6N/Er1W26L3fvcl5WUG9Mi9lxLIAABphVAWh9/H9GWybKpp0Z33r1dhToYevvcSzSzL9bokAACS\niunLOPw+U5BMNuY21bTorvvXqyQvUz++9xJNK87xuiQAAJKOUBYHI2Vjb/vhNr1/dTiQPfKRSzSl\niEAGAEhPTF/G4TdTPxuSj5nccDlVAAAcKklEQVQ9Dcd11+r1yssK6KF7LiaQAQDSGqEsjnBHf6+r\nmJhqjnXqffevl5npoXsu1oxS1pABANIboSwOPxuSj4kjrd167/0vqasvqB/es1JzKvK9LgkAAM8R\nyuJgQ/LR19zRqzvvf0nNHX168EMrtWhyodclAQCQEljoH4ePvS9HVVdvUHc/8Ipqmrv04IdWatmM\nYq9LAgAgZSQ0UmZmpWb2lJntjlxG3aDQzH5jZi1m9svT7v8/M9tvZpsiX8sTqWe0MVI2evqDIf3F\njzdqY02L/vPdy3XJnDKvSwIAIKUkOn35aUlrnXPzJa2N3I7m3yTdFeOxv3HOLY98bUqwnlHFNkuj\nwzmnf1izTU/vqNfnb12im86d4nVJAACknERD2e2SHohcf0DSHdEOcs6tldSe4GslHaFsdHzrmT36\n0fqD+tOr52rVZVVelwMAQEpKNJRNcs7VSVLksvIsvseXzew1M/uGmWXFOsjM7jWzDWa2obGx8Wzr\nHZFwR39CWSJ+sqFGX39ql956/jT97Y0LvS4HAICUNWQoM7OnzWxrlK/bR+H1PyNpkaSLJJVK+lSs\nA51z9znnVjjnVlRUVIzCSw/NZ3T0T8T6fUf12Z9v0RXzyvWVt50nM/O6JAAAUtaQZ186566L9ZiZ\n1ZvZFOdcnZlNkdQwkhcfGGWT1GNm35f0yZE8f6wxUnb2Dh7t1J/88FXNKM3Vt997gTIDdF8BACCe\nRP9SrpG0KnJ9laTHR/LkSJCThYdQ7pC0NcF6RpXP6Oh/Ntq7+3T3A68o5KTVqy5SUW6G1yUBAJDy\nEg1lX5F0vZntlnR95LbMbIWZ3T9wkJm9IOknkq41s1ozuzHy0ENmtkXSFknlkr6UYD2jyu8T05cj\nFAw5/eWPN2p/U4e+874LNLs8z+uSAAAYFxJqHuucOyrp2ij3b5B0z0m33xDj+dck8vpjLeDzsSH5\nCP3Lr3fo2epGffkt5+qyueVelwMAwLjBQp84wh39va5i/Pj5xlrd/+J+feCyKt158SyvywEAYFwh\nlMXhN9GnbJh21LXpMz/bootnl+pzt5zjdTkAAIw7hLI4fDSPHZbWrj79yQ9fVVFOhr713gsU8PO/\nFQAAI8WG5HH4jQ3JhxIKOX3i0U061NylRz5yiSoKYvb/BQAAcTCkEQfbLA3tv9ft0dM7GvT3b16s\nC2eVel0OAADjFqEsDr+PkbJ4XtzdpK8/tUt3LJ+q91/Kwn4AABJBKIvD7zP1M1IWVdPxHv31o5s0\ntyJf//zWpWyhBABAglhTFke4oz+h7HThdWSb1dbVpx/cvVK5mfxvBABAohgpi8PvY0PyaO5/cZ+e\n29Woz715sRZNLvS6HAAAJgRCWRxsSH6mzTUt+tffVOvGJZP0votnel0OAAATBqEsDp+Z2GXphPbu\nPv3lwxtVWZClr77tPNaRAQAwilgMFIffJ0bKTvKF/7ddNcc69chHLlVxbqbX5QAAMKEwUhaH3+dT\nMOTkCGZ6anu9Hnu1Vn969VxdVEU/MgAARhuhLA5/ZHou3df6Hz3eo8/87DWdM6VQf3XtAq/LAQBg\nQiKUxRHwh0NZfxovLHPO6XO/2Kq2rn59413LlBngfxkAAMYCf2Hj8PsiI2Xpm8m0ZvNhPbH1iP76\n+gW0vwAAYAwRyuIIREJZX5qmsiOt3fr7X2zVhbNKdO+Vc7wuBwCACY1QFsdAKAsG029R2cC0ZW8w\npK+/Y9ngqCEAABgbhLI4/P7wjycdR8qe2HpET++o1yeuX6iq8jyvywEAYMIjlMWRMTBSlmanX7Z2\n9ukfHt+mpdOK9MHLq7wuBwCAtEDz2DgGpuz602z68l+e2KHmzl498KGLFPCT2wEASAb+4sZxoiVG\n+oSyP+w9qodfqdGH3zBHS6YWeV0OAABpg1AWR8AX/vEE02RNWXdfUJ/9+RbNKsvVx66b73U5AACk\nFaYv4xhsiZEm05ffWbdX+5s69NA9Fys7w+91OQAApBVGyuIYWE+VDgv9Dx7t1Hee26vblk3V5fPK\nvS4HAIC0QyiL48RI2cSfvvynX25Ths/0d7ec43UpAACkJUJZHP40aYnxzM56Pb2jQX913XxNKsz2\nuhwAANISoSyOdDj7srsvqC/8v+2aW5GnD1w22+tyAABIWyz0j2Pg7MuJ3Kfsu8/v04Gjnfrh3Rcr\nM0BGBwDAK/wVjuPESNnEXFN2uKVL3163R29aOllXzGdxPwAAXiKUxRGY4B39v/bbaoWc9Nk3sbgf\nAACvEcriGNxmaQKuKdt6qFU/++Mhfejy2Zpekut1OQAApD1CWRwZE7RPmXNOX/rVdpXmZerP3jjX\n63IAAIAIZXGdGCmbWGvKntnZoJf2HdPHrpuvwuwMr8sBAAAilMWVMQHPvuwLhvTPv96hOeV5es/K\nmV6XAwAAIghlcfgn4NmXD79So72NHfr0zYsGp2cBAID3+KscR8YEW+jf2duv/3h6t1bOLtX1iyd5\nXQ4AADgJoSwO/wRrifHA7w+o6XiP/vbGhTIzr8sBAAAnIZTFMdjRfwKMlLV19+l/nturqxdWaEVV\nqdflAACA0xDK4hjo6B+cAGvKVr+wX61dffrkDQu9LgUAAERBKItjYPqyb5xPXzZ39Gr1i/t187mT\nde60Iq/LAQAAURDK4pgozWP/57m96ujt18evX+B1KQAAIAZCWRyRgTL1B8fv9GVDW7ce+MPresvy\naZo/qcDrcgAAQAyEsjjMTAGfjeuF/vc9v099Qae/um6+16UAAIA4CGVDCPht3E5fHj3eo4fWH9Tt\ny6dqVlme1+UAAIA4CGVDCPh843ah//d+t1/d/UH92dXzvC4FAAAMgVA2hPBI2fhbU9ba1acHf39A\nbzp3iuZV5ntdDgAAGAKhbAgBn6lvHE5fPvj719Xe068/e+Ncr0sBAADDQCgbQsDnU3CcTV929PRr\n9e/269pFlVoylb5kAACMB4SyIfh9pr5xNn350PoDauns00evYS0ZAADjBaFsCOPt7Mue/qDuf2G/\nLptbpgtmlnhdDgAAGCZC2RDGW5+yNZsOq6G9Rx+5irVkAACMJ4SyIQR8vnHT0d85p9Uv7tfCSQW6\ncn651+UAAIARIJQNYTxNX76wu0k7j7Tr7jfMlpl5XQ4AABgBQtkQAj4bN81jv/vCPlUUZOn25VO9\nLgUAAIwQoWwIAb9P/ePg7MudR9r0wu4mrbp0lrICfq/LAQAAI0QoG0KGf3yMlN3/wn5lZ/h058Wz\nvC4FAACcBULZEDL8PvWl+EL/hrZuPb7pkN5x4QyV5GV6XQ4AADgLhLIhZI6DUPbQ+oPqCzp96IrZ\nXpcCAADOEqFsCBl+n3r7UzeU9QVD+vHLB3XVggrNLs/zuhwAAHCWCGVDyAj4UnpN2ZPbjqihvUfv\nv5S1ZAAAjGeEsiFk+C2lR8oe/MMBTS/J0dULK70uBQAAJIBQNoSsQOquKdt5pE0v7z+m910yS34f\nzWIBABjPCGVDSOWzL3/whwPKDPj0zhUzvC4FAAAkiFA2hHAoS701ZW3dffr5xkO69bypKqUNBgAA\n4x6hbAgZfp96U3Ck7Gev1qqzN6hVl7HAHwCAiYBQNoRMv6kvGJJzqTNa5pzTj1+u0XnTi3Te9GKv\nywEAAKOAUDaEDL9Pzkn9odQJZZtrW1Vd3653XzTT61IAAMAoSSiUmVmpmT1lZrsjlyVRjlluZn8w\ns21m9pqZveukx2ab2frI8x8xs5RbHJURCP+IUmmx/yOvHFROhl+3LpvidSkAAGCUJDpS9mlJa51z\n8yWtjdw+Xaek9zvnlki6SdI3zWxgzu2rkr4ReX6zpLsTrGfUZfgjoaw/NUbKOnr6tWbTYd1y3hQV\nZGd4XQ4AABgliYay2yU9ELn+gKQ7Tj/AObfLObc7cv2wpAZJFWZmkq6R9Fi853stMzJSliqL/X+1\npU4dvUG96yLaYAAAMJEkGsomOefqJClyGbetvJmtlJQpaa+kMkktzrn+yMO1kqbFee69ZrbBzDY0\nNjYmWPbwZfrDTVlTZfry0VdqNKciTytmnTFTDAAAxrHAUAeY2dOSJkd56O9G8kJmNkXSDyStcs6F\nIiNlp4s5R+icu0/SfZK0YsWKpM0lDk5fpkAo29PQrg0HmvWZmxcp+o8PAACMV0OGMufcdbEeM7N6\nM5vinKuLhK6GGMcVSvqVpM85516K3N0kqdjMApHRsumSDo/4v2CMpVIoe3RDrQI+01svmO51KQAA\nYJQlOn25RtKqyPVVkh4//YDIGZU/l/Sgc+4nA/e7cOOvZyW9Pd7zvTYQyno83pS8LxjSz/5Yq2sW\nVaqiIMvTWgAAwOhLNJR9RdL1ZrZb0vWR2zKzFWZ2f+SYd0q6UtIHzGxT5Gt55LFPSfq4me1ReI3Z\n6gTrGXWZgYE1Zd6effni7iY1He/V2y9klAwAgIloyOnLeJxzRyVdG+X+DZLuiVz/oaQfxnj+Pkkr\nE6lhrKXK9OXPNh5ScW6Grl4Y91wKAAAwTtHRfwiZg33KvAtl7d19+u22I3rzeVMGW3QAAICJhb/w\nQ8hIgT5lv9l6RD39Ib3lfKYuAQCYqAhlQxgcKfNwTdkvNh3SrLJcXTCTzccBAJioCGVD8HpN2ZHW\nbv1+71HdsXwavckAAJjACGVDyPC4o//jmw7JOemO82NudgAAACYAQtkQvO5T9vONh3T+zGLNLs/z\n5PUBAEByEMqGMHC2oxcjZTuPtGnnkXa9hVEyAAAmPELZEDI8bInxy8118pl0y9IpSX9tAACQXISy\nIZwYKUvu2ZfOOf1qS50um1uusny2VQIAYKIjlA1hYKF/svuUba9r0/6mDt1yHqNkAACkA0LZEDI9\nWuj/q9fq5PeZblwyOamvCwAAvEEoG4KZKTPgU09/MGmveWLqskyleZlJe10AAOAdQtkwZAd86ulL\n3kjZtsNtOnC0U29m6hIAgLRBKBuG7Ay/uvuSN1L2y9fqFPCZbljM1CUAAOmCUDYMWRm+pK0pC09d\nHtbl88pVwtQlAABpg1A2DNmB5I2UbTnUqppjXZx1CQBAmiGUDUMyR8p+tSU8dXkjU5cAAKQVQtkw\nJGukzDmn326r16Vzy1SUmzHmrwcAAFIHoWwYkrXQf2/jce1v6tAN9CYDACDtEMqGISuQnOnLJ7fV\nS5KuP2fSmL8WAABILYSyYUjWSNlvtx3RshnFmlyUPeavBQAAUguhbBiyMnzqHuPmsXWtXdpc26ob\nFjNKBgBAOiKUDUNWwD/m05dPbw9PXbLXJQAA6YlQNgzZGT71jPH05W+312tORZ7mVeaP6esAAIDU\nRCgbhqyAX91juCF5a1ef/rD3KNsqAQCQxghlw5Cd4VNf0CkYcmPy/ddVN6g/5HTDEtaTAQCQrghl\nw5Cd4Zck9YzRaNlvt9WroiBLy6cXj8n3BwAAqY9QNgxZgfCPqWcMzsDs6Q9qXXWDrjtnknw+G/Xv\nDwAAxgdC2TAMjJSNxbqyV19vVkdvUNcuqhz17w0AAMYPQtkwZGeEf0xj0ats3a5GZfp9unRu2ah/\nbwAAMH4QyoYhKzB2a8rWVTfootklyssKjPr3BgAA4wehbBjGaqTscEuXdtUf11ULKkb1+wIAgPGH\nUDYM2ZGRstHe//K5XY2SpKsXsp4MAIB0RygbhpzMcCjr6h3dULauukFTi7I1ny7+AACkPULZMAys\n9+ro7R+179kXDOl3e47qqoUVMqMVBgAA6Y5QNgy5kZGyzp7RGyl79UCzjvf066oFTF0CAABC2bDk\nZY7+SNlzuxoV8Jkun0crDAAAQCgbltysyEjZKK4pW1fdqAtnlaggO2PUvicAABi/CGXDkOn3KeAz\ndfSMzkhZfVu3dtS1cdYlAAAYRCgbBjNTbqZ/1EbKTrTCoD8ZAAAII5QNU15WYNRGyp6rbtSkwiwt\nmlwwKt8PAACMf4SyYcrLCozKSFl/MKQXdjfqqgW0wgAAACcQyoYpL9Ov46MwUrappkVt3bTCAAAA\npyKUDVNuZkCdo9ASY111o/w+0xXzy0ehKgAAMFEQyoYpL8uvjlFoHvvcrkZdMLNYRTm0wgAAACcQ\nyoZpNEbKGtt7tOVQq65awFmXAADgVISyYcrL8qsjwYX+zw+2wmA9GQAAOBWhbJhyMwPqTHCh/3O7\nGlWen6XFUwpHqSoAADBREMqGKS8zPFIWCrmzen4w5PT87kZduaBcPh+tMAAAwKkIZcOUnx3elPz4\nWa4r21zbopbOPqYuAQBAVISyYRo4W7K1s++snv9cdaN8Jr1hHq0wAADAmQhlwzQQytq6zy6UrdvV\nqGUzilWSlzmaZQEAgAmCUDZMhQMjZV0jD2XHOnr1Wm2LrqaLPwAAiIFQNkyDI2VnEcpe2N0o56Sr\nFtKfDAAAREcoG6aiBEbK1lU3qjQvU+dNKxrtsgAAwARBKBumsw1loZDT87sadeV8WmEAAIDYCGXD\nlJ8VkM9GHsq2Hm7V0Y5epi4BAEBchLJhMjMV5mSorWtkfcrWVTfKTLpyPqEMAADERigbgaKcjBGP\nlD23q1HnTStSWX7WGFUFAAAmAkLZCIw0lLV09mrjwWZdtYBRMgAAEB+hbASKcjLU0tk77ONf2N2k\nEK0wAADAMBDKRqCiIEtNx4cfytZVN6o4N0PLZ5SMYVUAAGAiIJSNQEV+lhqP98g5N+SxoZDTc7sa\ndNWCCvlphQEAAIZAKBuBioIs9faH1NY99BmYWw+3qul4r964kK2VAADA0AhlI1BRED6DsrG9Z8hj\nn9nZEG6FwSJ/AAAwDISyEajIH34oe7a6UctnFKs0L3OsywIAABMAoWwEBkfKjscPZU3He/RabQtT\nlwAAYNgIZSNQPsyRsud3Nco5EcoAAMCwJRTKzKzUzJ4ys92RyzN6P5jZcjP7g5ltM7PXzOxdJz32\nf2a238w2Rb6WJ1LPWCvOzVB2hk91LV1xj3tmZ4PK87O0ZGphkioDAADjXaIjZZ+WtNY5N1/S2sjt\n03VKer9zbomkmyR908yKT3r8b5xzyyNfmxKsZ0yZmaaX5KqmuTPmMT39Qa2rbtS1iyrloxUGAAAY\npkRD2e2SHohcf0DSHacf4Jzb5ZzbHbl+WFKDpHF7SuLM0lwdPBZ7pOz3e47qeE+/blo6OYlVAQCA\n8S7RUDbJOVcnSZHLuIuozGylpExJe0+6+8uRac1vmFnMXbvN7F4z22BmGxobGxMs++zNLM1V7bHO\nmA1kn9hap4KsgC6bW5bkygAAwHg2ZCgzs6fNbGuUr9tH8kJmNkXSDyR90DkXitz9GUmLJF0kqVTS\np2I93zl3n3NuhXNuRUWFdwNt00ty1N7Tr5bOMzcm7w+G9NT2el1zTqWyAn4PqgMAAONVYKgDnHPX\nxXrMzOrNbIpzri4SuhpiHFco6VeSPuece+mk710XudpjZt+X9MkRVe+BWWV5kqR9TR268LQeZOv3\nH1NzZ59uPpepSwAAMDKJTl+ukbQqcn2VpMdPP8DMMiX9XNKDzrmfnPbYlMilKbwebWuC9Yy5c6YU\nSJK217Wd8dhP/1irgqyArlpAKwwAADAyiYayr0i63sx2S7o+cltmtsLM7o8c805JV0r6QJTWFw+Z\n2RZJWySVS/pSgvWMuWnFOSrKydD2w62n3N/R06/fbD2iNy+bopxMpi4BAMDIDDl9GY9z7qika6Pc\nv0HSPZHrP5T0wxjPvyaR1/eCmWnJ1EJtPXTqSNmvttSpszeot1843aPKAADAeEZH/7OwoqpU2w63\nqrmjV5LknNPqF/Zr4aQCXTDzjP65AAAAQyKUnYU3LqxQyEnP7Qq35li7o0HV9e368JVzFF4eBwAA\nMDKEsrNw3vRiTS/J0YN/eF3t3X364q+2a055nm5bNtXr0gAAwDhFKDsLfp/pz66epz8ebNEVX31W\ntc1d+vJbliozwI8TAACcnYQW+qez96ycoebOXr2076g+cFmVLqWDPwAASACh7CyZmT76xnn66Bvn\neV0KAACYAJhvAwAASAGEMgAAgBRAKAMAAEgBhDIAAIAUQCgDAABIAYQyAACAFEAoAwAASAGEMgAA\ngBRAKAMAAEgBhDIAAIAUQCgDAABIAYQyAACAFEAoAwAASAGEMgAAgBRAKAMAAEgBhDIAAIAUQCgD\nAABIAYQyAACAFGDOOa9rGDEza5R0YIxfplxS0xi/BkaO9yX18J6kJt6X1MN7knqS9Z7Mcs5VDHXQ\nuAxlyWBmG5xzK7yuA6fifUk9vCepifcl9fCepJ5Ue0+YvgQAAEgBhDIAAIAUQCiL7T6vC0BUvC+p\nh/ckNfG+pB7ek9STUu8Ja8oAAABSACNlAAAAKYBQBgAAkAIIZVGY2U1mVm1me8zs017Xk47MbIaZ\nPWtmO8xsm5n9VeT+UjN7ysx2Ry5LvK413ZiZ38w2mtkvI7dnm9n6yHvyiJllel1jujGzYjN7zMx2\nRj4zl/JZ8ZaZ/XXkd9dWM/uxmWXzWUk+M/uemTWY2daT7ov62bCw/4z87X/NzC5Idr2EstOYmV/S\ntyXdLGmxpPeY2WJvq0pL/ZI+4Zw7R9Ilkj4aeR8+LWmtc26+pLWR20iuv5K046TbX5X0jch70izp\nbk+qSm//Iek3zrlFkpYp/P7wWfGImU2T9JeSVjjnzpXkl/Ru8Vnxwv9Juum0+2J9Nm6WND/yda+k\n7ySpxkGEsjOtlLTHObfPOdcr6WFJt3tcU9pxztU55/4Yud6u8B+ZaQq/Fw9EDntA0h3eVJiezGy6\npFsk3R+5bZKukfRY5BDekyQzs0JJV0paLUnOuV7nXIv4rHgtICnHzAKSciXVic9K0jnnnpd07LS7\nY302bpf0oAt7SVKxmU1JTqVhhLIzTZNUc9Lt2sh98IiZVUk6X9J6SZOcc3VSOLhJqvSusrT0TUl/\nKykUuV0mqcU51x+5zecl+eZIapT0/ci08v1mlic+K55xzh2S9DVJBxUOY62SXhWflVQR67Ph+d9/\nQtmZLMp99A3xiJnlS/qppI8559q8riedmdmbJTU45149+e4oh/J5Sa6ApAskfcc5d76kDjFV6anI\nGqXbJc2WNFVSnsJTY6fjs5JaPP99Rig7U62kGSfdni7psEe1pDUzy1A4kD3knPtZ5O76geHkyGWD\nV/Wlocsl3WZmrys8rX+NwiNnxZEpGonPixdqJdU659ZHbj+mcEjjs+Kd6yTtd841Ouf6JP1M0mXi\ns5IqYn02PP/7Tyg70yuS5kfOkslUeHHmGo9rSjuRtUqrJe1wzv37SQ+tkbQqcn2VpMeTXVu6cs59\nxjk33TlXpfDn4hnn3J2SnpX09shhvCdJ5pw7IqnGzBZG7rpW0nbxWfHSQUmXmFlu5HfZwHvCZyU1\nxPpsrJH0/shZmJdIah2Y5kwWOvpHYWZvUngEwC/pe865L3tcUtoxsyskvSBpi06sX/qswuvKHpU0\nU+FffO9wzp2+iBNjzMyulvRJ59ybzWyOwiNnpZI2Snqfc67Hy/rSjZktV/jki0xJ+yR9UOF/dPNZ\n8YiZfUHSuxQ+k3yjpHsUXp/EZyWJzOzHkq6WVC6pXtI/SvqFonw2IgH6Wwqfrdkp6YPOuQ1JrZdQ\nBgAA4D2mLwEAAFIAoQwAACAFEMoAAABSAKEMAAAgBRDKAAAAUgChDAAAIAUQygAAAFLA/wdmTcEL\nTRTGfAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fg = plt.figure(figsize=(10,8))\n",
    "ax = plt.gca()\n",
    "#\n",
    "X = numpy.linspace(0., 100., 10000)\n",
    "\n",
    "def f(r, r0=1., q=1.5):\n",
    "    return -r/(q*(r0 + r)**q)\n",
    "\n",
    "ax.plot(X, f(X, r0=1., q=1.5), '-')\n",
    "\n",
    "print('*** ', f(1000.), f(1e5), f(1e7))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}