{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import relevant modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import kn_calc_functions as knf\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic example: fluid Love numbers $k_n$ for a constant-density planet\n",
    "Create a constant density (i.e., single-layer) planet by generating a random number for the radius and one for the density.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n   k_n\n",
      "02  1.5000\n",
      "03  0.7500\n",
      "04  0.5000\n",
      "05  0.3750\n",
      "06  0.3000\n",
      "07  0.2500\n",
      "08  0.2143\n",
      "09  0.1875\n",
      "10  0.1667\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAESCAYAAADTx4MfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xt0VfWd9/H3N6SIISJXQYUkjAjE\nWlGjrKpgAaMgKnCi9TIZa7WajtpWO7ZemjVabaPzzGrVzmi1jNRLjWJRJIAKiFzapxUQBhzAC9HI\ntRZGUWISEYjf548c8iQkQDgmZ/9O8nmtddY5Z5+dcz4JkA/7t/f+bXN3REREDlVa1AFERCQ1qUBE\nRCQhKhAREUmICkRERBKiAhERkYSoQEREJCEqEBERSYgKREREEqICERGRhKRHHaAt9e7d23NycqKO\nQXV1NV27du3wGULJEUKGUHKEkCGUHCFkCCXHihUrPnL3Pgdd0d3b7S0vL89DsHDhwqgjBJHBPYwc\nIWRwDyNHCBncw8gRQgb3MHIAy70Fv2M1hCUiIglRgYiISEKCKBAz+72ZbTOzNft5fZSZ7TCzVfHb\nncnOmIjS0lJycnIYM2YMOTk5lJaWRh1JRKTVhLIT/QngIeCpA6zzZ3e/MDlxvrrS0lKKioqoqakB\nYMOGDRQVFQFQWFgYZTQRkVYRxBaIu/8J2B51jtZUXFxcXx571dTUUFxcHFEiEZHWZR7IBaXMLAeY\n7e4nNvPaKOAFYDPwN+An7r52P+9TBBQB9O3bN2/q1KltlPjAxowZQ3M/WzNjwYIFSc9TVVVFZmZm\n0j83xBwhZAglRwgZQskRQoZQcowePXqFu5920BVbcqhWMm5ADrBmP691AzLjj8cD5S15zygP483O\nznagyS07OzuSPCEcGugeRo4QMriHkSOEDO5h5Aghg3sYOWhPh/G6e6W7V8Ufvwx8zcx6RxzrgEpK\nSsjIyGi0LCMjg5KSkogSiYi0rpQoEDPrZ2YWfzycutwfR5vqwAoLC5k8eTLHHnssAN27d2fy5Mna\ngS4i7UYQR2GZ2bPAKKC3mW0G7gK+BuDujwKXANeb2R7gc+Dy+GZW0AoLCyksLGTIkCH07NlT5SEi\n7UoQBeLuVxzk9YeoO8w3JY0cOZIpU6awZcuW+i0SEZFUlxJDWKlu5MiRAJSVlUWcRESk9ahAkiA7\nO5shQ4Ywffr0qKOIiLQaFUiSFBQUsGjRIrZvb1fnS4pIB6YCSZJYLEZtbS2zZs2KOoqISKtQgSTJ\naaedRv/+/XnxxRejjiIi0ipUIEliZsRiMebOnUt1dXXUcUREvjIVSBIVFBSwc+dO5syZE3UUEZGv\nTAWSRCNGjKBXr146GktE2gUVSBKlp6czceJEXnrpJXbt2hV1HBGRr0QFkmSxWIwdO3awcOHCqKOI\niHwlKpAky8/PJzMzU8NYIpLyVCBJ1qVLF8aPH09ZWRm1tbVRxxERSZgKJAKxWIytW7fy+uuvRx1F\nRCRhKpAIjB8/ns6dO+ukQhFJaSqQCHTr1o38/HymT5/e7HXTRURSgQokIgUFBaxfv54333wz6igi\nIglRgURkwoQJpKWl6WgsEUlZKpCI9OnTh5EjR2o/iIikLBVIhGKxGGvWrKG8vDzqKCIih0wFEqFY\nLAagrRARSUkqkAhlZWWRl5en/SAikpJUIBErKChg6dKlbNmyJeooIiKHRAUSsb3DWGVlZREnERE5\nNCqQiOXm5jJ06FANY4lIylGBBCAWi7Fo0SK2b98edRQRkRZTgQQgFotRW1vLrFmzoo4iItJiKpAA\nnHbaafTv31+H84pISlGBBMDMiMVizJ07l+rq6qjjiIi0iAokEAUFBezcuZM5c+ZEHUVEpEVUIIEY\nMWIEvXr10tFYIpIyVCCBSE9PZ+LEicyePZtdu3ZFHUdE5KBUIAGJxWJUVlayYMGCqKOIiByUCiQg\n+fn5ZGZm6mgsEUkJQRSImf3ezLaZ2Zr9vF5oZv8Tv/3VzIYlO2MydOnShfHjx1NWVkZtbW3UcURE\nDiiIAgGeAMYd4PUPgG+5+0nAL4DJyQgVhYKCArZu3crrr78edRQRkQMKokDc/U/AfufxcPe/uvsn\n8adLgP5JCRaB888/n86dO2sYS0SCZ+4edQYAzCwHmO3uJx5kvZ8AQ9392v28XgQUAfTt2zdv6tSp\nrZz00FVVVZGZmdni9W+//XY2bNjAM888g5lFkqGthJAjhAyh5AghQyg5QsgQSo7Ro0evcPfTDrqi\nuwdxA3KANQdZZzTwNtCrJe+Zl5fnIVi4cOEhrf/YY4854CtXrowsQ1sJIUcIGdzDyBFCBvcwcoSQ\nwT2MHMByb8Hv2CCGsFrCzE4CHgMmuvvHUedpSxMmTCAtLU0nFYpI0FKiQMwsC5gOXOnu66LO09b6\n9OnDyJEjtR9ERIIWRIGY2bPA68AQM9tsZt8zs382s3+Or3In0Av4rZmtMrPlkYVNklgsxpo1aygv\nL486iohIs4IoEHe/wt2PdvevuXt/d5/i7o+6+6Px16919x7ufnL8dvCdOylu76VutRUiIqEKokCk\nqaysLPLy8rQfRESCpQIJWEFBAUuXLmXLli1RRxERaUIFErC9w1hlZWURJxERaUoFErDc3FyGDh2q\nYSwRCZIKJHCxWIxFixaxfft+Z3oREYmECiRwBQUF1NbWMmvWrKijiIg0ogIJXF5eHgMGDNDhvCIS\nHBVI4MyMSZMmMXfuXKqrq6OOIyJSTwWSAgoKCti5cydz5syJOoqISD0VSAoYMWIEvXr10tFYIhIU\nFUgKSE9PZ+LEicyePZtdu3ZFHUdEBFCBpIxYLEZlZSULFiyIOoqICKACSRn5+flkZmbqaCwRCYYK\nJEV06dKF8ePHU1ZWRm1tbdRxRERUIKmkoKCArVu38vrrr0cdRUREBZJKzj//fDp37qxhLBEJggok\nhXTr1o1zzz2X6dOnU3fdexGR6KhAUkwsFmP9+vW8+eabUUcRkQ5OBZJiJkyYQFpamk4qFJHIqUBS\nTJ8+fRg5cqT2g4hI5FQgKSgWi7FmzRrKy8ujjiIiHZgKJAXtvdSttkJEJEoqkBSUlZVFXl6e9oOI\nSKRUICmqoKCApUuXsmXLlqijiEgHpQJJUXuHsWbMmBFxEhHpqFQgKSo3N5ehQ4dqP4iIREYFksJi\nsRiLFi1i+/btUUcRkQ5IBZLCCgoKqK2tZdasWVFHEZEOSAWSwvLy8hgwYICGsUQkEiqQFGZmxGIx\n5s6dS3V1ddRxRKSDUYGkuFgsxs6dO5kzZ07UUUSkg1GBpLgRI0bQu3dvnVQoIkmnAklx6enpTJgw\ngdmzZ7Nr166o44hIBxJMgZjZODN718zeM7Pbm3k928xeM7P/MbNFZtY/ipwhisViVFZWsmDBgqij\niEgHEkSBmFkn4GHgfOAE4AozO2Gf1X4FPOXuJwH3APclN2W48vPzyczM1NFYIpJUQRQIMBx4z90r\n3H0XMBWYuM86JwCvxR8vbOb1DqtLly6MHz+eGTNmUFtbG3UcEekgQimQY4FNDZ5vji9r6E3g4vjj\nGHCEmfVKQraUUFBQwLZt23j99dejjiIiHYS5e9QZMLNvA2Pd/dr48yuB4e7+wwbrHAM8BAwE/kRd\nmXzd3Xfs815FQBFA375986ZOnZqcb+IAqqqqyMzMbNPPqK6uJhaLMWnSJG644YZIMrRECDlCyBBK\njhAyhJIjhAyh5Bg9evQKdz/toCu6e+Q34AxgboPndwB3HGD9TGDzwd43Ly/PQ7Bw4cKkfM4FF1zg\nOTk5/uWXX0aW4WBCyBFCBvcwcoSQwT2MHCFkcA8jB7DcW/C7O5QhrDeA481soJl1Bi4HZjZcwcx6\nm9nevHcAv09yxuDFYjHWr1/Pm2++GXUUEekAgigQd98D/ACYC7wN/NHd15rZPWY2Ib7aKOBdM1sH\n9AVKIgkbsAkTJpCWlqaTCkUkKYIoEAB3f9ndB7v7ce5eEl92p7vPjD9+3t2Pj69zrbt/EW3i8PTp\n04eRI0fqcF4RSYpgCkRaR0FBAWvWrKG8vDzqKCLSzqlA2plJkyYBaCtERNqcCqSdycrKIi8vT/tB\nRKTNqUDaoYKCApYuXcqWLVuijiIi7VirFoiZdTOzbq35nnLoYrEYADNmzIg4iYi0Z61WIGZWAnwK\nfGJmG82szMzuanAYriRJbm4uQ4cO1X4QEWlTrbkFciNwFtAD+CdgAfAPwC9a8TOkhWKxGIsWLWL7\n9u1RRxGRdqo1C+QTYJm7V7r7n9z9N+5+lbsPa8XPkBYqKCigtraWWbNmRR1FRNqp1iyQycA/tuL7\nyVeQl5fHgAEDNIwlIm2mNQvkMuBRM/sPMxtjZt1b8b3lEJkZsViMuXPnUl1dHXUcEWmHWrNAfk7d\nVQNzgCeBj83sfTOb1oqfIYcgFouxc+dO5syZE3UUEWmHDrlAzOza5pa7+wx3v8vdJ7j7AOomPLye\nupl2JQIjRoygd+/eOqlQRNpEIlsg/2FmJx9oBTM70t0/cvd57v7vCWaTryg9PZ0JEyYwe/Zsdu/e\nHXUcEWlnEimQecAL+9vHYWY5wF+/QiZpRT169KCyspLzzjuPnJwcSktLo44kIu1EIgVyVfz+qX1f\nMLPTgSVAz68SSlpHaWkpjzzySP3zDRs2UFRUpBIRkVZxyAXiddcgvxQ418yK9y43s0nAQuBj6i5R\nKxErLi6mpqam0bKamhqKi4v38xUiIi2XnsgXufsKM/sX6vaHvA6cRN0RWIuAi+MlIxHbuHHjIS0X\nETkUBy0QM3sTWAnsvV/l7p+6+yNmNgKYDRwGPAF8P355WglAVlYWGzZsaLL86KOPjiCNiLQ3LRnC\n2kPdSYK/Bl6j7vyO9Wb2IvA3oBPwW3f/nsojLCUlJWRkZDRZvmvXLjZv3hxBIhFpTw5aIO6eB2QC\nw4DvAr8BKoBvAbcAXwNuMLOtZjbHzO4zs2+3XWRpqcLCQiZPnkx2djZmRnZ2Nvfccw9ffPEF+fn5\nbN26NeqIIpLCWrQPxN1rgdXx2x/2Lo8fsnsycEqD+/MAB3QGegAKCwspLCxk0aJFjBo1CoDRo0cz\nduxYzjvvPBYuXEjPnjpoTkQO3VeaysTd1zc4A32iu2cBvakrEQnUiBEjKCsr45133mHcuHFUVlZG\nHUlEUlCrX9LW3be7+2ut/b7SuvLz83n++edZuXIlF154YZPDfUVEDkbXRO/ALrroIp5++mn+8pe/\nEIvF+OKLL6KOJCIpRAXSwV122WU89thjzJs3j8suu0xzZolIi6lAhKuvvpqHHnqIsrIyrrrqKmpr\na6OOJCIpIKEz0aX9ufHGG6murua2224jIyODyZMnk5am/1+IyP6pQKTerbfeSlVVFb/4xS/o2rUr\nDz74IGYWdSwRCZQKRBq5++67qaqq4oEHHiAzM5OSkpKoI4lIoFQg0oiZ8etf/5qamhruvfdeunbt\nys9+9rOoY4lIgFQg0oSZ8dvf/pbq6mqKi4vp2rUrN910U9SxRCQwKhBpVlpaGo8//jg1NTXcfPPN\ndO3alWuvvTbqWCISEB1mI/uVnp7Os88+y/nnn09RURHPPPNM1JFEJCDBFIiZjTOzd83sPTO7fT/r\nXGpmb5nZWjPTb7Mk6Ny5My+88ALf+ta3+M53vsOLL74YdSQRCUQQBWJmnYCHgfOBE4ArzOyEfdY5\nHrgDOMvdvw7cnPSgHdThhx/OzJkzOf3007nsssuYM2dO1JFEJABBFAgwHHjP3SvcfRcwFZi4zzrX\nAQ+7+ycA7r4tyRk7tCOOOIJXXnmFr3/968RiMRYvXhx1JBGJWCgFciywqcHzzfFlDQ0GBpvZX8xs\niZmNS1o6AaB79+7MmzePgQMHcuGFF7J06dKoI4lIhMzdo85A/AqGY9392vjzK4Hh7v7DBuvMBnYD\nlwL9gT8DJ7r7p/u8VxFQBNC3b9+8qVOnJuebOICqqioyMzPbTYaPPvqIm266icrKSh544AEGDRoU\nSY5EhZAhlBwhZAglRwgZQskxevToFe5+2kFXdPfIb8AZwNwGz+8A7thnnUeB7zZ4/hpw+oHeNy8v\nz0OwcOHCqCO0eoYPPvjABwwY4H369PG33norshyJCCGDexg5QsjgHkaOEDK4h5EDWO4t+N0dyhDW\nG8DxZjbQzDoDlwMz91lnBjAawMx6UzekVZHUlFIvJyeH+fPnk5aWRn5+PhUV+qMQ6WiCKBB33wP8\nAJgLvA380d3Xmtk9ZjYhvtpc4GMzewtYCPzU3T+OJrEADB48mPnz57Nz507OOeccNm/eHHUkEUmi\nIAoEwN1fdvfB7n6cu5fEl93p7jPjj93d/8XdT3D3b7h79Ds3hBNPPJF58+axfft2zjnnHLZu3Rp1\nJBFJkmAKRFJXXl4eL730Eps3b+bcc89l+/btUUcSkSRQgUirGDFiBGVlZaxbt45x48ZRWVkZdSQR\naWMqEGk1+fn5TJs2jZUrV3LBBRdQXV0ddSQRaUMqEGlVF110EaWlpfz1r38lFovxxRdfRB1JRNqI\nCkRa3aWXXsqUKVN49dVXufTSS9m9e3fUkUSkDahApE1897vf5aGHHmLmzJl85zvfoba2NupIItLK\ndEEpaTM33ngj1dXV3HbbbXz44Yd88MEHbNq0iaysLEpKSigsLIw6ooh8BSoQaVO33norS5YsaXQd\nkQ0bNlBUVASgEhFJYRrCkja3YsWKJstqamooLi6OII2ItBYViLS5TZs2Nbt848aNSU4iIq1JBSJt\nLisrq9nl7k5hYSHvvfdekhOJSGtQgUibKykpISMjo9Gyww8/nIsuuogXX3yRoUOH8v3vf1+TMYqk\nGBWItLnCwkImT55MdnY2ZkZ2djb/9V//xcyZM3n//fe5/vrrefzxxxk0aBC33HILH330UdSRRaQF\nVCCSFIWFhaxfv54FCxawfv36+qOvjj76aP7zP/+TdevWccUVV/Dggw8ycOBAfv7zn2s+LZHAqUAk\nCDk5OTz++OOsWbOGcePGcffddzNw4EB+9atf8fnnn0cdT0SaoQKRoOTm5jJt2jSWL1/O8OHD+elP\nf8qgQYN49NFH2bVrV9TxRKQBFYgEKS8vj1deeYXFixczcOBArr/+enJzc3n66ac1LYpIIFQgErSz\nzz6bP//5z7z00kt069aNK6+8kmHDhjFjxgzcPep4Ih2aCkSCZ2aMHz+eFStW8Nxzz7F7925isRjf\n/OY3mT9/ftTxRDosFYikjLS0NC699FLWrl3LlClT+PDDDzn33HM555xzWLJkSdTxRDocFYiknPT0\ndK655hrKy8v5zW9+w+rVqznjjDOYOHEiq1evjjqeSIehApGUddhhh/GjH/2IiooKfvnLX7J48WKG\nDRum6VFEkkQFIikvMzOT4uJiKioquO222zQ9ikiSqECk3ejZsyf33XcfFRUVzU6PUlpaSk5ODmPG\njCEnJ4fS0tKoI4ukNBWItDv9+vVrMj1K//79ufrqq9mwYQPuXn9RK5WISOJUINJuNZwepVOnTuze\nvbvR67qolchXowKRdi83N3e/82lt2LCBJUuW8OWXXyY5lUjqU4FIh7C/i1oBnHHGGfTv35/rr7+e\nefPmac4tkRZSgUiH0NxFrTIyMvjd737H008/zZlnnslTTz3F2LFjOeqooygsLGTatGl89tlnESUW\nCV961AFEkmHv9UeKi4vZuHEjWVlZlJSU1C8vLCzk888/Z/78+cyYMYOZM2fyzDPPcNhhh5Gfn8+k\nSZOYMGECRx11VJTfhkhQtAUiHcb+Lmq1197L7E6ZMoW///3vLF68mBtuuIG1a9dy3XXX0a9fP0aO\nHMn9999PRUVFRN+FSDhUICLN6NSpE2effXZ9WaxatYo777yTzz77jFtuuYXjjjuOYcOGcdddd7Fq\n1SrNDCwdkgpE5CDMjGHDhvHzn/+cVatWUVFRwf3330/37t355S9/ySmnnMLAgQO5+eabWbx4MXv2\n7Ik6skhSBFMgZjbOzN41s/fM7PZmXv9nM1ttZqvM7P+a2QlR5BQZOHAgP/7xj1m8eDEffvghU6ZM\n4Rvf+AaPPvooo0aNol+/flxzzTXMnDlTl+OVdi2IAjGzTsDDwPnACcAVzRTEM+7+DXc/Gfh34P4k\nxxRp4qijjuKaa65h1qxZfPTRR0ybNo1x48Yxffp0Jk6cSO/evbn44ov5wx/+wCeffAKgKVWk3Qjl\nKKzhwHvuXgFgZlOBicBbe1dw98oG63cFNOgsQcnMzOSSSy7hkksuYdeuXSxevJgZM2YwY8YMpk+f\nTnp6OoMHD6a8vLz+rPi9U6oATXbqi4QuiC0Q4FhgU4Pnm+PLGjGzG83sfeq2QH6UpGwih6xz586c\ne+65PPzww2zatImlS5fyk5/8hHXr1jU7pcott9xCdXV1RGlFEmMhHD1iZt8Gxrr7tfHnVwLD3f2H\n+1n/H+PrX9XMa0VAEUDfvn3zpk6d2nbBW6iqqorMzMwOnyGUHFFmGDNmzH6P2EpLSyMnJ4fc3FyG\nDBlCbm4uOTk5pKe33UBBCH8eoeQIIUMoOUaPHr3C3U876IruHvkNOAOY2+D5HcAdB1g/DdhxsPfN\ny8vzECxcuDDqCEFkcA8jR5QZsrOznbrh10a3Pn36+L/+67/62LFjvWfPnvXLu3Tp4meeeabffPPN\nXlpa6uXl5f7ll1+2Wp4Q/jzcw8gRQgb3MHIAy70Fv7tD2QfyBnC8mQ0EtgCXA//YcAUzO97dy+NP\nLwDKEUkxJSUlFBUVUVNTU78sIyODBx54oH4fiLtTUVHBsmXLeOONN1i2bBm/+93vePDBBwHo0aMH\np59+OsOHD2f48OGcfvrp9OvXL5LvRzq2IArE3feY2Q+AuUAn4PfuvtbM7qGuCWcCPzCzfGA38AnQ\nZPhKJHQHm1IF6s47Oe644zjuuOO44oorANizZw9r165tVCr33XcftbW1AAwYMKC+TIYPH05eXh7d\nunVL/jcoHUoQBQLg7i8DL++z7M4Gj29KeiiRNlBYWEhhYSGLFi1i1KhRLfqa9PR0hg0bxrBhw7ju\nuuuAup3vK1eubFQqL7zwAlBXQkOHDm1UKieddBKHHXYYUHco8YFKTKQlgikQETk0GRkZnHXWWZx1\n1ln1yz7++GOWL1/OsmXLWLZsGa+88gpPPvkkUHdk2LBhw+jevTuLFy+un7ZehxJLolQgIu1Ir169\nGDt2LGPHjgXq9qds2rSpfgtl2bJlzJ8/v8mRYDU1Ndxwww3U1tYyePBgBg8eTM+ePaP4FiSFqEBE\n2jEzIysri6ysLC6++GKg7nDh5lRWVnLVVf9/12Lv3r3ry2TIkCH1jwcNGkSXLl2Skl/CpgIR6WCy\nsrLYsGFDs8tfffVV1q1bx7vvvsu6detYt24dc+fO5Yknnqhfz8zIzs5utlyysrL2W1DS/qhARDqY\n/R1KfO+999YXwYUXXtjoaz777DPKy8sbFcu7777Lk08+2eiqjYcddhjHH398s+XSu3fvJlm0Mz+1\nqUBEOpiWHEq8ryOOOIJTTz2VU089tdFyd2fr1q1NtlrWrl3LzJkzG01t37Nnz0bFsnXrViZPnszO\nnTsB7cxPRSoQkQ4okUOJm2Nm9OvXj379+nH22Wc3em3Pnj2sX7++yVbLa6+9xlNPPdXs+9XU1HD9\n9dezefNmBgwYUH875phj6Ny5c8I5pW2oQESkTaSnpzNo0CAGDRrEBRdc0Oi1qqoqunXr1uy8YJ99\n9hm33974kkBmRt++fRkwYAD9+/dvVC57nx9zzDEJzRumYbTEqUBEJOkyMzP3uzM/Ozub1atXs3nz\nZjZt2sSmTZsaPX7nnXd49dVXqaqqavR1aWlp9OvXr9ly2Xvr168fnTp1qv+a0tLSRvuDNIx2aFQg\nIhKJ/e3MLykp4YgjjiA3N5fc3Nz9fv2OHTualMve56tXr+bll19u9N5Qd637Y445pr5cXnnllSbr\n1NTUUFxcrAJpARWIiEQikZ35DR155JEceeSRnHjiic2+7u58+umnTcpl7+OVK1c2OoKsoQ0bNjBw\n4ED69u3LUUcdVX/f3OOePXs22qpJVCoOpalARCQyrbUzvzlmRo8ePejRowcnnXRSs+tkZ2ezcePG\nJsuPOOIIRowYwbZt29i4cSPLly9n27Zt9ZNXNpSWlkbv3r0PWDIN7w8//PAm75GqQ2kqEBHpsO69\n995mh9EeeeSRJr+4v/zySz755BO2bdvGtm3b2Lp1a6P7vY+XLVvGtm3b9rt1k5mZ2aRgnnvuuZQc\nSlOBiEiHdSjDaGlpafTq1YtevXodcN/MXp9//nmTstn3cUVFBUuWLGHHjh3NvkdzW0chUYGISIfW\nVsNohx9+ONnZ2WRnZx903f0NpWVlZbVanragSWtERCJ27733kpGR0WjZ3iPSQqYCERGJWGFhIZMn\nTyY7O7t+ssrJkycHvf8DNIQlIhKEtjwira1oC0RERBKiAhERkYSoQEREJCEqEBERSYgKREREEqIC\nERGRhKhAREQkISoQERFJiApEREQSogIREZGEqEBERCQhKhAREUmICkRERBKiAhERkYSoQEREJCHB\nFIiZjTOzd83sPTO7vZnXDzOz5+KvLzWznOSnFBGRvYIoEDPrBDwMnA+cAFxhZifss9r3gE/cfRDw\nAPB/kptSREQaCqJAgOHAe+5e4e67gKnAxH3WmQg8GX/8PHCOmVkSM4qISAOhFMixwKYGzzfHlzW7\njrvvAXYAvZKSTkREmgjlmujNbUl4AutgZkVAUfxplZm9+xWztYbewEfKAISRI4QMEEaOEDJAGDlC\nyABh5MhuyUqhFMhmYECD5/2Bv+1nnc1mlg4cCWzf943cfTIwuY1yJsTMlrv7aR09Qyg5QsgQSo4Q\nMoSSI4QMIeVoiVCGsN4AjjezgWbWGbgcmLnPOjOBq+KPLwEWuHuTLRAREUmOILZA3H2Pmf0AmAt0\nAn7v7mvN7B5gubvPBKYAfzCz96jb8rg8usQiIhJEgQC4+8vAy/ssu7PB453At5Odq5WEMKQWQgYI\nI0cIGSCMHCFkgDByhJABwslxUKZRIBERSUQo+0BERCTFqEDaiJkNMLOFZva2ma01s5siytHFzJaZ\n2ZvxHHdHkSOepZOZrTSz2RGeVS1SAAAE70lEQVRmWG9mq81slZktjyhDdzN73szeif/9OCOCDEPi\nP4O9t0ozuzmCHD+O/71cY2bPmlmXCDLcFP/8tcn8GZjZ781sm5mtabCsp5m9ambl8fseycqTCBVI\n29kD3OLuucA3gRubmZ4lGb4Axrj7MOBkYJyZfTOCHAA3AW9H9NkNjXb3kyM8VPI3wBx3HwoMI4Kf\nibu/G/8ZnAzkATXAi8nMYGbHAj8CTnP3E6k7gCapB8eY2YnAddTNhjEMuNDMjk/Sxz8BjNtn2e3A\na+5+PPBa/HmwVCBtxN0/dPf/jj/+jLpfEvueXZ+MHO7uVfGnX4vfkr7jy8z6AxcAjyX7s0NiZt2A\ns6k7qhB33+Xun0abinOA9919QwSfnQ4cHj+3K4Om53+1tVxgibvXxGe4WAzEkvHB7v4nmp7L1nDK\npieBScnIkigVSBLEZw4+BVga0ed3MrNVwDbgVXePIseDwK3AlxF8dkMOzDOzFfFZC5LtH4D/BR6P\nD+c9ZmZdI8jR0OXAs8n+UHffAvwK2Ah8COxw93lJjrEGONvMeplZBjCexic1J1tfd/8Q6v4TChwV\nYZaDUoG0MTPLBF4Abnb3yigyuHttfKiiPzA8vtmeNGZ2IbDN3Vck83P34yx3P5W6mZ9vNLOzk/z5\n6cCpwCPufgpQTYTDFPETdycA0yL47B7U/Y97IHAM0NXM/imZGdz9bepm9n4VmAO8Sd3ws7SACqQN\nmdnXqCuPUnefHnWe+FDJIpqOu7a1s4AJZraeupmWx5jZ00nOAIC7/y1+v426Mf/hSY6wGdjcYCvw\neeoKJSrnA//t7lsj+Ox84AN3/1933w1MB85Mdgh3n+Lup7r72dQNKZUnO0MDW83saID4/bYIsxyU\nCqSNxKeanwK87e73R5ijj5l1jz8+nLp/tO8kM4O73+Hu/d09h7rhkgXuntT/aQKYWVczO2LvY+A8\n6oYwksbd/w5sMrMh8UXnAG8lM8M+riCC4au4jcA3zSwj/u/lHCI4oMDMjorfZwEFRPfzgMZTNl0F\nlEWY5aCCORO9HToLuBJYHd//APCz+Bn3yXQ08GT8ol1pwB/dPbLDaCPWF3gxfhmZdOAZd58TQY4f\nAqXx4aMK4OoIMhAf8z8X+H4Un+/uS83seeC/qRs2Wkk0Z2G/YGa9gN3Aje7+STI+1MyeBUYBvc1s\nM3AX8G/AH83se9QVbNCzb+hMdBERSYiGsEREJCEqEBERSYgKREREEqICERGRhKhAREQkISoQERFJ\niApEREQSogIREZGEqEBERCQhKhCRJDKzX5uZm1mWmf2bmX1gZp/Hp5cfEXU+kUOhqUxEksjMXqPu\nCoBbqJtEcRHQB/gJdVO794/PTCsSPE2mKJJcJwNHAj9y96f2Loxfka8YyCHa6cRFWkxDWCJJYmbZ\nQE/gpYblEfdF/P7z5KYSSZwKRCR5TonfT23mtROBz6gb2hJJCSoQkeTZWyBLmnktD1jp2ikpKUQF\nIpI8JwM7gPcbLoxfMfIfqLuwkkjKUIGIJM8p1F1/fN+tjFMBQwUiKUYFIpIE8UumDgBWNPPyqfF7\nFYikFBWISHLs3f/RXEnkATXAO8mLI/LV6URCERFJiLZAREQkISoQERFJiApEREQSogIREZGEqEBE\nRCQhKhAREUmICkRERBKiAhERkYSoQEREJCEqEBERScj/A0rSf0KgIfYbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110e91e10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rad, den = np.random.rand(1)+1., np.random.rand(1)+1.\n",
    "\n",
    "# Initialize the variables for the Love numbers and the degree\n",
    "kn = []\n",
    "degree = range(2,11)\n",
    "\n",
    "# Calculate the Love numbers\n",
    "for n in degree:\n",
    "    kn.append(knf.calc_kn(n,rad,den))\n",
    "\n",
    "# Print out results\n",
    "print(\"n   k_n\")\n",
    "for n,k in zip(degree,kn): print(\"%02i  %5.4f\"%(n,float(k)))\n",
    "\n",
    "# Create figure\n",
    "nticks = [i for i in degree]\n",
    "nlabel = [str(i) for i in nticks]\n",
    "kticks = [i*0.300 for i in range(6)]\n",
    "klabel = [\"%2.1f\"%i for i in kticks]\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "ax.plot(degree,kn,'k-o')\n",
    "ax.set_ylim([0,1.55])\n",
    "ax.set_xlim([min(degree)-1,max(degree)+1])\n",
    "ax.set_xlabel(r\"$n$\",fontsize = 18)\n",
    "ax.set_ylabel(r\"$k_n$\",fontsize = 18)\n",
    "ax.xaxis.set_ticks(nticks)\n",
    "ax.xaxis.set_ticklabels(nlabel)\n",
    "ax.yaxis.set_ticks(kticks)\n",
    "ax.yaxis.set_ticklabels(klabel)\n",
    "plt.grid()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Earth Fluid Love number $k_2$ from PREM density profile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PREM data from the IRIS database\n",
    "The data can be loaded directly from the IRIS website using pandas. However, pandas is not (yet) part of every python installation. Thus, the file, located @http://ds.iris.edu/files/products/emc/data/PREM/PREM_1s.csv is included in the repository. The pandas snippet is commented out below.\n",
    "\n",
    "The PREM contains surface-to-center arrays. Before calculating the Love numbers one needs to flip them and remove the central point, given that the radial vector should contain the outer boundary of each shell. (See also the README.md in the repository.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "## ---------------------------- DOWNLOAD PREM USING PANDAS --------------------------------------------\n",
    "## Since the file is provided without header, the column names are supplied using information available\n",
    "## in the relevan IRIS webpage (http://ds.iris.edu/ds/products/emc-prem/)\n",
    "\n",
    "#import pandas as pd\n",
    "#prem_url = \"http://ds.iris.edu/files/products/emc/data/PREM/PREM_1s.csv\"\n",
    "#prem_data = pd.read_csv(prem_url,header=None)\n",
    "#prem_data.columns = ['radius','depth','density','Vpv','Vph','Vsv','Vsh','eta','Q-mu','Q-kappa']\n",
    "#rad = np.flip(np.asarray(prem_data.radius),0)[1:]\n",
    "#den = np.flip(np.asarray(prem_data.density),0)[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "## ---------------------- READ PREM FROM THE FILE IN THE REPOSITORY -------------------------------------\n",
    "prem_file = 'PREM_1s.csv'\n",
    "rad,den = np.loadtxt(prem_file,unpack=True,usecols=(0,2),delimiter=',')\n",
    "rad = np.flip(rad,0)[1:]\n",
    "den = np.flip(den,0)[1:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute and print values for $k_2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Earth fluid Love number k_2: 0.9346387\n"
     ]
    }
   ],
   "source": [
    "print(\"Earth fluid Love number k_2: %8.7f\"%(float(knf.calc_kn(2,rad,den))))"
   ]
  },
  {
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}