diff --git a/Spark4Physicists.ipynb b/Spark4Physicists.ipynb
new file mode 100644
index 0000000..4c79200
--- /dev/null
+++ b/Spark4Physicists.ipynb
@@ -0,0 +1,745 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Analysing billion-objects catalogue interactively: Spark for physicists\n",
+    "<br> Authors: **S. Plaszczynski, J. Peloton, C. Arnault and J.E. Campagne **\n",
+    "<br> ***LAL, Univ. Paris-Sud, CNRS/IN2P3, Universite Paris-Saclay, Orsay, France***\n",
+    "<br> Last Verified to Run: **2019-02-07**\n",
+    "\n",
+    "\n",
+    "\n",
+    "This jupyter notebook describes the analysis presented in  https://arxiv.org/abs/1807.03078\n",
+    "The results presented here were run on the Universite Paris-Sud Spark cluster: https://www.informatique-scientifique.u-psud.fr/services/spark.html\n",
+    "\n",
+    "In order to run this notebook you need:\n",
+    "- a cluster running `pypark` (spark version >=2.2) on which the following python libraries have been installed (on each executor):\n",
+    "`matplotlib`,`numpy`,`pandas`,`pyarrow`,`healpy`\n",
+    "- a spark jupyter kernel (NERSC users: https://github.com/astrolabsoftware/spark-kernel-nersc)\n",
+    "- the `SparkFITS` library : https://astrolabsoftware.github.io/spark-fits/\n",
+    "- a set of FITS files, as produced for instance by the `CoLoRe` simulation : https://github.com/damonge/CoLoRe\n",
+    "The directory that points to these files should be put in the FITSDIR environment variable. Any file(s) with (RA,DEC,redshift) coordinates may also be used with little adaptation to match the exact key names.\n",
+    "\n",
+    "For testing, Spark (and these dependencies) can also be installed on a personal computer https://spark.apache.org/docs/latest/.\n",
+    "\n",
+    "For LSST/DESC collaborators with a NERSC account see https://github.com/LSSTDESC/desc-spark (FITSDIR=/global/cscratch1/sd/plaszczy/LSST10Y)\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Initialisation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Spark session started\n"
+     ]
+    }
+   ],
+   "source": [
+    "from pyspark.sql import SparkSession\n",
+    "from pyspark import StorageLevel\n",
+    "from pyspark.sql import functions as F\n",
+    "from pyspark.sql.functions import randn\n",
+    "from pyspark.sql.types import IntegerType,FloatType\n",
+    "from pyspark.sql.functions import pandas_udf, PandasUDFType\n",
+    "\n",
+    "spark = SparkSession.builder.getOrCreate()\n",
+    "sc=spark.sparkContext\n",
+    "print(\"Spark session started\")\n",
+    "\n",
+    "#a class for timing\n",
+    "from time import time\n",
+    "class Timer:\n",
+    "    \"\"\"\n",
+    "    a simple class for printing time (s) since last call\n",
+    "    \"\"\"\n",
+    "    def __init__(self):\n",
+    "        self.t0=time()\n",
+    "    \n",
+    "    def start(self):\n",
+    "        self.t0=time()\n",
+    "        \n",
+    "    def stop(self):\n",
+    "        t1=time()\n",
+    "        print(\"{:2.1f}s\".format(t1-self.t0))\n",
+    "\n",
+    "timer=Timer()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Reading the data \n",
+    "\n",
+    "all CoLoRe FITS files assumed to be in a directory pointed by the FITSDIR env variable\n",
+    "Le us build the **gal** datafreame selecting only the \"RA\" and \"dec\" columns and building a new redshift column called \"z\" on the flight."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "6.1s\n"
+     ]
+    }
+   ],
+   "source": [
+    "import os\n",
+    "dirfits=os.environ['FITSDIR']\n",
+    "timer.start()\n",
+    "gal=spark.read.format(\"fits\").option(\"hdu\",1)\\\n",
+    "         .load(dirfits)\\\n",
+    "         .select(F.col(\"RA\"), F.col(\"Dec\"), (F.col(\"Z_COSMO\")+F.col(\"DZ_RSD\")).alias(\"z\"))\n",
+    "timer.stop()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Pretty fast isn't it? This is because of lazy evaluation : data have not been actaully loaded (yet). The pipeline (Direct Acyclic Graph, DAG) has been updated but its execution has not been triggered since read/load/slect are not actions but transformations.\n",
+    "At this level the schema of the data is however available so you can looks at it"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "root\n",
+      " |-- RA: float (nullable = true)\n",
+      " |-- Dec: float (nullable = true)\n",
+      " |-- z: float (nullable = true)\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "gal.printSchema()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "you may want to look at a few values to check eveything is OK. Let's print the first 5"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "+---------+---------+---------+\n",
+      "|       RA|      Dec|        z|\n",
+      "+---------+---------+---------+\n",
+      "|225.80168|18.519966|2.4199903|\n",
+      "|225.73839|18.588171|2.4056022|\n",
+      "|225.79999|18.635067| 2.396816|\n",
+      "|225.49783|18.570776|2.4139786|\n",
+      "|225.57983|18.638515|2.3995044|\n",
+      "+---------+---------+---------+\n",
+      "only showing top 5 rows\n",
+      "\n",
+      "14.1s\n"
+     ]
+    }
+   ],
+   "source": [
+    "timer.start()\n",
+    "gal.show(5)\n",
+    "timer.stop()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Still pretty fast. Here show is an action and has triggered the data access.But not the full data set : since we required only the first 5 rows, the DAG was (automatically) analyzed and only the first block was actually read"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now we add an extra column that performs gaussian smearing on the \"z\" column. We use the Spark \"rand\" function that is highly optimized"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "+---------+---------+---------+---------+\n",
+      "|       RA|      Dec|        z|     zrec|\n",
+      "+---------+---------+---------+---------+\n",
+      "|225.80168|18.519966|2.4199903|2.3029215|\n",
+      "|225.73839|18.588171|2.4056022|2.6602228|\n",
+      "|225.79999|18.635067| 2.396816|2.3686476|\n",
+      "|225.49783|18.570776|2.4139786|2.4345846|\n",
+      "|225.57983|18.638515|2.3995044|2.3599765|\n",
+      "+---------+---------+---------+---------+\n",
+      "only showing top 5 rows\n",
+      "\n",
+      "8.3s\n"
+     ]
+    }
+   ],
+   "source": [
+    "timer.start()\n",
+    "gal=gal.withColumn(\"zrec\",(gal.z+0.03*(1+gal.z)*randn()).astype('float'))\n",
+    "gal.show(5)\n",
+    "timer.stop()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now we have all our columns ready, we put data in cache. Recall that the amount of memory available for the cache is about 60% of your total executors memory. You may simply use the cache() function. Here in case you do not have enough memory we use a slighlty more evolved form that fills the cache and spills the rest over disk if necessary.\n",
+    "Remember to finalize your command with an actiom as count() here.\n",
+    "This is the lenghty part but you need to perform it only once."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "#galaxies=5926764520\n",
+      "97.6s\n"
+     ]
+    }
+   ],
+   "source": [
+    "timer.start()\n",
+    "print(\"#galaxies={}\".format(gal.persist(StorageLevel.MEMORY_AND_DISK).count()))\n",
+    "timer.stop()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We put about 6 billions of galaxies in cache, just as if we had a hge RAM. All further computattions will be boosted."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Basic statistics"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Le us first look at some simple statistics on a single column"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "+-------+-------------------+\n",
+      "|summary|                  z|\n",
+      "+-------+-------------------+\n",
+      "|  count|         5926764520|\n",
+      "|   mean|  0.875229446764415|\n",
+      "| stddev|0.47360539212211805|\n",
+      "|    min|        -5.93947E-4|\n",
+      "|    max|          2.4352543|\n",
+      "+-------+-------------------+\n",
+      "\n",
+      "5.6s\n"
+     ]
+    }
+   ],
+   "source": [
+    "timer.start()\n",
+    "gal.describe(['z']).show()\n",
+    "timer.stop()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "On all the columns:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "+-------+------------------+--------------------+------------------+-------------------+\n",
+      "|summary|                RA|                 Dec|                 z|               zrec|\n",
+      "+-------+------------------+--------------------+------------------+-------------------+\n",
+      "|  count|        5926764520|          5926764520|        5926764520|         5926764520|\n",
+      "|   mean|179.99419917017775|-0.00314190695924...|0.8752294467644149| 0.8752288274953041|\n",
+      "| stddev|103.92702755739798|   39.17648956108155|0.4736053921221175|0.47714711978268487|\n",
+      "|    min|      4.9597503E-8|          -89.997185|       -5.93947E-4|       -0.124362685|\n",
+      "|    max|             360.0|             89.9986|         2.4352543|            2.95336|\n",
+      "+-------+------------------+--------------------+------------------+-------------------+\n",
+      "\n",
+      "9.2s\n"
+     ]
+    }
+   ],
+   "source": [
+    "timer.start()\n",
+    "gal.describe().show()\n",
+    "timer.stop()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If you just need some specific values it will be more efficient to call them directly"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "min=-0.0005939470138400793 max=2.4352543354034424\n",
+      "4.5s\n"
+     ]
+    }
+   ],
+   "source": [
+    "timer.start()\n",
+    "minmax=gal.select(F.min(\"z\"),F.max(\"z\")).first()\n",
+    "print(\"min={} max={}\".format(minmax[0],minmax[1]))\n",
+    "timer.stop()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Histograms"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The distributed way of building histograms is:\n",
+    "- add a column of bin values\n",
+    "- groupBy this column and count the number of elements in each group\n",
+    "- you may then drop the bin and add the center of the bin locations."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## with Spark native functions\n",
+    "Here is how to do it using Spark dataframe functions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#histogram caracteristics\n",
+    "zmin=minmax[0]\n",
+    "zmax=minmax[1]\n",
+    "Nbins=100\n",
+    "dz=(zmax-zmin)/Nbins\n",
+    "\n",
+    "#create the bin number column\n",
+    "zbin=gal.select(gal.z,((gal['z']-zmin-dz/2)/dz).cast(IntegerType()).alias('bin'))\n",
+    "# groupBy it + count + sort by bin values\n",
+    "h=zbin.groupBy(\"bin\").count().orderBy(F.asc(\"bin\"))\n",
+    "# add the bin locations (zbin) and drop the bin numer\n",
+    "h=h.select(\"bin\",(zmin+dz/2+h['bin']*dz).alias('loc'),\"count\").drop(\"bin\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now data is reduced to a simple histogram, we go back to the standard python world by converting the dataframe to a pandas object (this is the action so triggers the DAG materialization)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Index(['loc', 'count'], dtype='object') size=100\n",
+      "13.0s\n"
+     ]
+    }
+   ],
+   "source": [
+    "timer.start()\n",
+    "p=h.toPandas()\n",
+    "print(p.columns,\"size={}\".format(p.index.size))\n",
+    "timer.stop()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let us plot the result with matplotlib"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAERCAYAAACAbee5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEuZJREFUeJzt3W2MZOdZ5vH/xcR2JCcEwQwi8oynDRjBECU4aTmBwGJIFtn+YH8goLEIYOQw4sXhJQjkBeRE5hNBu2jZdchOwDJEEBMCggFNMIIYBQGO3CHB2I4MgzONW05w581ZZHDw7s2HrrbL7eqpM91Vdaqe+v+kVqrqnKm5z9Tk8j33OeepVBWSpLZ8Ud8FSJImz3CXpAYZ7pLUIMNdkhpkuEtSgwx3SWpQr+Ge5I4kjyd5oMO+lya5J8lHktyf5NpZ1ChJi6jvzv1O4OqO+/488N6qugI4DrxjWkVJ0qLrNdyr6oPAZ4ZfS/JVSf4kyYeT/GWSr93eHfjiweOXAI/NsFRJWigv6LuAEU4CP1RV/5jk1Wx16N8OvA340yRvBi4GXt9fiZI03+Yq3JO8CPgm4HeTbL980eB/bwDurKr/nuQbgXcneVlV/f8eSpWkuTZX4c7WmOhzVfUNI7bdxGA+X1V/k+SFwEHg8RnWJ0kLoe8Tqs9RVZ8HPp7kuwCy5RWDzf8MvG7w+tcBLwQ2eylUkuZc+lwVMsl7gKvY6sD/BXgr8AHgV4GXAhcAd1XVbUmOAe8CXsTWydWfqao/7aNuSZp3vYa7JGk65mosI0majN5OqB48eLBWVlb6+u0laSF9+MMf/lRVHRq3X2/hvrKywtraWl+/vSQtpCTrXfZzLCNJDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3BfcysoKSUjCxRdfPPLxzucu+yC1b96+rEMdrKyssL7+7B3I2yt7Jhn5eNQ2SW2zc19A6+vrVBV7Xa75oosusqOXGmfnvoSeeuqpXbt9SW2wc9czhjt6u3hpsdm5L4idc/Zp2NnRS1pcdu4LYr9zdknLxXCXpAYZ7nNs+Br2WXP+Li02Z+5zbHsUA7OfgTt/lxabnbskNchw11iOaKTF41hmzsziksfz5YhGWjx27nPGSx4lTYLhLkkNGhvuSe5I8niSB3bZ/j1J7h/8/HWSV0y+TM0L5+/SYujSud8JXH2O7R8HvrWqXg78AnByAnVpTm3P36tq7s4NSHrW2BOqVfXBJCvn2P7XQ0/vBQ7vvyxJ0n5MeuZ+E/D+3TYmOZFkLcna5ubmhH/rxdXnnaj74YhGml8TuxQyybexFe7fvNs+VXWSwdhmdXXVy0EG+rwTdT+8RFKaXxMJ9yQvB34NuKaqPj2J95Qk7d2+xzJJLgV+H/jeqvqH/ZckSdqvsZ17kvcAVwEHk2wAbwUuAKiqdwK3Al8GvGPwT/Onq2p1WgVLksbrcrXMDWO2vwl408QqWhLzuMzAfmyfXAU4evQoZ8+e7bcgacm5tkxPFvUk6m48uSrNF5cfkKQGGe6S1CDDXRPnzU1S/5y5a+Kcv0v9s3OfoUVdZkDS4jHcZ8gv4pA0K4a7psr5u9QPZ+6aKufvUj/s3CWpQYa7JDXIcNfMOH+XZseZu2bG+bs0O3buktQgw12SGmS4S1KDDPcpc8mB0Ty5Kk2XJ1SnrLUv5ZgUT65K02XnLkkNMtwlqUGGu3rn/F2aPGfu6p3zd2nyxnbuSe5I8niSB3bZniS/kuRMkvuTvHLyZUqSzkeXscydwNXn2H4NcPng5wTwq/svS5K0H2PDvao+CHzmHLtcD/xmbbkX+JIkL51UgYvIa9sl9W0SJ1QvAR4der4xeG1p+XV6kvo2iXAf1Z6OTLUkJ5KsJVnb3NycwG+t1njljDQZk7haZgM4MvT8MPDYqB2r6iRwEmB1ddW2Vs/jlTPSZEyicz8FfN/gqpnXAE9U1Scm8L6SpD0a27kneQ9wFXAwyQbwVuACgKp6J3AauBY4AzwJ/MC0ipUkdTM23KvqhjHbC/jRiVUkDWzP3wGOHj3K2bNn+y1IWiDeoaq55fxd2jvXlpGkBhnuktQgw12SGmS4ayF4c5N0fjyhOiErKyusr6/3XUazPLkqnR879wlxPRlJ88Rwl6QGGe6S1CDDXQvHk6vSeJ5Q1cLx5Ko0np27JDXIcJekBhnuktQgw10LzZOr0mieUN0H70rtnydXpdHs3PfBu1IlzSvDXZIaZLhLUoMMd0lqkOGuZnjljPQsr5ZRM7xyRnqWnbskNahTuCe5OsnDSc4kuWXE9kuT3JPkI0nuT3Lt5EuVJHU1NtyTHABuB64BjgE3JDm2Y7efB95bVVcAx4F3TLpQ6Xw4f9ey6zJzvxI4U1WPACS5C7geeGhonwK+ePD4JcBjkyxSOl/O37XsuoxlLgEeHXq+MXht2NuANybZAE4Dbx71RklOJFlLsra5ubmHcvu3srLyTEcoSfOqS7iPSrGd99vfANxZVYeBa4F3J3nee1fVyapararVQ4cOnX+1c8AlByQtgi7hvgEcGXp+mOePXW4C3gtQVX8DvBA4OIkCJUnnr0u43wdcnuSyJBeydcL01I59/hl4HUCSr2Mr3Bdz7qLmeHJVy2jsCdWqejrJzcDdwAHgjqp6MMltwFpVnQJ+CnhXkp9ka2RzYzm30Jzw5KqWUac7VKvqNFsnSodfu3Xo8UPAaydbmiRpr7xDVZIaZLhLUoMMdy2V4ZOrnmBVy1wVUktl+OQqeIJV7bJzl6QGGe6S1CDDXUvNG5zUKmfuWmre4KRW2bl34EqQkhaN4d6BK0FKWjSGuyQ1yHCXBjy5qpZ4QlUa8OSqWmLnLkkNMtwlqUGGuzSC83ctOmfu0gjO37Xo7NwlqUGGuyQ1yHCXpAYZ7tIYnlzVIvKEqjSGJ1e1iDp17kmuTvJwkjNJbtlln+9O8lCSB5P89mTLnD1XgpS0yMZ27kkOALcD/xXYAO5LcqqqHhra53LgvwGvrarPJvnyaRU8K9srQYLdmqTF06VzvxI4U1WPVNUXgLuA63fs84PA7VX1WYCqenyyZUrzwfm7FkWXmfslwKNDzzeAV+/Y52sAkvwVcAB4W1X9yUQqlOaI83ctii7hPupv8M5vrXgBcDlwFXAY+MskL6uqzz3njZITwAmASy+99LyLlSR102UsswEcGXp+GHhsxD5/WFX/UVUfBx5mK+yfo6pOVtVqVa0eOnRorzVLksboEu73AZcnuSzJhcBx4NSOff4A+DaAJAfZGtM8MslCpXnj/F3zbOxYpqqeTnIzcDdb8/Q7qurBJLcBa1V1arDtO5I8BPw/4Ker6tPTLFzqm/N3zbNONzFV1Wng9I7Xbh16XMBbBj+SpJ65/IAkNchwlybA+bvmjWvLSBPg/F3zxs5dkhpkuEtSgwz3Ia4EqUlw/q554Mx9iCtBahKcv2se2LlLUoMMd2mKHNGoL45lpClyRKO+2LlLUoMMd0lqkOEuzYjzd82SM3dpRpy/a5bs3CWpQYa7JDXIcJd64Pxd0+bMXeqB83dNm527JDVo6cPdlSDVN0c0moalH8u4EqT65ohG07D0nbsktchwl+aIIxpNSqdwT3J1koeTnElyyzn2e0OSSrI6uRKl5bE9oqkq1tfX+y5HC2xsuCc5ANwOXAMcA25IcmzEfi8Gfgz40KSLlCSdny6d+5XAmap6pKq+ANwFXD9iv18A3g78+wTrkyTtQZdwvwR4dOj5xuC1ZyS5AjhSVX98rjdKciLJWpK1zc3N8y5WWibO37UfXS6FHHVtVj2zMfki4JeBG8e9UVWdBE4CrK6u1pjdpaXmJZLajy6d+wZwZOj5YeCxoecvBl4G/EWSs8BrgFOeVJWk/nQJ9/uAy5NcluRC4DhwantjVT1RVQeraqWqVoB7geuqam0qFUtLyBGNztfYsUxVPZ3kZuBu4ABwR1U9mOQ2YK2qTp37HSTtlyMana9Oyw9U1Wng9I7Xbt1l36v2X5ak3Wx38QBHjx7l7Nmz/RakubT0a8tIi8YuXl0s5fIDrgQpqXVLGe7bK0Fudz/SovJEq3bjWEZaYI5otJul7NwlqXWGu9QIRzQa5lhGaoQjGg2zc5ekBhnuUoMc0cixjNQgRzSyc5caZxe/nOzcpcbZxS8nO3dJapDhLi0RRzTLw7GMtEQc0SyPpencXQlSei67+LYtTee+vRIk2LFIYBffuqXp3CXtbriLt5Nvw9J07pJ2N9zFg518C+zcJT2P8/jFZ+cu6Xmcxy8+O3dJ52QXv5g6hXuSq5M8nORMkltGbH9LkoeS3J/kz5McnXypkvqw3cVXFevr632Xo47GhnuSA8DtwDXAMeCGJMd27PYRYLWqXg68D3j7pAuV1D+7+MXRZeZ+JXCmqh4BSHIXcD3w0PYOVXXP0P73Am+cZJGS5oOz+MXRZSxzCfDo0PONwWu7uQl4/6gNSU4kWUuytrm52b1KSXPHLn6+dQn3Uf95rhGvkeSNwCrwS6O2V9XJqlqtqtVDhw51r1LS3BmexX/yk5806OdMl7HMBnBk6Plh4LGdOyV5PfBzwLdW1VOTKU/SInBcM3+6dO73AZcnuSzJhcBx4NTwDkmuAP4PcF1VPT75MvfGxcKk2XNcMx/GhntVPQ3cDNwNfAx4b1U9mOS2JNcNdvsl4EXA7yb5aJJTu7zdTG0vFjZ8W7Wk6XJcMx863aFaVaeB0zteu3Xo8esnXJekBjiu6Y93qEqaCcc1s2W4S5oJxzWz5cJhkmbOcc302blL6pXjmukw3CX1ynHNdBjukuaGQT85hrukuWTQ709z4e5dqVJ7DPrz11y4e1eq1DaDvpvmwl3S8jDod2e4S2qCQf9chruk5hj0hrukxi1r0BvukpbGMgW94S5pKe0W9BdffHEToW+4S1p6w0H/5JNPjgz9RQt7w12SdjEc+ovW4S98uA/fkepdqZKmqUuHPy+hv/DhPnxHqnelSupDl9CfddD7ZR2SNCV9finJwnfukqTnM9wlqUGGuyQ1qFO4J7k6ycNJziS5ZcT2i5L8zmD7h5KsTLpQSVJ3Y8M9yQHgduAa4BhwQ5JjO3a7CfhsVX018MvAL0660GF+IYcknVuXzv1K4ExVPVJVXwDuAq7fsc/1wG8MHr8PeF2mmLx+IYekRXPRRRfN9LLILpdCXgI8OvR8A3j1bvtU1dNJngC+DPjU8E5JTgAnBk//NcnDeykaOJjkmffe+d+R4ecNPj4IfOp8f+0c1L3fx7t+5vt5PCfHNpXPvIE/l6l85n3+uWxbX18fN3k4yI78HHL0XL9wW5dwH1XBzpa5yz5U1UngZIff89wFJWtVtbrf91lEy3rsy3rcsLzHvqzHDZM59i5jmQ3gyNDzw8Bju+2T5AXAS4DP7KcwSdLedQn3+4DLk1yW5ELgOHBqxz6ngO8fPH4D8IFyIC5JvRk7lhnM0G8G7gYOAHdU1YNJbgPWquoU8OvAu5OcYatjPz7NopnAaGeBLeuxL+txw/Ie+7IeN0xifG2DLUnt8Q5VSWqQ4S5JDZrrcM8SL3vQ4dhvTLKZ5KODnzf1UeekJbkjyeNJHthle5L8yuDP5f4kr5x1jdPQ4bivSvLE0Od966xrnIYkR5Lck+RjSR5M8uMj9mn1M+9y7Hv/3Ie/6GKeftg6eftPwFcCFwJ/Bxzbsc+PAO8cPD4O/E7fdc/w2G8E/nfftU7h2P8L8ErggV22Xwu8n617K14DfKjvmmd03FcBf9x3nVM47pcCrxw8fjHwDyP+rrf6mXc59j1/7vPcuc/dsgcz1OXYm1RVH+Tc90hcD/xmbbkX+JIkL51NddPT4bibVFWfqKq/HTz+v8DH2LrjfVirn3mXY9+zeQ73Ucse7Dzw5yx7AGwve7Douhw7wHcO/pn6viRHRmxvUdc/mxZ9Y5K/S/L+JF/fdzGTNhirXgF8aMem5j/zcxw77PFzn+dwn9iyBwuoy3H9EbBSVS8H/oxn/wXTulY/83H+FjhaVa8A/hfwBz3XM1FJXgT8HvATVfX5nZtH/JJmPvMxx77nz32ew32Zlz0Ye+xV9emqemrw9F3Aq2ZUW9+6/L1oTlV9vqr+dfD4NHBBkoM9lzURSS5gK9x+q6p+f8QuzX7m4459P5/7PIf7Mi97MPbYd8wcr2NrXrcMTgHfN7iC4jXAE1X1ib6LmrYkX7F9PinJlWz9f/fT/Va1f4Nj+nXgY1X1P3bZrcnPvMux7+dz77IqZC9qPpc9mImOx/5jSa4Dnmbr2G/sreAJSvIetq4QOJhkA3grcAFAVb0TOM3W1RNngCeBH+in0snqcNxvAH44ydPAvwHHG2lkXgt8L/D3ST46eO1ngUuh7c+cbse+58/d5QckqUHzPJaRJO2R4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLA0l+aGjd7I8nuafvmqS98iYmaYfBeh8fAN5eVX/Udz3SXti5S8/3P9lap8hg18Ka27VlpD4kuRE4CtzccynSvjiWkQaSvIqtdfG/pao+23c90n44lpGedTPwpcA9g5Oqv9Z3QdJe2blLUoPs3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJatB/Aq3INiHIhMzLAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "%matplotlib inline\n",
+    "import matplotlib.pyplot as plt\n",
+    "plt.bar(p['loc'].values,p['count'].values,dz,color='white',edgecolor='black')\n",
+    "plt.xlabel(\"z\");"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## With a User Defined Function (udf)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Although absolutely unnecessary here, we illustrate how one could have plugged an external UDF (to compute the bin number. Performances are extrememly dowgraded so you should NOT do it this way (see next part)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "111.3s\n"
+     ]
+    }
+   ],
+   "source": [
+    "#define your UDF (here just binning)\n",
+    "binNumber_udf=F.udf(lambda z: int((z-zmin)/dz))\n",
+    "timer.start()\n",
+    "#create the bin number column\n",
+    "p_udf=gal.select(gal.z,binNumber_udf(gal.z).alias('bin')).groupBy(\"bin\").count().orderBy(F.asc(\"bin\")).toPandas()\n",
+    "timer.stop()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "That's not great.Let see how to be more efficient"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## With pandas_udf"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The recommended way to bind Spark to external python functions is through pandas_udf which achieve some level of vectorization the previous method miss."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "46.1s\n"
+     ]
+    }
+   ],
+   "source": [
+    "import pandas as pd\n",
+    "@pandas_udf(\"float\", PandasUDFType.SCALAR)\n",
+    "def binFloat(z):\n",
+    "    return pd.Series((z-zmin)/dz)\n",
+    "timer.start()\n",
+    "p_udf=gal.select(gal.z,binFloat(\"z\").astype('int').alias('bin')).groupBy(\"bin\").count().orderBy(F.asc(\"bin\")).toPandas()\n",
+    "timer.stop()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Tomography\n",
+    "This is a popular method in cosmology when one builds \"shells\" in some redshift range of over-densities and compute their 2D power-spectrum (using some spherical pixelization as HealPix). As the redshift decreases, matter clusters more and more and it is expected that these spectra increase in amplitude. Some cross-correlations between nearby bins also appear because of the smearing between the true redhifts (\"z\") and the observed ones (\"zrec\")."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# define a function that filters galaxies within a redshift range and prject them onto a healpix map\n",
+    "# that is returned\n",
+    "\n",
+    "import numpy as np\n",
+    "import healpy as hp\n",
+    "\n",
+    "nside=512\n",
+    "\n",
+    "def redshift_bin(gal,z1,z2):  \n",
+    "    timer.start()\n",
+    "    #filtering\n",
+    "    shell=gal.filter(gal['zrec'].between(z1,z2))\n",
+    "    print(\"shell=[{},{}] N={}M\".format(z1,z2,shell.count()/1e6))\n",
+    "    \n",
+    "    #projection into healpixels: build a pandas_udf that transforms (ra,dec) into (theta,phi) \n",
+    "    # and calls healpy ang2pix function          \n",
+    "    @pandas_udf('int', PandasUDFType.SCALAR)\n",
+    "    def Ang2Pix(ra,dec):\n",
+    "        return pd.Series(hp.ang2pix(nside,np.radians(90-dec),np.radians(ra)))\n",
+    "    map=shell.select(Ang2Pix(\"RA\",\"Dec\").alias(\"ipix\")).groupBy(\"ipix\").count().toPandas()\n",
+    " \n",
+    "    #back to python world\n",
+    "    myMap = np.zeros(hp.nside2npix(nside))\n",
+    "    myMap[map['ipix'].values]=map['count'].values\n",
+    "    timer.stop()\n",
+    "    return myMap"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Do a loop over redshift bins  (choosen here so that there are some overlap between adajacent ones for the true redshift) and store the Healpix maps with the number of galaxies in each pixel in a dictionary for later use."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 54,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "shell=[0.0,0.13] N=68.557105M\n",
+      "27.8s\n",
+      "shell=[0.13,0.27] N=322.103946M\n",
+      "24.7s\n",
+      "shell=[0.27,0.43] N=669.646368M\n",
+      "30.3s\n",
+      "shell=[0.43,0.63] N=1049.274954M\n",
+      "35.3s\n",
+      "shell=[0.63,0.82] N=980.216568M\n",
+      "37.5s\n",
+      "shell=[0.82,1.05] N=984.086242M\n",
+      "33.9s\n",
+      "shell=[1.05,1.32] N=812.257315M\n",
+      "30.7s\n",
+      "shell=[1.32,1.61] N=529.298439M\n",
+      "30.5s\n",
+      "shell=[1.61,1.95] N=325.212359M\n",
+      "27.8s\n",
+      "shell=[1.95,2.32] N=156.932408M\n",
+      "26.4s\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "<Figure size 1080x720 with 0 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "#redshift bins\n",
+    "zshell=[0.0,0.13,0.27,0.43,0.63,0.82,1.05,1.32,1.61,1.95,2.32]\n",
+    "tmap=dict()\n",
+    "\n",
+    "#loop over bins and store healpix maps in a dictionary indexed by z1\n",
+    "plt.figure(figsize=(15,10))\n",
+    "for (z1,z2) in zip(zshell[:-1],np.roll(zshell,-1)): \n",
+    "    tmap[z1]=redshift_bin(gal,z1,z2)    "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now do a standard analysis as the density contrast power spectra with (standard) healpix functions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 59,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "bin=[0.0,0.13] mean #gals=21.7937167485555\n",
+      "bin=[0.13,0.27] mean #gals=102.39408683776855\n",
+      "bin=[0.27,0.43] mean #gals=212.87484741210938\n",
+      "bin=[0.43,0.63] mean #gals=333.5555248260498\n",
+      "bin=[0.63,0.82] mean #gals=311.60245513916016\n",
+      "bin=[0.82,1.05] mean #gals=312.83259137471515\n",
+      "bin=[1.05,1.32] mean #gals=258.209646542867\n",
+      "bin=[1.32,1.61] mean #gals=168.25944232940674\n",
+      "bin=[1.61,1.95] mean #gals=103.38222471872966\n",
+      "bin=[1.95,2.32] mean #gals=49.88746897379557\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "Text(0,0.5,'power spectrum')"
+      ]
+     },
+     "execution_count": 59,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEKCAYAAAAiizNaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXl8lOW1x78nkxWSsCZACPuWEIRIWNyKCkVxpRT0yqUtVVzqpa2t9V6r7a3XWlutF+tu61brilr0ShFxi4oiGIOAhNUIEcIalhBIyDZ57h/vzGQymeXNMpmQnO/nMx8mz/tsk2hOznnOc35ijEFRFEVR2oqoSG9AURRF6Vyo4VEURVHaFDU8iqIoSpuihkdRFEVpU9TwKIqiKG2KGh5FURSlTVHDoyiKorQpangURVGUNkUNj6IoitKmREd6A5Gkd+/eZvDgwZHehqIoyinF2rVrDxljUpo7vlMbnsGDB5Ofnx/pbSiKopxSiMi3LRkf1lCbiMwQkW0iUigiv/bzPE5EXnE9/1xEBns9u83Vvk1ELvRqf0ZEDopIgc9c94nIVhH5SkTeEJHu4fxsiqIoSvMIm+EREQfwKHARMBqYKyKjfbotAI4aY4YDfwHudY0dDVwFZAEzgMdc8wE862rz5T1gjDFmLLAduK1VP5CiKIrSKoTT45kEFBpjdhhjqoHFwEyfPjOBf7je/xOYJiLial9sjKkyxuwECl3zYYxZCRzxXcwY864xptb15RogvbU/kKIoitJywnnG0x/Y7fV1MTA5UB9jTK2IHAN6udrX+Izt34S1rwFe8fdARK4HrgcYOHBgE6ZUFKWp1NTUUFxcTGVlZaS3ojSD+Ph40tPTiYmJadV5w2l4xE+br/hPoD52xvpfVOQ3QC3wor/nxpgngCcAJkyYoGJEihJGiouLSUpKYvDgwVjBDOVUwRjD4cOHKS4uZsiQIa06dzhDbcXAAK+v04G9gfqISDTQDSuMZmdsI0RkPnApMM+owp2iRJzKykp69eqlRucURETo1atXWLzVcBqeL4ARIjJERGKxkgWW+vRZCsx3vZ8D5LoMxlLgKlfW2xBgBJAXbDERmQHcClxujKloxc+hKEoLUKNz6hKun13YDI/roP+nwDvAFuBVY8wmEfm9iFzu6vY00EtECoGbgV+7xm4CXgU2AyuAhcYYJ4CIvAysBkaJSLGILHDN9QiQBLwnIutF5K+h9lhZ42ylT6soiqLYRTpzRKr34ExzqGhLpLehKB2WLVu2kJmZGeltKC3A389QRNYaYyY0d85OXavN2YmNrqIoSqTo3IanTg2PoiihKSoqIiEhgezsbE/bihUrGDVqFMOHD+eee+4JONZuv6b0veaaa0hNTWXMmDGetsrKSiZNmsS4cePIysrijjvu8Dw7efIk2dnZxMbGcujQITsfObwYYzrtq0vaCKMoSvjYvHlzpLdgiyeffNJMmDDB5OTkeF6XXXaZ5/nOnTtNVlaW5+va2lozdOhQ880335iqqiozduxYs2nTpkbz2u3X1L4ff/yxWbt2bYM91dXVmePHjxtjjKmurjaTJk0yq1evbjBu0KBBpqSkxP43xvj/GQL5pgW/ezt1kVANtSlK52Dq1KkcOWIVPNm6dSvPP/88V1xxBQCLFi1i48aN5ObmkpSUZGu+vLw8hg8fztChQwG46qqrePPNNxk9enSz+jW175QpUygqKmrQJiIkJiYC1sXdmpqadptR2KkNjzFQVeskLtoRurOiKC3izn9tYvPesladc3RaMndclhWyX25uLgCPP/44H374Id///vcBcDqdvPTSS+Tl5eFw2P89sGfPHgYMqL9qmJ6ezueff97sfk3tGwin00lOTg6FhYUsXLiQyZN9i8W0Dzq14QE4XllLXKIaHkXp6Dz33HO8/fbbLFmyxGNkSkpKKCwsJCcnp1H/tLQ0li9f7ncu4yda4s+7sNuvqX0D4XA4WL9+PaWlpcyaNYuCgoIG50DthU5veMpO1tA7MS7S21CUDo8dzyRcvPbaa7z44ou8+eabDeqOpaSkMGTIENauXdskjyc9PZ3du+tLURYXF5OWltbsfk3tG4ru3btz3nnnsWLFinZpeDp1VhtYHo+iKB2XZcuW8dhjj/H6668THx/f4JnD4eCHP/wh1157LcePH7c958SJE/n666/ZuXMn1dXVLF68mMsvv7xJ/aZNm8aePXuaPGcgSkpKKC0tBawstvfff5+MjAzb49uSTm94yiprIr0FRVHCyPz58ykuLubss88mOzubp59+usHzX/3qV5x11llMnTqVCRMmeF7BfulHR0fzyCOPcOGFF5KZmcmVV15JVla9R3fxxRezd+/egP3q6uooLCykZ8+eTZ4TYO7cuZx55pls27aN9PR0nn76afbt28f555/P2LFjmThxItOnT+fSSy9trW9jq9KpKxfE9RthXn9nJZeM7RfprShKh6SjVC4oKiri0ksvpaCgIHRnGxQUFPDMM89w//33t8p8dhk8eDD5+fn07t3b9hitXBAG1ONRFCUUDoeDY8eONbhA2hLGjBnTpkbHfYG0pqaGqKjI/9rv9MkFx30Mz9cHjjM8NbHd5r8ritL2DBgwoMHB/6lGQkIC69evj/Q2PETe9EWYspP1yQU7D5Uz/S8rWfl1OygpoSiK0kFRw+Pl8RQftWR8SiuqI7UdRVGUDo8anpP1hufwCcvg1Dg7b8KFoihKuOn0hsf7Hs+hE1UA1DjrIrUdRVGUDk+nNzzeobbD5W6PRw2PoihKuFDD45VccOi42+PRUJuiKEq4UMOjHo+iKEqb0qkNj+CbXODyeGrV8CiKUk9LFEibqhYaiFDr7d69m/PPP5/MzEyysrJ48MEHPc+2bdtGdna255WcnMwDDzwQOWXSlqjIneqvrmkjzaBbl5maWqcxxpiz/vSBGXTrMrPo3W1BFfkURbFHZ1cgNab5aqHe2Flv7969Zu3atcYYY8rKysyIESMCqqL26dPHFBUVedqCKZOqAmkr44iyqhOcqKqlW0IMh8s1q01ROiKRUiCF1lELtbNev3796NfPqjuZlJREZmYme/bsabSnDz74gGHDhjFo0CBbnzUcdGrDE+2wftBlJ2uJcURRWWMZHA21KUrHIlIKpMFoilpoU9crKipi3bp1fudcvHgxc+fObdJeW5tObXjcHo9vodDaOs1qU5RW5+1fw/6NrTtn39PgosDnK95EQoE0GE1RC23KeidOnGD27Nk88MADJCcnN3hWXV3N0qVL+dOf/tSkvbY2YTU8IjIDeBBwAE8ZY+7xeR4HPAfkAIeBfzPGFLme3QYsAJzAz40x77janwEuBQ4aY8Z4zdUTeAUYDBQBVxpjjgbbn7fhqfLycqo11KYoHYpIKZDawY5aqN31ampqmD17NvPmzfN4dd68/fbbjB8/nj59+jRrr61F2AyPiDiAR4HpQDHwhYgsNcZs9uq2ADhqjBkuIlcB9wL/JiKjgauALCANeF9ERhpjnMCzwCNYBsubXwMfGGPuEZFfu76+Ndgeo6PqQ21RUn+fR0NtihIGbHomrY1bgXTZsmVBFUgfeugh22c83mqh/fv3Z/Hixbz00ku291RSUkJMTAzdu3f3qIXeemv9r6tp06bx3HPP0b9/f9vrGWNYsGABmZmZ3HzzzX7XffnllyMeZoPwplNPAgqNMTuMMdXAYmCmT5+ZwD9c7/8JTBPLf5wJLDbGVBljdgKFrvkwxqwEjvhZz3uufwDfC7VBb4/HfYcn1hGloTZF6UBEUoEUmq4W2lR1Uvdaq1at4vnnnyc3N9eTNu0dKqyoqOC9997z6wm1NeEMtfUHvAUsigHfky5PH2NMrYgcA3q52tf4jO0fYr0+xph9rrn2iUhqqA1GRwm1WHd5KmucAKQmx2moTVE6EIcPHw7Z57rrruO6665r0rwXX3wxF198sd9n3r/wX375Zb991q1b57d98+bNzJ49m4SEBFvruddKS0vzexbkpkuXLra+F21BOD0efydfvt+VQH3sjG0WInK9iOSLSP6J48cBq1DooRPVJMVFkxgXTa0aHkVRvGhtBdJgtKU6aaSUScPp8RQDA7y+Tgf2BuhTLCLRQDesMJqdsb4cEJF+Lm+nH3DQXydjzBPAEwAZp2WbmLhoT6itd1IcMY4ordWmKEoDTnUF0kBESpk0nCbuC2CEiAwRkVisZIGlPn2WAvNd7+cAua5bsUuBq0QkTkSGACOAvBDrec81H3jTziaTE2IoO1nL4RNV9OoaS7RD9AKpoihKGAmb4THG1AI/Bd4BtgCvGmM2icjvRcR9avc00EtECoGbsTLRMMZsAl4FNgMrgIWujDZE5GVgNTBKRIpFZIFrrnuA6SLyNVYmXcgUGgGS4l0ez4lqeiXGujweNTyKoijhIqz3eIwxy4HlPm2/83pfCVwRYOzdwN1+2v3mAhpjDgPTmrI/EUiOj+F4ZQ2HTlSRM7gHJ6pqPRUMFEVRlNanU1enBiE5IZrSihqOVFTTu6t6PIqiKOGmUxset8ez60gFxkCvRE0uUBRFCTed2/BgJRdUVFt3eHonxhGjyQWKoihhpVMbHhCS4uuPuTS5QFEUJfx0asPjDrW56e0yPLUaalMURQkbndvwAMkJXh5PVyvUpiVzFEXxprnS14HkqANJUQfDrtS20+nk9NNP99R+g8BS25GSvu7UejxQ7/E4ooRuCTEaalOUTshTTz3F3/72twa1ztLS0li6tP7O+7Bhwzy3/J1OJwsXLuS9994jPT2diRMncvnllzdS+4yOjmbRokWMHz+e48ePk5OTw/Tp0xk9enSDufr378+sWbMC7s/uegAPPvggmZmZlJWVedri4uLIzc0lMTGRmpoazjnnHC666CLOOOMM1q9fz+DBg5v1fWsundvjEUhyGZ6eXWOJihINtSlKB2Tq1Kke7yI+Pp7XXnvN82zRokV8+umn5Obmkp+f73l5Gx1fvKWoY2NjPVLUvvTr14/x48cDDeWovbEjRW13veLiYt566y2uvfbaBu1NldoON53a4xHXPR6AXl1jAUsOW0NtitL63Jt3L1uPbG3VOTN6ZnDrpKCyW0D7kL4OJEdtR4ra7nq/+MUv+POf/8xxVwFkb5oitR1uOrXhwSu5ICUpDrD0eDTUpigdj0hKXweSo7YrRW1nvWXLlpGamkpOTg4fffRRo/5NkdoON53a8Ljv8UC9xxPjiMIYcNYZj1Ccoigtx45nEi4iKX0dTI7arhS1nfVWrVrF0qVLWb58OZWVlZSVlfGDH/yAF154oUE/O1Lb4UbPeFz3eHolWh5PtMMyNur1KErHwC19/frrrweVvvYXngqEtxR1dXU1ixcv9qtYGkqOOpAU9bRp0xqcBdlZ709/+hPFxcUUFRWxePFipk6d6jE6JSUllJaWAniktjMyMmx/3tamUxsesJIJbp4+ku9lWwKnsQ7rW6LnPIrSMYik9HUwOepAUtTNlb4ORjCp7UggwaRSOzo5ORPM2rX5Ddr+8VkRdyzdxJf/PZ2ervCboijNY8uWLWRmZkZ6Gy2mqKiISy+9lIKCgrCvVVBQwDPPPNNmKqQAgwcPJj8/n969ezd65u9nKCJrjTETmrtep/Z4/J0FaqhNURRfVPq6denUyQX+iHGH2mrV8CiKYqHS161Lp/Z4/OE+46mt67whSEVRlHCihscHDbUpiqKEFzU8PmioTVEUJbyo4fFBQ22KoijhRQ2PDxpqUxRFCS9qeHxwh9rU8CiKooQHNTw+1BseDbUpiqKEAzU8PsS4Q22aXKAoihIW1PD4EONJLlDDoyiKRXOlr934ylEHkqIOhp31SktLmTNnDhkZGWRmZrJ69eqg60VK+jqshkdEZojINhEpFJFf+3keJyKvuJ5/LiKDvZ7d5mrfJiIXhppTRKaJyJcisl5EPhWR4c3ZsyedWkNtitJpeOqpp5g4cWLQIqH+pK/ffvttNm/ezMsvv8zmzZsDzu+Wo3bjlqLesGED69evZ8WKFaxZsybgeLvr3XTTTcyYMYOtW7eyYcMGz5qB1nNXLggk6RAubBseEUkWkZ7ul43+DuBR4CJgNDBXRHwFwhcAR40xw4G/APe6xo4GrgKygBnAYyLiCDHn48A8Y0w28BLwW7ufzRsNtSlKxyNS0tfgX466qVLUdtYrKytj5cqVLFiwAIDY2Fi6d+/erPXCTchabSJyA/B74CTgdgMMMDTE0ElAoTFmh2uexcBMwNtMzwT+x/X+n8AjYn03ZgKLjTFVwE4RKXTNR5A5DeCW9usGBK8THoDmhNqqa+u48YW1/OqCUYxOSw49QFE6Ifv/+EeqtrSu9HVcZgZ9b789ZL9ISl8HkqNuihS1nfV27NhBSkoKV199NRs2bCAnJ4cHH3yQrl27Nnm9cGPH47kFyDLGDDbGDHG9QhkdgP6Ad1W9Yleb3z7GmFrgGNAryNhgc14LLBeRYuCHQPCgawCaE2o7UFbJB1sP8mlhSXOWVBSlDXBLX7/44ot+pa/dHpH7dfHFFwecy670tbcctS9uKeri4mLy8vKCSi7YWa+2tpYvv/ySG2+8kXXr1tG1a9cGZ0FNWS/c2KlO/Q1Q0Yy5/flxvt+9QH0CtfszlO45fwlcbIz5XET+E7gfyxg1XFDkeuB6gIEDBzaarDmhtipX30Mnqm2PUZTOhh3PJFxESvrajhy1HSlqO+ulp6eTnp7u8WTmzJnjNwnhVJG+vg34TET+JiIPuV82xhUDA7y+Tqdx+MvTR0SisUJkR4KM9dsuIinAOGOM2/d8BTjL36aMMU8YYyYYYyakpKQ0et6cUFtVrROAkuNVtscoitI2RFL6OpAcdSgp6uZIX/ft25cBAwawbds2AD744ANGj7aOwNub9LUdj+dvQC6wEWjKifsXwAgRGQLswUoW+HefPkuB+cBqYA6Qa4wxIrIUeElE7gfSgBFAHpYn5G/Oo0A3ERlpjNkOTAe2NGGvHppzgdTt8ajhUZT2x/z58+nZsydnn302AD/72c88B/BgSV8/+eSTTJ06tUFIKy0tLWCCgbcUtdPp5Jprrmkkff3UU08FzBbbt28f8+fPx+l0UldXx5VXXulJtQ4lfe27nvdaDz/8MPPmzaO6upqhQ4fy97//PeR6EcEYE/QFfBaqT5CxFwPbscJ1v3G1/R643PU+HngNKMQyLEO9xv7GNW4bcFGwOV3ts7CM4wbgI++5Ar1ycnKML3V1dWbQrcvM/e9ua/QsEKu/OWQG3brMXPiXj22PUZTOwObNmyO9hVZh586dJisrq03W2rhxo/nlL3/ZJmu5GTRokCkpKfH7zN/PEMg3zbQLxhhbHs+HrnORfwGeP+mNMUdsGLXlwHKftt95va8Erggw9m7gbjtzutrfAN4ItadQiAjRUdLEUJt6PIrSkfGWvg63YmdbS1+feeaZ7VL62h0eu82rzU469SlLjCOqaaG2GuuM50hFNTXOOk+4TlGUjoFKX7cuIQ2PMWZIW2ykPRHtkCYJwVW7KlkbA0fKq+mTHB9ihKIoSufFzgXSH/lrN8Y81/rbaR/EOqKaFmqrqe9bcrxKDY+iKEoQ7ITaJnq9jwemAV8CHdbwxDiiqKltelYbQMkJPedRFEUJhp1Q28+8vxaRbsDzYdtROyDaIU0SgnPf4wFNMFAURQlFc07BK7Du1XRYYh1R1NTZ93i8z4MOqcejKIoSFDtnPP+ivixNFFZV6FfDualIY4Xamp5OnRDjUI9HURQlBHbOeP7X630t8K0xpjhM+2kXNCfU5ogS+naLV8OjKIoSgqCGx6V/89/GmO+20X7aBTE+obbSimr2llYGlDyorq0j1hFFSmKcGh5FUZQQBD3jMcY4gQpXQkGnIdYn1PbEyh1c/sinfH3AfxHBqto64mKiSEmK0zMeRemAtET62p8cdbikr//yl7+QlZXFmDFjmDt3LpWVlQDs3r2b888/n8zMTLKysnjwwQeB9i19XQlsFJGnm1id+pTFN9R2tKKa2jrD797c5FcXo6qmjrjoKHonxqrHoyinIOGUvvYnRx0O6es9e/bw0EMPkZ+fT0FBAU6nk8WLFwNWkdFFixaxZcsW1qxZw6OPPsrmzZsjJn1t54znLdfLG/spX6cgMY4oyqvrU6QrXO9X7zjM0g17mZndUM+uqtZJXLSDlKQ4yiprqaxxEh9jX9tDUZTwMnXqVI4cscpLbt26leeff54rrrDKRC5atIiNGzeSm5tLUlKSrfm8pagBjxS1W4bAjVuO+tlnnwUsOerY2FiAZktfB1uvtraWkydPEhMTQ0VFhceg9OvXj379+gGQlJREZmYme/bsaTS+rbBjeLobYx70bhCRm8K0n3ZBjEMahNoqqp2M7JNIXLSDu9/awtSMVJLi68Wkqp11xEZboTawUqrTe3Rp830rSnvmk1e3c2j3iVads/eARL5z5ciQ/SIlfR1Mjrq1pa/79+/PLbfcwsCBA0lISOCCCy7gggsuaDRXUVER69ata/fS1/P9tP24lffRrojxKZlTUV1LYlw0d31vDCUnqngkt7BB//pQm9vwqBKporQ3IiF9HUyOurWlr48ePcqbb77Jzp072bt3L+Xl5Q2UTgFOnDjB7NmzeeCBB0hO9p8s1RYE9HhEZC5WZeohLmE2N0nA4XBvLJL4VqeuqHaSGBdN9oDunDGkF2t2NlSEqKq1DI/b49FzHkVpjB3PJFxESvrajhx1a0lfv//++wwZMgS3svL3v/99PvvsM37wgx8AVkhv9uzZzJs3z+PxRYpgHs9nwCJgq+tf9+tXwIzwby1y+FanrqhykuA6s0mKj/bIILiprm0YalPDoyjth0hKXweSow6H9PXAgQNZs2YNFRUVGGP44IMPyMzMBCyPacGCBWRmZnLzzTfb/pzhIqDHY4z5FvhWROYBe12ibYhIApAOFLXJDiOAb3XqippausRahich1sFJH8NTVeuke5dYenWtP+NRFKV9EGnpa39y1Lt372516evJkyczZ84cxo8fT3R0NKeffjrXX389AKtWreL555/ntNNO86SE//GPfwwaTgwndpILXgXO8vraiSVXPdF/91Mf31DbyWonXeKsb1VCjIOT1b6Gxwq1xUZH0b1LjHo8itKOOHw49MnAddddx3XXXdekeS+++OKAv7iXL68XSc7OziY/P7/B8x49erBu3Tq/Yzdv3szs2bNJSEiwtZ73WnfeeSd33nlnoz7nnHOO33OiSGEnuSDaGOM5LXe9jw3fliJPtE9WW3mVky6uUFt8jIPKRh5PHXGu51q9QFE6Ht7S1+GmraWvs7Oz26X0dYmIXG6MWQogIjOBtrviGgGs6tSW4amrM5yscTYItVXWNKzj5i6ZA9A7MU41eRSlg6HS162LHcPzE+BFEXkU6+JoMeBXlbSj4B1qq3Rp7XiH2qqdddQ664h2GZuqWidxMdb7lKQ4NhSXRmDXiqIopwZ2hOC+Ac4QkURAjDH2Uz9OUaIdgrPOUFdnKK9yGZ5Yd6jNMjCVtXUkug2P6x4PWIZHQ22KoiiBCRnUE5E+IvI08Jox5riIjBaRBaHGncrEuAxKTV2dJ5GgS2y9xwM0SDCoclUuACvUVlHtpLyqlt1HKti091hbbl1RFKXdY+c06VngHcB9W2k78Itwbag94D6vqXEayqtrAW+Px/rXnWBgjKG6to64aFdygesuzyUPfcJ3/vwhsx77rIE0tqIoSmfHjuHpbYx5FagDMMbUYqVUd1iiHVYpipraOk+B0ASv5AKoNzxu9VF3qC0rLZkusQ76JMczNSOV6to6jp2sadP9K4qitGfsGJ5yEemFqyK1iJwB2IoficgMEdkmIoUi8ms/z+NE5BXX889FZLDXs9tc7dtE5MJQc4rF3SKyXUS2iMjP7ezRH/5CbV19Q20BDE9mv2Q23Xkhr9xwJjOzLSexTA2PoiiKBztZbTcDS4FhIrIKSAHmhBrkUi99FJiOlQn3hYgsNcZ4i0gsAI4aY4aLyFXAvcC/icho4CogCyvE976IuAs9BZrzx8AAIMMYUyciqTY+m1+Chdp8z3iqfQyP67MD0C3BqgmlHo+iKEo9drLavhSRc4FRgADbjDF2fpNOAgqNMTsARGQxMBPwNjwzgf9xvf8n8IhYv7VnAouNMVXAThEpdM1HkDlvBP7dGOMOCR60sUe/eIfa6pMLLIMT18jjsf51n/F4o4ZHURSlMXay2uKBnwN3AXcCC11toegPeN+4Kna1+e3jOjs6BvQKMjbYnMOwvKV8EXlbREbY2KNf3KG22rr6Mx7frLZGZzwxjb+VangUpWPQEulrf3LUgaSog3HNNdeQmpoasIJ1qD0NHjzYU6ttwoQJQPuWvn4OK+T1MPAIMBp43sY4f3J6vsWCAvVpajtAHFBpjJkAPAk843dTIte7jFN+SUmJ3427DU91raHCFWprnFxQ5+pj/esOz3njMTwVangUpT0TLunrQHLUgaSog/HjH/+YFStWBHxuZ08ffvgh69ev99SOa8/S16OMMeO8vv5QRDbYGFeMdebiJh3YG6BPsYhEA92AIyHGBmovBpa43r8B/N3fpowxTwBPAEyYMMFv1bwYd6jN6e3x+Jzx2PB4kj0eT62/ZRRFaSMiJX0N/uWomyNFPWXKFIqKilplT5HGjuFZJyJnGGPWAIjIZGCVjXFfACNEZAiwBytZ4N99+izFUjhdjZWwkGuMMS7huZdE5H6s5IIRQB6WxxNozv8DpmJ5Oudi3TdqFt6htvLqWmIdUZ423+QCtzaPvzOeGEcUXWMdGmpTFODDZ5/g4Lc7WnXO1EFDOf/H14fsFynpazty1K0lRR1qTyLCBRdcgIhwww03eCQTIoEdwzMZ+JGI7HJ9PRDYIiIbAWOMGetvkDGmVkR+inX51AE8Y4zZJCK/B/JdRUefBp53JQ8cwTIkuPq9ipU0UAssNMY4AfzN6VryHqyacr8ETgDXNuk74YV3qO1ktdMTXgOIj7We+Xo8sdH+o5bdEmLU8ChKO8Atfb1kyRK/0te+pKWlNZAc8Mau9LW3HHX37t254ooreOGFFzyqoK0pRR1qT6tWrSItLY2DBw8yffp0MjIymDJlSovWbC5awns+AAAgAElEQVR2DE+z1UaNMcuB5T5tv/N6XwlcEWDs3cDdduZ0tZcClzR3r974htq6ehmeWEcUIvXJBf7Sqb1JVsOjKAC2PJNwESnp62By1K0tRR1qT+73qampzJo1i7y8vIgZHlt6PMB+lyLpEKz05WPGmG9dbR2OhllttQ08HhFpIAZXf4HU/3+03RJi9AKpokSQSEpfB5KjDiVF7St93dI9lZeXez5feXk57777btDsuHBjx/AsAZwiMhwrNDYEeCmsu4owDbPanHSNa+gYJsQ4PHIJ7ns8GmpTlPbJ/PnzKS4u5uyzzyY7O5unn366wfNf/epXnHXWWUydOjVoVps33lLUmZmZXHnllY2kr/fu3dtAjvq0006jrq6O66+/3iNFnZubS3Z2NtnZ2Z6wnj/pa4C5c+dy5plnsm3bNtLT0z2fw71WsD0dOHCAc845h3HjxjFp0iQuueQSZsxodjCrxdgJtdW5zmu+DzxgjHlYRPxrtnYQGoTaqpyehAI38TEOTlY3TKcOFGpTw6MokSXS0tf+5KiDSVEHkr5++eWXQ64VaE9Dhw5lwwY7ychtgx2Pp0ZE5mKJvy1ztcUE6X/K0yDUVlPrSaV2Y6mQ+q/V5osaHkU59VHp69bFjsdzNZYK6d3GmJ2uVOYXwrutyBLjMiI1rlBbFz+htkYlc2ICn/GcrHFa8tgBjJOiKO0blb5uXezUatuMVTLH/fVOrNTlDktMlCvUVmeF2ro0CrVFNSoS6q9yAUC3LpZzWFZZQ+/EuHBtWVEU5ZRB/wT3g0cWodbKavMNtcU38HjqEKk/F/JF67UpiqI0RA2PHzyhNqfhZI3/UJv3GU9cdJTfy2PgXTZHDY+iKAqEMDwi4hCR+9pqM+2FaFeoraLaSY3TNAq1NUguqHEGDLOBejyKoii+BDU8rjI1ORLoz/kOijvU5jYWwZILqp11ARMLoN7w6CVSRVEUC1tFQoE3ReQ1oNzdaIx5PWy7ijCOKMERJfWGx98Zj6dIaF3AVGqA5Hj1eBRFUbyxY3h6AoexKj+7MUCHNTxghduCGR63Ho/7jCcQqsmjKIrSkJDJBcaYq/28rmmLzUWSWEeUJzzmVh91kxDjoNpZh7POUFVbR2yAOm1gldJJiFFpBEU5lfGnQGpHERQCq4L6UwQNRqj1KisrmTRpEuPGjSMrK4s77rjD8yyQ4mm7VSAVkZEi8oGIFLi+Hisivw3/1iJLTHRUQI8nwSWNUFnjpKrWGdTjAa1eoCjtnaYqkEJoRVAIrQrqqwgajFDrxcXFkZuby4YNG1i/fj0rVqxgzZo1AAEVT9uzAumTwH8CfwMwxnwlIi8BfwjnxiKNd6gtwdfweKmQhgq1gRoeRYk0ra1ACqEVQaF1VUFDrSciJCYmAlBTU0NNTY3nmkdzFE/DiZ17PF2MMXk+bR1eyznGUe/xdPUJtcV7qZBW2SiFo4ZHUSJLbm4u69ev54YbbuDyyy9vpED69NNPN8no2MWfKqhb7sCtCJqTk8MTTzzRKus5nU6ys7NJTU1l+vTpflVNW0vxtCXY8XgOicgwrIQCRGQOsC+su2oHxEZHeVKm/RUJBSvUVl1bF1CLx01yQgzFRyvCs1FFOUUo/dc3VO8tD92xCcSmdaX7ZcNs9W1NBVK7BFMFDYciqMPhYP369ZSWljJr1iwKCgoanAm1puJpS7BjeBYCTwAZIrIH2AnMC+uu2gHuS6TgJ6st2jvU5iQuJrTHs3mvejyKEilaW4HULsFUQcOpCNq9e3fOO+88VqxY4TE8ra142hLsFAndAXxXRLoCUcYY+zJ9pzAxXtUIGmW1eTyeOuseT5DKBaChNkUBbHsmrY1bgXTZsmVBFUgfeuihVg+3eauC9u/fn8WLF/PSSy9RXl5OXV0dSUlJHkXQ3/3ud55x06ZN47nnnqN///621yopKSEmJobu3btz8uRJ3n//fW699VaAkIqnbY2drLZvRORF4IfAgFD9Owruem0iVjVqb+K9kgusygWhDU95tZMaZ114NqsoSkDCoUAKgRVBIbQqaDBF0OYqkO7bt4/zzz+fsWPHMnHiRKZPn86ll14KEFTxNBLYCbWNBiYD3wH+V0QygA3GmFlh3VmEcUsjdIlxNCoAmuCdXFDjDHnG0y3B+jaXnayhV2Ictc46RKzqCIqihJdwKZAGUgSF0KqgwRRBm6tAmpaWxrp1/sWhgymeRgI7WW1OoMb1bx1wADgYzk21B9yhtoTYxrbZO7nAVjp1l4Zlc2b/dTW//b+NrbldRVHCiCqQti52PJ4yYCNwP/CkMSb0nw8dAHeozTexABrf47GTTg2W4dlTepINu0vZsq+M/7owgx5dY1t554qitDaqQNq62DFxc4GVwH8Ai0XkThGZFt5tRR5PqM2P4XGf+RyvtDwYOxdIwTI8n2wvASzl0iVfFrfafhVFUU4V7NRqe9MY85/ADcBy4MfAsjDvK+K4Q23+DY/VVlrhNjyhzni8DM/Xh+iTHMfpA7vzct6udhV3VRRFaQvsZLUtEZFvgAeBROBHQA87k4vIDBHZJiKFIvJrP8/jROQV1/PPRWSw17PbXO3bROTCJsz5sIicsLO/YLhDbV3jGkcjLcVRKHWd2YQKtblVSI+WV/Np4SGmjEjh3ycN5JuScvJ2HmnpVhVFUU4p7ITa7gFGGmMuNMbcZYz52BhTGWqQiDiAR4GLsDLj5oqIb2GgBcBRY8xw4C/Ava6xo4GrgCxgBvCYSw016JwiMgHobuMzhcQdakvwI/ImIg0qTtsNtX1aeJhjJ2v4zsgULh2bRlJ8NC/n7WqN7SqKopwy2DE864GFIvJP1+tnIhITchRMAgqNMTuMMdXAYmCmT5+ZwD9c7/8JTHOpnc4EFhtjqowxO4FC13wB53QZpfuA/7Kxt5AEC7WBZZDcGjuh7vHERTuIj4li5fYSROCc4b1JiHXw/dP7s7xgP0fLq1tjy4qiKKcEdgzP40AO8JjrNd7VFor+gHcaSLGrzW8fY0wtcAzoFWRssDl/Ciw1xgStIyci14tIvojkl5SUBOwX7XAlF/gJtYF1zlN60jIYoc54wPJ6qp11nNa/Gz1dmWxzJw+kuraOZV/tDTleURSlo2DH8Ew0xsw3xuS6XlcDE22M83c70vckPVCfJrWLSBpwBfBwqE0ZY54wxkwwxkxISUkJ2M/j8fgJtYGV2eYOtcWGKJkD9eG274zo7WnL6JtMSlIc63cfCzleURSlo2DrAqmrOjUAIjIU6zJpKIppWGInHfD9097TR0SigW7AkSBjA7WfDgwHCkWkCOgiIoU29hiQ2CD3eMC6RGo31Abehqehscvsl8yWfWUt2aqiKMophR3D85/AhyLykYh8DOQCv7Ix7gtghIgMEZFYrGSBpT59lgLzXe/nALnGyi9eClzlynobAowA8gLNaYx5yxjT1xgz2BgzGKhwJSw0G3d16kChtoQYB2WVliyR3VBb11gH4wc2TAjM7JdE4cETVNdqHTdFaa80V/o6kBx1ICnqYNhZr7S0lDlz5pCRkUFmZiarV68OOr7dSl8bYz7A+sX/c9drlDHmQxvjarHOXd4BtgCvGmM2icjvRcRdfe9poJfLO7kZ+LVr7CbgVWAzsAJYaIxxBpqzKR/YLqGSC+K9QnCh0qkBrjl7CHd9b0yjvqP7JVPtrOObkhZngCuK0kzCJX0dSI46kBR1MOysd9NNNzFjxgy2bt3Khg0byMzMDDq+3Upfi0g8VtWCc7DOWT4Rkb/aSak2xizHunTq3fY7r/eVWGcz/sbeDdxtZ04/fRJD7S0U9aG2wB6Pm1Dp1ABnDe/tt310P0uMacu+MjL7RU6YSVE6MpGSvg4kR90cKepQ65WVlbFy5UqeffZZAGJjY4mNrS/JZWe/bYWdWm3PAcepP7ifCzxPAIPRUYgOUjIH6guFgj3DE4ghvbsSGx2l5zxKh+ftt99m//79rTpn3759ueiii0L2y83NBeDxxx/nww8/bCR9nZeXFxYhOPcaOTk5FBYWsnDhwkaS060lRb1jxw5SUlK4+uqr2bBhAzk5OTz44IN07dq1RfOGAzu/MUcZYxYYYz50va4HRoZ7Y5Gmvjp1gFBbdNNCbYGIdkSR0TeJzWp4FCWsuKWvX3zxRb/S126dGvfLV8qgubjlqIuLi8nLy6OgoMDzrDWlqGtra/nyyy+58cYbWbduHV27duWee+5p6fbDgh2PZ52InGGMWQMgIpOBVeHdVuTxlMwJFGpr4PG07C+lzL7JvLflAMaYRto/itJRsOOZhItISV974ytH3dpS1Onp6aSnp3s8pzlz5rRbw2PnT/XJwGciUuRKVV4NnCsiG0Xkq7DuLoIEq04NDZML7KRTB2N0WjJHyqs5eLyqRfMoitIYt/T166+/HlT6+vjx462+dklJCaWlpQAeOeqMjIyQUtTTpk1jz549TVqrb9++DBgwgG3btgHwwQcfBD0ziiR2fmPOAIYA57peQ4CLgUuBy8K3tcjiTqNOjvdfHaipyQXBcCcVbN6r4TZFaW0iKX0dSI46mBR1c6WvAR5++GHmzZvH2LFjWb9+Pbfffrut/bY1IUNtxphv22Ij7Y0LRvfhiR/mMLBXF7/PE2LrjY2dygXByOhnZdJs3lfG+RmpLZpLUZSGRFL6OpAcdTAp6uZKXwNkZ2eTn5/f5P22NW2ndXqKER/j4IKsvgGfuz2e2OioFp/LJMfHMKBngiYYKEo7RaWvWxc7yQWKH+JchqelYTY3mX3rS+dU1To5WFbFgJ7+vS1FUdoWlb5uXdTjaSYJrWx4Rqcls/NQOU99soNz//wRUxd9xP5jIe/oNqDGWceqwrYre6EoitIc1PA0k3rD0zopmJn9kjEG/vDWFnp2jaXGaVi9o2lG5K2v9jHvqc8p2KPVrhVFab+o4Wkm7ns8reXxTBmRwg1ThrL4+jNY9rNz6JYQw+pvQh+KeuOu9/a5ymkritKOUcPTTOK9kgtag4RYB7ddnMkZQ3sRFSVMGtKTNTuaZkCKDlcAsPZbNTyKorRf1PA0k9Y+4/HlzKG92HWkgj2lJ22P2XW4HID8oqMBUzUVRVEijRqeZhLvqlbQWmc8vpwxtBcAa5oQbis6XEGXWAcHj1ex+4h9g6UoitKWqOFpJp4znhaWywlERt8kenSJYfUOe4antKKaYydruGiMVWo9X8NtiqK0U9TwNBPPBdIWVi0IRFSUMHlIL9bYNDzu850LsvqQFBdN/rdHw7IvRVGUlqKGp5m4kwvC5fEAnDG0J8VHT7L7SEXIvt+6zneG9u5K9sDurC1Sw6MorUVzpa8hsBy13fFNWe/BBx9kzJgxZGVl8cADDzR4NnjwYE477TSys7OZMGEC0I6lrxX/xEVHIRK+Mx6AM4dZqqV2vJ5vD1cgAgN6dmHi4J5sP3icYydrwrY3RelIhEv6GgLLUdsdb3e9goICnnzySfLy8tiwYQPLli3j66+/btDnww8/ZP369Z56bpGSvlbD00xEhIQYR9hCbQAjUhPp2TW2UVq1MYaFL37JioJ6Nceiw+X0TY4nPsbBhEE9MAa+3KVej6KAJX3trgIdHx/Pa6+95nm2aNEiPv30U3Jzc8nPz/e8li5dGnTOKVOmNKog7YtbjnrBggWAJUfdvXt32+Obst6WLVs444wz6NKlC9HR0Zx77rm88cYbtudvS7RWWwsY1KtrwOrVrUFUlHDG0J6NPJ79ZZW8tXEfJ6pqmTHGKmT67eEKBrn2kj2wO44oIb/oCOeP0mrXSvtg+/a7OH5iS6vOmZSYyciR/x2yX6Skr9tSjnrMmDH85je/4fDhwyQkJLB8+XJPSA2sP5YvuOACRIQbbriB66+/vtX3YBc1PC3gXz89m6gwK4aePqAHyzfu59CJKnonxgGwaY9VTDRv5xGqap3ERTv49nA5383sA0CX2GhG90smX895FMWDW/p6yZIlfqWvfUlLS2sgOdAc3HLUDz/8MJMnT+amm27innvu4a677mrRvP7IzMzk1ltvZfr06SQmJjJu3Diio+t/xa9atYq0tDQOHjzI9OnTycjIYMqUKa2+Dzuo4WkB0WEMs7kZm94NgK+KS5maYRmWgr1WLbaTNU7W7SplTP9uHDpR3cD7yhnUg8Vf7KLWWdcm+1SUUNjxTMJFpKSv21qOesGCBZ6w3u233056errnmfscJzU1lVmzZpGXlxcxw6O/kdo5Y/p3I0rgq+L6wp8Fe8ro1y2eKIHPCg95MtoG96p338cN6EZlTR3flJS3+Z4VpT0RSenr5spRN0f6GuDgwYMA7Nq1i9dff525c+cCUF5e7vl85eXlvPvuu7az6cKBGp52Tte4aIanJjYwPJv2HmPSkJ6MG9CdTwsP8a3rDs8gL49nTJrlKWmlaqWzE0npawgsRx1ofEukr2fPns3o0aO57LLLePTRR+nRowcABw4c4JxzzmHcuHFMmjSJSy65hBkzZjT3W9piNNR2CnBa/+58vP0gxhgOl1ez71glY9K6MaBHFx7/+BuPcRnk5fEMTUkkIcbBxj3HmJ2THmhqRenwRFL6GgLLUQca3xLp608++cRvn6FDh7Jhw4aA+21rwurxiMgMEdkmIoUi8ms/z+NE5BXX889FZLDXs9tc7dtE5MJQc4rIi672AhF5RkRiCDd5T8LuL8K+zLgB1hnO3mOVbNprJRZk9U/m7OG9cdYZlnxZTO/EWBLj6v+OcEQJo9OS2bRXPR5FaSkqfd26hG0lEXEAjwIXAaOBuSLiG9xcABw1xgwH/gLc6xo7GrgKyAJmAI+JiCPEnC8CGcBpQAJwbbg+GwA1J+HtW2HNo2FdBmBsupX3/9XuUo8hyUrrxvhB3YmPieJAWVUDb8fNmLRkNu0to65OK1UrSktwS19HQiY6nLgvkO7Zs6dJd4paSjhN3CSg0BizwxhTDSwGZvr0mQn8w/X+n8A0ERFX+2JjTJUxZidQ6Jov4JzGmOXGBZAHhDe+dHAzGCccbN17Cf7I7JdEjEP4as8xNu0pY2DPLnRLiCEu2sHEwdZ/LIN6Nr5PNKZ/Nyqqnew4pAkGiqK0H8JpePoDu72+Lna1+e1jjKkFjgG9gowNOacrxPZDwG9tCRG5XkTyRSS/pKSkiR/Ji31fWf8e+hpqq5o/jw3ioh2M6pvEV8WlFOw9RlZasufZOcOtsjp+PZ7+VoKBv3DbkfJqxt/1HisK9oVp14qiKP4Jp+Hxd7PSN+YTqE9T2715DFhpjPF7ymaMecIYM8EYMyElJcVfF3vscx3UGadlfMLM2PTurNtVyreHKzwGBeDcUdZnGNU3qdGY4amJxEZH+c1se2/zfo6UV/OvDWp4FEVpW8JpeIqBAV5fpwN7A/URkWigG3AkyNigc4rIHUAKcHOrfIJg7P8KEq0LnW0RbhuXboXNgAYeT0bfZN6/eQoXjO7TaEyMI4rMfsls9GN43tl0AIBPvi6h1lkXpl0riqI0JpyG5wtghIgMEZFYrGQB36p7S4H5rvdzgFzXGc1S4CpX1tsQYATWuU3AOUXkWuBCYK4xJry/SZ21cGATjP4eREXDwU1hXQ6slGo3WWndGjwbnppEVJT/0j1j0pLZtKdhgsHxyho+/foQg3p1oayylg3FpeHZtKIoih/CZnhcZzY/Bd4BtgCvGmM2icjvRcR9M+tpoJeIFGJ5Kb92jd0EvApsxjqrWWiMcQaa0zXXX4E+wGoRWS8ivwvXZ+PQdqithP450GtEm3g8I/skEh8TRd/keFKS4myPO61/N45X1bL7aL2mz0fbSqh21vHbS0YTJfDx9rbT4VAURQnrBVJjzHJguU/b77zeVwJXBBh7N3C3nTld7W13GXa/K7Gg3zjoMxqKA9zlqauDmnKIa3z+0lSiHVGcMzyF7l2adj3JfR60cc8xTwLCik376Z0Yy9SMVLIHdOfj7SXcPH1ki/eoKIpiBy2Z0xz2fQXRCdB7BKRmQukuqPJT5yn/abhnECz9ORxret0lX578UQ73zRnbpDEj+iQS4xAKXBWtK2ucfLT1INNH98ERJZw7MpWviks5Ul7dor2drHby8fYWZAkqitJpUMPTHPZtgD5ZEOWA1Cyr7eDWxv2+XQXRcbD+JXjodFjdssumIoI0UYbBnYr9zqb9bNt/nFWFhyivdnJBlqXjc+6oFIyxkgyC8f7mA6zfHfgs6KlPdjD/mTytDad0SFoifR1IjtqfFHUgdu/ezfnnn09mZiZZWVk8+OCDjfpUVlYyadIkxo0bR1ZWFnfccUfI8Sp9fapgDOzfaIXZwPJ4wLpQ6sv+Ahg2FX62FgadCe/fCZVt/4v55ukjOVpRzSUPfcJdyzaTFBfNWcN6AdYZUI8uMax0nfPsPFTOO5v2Y+V4WKz+5jDXP5/PD5/+3FMJ25e3Nlpp2f/6yjdxUVHaP+GSvg4lR+0rRR2I6OhoFi1axJYtW1izZg2PPvoomzc3/J0TFxdHbm4uGzZsYP369axYsYI1a9YEHa/S16cKR4ug6hj0c4W8ug+CmK6NDU91ORwuhD5joMcgOP+34KyCLcvafMtTM/qQ+6vzmJOTTtHhCqaP7kNctKU94ogSzhmRwsfbD3LzK+uZtugjbnh+LXe/tcUqSnqiil+8so6BPbsQJcLCl76kssbZYP6dh8rZuv84MQ7hra/2NTBaitIeiJT0dWvJUffr14/x48cDkJSURGZmZiPZBBEhMTERgJqaGmpqajwREjvj2xKtTm2HyjI4uhP6jq2/ONrXZXiioiA1o7HhObgFMNDX5YanT7CM1MbX4PR5bbZ1Nz27xnLP7LHccO4weiXGNnh27sgU/rVhL8sL9rHgnCFU1tTx1Kc7OVnjZN+xSo5W1PDMf0xk/7FKFvwjnz+8tZk/fO80z/i3XdUPfj51BIve28763aWcPrBHm34+pf3z318XU3DiZKvOOSYxgbtGhK6OFSnp62By1M2Voi4qKmLdunUecTlvnE4nOTk5FBYWsnDhQr99go1vK9TweFNdDivvgzMWQqJXVYMP7oQvnrKMTZdeIA5I9ap3mpoJ299pONf+jda/fV2/oEXgtCvg0/vhxEFITA3vZwnAkN6NS+vMzLbc7PNGpdA7MQ5jDF3jovnrx98A8PuZWWSldSMrrRs3TBnK31buYPKQXlw2zhr39sb9jBvQnflnD+bh3EKWfbVPDY/S7oiE9HUwOermSFGfOHGC2bNn88ADD5CcnNzoucPhYP369ZSWljJr1iwKCgoanEGFGt9WqOHx5rNH4NO/WKGzc//TajMGtr1thcxqq2DHh9DnNIjxUjJMzYJ1L8CJknqDdaAA4pItL8fNaVfAJ/8Lm96AyTe03ecKQYwjijlemj0iwq0zRpGSFMfBskp+eEb9Z7jlwlHkFR3hN29sZMLgHtQ6DRv3HOO2izJIjo9hysgU3vpqH7+5ODPgpValc2LHMwkXkZK+hsBy1E2Voq6pqWH27NnMmzfP47EFonv37px33nmsWLHCY3iaMj7c6BmPm+MHYJUrU2Sr1znMgQIo2wOTfwL/sQZ+9CbM+mvDsf4SDPYXWJlv3lloqRmW0dr4Go0wBoo+teQWAnFkJ5S1zeG9iLDgnCHcdnFmg0y6GEcUf7kym9o6wy2vbWC5K6ngojH9ALhsXD/2l1WydtfRNtmnooQiktLX4F+OOpQUta/0tTGGBQsWkJmZyc03+68IVlJSQmmplXl68uRJ3n//fTIyMmyPb0vU8Lj56E/W4f/4H8G+9VDqKoLtDqGNuMA6zxl6Xv25jZu+pwFipU+DdXH0wCbLS/LltDnWhdMjOxu2f/EUPHsJPDAWPnvYCvt5c/RbeOJceGYG1FTWtxtjidHVNTzwDyeDe3flvy8dzarCwzzw/tdkpSUz0CW7PS2zD3HRUSzboNltSvsg0tLX/uSog0lR+5O+XrVqFc8//zy5ubmeJAl3GNC91r59+zj//PMZO3YsEydOZPr06Vx66aUhx0cEY0ynfeXk5FgCPge3GvM/PYx56xZjDhUac0eyMasft549+V1j/jrFhOQflxvzlzHGOJ3GHP7GmiP/7437Hd1lPfvgrvq20mJj7u5vzFMXWPPckWzMfSOMKfrMel5TZcwTU425q4/17OM/149d81er7bNHG6+zdbkJF3V1dWbBs3lm0K3LzCO5Xzd49pPn8032ne+Yo+VVIec5Wl5lXsnbZb7YedhU1tSGa7tKhNi8eXOkt9Aq7Ny502RlZbXJWhs3bjS//OUv22QtN4MGDTIlJSV+n/n7GQL5pgW/ezu3x1NeAp/cD28uhNiucO6t0GsYpGRa4bbyw5Z3MvLC0HNlz7MqGHy7ygqzgRVW86X7AMi8DFb+L3z5nOWxLL8F6mqtEN6P3oRr3oXYRHjucuvy6Qd3wp58mPW4NfaT++FYMexdD+/+1pp37d+tudy8dTO8fBWUbGv598kPIsI9s8fybxMGcMWEhrH7n00dwfHKWn6/zM/dJi8KD57ge4+u4r+WfMWcv67mtP95l5+/vE4VU5V2h0pfty6dO7ngWLH1Sz0qGi6+D7paompkXIL59H7K1j9HNwyMsGF4Mi6F2CTY8DJ0SweJqj/78eX7T8IrP4ClP4Odn8C25TD9Lug5xHo+cDJc+z68Nh/+70arbeK1kDUL0sbD1+/B8v+Cki3QpbeVqPD+HbBrjXVR9eAW+Ppda9ynD1gGKwz0TozjXj8lfEanJfOTc4fxyIeFXDYujfNHpXKy2smLn39LVW0dw1MTqXHWcdvrG4mLjuLZqydSVVvHO5v28/qXe5iTk86UkS3QSlKUVsYtfd3RcF8gbWs6t+FJ7AO/3QmO2AZJAOv6jeR/+6aw5eun+WdyKkPTTg89V2wXGDMLNi6x7uz0Gm61+SMmAa56CV6dDxtftdK0z/iPhn269IQfvG55NEd2wAWueqk9BsHZv4CP77GM21dadu8AACAASURBVI/fsqoorPxfWPusZXhWP2LVkht9uTX/+bdB94HN+x41k59NG86KTfv5zesb+d1lWfxx+RZ2Halo0CejbxJPzZ9Aeg/r+3TeqBQ+3lbCc6u/VcOjKB2Yzh1qi0mwaqm5jE6Ns4ZbPr6FH33xB/bFxmKA1/sNs5IK7JA9z6pGvfNj/4kF3kTHwZXPwXfvhCueBYefvwEcMXDRvTDvtYbp22ffBIPOhgv/CIPOssKEY6+Azf9nhda+etW6pDrtd4BYyQotZfs7sHed7e5x0Q7unT2WfWWV/OSFtTiihJeum8zG/7mA/1t4Nn/9QQ5LbjzLY3TcY66aNIDcrQcoPloRZHb/GGN4d9N+lqwtbvJYRVHajs5teKITGny5et9q3il6h6uzrmZZ6gWcV3GSpXVHqXHW2JtvwGToOdR675v55nf9WDjnF9a5UlOI7QJXL4czbqxvy/mxpRH04hXgrLE8qG7pMPbfrLOk4wdg61vwj8tg9WNNW293nnVetHhew4y6EOQM6sFdM8fwnxeO4u2bvsNZw3qTFB9D9oDuzBjTl65xjY3t3EmWZ/bS57uatMXtB44z76nPuf75tfzqtQ0seneblu5pJ+jP4dQlXD+7zh1qi24oqJa3L4+YqBj+I/s/iD+6i1ll3/B+9Q4+Kv6I6YOmh55PBMb9O3z4B/+JBUHYVbaL+764j5vG38TwHsObNBawwm39sq1U8MzL6o3ZOb+A9S/CIxOtGnMxXaxzpZRRMHya1afOaSUq9B4B8T63mavL4Y0brMuwZXsg/xk40ycsGIQfeF0+Zf9Ga52TR6z7SuOugh6DG/RP79GFqRl9eOWL3dz03RFU19bx91VFrNt1lF1HKth3rJIx/btxyWn9OHNYL9Z+e5T3Nh/g4+0lJMZF8/uZWWzeW8bDuYVU1ji5/eJMDp2oZteRCkb1TSLRj7FTwkd8fDyHDx+mV69eTa6srkQWYwyHDx9udPepNejc/xf6/I+Qtz+PcSnjiI+Oh5SRnP1vr5O65AKWfL0kqOE5dPIQ1717HZcPu5wfT1yAOKthSPDSF95UO6u55eNb2HJkC98e/5bFlyymS0yA86FgTLwWlv4Uzrqpvq33CJhwjZWdd9bPYOQMeOZCWHIt3LDSknZYcq2VjRcVY50RjbzISmRI7gfv/rd15+jHy+DjP8Mni6y7TnGJDdfeuRI2L7UyAIdNteb1Ztfn8PeLwHjdN/rsESupY9xVDX4WPzpzEO9vOcBtr29k5fYSDp2oJrNfMiNSkzhrWG/W7DjMHUvr5cbTeySw4Jwh/OTcYfTsGktdnSEuOoonP9nJS5/vorzaWrNHlxh+cu4wfnTmYBJiw3NLXWlIeno6xcXFlJSoVtOpSHx8vKfSQmsindkNnjBhgnGXIz9WdYzvLP4ON467kRuz60NYD697mCe/epJ357xL3659/c5z3xf38dzm5wC4cuSV3Db5NqKj7Nv0e/Pu5YUtL3D1mKt5tuBZLht2GXef00h8NTTGWNWz3dlxgThUCE+cZ6V2H99vlQKa+hs4cQC2v2tly0kUpE+C3WvgzJ/ChXdDcT48NQ2m/hamuEoK1Zy05B4+f9waY+ogub9lBM/6mXVOVX4Y/vYd6/28JZDUByqOwBs/gV2fWUZuxr1WO1BXZ5i66COKDlcwaXBPfnNJJuMGdG/wEb7ed5SiDSsZMO5cRvXt1uivaWMMz35WxDclJxiWkkif5HgWf7GbldtLSEmK46IxfZkyIoWJg3uCQI2zjvgYRwOPyBhD0eEK+ndPIDa6c0elFcUbEVlrjAkuIhSEzu3xeJF/IB+DYVK/SQ3aZw2fxRNfPcEbhW9w47gbG407dPIQr257lUuHXkpql1SeKXiGfeX7uH3y7aQnhf5LYWXxSl7Y8gJzM+Zyc87NxDvieXzD40zqO4mZw2c27UOIhDY6AL2Hw/ceg1d/aCVBXPGs5RkBTP89HPraSlDY+Br0G4fzvNtZ+vUb9E/sz6RRl8Cqh6HbACt0tm25lXU3+Sdw3m2w4yP48h9WmvqWpTDrCXjndqrKS7hl/EVMPbqRWb2HW3LgP14Gqx6AD/8EhR9Y96gm30CUI4bH5uVw4Hgl541MaRyiqalkxEcLGbHtLUj4DfT7Lz/fCuHqsxt+Ly4+rR95O4/wxMod/HNtMc+t/rbRt29UnyQmDO5BZU0dH28voeR4FePSu/HCtZNJim+a7LiiKP5Rj8fl8dyTdw9Lti9h1dxVxDoaygZc++617Dy2k9cvf51ucd0aPLvvi/t4YcsLLP3eUgYlD+LVba/yp7w/UWfqmDZwGvMy5zE+dXyjX56HTh5i8dbFvLjlRfon9ufFS14kzhGHs87Jde9d9//tnXm4HUd14H+nuvvu7963a3uSJdnajG1Z2NjBYBsvgO0QyEAGjE3CBL4EAvnI8EHGQJgsJEwyM86MIUOAhGRCDIHEC2AHL3iceMU23hdZtiTb2pcnvf29u3R31Zk/+sp+1mIJZJ4sUb/v6+/e7q7qrtPV3afrVNU5PLjjQfqKfcytzKW/1E81V6Uj18Hi2mLOHjib3mIvQ40hrll7DbdtvI2VfSt5z5L3sKJnBfdvv5/vr/s+G8Y3cNqs0zhzzpmc3Hsy3YVuRISR5gi3b7qdBzfczruWX8pZC8494DVaPbSaL9z3BZ4eehojhiuW/waX3fSngEKQz/qWzv+DzJXQdJ7+Adz4e1ngO3Vcdfq7+buhh8iZHP/yK//C8Z3TBlTsXg+3fAbW35YNzjj9Q7Dy/S/Nq5pOcxy+exlsuDsbhr7zKfjgjbDwzS+lcQ4GV2d9WflKZl6c7g1clZZ1PLxxhKe2jhEYQy4QhqcSHto4zKObRgmMcPaSXpb0d/BX/7aOVQs6+eaHzqCU899qHs/htni84mkrnnff8G56Cj387dv+dp90jw0+xm/e+pus6l/F1y/8OlGQffnubuzm4usu5m0L3/Yy09jOqZ1855nvcM3aaxiPx5nfMZ93LH4HJ3SewHNjz/Hs8LPcteUuUpdy7vxzueINV7ysdTTaHOWatdeweWIz2ya3MdgYZDKeZDwep2VbACzrWsYLYy8Qu5hT+k5h7fBamrZJR9TBRDJBNVdladdSntz95It58kGe/lI/2ya3YdVSDIvENuYP3/iHvHtJ5q02tjGP73qcx3c9zmODj3H31rvpLnTzydM+yW0bbuOOLXdw+cBb+fTyXyecfVJmPtuL4eYwQ40hlgQVuOnTrC6UuHz8Qc5fcD4P7XiIuZW5XH3J1URmr7zP3pKFjdj8QNbfNHdVpjj29HfZOGuNjW2GX/0aLLsoMxm2JuF37oX6EDzw9Uzp1aeH8RUYeEM2UXh0Y2ZeLPVkpsa+5ZnZcFpL0TmFtInZ/QzsXM2TG3fyrQd3MDCrhxPOvIRx00mcOgJjiAKhryPPm07oJQpebo5TVaZiy8hUzJxagXDX6syL+VmfgNo8AFqp5YdPbKeZON57+gBh4E16ntc+XvEcBnsUz3BzmHP/+Vw+seoT/NYpv7XftDc+dyOfu+dzvL/rrXz6tE+R9ndx5UNXct26615s7exNPalz28bbuPH5G/nJ9p+gKIIwrzKPN817Ex9Y8QEW1hYecnlVlbUja7lzy53ct+0+ju88nstWXMbi2mLG43Fufv5mHtv1GOcOnMt5C84jH+Rp2RaPDz7OutF1bJ/czvap7SyoLuDtC9/OQGWAT9/5ae7ddi/vXfpexuIx7t5yN/U0m0OzsLqQcwbO4SMrP0I1V8U6y5UPXcm31nyL/mI/Fy26iEsWX8LyruUEJsA6yzVrr+HLj3yZiWSCdx7/Tn731N/lY7d/jPF4nO+/6/vct+0+PnXnp/j4qR/noys/+jL5Rpuj1NM6c+tj8MjVsPNJiOtokjlMlSCfzVl68ydhSXuwx/Yn4BsXQrELJndkrbAV74ATLswGeDRG4JmbstZUkMsm0nbMhqndmeeKLW0Hq+d9Dk69LFNaT/xz1p+l+zpeTTTgdvd6brWnM0928zqzgS6ZZF24hMoJZxEufCM/3hHw8MZhNg3XaSaOCnWuKHyPy7kZg2OqspAHzvkmT46Xufr+jeyezD4MThmo8Zf/cSVLZnWAs6SbHiDsPWG/sZtUlXvXD7FxeIp3rpzrzYCeGcUrnsNgj+K5dcOtfPrOT3P1xVdzav+BfTH907c/x7L/8T1sFPD7Hw4ZLlneveTd/MlZf3LQc+2c2snu5m4WVRf9TCPW1Fp2fenLmFKRno985FUbmpq4hC/e/0WuW3cd3YVuLph9DufJcl63/By6++fvN88dm+/g+nXXc/fWu0ldSjEssqxrGfW0ztqRtZw550xWdK/gW2u+hapi1fKVC77COQPZSL8r7rqCH234EW+Z/xa6C90YMTw6+ChrR9aiKOfPP5+Prvwo3YVurl13LdeuvZax1hid+U66Cl1ceNyFXLb8Mmr5Gs20yU13/leefe5WLjru7Zz6pv+CVPqygQHjGxCEBdUFGDlAS2JsK/zwU7D25pe29a2AZRdnZsTZJ2f9UWmTXYPbCFdfT3XttQSNoez61RZRDzooDa8hIpvv9QJz2Vg5lWoxx+x4E71T6wjTSa7lQm6KT+X/RH/FDu3m/fHned2yJXz4zYsYb6R8/ntPUE0G+Ujng5w3eRNz2EWTPBsWvY8Fv/JZSt1zaSaWRzYM8bXbHmf9pm1Upc7cfINfXV7mzIU1OjpqFMtVpH9F5v3ip8A59TGUPIeEVzyHwR7F82f3/xk3PHcD977/3n3NP23Gb7qJbVd8hvHeIsVdE4ycNJ/wyj/kjLln/lQj2ADs6ChSLGLy+YMnBjSO2XrFFUzcfAsAtV97D3P++I+RMCTZuZOJW26hcPLJFFetQkRoPfccu666isbq1VTOOYfqRRdhSiUm77iTyXvuwdWnCKo1gs5O8suWUly1iomuPObWuxm/7nrsSBZLx9RqBJ01SC1qLfmlS+i67DIqZ58N1rLr/rtZe9cP2OZGeD4YYbAYc8kZv855r38PQanEhrENXPXIVcwvz+NjpUuoP/ggrefW09jwPIObnuWZxTluXSU835tySu8pvGH2G0g15dtrvs1EPEEgAU4dZw+czfGdxzPWGmPrxFYe2PEA5ajMefPP456t9zDaGiWUkFRTlnQtYUnnEh7a8RCDjSwOSikssbRrKeWoTOISBGFl/0reMvAWlvcsZ+PYBlY/+W22Dz1L2H8iQXUeu5q7WTu8lnWj6ygEBQY6BhjoGKC70E1n1EF33GBe/ynM71lOb7EXsTHD63+C3XAfhd0PcPfux3AiHF9ZwHG9y9i14pdZE+V4dPtz9NVHOPWRv6Y/rJDrXw6FTgpJg8rWxwgndwCwtnwaq2ddRGHzj3lrfAcWw6agwIa8MBLA0iRmWZxQPMDzazGsK5zCtjkX0Fr0VnrmL6O/I89wPWZwrE7SrDO3v5eFPSU2jzT4/v1reOHJ+1jVMcoHXlegV8ag2IWrzWcbfQQ9i+mafRy5MGRwosWWkTr5MOCkeVU/P+cXkNe04hGRi4AvAQHwDVX9i73254F/BE4DhoD3qeqG9r7PAh8GLPAJVb31lY4pIouA7wLdwCPAr6tq/ErlO+200/W2e3/E5T+8nBVmLp8fP4fJu+4mt3AhpTPPIL94MY3HHmPqxz9m/KabKZ5+GvO/8hXGbriRnV/8IrM+/3m6P3A5LzyynQeuXUOHG2X+6MOUNz9OMGs29dnLSaqzKBcd1ZKlvmWQTWuG2Zn0EkSGgRP7OP6958H4CKN33cvkk89iymXCvj6kp59GqY+6dDD50MMUVt/DossuppCMMfTVr1K58AJc12w23fU0o+XjiHMd2GovVGqE29ZTsmN0zqsRPnEPxZFNGHVgDIWVK4l6e7Hj49jhIVrPPZ91xgMYQ8cF51O54ALs0DDxls248QkkDMEYJu+9l2TXboI5c9GxMXRqApWQiY75jNUWM1WaAzhEHREpFRmnGkwhO7dgG01UAqSrh2DOXEy5g9bjj2BaUxROPBE541zSBSfiyjUKXY77d/2QcNsWTtnVj31hiGB4J8WxzVCfwlYKbC45Xogcs3OwuDhAhymzsznI5vpWWmmLPlOjU8okPT08P1BjdZcyUhsniFJsHDO4dS21SUf3lFCuK5WmkgTw1HHCswNCmC9yUm4hpzb6MFMNRscGmZwYQsam6Gh2ErkqTyzYypMLUwq5Mid0ncDSrqWMNkZ4eP1dFCdjRsrQKBz8pSyqiIIxhu6wBGGBsWTyxb45g8GoksrLn9UA4YTiLE6uHc/cwgmk9Spxo0HSmKI2vIalk49xnN3G1ijkgbCXx8NO0miKZtjEiXJSw7KoUaSslq3FKR4p5IlF6LOWbid0pTEdzr24FKxhl3SwOh+ypiA0jLC8VeGs4iKW1HoZclMMa51Kq8H8xjhzGuOk3SuI558Fc05lpLmNbRPPM56O09HRT6HST3+uxuJWTGvn84S2QSE0Waur1Jt53uiY03ZtVcg8uI9tyZZkKpvUnK9mHjumdjEyuY18vkapNh/K/S/NJTNhZl6tzNpvnySQPQNpM+tHzFcP3U3WLyivWcUjIgGwFngrsAV4EHi/qj49Lc3HgFNU9aMicinwH1T1fSJyIvAd4AxgLvD/gKXtbPs9poj8C3C9qn5XRL4GPK6qr+iWefacZfre3/oCk7lRTl6/htOfvJfcnD4aE8pIeRlTHYvpSSboy0NxwTyG5nWxcevTOE2pDe6iZ+1TDK98H5vdIgqNYfLFThaXisyJYNxadjW3s6u5iXEtkZjZiOlGTEAhTHGpJWbfFo9qiqSbyJtRIqoYKaHqaKaDNNNt5IzQneuhGtSYsDE7WsM4EXJRhSgNCWJolUu00iZoEzCICFEg9OS6mRN1YbXJzuYGBhvrCAq9lCvHU4z6aUpKc3wrLhkn0QCRAkiIOoc6h7MTqBtC7ShiKkg4BxP0AwKaEEUJzk7hbAPnFDGdBKaGISFwQwQ6hiKoFNtLGaSMSpHU1bFuArTZ7guDQEqYsB8T9GF0isBuJdJBDDkwHRhTwjlLvjVMJdmJ2AbqWiSmRKO4iPH8UlqmiLopxE1iUDqaw1TrgzgXkyDEJiTJ1UjyPaBCZXw9leYmJFSmKFAvVFHTgTGdaNRLs9CPC0LQFElHKE2tQ8PdaJpA6gisEqVKZC05p0RdNZoDnRgJMRMtkkYLS4QzeazLkSaWtB0GIh+C7XIkJSFMQGKLqSeYqRTTtES0KOct5Uqe8VzIqBEm4xSpK1GSB5MjDULSMMSooxhPUmyM0gxH2F0eZqQ6jimViYo9OIo0xkbomEjIJYoLDR2FKq5UYrPGjEQt8pqnnNQopGXUDYPdAW4EF5aohl0YzTFuE/KJIzUNhku7Ga40SQ3kXIW8rYJO4mQcUci1AqqtLvK2gnEGo4YkSGnlpqhFY6SRMBR1MBaWiJylksaUXAsNW7hckzSf0Mgr9ZwgQGfsqLYc9TRgp41IE4MKiFGKRtGcw+UUCZXOWKm1HLU4pK+hdLUgBrZEwtZ8DiMxvZLSjcUYgyv04aJZVKIWnWETwbHZxWy1KVZD+rVGp/SSz9eIygVcLiQJLXGUYE2LqDFJvj5JVG9g0xibpuRNwPxqPwu6j6NQ7WM8DBg30HAJsU1INKVgclTCAqWgQDEsU4pKSFikUahSz5dBchQSIWo5Gnlll5lkd3OIlo1JURBDf3k2czoW0FOehbgUcSkOJQZaAuOtMUamdjLaGKIUlekp9dFV6qMQlckFeUITIgBJPZuXl6uwR0vsad2+lhXPG4E/VtW3t9c/C6Cqfz4tza3tNPeJSAjsAPqAz0xPuyddO9s+xwT+AtgFzFbVdO9zH4ho2Yna87V/enG9nDZYNbmevniSYirkrdBCmVQllZRc0CQKmuQ0JWyBNFOaQZ7JYgfj+TITYY2xsMpkWKRoE2qppStN6IrrVFuTBGmDySiiHoaEJqZDxukIxohSReI81kYMF4ShYp5mGNHdarBgCnrigNFiwnjBkRpH6BKMS0ldSKx5mhrRCvI0gwKpCSnaOqV4ElHLSLGLwUIPqQnp0jFqjDImVTYHA2wPe6nZKQbiQebGu4lcjEtjrArNsEgrKhObCMFi1BKqJecsBeeInBJYCFRoBiH1MEfLBISmQY4pIlqQBmgagAsQURAIsYQaE2mCC4U0EmxoCJySc5YARxoqaQhWBOdyqAsxzpFzMTkXk5CjIQVakgckay2oUrYxFdcilBZp3pJEjpbkaLoyTS2RuByphjgNCV1KziXkXULJxZRtQk4SXJhio4SmiZjQChNSJHBKSVuUXEJKQExIU/JMmiJTpkRqAorUKTNFJDFOAxwBkbOU0xaFNCExAfUgRzMIybmUgovJaYLVEKsRikEytUwoKbmwSS5okrgCU7ZGQ8vExmADRQXyNqGYJuSsfdHpg4pgJcASEAcBLROQmIC8tii6OgXXQI0hNSFWDOIcoorBEhlLSAoYUhvhXERKli7FkBKSSEgqAWSvJQwQYYlICTTFqGsv2b0hzqBicMaQGoMzihVFxJF3MQWNMU5xhDgCRCHAYnBYhBQhwYAIarKzRpqQcymRpgROs7QS0JIcsUSggnEQOM3KhAVRmlKgIQVSCcmpJacODDhRnMnun1xqyVmHOABB1dAKDc3QYMWQ05S8TQmdItaAFQTFqEWwOAMaZMcTNaAGcYIgiCooWFFsu8KMmkxmVcL2M6ZOsRacGlxo0DCAQAmcJXTZNc5OJG1loKgoKoKKvOQBZM97vX0vGjSrJ7EoitWsXq0YUglwIuStJZdawlRJg6zOrJGsTtQi4kAURPnib//Ra3YC6TxgegCLLcCZB0rTVhhjQE97+/175Z3X/r+/Y/YAo6qa7if9AZmt2/kL/ThVxnmGE7k3OJuf1N5ALIfW9zKdQFO6GaKbIRbyAg1KTFBhC52MMBeVn85FS6QxSUcO9jOV5UCIZi+B6eUXtcxmBwVtsl0WMUaNMlMscBtZro8yGnSyubSAx0tn7lPGik5QoIHDYAlIiEjIkUhu71OT1wZ5WsTkaFFAD9SZ/yoi6pB2++iVzifqyNMiJCXAIigxORoUDloveW1iCUjl5SaaQFPKTFJhghwxw8yiTpmEiACH4Ij3c63y2iQm91Nfn5y2iIgJsYDSpEhLDuxDK9IWeVrkSGiSpy4vd3EUaIqVQ3v8jVoiYnLEhKQvKkjXVkgxeSyZwjuQXKLZXWTa9bX39fQcbfzRYeX+eSqe/Rm3925eHSjNgbbv765+pfT7Fkrkt4HfBlgwMJ+zer5ElNvFm7o7+XjHHKKom0YaM9wYZ7w1QT4CXOaRWcIeJOomVWGiPsbY9kE6ImFOf5XeShG0G+cGsK6JOmg0WrRaLSSqM2pyxAi10FAJQMIqDdPLhBZo2oR6PEqcjjNQ6mJRxyxKQcjOVoPVQ1sYnJhgXnUWc6vd5MOQhnU0XPZ1maNFoA1qYUQ1ijASMp402N5s0HKOpR09VHIvTWB1mpmxRF7yJaeqWOeo25jxeJzUtujNlylGFYzJvZg3C1ubkrqYlm3RSmNaNqYjiigEEcbkCYLsBRerktiUejyMdU2MyRFISIqhYaGlQjkqUw1DCoEhcUrTKVaVghEKgSESwbTPHTvHpHVMWUfBCJUgoGDkxbJZVSZSy1hqiZ1SCgzFwFAyhvy0dADOtXAuRdXSsI66UyadkGhIIYiIjFA0hmoYELVHebWcYzJ1L+4LBVQTnMv6YYKgst9O9oZ1jCYpxcDQEQYEIqgqdetoOCVA21/5Nru9RXCEJAqxOgrGUDYGo1MYk8OYlz4qmmmTqbSZmfEkIJCASCAUxZgImaYEUqdMWUvOmJddN6dKokqqSuoUESHInpPs+qOYaX1LqtrOmy17y+xUadmYRjJFK62TCwKKQYFckCMICkhb0SdOqVtLohDJSw+2be8LjRCQtWhCExIZgyrUnaNuHS3ncAoWJRKhIELOZPdnbGNiG6MSkFiDWqUcOUomIRJLSkA9dkBAKVciFxWIbcpofYLxxhRqHCYUMErFRBQ1JLBCS5SGdTRti9S2iJM6DkUlQE1AoVAmnysQSgRisWlCmrZoJQlxHKPOUsoXyYcREghWLYnNLAyJGmKnkCa4pImNYyqlMqVymSCIiC00Y4t1FieK1Ti79lawqSWUgADAOiQIcKqkalGxpM5iVUnIkbogM2OnMWJb5KMIg0GdpelSJtOElrXkRAgBo+AwOCeoBKQEpDbrfD8cfp6KZwswfTzuALDtAGm2tE1tNWD4IHn3t3030CkiYbvVs79zAaCqfwP8DWSj2hat3Le1GEVQLR6kqVHtgtkLXzFJV9crH+Ll1PbZMrtQYva8pftJO53KPlt6oio9BxixbfbzchQRwiCgGhSp5or7yfVSOpGInInIhWU6XqFhmBchb3JUov37t9snvdmfJC+RM4ZuY+g+wIdyIEJnFNIZHfyWNib/4gs8iqB6kPRZ+Qz53Mu/e0RyGLNv6286xcBQ3MsThohQDgPKL245lK//jn22FMJC5tD2EAiNUNvP6EsjktUVZMN19uHl98vBBrAZEYphnmKYJxvns3+iA5TnFRGomoBq+Eqt1DxMu7IHorbX81GM8tQKB8/nyThcxfPztIc8CCwRkUUikgMuBW7YK80NwAfb/38N+DfNOp1uAC4VkXx7tNoS4CcHOmY7z7+3j0H7mD/4Ocrm8Xg8np+Rn1uLp91n87vArWTfUn+vqqtF5AvAQ6p6A/B3wNUisp6spXNpO+/q9ii1p4EU+LhqNo18f8dsn/IK4Lsi8mfAo+1jezwej+c1hp9A2vbV5vF4PJ5D43CHU/tZUh6Px+OZUbzi8Xg8Hs+M4hWPx+PxeGYUr3g8Ho/HM6N4xePxeDyeGeUXelSbiOwCNh7pcrzK9JJNqD2WOdZl9PId9e2o8AAABgFJREFU3Rzr8gEsU9V9ZzUfIr/QAeRVte9Il+HVRkQeOpxhjkcDx7qMXr6jm2NdPshkPJz83tTm8Xg8nhnFKx6Px+PxzChe8Rx7/M2RLsAMcKzL6OU7ujnW5YPDlPEXenCBx+PxeGYe3+LxeDwez4ziFc9RjIjMF5F/F5E1IrJaRH6vvb1bRG4TkXXt358qMtBrDREJRORREfnX9voiEXmgLd8/t0NkHJWISKeIXCsiz7Tr8Y3HYP19sn1/PiUi3xGRwtFchyLy9yIyKCJPTdu23zqTjC+LyHoReUJEXn/kSn5oHEC+/9m+R58Qke+JSOe0fZ9ty/esiLz9UM7hFc/RTQp8SlVXAL8EfFxETgQ+A9yuqkuA29vrRzO/B6yZtv7fgf/dlm+Ew49LdST5EnCLqi4HVpLJeczUn4jMAz4BnK6qJ5GFM7mUo7sO/wG4aK9tB6qzi8niiS0hi3z81Rkq4+HwD+wr323ASap6CrAW+CxA+31zKfC6dp6/FjlIPHm84jmqUdXtqvpI+/8E2UtrHvAu4JvtZN8EfvXIlPDwEZEB4JeBb7TXBTgfuLad5KiVT0SqwDm0Y0epaqyqoxxD9dcmBIrtKMMlYDtHcR2q6l1k8cOmc6A6exfwj5pxP1mk5DkzU9Kfjf3Jp6o/akd3BrifLMozZPJ9V1VbqvoCsB4442Dn8IrnGEFEFgKrgAeAWaq6HTLlBPQfuZIdNlcB/wVw7fUeYHTaQ7CFTNkejSwGdgH/t21K/IaIlDmG6k9VtwJXApvIFM4Y8DDHTh3u4UB1Ng/YPC3dsSDrh4Cb2/9/Jvm84jkGEJEKcB3wn1V1/EiX59VCRN4BDKrqw9M37yfp0To0MwReD3xVVVcBUxzFZrX90e7reBewCJgLlMnMT3tztNbhwTiW7ldE5A/ITPzf3rNpP8kOKp9XPEc5IhKRKZ1vq+r17c079zTn27+DR6p8h8mbgHeKyAbgu2TmmavIzBV73D0NANuOTPEOmy3AFlV9oL1+LZkiOlbqD+BC4AVV3aWqCXA9cBbHTh3u4UB1tgWYPy3dUSuriHwQeAdwub40D+dnks8rnqOYdn/H3wFrVPV/Tdt1A/DB9v8PAj+Y6bK9GqjqZ1V1QFUXknVg/puqXg78O/Br7WRHs3w7gM0isqy96QLgaY6R+muzCfglESm179c9Mh4TdTiNA9XZDcBvtEe3/RIwtsckdzQhIhcBVwDvVNX6tF03AJeKSF5EFpENovjJQQ+oqn45ShfgzWTN2ieAx9rLJWT9ILcD69q/3Ue6rK+CrG8B/rX9f3H75l4PXAPkj3T5DkOuU4GH2nX4faDrWKs/4E+AZ4CngKuB/NFch8B3yPqrErIv/g8fqM7ITFFfAZ4DniQb3XfEZfgZ5FtP1pez5z3ztWnp/6At37PAxYdyDu+5wOPxeDwzije1eTwej2dG8YrH4/F4PDOKVzwej8fjmVG84vF4PB7PjOIVj8fj8XhmFK94PB6PxzOjeMXj8RwBRGShiPyn9v8PiMgVh5hvclr+pw6W3uN5LeIVj8czw4jI7wC3An8qIncAdwKHFMfE4zkWCA+exOPxvFqISAfZTP5fAVYAdwCjQEFEqtp28ioivw+8l2yW//dU9Y+OTIk9nlcf3+LxeGYWB+SAKoCqbtAsltLtZH7MEJG3kfm8OoPMpc5pInLOkSmux/Pq4xWPxzODqOoU8BvAfyMztV0pIiWy+CZ7oj6+rb08CjwCLCdTRB7PMYE3tXk8M4yq3iAiT5CZ204HPkWmiP66nUSAP1fVrx+hIno8P1d8i8fjmUFEpCIix7VX94Qr71BVCzzbjmF/K/ChdoA/RGSeiBy1UUg9nr3xLR6PZ2aJgK8DvWSu9DcBl7X33ULmVv4vRWQFcF8WwoZJ4AMc3QHhPJ4X8WERPJ4jgIgsBN6iqv8wbdsc4Juq+rYjVCyPZ0bwisfjOQKISCewUFUfO9Jl8XhmGq94PB6PxzOj+MEFHo/H45lRvOLxeDwez4ziFY/H4/F4ZhSveDwej8czo3jF4/F4PJ4Z5f8DEb/N/XiJb+YAAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "for (z1,z2) in zip(zshell[:-1],np.roll(zshell,-1)): \n",
+    "    galmap=tmap[z1]\n",
+    "    #compute mean\n",
+    "    Nmean=np.mean(galmap)\n",
+    "    print(\"bin=[{},{}] mean #gals={}\".format(z1,z2,Nmean))\n",
+    "    #and do fractional difference\n",
+    "    dens_map=galmap/Nmean-1\n",
+    "    #compute power-spectrum on the spehere using a standrd healpix function\n",
+    "    cl=hp.anafast(dens_map)\n",
+    "    plt.plot(cl,label=r\"$z \\in [{},{}]$\".format(z1,z2))\n",
+    "\n",
+    "plt.xlim(1,120)\n",
+    "plt.legend()\n",
+    "plt.xlabel(r\"$\\ell$\")\n",
+    "plt.ylabel(\"power spectrum\")"
+   ]
+  }
+ ],
+ "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.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}